CDC课件-SAS讲义_协方差分析.doc

协方差分析用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 研究牡蛎在不同温度的水中不同位置上的生长情况。有人做了如下试验：分别在通向发电站的入口处（温度较低）不同位置（底部和表层）和出口处（温度较高）不同位置（底部和表层）及电站附近的深水处（底部和表层的中间）总共5个不同位置点上，随机地各放4袋牡蛎（每袋中有10个），共54=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（出口顶部）29.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 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;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。程序输出的主要结果见表1.21.31.4所示。The SAS System Analysis of Variance ProcedureDependent Variable: YSource DF Sum of Squares Mean Square F Value Pr 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 ProcedureDependent 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 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 EstimateINTERCEPT 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 XX matrix has been 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 |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.3975719 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.557715 -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+31)/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检验，获得其他多重比较的检验结果。从最后的55修正均值比较结果表中，可得到（）中的任何一个与（）中的任何一个之间有差别有统计学意义。表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.1531125 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语句通过其后的参数项设置，用来假设检验我们自己计划的原假设，结果显示差别有统计学意义（0.0001 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 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.78952901 9.00 0.0059 T for H0: Pr |T| Std Error ofParameter Estimate Parameter=0 EstimateINTERCEPT -1118.730592 B -2.25 0.0331 497.2296650SEX female -136.828607 B -1.85 0.0762 74.0867551 male 0.000000 B . . .X1 54.477217 4.78 0.0001 11.3