R语言学习系列28-协方差分析

23. 协方差分析

一、基本原理

1. 基本思想

在实际问题中，有些随机因素是很难人为控制的，但它们又会对结果产生显著影响。如果忽略这些因素的影响，则有可能得到不正确的结论。这种影响的变量称为协变量（一般是连续变量）。

例如，研究3种不同的教学方法的教学效果的好坏。检查教学效果是通过学生的考试成绩来反映的，而学生现在考试成绩是受到他们自身知识基础的影响，在考察的时候必须排除这种影响。

协方差分析将那些难以控制的随机变量作为协变量，在分析中将其排除，然后再分析控制变量对于观察变量的影响，从而实现对控制变量效果的准确评价。

协方差分析要求协变量应是连续数值型，多个协变量间互相独立，且与控制变量之间没有交互影响。前面单因素方差分析和多因素方差分析中的控制变量都是一些定性变量，而协方差分析中既包含了定性变量（控制变量），又包含了定量变量（协变量）。

协方差分析在扣除协变量的影响后再对修正后的主效应进行方差分析，是一种把直线回归或多元线性回归与方差分析结合起来的方法，其中的协变量一般是连续性变量，并假设协变量与因变量间存在线性关系，且这种线性关系在各组一致，即各组协变量与因变量所建立的回归直线基本平行。

当有一个协变量时，称为一元协方差分析，当有两个或两个以上

的协变量时，称为多元协方差分析。

2. 协方差分析需要满足的条件

（1）自变量是分类变量，协变量是定距变量，因变量是连续变量；对连续变量或定距变量的协变量的测量不能有误差；

（2）协变量与因变量之间的关系是线性关系，可以用协变量和因变量的散点图来检验是否违背这一假设；协变量的回归系数（即各回归线的斜率）是相同的，且不等于0，即各组的回归线是非水平的平行线。否则，就有可能犯第一类错误，即错误地接受虚无假设；

（3） 自变量与协变量相互独立，若协方差受自变量的影响，那么协方差分析在检验自变量的效应之前对因变量所作的控制调整将是偏倚的，自变量对因变量的间接效应就会被排除；

（4）各样本来自具有相同方差σ2的正态分布总体，即要求各组方差齐性。

二、协方差理论

1. 观测值=均值+分组变量影响+协变量影响+随机误差. 即

()ij i ij ij y u t x x βε=++-+ （1） 其中，X 为所有协变量的平均值。 注：在方差分析中，协变量影响是包含在随机误差中的，在协方差分析中需要分离出来。

用协变量进行修正，得到修正后的y ij (adj)为

(adj)()ij ij ij i ij y y x x u t βε=--=++

就可以对y ij (adj)做方差分析了。关键问题是求出回归系数β.

2. 总离差=分组变量离差+协变量离差+随机误差，

（1）计算总离差平方和时，记

11

()()k n

xy ij ij i j T x x y y ===--∑∑

211

()k n

xx ij i j T x x ===-∑∑ 总离差平方和：211()k n

yy ij i j T y y ===-∑∑

最终要检验分组自变量对因变量有无显著作用。原假设H 0：无显著作用。假设检验是在H 0为真条件下进行，可认为t i =0，则

()ij T ij ij y u x x βε=+-+

按最小二乘法原理线性回归可得到β的估计值

?xy T xx

T T β= 记修正的总离差平方和（残差平方和）为T yy(adj)，则

22(adj)?T xy yy yy xx yy xx T T T T T T β=-=-，自由度为n -2

注：2?T xx T β为回归平方和，若?0T

β=（回归线为水平线），表示协变量x 对y 无作用，用方差分析就可以解决了。

（2）计算组内离差平方和时，记

11

()()k n

xy ij i ij i i j E x x y y ===--∑∑

211()k n

xx ij i i j E x x ===-∑∑

组内总离差平方和：211()k n

yy ij i i j E y y ===-∑∑

根据协方差分析的基本假设：各组内回归系数相等（做协方差分析时需要检验这一点），得到组内回归系数βw 的估计值

?xy w xx

E E β= 记修正的组内总离差平方和（组内残差平方和）为E yy(adj), 则

22(adj)?xy yy yy w xx yy xx E E E E E E β=-=-， 自由度为n -k -1

其中，2?w xx E β为组内回归平方和，当1??w wk

