supp

R语⾔实战笔记--第九章⽅差分析

R语⾔实战笔记–第九章⽅差分析

标签（空格分隔）：R语⾔⽅差分析

术语

组间因⼦，组内因⼦，⽔平：组间因⼦和组同因⼦的区别是，组间因⼦对所有测试对象进⾏分组，⽽组内因⼦则把所有测试对象归为同

⼀组，⽔平则是因⼦的分类值

单因素⽅差分析，多因素⽅差分析，协⽅差分析，多元⽅差分析，协变量：单因素，多因素都是⼀元⽅差分析，只有⼀个因变量(y),协

⽅差分析也是，多元就是有多个因变量，协变量的意思其实就是不感兴趣或不能控制的变量，把它从⾃变量（可控制变量）中剔除出去的变

量，它代表着每个测试对象的某些初始状态。

均衡设计，⾮均衡设计：分组时，各组的观测数若相同，则为均衡设计，否则为⾮均衡设计。

下⾯看两个图表，代表的是因⼦数、协变量、因变量的数⽬不同时，⽅差的叫法不同，以及⼀个书上的例⼦。

ANOVA模型拟合

模型拟合的函数⽅法是aov(formula,data=dataframe),其中formula的公式与回归拟合中的格式⼀样，只是少了⼀些幂级及变量替换

的数据。

另外，需要⼗分注意的是，在⽅差分析中，formula公式的⾃变量（含协变量）顺序很重要，顺序很重要，顺序很重要！R中的计算效

应的顺序为序贯型，即如公式：y~A+B+A:B，R将评价1）A对y的影响，2）控制A，B对于y的影响，3）控制A和B的主效应，A与B的交

互效应。样本⼤⼩越不平衡，效应项的顺序对结果的影响就越⼤。越基础性的效应越需要放在前⾯，具体来说，就是协变量，然后是主效

应，再然后是双因素交互效应，再然后是三因素交互效应，再然后是四因素……基础性，⽬前我的理解就是变量的⽔平越简单，⽐如性别

（只有两个，三个也⾏）。直接引⽤R语⾔实战中的补充内容：《顺序很重要！》来解释⼀下顺序问题。

当⾃变量与其他⾃变量或者协变量相关时，没有明确的⽅法可以评价⾃变量对因变量的贡献。例如，含因⼦A、B和因变量y的双因素不平衡因⼦设计，有三

种效应：A和B的主效应，A和B的交互效应。假设你正使⽤如下表达式对数据进⾏建模：

Y~A+B+A:B

有三种类型的⽅法可以分解等式右边各效应对y所解释的⽅差。

类型I（序贯型）

效应根据表达式中先出现的效应做调整。A不做调整，B根据A调整，A:B交互项根据A和B调整。

类型II（分层型）

效应根据同⽔平或低⽔平的效应做调整。A根据B调整，B依据A调整，A:B交互项同时根据A和B调整。

类型III（边界型）

每个效应根据模型其他各效应做相应调整。A根据B和A:B做调整，A:B交互项根据A和B调整。

R默认调⽤类型I⽅法，其他软件（⽐如SAS和SPSS）默认调⽤类型III⽅法。

单因素⽅差分析

⾸先，我们要知道我们的数据结构，才可以使得aov来进⾏分析，以书中例⼦来看，它应该是属于我们第⼀章所说的融合后，即在数据

框中只有⼀个变量存放观测结果，其它变量均为因⼦向量，它们的组合唯⼀确定观测结果的值。

其次，使⽤aggregate(fit,by,FUN)来对数据集进⾏均值、⽅差等函数来进⾏初步的统计描述，得出初步结论。

第三，使⽤aov(formula,data)来进⾏⽅差分析，检验各个⽔平间是否显著差异，若p值⼩于显著⽔平（⼀般取0.05），则为各⽔平间

有显著差异，但是，⽅差分析函数aov并没有给出各个⽔平间的差异是否显著，所以需要继续分解。

第四，使⽤TukeyHSD(fit)函数（包含在基础包stats中）对数据进⾏多重⽐较，可以由结果直接得知，两两⽔平之间的显著差异

第五，作图，可以使⽤plot(TukeyHSD(fit))直接得出图表，微调标签这些这⾥不说了。得出的图中，若置信区间包含了0，则表⽰这

对数据不显著，书上也介绍了multcomp包中的glht函数，⽤法为glht(fit,linfct=mcp(trt=”Tukey”))，得到⼀个glht模型，使⽤

summary(glht)的话，可以得到⼀个与TukeyHSD⼀致的两两对⽐p值（参数linfct为线性假设的规范，mcp为多重对⽐的函数，返回⼀个

