方差分析

使用的包有car、gplots、HH、rrcov和mvoutlier。

方差分析的基本原理是认为不同处理组的均数间的差别基本来源有两个：

1. 实验条件，即不同的处理造成的差异，称为组间差异。用变量在各组的均值与总均值之偏差平方和的总和表示，记作SSb，组间自由度dfb。

2. 随机误差，如测量误差造成的差异或个体间的差异，称为组内差异，用变量在各组的均值与该组内变量值之偏差平方和的总和表示， 记作SSw，组内自由度dfw。

总偏差平方和 SSt = SSb + SSw。

组内SSw、组间SSb除以各自的自由度(组内dfw=n-m，组间dfb=m-1，其中n为样本总数，m为组数)，得到其均方MSw和MSb，一种情况是处理没有作用，即各组样本均来自同一总体，MSb/MSw≈1。另一种情况是处理确实有作用，组间均方是由于误差与不同处理共同导致的结果，即各样本来自不同总体。那么，MSb>>MSw(远远大于)。 MSb/MSw比值构成F分布。用F值与其临界值比较，推断各样本是否来自相同的总体。

F分布表

术语速成

单因素方差分析

多因素方差分析

因素（组内、组间、混淆）

协方差分析（单、多）

以焦虑症治疗为例，现有两种治疗方案：认知行为疗法（CBT）和眼动脱敏再加工法（EMDR）。我们招募10个焦虑症患者作为志愿者，随机分配一半的人接受为期五周的CBT，另外一半接受为期五周的EMDR。在治疗结束时，要求每个患者都填写状态特质焦虑问卷（STAI），也就是一份焦虑度测量的自我评测报告。

单因素组间方差分析

方差分析主要通过F检验来进行效果评测，若治疗方案的F检验显著，则说明五个星期后两种疗法的STAI得分均值不同。

单因素组内方差分析

当时间的F检验显著时，说明患者的STAI得分均值在五周和六个月间 发生了改变。

含组间和组内因子的双因素方差分析

疗法（therapy）和时间（time）都作为因子时，我们既可分析疗法的影响（时间跨度上的平 均）和时间的影响（疗法类型跨度上的平均），又可分析疗法和时间的交互影响。前两个称作主效应，交互部分称作交互效应。

本例中，你将做三次F检验：疗法因素一次，时间因素一次，两者交互因素一次。若疗法结 果显著，说明CBT和EMDR对焦虑症的治疗效果不同；若时间结果显著，说明焦虑度从五周到六 个月发生了变化；若两者交互效应显著，说明两种疗法随着时间变化对焦虑症治疗影响不同（也就是说，焦虑度从五周到六个月的改变程度在两种疗法间是不同的）。

ANOVA 模型拟合

aov()函数

aov()函数的语法为aov(formula, data = dataframe)

R表达式中的特殊符号

小写字母表示定量变量，大写字母表示组别因子，Subject是对被试者独有的标识变量。

常见研究设计的表达式

随机化区组：依据数学上概率的原理，将被试材料按相等机会原则分组。理论上可使不同组的被试材料除实验处理之外，其他无关变量保持相等，可弥补配对法顾此失彼的特点，是控制无关变量较好的方法。

表达式中各项的顺序

表达式中效应的顺序在两种情况下会造成影响：

(a)因子不止一个，并且是非平衡设计；

(b)存在协变量。出现任意一种情况时，等式右边的变量都与其他每个变量相关。

此时，我们无法清晰地划分它们对因变量的影响。例如，对于双因素方差分析，若不同处理方式中的观测数不同，那么模型y~ A∗BA∗B与模型y ~ B∗AB∗A的结果不同。 R默认类型I（序贯型）方法计算ANOVA效应。第一个模型可以这样写：y ~ A + B + A:B。

R中的ANOVA表的结果将评价：

A对y的影响；

控制A时，B对y的影响；

控制A和B的主效应时，A与B的交互效应。

R表达式顺序

