参考学习资料：《R语言实战（第二版）》第九章

前面记录了通过量化的预测变量来预测量化的响应变量的回归模型。这部分内容能将名义型或有序型因子作为预测变量进行建模。当包含的因子是解释型变量时，我们关注的重点通常会从预测转向组别差异的分析，这种分析方法就是方差分析（ANOVA）。

术语速成

实验设计和方差分析都有自己相应的语言。在讨论实验设计分析前，先快速回顾一些重要的术语，并通过对一系列复杂度逐步增加的实验设计的学习，引入模型最核心的思想。
以焦虑症治疗为例：现在有两种治疗方案：认知行为疗法（CBT）和眼动脱敏再加工法（EMDR）。我们招募10位焦虑症患者作为志愿者，随机分配一半的人接受为期5周的CBT，另一半人接受EMDR，设计方案如表1所示。在治疗结束时，要求每位患者都填写状态特质焦虑问卷（STAI），也就是一份焦虑度测量的自我评测报告。

表1 单因素组间方差分析

在这个实验中，治疗方案CRT,EMDR是组间因子。就是两组的人员不会有重叠。表中的字母s就代表每位受试者。STAI是因变量，治疗方案是自变量。由于每组观察数相等，这个设计也称平衡设计（balanced design）;若观测数不同，则称作非均衡设计（unbalanced design）。
因为只有一个类别型变量，表1的统计设计又称为单因素方差分析（one-way ANOVA）,或者进一步称为单因素组间方差分析。方差分析主要通过F检验来进行效果评测，若治疗方案的F检验显著，则说明5周后两种疗法的STAI得分均值不同。
假设你只对CBT的效果感兴趣，则需要将10个受试者都放在CBT组，然后在治疗5周和六个月后分别评价疗效，设计方案如表2所示

表2 单因素组内方差分析

此时，时间（time）是2个水平（5周，6个月）的组内因子。因为每位患者所在的水平下都进行了测量，所有这种设计成为单因素组内方差分析；又由于每个受试者都不止一次呗测量，也称作重复测量方差分析。当时间的F检验显著时，说明STAI得分均值在这2个时间点发生了改变。
现在假设你对治疗方案差异和它随时间的改变都感兴趣，则将两个设计结合起来看即可：随机分配5位患者到CBT和EMDR，在5周和6个月后分别评价他们的STAI结果，如表3：

疗法和时间都作为引子是，我们既可以分析疗法的影响和时间的影响，又可以分析疗法和时间的交互影响。前两个称作主效应，交互部分称作交互效应。
当设计包含2个以上的因子的实验，就是因素方差分析设计，比如2个因子就是双因素方差分析，3个因子就是三因素方差分析，以此类推。如果因子的设计同时包括组间和组内的因子，又称作呼和模型方差分析。如表3即为典型的双因素混合模型方差分析。
这个例子需要做3次F检验。疗法因素一次，时间因素一次，两者交互因素一次。若疗法结果显著，说明CBT,EMDR对焦虑症的治疗效果有差异，若时间结果显著，则表示焦虑度从5周到6个月发生了显著变化，若两者交互效应显著，说明2种疗法随时间的变化对焦虑症治疗影响有差异。换句话说，焦虑度从5周到6个月的改变程度在两种疗法间是有差异的。

下面继续拓展。众所周知，抑郁症对病症治疗有影响，而且抑郁症和焦虑症常同时出现。即使受试者被随机分配到不同的治疗方案中，在研究开始时，两组的疗法中的患者抑郁水平就可能不同，任何治疗后的差异都有可能是最初的抑郁水平不同导致的，而不是实验中的操作问题。抑郁症也可以解释因变量的组间差异，因此它常被成为混淆因素(confounding factor)。由于你对抑郁症不感兴趣，它也被称作干扰变数(nuisance variable)。
假设招募患者时使用抑郁症的自我评测报告，比如白氏抑郁症量表（BDI），记录了他们的抑郁水平，那么你可以在评测两发类型的影响前，对任何抑郁水平的组间差异进行统计性调整。那么这个BDI就成为协变量，该设计就成为协方差分析（ANCOVA）。
以上设计只记录了单个因变量情况（STAI），为增强研究的有效性，可以对焦虑症进行其他的测量（比如家庭评分、医师评分，以及焦虑症对日常行为的影响评价）。当因变量不止一个时，设计被称作多元方差分析（MANOVA），若也存在协变量，那么就叫多元协方差分析（MANCOVA）。