对照组，可以由特殊字符，如本例的“Tukey”，也可以由rbind(“rowName=c(1,0,-1))的组合给出，⽐如有c(L，M，H)，则c(-1,1,0)

为M-L，c(1,0,-1)为L-H，c(2,-1,-1)为L-M-H，向量和需为0），再使⽤cld返回⼀个cld模型，使⽤它作图函数为

plot(cld(glht,level=0.05)），可以作出⼀个头顶带着⼏个字母的箱线图，头顶的⼏个字母的意思是，若⼀对变量对应的字母包含相同字

母，则这⼀对变量的均值不显著差异。

第六，使⽤QQ图来验证因变量的正态性（注意QQ图需要使⽤回归模型），以及使⽤Bartlett检验

（(formula,data=datafram)）其⽅差齐性，检验的零假设为⽅差齐性（即若p值⼩于显著⽔平，则⽅差不齐）。

单因素协⽅差分析

与单因素⽅差分析步骤基本⼀致，只是⽐它多了⼀些协变量的存在。

第⼀，数据结构调整，观测数、协变量观测值均由其它标识列唯⼀确定，为⼀组数据。

第⼆，使⽤aggregate函数求均值，但其均值由于受协变量影响，此处均值的意义不⼤，此处可省略。

第三，使⽤aov(formula,data)进⾏协⽅差分析，把所有需要的协变量放到实际⽔平（或叫分组）变量前⾯，如

aov(y~coval1+coval2+…+covaln+val,data)。书中例⼦得到结果如下：

>fit<-aov(weight~gesttime+do,data=multcomp::litter)

>summary(fit)

DfSumSqMeanSqFvaluePr(>F)

gesttime1134.3134.308.0490.00597**

do3137.145.712.7390.04988*

Residuals691151.316.69

结果表明，gesttime与weight相关（p值显著），控制gesttime，do与weight相关，可以解释的是，我们没有办法⼈⼯控制

gesttime（在例⼦中，它是怀孕时间，没办法⼈⼯控制，所以把它作为⼀个协变量），但我们可以测量gesttime，以此作为协变量得到统

计上的控制，把gesttime作为⼀个统计控制量（协变量）来矫正其do分组对weight的均值影响，此时，其均值可以使⽤effects包中的

effect(do,fit)函数（在这⾥直接使⽤例⼦的do变量及得到的fit模型，帮助更好理解）。可以得到在⼈⼯控制变量do的各⽔平下，控

制了协变量之后的各组间均值。

第四，多重⽐较使⽤multcomp包中的glht函数，mcp参数由rbind(c(3,-1,-1,-1))给出，意思为⽤量为0的和其它三组⽤量不为0的作

⽐较，也可以使得c(1,-1,0,0)得到⽤量0和⽤量5的做⽐较。代码及结果如下：

>library(multcomp)

>contrast<-rbind("nodrugvsdrug"=c(3,-1,-1,-1))

>summary(glht(fit,linfct=mcp(do=contrast)))

MultipleComparisonsofMeans:Ur-definedContrasts

Fit:aov(formula=weight~gesttime+do)

LinearHypothes:

valuePr(>|t|)

nodrugvsdrug==08.2843.2092.5810.012*

上⾯的结果表明，⽤和不⽤药剂，其结果是有明显差异的。

第五，作图，可以使⽤HH包中的ancova（）函数来作图。函数式为：

ancova(weight~gesttime+do,data=litter)

需要⾃⾏加载HH包。图形如下，此项表明，在不同的剂量下，重量均随着怀孕时间的增加⽽增加，且各剂量的线相互平⾏（注意：各

直线平⾏是由于公式的设⽽⾮结果为平⾏）。

第六，假设验证，与单因素⽅差分析⼀样，均需要要求因变量的正态性及各⽔平的⽅差齐性，使⽤QQ图及来验证，另

外，协⽅差分析还需要对协变量与因变量的回归斜率相同（⽬的是为了验证协变量对因变量的影响是不受其它因素的变化⽽影响的），此项

可以直接使⽤⽅差分析中的交互项来判定，若交互项不显著，可以认为协变量与因变量是回归斜率相同，若交互项显著，则认为因变量与协

变量的关系依赖于⾃变量的⽔平：

>qqPlot(lm(weight~gesttime+do,data=litter))

>(weight~do,data=litter)

Bartlett'sK-squared=9.6497,df=3,p-value=0.02179

>summary(aov(weight~gesttime*do,data=litter))

DfSumSqMeanSqFvaluePr(>F)

gesttime1134.3134.308.2890.00537**

do3137.145.712.8210.04556*

gesttime:do381.927.291.6840.17889

Residuals661069.416.20

从结果上看，正态性及⽅差齐性均符合假设，回归斜率检验中，交互项不显著，可以说明回归斜率相同。

双因素⽅差分析

与单因素协⽅差分析类似，双因素⽅差分析的区别就在于它们有交互效应，其实协⽅差分析中后⾯的回归斜率检验就是双因素分析，⼀

个协变量的单因素协⽅差分析与双因素⽅差分析差异就在于，交互效应是否显著，当然从本质上不是这样的，不要混淆了。

>fit<-aov(len~supp*do,data=ToothGrowth)

>summary(fit)

DfSumSqMeanSqFvaluePr(>F)

supp1205.4205.412.3170.000894***

do12224.32224.3133.415<2e-16***

supp:do188.988.95.3330.024631*

Residuals56933.616.7

由此可以看到，使⽤的药物、剂量都对⽛齿的长度有影响，⽽且，药物和剂量的交互效应也⼗分显著。

R语⾔实战对双因素⽅差分析的交互效应绘图给出三个⽅法：

>#()函数，在基础包{stats}中

>(do,supp,len,type="b",col=c("blue","red"),pch=c(16,18),main="Interactionbetweendoandsupplementtype")

>#plotmeans()函数，在包gplots中

>plotmeans(len~interaction(supp,do,p=""),connect=list(c(1,3,5),c(2,4,6)))

>#interaction2wt()函数，在HH包中

>interaction2wt(len~supp*do)

对应的图如下：

重复测量⽅差分析

这⾥稍微说明⼀下什么叫重复测量⽅差分析，重复测量，测试对象会被重复测量其结果，打个⽐⽅就很容易明⽩了，⽐如，两个⼈在同

⼀年中各⽉的⾝⾼变化，这⾥就包含了⼀个因变量（⾝⾼），⼀个组间因⼦（⼈，有两个，所以是两个⽔平），⼀个组内因⼦（时间，12

个⽉，就是12个⽔平），此时，对这两个⼈来说（测试对象），他们⽣个⼈要测试12次，即组内因⼦的⽔平数，所以叫重复测量⽅差分

析。

但书中在分析关于CO2的例⼦时，使⽤了传统的重复测量⽅差分析。该⽅法假设任意组内因⼦的协⽅差矩阵为球形，并且任意组内因⼦

两⽔平间的⽅差之差都相等。但在现实中这种假设不可能满⾜，所以，本章只是说明⼀下理想情况下的分析，其它分析在书中未给出，后续

再介绍。

了解之后，回归书本上的例⼦，植被类型与⼆氧化碳浓度对⼆氧化碳吸收量的影响。若有两个或以上组内因⼦（假设为A,B,C,…），则

Error中代码为Error(Subject/(A+B+C+…))代码如下：

>w1b1<-subt(CO2,Treatment=="chilled")

>fit<-aov(uptake~Type*conc+Error(Plant/conc),data=w1b1)

>summary(fit)

Error:Plant

DfSumSqMeanSqFvaluePr(>F)

Type12667.22667.260.410.00148**

Residuals4176.644.1

Error:Plant:conc

DfSumSqMeanSqFvaluePr(>F)

conc1888.6888.6215.460.000125***

Type:conc1239.2239.258.010.001595**

Residuals416.54.1

Error:Within

DfSumSqMeanSqFvaluePr(>F)

Residuals3086928.97

>with(w1b1,(conc,Type,uptake,type="b",pch=c(16,18)))

>boxplot(uptake~Type*conc,data=w1b1,las=2)

对重复测量⽅差分析（混合模型设计，即既有组间因⼦，也有组内因⼦），书中给出了以下⼏种⽅法：

使⽤lme4包中的lmer()函数拟合线性混合模型（Bates，2005）

使⽤car包中的Anova()函数调整传统检验统计量以弥补球形假设的不满⾜（例如Geisr-Greenhou校正）

使⽤nlme包中的gls()函数拟合给定⽅差-协⽅差结构的⼴义最⼩⼆乘模型（UCLA，2009）

⽤多元⽅差分析对重复测量数据进⾏建模（Hand，1987）。

多元⽅差分析

使⽤manova函数，⾥⾯最主要的是，因变量需要预先使⽤cbind合并，并不是像formula⼀样直接使⽤y1+y2+y3~x这样的公式！正

确写法应该是manova(cbind(y1,y2,y3)~x1+x2,data=dataframe)，summary得到的结果为⾃变量是否对各个因变量有显著影响。若显

著，则可以使⽤(fit)进⾏进⼀步的分析，它可以显⽰各因变量与⾃变量的⽅差分析结果，即变为⼀元⽅差分析。

多元⽅差分析要求两个假设前提，多元正态性及⽅差-协⽅差矩阵同质性。多元正态性可以由QQ图检验，理论补充及代码如下：

#理论补充

#若有⼀个p*1的多元正态随机向量x，均值为µ，协⽅差矩阵为Σ，那么x与µ的马⽒距离的平⽅服从⾃由度为p的卡⽅分布。Q-Q图展⽰卡⽅分布的分位数，横纵

坐标分别是样本量与马⽒距离平⽅值。如果点全部落在斜率为1、截距项为0的直线上，则表明数据服从多元正态分布。

>center<-colMeans(y)

>n<-nrow(y)

>p<-ncol(y)

>cov<-cov(y)

>d<-mahalanobis(y,center,cov)

>coord<-qqplot(qchisq(ppoints(n),df=p),d,main="Q-QPlot",ylab="MaalanobisD2")

>abline(a=0,b=1)

>identify(coord$x,coord$y,labels=(UScereal))

⽅差-协⽅差矩阵同质性即指各组的协⽅差矩阵相同，R语⾔实战建议使⽤Box’sM检验来进⾏，但是R语⾔中却没有这个函数，⽽我

找资料的时候，看到了另⼀个⼈的说法，说Box’sM检验的准确性⾮常不好，使⽤的是F检验量，当样本量较⼤的时候，⽆论原数据怎么

样它都会显著，但也没给出别的检验办法，⽹上也暂时没找到合适的⽅法来进⾏检验，若后续找到了再补充。

多元离群点的检测，使⽤mvoutlier包中的(cbind(y1,y2,y3……))来进⾏，得出四张图，可以看到它已标记出对应的离群点。

稳健多元⽅差分析

在多元正态性或⽅差-协⽅差同质性或担⼼离群点时，可以使⽤稳健多元⽅差分析或⾮参数多元⽅差分析，稳健多元单因素⽅差分析可

以通过rrcov包中的(ame,x,method=”mcp”)，另外，vegan包中的adonis()函数则提供了⾮参数MANOVA的等同

形式，代码如下：

library(MASS)

attach(UScereal)

#稳健多元单因素⽅差分析

library(rrcov)

y<-cbind(calories,fat,sugars)

(y,shelf,method="mcd")#windows下运⾏超级超级超级慢!

#⾮参数多元⽅差分析

library(vegan)

adonis(y~shelf)

可以看出，这两个⽅差分析的结果都是显著的。即表明，货架的不同层对⾕物的三种物质的含量是有影响的。

⽤回归来做ANOVA

这章看得不是很明⽩，它的做法是，回归函数lm若碰到了⾃变量是因⼦（原要求因变量和⾃变量都是数值型）的时候，若有k个⽔平，

将会⾃动⽣成k-1个对照变量来代替⾃变量，然后进⾏回归拟合，代码如下：

>library(multcomp)

>levels(cholesterol$trt)

[1]"1time""2times""4times""drugD""drugE"

><-aov(respon~trt,data=cholesterol)

>summary()

DfSumSqMeanSqFvaluePr(>F)

trt41351.4337.832.439.82e-13***

Residuals45468.810.4

---

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

><-lm(respon~trt,data=cholesterol)

>summary()

Call:

lm(formula=respon~trt,data=cholesterol)

Residuals:

Min1QMedian3QMax

-6.5418-1.9672-0.00161.89016.6008

Coefficients:

valuePr(>|t|)

(Intercept)5.7821.0215.6659.78e-07***

trt2times3.4431.4432.3850.0213*

trt4times6.5931.4434.5683.82e-05***

trtdrugD9.5791.4436.6373.53e-08***

trtdrugE15.1661.44310.5071.08e-13***

---

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

Residualstandarderror:3.227on45degreesoffreedom

MultipleR-squared:0.7425,AdjustedR-squared:0.7196

F-statistic:32.43on4and45DF,p-value:9.819e-13

看到了么，⽅差分析所得的F检验的P值与回归得到的P值是⼀样的。⾄于原理，为什么会这样，我还不是很明⽩，也在找资料找⼈问当

中。

R中创建对照变量有五种⽅法，如下表：

创建⽅法描述

t第⼆个⽔平对照第⼀个⽔平，第三个⽔平对照前两个的均值，第四个⽔平对照前三个的均值，以此类推

基于正交多项式的对照，⽤于趋势分析（线性、⼆次、三次等）和等距⽔平的有序因⼦

对照变量之和限制为0。也称作偏差找对，对各⽔平的均值与所有⽔平的均值进⾏⽐较

ent各⽔平对照基线⽔平（默认第⼀个⽔平）。也称作虚拟编码

类似于ent，只是基线⽔平变成了最后⼀个⽔平。⽣成的系数类似于⼤部分SAS过程中使⽤的对照变量

使⽤<-lm(respon~trt,data=cholesterol,contrasts=””)的contrasts参数来改变创建⽅法，默认的创建⽅

法为ent)

以上是⽅差分析。