ββ==L 时，组内总离差平方和认为完全是由随机因素引起的，E yy(adj)就是随机为误差。这里的?w β是1??,,w wk

ββL 的加权平均值。 （3）计算分组变量离差平方和B yy(adj)，它反映的是各个水平之间的差异。

2(adj)(adj)(adj)(adj)?T yy yy yy yy xx yy B T E T T E β=-=--

即，分组变量离差=总离差-协变量离差-随机误差。

于是，就可以进行组间无差异检验了：

(adj)(adj)/1

/1yy yy B k F E n k -=--

3. 因此，在做协方差分析前，需要依次做两个假设检验：

（1）协变量对因变量的影响对与各组来说都是相同的，即各组

回归系数相等：1???:w wk w

βββ===L ; 步骤：

① 先按回归系数相等和不相等分别表示模型

()ij i w ij ij y u t x x βε=++-+

()ij i wi ij ij y u t x x βε=++-+

并计算出误差平方和

2(adj)yy yy w xx E E E β=-

211

i k

yy wi xx i S E E β==-∑ 其中，1

i k yy yy i E E ==∑. ② 计算F 值

(adj)11/1

/2yy E S k F S n k --=-

若F 值小于临界值F α，则说明各组回归系数无显著差异（相等）。

（2）这些相等的回归系数?0w

β≠. 即采用一元线性回归的显著性检验，

2(adj)/1=//(1)

w xx yy E F E n k β=--回归平方和/自由度残差平方和自由度 222

2

/(1)(/)/(1)xy xx

xy yy xy xx yy xx xy E E E n k E E E n k E E E --==----

4. 协方差分析的步骤

（1）检验数据是否满足假设条件：正态分布性、方差齐性、各分组通过协变量预测因变量的回归斜率相同；

（2）检验效应因子的显著性；

（3）估计校正的组均值；

（4）检验校正的组均值之间的差异。

三、R语言实现

协方差分析要求数据满足：正态性、方差齐性、各分组通过协变量预测因变量的回归斜率相同。

R语言用aov()函数进行协方差分析，基本格式为：

aov(formula, data, ...)

其中，data为数据框；

formula为协方差公式形式，形如y~x+A, x为连续型协变量，A 为组别因子。

例1研究分别接受了3种不同的教学方法的3组学生，在数学成绩上是否有显著差异，数据文件“ex28_cov.Rdata”。

先不考虑数学入学成绩，只以“教学方法”为分组变量，“后测成绩”为因变量进行单因素方差分析：

setwd("E:/办公资料/R语言/R语言学习系列/codes")

load("ex28_cov.Rdata")

head(scores)

before after teach

1 39 68 1

2 38 6

3 1

3 51 65 1

4 56 68 1

5 74 74 1

6 40 60 1

attach(scores)

table(teach) #各组的样本数

teach

1 2 3

30 32 33

aggregate(after, by=list(teach), mean) #各组均值 Group.1 x

1 1 62.88333

2 2 72.67188

3 3 65.06061

shapiro.test(after) #正态性检验

Shapiro-Wilk normality test

data: after

W = 0.99105, p-value = 0.7772

bartlett.test(after~teach,data=scores) #方差齐性检验Bartlett test of homogeneity of variances

data: after by teach

Bartlett's K-squared = 0.69854, df = 2, p-value = 0.7052

fit.aov<-aov(after~teach,data=scores)

summary(fit.aov)

Df Sum Sq Mean Sq F value Pr(>F)

teach 2 1662 830.8 10.44 8.23e-05 ***

Residuals 92 7325 79.6

---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

说明：单因素方差分析的p值=8.23e-05, 远小于0.05, 表明，两种教学方法有非常显著的差异。

但是，后测成绩肯定会受到前测成绩（连续型）的影响，假定前测成绩与教学方法（即组别，是控制变量）不存在交互影响。因此，将后测成绩作为因变量；教学方法作为控制变量；前测成绩作为协变量进行协方差分析。

回归斜率相同检验，即前测成绩与后测成绩的回归线是否平行：

scores1<-subset(scores,teach==1)

scores2<-subset(scores,teach==2)

scores3<-subset(scores,teach==3)

par(mfrow=c(1,3))

plot(scores1$before,scores1$after,xlab="before",ylab= "after",main="teach=1")

abline(lm(after~before,data=scores1))

