R语言执行多元方差分析
前文已分别简介了单因素方差分析（单因素ANOVA）、单因素协方差分析（单因素ANCOVA）、双因素方差分析（双因素ANOVA）、重复测量方差分析（重复测量ANOVA）在R语言中的实现方法。通过前文几种常见的方差分析示例，我们可知当因子变量只有一组时，称为单因素方差分析，因子变量有两组时，称为双因素方差分析，当因子变量存在多组时，即为多因素方差分析（因子变量越多，解释起来也就越复杂，所以一般不会涉及更多因素）；存在协变量时，称为协方差分析。
数据预处理
示例数据说明
我们首先将示例数据读到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')
#选择数据，并将分组列转换为因子变量
muti <- soil[ ,c('sample', 'treat', 'chao1', 'pH', 'NR')]
muti$treat <- factor(muti$treat)
str(muti)
head(muti)
假设存在这么一个关注的问题：
我们采集了来源于同一环境中的土壤（土壤类型一致），分为了3组，分别添加了三种类型的化学物质（a、b、c），并将土壤孵育了一定的时间（时间相同）。最后取样后，通过16S测序，获得了土壤细菌群落的Alpha多样性指数；通过土壤理化测定，获得了土壤中多种理化指标；通过土壤酶活性测定，获得了主要的几种土壤酶活性数据。通过这些数据，我们想要得知土壤细菌群落的Alpha指数、土壤理化、以及土壤酶活性是否因所添加化学物质类型不同而显著改变。对应于上述挑选出的测试数据“muti”：sample，试验样本名称，土壤类型完全相同，土壤孵育时间完全一致；treat，在土壤中添加的三种化学物质（a、b、c），这列作为分组列，需要转换为因子变量类型，各组处理间相互独立；chao1，Alpha多样性指数中的Chao1指数，数值变量；pH，土壤理化数据中的pH值数据，数值变量；NR，土壤硝酸还原酶活性，数值变量。这里存在3组因变量：“土壤菌群Chao1指数”、“土壤pH值”、“土壤硝酸还原酶（NR）活性”，对应于1组因子变量“化学物质类型”，因此我们使用单因素多元方差分析（单因素MANOVA），探究土壤细菌群落Chao1指数、土壤pH、土壤硝酸还原酶（NR）活性是否因所添加化学物质类型不同而发生显著改变。评估检验的假设条件
单因素MANOVA有两个前提假设，一是多元正态性，二是方差-协方差矩阵同质性。
多元正态性验证
所谓多元正态性，即指因变量组合成的向量服从一个多元正态分布，在R中同样可使用Q-Q图验证多元正态性。#QQ-plot 检验多元正态性y <- cbind(muti$NR, muti$pH, muti$chao1)
coord <- qqplot(qchisq(ppoints(nrow(y)), df = ncol(y)), mahalanobis(y, colMeans(y), cov(y)))
abline(a = 0, b = 1)
若数据服从多元正态性，则点将落在直线上。根据结果可知，我们的数据基本服从多元正态性。
不过我们看到似乎有两个离群点。
这时可以继续使用identify()，交互式地在图中点击这两个点，查看它们是那些样本（本示例中是C4_2和C4_3两个样本）。若有必要，可以将这两个样本剔除，然后再继续下一步（本示例演示不再将它们删除，直接进入下一步分析）。
#可以交互式展示样本位置，可用于观测离群点identify(coord$x, coord$y, labels = muti$sample)
方差-协方差矩阵同质性即指各组的协方差矩阵相同，通常可使用Box’s M检验来评估该假设。注：Box’s M检验对正态性假设很敏感。
这里使用biotools包中的boxM()函数来实现，详情可使用?boxM参阅帮助文档。#Box's M 检验验证方差-协方差矩阵同质性（p 值大于 0.05 即说明各组的协方差矩阵相同）
library(biotools)
boxM(muti[ ,c('chao1', 'pH', 'NR')], muti[ ,'treat'])很遗憾，我们的数据并未通过前提假设。
单因素多元方差分析（单因素MANOVA）
单因素MANOVA
一般来讲，相较于一元的方差分析，MANOVA的前提假设非常地苛刻，大多数情况下都是直接拒绝的。我们的数据就未通过前提假设，理论上不能再执行单因素MANOVA了。但是请允许我继续使用单因素MANOVA来做，仅仅用来作个方法的演示。
前提假设未满足的前提下，可以尝试使用稳健多元方差分析（稳健MANOVA，如rrcov包Wilks.test()，参见下文）；或者更换非参数MANOVA，如置换多元方差分析（PERMANOVA）（vegan包adonis()）。
#假设通过，执行单因素 MANOVA，详情使用 ?manova 查看帮助。fit <- manova(cbind(NR, pH, chao1)~treat, data = muti)
summary(fit) #查看整体结果
summary.aov(fit) #对每一个变量做单因素方差分析
结果中，首先整体显著，其次对各因变量的结果也显著，即土壤细菌群落Chao1指数、土壤pH、土壤硝酸还原酶（NR）活性均因所添加化学物质类型不同而发生显著改变。
对于3个子单因素方差分析（单因素ANOVA）结果，我们还可继续使用多重比较（Tukey HSD检验），探究对于每个因变量，3种化学物质的处理结果之间更具体的差异是怎样的。可参考前文单因素ANOVA。
稳健单因素MANOVA
#通过 rrcov 包中的 Wilks.test() 函数实现，详情可使用 ?Wilks.test 查看帮助
library(rrcov)
muti$treat <- factor(muti$treat)
Wilks.test(treat~., data = muti[c('treat', 'NR', 'pH', 'chao1')], method = 'c')
结果表明，土壤细菌群落Chao1指数、土壤pH、土壤硝酸还原酶（NR）活性均因所添加化学物质类型不同而发生显著改变。