下面就开始用R拟合ANOVA/ANCOVA/MANOVA模型了。

ANOVA模型拟合

虽然ANOVA和回归方法都是独立发展而来，但是从函数形式上来看，他们都是广义线性模型的特例。用之前提到的lm()函数也能分析ANOVA模型，不过，今天主要使用avo()函数。两个函数的结果是相同的，只是函数的展示结果的格式不同，后面可以比较一下。

1. aov()函数

aov()函数需要的数据是一个dataframe，表4列举了表达式中可以使用的特殊符号，其中y是因变量，字母A 、B、C代表因子。

表4 R表达式中的特殊符号

表5就列举了一些常见的研究设计表达式。在表5中，小写字母表示定量变量，大写字母表示组别因子，Subject 是对受试者独有的标识变量。

表5 常见研究设计的表达式

表5就列举了一些常见的研究设计表达式。在表5中，小写字母表示定量变量，大写字母表示组别因子，Subject 是对受试者独有的标识变量。
后面将深入讨论几个实验设计的例子

2.表达式中各项的顺序

表达式中效应的顺序在2中情况下会造成影响：（a）因子不止一个，并且是非平衡设计；（b）存在协变量。出现任意一种，等式右边的变量都与其他变量相关。那么就无法清晰地划分他们对因变量的影响。例如，对于双因素方差分析，若不同的处理方式中的观测数不同，那么模型y~A*B与模型y~B*A的结果不同。
R默认类型I（序惯型）方法计算ANOVA效应（参考补充内容“顺序很重要！”）第一个模型可以这样写：y~A+B+A:B。R中的ANOVA表的结果将评价：

• A对y的影响；
• 控制A时，B对y的影响；
• 控制A,B的主效应时，A与B的交互效应。

顺序很重要

• 当自变量与其他自变量或者协变量相关时，没有明确的方法可以评价自变量对因变量的贡献。例如，含因子A，B和因变量y的双因素不平衡因子的设计，有三种效应：A和B的主效应，A和B的交互效应。假设你正使用如下表达式对数据进行建模：
  Y~ A + 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：B交互项根据A和B调整。
  R默认调用类型I方法，其他软件（比如SAS和SPSS）默认调用类型III方法。

样本大小越不平衡，效应项的顺序对结果的影响越大。一般来说，越基础性的效应越需要放在表达式前面。具体来讲，首先是协变量，然后是主效应，接着是双因素的交互项，再接着是三因素的交互项，以此类推。对于主效应，月基础性的变量月应放在表达式的前面，因此性别要放在处理方式之前。有一个基本的准则：若研究设计不是正交的（也就是说，因子和/协变量相关），一定要谨慎设置效应的顺序。
在讲解具体的例子前，请注意car包中的Anova()函数（不是标准的anova()函数，R中严格区分大小写）提供了使用类型II或类型III方法的选项，而aov()函数使用的是类型I方法，若想使结果与其他软件(如SAS和SPSS)提供的结果保持一致，可以使用Anova()函数,具体细节可参考help(Anova,package="car")。

单因素方差分析

一般感兴趣的是比较分类因子定义的两个或多个组别中的因变量均值。以multcomp包中的cholesterol数据集为例，50个患者均接受降低胆固醇的药物治疗（trt）五种疗法汇总的一种疗法。其中3种治疗条件使用药物相同，分别是20mg qd,10mg bid,5mg qid。剩下的两种方式（drugD和drugE）代表候选药物。哪种药物疗法降低胆固醇（响应变量）最多呢？分析过程如下：

