引言：上一章我们学习了如何通过量化的预测变量来预测量化的响应变量的回归模型。本期我们将一起学习如何对类别型预测变量建立合适的统计模型进行分析。

当包含的因子是解释变量时，关注的重点通常会从预测转向组别差异的分析，这种分析法称作方差分析(ANOVA)。

9.1 ANOVA模型拟合

lm()函数也能分析ANOVA模型，但基本都使用aov()函数，两个函数的结果是等同的。

1. aov()函数

aov()函数的语法为 ：

aov(formula, data=dataframe)

表9-1 常见研究设计的表达式

设计

表达式

单因素ANOVA

y ~ A

含单个协变量的单因素 ANCOVA

y ~ x + A

双因素 ANOVA

y ~ A* B

含两个协变量的双因素 ANCOVA

y ~ x1 + x2 + A * B

随机化区组

y ~ B + A(B是区组因子)

单因素组内 ANOVA

y ~ A + Error(Subject/A)

含单个组内因子(W)和单个组间因子(B)的重复测量 ANOVA

y ~ B * W + Error(Subject/W)

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

2. 表达式中各项的顺序

当出现：(a)因子不止一个，并且是非平衡设计；(b)存在协变量

以上任意一种情况时，等式右边的变量都与其他每个变量相关，无法清晰地划分它们对因变量的影响。对于双因素方差分析，若不同处理方式中的观测数不同，那么模型y ~ A * B与模型y ~ B * A的结果不同。

R默认类型I(序贯型)方法计算ANOVA效应。R中的ANOVA表的结果将评价：

A对y的影响；

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

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

一般来说，越基础性的效应越需要放在表达式前面。

首先是协变量，然后是主效应，接着是双因素的交互项，再接着是三因素的交互项，以此类推。

对于主效应，越基础性的变量越应放在表达式前面，因此性别要放在处理方式之前。

基本准则：若研究设计不是正交的(也就是说，因子和/或协变量相关)，一定要谨慎设置效应的顺序。

9.2 单因素方差分析

单因素方差分析(one-way ANOVA)比较分类因子定义的两个或多个组别中的因变量均值。

例：multcomp包中的cholesterol数据集50个患者均接受降低胆固醇药物治疗(trt)五种疗法中的一种疗法。五种疗法为：相同药物分别20mg一天一次(1time)、10mg一天两次(2times)和5mg一天四次(4times)治疗及(drugD和drugE)候选药物治疗。

分析哪种药物疗法降低胆固醇(响应变量)最多：

library(multcomp)attach(cholesterol)#列出各组样本大小table(trt)trt1time 2times 4times  drugD  drugE10     10     10     10     10#每10个患者接受其中一个药物疗法#计算各组均值aggregate(response, by=list(trt), FUN=mean)Group.1        x1   1time  5.781972 2times 9.224973 4times 12.374784 drugD 15.361175 drugE 20.94752#均值显示drugE降低胆固醇最多，而1time降低胆固醇最少#计算各组标准差aggregate(response, by=list(trt), FUN=sd)Group.1        x1   1time 2.8781132  2times 3.4830543  4times 2.9231194 drugD 3.4546365 drugE 3.345003#各组的标准差相对恒定，在2.88到3.48间浮动#检验组间差异(ANOVA)fit 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

ANOVA对治疗方式(trt)的F检验非常显著(p<0.0001)，说明五种疗法的效果不同。

gplots包中的plotmeans()绘制带有置信区间的组均值图形：