样本大小越不平衡，效应项的顺序对结果的影响越大。一般来说，越基础性的效应越需要放在 表达式前面。具体来讲，首先是协变量，然后是主效应，接着是双因素的交互项，再接着是三因素的交互项，以此类推。对于主效应，越基础性的变量越应放在表达式前面。

car包中的Anova()函数（不要与标准anova()函数混淆）提供了使用类型II或类型III方法的选项，而aov()函数使用的是类型I方法。

单因素方差分析

单因素方差分析中，你感兴趣的是比较分类因子定义的两个或多个组别中的因变量均值。以multcomp包中的cholesterol数据集为例（取自Westfall，Tobia，Rom，Hochberg，1999），50个患者均接受降低胆固醇药物治疗（trt）五种疗法中的一种疗法。其中三种治疗条件使用药物相同，分别是20 mg一天一次（1time）、10mg一天两次（2times）和5mg一天四次（4times）。剩下的两种方式（drugD和drugE）代表候选药物。哪种药物疗法降低胆固醇（响应变量）最多呢？

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)#检验组间差异（anova）

##             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

library(gplots)

## ## Attaching package: 'gplots'

## The following object is masked from 'package:stats':## ##     lowess

plotmeans(response~trt,xlabl="Treatment",ylab = "Response",main="Mean Plot\nwith 95% CI")#绘制各组均值及其置信区间的图形

## Warning in plot.window(...): "xlabl"不是图形参数

## Warning in plot.xy(xy, type, ...): "xlabl"不是图形参数

## Warning in axis(side = side, at = at, labels = labels, ...): "xlabl"不是图## 形参数## Warning in axis(side = side, at = at, labels = labels, ...): "xlabl"不是图## 形参数

## Warning in box(...): "xlabl"不是图形参数

## Warning in title(...): "xlabl"不是图形参数

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "xlabl"不是图形参数

## Warning in axis(1, at = 1:length(means), labels = legends, ...): "xlabl"不## 是图形参数

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "xlabl"不是图形参数

detach(cholesterol)

从输出结果可以看到，每10个患者接受其中一个药物疗法。均值显示drugE降低胆固醇最 多，而1time降低胆固醇最少，各组的标准差相对恒定，在2.88到3.48间浮动。ANOVA对治 疗方式（trt）的F检验非常显著（p<0.0001），说明五种疗法的效果不同。

多重比较

虽然ANOVA对各疗法的F检验表明五种药物疗法效果不同，但是并没有告诉你哪种疗法与其他疗法不同，多重比较可以解决这个问题。 TukeyHSD()函数提供了对各组均值差异的成对检验。但要注意TukeyHSD()函数与本章使用的HH包存在兼容性问题：若载入HH包，TukeyHSD()函数将会失效。对于上例，使用detach(“package::HH”)将它从搜寻路径中删除，然后再调用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

原假设应为差异不显著，置信区间为95%，可以看到，1time和2times的均值差异不显著（p=0.138），而1time和4times间的差异非 常显著（p<0.001）。

par(last=2)

## Warning in par(last = 2): "last"不是图形参数

plot(TukeyHSD(fit))

成对比较图形如图所示。第一个par语句用来旋转轴标签，第二个用来增大左边界的面积，可使标签摆放更美观（par选项参见第3章）。原假设为差异显著，置信区间为95%，图形中置信区间包含0的疗法说明差异不显著（p>0.05）。

multcomp包中的glht()函数提供了多重均值比较更为全面的方法，既适用于线性模型，也适用于广义线性模型。

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

有相同字母的组（用箱线图表示）说明均值差异不显著。可以看到，1time和2times差异不显 著（有相同的字母a），2times和4times差异也不显著（有相同的字母b），而1time和4times差异显著（它们没有共同的字母）。

从结果来看，使用降低胆固醇的药物时，一天四次5 mg剂量比一天一次20 mg剂量效果更佳， 也优于候选药物drugD，但药物drugE比其他所有药物和疗法都更优。

评估检验的假设条件

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

library(car)qqPlot(lm(response~trt,data = cholesterol),simulate=TRUE,main="Q-QPlot",labels=FALSE)