library(multcomp)
attach(cholesterol)
table(trt) #各组样本大小
aggregate(response,by =list(trt),FUN =mean) #各组均值
aggregate(response,by =list(trt),FUN =sd) #各组标准差
fit <- aov(response ~ trt) #检验组间差异（ANOVA）
summary(fit) 
library(gplots) #绘制各组均值及置信区间的图形
plotmeans(response ~ trt, xlab = "Treatment", ylab = "Response",
          main="Mean Plot\nwith 95 % CI")

运行结果如下

> library(multcomp)
Loading required package: mvtnorm
Loading required package: survival
Loading required package: TH.data
Loading required package: MASS

Attaching package: ‘TH.data’

The following object is masked from ‘package:MASS’:

    geyser

> attach(cholesterol)
> table(trt)
trt
 1time 2times 4times  drugD  drugE 
    10     10     10     10     10 
> aggregate(response,by =list(trt),FUN =mean)
  Group.1        x
1   1time  5.78197
2  2times  9.22497
3  4times 12.37478
4   drugD 15.36117
5   drugE 20.94752
> aggregate(response,by =list(trt),FUN =sd)
  Group.1        x
1   1time 2.878113
2  2times 3.483054
3  4times 2.923119
4   drugD 3.454636
5   drugE 3.345003
> fit <- aov(response ~ trt)
> summary(fit)
            Df Sum Sq Mean Sq F value   Pr(>F)    
trt          4 1351.4   337.8   32.43 9.82e-13 ***
Residuals   45  468.8    10.4                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

图1，5种降低胆固醇药物疗法的均值，含95%的置信区间

多重比较

虽然ANOVA对各种疗法的F检验表明5种药物疗法效果不同，但是并没有告诉你哪种疗法与其他疗法不同。多重比较可以解决这个问题。例如，TukeyHSD()函数提供了各组均值差异的成对检验：

> TukeyHSD(fit)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = response ~ trt)

$trt
                  diff        lwr       upr     p adj
2times-1time   3.44300 -0.6582817  7.544282 0.1380949
4times-1time   6.59281  2.4915283 10.694092 0.0003542
drugD-1time    9.57920  5.4779183 13.680482 0.0000003
drugE-1time   15.16555 11.0642683 19.266832 0.0000000
4times-2times  3.14981 -0.9514717  7.251092 0.2050382
drugD-2times   6.13620  2.0349183 10.237482 0.0009611
drugE-2times  11.72255  7.6212683 15.823832 0.0000000
drugD-4times   2.98639 -1.1148917  7.087672 0.2512446
drugE-4times   8.57274  4.4714583 12.674022 0.0000037
drugE-drugD    5.58635  1.4850683  9.687632 0.0030633

> par(las=2)
> par(las=2)
> par(mar=c(5,8,4,2))
> plot(TukeyHSD(fit))

可以看到，1time和2time的均值差异不显著(p=0.138),而1time和4time的差异非常显著(p<0.001)。
成对比较图形如图2所示。第一个par预计用来旋转轴标签，第二个用来增大左边界的面积，可使标签摆放更美观，具体参数设置可查阅帮助文档。图形中执行区间包含0的疗法说明差异不显著(p>0.5)。

图2 Tukey HSD均值成对比较图

multcomp包中的glht()函数提供了多重均值比较更为全面的方法，既适用于线性模型，也适用于广义线性模型。下面的代码重现了Tukey HSD检验，并用一个不同的图形对结果进行展示:

library(multcomp)
par(mar=c(5,4,6,2))
tuk <- glht(fit,linfct=mcp(trt="Tukey"))
plot(cld(tuk,level = .05,col="lightgrey"))

图3 multcomp包中的Tukey HSD检验

评估检验的假设条件

我们对于结果的信心依赖于做统计检验时数据满足假设条件的程度。单因素方差分析中，我们假设因变量服从正态分布，各组方差相等。可以使用Q-Q图来检验正态性假设：

> library(car)
> qqPlot(lm(response ~ trt, data =cholesterol),
+        simulate=TRUE, main="Q-Q Plot",labels=FALSE)
[1] 19 38

结果如图4所示：

图4 正态性检验

> bartlett.test(response ~ trt,data = cholesterol)

    Bartlett test of homogeneity of variances