library(gplots)plotmeans(response ~ trt, xlab="Treatment", ylab="Response",   #绘制各组均值及其置信区间的图形main="Mean Plotnwith 95% CI")detach(cholesterol)

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

1. 多重比较

ANOVA对各疗法的F检验表明五种药物疗法效果不同，多重比较可以知道哪种疗法与其他疗法不同。

TukeyHSD()函数提供了对各组均值差异的成对检验：

TukeyHSD(fit)Tukey multiple comparisons of means95% family-wise confidence levelFit: aov(formula = response ~ trt)$trtdiff        lwr       upr     p adj2times-1time   3.44300 -0.6582817  7.544282 0.13809494times-1time   6.59281  2.4915283 10.694092 0.0003542drugD-1time    9.57920  5.4779183 13.680482 0.0000003drugE-1time   15.16555 11.0642683 19.266832 0.00000004times-2times  3.14981 -0.9514717  7.251092 0.2050382drugD-2times   6.13620  2.0349183 10.237482 0.0009611drugE-2times  11.72255  7.6212683 15.823832 0.0000000drugD-4times   2.98639 -1.1148917  7.087672 0.2512446drugE-4times   8.57274  4.4714583 12.674022 0.0000037drugE-drugD    5.58635  1.4850683  9.687632 0.0030633

1time和2times的均值差异不显著(p=0.138)；1time和4times间的差异非常显著 (p<0.001)。

par(las=2)           #旋转轴标签par(mar=c(5,8,4,2))  #增大左边界的面积,可使标签摆放更美观plot(TukeyHSD(fit))

图9-2 Tukey HSD均值成对比较图

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

library(multcomp)par(las=1)par(mar=c(5,4,6,2))  #增大顶部边界面积,适合字母阵列摆放tuk

图9-3 multcomp包中的TukeyHSD检验

有相同字母的组(用箱线图表示)说明均值差异不显著。

2. 评估检验的假设条件

单因素方差分析中，假设因变量服从正态分布，各组方差相等。可以使用Q-Q图来检验正态性假设：

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

图9-4 正态性检验

注意qqPlot()要求用lm()拟合。

方差齐性检验：

可以通过如下代码来做Bartlett检验：

bartlett.test(response ~ trt, data=cholesterol)Bartlett test of homogeneity of variancesdata:  response by trtBartlett'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 Bonferroni p

没有证据说明胆固醇数据中含有离群点(当p>1时将产生NA)。根据Q-Q图、Bartlett检验和离群点检验，该数据可以用ANOVA模型拟合得很好。

9.3 单因素协方差分析

单因素协方差分析(ANCOVA)包含一个或多个定量的协变量。

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

library(multcomp)data(litter, package="multcomp")attach(litter)table(dose)dose0   5  50 50020  19  18  17#每种剂量下所产的幼崽数并不相同aggregate(weight, by=list(dose), FUN=mean)Group.1 x1       0 32.308502       5 29.308423      50 29.866114     500 29.64647#获得各组均值，可以发现未用药组幼崽体重均值最高(32.3)fit 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

ANCOVA的F检验表明：(a)怀孕时间与幼崽出生体重相关；(b)控制怀孕时间，药物剂量与出生体重相关。

可使用effects包中的effects()函数获取调整的组均值，即去除协变量效应后的组均值：

library(effects)effect("dose", fit)dose effectdose0        5       50      50032.35367 28.87672 30.56614 29.33460

同样可使用multcomp包来对所有均值进行成对比较，还可以用来检验用户自定义的均值假设：

library(multcomp)contrast |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)

假设检验的t统计量(2.581) 在p<0.05水平下显著，可以得出未用药组比其他用药条件下的出生体重高的结论。

1. 评估检验的假设条件

ANCOVA与ANOVA相同，都需要正态性和同方差性假设，可用上文提到的Q-Q图、Bartlett检验和离群点检验。

ANCOVA还假定回归斜率相同。

例：multcomp包中的litter数据集假定四个处理组通过怀孕时间来预测出生体重的回归斜率都相同。ANCOVA模型包含怀孕时间×剂量的交互项时，可对回归斜率的同质性进行检验，交互效应若显著，则意味着时间和幼崽出生体重间的关系依赖于药物剂量的水平。

检验回归斜率的同质性：

library(multcomp)fit2 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.17889Residuals     66 1069.4   16.20---Signif. codes:0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

可以看到交互效应不显著，支持了斜率相等的假设。若假设不成立，可以尝试变换协变量或因变量，或使用能对每个斜率独立解释的模型，或使用不需要假设回归斜率同质性的非参数 ANCOVA方法。

2. 结果可视化

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

library(HH)ancova(weight ~ gesttime + dose, data=litter)Analysis of Variance TableResponse: weightDf  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

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

用怀孕时间来预测出生体重的回归线相互平行，只是截距项不同。

若用ancova(weight ~ gesttime*dose)，生成的图形将允许斜率和截距项依据组别而发生变化，可以可视化违背回归斜率同质性的实例。

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

图9-6 四种药物处理组的怀孕时间和出生体重的关系图(斜率和截距项可改变)

9.4 双因素方差分析

在双因素方差分析中，受试者被分配到两因子的交叉类别组中。

例：基础安装中的ToothGrowth数据集随机分配60只豚鼠，分别采用两种喂食方法(橙汁或维生素C)，各喂食方法中抗坏血酸含量有三种水平(0.5mg、1mg或2mg)，每种处理方式组合都被分配10只豚鼠。