plot(scores2$before,scores2$after,xlab="before",ylab= "after",main="teach=2")

abline(lm(after~before,data=scores2))

plot(scores3$before,scores3$after,xlab="before",ylab=" after",main="teach=3")

abline(lm(after~before,data=scores3))

可见两组的直线趋势的斜率比较接近（平行），基本符合协方差假定。

除了图形判断外，还可以通过交互作用是否显著，来判断斜率是否相同。因为若前验成绩与教学方法的交互作用显著，则说明前验成绩与后验成绩的关系，依赖于教学方法。

library(multcomp)

fit2<-aov(after~before*teach,data=scores)

summary(fit2)

Df Sum Sq Mean Sq F value Pr(>F) before 1 2432 2432.0 35.391 5.22e-08 ***

teach 2 362 180.9 2.633 0.0775 .

before:teach 2 76 38.2 0.556 0.5752

Residuals 89 6116 68.7

---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1说明：交互项的p值=0.5752>0.05, 故不显著，支持了斜率相同的假设。

fit<-aov(after~before+teach,data=scores) #协方差分析

summary(fit)

Df Sum Sq Mean Sq F value Pr(>F)

before 1 2432 2432.0 35.739 4.35e-08 ***

teach 2 362 180.9 2.659 0.0755 .

Residuals 91 6192 68.0

---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1说明：协方差分析的结果表明：前测成绩的p值=4.35e-08远小于0.05, 说明“前测成绩”对“后测成绩产生了非常显著的影响；“教学方法”的p值=0.0755>0.05, 说明“教学方法”对“后测成绩”的影响不显著。

由于受到协变量的影响，我们希望获取调整后的各组均值——即去除协变量效应后的各组均值。可使用effects包中的effects()函数来计算调整的均值：

library(effects)

effect("teach", fit)

teach effect

teach

1 2 3

65.05722 70.04958 65.62718

与方差分析时一样，要想得到教学方法两两之间有无差异，可以均值的成对比较（略）。

下面讲一下自定义比较（使用multcomp包可以实现），例如，分组变量有4个水平ABCD，要比较A与D时，比较矩阵=[1 0 0 -1]T, 有

[A B C D]×[1 0 0 -1]T=0 等价于A=D

要想将A与D合并再与B比较有无差异，则可以指定L矩阵=[1 -2 0

1]T, 则

[A B C D]×[1 -2 0 1]T =0 等价于(A+D)/2 = B

注意：是从(A+D)/2 = B倒推比较矩阵，该式即A-2B+0C+D=0.

根据调整后的各组均值，教学方法1和3基本相同，虽然总体上3种差异不显著，教学方法2与1、3是否有显著差异呢？那么就需要自定义比较。

library(multcomp)

contrast<-rbind("2 vs. 13"=c(-1,2,-1))

res.vs<-glht(fit,linfct=mcp(teach=contrast))

summary(res.vs)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts

Fit: aov(formula = after ~ before + teach, data = scores)

Linear Hypotheses:

Estimate Std. Error t value Pr(>|t|)

2 vs. 1

3 == 0 9.415 4.083 2.306 0.023

4 *

---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Adjusted p values reported -- single-step method)

说明：自定义比较：教学方法2与教学方法1,3是否有差异，设置比较矩阵为[-1,2,-1], 结果p值=0.0234<0.05, 拒绝原假设，即有显著差异。

另外，HH包中的ancova()函数，也是用来做协方差分析的，还能将结果可视化。基本格式为：

ancova(formula, data=, x, groups, ...)

其中，若formula不包括，x和groups为作图时需指明协变量和因子。

library(HH)

ancova(after~before+teach, data=scores)

Analysis of Variance Table

Response: after

Df Sum Sq Mean Sq F value Pr(>F)

before 1 2432.0 2431.96 35.7392 4.354e-08 ***

teach 2 361.9 180.93 2.6589 0.07546 .

Residuals 91 6192.3 68.05

---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

从图中可以看出，用前验成绩来预测后验成绩的回归线相互平行，只是截距项不同。教学方法13基本相同，教学方法2明显好于13。

上述代码会让直线保持平行，若用

ancova(after~before*teach, data=scores)

则生成的图形将允许斜率和截距项依据组别而发生变化，这对可视化那些违背回归斜率同质性的实例非常有用。