CDC课件SAS讲义协方差分析

1、协方差分析用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. 一元协方差分析例1 研究牡蛎在不同温度的水中不同位置上的生长情况。有人做了如下试验：分别在通向发电站的入口处（温度较低）不同位置（底部和表层）和出口处（温度较高）不同位置（底部和表层）及电站附近的深水处

2、（底部和表层的中间）总共5个不同位置点上，随机地各放4袋牡蛎（每袋中有10个），共5×4=20袋。在将每袋牡蛎放入位置点之前，先洗干净称出每袋的初始体重，放在5个不同点一个月后再称出最后体重。试验结果数据见表26.1所示。表1.1 牡蛎在不同温度和位置上的生长数据位置trt重复数rep（x为初始体重，y为最后体重）1234xyxyxyxy1（入口底部）27.232.632.036.633.037.726.831.02（入口顶部）28.633.826.831.726.530.726.830.43（出口底部）28.635.222.429.123.228.924.430.24（出口顶部）2

3、9.335.021.827.030.336.424.330.55（附近中部）20.424.619.623.425.130.318.121.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.028.6 33.8 26.8 31.7 26.5 30.7 26.8 30.428.6 35.2 22.4 29.1 23.2 28.9 24.4 30.229.3 35.0 21.8 27.0 30.3 36.4 24.3

4、30.520.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;est

5、imate '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

6、 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的各修正均

7、值之间两两比较所对应的t值和p值。Contrast语句是用来比较入口处底部和顶部均值之和与出口处底部和顶部均值之和是否相等。前三条estimate语句是用来估计入口处底部和顶部调整后的均值及它们之差，并假设检验是否为0，后三条estimate语句是用来估计入口处底部和顶部未调整的均值及它们之差，并假设检验是否为0。程序输出的主要结果见表1.21.31.4所示。The SAS System Analysis of Variance ProcedureDependent Variable: YSource DF Sum of Squares Mean Square F Value Pr >

8、FModel 4 198.40700000 49.60175000 4.64 0.0122Error 15 160.26250000 10.68416667Corrected Total 19 358.66950000 R-Square C.V. Root MSE Y Mean 0.553175 10.59706 3.26866436 30.84500000Source DF Anova SS Mean Square F Value Pr > FTRT 4 198.40700000 49.60175000 4.64 0.0122General Linear Models Procedur

9、eDependent Variable: YSource DF Sum of Squares Mean Square F Value Pr > FModel 5 354.44717675 70.88943535 235.05 0.0001Error 14 4.22232325 0.30159452Corrected Total 19 358.66950000 R-Square C.V. Root MSE Y Mean 0.988228 1.780438 0.54917622 30.84500000Source DF Type I SS Mean Square F Value Pr >

10、; FTRT 4 198.40700000 49.60175000 164.47 0.0001X 1 156.04017675 156.04017675 517.38 0.0001Source DF Type III SS Mean Square F Value Pr > FTRT 4 12.08935928 3.02233982 10.02 0.0005X 1 156.04017675 156.04017675 517.38 0.0001 T for H0: Pr > |T| Std Error ofParameter Estimate Parameter=0 EstimateI

11、NTERCEPT 2.494859769 B 2.43 0.0293 1.02786287TRT 1 -0.244459378 B -0.42 0.6780 0.57658196 2 -0.280271345 B -0.57 0.5786 0.49290825 3 1.654757698 B 3.85 0.0018 0.42943036 4 1.107113519 B 2.35 0.0342 0.47175112 5 0.000000000 B . . .X 1.083179819 22.75 0.0001 0.04762051NOTE: The X'X matrix has been

12、 found to be singular and a generalized inverse was used to solve the normal equations. Estimates followed by the letter 'B' are biased, and are not unique estimators of the parameters.表1.2 单因素trt一元x的协方差分析表1.2中结果分析：对分组变量trt的方差分析表明，即使当初始体重x不考虑，各分组最后体重均值的区别也统计学意义（0.0122<0.05），其中分组变量trt的平方和为