attach(ToothGrowth)table(supp, dose)dosesupp 0.5  1  2OJ  10 10 10VC  10 10 10#获得各单元的均值aggregate(len, by=list(supp, dose), FUN=mean)Group.1 Group.2     x1      OJ     0.5 13.232      VC     0.5  7.983      OJ     1.0 22.704      VC     1.0 16.775      OJ     2.0 26.066      VC     2.0 26.14#获得各单元的标准差aggregate(len, by=list(supp, dose), FUN=sd)Group.1 Group.2        x1      OJ     0.5 4.4597092      VC     0.5 2.7466343      OJ     1.0 3.9109534      VC     1.0 2.5153095      OJ     2.0 2.6550586      VC     2.0 4.797731#dose变量被转换为因子变量dose F)supp         1  205.4   205.4  15.572 0.000231 ***dose         2 2426.4  1213.2  92.000

用summary()函数得到方差分析表，可以看到主效应(supp和dose)和交互效应都非常显著。

结果可视化

用interaction.plot()函数展示双因素方差分析的交互效应：

attach(ToothGrowth)interaction.plot(dose, supp, len, type="b",col=c("red","blue"), pch=c(16, 18),main = "Interaction between Dose and Supplement Type")detach(ToothGrowth)

图9-7 喂食方法和剂量对牙齿生长的交互作用

用interaction.plot()函数绘制了牙齿长度的均值。

用gplots包中的plotmeans()函数展示交互效应：

library(gplots)attach(ToothGrowth)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")detach(ToothGrowth)

图9-8 喂食方法和剂量对牙齿生长的交互作用

用plotmeans()函数绘制的95%的置信区间的牙齿长度均值

用HH包中的interaction2wt()函数来可视化结果：

library(HH)attach(ToothGrowth)interaction2wt(len ~ supp * dose)detach(ToothGrowth)

图9-9 ToothGrowth数据集的主效应和交互效应

推荐HH包中的interaction2wt()函数，能展示任意复杂度设计(双因素方差分析、三因素方差分析等)的主效应(箱线图)和交互效应。

9.5 重复测量方差分析

重复测量方差分析(单因素组内方差分析)，即受试者被测量不止一次。

例：基础安装包中的CO2数据集包含了北方和南方牧草类植物Echinochloa crus-galli的寒冷容忍度研究结果。比较某浓度二氧化碳的环境中寒带植物与非寒带植物的光合作用率。自变量是植物类型Type(组间因子)和七种水平的二氧化碳浓度conc(组内因子)。

含一个组间因子和一个组内因子的重复测量方差分析：

CO2$conc F)Type       1 2667.2  2667.2   60.41 0.00148 **Residuals  4  176.6    44.1---Signif. codes:0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Error:Plant:concDf Sum Sq Mean Sq F value   Pr(>F)conc       6 1472.4  245.40   52.52 1.26e-12 ***conc:Type 6 428.8 71.47 15.30 3.75e-7 ***Residuals 24 112.1 4.67---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

方差分析表表明在0.01的水平下，主效应类型和浓度以及交叉效应类型×浓度都非常显著。

通过interaction.plot()函数展示了交互效应：

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

图9-10 interaction.plot()函数绘制CO2浓度和植物类型对CO2吸收的交互影响

可以使用boxplot()函数对相同的数据画图，展示交互效应其他不同的侧面：

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-11 boxplot()函数绘制CO2浓度和植物类型对CO2吸收的交互效应

魁北克省的植物比密西西比州的植物二氧化碳吸收率高，而且随着CO2浓度的升高，差异越来越明显。

9.6 多元方差分析

当因变量(结果变量)不止一个时，可用多元方差分析(MANOVA)对它们同时进行分析。

例：MASS包中的UScereal数据集。

研究美国谷物中的卡路里、脂肪和糖含量是否会因为储存架位置的不同而发生变化。

单因素多元方差分析：

library(MASS)attach(UScereal)#将shelf变量转换为因子变量，从而使它在后续分析中能作为分组变量shelf 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

F值显著，说明三个组的营养成分测量值不同。

#summary.aov()函数对每一个变量做单因素方差分析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 ‘ ’ 1Response 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 ‘ ’ 1Response 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

三组中每种营养成分的测量值都是不同的。

1. 评估假设检验

单因素多元方差分析有两个前提假设：