data:  response by trt
Bartlett's K-squared = 0.57975, df = 4, p-value = 0.9653

Barlett检验表明5组的方差并没有显著不同（p=0.97）。其他检验如Fligner-Killeen检验（fligner.test()）函数和Brown-Forsythe检验（HH包中的hov()函数）此处没有做演示，但他们获得的结果与Barlett检验相同。
不过，方差齐性分析对离群点非常敏感。可利用car包中的outlierTest()函数来检测离群点：

> outlierTest(fit)
No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:
   rstudent unadjusted p-value Bonferroni p
19 2.251149           0.029422           NA

从输出结果来看，并没有证据说明胆固醇数据中含有离群点（当p>1时将产生NA）。因此根据Q-Q图、Barlett检验和离群点检验，该数据似乎可以用ANOVA模型拟合的很好。这些方法反过来增强了我们对于所得结果的信心。

单因素协方差分析

单因素协方差分析（ANCOVA），包含一个或多个定量的协变量。下面的例子来自multcomp包中的litter数据集。怀孕的小鼠被分为4组，每组接受不同剂量（0，5，50，或500）的药物处理。产下的幼崽的体重均值为因变量，怀孕的时间为协变量。如下所示

> data("litter",package = "multcomp")
> attach(litter)
> table(dose)
dose
  0   5  50 500 
 20  19  18  17 
> aggregate(weight, by=list(dose),FUN=mean)
  Group.1        x
1       0 32.30850
2       5 29.30842
3      50 29.86611
4     500 29.64647
> fit <- aov(weight ~ gesttime + dose)
> summary(fit)
            Df Sum Sq Mean Sq F value  Pr(>F)   
gesttime     1  134.3  134.30   8.049 0.00597 **
dose         3  137.1   45.71   2.739 0.04988 * 
Residuals   69 1151.3   16.69                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

可以看到，每种剂量下所产的幼崽数并不相同：0剂量时（未用药）产仔20个，500剂量时产仔17个。再看均值，发现未用药组体重均值最高32.3，ANCOVA的F检验表明：（a）怀孕时间与幼崽出生体重相关；（b）控制怀孕时间，药物剂量与出生体重相关。控制怀孕时间，确实发现每种药物剂量下幼崽出生体重均值不同。由于使用了协变量，你可能想要获取调整的组均值，即去除协变量效应后的组均值。可使effects包中的effects()函数来计算调整的均值：

> library(effects)
Registered S3 methods overwritten by 'lme4':
  method                          from
  cooks.distance.influence.merMod car 
  influence.merMod                car 
  dfbeta.influence.merMod         car 
  dfbetas.influence.merMod        car 
lattice theme set by effectsTheme()
See ?effectsTheme for details.
> effect("dose",fit)

 dose effect
dose
       0        5       50      500 
32.35367 28.87672 30.56614 29.33460 

本例中，调整的均值与aggregate()函数得出的为调整的均值类似，但并非所有的情况都是如此。总之effects包为复杂的研究设计提供了强大的计算调整均值的方法，并能将结果可视化，更多细节参考帮助文档。
同前面的ANOVA例子一样，剂量的F检验虽然表明了不同处理方式幼崽的体重均值不同，但无法告知我们哪种处理方式与其他方式不同。同样需要进行假设检验。
假设对未用药条件与其他三种用药条件影响是否不同感兴趣，检验结果如下：

> library(multcomp)
> contrasts <- rbind("no drug vs. drug"=c(3,-1,-1,-1))
> contrast <- rbind("no drug vs. drug"=c(3,-1,-1,-1))
> summary(glht(fit,linfct=mcp(dose=contrast)))

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: aov(formula = weight ~ gesttime + dose)

Linear Hypotheses:
                      Estimate Std. Error t value Pr(>|t|)  
no drug vs. drug == 0    8.284      3.209   2.581    0.012 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

对照c(3,-1,-1,-1)设定第一组和其他三组的均值进行比较。假设检验的t统计量（2.581）在p<0.05水平下显著，因此，可以得出未用药组比其他用药条件下的出生体重高的结论。其他对照可用rbind()函数添加。