13、198.40700000。而在协方差分析中，分组变量trt的类型1的平方和等于方差分析中的平方和198.40700000，分组变量trt的类型3的平方和为12.08935928，大大小于类型1的平方和，是因为类型3的平方和反映了经过共同的协变量x调整后的平方和，减去了协变量的影响，所以平方和大幅减小。类型1是一种未经过调整的平方和，因为它的优先级高于协变量的调整。更进一步分析，我们注意到方差分析中均方误差为10.68416667，而协方差分析中却缩小到0.30159452，相应地分组变量trt的F统计量从4.64增加到10.02，说明包含了协变量后分组的区别有统计学意义的概率增加，原因是简单方

14、差分析中，大多数的误差是由于初始体重x的变异造成的。表中的最后一部分是选择项solution的输出结果，对模型中的截距、各分组变量和协变量的回归系数进行估计和检验，在这个单因素trt的情况下，估计是以最后一个水平trt5（trt=5）为对照组，并且设置它的系数为0，因此截距intercept的估计值是分组trt5的估计值。其他四个分组trt的系数估计是每一个与trt5进行比较而得到的。注意，出口处的trt3和trt4分组是不同与trt5分组。协变量x的系数是合并各组内y和x所得到的回归系数，即是由5个独立的trt分组，分别回归y和x后得到回归系数然后加权平均。协变量x的系数估计值表明，初始体重

15、变动1个单位最后体重y相关地要变动1.083179819单位。表1.3 未调整均值和调整均值及均值之间的比较The SAS SystemGeneral 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.

16、29757684 26.4250000 4.04917687 5 4 25.0250000 3.69898635 20.8000000 3.02103735Least 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 5 30.397

17、5719 0.3621988 0.0001 5 T for H0: LSMEAN(i)=LSMEAN(j) / Pr > |T| i/j 1 2 3 4 5 1 . 0.087941 -4.1466 -3.22289 -0.42398 0.9312 0.0010 0.0061 0.6780 2 -0.08794 . -4.76003 -3.55771 -0.56861 0.9312 0.0003 0.0032 0.5786 3 4.146599 4.76003 . 1.378002 3.853378 0.0010 0.0003 0.1898 0.0018 4 3.222892 3.557

18、715 -1.378 . 2.346817 0.0061 0.0032 0.1898 0.0342 5 0.42398 0.568608 -3.85338 -2.34682 . 0.6780 0.5786 0.0018 0.0342NOTE: To ensure overall protection level, only probabilities associated with pre-planned comparisons should be used.表1.3中结果分析：means语句要求计算按trt每个水平分组的未调整的y和x均值。如=34.475=(32.6+36.6+37.7+3

19、1)/4，=29.75=(27.2+32+33+26.8)/4。Lsmeans语句要求计算调整后的y的均值，或称最小二乘均值估计，我们可以由公式求分组平均得到：再由公式求分组平均代入上式：例如，初始体重的整体平均值为=（29.750+27.175+24.650+26.425+20.800）/5=25.76，以trt1分组为例，调整后=30.1531125=34.4751.083179819×（29.7525.76）。tdiff选择项要求对已调整均值的两两比较采用lsd检验，可以使用adjust= duncan/waller等选项替代lsd检验，获得其他多重比较的检验结果。从最后的5&

20、#215;5修正均值比较结果表中，可得到（）中的任何一个与（）中的任何一个之间有差别有统计学意义。表1.4 有计划的均值对比和参数估计The SAS SystemDependent Variable: YContrast DF Contrast SS Mean Square F Value Pr > Ftrt12 vs trt34 1 8.59108077 8.59108077 28.49 0.0001 T for H0: Pr > |T| Std Error ofParameter Estimate Parameter=0 Estimatetrt1 adj mean 30.153

