R语言执行单因素方差分析及多重比较
对于两组数据间的差异分析，最常见的方法就是使用T检验比较两组均值是否存在显著不同。当拓展到多组（三组及以上）时，使用T检验逐一两两比较的方法无疑是低效的，不仅仅由于需要的检验次数增多，而且发生I型错误（拒绝真）的概率也会增大。Fisher提出一种广义T检验的方法来比较三组及以上总体的均值，称为方差分析（ANOVA）。
说到ANOVA，相信大家也并不陌生，这也是在统计学中最常见的统计推断方法之一。几种常见的ANOVA包含单因素方差分析（单因素ANOVA）、单因素协方差分析（ANCOVA）、双因素方差分析（双因素ANOVA）、重复测量方差分析（重复测量ANOVA）、多元方差分析（MANOVA）等。本篇首先介绍其中最常涉及的单因素ANOVA在R语言中的实现过程，一组因子变量对应一组因变量；其它几种类型，会在后续的几篇文章中再一一阐述。
本文使用的作图数据的网盘链接（提取码z4w4）：
https://pan.baidu.com/s/1J-9GsmoHuQ_CEpxeWyEQsA数据预处理
示例数据说明
我们首先将示例数据读到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')
#以 chao1 指数为例，同时将分组列转换为因子变量
chao1 <- soil[ ,c('sample', 'site', 'chao1')]
chao1$site <- factor(chao1$site)
str(chao1)
head(chao1)
假设存在这么一个研究：
我们在3个地域（A、B、C）分别采集了土壤样本，即获得了3种类型的土壤，并通过16S测序，获得了每种类型土壤中细菌群落的Alpha多样性指数。我们想要得知，三种土壤环境下的细菌群落的Alpha多样性指数是否存在显著不同。
对应于上述挑选出的测试数据“chao1”：sample，采集的土壤样本名称；site，土壤样本来源的环境（A、B、C），这列作为分组列，需要转换为因子变量类型，各组之间相互独立；chao1，Alpha多样性指数中的Chao1指数，数值变量。
由于此处只存在“土壤采集环境”这么一组分组变量，对应于单因素，接下来我们考虑使用单因素ANOVA来对3种土壤环境下的细菌群落的Chao1指数进行比较。
评估检验的假设条件
与T检验相似，ANOVA同样要求数据服从正态分布；此外，ANOVA还建立在各组方差相等的基础上。因此，在执行单因素ANOVA之前，我们首先应当对数据进行正态性分布验证，以及方差齐性检验。
正态性检验
首先是正态性检验，这里使用Q-Q图来检验正态性假设。除了Q-Q图，其它的常用方法还有如Shapiro-Wilk检验等。#QQ-plot 检查数据是否符合正态分布（所有的点都离直线很近，落在置信区间内说明正态性良好）library(car)
qqPlot(lm(chao1~site, data = chao1), simulate = TRUE, main = 'QQ Plot', labels = FALSE)qqPlot()提供了精确的正态假设检验方法，它画出了在n-p-1个自由度的t分布下的学生化残差（studentized residual，也称学生化删除残差或折叠化残差）图形，其中n是样本大小，p是回归参数的数目（包括截距项）。图中横坐标是标准的正态分布值，纵坐标是我们数据的值。如果两者基本相等，或者说所有的点都离直线很近，落在置信区间内（图中虚线部分，默认展示95%置信区间），即表明正态性假设符合得很好。由图可知，我们的数据符合正态分布模型。
方差齐性检验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)
结果显示，我们的数据各组方差相等。（尽管示例数据经有偏离的趋势了，凑合用吧）
单因素方差分析（单因素ANOVA）
单因素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)
#若想查看各组均值及标准差，可使用 aggregate()
chao1_mean <- aggregate(chao1$chao1, by = list(chao1$site), FUN = mean)
chao1_sd <- aggregate(chao1$chao1, by = list(chao1$site), FUN = sd)
单因素ANOVA结果表明，3种土壤环境下的细菌群落的Chao1指数具有显著差异，p值远低于0.05水平。
多重比较
#Tukey HSD 检验
tuk <- TukeyHSD(fit, conf.level = 0.95)
plot(tuk)
显著水平默认为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二者无差异。
这里顺便再提一个可能存在的误区。上述标注显著性abc时，由最大值从a开始，逐渐往小值标注b、c等；而有些图中，由最小值从a开始，逐渐往大值标注b、c等（上述cld()参数中，你使用decreasing = FALSE就反过来了）。事实上这两种方法都是可以的，只是更普遍的习惯可能是由大值逐渐往小值标注，所以可能好多同学误以为反过来是错误的。
ggplot2柱状图示例
通过上述各步，我们初步获得了各组间差异分析结果。在文献中，常能见到以均值±误差棒（常用标准差或标准误差）的柱状图，对ANOVA的结果可视化呈现，组间差异水平高低一目了然。
这里根据上述统计结果，简单地使用ggplot2绘制柱状图，以展示3种土壤环境下的细菌群落的Chao1指数的差异水平。
#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')
友情链接