SAS讲义 第二十六课协方差分析

第二十六课 协方差分析
当定量的影响因素对观察结果有难以控制的影响，甚至还有交互作用时，采用协方差分析，这些影响变量称为协变量，扣除（或消除）协变量的影响，可以得到修正后的均值估计。

一、 协方差分析概述
1. 协方差分析概念
协方差分析（analysis of covariance ）又称带有协变量的方差分析（analysis of variance with covariates ），是将回归分析与方差分析结合起来使用的一种分析方法。

在各种试验设计中，对主要变量y 研究时，常常希望其他可能影响和干扰y 的变量保持一致以到达均衡或可比，使试验误差的估计降到最低限度，从而可以准确地获得处理因素的试验效应。

但是有时，这些变量难以控制，或者根本不能控制。

为此需要在试验中同时记录这些变量的值，把这些变量看作自变量，或称协变量（covariate ），建立因变量y 随协变量变化的回归方程，这样就可以利用回归分析把因变量y 中受协变量影响的因素扣除掉，从而，能够较合理地比较定性的影响因素处在不同水平下，经过回归分析手段修正以后的因变量的总体均值之间是否有显著性的差别。

简单地说，协方差分析是扣除协变量的影响，或者将这些协变量处理成相等，再对修正的y 的均值作方差分析。

2. 协方差分析的假定
协方差分析需要满足的假定为：
①各样本来自具有相同方差2
σ的正态分布总体，即要求各组方差齐性。

②协变量与主要变量y 间的总体回归系数不等于0。

③各组的回归线平等，即回归系数 ==21ββ。

如果上述的假定满足，就作协方差分析。

前述的各种试验设计，如完全随机化设计、随机区组设计、析因设计、拉丁方设计等，都可以带一个或多个协变量，按设计方案扣除协变量的影响后，对主要变量y 的修正均值作比较，得出统计结论。

