R语言执行双因素方差分析

双因素方差分析（双因素ANOVA）在R中实现

前文简介了单因素方差分析（单因素ANOVA，一组因子变量对应一组因变量）、单因素协方差分析（单因素ANCOVA，在单因素ANOVA的基础上，加入协变量，即一组因子变量和协变量对应一组因变量）在R语言中的实现方法，本篇继续简介双因素方差分析（双因素ANOVA）。双因素方差分析，顾名思义，两组因子变量对应一组因变量。

本文使用的作图数据的网盘链接（提取码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', 'treat', 'times', 'chao1')]
chao1$treat <- factor(chao1$treat)
chao1$times <- factor(chao1$times)
str(chao1)
head(chao1)

假设存在这么一个关注的问题：

我们采集了来源于同一环境中的土壤（土壤类型一致），分为了3组，分别添加了三种类型的化学物质（a、b、c），并将土壤孵育了不同的时间（10、15、20、25 天）。在不同时间段收集样本后，通过16S测序，获得了土壤细菌群落的Alpha多样性指数，期望关注化学物质类型以及处理时间对土壤细菌群落的影响（关注交互作用）。对应于上述挑选出的测试数据“chao1”：sample，试验样本名称，土壤类型完全相同；treat，在土壤中添加的三种化学物质（a、b、c），这列作为分组列，需要转换为因子变量类型，各组处理间相互独立；times，土壤孵育时间（天数），这列作为分组列，需要转换为因子变量类型（尽管它本来是数值类型的，但必需转化为因子类型后，才可作为分组变量用于方差分析）；chao1，Alpha多样性指数中的Chao1指数，数值变量。此处存在“化学物质类型”以及“土壤孵育时间”两组分组变量，对应于双因素，同时还需要关注二者的交互作用，接下来我们考虑使用双因素方差分析（双因素ANOVA）来探究化学物质类型以及处理时间是否对土壤细菌群落产生了显著影响。

评估检验的假设条件

双因素ANOVA建立在两组数据必须均服从正态分布，以及各组方差相等的基础上才能执行。

正态性检验

同单因素ANOVA的方法，可使用Q-Q图来检验正态性假设。这里需要检查两组（化学物质类型、处理时间）是否均满足。#QQ-plot 检查数据是否符合正态分布
par(mfrow = c(1, 2))
qqPlot(lm(chao1~treat, data = chao1), simulate = TRUE, main = 'QQ Plot', labels = FALSE)
qqPlot(lm(chao1~times, data = chao1), simulate = TRUE, main = 'QQ Plot', labels = FALSE)

使用car包中的qqPlot()来绘制Q-Q图，由图可知，我们的数据符合正态分布模型。

方差齐性检验
同单因素ANOVA的方法，可使用Bartlett检验进行方差齐性检验。这里需要检查两组（化学物质类型、处理时间）是否均满足。#使用 Bartlett 检验进行方差齐性检验（p 值大于 0.05 说明方差齐整）
bartlett.test(chao1~treat, data = chao1)
bartlett.test(chao1~times, data = chao1)结果显示，我们的数据各组方差相等。

双因素方差分析（双因素ANOVA）

双因素ANOVA

我们的数据通过了正态性检验和方差齐性检验，接下来进行双因素ANOVA。

如果不满足上述前提假设，一是可以考虑转化数据（当然，我们需要确保转换后的数据能够被合理解释，否则将无意义）；二是可以考虑使用非参数的检验方法，对于双因素方差分析的非参数替代方法，常使用Scheirer-Ray-Hare检验（rcompanion包scheirerRayHare()）。双因素ANOVA的aov()函数书写为aov(y~A*B)的样式，表示考虑所有可能的交互项：A、B以及A和B的交互（A:B），其中A、B分别为两组因子变量。因此，双因素ANOVA的aov()函数也可书写为aov(y~A+B+A:B)的样式。此外，表达式中各效应存在主次顺序，即y~A*B与y~B*A的处理方式是不同的，一般情况下，越基础性的变量越应放在表达式前面。在这里，主因素（treat）在前，次因素（times）在后。#满足假设，双因素方差分析，详情使用 ?aov 查看帮助
fit <- aov(chao1~treat*times, data = chao1)
summary(fit)

#上式等同于
fit <- aov(chao1~treat+times+treat:times, data = chao1)
summary(fit)

#查看各组均值及标准差
aggregate(chao1$chao1, by = list(chao1$treat, chao1$times), FUN = mean)
aggregate(chao1$chao1, by = list(chao1$treat, chao1$times), FUN = sd)

双因素ANOVA结果表明，化学物质类型以及处理时间均对土壤细菌群落产生了显著影响，并且二者交互作用也显著。

交互效应展示

若想展示双因素ANOVA的交互效应，以查看数据分布，有多种方法可供选择。以下提供几种参考方法。#例如，interaction.plot() 函数，展示各组均值趋势
interaction.plot(chao1$times, chao1$treat, chao1$chao1)#再例如，boxplot() 函数，以箱线图展示各组数据分布
boxplot(chao1~treat*times, data = chao1, col = c('#B3DE69', '#FDB462', '#80B1D3'))
#再例如，gplots 包 plotmeans() 函数，展示了均值和误差棒（95% 置信区间，以及各组样本量大小
library(gplots)
plotmeans(chao1~interaction(treat, times), data = chao1, connect = list(c(1, 4, 7, 10), c(2, 5, 8, 11), c(3, 6, 9, 12)))#再例如，HH 包 interaction2wt() 函数
library(HH)
interaction2wt(chao1~treat*times, data = chao1)
友情链接