多元正态性，指因变量组合成的向量服从一个多元正态分布，可用Q-Q图来检验。

方差-协方差矩阵同质性，R中还没有好的方法检验。

检验多元正态性：

center

图9-12 检验多元正态性的Q-Q图

观测点“Wheaties Honey Gold”和“Wheaties”异常，数据集似乎违反了多元正态性，可以删除这两个点再重新分析。

可以使用mvoutlier包中的ap.plot()函数来检验多元离群点：

library(mvoutlier)outliers

图9-13 检验多元离群点

2. 稳健多元方差分析

如果多元正态性或者方差-协方差均值假设都不满足，或者担心多元离群点，那么可以考虑用稳健或非参数版本的MANOVA检验。

稳健单因素MANOVA可通过rrcov包中的Wilks.test()函数实现：

library(rrcov)Wilks.test(y, shelf, method="mcd")Robust One-way MANOVA (Bartlett Chi2)data:xWilks' Lambda=0.511, Chi2-Value = 23.96, DF = 4.98, p-value = 0.0002167sample estimates:calories    fat  sugars1 120 0.701 5.662      128  1.185   12.543 161 1.652 10.35

稳健检验对离群点和违反MANOVA假设的情况不敏感，验证了存储在货架顶部、中部和底部的谷物营养成分含量不同。

9.7 用回归来做ANOVA

例：multcomp包中的cholesterol数据集比较五种降低胆固醇药物疗法(trt)的影响：

library(multcomp)levels(cholesterol$trt)[1] "1time" "2times" "4times" "drugD" "drugE"fit.aov F)trt          4  1351.37   337.84   32.433  9.819e-13 ***Residuals   45   468.75    10.42---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

用lm()函数拟合同样的模型：

fit.lm |t|)(Intercept) 5.782 1.021 5.665 9.78e-07 ***trt2times      3.443      1.443   2.385   0.0213 *    #变量trt2times表示水平1time和2times的一个对照trt4times      6.593      1.443   4.568 3.82e-05 ***  #trt4times是1time和4times的一个对照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 ‘ ’ 1Residual standard error: 3.227 on 45 degrees of freedomMultiple R-squared:0.7425,  Adjusted R-squared:  0.7196F-statistic:32.43 on 4 and 45 DF,  p-value: 9.819e-13

从输出的概率值来看，各药物条件与第一组(1time)显著不同。

因为线性模型要求预测变量是数值型，lm()函数会用一系列与因子水平相对应的数值型对照变量来代替因子。默认情况下，对照处理用于无序因子，正交多项式用于有序因子。

表9-5 R提供的五种创建对照变量的内置方法

对照变量创建方法

描述

contr.helmert

第二个水平对照第一个水平，第三个水平对照前两个的均值，第四个水平对照前三个的均值，以此类推

contr.poly

基于正交多项式的对照，用于趋势分析(线性、二次、三次等)和等距水平的有序因子

contr.sum

对照变量之和限制为0。也称作偏差找对，对各水平的均值与所有水平的均值进行比较

contr.treatment

各水平对照基线水平(默认第一个水平)。也称作虚拟编码

contr.SAS

类似于contr.treatment，只是基线水平变成了最后一个水平。生成的系数类似于大部分SAS过程中使用的对照变量

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

contrasts(cholesterol$trt)2times 4times drugD drugE1time 0 0 0 02times      1      0     0     04times 0 1 0 0drugD 0 0 1 0drugE 0 0 0 1

若患者处于drugD条件下，变量drugD等于1，其他变量2times、4times和drugE都等于0；无需列出第一组的变量值，因为其他四个变量都为0，这已经说明患者处于1time条件。

通过设定contrasts选项，可以修改lm()中默认的对照方法：

fit.lm

还能通过options()函数修改R会话中的默认对照方法：

options(contrasts = c("contr.SAS", "contr.helmert"))

设定无序因子的默认对比方法为contr.SAS，有序因子的默认对比方法为contr.helmert。

再查看它的编码过程：

contrasts(cholesterol$trt)1time 2times 4times drugD1time 1 0 0 02times 0 1 0 04times 0 0 1 0drugD 0 0 0 1drugE 0 0 0 0

9.8 小结

本章我们通过组内和组间设计的示例介绍了方差分析模型以及假设检验。学习完回归和方差分析这两种常用的统计方法后，下一章我们将学习对研究设计非常重要的功效分析。敬请期待～