2021-06-19 R语言执行单因素方差分析（单因素ANOVA）及多重比较

对于两组数据间的差异分析，最常见的方法就是使用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')