注意qqPlot()要求用lm()拟合。数据落在95%的置信区间范围内，说明满足正态性假设。

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

Bartlett检验表明五组的方差并没有显著不同（p=0.97）。

但方差齐性分析对离群点非常敏感。可利用car包中的outlierTest()函数来检测离群点：

library(car)outlierTest(fit)

## ## No Studentized residuals with Bonferonni p < 0.05## Largest |rstudent|:##    rstudent unadjusted p-value Bonferonni p## 19 2.251149           0.029422           NA

从输出结果来看，并没有证据说明胆固醇数据中含有离群点（p<1）。因此根 据Q-Q图、Bartlett检验和离群点检验，可以验证数据可以用ANOVA模型拟合得很好。

单因素协方差分析

下面的例子来自于multcomp包中的litter数据集。怀孕小鼠被分为四个小组，每个小组接受不同剂量（0、5、50或500）的药物处理。产下幼崽的体重值为因变量，怀孕时间为协变量。

data("litter",package = "multcomp")attach(litter)table(dose)#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

利用table()函数，可以看到每种剂量下所产的幼崽数并不相同：0剂量时（未用药）产崽 20个，500剂量时产崽17个。再用aggregate()函数获得各组均值，可以发现未用药组幼崽体重均值最高（32.3）。

ANCOVA的F检验表明：

(a)怀孕时间与幼崽出生体重相关；

(b)控制怀孕时间，药物剂量与出生体重相关。控制怀孕时间，确实发现每种药物剂量下幼崽出生体重均值不同。

由于使用了协变量，可使用effects包中的effects()函数来获得去除协变量效应后的组均值：

library(effects)

## ## Attaching package: 'effects'

## The following object is masked from 'package:car':## ##     Prestige

effect("dose",fit)

## ##  dose effect## dose##        0        5       50      500 ## 32.35367 28.87672 30.56614 29.33460

和上一节的单因素ANOVA例子一样，剂量的F检验虽然表明了不同的处理方式幼崽的体重均值不同，但无法告知我们哪种处理方式与其他方式不同。同样，我们使用multcomp包来对所有均值进行成对比较。

假定你对未用药条件与其他三种用药条件影响是否不同感兴趣。

library(multcomp)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水平下显著，因此，可以得出未用药条件与其他三种用药条件影响不同的结论。

评估检验的假设条件

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

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

可以看到交互效应不显著，支持了斜率相等的假设。若假设不成立，可以尝试变换协变量或因变量，或使用能对每个斜率独立解释的模型。

结果可视化

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

library(HH)

ancova(weight~gesttime+dose,data=litter)

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

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

双因素方差分析

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

attach(ToothGrowth)

## 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.0 22.70## 4      VC     1.0 16.77## 5      OJ     2.0 26.06## 6      VC     2.0 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.0 3.910953## 4      VC     1.0 2.515309## 5      OJ     2.0 2.655058## 6      VC     2.0 4.797731

fit<-aov(len~supp*dose)summary(fit)

##             Df Sum Sq Mean Sq F value   Pr(>F)    ## supp         1  205.3   205.3  12.317 0.000894 ***## dose         1 2224.3  2224.3 133.415  < 2e-16 ***## supp:dose    1   88.9    88.9   5.333 0.024631 *  ## Residuals   56  933.6    16.7                     ## ---## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

可以看到主效应（supp和dose）和交互效应都非常显著。

有多种方式对结果进行可视化处理。此处可用interaction.plot()函数来展示双因素方差分析的交互效应。

interaction.plot(dose,supp,len,type = "b",col=c("red","blue"),pch = c(16,18),main="Interaction")

喂食方法和剂量对牙齿生长的交互作用。

还可以用gplots包中的plotmeans()函数来展示交互效应。

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")

以上图形都表明随着橙汁和维生素C中的抗坏血酸剂量的增加，牙齿长度变长。对于0.5 mg 和1 mg剂量，橙汁比维生素C更能促进牙齿生长；对于2mg剂量的抗坏血酸，两种喂食方法下牙齿长度增长相同。