SAS系统和数据分析协方差分析.doc

第二十六课 协方差分析当定量的影响因素对观察结果有难以控制的影响，甚至还有交互作用时，采用协方差分析，这些影响变量称为协变量，扣除（或消除）协变量的影响，可以得到修正后的均值估计。一、 协方差分析概述1. 协方差分析概念协方差分析（analysis of covariance）又称带有协变量的方差分析（analysis of variance with covariates），是将回归分析与方差分析结合起来使用的一种分析方法。在各种试验设计中，对主要变量y研究时，常常希望其他可能影响和干扰y的变量保持一致以到达均衡或可比，使试验误差的估计降到最低限度，从而可以准确地获得处理因素的试验效应。但是有时，这些变量难以控制，或者根本不能控制。为此需要在试验中同时记录这些变量的值，把这些变量看作自变量，或称协变量（covariate），建立因变量y随协变量变化的回归方程，这样就可以利用回归分析把因变量y中受协变量影响的因素扣除掉，从而，能够较合理地比较定性的影响因素处在不同水平下，经过回归分析手段修正以后的因变量的总体均值之间是否有显著性的差别。简单地说，协方差分析是扣除协变量的影响，或者将这些协变量处理成相等，再对修正的y的均值作方差分析。2. 协方差分析的假定协方差分析需要满足的假定为：各样本来自具有相同方差的正态分布总体，即要求各组方差齐性。协变量与主要变量y间的总体回归系数不等于0。各组的回归线平等，即回归系数如果上述的假定满足，就作协方差分析。前述的各种试验设计，如完全随机化设计、随机区组设计、析因设计、拉丁方设计等，都可以带一个或多个协变量，按设计方案扣除协变量的影响后，对主要变量y的修正均值作比较，得出统计结论。3. 协方差分析的模型最简单的单因素一元协方差分析的模型，是由单因素效应模型加上协变量的影响因素而得出：(26.1)其中，为协变量，为协变量在分类水平i和j上的记录值，为所有协变量的平均值，为相关的回归系数。设，为平均截距。上式可以化简成：(26.2)设，上式可以化简成：(26.3)很明显是第i组回归线的截距，等于回归线的平均截距加上本组的效应。这个式揭示了，观察值的模型可以表示成一组相似的回归线，且各组具有共同的回归系数，和各组自己的截距。用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个），共54=20袋。在将每袋牡蛎放入位置点之前，先洗干净称出每袋的初始体重，放在5个不同点一个月后再称出最后体重。试验结果数据如表26.1所示。表26.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。程序输出的主要结果如表26.2（a）、表26.2（b）、表26.2（c）所示。表26.2（a） 单因素trt一元x的协方差分析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.表26.2(a)中结果分析：对分组变量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.表26.2（b）中结果分析：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的均值，或称最小二乘均值估计，我们可以由公式（26.1）求分组平均得到：(26.4)再由公式（25.2）求分组平均代入上式：(26.5)例如，初始体重的整体平均值为=（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修正均值比较结果表中，可得到（）中的任何一个与（）中的任何一个之间有显著或非常显著性差别。表26.2（c） 有计划的均值对比和参数估计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表26.2（c）中结果分析：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.3853803X2 130.645108 3.00 0.0059 43.5387744NOTE: 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. Least Squares Means SEX Y Std Err Pr |T| T / Pr |T| H0: LSMEAN LSMEAN H0:LSMEAN=0 LSMEAN1=LSMEAN2 female 3876.62903 52.32694 0.0001 -1.84687 0.0762 male 4013.45764 52.32694 0.0001表26.4中结果分析：由类型3的平方和计算结果表明，身高、体重对体表面积都有非常显著性的影响（0.00010.05，0.00590.05）。由回归分析的结果可知道，与x1、x2相对应的公共偏回归系数为54.477217、130.645108，它们与0之间差别的检验结果为p=0.0001和p=0.0059。男、女两性体表面积的修正均值分别为52.32694和52.32694，两者之间无显著性差别（p=0.0762）。第二十七课 符号检验和Wilcoxon符号秩检验在统计推断和假设检验中，传统的检验统计量都叫做参数检验，因为它们都依赖于确定的概率分布，这个分布带有一组自由的参数。参数检验被认为是依赖于分布假定的。通常情况下，我们对数据进行分析时，总是假定误差项服从正态分布，这是人们易于接受的事实，因为正态分布的原始出发点就是来自于误差分布，至于当样本相当大时，数据的正态近似，这是由于大样本理论所保证的。但有些资料不一定满足上述要求，或不能测量具体数值，其观察结果往往只有程度上的区别，如颜色的深浅、反应的强弱等，此时就不适用参数检验的方法，而只能用非参数统计方法（non-parametric statistical analysis）来处理。这种方法对数据来自的总体不作任何假设或仅作极少的假设，因此在实用中颇有价值，适用面很广。三、 单样本的符号检验符号检验（sign test）是一种最简单的非参数检验方法。它是根据正、负号的个数来假设检验。首先需要将原始观察值按设定的规则，转换成正、负号，然后计数正、负号的个数作出检验。该检验可用于样本中位数和总体中位数的比较，数据的升降趋势的检验，特别适用于总体分布不服从正态分布或分布不明的配对资料，有时当配对比较的结果只能定性的表示，如试验前后比较结果为颜色从深变浅、程度从强变弱，成绩从一般变优秀，即不能获得具体数字，也可用符号检验，例如用正号表示颜色从深变浅，用负号表示颜色从浅变深。用于配对资料时，符号检验的计算步骤为：首先定义成对数据指定正号或负号的规则，然后计数正号的个数及负号的个数，由于在具体比较配对资料时，可能存在配对资料的前后没有变化，或等于假设中的中位数，此时仅需要将这些观察值从资料中剔除，当然样本大小也随之减少，故修正样本大小。当样本较小时，应使用二项分布确切概率计算法，当样本较大时，常利用二项分布的正态近似。1. 小样本时的二项分布概率计算当时，或的检验值由精确计算尺度二项分布的卷积获得。在比较配对资料试验前后有否变化，或增加或减小的假设检验时，如果我们定义试验后比试验前增加为正号，反之为负号，那么对于原假设：试验前后无变化来说，正号的个数和负号的个数可能性应当相等，即正号出现的概率=0.5，于是与均服从二项分布，对于太大的相应太小的，或者太大的相应太小的，都将拒绝接受原假设；对于原假设：试验后比试验前有增加来说，正号的个数大于负号的个数的可能性应该大，即正号出现的概率，对于太小的相应太大的，将拒绝接受原假设；对于原假设：试验后比试验前减小来说，正号的个数小于等于负号的个数的可能性应该大，即正号出现的概率，对于太大的相应太小的，将拒绝接受原假设。例27.1有一种提高学生某种素质的训练，有人说它是无效的，有人说它是有效的，那么真实情况究竟应该是怎样的呢？随机地选取15名学生作为试验样本，在训练开始前做了一次测验，每个学生的素质按优、良、中、及格、差打分，经过三个月训练后，再做一次测试对每个学生打分。数据如表27.1所示。我们将素质提高用正号表示，反之用负号表示，没有变化用0表示。显著性水平取0.1。表27.1 训练前后的素质比较学生编号训练之前训练之后差异符号1中优2及格良3良中4差中5良良06中优7差及格8良优9中差10差中11中优12及格良13中及格14中优15差中从表27.1中15名学生训练前后的差异分析可得出：有14名学生有差异，其中=11，=3。1名学生无差异（学生编号为5），应该从分析中去掉，所以=151=14。假设检验为： 即训练之后学生素质没有提高。 即训练之后学生素质有提高。由于试验的结果只有两种可能，正号或负号，对每一个学生试验出现正号的假定概率为=0.5，负号为1-=0.5，这样整个试验的概率是相同的，并且每一个试验是相互独立的。因此在=14次独立的试验中，正号出现的次数服从二项分布，如表27.2所示。表27.2 二项分布的概率和累计概率n=14,p=0.5正号出现的次数正号出现的概率累计概率00.00010.000110.00090.000920.00560.006530.02220.028740.06110.089850.12220.212060.18330.395370.20950.604780.18330.788090.12220.9102100.06110.9713110.02220.9935120.00560.9991130.00090.9999140.00011.0000从表27.2的累计概率列中我们看到，正号出现的次数大于10的概率为10.9713=0.0287，或者换一种方法计算为=0.0001+0.0009+0.0056+0.0222=0.0287，二者的微小差异是因为小数点后舍入问题造成的。而试验的结果：正号出现的次数为11，大于10，出现的概率不会超过0.0287，我们开始设定的显著性水平为0.1，由于0.02870.1，所以我们拒绝原假设，接受备选假设。如果我们的原假设为=0.5，既训练前后学生素质相等，那么就是双侧检验，应该加上正号出现的次数小于4的概率0.0287，即20.0287=0.0574时，应该修正为0.5；当时，应该修正为0.5。值加或减的0.5是连续性修正因子，目的是为了能将连续分布应用到近似的离散型分布。四、 配对资料的Wilcoxon符号秩检验当两组配对资料近似服从正态分布，它们差值的检验可以使用配对t检验法。如果配对资料的正态分布的假设不能成立，就可以使用Frank Wilcoxon（1945）符号秩检验，它是一种非参数检验方法，对配对资料的差值采用符号秩方法来检验。它的基本要求是差值数据设置为最小的序列等级和两组配对资料是相关的（配成对）。在两组配对资料的差异有具体数值的情况下，符号检验只利用大于0和小于0的信息，即正号和负号的信息，而对差异大小所包含的信息却未加利用，但Wilcoxon符号秩检验方法既考虑了正、负号，又利用了差值大小，故效率较符号检验法高。例27.2某制造商想要比较两种不同的生产方法所花费的生产时间是否有差异。随机地选取了11个工人，每一个工人都分别使用两种不同的生产方法