R繪圖應用實例:單因素方差分析ANOVA及繪圖

本文主要是利用日常實驗數據，嘗試用R進行單因素方差分析並繪製柱形圖。

ANOVA原理參考：單因素方差分析（One-way Anova）
實驗數據：
在隨機劃分的試驗田中，施加三種複合肥（B,C,D），飼料填充物做空白對照（A），一段時間後，測定試驗田內植株高度，比較數據有無不同，差異性是否具有統計學意義。
計算結果錄入Excel，如下：

圖片.png

（一）：數據準備

rm(list = ls())
# prepare the data
data <- data.frame(A=c(47,65,44,59,62,37,51),
B=c(68,55,49,62,70,59,63),
C=c(78,76,72,81,76,71,83),
D=c(85,65,81,98,75,92,79))
head(data)

# convert data from a wide format to a long format.
library(reshape2)
data_long <- melt(data, measure.vars= c("A","B","C","D"))
head(data_long)

（二）：方差分析（條件適用檢驗及方差分析）

# 單因素方差分析（One-way Anova）前提假設的檢驗：
# 1. 正態性檢驗(對4個水平下的每組數據都做一次正態檢驗)
shapiro.test(data$A) #或shapiro.test(data_long$value[1:7])
shapiro.test(data$B) #或shapiro.test(data_long$value[8:14])
shapiro.test(data$C) #或shapiro.test(data_long$value[15:21])
shapiro.test(data$D) #或shapiro.test(data_long$value[22:28])
############################################################################
#split:拆分為list數據集; #
#lapply:對數據集進行shapiro.test檢驗循環; #
#unlist():向量的原樣返回 #
group_data <- split(data_long[,2], data_long[,1]) #
group_data #
unlist(lapply(group_data, function(x){ #
shapiro.test(x)$p.value #
})) #
############################################################################
# 使用Q-Q圖來檢驗正態性
library(car)
qqPlot(group_data[[1]]) #或qqPlot(data$A)
qqPlot(group_data[[2]]) #或qqPlot(data$A)
qqPlot(group_data[[3]]) #或qqPlot(data$A)
qqPlot(group_data[[4]]) #或qqPlot(data$D)

# 2. 方差齊性檢驗，使用car包的leveneTest()
leveneTest(value~variable, data = data_long)

# 3. 離群點檢驗：
outlierTest(lm(value~variable, data = data_long))

Q-Q圖示例

# 4. 單因素方差分析ANOVA
aov1 <- aov(value~variable, data = data_long)
summary(aov1)
Df Sum Sq Mean Sq F value Pr(>F)
variable 3 4049 1349.7 18.2 2.25e-06 ***
Residuals 24 1780 74.2
---
Signif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1

# Show the means
model.tables(aov1, "means")

# Tukey HSD post-hoc test
TukeyHSD(aov1)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = value ~ variable, data = data_long)
$variable
diff lwr upr p adj
B-A 8.714286 -3.984450 21.41302 0.2574525
C-A 24.571429 11.872693 37.27016 0.0000986
D-A 30.000000 17.301265 42.69874 0.0000055
C-B 15.857143 3.158407 28.55588 0.0106071
D-B 21.285714 8.586979 33.98445 0.0005874
D-C 5.428571 -7.270164 18.12731 0.6453377

# Alternative: Multiple comparisons using multcomp package
library(multcomp)
e <- glht(aov1, linfct = mcp(variable = "Tukey"))
summary(e)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrast
Fit: aov(formula = value ~ variable, data = data_long)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
B - A == 0 8.714 4.603 1.893 0.2573
C - A == 0 24.571 4.603 5.338 <0.001 ***
D - A == 0 30.000 4.603 6.517 <0.001 ***
C - B == 0 15.857 4.603 3.445 0.0106 *
D - B == 0 21.286 4.603 4.624 <0.001 ***
D - C == 0 5.429 4.603 1.179 0.6453
---
Signif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1
(Adjusted p values reported -- single-step method)

（三）：統計分析

# summarySE 函數提供了標準差、標準誤以及95%的置信區間
library(Rmisc)
data_long_count <- summarySE(data_long, measurevar="value",
groupvars= "variable")
data_long_count

（四）：繪圖

繪圖應用1

# 繪製帶有顯著性標記的條形圖
marker <- c("", "", "***", "***")
library(ggplot2)
ggplot(data_long_count, aes(x=variable, y=value,fill=variable)) +
geom_bar(stat="identity", color="black",size=.3) +
geom_errorbar(aes(ymin=value-se, ymax=value+se),
size=.3, width=.2) +
geom_text(aes(y = value + 1.5 * se, label = marker),
size = 5, fontface = "bold") +
scale_color_brewer("Set1") +
ggtitle("Effects of Three Compound Fertilizers on Plant Height") +
xlab("Different compound fertilizer treatments") +
ylab("Plant height")

Rplot04.png

繪圖應用2

# 最終繪圖效果
library(ggsignif)
library(ggplot2)
ggplot(data_long_count, aes(x=variable, y=value,fill=variable)) +
geom_bar(stat="identity", color="black",size=.3) +
# 添加誤差棒
geom_errorbar(aes(ymin=value-se, ymax=value+se), size=.3, width=.2) +
# 手動添加顯著性
geom_signif(annotations = c("***","***"), y_position = c(83, 90),
xmin = c(1, 1), xmax = c(3, 4),
tip_length = c(c(0.65, 0.05),c(0.7, 0.05)), vjust = 0 ) +
# 設置色板
scale_color_brewer(palette = "Set1") +
# 更改主標題
ggtitle("Effects of Three Compound Fertilizers on Plant Height") +
xlab("Different compound fertilizer treatments") +
ylab("Plant height") +
# 更改軸標籤的外觀
theme(plot.title = element_text(color="red", size=14, face="bold.italic"),
axis.title.x = element_text(color="blue", size=14, face="bold"),
axis.title.y = element_text(color="#993333", size=14, face="bold")) +
# 更改圖例標籤
labs(fill="Treated")

Rplot05.png