评估检验的假设条件

ANCOVA与ANOVA相同，都需要正态性和同方差性假设，可以用上面的步骤来检验这些假设条件。另外ANCOVA还假定回归斜率相同。本例中，假定四个处理组通过怀孕时间来预测出生体重的回归斜率都相同。ANCOVA模型包含怀孕时间X剂量的交互项时，可对回归斜率的同质性进行检验。交互效应若显著，则意味着时间和幼崽出生体重间的关系依赖于药物剂量的水平。
检验斜率的同质性：

> library(multcomp)
> fit2 <- aov(weight ~ gesttime*dose,data = litter)
> summary(fit2)
              Df Sum Sq Mean Sq F value  Pr(>F)   
gesttime       1  134.3  134.30   8.289 0.00537 **
dose           3  137.1   45.71   2.821 0.04556 * 
gesttime:dose  3   81.9   27.29   1.684 0.17889   
Residuals     66 1069.4   16.20                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

上面结果提示交互效应不显著，支持了斜率相等的假设。若假设不成立，可以尝试变换协变量或因变量，或使用能对每个斜率独立解释的模型，或使用不需要假设斜率同质性的非参数ANCOVA方法。 sm包中的sm.ancova()函数为后者提供解释方案。

结果可视化

HH包中ancova()函数可以绘制因变量、协变量和因子直接的关系图。

> library(HH)
Loading required package: lattice
Loading required package: grid
Loading required package: latticeExtra
Loading required package: RColorBrewer
Loading required package: gridExtra

Attaching package: ‘HH’

The following objects are masked from ‘package:car’:

    logit, vif

The following object is masked from ‘package:gplots’:

    residplot

> ancova(weight ~ gesttime + dose, data=litter)
Analysis of Variance Table

Response: weight
          Df  Sum Sq Mean Sq F value   Pr(>F)   
gesttime   1  134.30 134.304  8.0493 0.005971 **
dose       3  137.12  45.708  2.7394 0.049883 * 
Residuals 69 1151.27  16.685                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

图5 四种药物处理组的怀孕时间和出生体重的关系图

图5提示，用怀孕时间来预测出生体重的回归线相互平行，知识截距不同。随着怀孕时间的增加，幼崽出生体重也会增加。另外，还可以看到0剂量组截距项最大，5剂量组截距项最小，由于你上面的设置，直线会保持平行，若用ancova(weight ~ gesttime*dose)，生成的图形将允许斜率和截距项依据组别而发生变化，这对可视化那些违背回归斜率同质性的实例非常有用。

双因素方差分析

在双因素方差分析中，受试者被分配到两因子的交叉类别组中。以基础安装中的ToothGrowth数据集为例，随机分配60只豚鼠，分别采用两种喂食方法（橙汁或维C），个喂食方法中抗坏血酸含量有三种水平（0.5mg、1mg或2mg）,每种处理方式组合都被分配10只豚鼠。牙齿长度为因变量，分析如下：

> attach(ToothGrowth)
The following object is masked _by_ .GlobalEnv:

    dose

The following object is masked from litter:

    dose

> table(supp,dose)
    dose
supp 0.5  1  2
  OJ  10 10 10
  VC  10 10 10
> aggregate(len,by=list(supp,dose),FUN=mean)
  Group.1 Group.2     x
1      OJ     0.5 13.23
2      VC     0.5  7.98
3      OJ       1 22.70
4      VC       1 16.77
5      OJ       2 26.06
6      VC       2 26.14
> aggregate(len,by=list(supp,dose),FUN=sd)
  Group.1 Group.2        x
1      OJ     0.5 4.459709
2      VC     0.5 2.746634
3      OJ       1 3.910953
4      VC       1 2.515309
5      OJ       2 2.655058
6      VC       2 4.797731
> dose <- factor(dose)
> fit <- aov(len ~ supp*dose)
> summary(fit)
            Df Sum Sq Mean Sq F value   Pr(>F)    
