对于两组数据间的差异分析，最常见的方法就是使用T检验比较两组均值是否存在显著不同。当拓展到多组（三组及以上）时，使用T检验逐一两两比较的方法无疑是低效的，不仅仅由于需要的检验次数增多，而且发生I型错误的概率也会增大。Fisher提出一种广义T检验的方法来比较三组及以上总体的均值，称为方差分析（ANOVA）。
说到ANOVA，相信大家也并不陌生，这也是在统计学中最常见的统计推断方法之一。本文就几种常见的ANOVA方法为例，包含单因素方差分析（单因素ANOVA）、单因素协方差分析（ANCOVA）、双因素方差分析（双因素ANOVA）、重复测量方差分析（重复测量ANOVA）、多元方差分析（MANOVA），简介怎样在R中进行ANOVA，以实现多组间数据总体均值的差异分析。
示例数据、R脚本等，已上传至百度盘（提取码z4w4）：
/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')
单因素方差分析（单因素ANOVA）
单因素方差分析是肯定是最常见的方差分析了，一组因子变量对应一组因变量。
示例数据说明
从总示例数据中挑选部分数据作为测试。
#以chao1指数为例，同时将分组列转换为因子变量chao1<-soil[,c('sample','site','chao1')]chao1$site<-factor(chao1$site)
假设存在这么一个研究：
我们在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检验，我在介绍T检验的博文中有提及，本文就不再讲述了。
library(car)#使用其中的qqPlot()用于绘制QQ图#QQ-plot检查数据是否符合正态分布qqPlot(lm(chao1~site,data=chao1),simulate=TRUE,main='QQPlot',labels=FALSE)使用car包中的qqPlot()来绘制Q-Q图。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。R语言执行方差分析的命令是aov()，详情使用?aov查看帮助，本文末的附录中也有简介。对于单因素方差分析，aov()函数书写为aov(y~A)的样式，A即为因子变量。
备注：如果不满足上述前提假设，可以尝试转换数据，例如，log转换等或许可以使非正态分布的原始数据转变为正态分布类型（当然，我们需要确保转换后的数据能够被合理解释，否则将无意义）；或者可以考虑使用非参数的检验方法，对于单因素的分析，可选的非参数替代方法例如Kruskal-Wallis检验（kruskal.test()）、Friedman检验（friedman.test()）等。
#满足假设，单因素方差分析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水平。
多重比较
上述单因素ANOVA告诉我们3种土壤环境下的细菌群落的Chao1指数具有显著差异，这种差异是在整体水平而言的，并没有告诉我们究竟谁和谁存在差异。如果我们想继续获知两两分组之间的差异，进行多重比较即可。
常用Tukey HSD检验，在ANOVA结果的基础上继续执行事后两两比较。不推荐使用T检验（T检验和Tukey检验是两回事），原因正如本文开始时所提，多次T检验容易提高I型错误的概率。
##方差分析后，多重比较，继续探寻两两分组间的差异#TukeyHSD检验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二者无差异。
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='TukeyHSDtest')示例结果如下。
除了柱状图，好的可视化方案还有很多。由于本文的目的并非可视化，因此不再举更多示例，还需大家在文献中获取灵感了。
单因素协方差分析（ANCOVA）
当方差分析中存在协变量时，即可称为协方差分析。其中单因素协方差分析是最常见的，在单因素方差分析中引入了协变量。
示例数据说明
从总示例数据中挑选部分数据作为测试。
#以shannon指数为例，同时将分组列转换为因子变量shannon<-soil[,c('sample','treat','shannon','days')]shannon$treat<-factor(shannon$treat)
假设存在这么一个研究：
我们使用来源于同一环境中的土壤进行盆栽试验（土壤类型一致），并种植了同种植株（植物类型一致）。我们将盆栽（1植株/1盆栽）分为了3组，分别在土中添加了三种化学物质（a、b、c）；然后等待植物到达花期后，收集每个植株的根际土，并通过16S测序，获得了植物根际细菌群落的Alpha多样性指数，意在探究不同类型的化学物质是如何影响植物根际细菌群落的。但由于各盆栽中植株开花期时间并非一致，不同个体间在开花时间上可能相隔数天（这个我们无法控制），而该植物开花时间又很短，可能早开花的植株还未等晚开花的植株开花就已经凋谢，因此无法保证所有植株均在同一天取样观察，于是我们就先开花的先取样后开花的后取样。尽管期间并未相隔很长时间（天），理论上单因素ANOVA就可以满足需求，但我们仍然想要将植物生长时间这个因素考虑在内（如果觉得植物生长时间相隔几天的差异可能会导致其根际菌群产生较大的变异时），即将它作为协变量处理（协变量是独立变量，实验者不能操纵，但仍影响实验结果），尝试使用单因素协方差分析（ANCOVA）。
对应于上述挑选出的测试数据“shannon”：sample，试验样本名称，每一个样本即对应一个盆栽，各盆栽中土壤类型、植物类型完全相同，而且均是1植株/1盆栽；treat，在土壤中添加的三种化学物质（a、b、c），这列作为分组列，需要转换为因子变量类型，各组处理间相互独立；shannon，Alpha多样性指数中的Shannon指数，数值变量；days，各盆栽中植株的开花时间（即生长天数），这里为数值变量。
这里，植物根际菌群的Shannon指数为因变量，植物生长天数（days）为协变量，考虑使用单因素协方差分析（ANCOVA），探究不同类型的化学物质处理下的植物根际细菌群落的Shannon指数是否存在显著不同。
评估检验的假设条件
同样地，ANCOVA要求数据服从正态分布，以及各组方差相等。此外，ANCOVA还假定回归斜率相同。
正态性检验
同上文的方法，使用Q-Q图来检验正态性假设。
#QQ-plot检查数据是否符合正态分布qqPlot(lm(shannon~treat,data=shannon),simulate=TRUE,main='QQPlot',labels=FALSE)使用car包中的qqPlot()来绘制Q-Q图，由图可知，我们的数据符合正态分布模型（尽管似乎有个离群点，这时候可以考虑删除这个样本后再继续，本示例演示不再将它删除，直接进入下一步分析）。
方差齐性检验
同上文，使用Bartlett检验进行方差齐性检验。
#使用Bartlett检验进行方差齐性检验（p值大于0.05说明方差齐整）bartlett.test(shannon~treat,data=shannon)结果显示，我们的数据各组方差相等。
回归斜率的同质性检验
ANCOVA包含植物生长时间×化学物质类型的交互项时，可对回归斜率的同质性进行检验，若交互效应显著，则意味着植物生长时间和植物根际菌群的Shannon指数的关系依赖于所添加的化学物质类型。
#检验回归斜率的同质性，交互效应不显著则支持斜率相等的假设（即p值大于0.05说明回归斜率相等）summary(aov(shannon~days*treat,data=shannon))根据aov()公式可知，这实际上是一个双因素方差分析（见下文），根据双因素方差分析中的交互项结果来判断回归斜率的同质性。
结果显示交互作用不显著，支持了斜率相等的假设。如果假设不成立，可以尝试变换协变量或因变量，或使用能对每个斜率独立解释的模型，或使用无需回归斜率相等的非参数ANCOVA方法（如sm包sm.ancova()）。
单因素协方差分析（ANCOVA）
上述检验均已通过，接下来进行ANCOVA。
备注：同样地，如果不满足上述前提假设，一是可以考虑转化数据（当然，我们需要确保转换后的数据能够被合理解释，否则将无意义），二是可以考虑使用非参数的检验方法。上述提及了一个对应单因素协方差分析的非参数替代方法，sm包sm.ancova()。
#满足假设，单因素协方差分析fit<-aov(shannon~days+treat,data=shannon)summary(fit)同样地，使用aov()执行，更改公式书写方式即可。对于带协变量的项，以单因素协方差为例，aov()函数书写为aov(y~x+A)的样式，其中x为协变量，A为因子变量，注意需要将协变量写在因子前面，如上所示，协变量植物生长时间（days）在前，因子化学物质类型（treat）在后，顺序不可颠倒。
ANCOVA结果表明：（1）植物生长时间相隔几天的差异并未导致其根际菌群产生较大的变异（p值不显著）；（2）控制植物生长时间，不同类型的化学物质处理下的植物根际细菌群落的Shannon指数存在显著不同（p值显著）。
在这种协变量并未对因变量产生显著影响的例子中，直接使用单因素ANOVA（忽略协变量的存在）其实就可以了。
对于各组均值的获得方式。
#查看各组均值及标准差aggregate(shannon$shannon,by=list(shannon$treat),FUN=mean)aggregate(shannon$shannon,by=list(shannon$treat),FUN=sd)#由于使用了协变量，若想获取去除协变量效应后的组均值（调整的组均值）library(effects)effect('treat',fit)
因变量、协变量和因子之间的关系图
我们可以使用HH包ancova()函数，绘制因变量、协变量和因子之间的关系图查看。详情使用?ancova查看帮助。
#HH包ancova()可绘制因变量、协变量和因子之间的关系图library(HH)ancova(shannon~days+treat,data=shannon)ancova(shannon~days*treat,data=shannon)ancova()函数既可根据输入的公式执行对应的方差分析，并在屏幕输出方差分析结果；同时又可生成一张关系图，由两部分组成，左侧面板取决于因子的水平groups。右侧面板是所有groups的叠加。
如上示例，通过“ancova(shannon~days+treat, data = shannon)”，我们又执行了一次ANCOVA，结果屏幕输出，和上文结果一致。同时通过关系图可知，3条回归线相互平行，只是截距项不同，b组截距项最大，c组截距项最小；回归线拟合效果并不理想，也对应了我们先前的结论，在“数天”这么一个短期时间内，植物根际菌群的Shannon指数没有发生剧烈的改变。
我们还通过“ancova(shannon~days*treat, data = shannon)”，执行了一次双因素ANOVA。不过在这里我们并不是想做个双因素ANOVA分析，而是在更改了函数公式后，意在可视化允许斜率和截距项依据组别而发生改变，从而体现那些违背回归斜率同质性的实例。（上文已知，回归斜率的同质性检验是通过双因素ANOVA中的交互作用判断的，本示例中的回归斜率的同质性检验已通过，大家可以尝试自行寻找一例回归斜率不相等的数据，然后使用ancova()作图查看其与回归斜率相等的数据的差异）
双因素方差分析（双因素ANOVA）
双因素方差分析，顾名思义，两组因子变量对应一组因变量。
示例数据说明
从总示例数据中挑选部分数据作为测试。
#同样以chao1指数为例，同时将分组列转换为因子变量chao1<-soil[,c('sample','treat','times','chao1')]chao1$treat<-factor(chao1$treat)chao1$times<-factor(chao1$times)
假设存在这么一个研究：
我们采集了来源于同一环境中的土壤（土壤类型一致），分为了3组，分别添加了三种类型的化学物质（a、b、c），并将土壤孵育了不同的时间（10、15、20、25 天）。在不同时间段收集样本后，通过16S测序，获得了土壤细菌群落的Alpha多样性指数，期望关注化学物质类型以及处理时间对土壤细菌群落的影响（关注交互作用）。
对应于上述挑选出的测试数据“chao1”：sample，试验样本名称，土壤类型完全相同；treat，在土壤中添加的三种化学物质（a、b、c），这列作为分组列，需要转换为因子变量类型，各组处理间相互独立；times，土壤孵育时间（天数），这列作为分组列，需要转换为因子变量类型（尽管它本来是数值类型的，但必需转化为因子类型后，才可作为分组变量用于方差分析）；chao1，Alpha多样性指数中的Chao1指数，数值变量。
此处存在“化学物质类型”以及“土壤孵育时间”两组分组变量，对应于双因素，同时还需要关注二者的交互作用，接下来我们考虑使用双因素方差分析（双因素ANOVA）来探究化学物质类型以及处理时间是否对土壤细菌群落产生了显著影响。
评估检验的假设条件
同样地，双因素ANOVA要求数据服从正态分布，以及各组方差相等。
正态性检验
同上文的方法，使用Q-Q图来检验正态性假设。
#QQ-plot检查数据是否符合正态分布par(mfrow=c(1,2))qqPlot(lm(chao1~treat,data=chao1),simulate=TRUE,main='QQPlot',labels=FALSE)qqPlot(lm(chao1~times,data=chao1),simulate=TRUE,main='QQPlot',labels=FALSE)这里需要检查两组是否均满足。使用car包中的qqPlot()来绘制Q-Q图，由图可知，我们的数据符合正态分布模型。
方差齐性检验
同上文，使用Bartlett检验进行方差齐性检验。
#使用Bartlett检验进行方差齐性检验（p值大于0.05说明方差齐整）bartlett.test(chao1~treat,data=chao1)bartlett.test(chao1~times,data=chao1)这里需要检查两组是否均满足。结果显示，我们的数据各组方差相等。
双因素方差分析（双因素ANOVA）
我们的数据通过了正态性检验和方差齐性检验，接下来进行双因素ANOVA。
备注：同样地，如果不满足上述前提假设，一是可以考虑转化数据（当然，我们需要确保转换后的数据能够被合理解释，否则将无意义）；二是可以考虑使用非参数的检验方法，对于双因素方差分析的非参数替代方法，常使用Scheirer-Ray-Hare检验（rcompanion包scheirerRayHare()）。
#满足假设，双因素方差分析fit<-aov(chao1~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)同样地，使用aov()执行，更改公式书写方式即可。双因素ANOVA的aov()函数书写为aov(y~A*B)的样式，表示考虑所有可能的交互项：A、B以及A和B的交互（A:B），其中A、B分别为两组因子变量。
双因素ANOVA结果表明，化学物质类型以及处理时间均对土壤细菌群落产生了显著影响，并且二者交互作用也显著。
交互效应展示示例
若想展示双因素ANOVA的交互效应，以查看数据分布，有多种方法可供选择。以下提供几种参考方法。
#例如，interaction.plot()函数，展示各组均值趋势interaction.plot(chao1$times,chao1$treat,chao1$chao1)
#再例如，boxplot()函数，以箱线图展示各组数据分布boxplot(chao1~treat*times,data=chao1,col=c('red','green','blue'))
#再例如，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)
重复测量方差分析（重复测量ANOVA）
所谓重复测量方差分析（重复测量ANOVA），即受试者被测量不止一次，含一个组内和一个组间因子的重复测量方差分析（这也是一个常见的设计）。
示例数据
示例数据就直接使用R中自带的CO2数据集，以下直接搬运《R语言实战（第二版）》214页中的原话。
CO2数据集，包含了北方和南方牧草类植物Echinochloa crus-galli的寒冷容忍度研究结果，在某浓度二氧化碳的环境中，对寒带植物与非寒带植物的光合作用率进行了比较。研究所用植物一半来自于加拿大的魁北克省（Quebec），另一半来自美国的密西西比州（Mississippi）。
我们只关注其中的寒带植物。因变量是二氧化碳吸收量（uptake），单位为ml/L，自变量是植物类型Type（魁北克VS.密西西比）和七种水平（95~1000 umol/m^2 sec）的二氧化碳浓度（conc）。Type是组间因子，conc是组内因子。
data(CO2)w1b1<-subset(CO2,Treatment=='chilled')#我们只关注其中的寒带植物#将分组转变为因子类型w1b1$Type<-factor(w1b1$Type)w1b1$conc<-factor(w1b1$conc)
评估检验的假设条件
重复测量ANOVA同样要求数据服从正态分布，以及各组方差相等。这里我就不再演示这一步了，暂且默认数据满足前提假设，直接进入方差分析过程。
重复测量方差分析
同样地，使用aov()执行，更改公式书写方式即可。含单个组内因子（W）和单个组间因子（B）的重复测量方差分析的aov()函数书写为aov(y~B*W+Error(Subject/W))的样式。
#重复测量方差分析fit<-aov(uptake~conc*Type+Error(Plant/conc),data=w1b1)summary(fit)结果显示，主效应类型和浓度以及交叉效应类型×浓度都非常显著。
交互效应展示示例
对于重复测量ANOVA中的交互效应，同样可以使用interaction.plot()、boxplot()等命令来查看，如下示例使用boxplot()，更多可参考上文双因素ANOVA中的交互效应展示示例。
#交互效应展示示例boxplot(uptake~Type*conc,data=w1b1,col=c('red','blue'),las=2)结果显示，Quebec的寒带植物比Mississippi的寒带植物的CO2吸收率高，且随着CO2浓度升高，差异越来越明显。
多元方差分析（MANOVA）
通过上文几种常见的方差分析示例，我们可知当因子变量只有一组时，称为单因素方差分析，因子变量有两组时，称为双因素方差分析，当因子变量存在多组时，即为多因素方差分析（因子变量越多，解释起来也就越复杂，所以一般不会涉及更多因素）；存在协变量时，称为协方差分析。
上文所示的方差分析，因变量只有一种，即一个因变量对应一个或多个因子变量或协变量。当因变量不止一个时，即一个或多个因子变量对应了多个因变量时，可使用多元方差分析（MANOVA）。同样地，当因子变量只有一组时，称为单因素多元方差分析，因子变量有多组时，称为多因素多元方差分析。下文将展示一例单因素多元方差分析的示例。
示例数据说明
从总示例数据中挑选部分数据作为测试。
#选择数据，并将分组列转换为因子变量muti<-soil[,c('sample','treat','chao1','pH','NR')]muti$treat<-factor(muti$treat)
假设存在这么一个研究：
我们采集了来源于同一环境中的土壤（土壤类型一致），分为了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'sM检验验证方差-协方差矩阵同质性（p值大于0.05即说明各组的协方差矩阵相同）library(biotools)boxM(muti[,c('chao1','pH','NR')],muti[,'treat'])很遗憾，我们的数据并未通过前提假设。
单因素多元方差分析（单因素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）
如果多元正态性或方差-协方差矩阵同质性的前提假设未能满足（我们的示例数据即是如此），或者比较担心多元离群点的影响，那么可以尝试稳健多元方差分析（稳健MANOVA）。
对于我们的示例数据，即对应了稳健单因素多元方差分析（稳健单因素MANOVA）。在R中，可使用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）活性均因所添加化学物质类型不同而发生显著改变。
回归和ANOVA
以上均使用方差分析aov()函数来实现方差分析过程。
事实上，从函数形式上看，ANOVA是广义线性模型的特例，尽管ANOVA和回归方法都是独立发展而来的。因此在R中，ANOVA也可使用回归函数lm()来完成，lm()所得结果将和方差分析函数aov()的结果一致。本文不详细阐述ANOVA和回归的关系，将来写到回归分析时再叙述。但是为了加深理解，以下将简单展示使用回归分析函数lm()实现一个单因素方差分析的例子。
##使用回归函数lm()实现方差分析#以单因素方差分析为例，数据同上文#假设样本采自三种土壤环境（A、B、C），我们比较三种土壤环境下的细菌群落的chao1指数是否存在显著差异#以chao1指数为例，同时将分组列转换为因子变量chao1<-soil[,c('sample','site','chao1')]chao1$site<-factor(chao1$site)#正态性、方差齐性检验，略#这个上文已经通过了检验#aov()实现单因素方差分析fit_aov<-aov(chao1~site,data=chao1)summary(fit_aov)#lm()实现单因素方差分析fit_lm<-lm(chao1~site,data=chao1)summary(fit_lm)结果如下，比较aov()和lm()的结果，发现二者的所得p值是一致的。
这里有一点需要注意。线性模型要求输入的预测变量（自变量）为数值类型，而在ANOVA中，分组变量（自变量）为因子类型。因此，当我们使用lm()执行ANOVA时，由于输入数据中的分组变量为因子类型，lm()会首先使用一系列与因子水平相对应的数值型对照变量来代替因子，然后再执行线性拟合。具体细节本文不再细说，待将来写到线性模型时再叙述。
附录，关于R语言方差分析函数aov()
以下参考自《R语言实战（第二版）》201-202页中的内容。
aov()表达式
R语言运行ANOVA的默认函数是aov()，aov()函数的基本语法为aov(formula, data=dataframe)。下表列举了一些常见的ANOVA表达式，各式中“~”左侧y为因变量，右侧大写字母表示组别因子（因子变量），小写字母表示定量变量（主要用作协变量），Subject是对被测试者独有的标识变量。
表达式中各项的顺序
表达式中效应的顺序在两种情况下会造成影响：（a）因子不止一个，并且是非平衡设计；（b）存在协变量。出现任意一种情况时，等式右边的变量都与其他每个变量相关。此时，我们无法清晰地划分它们对因变量的影响。例如，对于双因素方差分析，若不同处理方式中的观测数不同，那么y~A*B与y~B*A的结果将不同。
R默认类型I（序贯型，见下图）方法计算ANOVA效应。若某模型写作y~A*B（根据上表，可知该式等同于y ~ A + B + A:B），那么在R中ANOVA表的结果将评价：
R中的ANOVA表的结果将评价：
（1）A对y的影响；
（2）控制A时，B对y的影响；
（3）控制A和B的主效应时，A与B的交互效应。
样本大小越不平衡，效应项的顺序对结果的影响越大。一般来说，越基础性的效应越需要放在表达式前面。具体来讲，首先是协变量，然后是主效应，接着是双因素的交互项，再接着是三因素的交互项，以此类推。对于主效应，越基础性的变量越应放在表达式前面，例如对于我们上文展示的双因素ANOVA示例，化学物质类型即需要在前，而处理时间在后。有一个基本是准则：若研究设计不是正交的（因子和/或协变量相关），一定要谨慎设置效应的顺序。
参考资料
Robert I. Kabacoff. R语言实战（第二版）（王小宁 刘撷芯 黄俊文 等 译）. 人民邮电出版社, .
钱松. 环境与生态统计:R语言的应用（曾思育 译）. 高等教育出版社, .
