方差分析主要通过F检验来进行效果评测，若治疗方案的F检验显著，则说明检验样本组间均值不同。

ANOVA模型拟合

从函数形式上看，ANOVA和回归方法都是广义线性模型的特例。因此回归分析章节中提到的lm()函数也能分析ANOVA模型。不过，在这个章节中，我们基本使用aov()函数。最后，会提供了个lm()函数的例子。

aov()函数

aov()函数的语法为aov(formula, data=dataframe)。下表列举了表达式可以使用的特殊符号。

符号

用途

~

分隔符号，左边为响应变量(因变量)，右边为解释变量(自变量)

:

表示预测变量的交互项

*

表示所有可能交互项的简洁方式

^

表示交互项达到某个次数

.

表示包含除因变量外的所有变量

下面是常见研究设计的表达式

设计

表达式

单因素ANOVA

y ~ A

含单个协变量的单因素ANOVA

y ~ x + A

双因素ANOVA

y ~ A * B

含两个协变量的双因素ANOVA

y ~ x1 + x2 + A * B

随机化区组

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

单因素组内ANOVA

y ~ A + Error(subject/A)

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

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

表达式中各项的顺序

当

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

​ 存在协变量

两者之一时，等式右边的变量都与其他变量相关。此时，我们无法清晰地划分它们对因变量的影响。

例如，对于双因素方差分析，若不同处理方式中的观测数不同，那么模型y ~ A*B与模型y ~ B*A的结果不同。

R默认类型I(序贯型)方法计算ANOVA效应(类型II和III分别为分层和边界型，详见R实战(第2版)202页)。R中的ANOVA表的结果将评价：

A对y的影响

控制A时，B对y的影响

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

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

car包的Anova()函数提供了三种类型方法，若想与其他软件(如SAS SPSS)提供的结果保持一致，可以使用它，细节可参考helo(Anova, package="car")。

单因素方差分析

单因素方法分析中，你感兴趣的是比较分类因子定义的两个或多个组别中的因变量均值。以multcomp包中cholesterol数据集为例(包含50个患者接收5种降低胆固醇疗法的一种，前三种是同样的药物不同的用法，后二者是候选药物)。哪种药物疗法降低胆固醇最多呢？

> library(mvtnorm)

> library(survival)

> library(TH.data)

> library(MASS)

> library(multcomp)

> 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

> 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

> library(gplots)

载入程辑包：‘gplots’

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

lowess

> plotmeans(response ~ trt, xlab = "Treatment", ylab = "Response", main = "Mean Plot \nwith 95% CI")

> detach(cholesterol)

mean_plot.png

从结果可以看到，均值显示drugE降低胆固醇最多，各组标准差相对恒定。ANOVA对治疗方式的F检验非常显著，说明五种疗法的效果不同。

多重比较

虽然ANOVA对各种疗法的F检验表明五种药物的治疗效果不同，但是没有告诉你哪种疗法与其他疗法不同。多重比较可以解决这个问题。例如，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(mar=c(5,8,4,2))

> plot(TukeyHSD(fit))

第一个par语句用来旋转轴标签，第二个用来增大左边界的面积，可使标签摆放更美观。

成对比较图形如下图所示。

conf_level.png

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

> library(multcomp)

> par(mar=c(5,4,6,2))

> tur

> plot(cld(tur, level=0.05), col="lightgrey")

glht.png

par语句增大了顶部边界面积，cld()函数中的level选项设置了使用的显著水平。

有相同的字母的组说明均值差异不显著。

评估检验的假设条件

可以使用Q-Q图来检验正态性

> library(car)

> qqPlot(lm(response ~ trt, data = cholesterol), simulate=T, main="Q-Q Plot", labels=FALSE)

qqplot.png

注意qqPlot需要lm()拟合。

方差齐次性检验：

例如，可以通过如下代码做Bartlett检验

> 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

结果表明五组方差没有显著地不同。

注意，方差齐性分析对离群点非常敏感。可以利用car包outlierTest()检验。

单因素协方差分析

ANCOVA扩展了ANOVA，包含一个或多个定量的协变量。

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

> 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

> 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

由于使用了协变量，我们可以通过effect包中的effects()函数计算调整的组均值。

> library(effects)

载入程辑包：‘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

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

下面代码清单可以用来检验未用药和其他三种药条件影响是否不同。

> library(multcomp)

> contrast

> 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在p<0.05下显著。由此可以得出结论。详见help(glht)。

评估检验的假设条件

ANCOVA与ANOVA相同，都需要正态性和同方差假设，检验可以参考上一节。另外ANCOVA还假定回归斜率相同。

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

fit2

summary(fit2)

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

library(HH)

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

用回归来做ANOVA

同样是之前比较五种降低胆固醇药物疗法的影响的例子，我们分别用两种不同的方法来做(aov()和lm())。

> library(multcomp)

载入需要的程辑包：mvtnorm

载入需要的程辑包：survival

载入需要的程辑包：TH.data

载入需要的程辑包：MASS

载入程辑包：‘TH.data’

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

geyser

> levels(cholesterol)

NULL

> levels(cholesterol$trt)

[1] "1time" "2times" "4times" "drugD" "drugE"

> fit.aov

> 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

> 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