supp         1  205.4   205.4  15.572 0.000231 ***
dose         2 2426.4  1213.2  92.000  < 2e-16 ***
supp:dose    2  108.3    54.2   4.107 0.021860 *  
Residuals   54  712.1    13.2                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> detach(ToothGrowth)
> interaction.plot(dose,ToothGrowth$supp, ToothGrowth$len,type = "b",
+                  col = c("red","blue"),pch = c(16,18),
+                  main = "Interaction between Dose and Supplement Type")

可见主效应（supp,dose）和交互效应都非常显著。此处可以用 interaction.plot函数来可视化，结果如下图所示

图6 喂食方法和剂量对牙齿生长的交互作用。用interaction.plot()函数绘制了牙齿长度的均值

library(gplots)
plotmeans(len ~ interaction(supp, dose, sep = " "),
          connect = list(c(1,3,5),c(2,4,6)),
          col = c("red","darkgreen"),
          main ="Interaction Plot with 95% CIs",
          xlab = "Treatment and Dose Combination")

图7 喂食方法和剂量对牙齿生长的交互作用。用plotmeans()函数绘制95%的置信区间的牙齿长度均值

还可用用HH包的interaction2wt()函数来可视化，图形对任意顺序的因子设计的主效应和交互效应都会进行展示。

library(HH)
interaction2wt(len ~ supp*dose)

图8 ToothGrowth数据集的主效应和交互效应。图形由interaction2wt()函数创建

显然三种方法我认为第三种更适合解释任意复杂度设计。

重复测量方差分析

即受试者测量不止一次。重点关注一个常见的含组内和组间因子重复测量方差分析，示例数据来源于生态学领域，研究方向是生命系统的生理和生化过程如何响应环境因素的变异（此为应对全球变暖的一个非常重要的研究领域）。
基础安装包中的数据集包含了南方与北方的牧草类植物Echinochloa crus-galli的寒冷容忍度研究结果，在某浓度环境中，对寒带植物的光合作用率进行了比较。研究用植物一般来自加拿大的魁北克省与另一半来自美国的密西西比州。
首先，关注的是寒带植物，因变量为吸收量，单位为ml/L，自变量是植物类型，Type和7种水平（ sec）的浓度conc。另外，Type是组间因子，conc是组内因子。

CO2$conc <- factor(CO2$conc)
w1b1 <- subset(CO2,Treatment =='chilled')
fit <- aov(uptake ~ conc*Type +Error(Plant/(conc)),w1b1)
summary(fit)
par(las=2)
par(mar = c(10,4,4,2))
with(w1b1,interaction.plot(conc,Type,uptake,
                           type = "b",col = c("red","blue"),pch = c(16,18),
                           main ="Interaction Plot for Plant Type and Concentration"))
boxplot(uptake ~ Type*conc,data = w1b1,col=(c("gold","green")),
        main = "Chilled Quebec and Mississippi Plants",
        ylab = "Carbon dioxide uptake rate (umol/m^2 sec)")

图9 CO2浓度和植物类型对其吸收的交互影响

图10 交互效应

注意：通常处理的数据集是宽格式，即列是变量，行是观测值，而且一行一个授试对象，litter数据框就是一个很好的例子，不过在处理长格式数据才能拟合模型，在长格式中，因变量的每次测量都要放到他独有的行中，CO2数据集就是这种形式，可以用reshape包等进行格式转换。在分析自己的数据时候先把自己的数据做成类似的格式即可套用模型。

混合模型涉及到各种方法

在分析CO2的例子时，我们使用了传统的重复测量方差分析。该方法假设任一组内因子的协方差矩阵为球形，并且任意组内因子两水平间的方差之差都相等。但是现实中这种假设不可能满足，于是衍生了一系列备选方案：

• 使用lme4包的lmer()函数拟合线性混合模型
• 使用car包的Anova()函数调整传统检验统计量以弥补球形假设的不满足（例如Geisser-Greenhouse校正）
• 使用nlme 包中的gls()函数拟合给定方差-协方差结构的广义最小二乘模型。
• 用多元方差分析对重复测量数据进行建模
  以上方法还需进一步探索。

目前为止都是对单个因变量的情况进行分析的，多个结果变量的设计常用方法：多元方差分析。