21、1125 90.30 0.0001 0.33391743trt2 adj mean 30.1173006 106.52 0.0001 0.28273504adj trt diff 0.0358120 0.09 0.9312 0.40722674trt1 unadj mean 34.4750000 125.55 0.0001 0.27458811trt2 unadj mean 31.6500000 115.26 0.0001 0.27458811unadj trt diff 2.8250000 7.27 0.0001 0.38832623表1.4中结果分析：contrast语句通过其后的参数项设

22、置，用来假设检验我们自己计划的原假设，结果显示差别有统计学意义（0.0001<0.05），即入口处底部和顶部均值之和与出口处底部和顶部均值之和是有统计学意义的，说明水中的温度不同对牡蛎生长是不同的。本程序中的estimate语句，有计划地设计了对入口处的底部和顶部调整后均值进行估计，及它们之差是否为0的假设检验，结果无统计学意义。但如果对未调整均值之差是否为0进行假设检验，结果却有统计学意义。因此，我们可以看到使用调整后均值进行估计是必要的。2. 多元协方差分析（课上没讲）例2 研究男女儿童的体表面积是否相同。考虑到儿童的身高和体重对表面积可能有影响，在某地测量了男女各15名初生至3周岁

23、儿童的身高、体重和体表面积，得到测量数据见表2.1所示。表2.1 3周岁男女儿童的身高、体重和体表面积男（male）女（female）身高（x1）体重（x2）表面积（y）身高（x1）体重（x2）表面积（y）54.03.002446.254.03.002117.350.52.251928.453.02.252200.251.02.502094.551.52.501906.256.53.502506.751.03.001850.352.03.002121.051.03.001632.576.09.503845.977.07.503934.080.09.004380.877.010.04180.47

24、4.09.504314.277.09.504246.180.09.004078.474.09.003358.876.08.004134.573.07.503809.796.013.55830.291.012.05358.497.014.06013.691.013.05601.799.016.06410.694.015.06074.992.011.05283.392.012.05299.494.015.06101.691.012.55291.5程序如下：proc format;value sexname 1=male 2=female;data child;do i=1 to 15;do sex

25、=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.350.5 2.25 1928.4 53.0 2.25 2200.251.0 2.50 2094.5 51.5 2.50 1906.256.5 3.50 2506.7 51.0 3.00 1850.352.0 3.00 2121.0 51.0 3.00 1632.576.0 9.50 3845.9 77.0 7.50 3934.080.0 9.00 4380.8 77.0 10.0 4180.474.0

26、 9.50 4314.2 77.0 9.50 4246.180.0 9.00 4078.4 74.0 9.00 3358.876.0 8.00 4134.5 73.0 7.50 3809.796.0 13.5 5830.2 91.0 12.0 5358.497.0 14.0 6013.6 91.0 13.0 5601.799.0 16.0 6410.6 94.0 15.0 6074.992.0 11.0 5283.3 92.0 12.0 5299.494.0 15.0 6101.6 91.0 12.5 5291.5;proc glm data=child;class sex;model y=s

27、ex 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值

28、和p值。程序输出的主要结果见表2.2所示。The SAS SystemGeneral Linear Models ProcedureDependent Variable: YSource DF Sum of Squares Mean Square F Value Pr > FModel 3 68523072.11494280 22841024.03831420 557.41 0.0001Error 26 1065399.75872373 40976.91379707Corrected Total 29 69588471.87366650 R-Square C.V. Root MSE Y

29、Mean 0.984690 5.131187 202.42755197 3945.04333333Source DF Type I SS Mean Square F Value Pr > FSEX 1 714100.40833333 714100.40833333 17.43 0.0003X1 1 67440016.91708050 67440016.91708050 1645.81 0.0001X2 1 368954.78952901 368954.78952901 9.00 0.0059Source DF Type III SS Mean Square F Value Pr > FSEX 1 139769.33971381 139769.33971381 3.41 0.0762X1 1 938153.70360865 938153.70360865 22.89 0.0001X2 1 368954.78952901 368954.78