1、多因素的试验设计，有时需要分析因素间的交互影响（interaction），2个因素间的交互影响称为一级交互影响（AB）；3个因素间的交互影响称为二级交互影响（ABC）。当交互影响项呈现统计不显著时，表明各个因素独立，当呈现统计显著时，就需要列出这个交互影响项的效应，以助于作出正确的统计推断。举例解释上述概念：要考察焦虑症的治疗疗效，一个因素是治疗方案，有2种治疗方案，即该因素有2个水平；（治疗方案称为组间因子，因为每个患者只能被分配到一个组别中，没有患者同时接受两种治疗）；再考虑一个因素治疗时间，也有两个水平：治疗5周和治疗6个月，同一患者在5周和6个月不止一次地被测量（两次），称为重复测量（
2、治疗时间称为组内因子，因为每个患者在所有水平下都进行了测量）。建立方差分析模型时，既要考虑两个因素治疗方案和治疗时间（主效应），又要考虑治疗方案和时间的交互影响（交互效应），此时即两因素混合模型方差分析。当某个因素的各个水平下的因变量的均值呈现统计显著性差异时，必要时可作两两水平间的比较，称为均值间的两两比较。二、R语言实现方差分析对数据的要求：满足正态性（来自同一正态总体）和方差齐性（各组方差相等），在这两个条件下，若各组有差异，则只可能是来自影响因素的不同水平。用aov（）函数进行方差分析，基本格式为：aov（formula, data=NULL, projections=FALSE, q
3、r=TRUE,contrasts=NULL, .）其中，formula为方差分析公式；data为数据框；projection设置是否返回预测结果；qr设置是否返回QR分解结果；contrasts为公式中一些因子的列表。formula公式的表示：（y为因变量，ABC为分组因子）符号用法分隔符号，左边为响应变量，右边为解释变量eg：yA+B+C+分隔解释变量：表示变量的交互项yA+B+A:B*表示所有可能交互项yA*B*C可展开为：yA+B+C+A:B+A:C+B:C+A:B:C表示交互项达到次数y（A+B+C）2展开为：.表示包含除因变量外的所有变量若一个数据框包括变量y,A、B和C，代码y.可
4、展开为yA+B+C常见研究设计的表达式：（小写字母表示定量变量，大写字母表示组别因子，Subject是对被试者独有的标识变量）设计表达式单因素ANOVAyA含单个协变量的单因素ANCOVAyx+A双因素ANOVAyA*B含两个协变量的双因素ANCOVAyx1+x2+A*B随机化区组yB+A, B为区组因子单因素组内ANOVAyA+Error（Subject/A）含单个组内因子（W）和单个组间因子（B）的重复测量ANOVAyB*W+Error（Subject/W）注意：非均衡设计时或存在协变量时，效应项的顺序对结果影响较大，越基础的效应越需要放在表达式前面，首先是协变量、然后是主效应、接着是双因
5、素的交互项，再接着是三因素的交互项。若研究不是正交的，一定要谨慎设置效应的顺序。有三种类型的方法可以分解yA+B+A:B右边各效应对y所解释的方差：类型I（序贯型）效应根据表达式中先出现的效应做调整。A不做调整，B根据A调整，A:B交互项根据A和B调整。类型II（分层型）效应根据同水平或低水平的效应做调整。A根据B调整，B依据A调整，A:B交互项同时根据A和B调整。类型III（边界型）每个效应根据模型其他各效应做相应调整。A根据B和A:B做调整，A:R默认调用类型I方法，其他软件（比如SAS和SPSS）默认调用类型III方法。car包中的Anova（）函数（不要与标准anova（）函数混淆）提
6、供了使用类型II或类型III方法的选项，而aov（）函数使用的是类型I方法。若想使结果与其他软件（如SAS和SPSS）提供的结果保持一致，可以使用Anova（）函数。三、单因素方差分析1个因变量，1个影响因素：总差异Yij = 平均差异 + 因素差异i + 随机差异ij例1 比较4种品牌的胶合板的耐磨性，各抽取5个样品，相同转速磨损相同时间测得磨损深度（mm），比较4个品牌胶合板的耐磨性有无差异部分数据如下（）：setwd（E:/办公资料/R语言/R语言学习系列/codes）load（head（datas） wear brand1 A2 A3 A4 A5 A6 Battach（datas）ta
7、ble（brand） #各组的样本数brandA B C D 5 5 5 5 aggregate（wear,by=list（brand）,mean） #各组均值 x1 A 2 B 3 C 4 D aggregate（wear,by=list（brand）,sd） #各组标准差1 A 0.2 B 0.3 C 0.library（car）qqPlot（lm（wearbrand,data=datas）,simulate=TRUE） #用Q-Q图检验数据的正态性leveneTest（wear（brand）,data=datas） #方差齐性检验Levenes Test for Homogeneity
8、of Variance （center = median） Df F value Pr（F）group 3 16 fitF） brand 3 *Residuals 16 -Signif. codes: 0 * * * . 1 说明：方差齐性检验，原假设H0：方差齐，p值=, 故接受原假设，即方差齐。单因素方差分析结果，brand是因素，Residuals是残差，各列依次为自由度、平方和、均方和、F统计量，p值=, 拒绝原假设，即不同品牌的磨损（均值）有显著差别。library（gplots）plotmeans（wearbrand,xlab=品牌, ylab=磨损） #图形展示带95%置信区间的
9、各组均值通过前面的分析知道，不同品牌的磨损（均值）有显著差别，但并不知道哪个品牌与其它品牌有显著差别。TukeyHSD（）函数提供了对各组均值差异的成对检验。TukeyHSD（fit） Tukey multiple comparisons of means 95% family-wise confidence levelFit: aov（formula = wear brand, data = datas）$brand diff lwr upr p adjB-A -0. C-A -0. D-A C-B -0. D-B D-C 0. 说明：可以看出（H0：无差异），B与A的差异非常不显著，C与A
10、、C与B、D与C的差异非常显著。multcomp包中的glht（）函数提供了更为全面的多重均值比较方法。library（multcomp）tuk |t|） B - A = 0 C - A = 0 *D - A = 0 . C - B = 0 D - B = 0 . D - C = 0 （Adjusted p values reported - single-step method）plot（cld（tuk, level = , col = lightgrey标记相同字母（标记b）的品牌ABD认为是无显著差异，在同一亚组，而品牌C（标记a）与另外三个品牌有显著差异。另外，也可以进行多重t检验，使
11、用函数： g, 其中，x为因变量，g为因子型的分组变量；设置p值的修正方法，由于多次重复t检验会大大增加犯第一类错误的概率，为此要进行p值的修正，使用bonferroni法修正效果较好。bonferroni Pairwise comparisons using t tests with pooled SD data: wear and brand A B C B - - C - D P value adjustment method: bonferroni原假设H0: 无差异，可见A与B无差异，C与ABD有显著差异。最后，方差分析对离群点非常敏感，检验是否有离群点：outlierTest（fi
12、t）No Studentized residuals with Bonferonni p Largest |rstudent|: rstudent unadjusted p-value Bonferonni p9 经检验无离群点。三、两因素方差分析1个因变量，2个影响因素：总差异Yijk = 平均差异 + 因素1差异i + 因素2差异i+ 因素1,2交互作用差异ij + 随机差异ijk例2 研究60只豚鼠的牙齿生长数据，按2种喂食方法：橙汁、维生素C，各喂食方法中抗坏血酸含量都有3个水平：、1mg、2mg，分配为6组，每组各10只，牙齿长度为因变量。做两因素方差分析。attach（ToothG
13、rowth）head（ToothGrowth） len supp dose1 VC 2 VC 3 VC 4 VC 5 VC 6 VC table（supp, dose） #各组样本数相同，即为均衡设计 dosesupp 1 2 OJ 10 10 10 VC 10 10 10aggregate（len, by=list（supp, dose）, mean） #计算各组均值1 OJ 3 OJ 5 OJ aggregate（len, by=list（supp, dose）, sd） #计算各组标准差（lensupp,data=ToothGrowth） #关于因素supp的方差齐性检验 Bartlet
14、t test of homogeneity of variances len by suppBartletts K-squared = , df = 1, p-value = （lendose,data=ToothGrowth） #关于因素dose的方差齐性检验 len by doses K-squared = , df = 2, p-value = fit-aov（lensupp*dose,data=ToothGrowth） #做两因素方差分析，考虑全部效应supp 1 *dose 1 2e-16 *supp:dose 1 * Residuals 56 :可以看出，主效应supp和dose都
15、非常显著（p值都远小于），交互效应也显著（p值=）。若交互作用不显著，可以可以只做去掉交互效应的方差分析。图形化展示两因素方差分析的交互效应：par（mfrow=c（1,2）（dose, supp, len, type=b, col = c（red, blue）, pch = c（16, 18）, main=Interaction between Dose and Supp（supp, dose, len, type=, col=c（Interaction between Supp and Dose有一个图的线有交叉，说明有交互作用。可以看出随着橙汁和维生素C中的抗坏血酸剂量的增加，牙齿长度变
16、长；。对于 mg和1 mg剂量，橙汁比维生素C更能促进牙齿生长；对于2 mg剂量的抗坏血酸，两种喂食方法下牙齿长度增长相同。也可以用HH包中的interaction2wt（）函数（也适合三因素方差分析）来展示更全面的可视化结果：library（HH）interaction2wt（lensupp*dose）三、重复测量方差分析 重复测量方差分析，即受试者被测量不止一次。例3（1个组内1个组间因子的重复测量）在某浓度CO2的环境中，对寒带植物（来自魁北克）和非寒带植物的（来自密西西比）光合作用率进行比较。因变量uptake为CO2吸收量，自变量Type（组间因子）为植物类型，自变量conc（组内因
17、子）为七种水平的CO2浓度。attach（CO2）head（CO2） #注意CO2是长格式的数据 Plant Type Treatment conc uptake1 Qn1 Quebec nonchilled 95 2 Qn1 Quebec nonchilled 175 3 Qn1 Quebec nonchilled 250 4 Qn1 Quebec nonchilled 350 5 Qn1 Quebec nonchilled 500 6 Qn1 Quebec nonchilled 675 w1b1-subset（CO2, Treatment=chilled） #先只考虑寒带植物-aov（up
18、take（conc*Type）+Error（Plant/conc）, data=w1b1）Error: PlantType 1 *Residuals 4 Plant:concconc 1 *conc:Type 1 * WithinResiduals 30 869 在的显著水平下，主效应“类型”（p值=）和“浓度”（p值=）以及交叉效应“类型*浓度”（p值=）都非常显著。attach（w1b1）（conc, Type, uptake, type=）, pch=c（16, 18）, main=Interaction Plot for Plant Type and Concentrationboxp
19、lot（uptakeType*conc, data=w1b1, col=（c（goldgreen）, main=Chilled Quebec and Mississippi Plants, ylab = Carbon dioxide uptake rate （umol/m2 sec）detach（w1b1）可以看出，魁北克省的植物比密西西比州的植物二氧化碳吸收率高，而且随着CO2浓度的升高，差异越来越明显。注1：重复测量设计时，需要有长格式数据才能拟合模型，若是宽数据，需要用reshape2包中的melt（）函数转化为长数据；注2：这里是用的传统的重复测量方差分析，假设任意组内因子的协方差矩阵
20、为球形，并且任意组内因子两水平间的方差之差都相等。实际中，该假设一般不满足，可以尝试：使用lme4包中的lmer（）函数拟合线性混合模型（Bates，2005）； 使用car包中的Anova（）函数调整传统检验统计量以弥补球形假设的不满足（例如Geisser-Greenhouse校正）； 使用nlme包中的gls（）函数拟合给定方差-协方差结构的广义最小二乘模型（UCLA，2009）； 用多元方差分析对重复测量数据进行建模（Hand，1987）。四、多元方差分析1. 因变量不只一个的方差分析，称为多元方差分析。要求数据满足：多元正态性、方差协方差矩阵同质性（即指各组的协方差矩阵相同，通常可用Boxs M检验来评估该假设）。例4 使用MASS包中的UScereal数据集，研究美国谷物中的卡路里、脂肪、糖含量是否会因为储物架位置的不同而发生变化。自变量货架位置shelf有1, 2, 3三个水平。library（MASS）attach（UScereal）yshelf 1 3 61 *Residuals 63 （fit） #输出每个单变量的方差分析结果 Response calories :shelf 1 45313 45313 *Residuals 63 203982 3238 Response fat :shelf 1 *