多元方差分析

MANOVA，以MASS包中的UScereal数据集演示，研究美国谷物中的卡路里、脂肪和糖含量是否会因为储存架位置不同而发生变化；其中1代表底层货架，2代表中层货架，3代表顶层货架。卡路里、脂肪和糖含量是因变量，货架是3个水平的自变量。分析过程如下：

rm(list = ls())
library(MASS)
attach(UScereal)
shelf <- factor(shelf)
y <- cbind(calories, fat, sugars)
aggregate(y, by=list(shelf),FUN=mean)
cov(y)
fit <- manova(y ~ shelf)
summary(fit)
summary.aov(fit)

结果为：

> summary(fit)
          Df Pillai approx F num Df den Df    Pr(>F)    
shelf      2 0.4021   5.1167      6    122 0.0001015 ***
Residuals 62                                            
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> summary.aov(fit)
 Response calories :
            Df Sum Sq Mean Sq F value    Pr(>F)    
shelf        2  50435 25217.6  7.8623 0.0009054 ***
Residuals   62 198860  3207.4                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response fat :
            Df Sum Sq Mean Sq F value  Pr(>F)  
shelf        2  18.44  9.2199  3.6828 0.03081 *
Residuals   62 155.22  2.5035                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response sugars :
            Df  Sum Sq Mean Sq F value   Pr(>F)   
shelf        2  381.33 190.667  6.5752 0.002572 **
Residuals   62 1797.87  28.998                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

由于多元检验是显著的，可以使用summary.aov(fit)对每一个变量做单因素方差分析。从结果二月看到，三组中每种营养成分的测量值都不同，另外还可以用均值比较步骤（TukeyHSD）来判断对于每一个因变量，哪种货架与其他货架都是不同的。

评估假设检验

若有一个px1的多元正态随机向量,其均值为,协方差矩阵为反3，那么,,的马氏距离的平方服从自由度为p的卡方分布。Q-Q图展示卡方分布的分位数，横纵坐标分别是样本量与马氏距离平方值。如果点全部落在斜率为1，截距为0的直线上，则表明数据服从多元正态分布。
分析如下：

center <- colMeans(y)
n <- nrow(y)
p <- ncol(y)
cov <- cov(y)
d <- mahalanobis(y,center,cov)
coord <- qqplot(qchisq(ppoints(n),df=p),
                d, main = "Q-Q Plot Assessing Multivariate Normality",
                ylab = "Mahalanobis D2")
abline(a=0,b=1)
identify(coord$x,coord$y,labels = row.names(UScereal))

图11 检验多元正态性的Q-Q图

数据不完全符合多元正态分布，有2 个离群点，如果符合的情况下可以用identify()交互性的对图中的点进行鉴别，该数据不符合，也可以删掉这2个离群点重新分析
可以使用mvoutlier包中的ap.plot()函数来检验多元离群点

library(mvoutlier)
outliers <- aq.plot(y)
outliers

结果是：

> library(mvoutlier)
Loading required package: sgeostat
Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2
sROC 0.1-2 loaded
> outliers <- aq.plot(y)
Projection to the first and second robust principal components.
Proportion of total variation (explained variance): 0.9789888
> outliers
$outliers
 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE
[17] FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE
[33] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE
[49] FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[65] FALSE

image.png

稳健多元方差分析

如果多元正态性或者方差-协方差均值假设都不满足，或者你担心多元离群点，那么可以考虑用稳健或非参数版本的MANOVA检验。稳健单因素MANOVA可以通过rrcov包中的Wilks.test()函数实现。vegan包中的adonis()函数则提供了非参数MANOVA的同等形式。

> library(rrcov)
Loading required package: robustbase
Scalable Robust Estimators with High Breakdown Point (version 1.4-7)

> Wilks.test(y,shelf,methods="mcd")

    One-way MANOVA (Bartlett Chi2)

data:  x
Wilks' Lambda = 0.63794, Chi2-Value = 27.42, DF = 6.00, p-value = 0.0001208
sample estimates:
  calories       fat    sugars
