上海財經(jīng)大學經(jīng)濟信息管理系IS/SHUFEPagePAGE1ofNUMPAGES8協(xié)方差分析當定量的影響因素對觀察結果有難以控制的影響，甚至還有交互作用時，采用協(xié)方差分析，這些影響變量稱為協(xié)變量，扣除（或消除）協(xié)變量的影響，可以得到修正后的均值估計。協(xié)方差分析概述協(xié)方差分析概念協(xié)方差分析（analysisofcovariance）又稱帶有協(xié)變量的方差分析（analysisofvariancewithcovariates），是將回歸分析與方差分析結合起來使用的一種分析方法。在各種試驗設計中，對主要變量y研究時，常常希望其他可能影響和干擾y的變量保持一致以到達均衡或可比，使試驗誤差的估計降到最低限度，從而可以準確地獲得處理因素的試驗效應。但是有時，這些變量難以控制，或者根本不能控制。為此需要在試驗中同時記錄這些變量的值，把這些變量看作自變量，或稱協(xié)變量（covariate），建立因變量y隨協(xié)變量變化的回歸方程，這樣就可以利用回歸分析把因變量y中受協(xié)變量影響的因素扣除掉，從而，能夠較合理地比較定性的影響因素處在不同水平下，經(jīng)過回歸分析手段修正以后的因變量的總體均值之間是否有顯著性的差別。簡單地說，協(xié)方差分析是扣除協(xié)變量的影響，或者將這些協(xié)變量處理成相等，再對修正的y的均值作方差分析。協(xié)方差分析的假定協(xié)方差分析需要滿足的假定為：=1\*GB3①各樣本來自具有相同方差的正態(tài)分布總體，即要求各組方差齊性。=2\*GB3②協(xié)變量與主要變量y間的總體回歸系數(shù)不等于0。=3\*GB3③各組的回歸線平等，即回歸系數(shù)如果上述的假定滿足，就作協(xié)方差分析。前述的各種試驗設計，如完全隨機化設計、隨機區(qū)組設計、析因設計、拉丁方設計等，都可以帶一個或多個協(xié)變量，按設計方案扣除協(xié)變量的影響后，對主要變量y的修正均值作比較，得出統(tǒng)計結論。協(xié)方差分析的模型最簡單的單因素一元協(xié)方差分析的模型，是由單因素效應模型加上協(xié)變量的影響因素而得出：(26.1)其中，為協(xié)變量，為協(xié)變量在分類水平i和j上的記錄值，為所有協(xié)變量的平均值，為相關的回歸系數(shù)。設，為平均截距。上式可以化簡成：(26.2)設，上式可以化簡成：(26.3)很明顯是第i組回歸線的截距，等于回歸線的平均截距加上本組的效應。這個式揭示了，觀察值的模型可以表示成一組相似的回歸線，且各組具有共同的回歸系數(shù)，和各組自己的截距。用SAS中的glm過程進行協(xié)方差分析時，要注意不同試驗設計時class語句和model語句的寫法。設分類變量為A、B，協(xié)變量為X，觀察值為Y，則有：=1\*GB3①單因素k水平設計的協(xié)方差分析模型classA;modelXA;=2\*GB3②隨機區(qū)組設計的協(xié)方差分析模型classAB;modelXAB;=3\*GB3③兩因素析因設計的協(xié)方差分析模型classAB;modelXABA*B;實例分析一元協(xié)方差分析例26.1研究牡蠣在不同溫度的水中不同位置上的生長情況。有人做了如下試驗：分別在通向發(fā)電站的入口處（溫度較低）不同位置（底部和表層）和出口處（溫度較高）不同位置（底部和表層）及電站附近的深水處（底部和表層的中間）總共5個不同位置點上，隨機地各放4袋牡蠣（每袋中有10個），共5×4=20袋。在將每袋牡蠣放入位置點之前，先洗干凈稱出每袋的初始體重，放在5個不同點一個月后再稱出最后體重。試驗結果數(shù)據(jù)如表26.1所示。表26.1牡蠣在不同溫度和位置上的生長數(shù)據(jù)位置trt重復數(shù)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程序如下：datagrowth;dotrt=1to5;dorep=1to4;inputxy@@;output;end;end;cards;27.232.632.036.633.037.726.831.028.633.826.831.726.530.726.830.428.635.222.429.123.228.924.430.229.335.021.827.030.336.424.330.520.424.619.623.425.130.318.121.8;procanovadata=growth;classtrt;modely=trt;procglmdata=growth;classtrt;modely=trtx/solution;meanstrt;lsmeanstrt/stderrtdiff;contrast'trt12vstrt34'trt-1-1110;estimate'trt1adjmean'intercept1trt10000x25.76;estimate'trt2adjmean'intercept1trt01000x25.76;estimate'adjtrtdiff'trt1-1000;estimate'trt1unadjmean'intercept1trt10000x29.75;estimate'trt2unadjmean'intercept1trt01000x27.175;estimate'unadjtrtdiff'trt1-1000x2.575;run;程序說明：定性變量trt的5個不同位置點對y可能有較大的影響，因此class語句中分組變量為trt，先選用anova過程進行方差分析。然而，牡蠣的初始體重x對牡蠣的最后體重y可能也有一定的影響，故適合選用glm過程進行協(xié)方差分析，在model語句中不僅包括分組變量trt，而且應包括協(xié)變量x。選擇項solution要求輸出回歸系數(shù)的估計值及其標準誤差和假設檢驗等結果。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的協(xié)方差分析TheSASSystemTheSASSystemAnalysisofVarianceProcedureDependentVariable:YSourceDFSumofSquaresMeanSquareFValuePr>FModel4198.4070000049.601750004.640.0122Error15160.2625000010.68416667CorrectedTotal19358.66950000R-SquareC.V.RootMSEYMean0.55317510.597063.2686643630.84500000SourceDFAnovaSSMeanSquareFValuePr>FTRT4198.4070000049.601750004.640.0122GeneralLinearModelsProcedureDependentVariable:YSourceDFSumofSquaresMeanSquareFValuePr>FModel5354.4471767570.88943535235.050.0001Error144.222323250.30159452CorrectedTotal19358.66950000R-SquareC.V.RootMSEYMean0.9882281.7804380.5491762230.84500000SourceDFTypeISSMeanSquareFValuePr>FTRT4198.4070000049.60175000164.470.0001X1156.04017675156.04017675517.380.0001SourceDFTypeIIISSMeanSquareFValuePr>FTRT412.089359283.0223398210.020.0005X1156.04017675156.04017675517.380.0001TforH0:Pr>|T|StdErrorofParameterEstimateParameter=0EstimateINTERCEPT2.494859769B2.430.02931.02786287TRT1-0.244459378B-0.420.67800.576581962-0.280271345B-0.570.57860.4929082531.654757698B3.850.00180.4294303641.107113519B2.350.03420.4717511250.000000000B...X1.08317981922.750.00010.04762051NOTE:TheX'Xmatrixhasbeenfoundtobesingularandageneralizedinversewasusedtosolvethenormalequations.Estimatesfollowedbytheletter'B'arebiased,andarenotuniqueestimatorsoftheparameters.表26.2(a)中結果分析：對分組變量trt的方差分析表明，即使當初始體重x不考慮，各分組最后體重均值的區(qū)別也統(tǒng)計顯著（0.0122<0.05），其中分組變量trt的平方和為198.40700000。而在協(xié)方差分析中，分組變量trt的類型1的平方和等于方差分析中的平方和198.40700000，分組變量trt的類型3的平方和為12.08935928，大大小于類型1的平方和，是因為類型3的平方和反映了經(jīng)過共同的協(xié)變量x調整后的平方和，減去了協(xié)變量的影響，所以平方和大幅減小。類型1是一種未經(jīng)過調整的平方和，因為它的優(yōu)先級高于協(xié)變量的調整。更進一步分析，我們注意到方差分析中均方誤差為10.68416667，而協(xié)方差分析中卻縮小到0.30159452，相應地分組變量trt的F統(tǒng)計量從4.64增加到10.02，說明包含了協(xié)變量后分組的區(qū)別更加顯著，原因是簡單方差分析中，大多數(shù)的誤差是由于初始體重x的變異造成的。表中的最后一部分是選擇項solution的輸出結果，對模型中的截距、各分組變量和協(xié)變量的回歸系數(shù)進行估計和檢驗，在這個單因素trt的情況下，估計是以最后一個水平trt5（trt=5）為對照組，并且設置它的系數(shù)為0，因此截距intercept的估計值是分組trt5的估計值。其他四個分組trt的系數(shù)估計是每一個與trt5進行比較而得到的。注意，出口處的trt3和trt4分組不同于trt5分組。協(xié)變量x的系數(shù)是合并各組內y和x所得到的回歸系數(shù)，即由5個獨立的trt分組分別回歸y和x后得到回歸系數(shù)然后加權平均。協(xié)變量x的系數(shù)估計值表明，初始體重變動1個單位，最后，體重y相關地要變動1.083179819單位。表26.2（b）未調整均值和調整均值及均值之間的比較TheSASSystemTheSASSystemGeneralLinearModelsProcedureLevelofYXTRTNMeanSDMeanSD1434.47500003.1889130929.75000003.205724052431.65000001.5373136727.17500000.960468643430.85000002.9557852924.65000002.758622844432.22500004.2975768426.42500004.049176875425.02500003.6989863520.80000003.02103735LeastSquaresMeansTRTYStdErrPr>|T|LSMEANLSMEANLSMEANH0:LSMEAN=0Number130.15311250.33391740.00011230.11730060.28273500.00012332.05232960.27962950.00013431.50468540.27640820.00014530.39757190.36219880.00015TforH0:LSMEAN(i)=LSMEAN(j)/Pr>|T|i/j123451.0.087941-4.1466-3.22289-0.423980.93120.00100.00610.67802-0.08794.-4.76003-3.55771-0.568610.93120.00030.00320.578634.1465994.76003.1.3780023.8533780.00100.00030.18980.001843.2228923.557715-1.378.2.3468170.00610.00320.18980.034250.423980.568608-3.85338-2.34682.0.67800.57860.00180.0342NOTE:Toensureoverallprotectionlevel,onlyprobabilitiesassociatedwithpre-plannedcomparisonsshouldbeused.表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.475－1.083179819×（29.75－25.76）。tdiff選擇項要求對已調整均值的兩兩比較采用lsd檢驗，可以使用adjust=duncan/waller等選項替代lsd檢驗，獲得其他多重比較的檢驗結果。從最后的5×5修正均值比較結果表中，可得到（）中的任何一個與（）中的任何一個之間有顯著或非常顯著性差別。表26.2（c）有計劃的均值對比和參數(shù)估計TheSASSystemTheSASSystemDependentVariable:YContrastDFContrastSSMeanSquareFValuePr>Ftrt12vstrt3418.591080778.5910807728.490.0001TforH0:Pr>|T|StdErrorofParameterEstimateParameter=0Estimatetrt1adjmean30.153112590.300.00010.33391743trt2adjmean30.1173006106.520.00010.28273504adjtrtdiff0.03581200.090.93120.40722674trt1unadjmean34.4750000125.550.00010.27458811trt2unadjmean31.6500000115.260.00010.27458811unadjtrtdiff2.82500007.270.00010.38832623表26.2（c）中結果分析：contrast語句通過其后的參數(shù)項設置，用來假設檢驗我們自己計劃的原假設，結果顯示非常顯著（0.0001<0.05），即入口處底部和頂部均值之和與出口處底部和頂部均值之和是有顯著差異的，說明水中的溫度不同對牡蠣生長是不同的。本程序中的estimate語句，有計劃地設計了對入口處的底部和頂部調整后均值進行估計，及它們之差是否為0的假設檢驗，結果為不顯著。但如果對未調整均值之差是否為0進行假設檢驗，結果卻為非常顯著。因此，我們可以看到使用調整后均值進行估計是必要的。多元協(xié)方差分析例26.2研究男女兒童的體表面積是否相同?？紤]到兒童的身高和體重對表面積可能有影響，在某地測量了男女各15名初生至3周歲兒童的身高、體重和體表面積，得到測量數(shù)據(jù)如表26.3所示。表26.33周歲男女兒童的身高、體重和體表面積男（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.474.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程序如下：procformat;valuesexname1=’male’2=’female’;datachild;doi=1to15;dosex=1to2;inputx1x2y@@;formatsexsexname.;output;end;end;cards;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.474.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;procglmdata=child;classsex;modely=sexx1x2/solution;lsmeanssex/stderrtdiff;run;程序說明：本例為帶有兩個協(xié)變量x1和x2，一個分組變量sex的完全隨機化設計的多元協(xié)方差分析。data步中為了便于讀入數(shù)據(jù)，sex分組變量取值為1和2，但又為了顯示清楚，用format過程自定義了sexname.格式，用于sex變量的顯示格式。在class語句中只能有sex分組變量，而在model語句中應把觀察指標放在等號的左邊，分組變量和協(xié)變量放在等號的右邊，solution選項求回歸方程的系數(shù)估計。lsmeans語句求修正后均值，stderr選項求均值的標準誤差，tdiff選項求均值對比的t值和p值。程序輸出的主要結果如表26.4所示。表26.4單因素的多元協(xié)方差分析TheSASSystemTheSASSystemGeneralLinearModelsProcedureDependentVariable:YSourceDFSumofSquaresMeanSquareFValuePr>FModel368523072.1149428022841024.03831420557.410.0001Error261065399.7587237340976.91379707CorrectedTotal2969588471.87366650R-SquareC.V.RootMSEYMean0.9846905.131187202.427551973945.04333333SourceDFTypeISSMeanSquareFValuePr>FSEX1714100.40833333714100.4083333317.430.0003X1167440016.9170805067440016.917080501645.810.0001X21368954.78952901368954.789529019.000.0059SourceDFTypeIIISSMeanSquareFValuePr>FSEX1139769.33971381139769.339713813.410.0762X11938153.70360865938153.7036086522.890.0001X21368954.78952901368954.789529019.000.005