3. 协方差分析的模型
最简单的单因素一元协方差分析的模型，是由单因素效应模型ij i ij a y εμ++=加上协变量的影响因素)(x x ij -β而得出：
ij ij i ij x x a y εβμ+-++=)(
(26.1)
其中x 为协变量，ij x 为协变量在分类水平i 和j 上的记录值，x 为所有协变量的平均值，β为相关的回归系数。

设x βμβ-=0，为平均截距。

上式可以化简成
ij ij i ij x a y εββ+++=0
(26.2)
设i i a +=00ββ，上式可以化简成
ij ij i ij x y εββ++=0
(26.3)
很明显i 0β是第i 组回归线的截距，等于回归线的平均截距0β加上本组的效应i a 。

这个式揭示了，观察值ij y 的模型可以表示成一组相似的回归线，且各组具有共同的回归系数β，和各组自己的截距i i a +=00ββ。

用SAS 中的glm 过程进行协方差分析时，要注意不同试验设计时class 语句和model 语句的写法。

设分类变量为A 、B ，协变量为X ，观察值为Y ，则有：
①单因素k 水平设计的协方差分析模型
class A; model X A ;
②随机区组设计的协方差分析模型
class A B; model X A B ;
③两因素析因设计的协方差分析模型
class A B;
model X A B A*B;
二、 实例分析
1. 一元协方差分析
例26.1 研究牡蛎在不同温度的水中不同位置上的生长情况。

有人做了如下试验：分别在通向发电站的入口处（温度较低）不同位置（底部和表层）和出口处（温度较高）不同位置（底部和表层）及电站附近的深水处（底部和表层的中间）总共5个不同位置点上，随机地各放4袋牡蛎（每袋中有10个），共5×4=20袋。

在将每袋牡蛎放入位置点之前，先洗干净称出每袋的初始体重，放在5个不同点一个月后再称出最后体重。

试验结果数据见表26.1所示。

表26.1 牡蛎在不同温度和位置上的生长数据
位置 trt 重复数rep （x 为初始体重，y 为最后体重） 1 2
3
4
x y x y x y x y 1（入口底部） 27.2 32.6 32.0 36.6 33.0 37.7 26.8 31.0 2（入口顶部） 28.6 33.8 26.8 31.7 26.5 30.7 26.8 30.4 3（出口底部） 28.6 35.2 22.4 29.1 23.2 28.9 24.4 30.2 4（出口顶部） 29.3 35.0 21.8 27.0 30.3 36.4 24.3 30.5 5（附近中部） 20.4
24.6
19.6
23.4
25.1
30.3
18.1
21.8
程序如下：
data growth;
do trt=1 to 5; do rep=1 to 4;
input x y @@; output;
end;
end;
cards;
27.2 32.6 32.0 36.6 33.0 37.7 26.8 31.0
28.6 33.8 26.8 31.7 26.5 30.7 26.8 30.4
28.6 35.2 22.4 29.1 23.2 28.9 24.4 30.2
29.3 35.0 21.8 27.0 30.3 36.4 24.3 30.5
20.4 24.6 19.6 23.4 25.1 30.3 18.1 21.8
;
proc anova data=growth;
class trt;
model y=trt;
proc glm data=growth;
class trt;
model y=trt x /solution;
means trt;
lsmeans trt /stderr tdiff;
contrast 'trt12 vs trt34' trt -1 -1 1 1 0;
estimate 'trt1 adj mean' intercept 1 trt 1 0 0 0 0 x 25.76;
estimate 'trt2 adj mean' intercept 1 trt 0 1 0 0 0 x 25.76;
estimate 'adj trt diff' trt 1 -1 0 0 0;
estimate 'trt1 unadj mean' intercept 1 trt 1 0 0 0 0 x 29.75;
estimate 'trt2 unadj mean' intercept 1 trt 0 1 0 0 0 x 27.175;
estimate 'unadj trt diff' trt 1 -1 0 0 0 x 2.575;
run;
程序说明：定性变量trt的5个不同位置点对y可能有较大的影响，因此class语句中分组变量为trt，先选用anova过程进行方差分析。

然而，牡蛎的初始体重x对牡蛎的最后体重y可能也有一定的影响，故适合选用glm过程进行协方差分析，在model语句中不仅包括分组变量trt，而且应包括协变量x。

选择项solution要求输出回归系数的估计值及其标准误差和假设检验等结果。

means和lsmeans语句要求输出分组变量trt各水平下y的未修正均值和修正后的均值，选择项stderr要求输出y的修正均值的标准误差、各修正均值与0比较的假设检验结果；选择项tdiff要求输出y的各修正均值之间两两比较所对应的t值和p值。

Contrast语句是用来比较入口处底部和顶部均值之和与出口处底部和顶部均值之和是否相等。

前三条estimate语句是用来估计入口处底部和顶部调整后的均值及它们之差，并假设检验是否为0，后三条estimate语句是用来估计入口处底部和顶部未调整的均值及它们之差，并假设检验是否为0。

程序输出的主要结果见表26.2（a）（b）（c）所示。

表26.2（a）单因素trt一元x的协方差分析
表26.2(a)中结果分析：对分组变量trt的方差分析表明，即使当初始体重x不考虑，各分
The SAS System
Analysis of Variance Procedure
Dependent Variable: Y
Source DF Sum of Squares Mean Square F Value Pr > F Model 4 198.40700000 49.60175000 4.64 0.0122 Error 15 160.26250000 10.68416667
Corrected Total 19 358.66950000
R-Square C.V. Root MSE Y Mean
0.553175 10.59706 3.26866436 30.84500000
Source DF Anova SS Mean Square F Value Pr > F TRT 4 198.40700000 49.60175000 4.64 0.0122
General Linear Models Procedure
Dependent Variable: Y
Source DF Sum of Squares Mean Square F Value Pr > F Model 5 354.44717675 70.88943535 235.05 0.0001 Error 14 4.22232325 0.30159452
Corrected Total 19 358.66950000
R-Square C.V. Root MSE Y Mean
0.988228 1.780438 0.54917622 30.84500000
Source DF Type I SS Mean Square F Value Pr > F TRT 4 198.40700000 49.60175000 164.47 0.0001 X 1 156.04017675 156.04017675 517.38 0.0001 Source DF Type III SS Mean Square F Value Pr > F TRT 4 12.08935928 3.02233982 10.02 0.0005 X 1 156.04017675 156.04017675 517.38 0.0001
T for H0: Pr > |T| Std Error of Parameter Estimate Parameter=0 Estimate
INTERCEPT 2.494859769 B 2.43 0.0293 1.02786287
组最后体重均值的区别也统计显著（0.0122<0.05），其中分组变量trt的平方和为198.40700000。

而在协方差分析中，分组变量trt的类型1的平方和等于方差分析中的平方和198.40700000，分组变量trt的类型3的平方和为12.08935928，大大小于类型1的平方和，是因为类型3的平方和反映了经过共同的协变量x调整后的平方和，减去了协变量的影响，所以平方和大幅减小。

类型1是一种未经过调整的平方和，因为它的优先级高于协变量的调整。

更进一步分析，我们注意到方差分析中均方误差为10.68416667，而协方差分析中却缩小到0.30159452，相应地分组变量trt的F统计量从4.64增加到10.02，说明包含了协变量后分组的区别更加显著，原因是简单方差分析中，大多数的误差是由于初始体重x的变异造成的。

表中的最后一部分是选择项solution的输出结果，对模型中的截距、各分组变量和协变量的回归系数进行估计和检验，在这个单因素trt的情况下，估计是以最后一个水平trt5（trt=5）为对照组，并且设置它的系数为0，因此截距intercept的估计值是分组trt5的估计值。

其他四个分组trt的系数估计是每一个与trt5进行比较而得到的。

注意，出口处的trt3和trt4分组是不同与trt5分组。

协变量x的系数是合并各组内y和x所得到的回归系数，即是由5个独立
的trt 分组，分别回归y 和x 后得到回归系数然后加权平均。

协变量x 的系数估计值表明，初始体重变动1个单位最后体重y 相关地要变动1.083179819单位。

表26.2（b ） 未调整均值和调整均值及均值之间的比较
表26.2（b ）中结果分析：means 语句要求计算按trt 每个水平分组的未调整的y 和x 均值。

如•1y =34.475=(32.6+36.6+37.7+31)/4，•1x =29.75=(27.2+32+33+26.8)/4。

Lsmeans 语句要求计算调整后的y 的均值，或称最小二乘均值估计，我们可以由（26.1）公式求分组平均得到：
)
()(x x y y x x y a i i i ij
ij ij i --=+--=+•••βεβμ
(26.4)
再由（25.2）公式求分组平均代入上式：
x
a x x y x x y y i i i i i i βββββ
++=+-=--=•••••0)
( (26.5)
例如，初始体重的整体平均值为x =（29.750+27.175+24.650+26.425+20.800）/5=25.76，以trt 1分组为例，调整后•1y =30.1531125=34.475－1.083179819×（29.75－25.76）。

tdiff 选择项要求对已调整均值的两两比较采用lsd 检验，可以使用adjust = duncan/waller 等选项替代lsd 检验，获得其他多重比较的检验结果。

从最后的5×5修正均值比较结果表中，可得到（•••521,,y y y ）中的任何一个与（••43,y y ）中的任何一个之间有显著或非常显著性差别。

表26.2（c ） 有计划的均值对比和参数估计
The SAS System
General Linear Models Procedure
Level of --------------Y-------------- --------------X-------------- TRT N Mean SD Mean SD 1 4 34.4750000 3.18891309 29.7500000 3.20572405 2 4 31.6500000 1.53731367 27.1750000 0.96046864 3 4 30.8500000 2.95578529 24.6500000 2.75862284 4 4 32.2250000 4.29757684 26.4250000 4.04917687 5 4 25.0250000 3.69898635 20.8000000 3.02103735
Least Squares Means
TRT Y Std Err Pr > |T| LSMEAN LSMEAN LSMEAN H0:LSMEAN=0 Number 1 30.1531125 0.3339174 0.0001 1 2 30.1173006 0.2827350 0.0001 2 3 32.0523296 0.2796295 0.0001 3 4 31.5046854 0.2764082 0.0001 4
表26.2（c ）中结果分析：contrast 语句通过其后的参数项设置，用来假设检验我们自己计划的原假设••••+=+43210:y y y y H ，结果显示非常显著（0.0001<0.05），即入口处底部和顶部均值之和与出口处底部和顶部均值之和是有显著差异的，说明水中的温度不同对牡蛎生长是不同的。

本程序中的estimate 语句，有计划地设计了对入口处的底部和顶部调整后均值进行估计，及它们之差是否为0的假设检验，结果为不显著。

但如果对未调整均值之差是
否为0进行假设检验，结果却为非常显著。

因此，我们可以看到使用调整后均值进行估计是必要的。

2. 多元协方差分析
例26.2 研究男女儿童的体表面积是否相同。

考虑到儿童的身高和体重对表面积可能有影响，在某地测量了男女各15名初生至3周岁儿童的身高、体重和体表面积，得到测量数据见表26.3所示。

表26.3 3周岁男女儿童的身高、体重和体表面积
男（male ）
女（female ）
身高（x 1）
体重（x 2）
表面积（y ） 身高（x 1）
体重（x 2）
表面积（y ） 54.0 3.00 2446.2 54.0 3.00 2117.3 50.5 2.25 1928.4 53.0 2.25 2200.2 51.0 2.50 2094.5 51.5 2.50 1906.2 56.5 3.50 2506.7 51.0 3.00 1850.3 52.0 3.00 2121.0 51.0 3.00 1632.5 76.0 9.50 3845.9 77.0 7.50 3934.0 80.0 9.00 4380.8 77.0 10.0 4180.4 74.0 9.50 4314.2 77.0 9.50 4246.1 80.0 9.00 4078.4 74.0 9.00 3358.8 76.0 8.00 4134.5 73.0 7.50 3809.7 96.0 13.5 5830.2 91.0 12.0 5358.4 97.0 14.0 6013.6 91.0 13.0 5601.7 99.0 16.0 6410.6 94.0 15.0 6074.9 92.0 11.0 5283.3 92.0 12.0 5299.4 94.0 15.0
6101.6
91.0
12.5
5291.5
程序如下：
proc format;
The SAS System
Dependent Variable: Y
Contrast DF Contrast SS Mean Square F Value Pr > F trt12 vs trt34 1 8.59108077 8.59108077 28.49 0.0001
T for H0: Pr > |T| Std Error of Parameter Estimate Parameter=0 Estimate trt1 adj mean 30.1531125 90.30 0.0001 0.33391743 trt2 adj mean 30.1173006 106.52 0.0001 0.28273504 adj trt diff 0.0358120 0.09 0.9312 0.40722674
value sexname 1=’male’ 2=’female’;
data child;
do i=1 to 15;
do sex=1 to 2;
input x1 x2 y @@;
format sex sexname.;
output;
end;
end;
cards;
54.0 3.00 2446.2 54.0 3.00 2117.3
50.5 2.25 1928.4 53.0 2.25 2200.2
51.0 2.50 2094.5 51.5 2.50 1906.2
56.5 3.50 2506.7 51.0 3.00 1850.3
52.0 3.00 2121.0 51.0 3.00 1632.5
76.0 9.50 3845.9 77.0 7.50 3934.0
80.0 9.00 4380.8 77.0 10.0 4180.4
74.0 9.50 4314.2 77.0 9.50 4246.1
80.0 9.00 4078.4 74.0 9.00 3358.8
76.0 8.00 4134.5 73.0 7.50 3809.7
96.0 13.5 5830.2 91.0 12.0 5358.4
97.0 14.0 6013.6 91.0 13.0 5601.7
99.0 16.0 6410.6 94.0 15.0 6074.9
92.0 11.0 5283.3 92.0 12.0 5299.4
94.0 15.0 6101.6 91.0 12.5 5291.5
;
proc glm data=child;
class sex;
model y=sex x1 x2 /solution;
lsmeans sex /stderr tdiff;
run;
程序说明：本例为带有两个协变量x1和x2，一个分组变量sex的完全随机化设计的多元协方差分析。

data步中为了便于读人数据，sex分组变量取值为1和2，但又为了显示清楚，用format过程自定义了sexname.格式，用于sex变量的显示格式。

在class语句中只能有sex 分组变量，而在model语句中应把观察指标放在等号的左边，分组变量和协变量放在等号的右边，solution选项求回归方程的系数估计。

lsmeans语句求修正后均值，stderr选项求均值的标准误差，tdiff选项求均值对比的t值和p值。

程序输出的主要结果见表26.4所示。

表26.4 单因素的多元协方差分析
表26.4中结果分析：由类型3的平方和计算结果表明，身高、体重对体表面积都有非常显著性的影响（0.0001<0.05，0.0059<0.05），而男、女两性之间无显著性差别（0.0762>0.05）。

由回归分析的结果可知道，与x 1、x 2相对应的公共偏回归系数系数为=1β54.477217、
=2β130.645108，它们与0之间差别的检验结果为p =0.0001和p =0.0059。

男、女两性体表面
积的修正均值分别为52.32694和52.32694，两者之间无显著性差别（p =0.0762）。

The SAS System
General Linear Models Procedure
Dependent Variable: Y
Source DF Sum of Squares Mean Square F Value Pr > F Model 3 68523072.11494280 22841024.03831420 557.41 0.0001 Error 26 1065399.75872373 40976.91379707 Corrected Total 29 69588471.87366650
R-Square C.V. Root MSE Y Mean 0.984690 5.131187 202.42755197 3945.04333333 Source DF Type I SS Mean Square F Value Pr > F SEX 1 714100.40833333 714100.40833333 17.43 0.0003 X1 1 67440016.91708050 67440016.91708050 1645.81 0.0001 X2 1 368954.78952901 368954.78952901 9.00 0.0059 Source DF Type III SS Mean Square F Value Pr > F SEX 1 139769.33971381 139769.33971381 3.41 0.0762 X1 1 938153.70360865 938153.70360865 22.89 0.0001 X2 1 368954.78952901 368954.78952901 9.00 0.0059 T for H0: Pr > |T| Std Error of Parameter Estimate Parameter=0 Estimate INTERCEPT -1118.730592 B -2.25 0.0331 497.2296650 SEX female -136.828607 B -1.85 0.0762 74.0867551 male 0.000000 B . . . X1 54.477217 4.78 0.0001 11.3853803 X2 130.645108 3.00 0.0059 43.5387744 NOTE: The X'X matrix has been found to be singular and a generalized inverse was used to solve。