R语言学习系列28-协方差分析

23.协方差分析一、基本原理1.基本思想在实际问题中，有些随机因素是很难人为控制的，但它们又会对结果产生显著影响。如果忽略这些因素的影响，则有可能得到不正确的结论。这种影响的变量称为协变量（一般是连续变量）。例如，研究3种不同的教学方法的教学效果的好坏。检查教学效果是通过学生的考试成绩来反映的，而学生现在考试成绩是受到他们自身知识基础的影响，在考察的时候必须排除这种影响。协方差分析将那些难以控制的随机变量作为协变量，在分析中将其排除，然后再分析控制变量对于观察变量的影响，从而实现对控制变量效果的准确评价。协方差分析要求协变量应是连续数值型，多个协变量间互相独立，且与控制变量之间没有交互影响。前面单因素方差分析和多因素方差分析中的控制变量都是一些定性变量，而协方差分析中既包含了定性变量（控制变量），又包含了定量变量（协变量）。协方差分析在扣除协变量的影响后再对修正后的主效应进行方差分析，是一种把直线回归或多元线性回归与方差分析结合起来的方法，其中的协变量一般是连续性变量，并假设协变量与因变量间存在线性关系，且这种线性关系在各组一致，即各组协变量与因变量所建立的回归直线基本平行。当有一个协变量时，称为一元协方差分析，当有两个或两个以上的协变量时，称为多元协方差分析。2.协方差分析需要满足的条件（1）自变量是分类变量，协变量是定距变量，因变量是连续变量；对连续变量或定距变量的协变量的测量不能有误差；（2）协变量与因变量之间的关系是线性关系，可以用协变量和因变量的散点图来检验是否违背这一假设；协变量的回归系数（即各回归线的斜率）是相同的，且不等于0，即各组的回归线是非水平的平行线。否则，就有可能犯第一类错误，即错误地接受虚无假设；（3）自变量与协变量相互独立，若协方差受自变量的影响，那么协方差分析在检验自变量的效应之前对因变量所作的控制调整将是偏倚的，自变量对因变量的间接效应就会被排除；（4）各样本来自具有相同方差σ2的正态分布总体，即要求各组方差齐性。二、协方差理论1.观测值=均值+分组变量影响+协变量影响+随机误差.即()ijiijijyutxx（1）其中，X为所有协变量的平均值。注：在方差分析中，协变量影响是包含在随机误差中的，在协方差分析中需要分离出来。用协变量进行修正，得到修正后的yij(adj)为(adj)()ijijijiijyyxxut就可以对yij(adj)做方差分析了。关键问题是求出回归系数β.2.总离差=分组变量离差+协变量离差+随机误差，（1）计算总离差平方和时，记11()()knxyijijijTxxyy211()knxxijijTxx总离差平方和：211()knyyijijTyy最终要检验分组自变量对因变量有无显著作用。原假设H0：无显著作用。假设检验是在H0为真条件下进行，可认为ti=0，则()ijTijijyuxx按最小二乘法原理线性回归可得到β的估计值ˆxyTxxTT记修正的总离差平方和（残差平方和）为Tyy(adj)，则22(adj)ˆTxyyyyyxxyyxxTTTTTT，自由度为n-2注：2ˆTxxT为回归平方和，若ˆ0T（回归线为水平线），表示协变量x对y无作用，用方差分析就可以解决了。（2）计算组内离差平方和时，记11()()knxyijiijiijExxyy211()knxxijiijExx组内总离差平方和：211()knyyijiijEyy根据协方差分析的基本假设：各组内回归系数相等（做协方差分析时需要检验这一点），得到组内回归系数βw的估计值ˆxywxxEE记修正的组内总离差平方和（组内残差平方和）为Eyy(adj),则22(adj)ˆxyyyyywxxyyxxEEEEEE，自由度为n-k-1其中，2ˆwxxE为组内回归平方和，当1ˆˆwwk时，组内总离差平方和认为完全是由随机因素引起的，Eyy(adj)就是随机为误差。这里的ˆw是1ˆˆ,,wwk的加权平均值。（3）计算分组变量离差平方和Byy(adj)，它反映的是各个水平之间的差异。2(adj)(adj)(adj)(adj)ˆTyyyyyyyyxxyyBTETTE即，分组变量离差=总离差-协变量离差-随机误差。于是，就可以进行组间无差异检验了：(adj)(adj)/1/1yyyyBkFEnk3.因此，在做协方差分析前，需要依次做两个假设检验：（1）协变量对因变量的影响对与各组来说都是相同的，即各组回归系数相等：1ˆˆˆ:wwkw;步骤：①先按回归系数相等和不相等分别表示模型()ijiwijijyutxx()ijiwiijijyutxx并计算出误差平方和2(adj)yyyywxxEEE211ikyywixxiSEE其中，1ikyyyyiEE.②计算F值(adj)11/1/2yyESkFSnk若F值小于临界值Fα，则说明各组回归系数无显著差异（相等）。（2）这些相等的回归系数ˆ0w.即采用一元线性回归的显著性检验，2(adj)/1=//(1)wxxyyEFEnk回归平方和/自由度残差平方和自由度2222/(1)(/)/(1)xyxxxyyyxyxxyyxxxyEEEnkEEEnkEEE4.协方差分析的步骤（1）检验数据是否满足假设条件：正态分布性、方差齐性、各分组通过协变量预测因变量的回归斜率相同；（2）检验效应因子的显著性；（3）估计校正的组均值；（4）检验校正的组均值之间的差异。三、R语言实现协方差分析要求数据满足：正态性、方差齐性、各分组通过协变量预测因变量的回归斜率相同。R语言用aov()函数进行协方差分析，基本格式为：aov(formula,data,...)其中，data为数据框；formula为协方差公式形式，形如y~x+A,x为连续型协变量，A为组别因子。例1研究分别接受了3种不同的教学方法的3组学生，在数学成绩上是否有显著差异，数据文件“ex28_cov.Rdata”。先不考虑数学入学成绩，只以“教学方法”为分组变量，“后测成绩”为因变量进行单因素方差分析：setwd(E:/办公资料/R语言/R语言学习系列/codes)load(ex28_cov.Rdata)head(scores)beforeafterteach139681238631351651456681574741640601attach(scores)table(teach)#各组的样本数teach123303233aggregate(after,by=list(teach),mean)#各组均值Group.1x1162.883332272.671883365.06061shapiro.test(after)#正态性检验Shapiro-Wilknormalitytestdata:afterW=0.99105,p-value=0.7772bartlett.test(after~teach,data=scores)#方差齐性检验Bartletttestofhomogeneityofvariancesdata:afterbyteachBartlett'sK-squared=0.69854,df=2,p-value=0.7052fit.aov-aov(after~teach,data=scores)summary(fit.aov)DfSumSqMeanSqFvaluePr(F)teach21662830.810.448.23e-05***Residuals92732579.6---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1说明：单因素方差分析的p值=8.23e-05,远小于0.05,表明，两种教学方法有非常显著的差异。但是，后测成绩肯定会受到前测成绩（连续型）的影响，假定前测成绩与教学方法（即组别，是控制变量）不存在交互影响。因此，将后测成绩作为因变量；教学方法作为控制变量；前测成绩作为协变量进行协方差分析。回归斜率相同检验，即前测成绩与后测成绩的回归线是否平行：scores1-subset(scores,teach==1)scores2-subset(scores,teach==2)scores3-subset(scores,teach==3)par(mfrow=c(1,3))plot(scores1$before,scores1$after,xlab=before,ylab=after,main=teach=1)abline(lm(after~before,data=scores1))plot(scores2$before,scores2$after,xlab=before,ylab=after,main=teach=2)abline(lm(after~before,data=scores2))plot(scores3$before,scores3$after,xlab=before,ylab=after,main=teach=3)abline(lm(after~before,data=scores3))可见两组的直线趋势的斜率比较接近（平行），基本符合协方差假定。除了图形判断外，还可以通过交互作用是否显著，来判断斜率是否相同。因为若前验成绩与教学方法的交互作用显著，则说明前验成绩与后验成绩的关系，依赖于教学方法。library(multcomp)fit2-aov(after~before*teach,data=scores)summary(fit2)DfSumSqMeanSqFvaluePr(F)before124322432.035.3915.22e-08***teach2362180.92.6330.0775.before:teach27638.20.5560.5752Residuals89611668.7---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1说明：交互项的p值=0.57520.05,故不显著，支持了斜率相同的假设。fit-aov(after~before+teach,data=scores)#协方差分析summary(fit)DfSumSqMeanSqFvaluePr(F)before124322432.035.7394.35e-08***teach2362180.92.6590.0755.Residuals91619268.0---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1说明：协方差分析的结果表明：前测成绩的p值=4.35e-08远小于0.05,说明“前测成绩”对“后测成绩产生了非常显著的影响；“教学方法”的p值=0.07550.05,说明“教学方法”对“后测成绩”的影响不显著。由于受到协变量的影响，我们希望获取调整后的各组均值——即去除协变量效应后的各组均值。可使用effects包中的effects()函数来计算调整的均值：library(effects)effect(teach,fit)teacheffectteach12365.0572270.0495865.62718与方差分析时一样，要想得到教学方法两两之间有无差异，可以均值的成对比较（略）。下面讲一下自定义比较（使用multcomp包可以实现），例如，分组变量有4个水平ABCD，要比较A与D时，比较矩阵=[100-1]T,有[ABCD]×[100-1]T=0等价于A=D要想将A与D合并再与B比较有无差异，则可以指定L矩阵=[1-201]T,则[ABCD]×[1-201]T=0等价于(A+D)/2=B注意：是从(A+D)/2=B倒推比较矩阵，该式即A-2B+0C+D=0.根据调整后的各组均值，教学方法1