1 119.4774 0.6621338  6.295493
2 129.8162 1.3413488 12.507670
3 180.1466 1.9449071 10.856821

从结果来看，文件检验对离群点和违反MANOVA假设的情况不敏感，而且再一次验证了存储在货架顶部、中部和底部的谷物营养成分含量不同。

用回归来做ANOVA

ANOVA和回归都是广义线性模型的特例。因此，本次的所有设计都可以用lm()函数来分析。但是为了更好的理解输出的结果，需要弄明白在拟合模型时，R是如何处理类别变量的。
下面来比较5种降低胆固醇药物疗法（trt）的影响。

library(multcomp)
levels(cholesterol$trt)
fit.aov <- aov(response ~ trt,data = cholesterol)
summary(fit.aov)
fit.lm <- lm(response ~ trt, data = cholesterol)
summary(fit.lm)

结果为：

> summary(fit.aov)
            Df Sum Sq Mean Sq F value   Pr(>F)    
trt          4 1351.4   337.8   32.43 9.82e-13 ***
Residuals   45  468.8    10.4                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> fit.lm <- lm(response ~ trt, data = cholesterol)
> summary(fit.lm)

Call:
lm(formula = response ~ trt, data = cholesterol)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.5418 -1.9672 -0.0016  1.8901  6.6008 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    5.782      1.021   5.665 9.78e-07 ***
trt2times      3.443      1.443   2.385   0.0213 *  
trt4times      6.593      1.443   4.568 3.82e-05 ***
trtdrugD       9.579      1.443   6.637 3.53e-08 ***
trtdrugE      15.166      1.443  10.507 1.08e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.227 on 45 degrees of freedom
Multiple R-squared:  0.7425,    Adjusted R-squared:  0.7196 
F-statistic: 32.43 on 4 and 45 DF,  p-value: 9.819e-13

因为线性模型要求预测变量是数值型，当lm()函数碰到因子时，它会用一系列与因子水平相对应的数值型对照变量来代替因子。如果因子有个水平，将会创建个对照变量，R提供了5种创建对照变量的内置方法见表6，你也可以自己重新创建。默认情况下，对照处理用于无序因子，正交多项式用于有序因子。

表6 内置对照

以对照（treament contrast）为例，因子的第一个水平变成了参考组，随后的变量都以它为标准。可以通过contrasts()函数查看他的编码过程。

> contrasts(cholesterol$trt)
       2times 4times drugD drugE
1time       0      0     0     0
2times      1      0     0     0
4times      0      1     0     0
drugD       0      0     1     0
drugE       0      0     0     1

若患者处于drugD条件下，变量drugD等于1，其他变量2times,4times和drugE都等于0。无需列出第一组的变量值，因为其他4个变量都为0。这已经说明患者处于1time条件。
在前面一段代码中，变量trt2times表示水平1time和2times的一个对照。类似地，trt4times，是1time和4times的一个对照。其余以此类推。从输出的概率来看，各药物条件与第一组（1time）显著不同。
通过设定contrasts选项，可以修改lm()中默认的对照方法。例如，若使用Helmert对照：fit.lm <- lm(response ~ trt, data=cholesterol, contrasts = "contr.helmert")还能通过options()函数修改R会话中的默认对照方法，例如，options(contrasts =c("contr.SAS","contr.helmert"))
设定无序因子的默认对比方法为contr.SAS，有序因子的默认对比方法为contr.helmert。虽然我们一直都是在线性模型范围中讨论对照方法的使用，但是在R中，你完全可以将其应用到其他模型中，包括广义线性模型。

小结

本教程回顾了基本实验和准实验设计的分析方法，包括ANOVA/ANCOVA/MANOVA。然后通过组内和组间设计的示例介绍了基本方法的使用，如单因素ANOVA、单因素ANCOVA、双因素ANOVA、重复测量ANOVA和单因素MANOVA。
除了这些基本分析，我们还回顾了模型的假设检验，以及应用多重比较过程来进行综合检验的方法。最后对各种结果可视化方法也进行了探索。如果你对用R分析DOE（Design Of Experiment,实验设计）感兴趣，请参阅说明书。