对于两组数据间的差异分析，最常见的方法就是使用T检验比较两组均值是否存在显著不同。当拓展到多组（三组及以上）时，使用T检验逐一两两比较的方法无疑是低效的，不仅仅由于需要的检验次数增多，而且发生I型错误（拒绝真）的概率也会增大。Fisher提出一种广义T检验的方法来比较三组及以上总体的均值，称为方差分析（ANOVA）。
几种常见的ANOVA包含单因素方差分析（单因素ANOVA）、单因素协方差分析（ANCOVA）、双因素方差分析（双因素ANOVA）、重复测量方差分析（重复测量ANOVA）、多元方差分析（MANOVA）等。本篇首先介绍其中最常涉及的单因素ANOVA在R语言中的实现过程，一组因子变量对应一组因变量
数据预处理
读入文件
#读入文件
soil <- read.table('soil.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
group <- read.table('group.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
soil <- merge(soil, group, by = 'sample')
#假设样本采自三种土壤环境（A、B、C），我们比较三种土壤环境下的细菌群落的 chao1 指数是否存在显著差异
#以 chao1 指数为例，同时将分组列转换为因子变量
chao1 <- soil[ ,c('sample', 'site', 'chao1')]
chao1$site <- factor(chao1$site)
str(chao1)
head(chao1)
sample site chao1
1 A1_1 A 2711.889
2 A1_2 A 3235.087
3 A1_3 A 3262.864
4 A2_1 A 3720.159
5 A2_2 A 3021.062
6 A2_3 A 3435.125
评估检验的假设条件
与T检验相似，ANOVA同样要求数据服从正态分布；此外，ANOVA还建立在各组方差相等的基础上。因此，在执行单因素ANOVA之前，我们首先应当对数据进行正态性分布验证，以及方差齐性检验。
正态性检验
#QQ-plot 检查数据是否符合正态分布（所有的点都离直线很近，落在置信区间内说明正态性良好）
library(car)
qqPlot(lm(chao1~site, data = chao1), simulate = TRUE, main = 'QQ Plot', labels = FALSE)
##Shapiro-Wilk 检验，当且仅当三者 p 值均大于 0.05 时表明数据符合正态分布，从中看出B组p值小于0.05，不符合正态分布（但是由于学术不精，暂且就当它符合正态分布吧）
shapiro <- tapply(chao1$chao1, chao1$site, shapiro.test)
shapiro
$A
Shapiro-Wilk normality test
data: X[[i]]
W = 0.96512, p-value = 0.8536
$B
Shapiro-Wilk normality test
data: X[[i]]
W = 0.83571, p-value = 0.02456
$C
Shapiro-Wilk normality test
data: X[[i]]
W = 0.91612, p-value = 0.2554
方差齐性检验
R语言中提供了一些可用来做方差齐性检验的函数，例如Bartlett检验（bartlett.test）、Fligner-Killeen检验（fligner.test()）、Brown-Forsythe检验（HH包hov()）等。
对于已经通过正态性检验的数据，推荐使用Bartlett检验来进行方差齐性检验（它建立在数据分布正态性的前提下，如果数据服从正态分布，这是最好的检验方法）；Fligner-Killeen检验是一个非参数检验，通常在数据偏离正态性时使用（当然，如果数据已经偏离正态分布了，也没必要再继续了，所以Fligner-Killeen检验似乎并不能很好地适用在方差分析过程中）。
#使用 Bartlett 检验进行方差齐性检验（p 值大于 0.05 说明方差齐整）
bartlett.test(chao1~site, data = chao1)
Bartlett test of homogeneity of variances
data: chao1 by site
Bartlett's K-squared = 5.9425, df = 2, p-value = 0.05124
单因素方差分析
单因素ANOVA
我们的数据通过了正态性检验和方差齐性检验，接下来进行单因素ANOVA。R语言执行方差分析的命令是aov()，对于单因素方差分析，aov()函数书写为aov(y~A)的样式，A即为因子变量。
如果不满足上述前提假设，一是可以考虑转化数据（当然，我们需要确保转换后的数据能够被合理解释，否则将无意义），二是可以考虑使用非参数的检验方法，对于单因素的分析，可选的非参数替代方法例如Kruskal-Wallis检验（kruskal.test()）、Friedman检验（friedman.test()）等。
#满足假设，单因素方差分析，详情使用?aov查看帮助
fit <- aov(chao1~site, data = chao1)
summary(fit)
Df Sum Sq Mean Sq F value Pr(>F)
site 2 4635658 2317829 17.48 6.66e-06 ***
Residuals 33 4376131 132610
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#若想查看各组均值及标准差，可使用 aggregate()
chao1_mean <- aggregate(chao1$chao1, by = list(chao1$site), FUN = mean)
Group.1 x
1 A 3340.115
2 B 2608.779
3 C 2552.170
chao1_sd <- aggregate(chao1$chao1, by = list(chao1$site), FUN = sd)
Group.1 x
1 A 301.8213
2 B 497.8551
3 C 242.6404
多重比较
上述单因素ANOVA告诉我们3种土壤环境下的细菌群落的Chao1指数具有显著差异，这种差异是在整体水平而言的，并没有告诉我们究竟谁和谁存在差异。如果我们想继续获知两两分组之间的差异，进行多重比较即可。
常用Tukey HSD检验，在ANOVA结果的基础上继续执行事后两两比较。不推荐使用T检验（注意T检验和Tukey检验是两回事），原因正如本文开始时所提，多次T检验容易提高I型错误的概率。
##方差分析后，多重比较，继续探寻两两分组间的差异
#Tukey HSD 检验
tuk <- TukeyHSD(fit, conf.level = 0.95)
plot(tuk)
tuk
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = chao1 ~ site, data = chao1)
$site
diff lwr upr p adj
B-A -731.33642 -1096.1330 -366.5398 0.0000681
C-A -787.94475 -1152.7413 -423.1482 0.0000223
C-B -56.60833 -421.4049 308.1882 0.9233767
显著水平默认为0.05。Tukey检验显示，A组和B组、A组和C组存在显著差异，但B组和C组无差异。（根据文字部分p值判断；或者根据图片判断，未越过虚线则表示无差异）
multcomp包中提供了更直观的方法，展示Tukey检验的结果。
library(multcomp)
tuk <- glht(fit, alternative = 'two.sided', linfct = mcp(site = 'Tukey'))
plot(cld(tuk, level = 0.05, decreasing = TRUE))
同样地，显著水平默认为0.05。结果以箱线图的方式，直观地为我们展示出组间差异。从图中我们可以轻易得知，A组（A环境下的土壤细菌群落）的Chao1指数显著高于其它两组（B、C环境下的土壤细菌群落），同时B、C二者无差异。
ggplot2柱状图示例
dat <- merge(chao1_mean, chao1_sd, by = 'Group.1')
names(dat) <- c('group', 'mean', 'sd')
dat <- cbind(dat, sign = c('a', 'b', 'b'))
library(ggplot2)
ggplot(dat, aes(group, mean)) +
geom_col(aes(fill = group), width = 0.4, show.legend = FALSE) +
geom_errorbar(aes(ymax = mean + sd, ymin = mean - sd), width = 0.15, size = 0.5) +
geom_text(aes(label = sign, y = mean +sd + 200)) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), plot.title = element_text(hjust = 0.5)) +
labs(x = 'Group', y = 'Chao1', title = 'Tukey HSD test')
来源：小白鱼的生统笔记
总代码
##单因素方差分析（单因素 ANOVA）
#读入文件
soil <- read.table('soil.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
group <- read.table('group.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
soil <- merge(soil, group, by = 'sample')
#假设样本采自三种土壤环境（A、B、C），我们比较三种土壤环境下的细菌群落的 chao1 指数是否存在显著差异
#以 chao1 指数为例，同时将分组列转换为因子变量
chao1 <- soil[ ,c('sample', 'site', 'chao1')]
chao1$site <- factor(chao1$site)
str(chao1)
head(chao1)
#QQ-plot 检查数据是否符合正态分布（所有的点都离直线很近，落在置信区间内说明正态性良好）
library(car)
qqPlot(lm(chao1~site, data = chao1), simulate = TRUE, main = 'QQ Plot', labels = FALSE)
##Shapiro-Wilk 检验，当且仅当两者 p 值均大于 0.05 时表明数据符合正态分布
shapiro <- tapply(chao1$chao1, chao1$site, shapiro.test)
shapiro
#使用 Bartlett 检验进行方差齐性检验（p 值大于 0.05 说明方差齐整）
bartlett.test(chao1~site, data = chao1)
#满足假设，单因素方差分析，详情使用?aov查看帮助
fit <- aov(chao1~site, data = chao1)
summary(fit)
#不满足假设，可使用非参数方法，例如 Kruskal-Wallis Test、Friedman Test
#kruskal.test()、friedman.test()
#若想查看各组均值及标准差，可使用 aggregate()
chao1_mean <- aggregate(chao1$chao1, by = list(chao1$site), FUN = mean)
chao1_sd <- aggregate(chao1$chao1, by = list(chao1$site), FUN = sd)
##方差分析后，多重比较，继续探寻两两分组间的差异
#Tukey HSD 检验
tuk <- TukeyHSD(fit, conf.level = 0.95)
plot(tuk)
tuk
#或者更直观地展示 Tukey HSD 检验结果
library(multcomp)
tuk <- glht(fit, alternative = 'two.sided', linfct = mcp(site = 'Tukey'))
plot(cld(tuk, level = 0.05, decreasing = TRUE))
#不妨自己作个图展示下（ggplot2 柱状图示例）
dat <- merge(chao1_mean, chao1_sd, by = 'Group.1')
names(dat) <- c('group', 'mean', 'sd')
dat <- cbind(dat, sign = c('a', 'b', 'b'))
library(ggplot2)
ggplot(dat, aes(group, mean)) +
geom_col(aes(fill = group), width = 0.4, show.legend = FALSE) +
geom_errorbar(aes(ymax = mean + sd, ymin = mean - sd), width = 0.15, size = 0.5) +
geom_text(aes(label = sign, y = mean +sd + 200)) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), plot.title = element_text(hjust = 0.5)) +
labs(x = 'Group', y = 'Chao1', title = 'Tukey HSD test')
网友评论