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

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

《R语言学习系列28协方差分析.docx》由会员分享，可在线阅读，更多相关《R语言学习系列28协方差分析.docx（11页珍藏版）》请在冰豆网上搜索。

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

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

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

组内总离差平方和：

根据协方差分析的基本假设：

各组内回归系数相等（做协方差分析时需要检验这一点），得到组内回归系数βw的估计值

记修正的组内总离差平方和（组内残差平方和）为Eyy（adj）,则

，自由度为n-k-1

其中，为组内回归平方和，当时，组内总离差平方和认为完全是由随机因素引起的，Eyy（adj）就是随机为误差。

这里的是的加权平均值。

（3）计算分组变量离差平方和Byy（adj），它反映的是各个水平之间的差异。

即，分组变量离差=总离差-协变量离差-随机误差。

于是，就可以进行组间无差异检验了：

3.因此，在做协方差分析前，需要依次做两个假设检验：

（1）协变量对因变量的影响对与各组来说都是相同的，即各组回归系数相等：

;

步骤：

①先按回归系数相等和不相等分别表示模型

并计算出误差平方和

其中，.

②计算F值

若F值小于临界值Fα，则说明各组回归系数无显著差异（相等）。

（2）这些相等的回归系数.

即采用一元线性回归的显著性检验，

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）

beforeafterteach

139681

238631

351651

456681

574741

640601

attach（scores）

table（teach）#各组的样本数

teach

123

303233

aggregate（after,by=list（teach）,mean）#各组均值

Group.1x

1162.88333

2272.67188

3365.06061

shapiro.test（after）#正态性检验

Shapiro-Wilknormalitytest

data:

after

W=0.99105,p-value=0.7772

bartlett.test（after~teach,data=scores）#方差齐性检验

Bartletttestofhomogeneityofvariances

data:

afterbyteach

Bartlett'sK-squared=0.69854,df=2,p-value=0.7052

fit.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.5752

Residuals89611668.7

---

Signif.codes:

0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

说明：

交互项的p值=0.5752>0.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.0755>0.05,说明“教学方法”对“后测成绩”的影响不显著。

由于受到协变量的影响，我们希望获取调整后的各组均值——即去除协变量效应后的各组均值。

可使用effects包中的effects（）函数来计算调整的均值：

library（effects）

effect（"teach",fit）

teacheffect

teach

123

65.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和3基本相同，虽然总体上3种差异不显著，教学方法2与1、3是否有显著差异呢？

那么就需要自定义比较。

library（multcomp）

contrast<-rbind（"2vs.13"=c（-1,2,-1））

res.vs<-glht（fit,linfct=mcp（teach=contrast））

summary（res.vs）

SimultaneousTestsforGeneralLinearHypotheses

MultipleComparisonsofMeans:

User-definedContrasts

Fit:

aov（formula=after~before+teach,data=scores）

LinearHypotheses:

EstimateStd.ErrortvaluePr（>|t|）

2vs.13==09.4154.0832.3060.0234*

---

Signif.codes:

0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

（Adjustedpvaluesreported--single-stepmethod）

说明：

自定义比较：

教学方法2与教学方法1,3是否有差异，设置比较矩阵为[-1,2,-1],结果p值=0.0234<0.05,拒绝原假设，即有显著差异。

另外，HH包中的ancova（）函数，也是用来做协方差分析的，还能将结果可视化。

基本格式为：

ancova（formula,data=,x,groups,...）

其中，若formula不包括，x和groups为作图时需指明协变量和因子。

library（HH）

ancova（after~before+teach,data=scores）

AnalysisofVarianceTable

Response:

after

DfSumSqMeanSqFvaluePr（>F）

before12432.02431.9635.73924.354e-08***

teach2361.9180.932.65890.07546.

Residuals916192.368.05

---

Signif.codes:

0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

从图中可以看出，用前验成绩来预测后验成绩的回归线相互平行，只是截距项不同。

教学方法13基本相同，教学方法2明显好于13。

上述代码会让直线保持平行，若用

ancova（after~before*teach,data=scores）

则生成的图形将允许斜率和截距项依据组别而发生变化，这对可视化那些违背回归斜率同质性的实例非常有用。