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检验（bartlett.test(formula,data=datafram)）其方差齐性，检验的零假设为方差齐性（即若p值小于显著水平，则方差不齐）。
　　
单因素协方差分析
　　与单因素方差分析步骤基本一致，只是比它多了一些协变量的存在。
　　第一，数据结构调整，观测数、协变量观测值均由其它标识列唯一确定，为一组数据。
　　第二，使用aggregate函数求均值，但其均值由于受协变量影响，此处均值的意义不大，此处可省略。
　　第三，使用aov(formula,data)进行协方差分析，把所有需要的协变量放到实际水平（或叫分组）变量前面，如aov(y~coval1+coval2+…+covaln+val,data)。书中例子得到结果如下：

> fit <- aov(weight~gesttime+dose,data=multcomp::litter)
> summary(fit)
Df Sum Sq Mean Sq F value  Pr(>F)
gesttime     1  134.3  134.30   8.049 0.00597 **
dose         3  137.1   45.71   2.739 0.04988 *
Residuals   69 1151.3   16.69

　　结果表明，gesttime与weight相关（p值显著），控制gesttime，dose与weight相关，可以解释的是，我们没有办法人工控制gesttime（在例子中，它是怀孕时间，没办法人工控制，所以把它作为一个协变量），但我们可以测量gesttime，以此作为协变量得到统计上的控制，把gesttime作为一个统计控制量（协变量）来矫正其dose分组对weight的均值影响，此时，其均值可以使用effects包中的effect(dose,fit)函数（在这里直接使用例子的dose变量及得到的fit模型，帮助更好理解）。可以得到在人工控制变量dose的各水平下，控制了协变量之后的各组间均值。
　　第四，多重比较使用multcomp包中的glht函数，mcp参数由rbind(c(3,-1,-1,-1))给出，意思为用量为0的和其它三组用量不为0的作比较，也可以使得c(1,-1,0,0)得到用量0和用量5的做比较。代码及结果如下：

> library(multcomp)
> contrast<-rbind("no drug vs drug"=c(3,-1,-1,-1))
> summary(glht(fit,linfct=mcp(dose=contrast)))
Multiple Comparisons of Means: User-defined Contrasts
Fit: aov(formula = weight ~ gesttime + dose)
Linear Hypotheses:Estimate Std. Error t value Pr(>|t|)
no drug vs drug == 0    8.284      3.209   2.581    0.012 *

　　上面的结果表明，用和不用药剂，其结果是有明显差异的。
　　第五，作图，可以使用HH包中的ancova（）函数来作图。函数式为：

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

　　需要自行加载HH包。图形如下，此项表明，在不同的剂量下，重量均随着怀孕时间的增加而增加，且各剂量的线相互平行（注意：各直线平行是由于公式的设而非结果为平行）。
　　
　　第六，假设验证，与单因素方差分析一样，均需要要求因变量的正态性及各水平的方差齐性，使用QQ图及bartlett.Test来验证，另外，协方差分析还需要对协变量与因变量的回归斜率相同（目的是为了验证协变量对因变量的影响是不受其它因素的变化而影响的），此项可以直接使用方差分析中的交互项来判定，若交互项不显著，可以认为协变量与因变量是回归斜率相同，若交互项显著，则认为因变量与协变量的关系依赖于自变量的水平：

> qqPlot(lm(weight~gesttime+dose,data=litter))
> bartlett.test(weight~dose,data=litter)
Bartlett's K-squared = 9.6497, df = 3, p-value = 0.02179
> summary(aov(weight~gesttime*dose,data=litter))Df Sum Sq Mean Sq F value  Pr(>F)
gesttime       1  134.3  134.30   8.289 0.00537 **
dose           3  137.1   45.71   2.821 0.04556 *
gesttime:dose  3   81.9   27.29   1.684 0.17889
Residuals     66 1069.4   16.20

　　从结果上看，正态性及方差齐性均符合假设，回归斜率检验中，交互项不显著，可以说明回归斜率相同。
　　
双因素方差分析
　　与单因素协方差分析类似，双因素方差分析的区别就在于它们有交互效应，其实协方差分析中后面的回归斜率检验就是双因素分析，一个协变量的单因素协方差分析与双因素方差分析差异就在于，交互效应是否显著，当然从本质上不是这样的，不要混淆了。

> fit<-aov(len~supp*dose,data=ToothGrowth)
> summary(fit)Df Sum Sq Mean Sq F value   Pr(>F)
supp         1  205.4   205.4  12.317 0.000894 ***
dose         1 2224.3  2224.3 133.415  < 2e-16 ***
supp:dose    1   88.9    88.9   5.333 0.024631 *
Residuals   56  933.6    16.7

　　由此可以看到，使用的药物、剂量都对牙齿的长度有影响，而且，药物和剂量的交互效应也十分显著。
　　R语言实战对双因素方差分析的交互效应绘图给出三个方法：

> #interaction.plot()函数，在基础包{stats}中
> interaction.plot(dose,supp,len,type="b",col=c("blue","red"),pch=c(16,18),main="Interaction between dose and supplement type")
> #plotmeans()函数，在包gplots中
> plotmeans(len~interaction(supp,dose,sep=" "),connect=list(c(1,3,5),c(2,4,6)))
> #interaction2wt()函数，在HH包中
> interaction2wt(len~supp*dose)

　　对应的图如下：
　　
　　

重复测量方差分析
　　这里稍微说明一下什么叫重复测量方差分析，重复测量，测试对象会被重复测量其结果，打个比方就很容易明白了，比如，两个人在同一年中各月的身高变化，这里就包含了一个因变量（身高），一个组间因子（人，有两个，所以是两个水平），一个组内因子（时间，12个月，就是12个水平），此时，对这两个人来说（测试对象），他们生个人要测试12次，即组内因子的水平数，所以叫重复测量方差分析。
　　但书中在分析关于CO2的例子时，使用了传统的重复测量方差分析。该方法假设任意组内因子的协方差矩阵为球形，并且任意组内因子两水平间的方差之差都相等。但在现实中这种假设不可能满足，所以，本章只是说明一下理想情况下的分析，其它分析在书中未给出，后续再介绍。
　　了解之后，回归书本上的例子，植被类型与二氧化碳浓度对二氧化碳吸收量的影响。若有两个或以上组内因子（假设为A,B,C,…），则Error中代码为Error(Subject/(A+B+C+…))代码如下：

> w1b1<-subset(CO2,Treatment=="chilled")
> fit<-aov(uptake~Type*conc+Error(Plant/conc),data=w1b1)
> summary(fit)
Error: PlantDf Sum Sq Mean Sq F value  Pr(>F)
Type       1 2667.2  2667.2   60.41 0.00148 **
Residuals  4  176.6    44.1
Error: Plant:concDf Sum Sq Mean Sq F value   Pr(>F)
conc       1  888.6   888.6  215.46 0.000125 ***
Type:conc  1  239.2   239.2   58.01 0.001595 **
Residuals  4   16.5     4.1
Error: WithinDf Sum Sq Mean Sq F value Pr(>F)
Residuals 30    869   28.97
> with(w1b1,interaction.plot(conc,Type,uptake,type="b",pch=c(16,18)))
> boxplot(uptake~Type*conc,data=w1b1,las=2)

　　
　　对重复测量方差分析（混合模型设计，即既有组间因子，也有组内因子），书中给出了以下几种方法：
　　使用lme4包中的lmer()函数拟合线性混合模型（Bates，2005）
　　使用car包中的Anova()函数调整传统检验统计量以弥补球形假设的不满足（例如Geisser-Greenhouse校正）
　　使用nlme包中的gls()函数拟合给定方差-协方差结构的广义最小二乘模型（UCLA，2009）
　　用多元方差分析对重复测量数据进行建模（Hand，1987）。
　　
多元方差分析
　　使用manova函数，里面最主要的是，因变量需要预先使用cbind合并，并不是像formula一样直接使用y1+y2+y3~x这样的公式！正确写法应该是manova(cbind(y1,y2,y3)~x1+x2,data=dataframe)，summary得到的结果为自变量是否对各个因变量有显著影响。若显著，则可以使用summary.aov(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-Q Plot",ylab="Maalanobis D2")
> abline(a=0,b=1)
> identify(coord$x,coord$y,labels=row.names(UScereal))

　　方差-协方差矩阵同质性即指各组的协方差矩阵相同，R语言实战建议使用Box’s M检验来进行，但是R语言中却没有这个函数，而我找资料的时候，看到了另一个人的说法，说Box’s M检验的准确性非常不好，使用的是F检验量，当样本量较大的时候，无论原数据怎么样它都会显著，但也没给出别的检验办法，网上也暂时没找到合适的方法来进行检验，若后续找到了再补充。
　　多元离群点的检测，使用mvoutlier包中的aq.plot(cbind(y1,y2,y3……))来进行，得出四张图，可以看到它已标记出对应的离群点。
　　
稳健多元方差分析
　　在多元正态性或方差-协方差同质性或担心离群点时，可以使用稳健多元方差分析或非参数多元方差分析，稳健多元单因素方差分析可以通过rrcov包中的Wilks.test(y.dataframe,x,method=”mcp”)，另外，vegan包中的adonis()函数则提供了非参数MANOVA的等同形式，代码如下：

library(MASS)
attach(UScereal)
#稳健多元单因素方差分析
library(rrcov)
y<-cbind(calories,fat,sugars)
Wilks.test(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"
> fit.aov <- aov(response ~ trt, data=cholesterol)
> summary(fit.aov)Df Sum Sq Mean Sq F value   Pr(>F)
trt          4 1351.4   337.8   32.43 9.82e-13 ***
Residuals   45  468.8    10.4
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> fit.lm <- lm(response ~ trt, data=cholesterol)
> summary(fit.lm)
Call:
lm(formula = response ~ trt, data = cholesterol)
Residuals:
Min      1Q  Median      3Q     Max
-6.5418 -1.9672 -0.0016  1.8901  6.6008
Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept)    5.782      1.021   5.665 9.78e-07 ***
trt2times      3.443      1.443   2.385   0.0213 *
trt4times      6.593      1.443   4.568 3.82e-05 ***
trtdrugD       9.579      1.443   6.637 3.53e-08 ***
trtdrugE      15.166      1.443  10.507 1.08e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.227 on 45 degrees of freedom
Multiple R-squared:  0.7425,    Adjusted R-squared:  0.7196
F-statistic: 32.43 on 4 and 45 DF,  p-value: 9.819e-13

　　看到了么，方差分析所得的F检验的P值与回归得到的P值是一样的。至于原理，为什么会这样，我还不是很明白，也在找资料找人问当中。
　　R中创建对照变量有五种方法，如下表：

创建方法	描述
contr.helmert	第二个水平对照第一个水平，第三个水平对照前两个的均值，第四个水平对照前三个的均值，以此类推
contr.poly	基于正交多项式的对照，用于趋势分析（线性、二次、三次等）和等距水平的有序因子
contr.sum	对照变量之和限制为0。也称作偏差找对，对各水平的均值与所有水平的均值进行比较
contr.treatment	各水平对照基线水平（默认第一个水平）。也称作虚拟编码
contr.SAS	类似于contr.treatment，只是基线水平变成了最后一个水平。生成的系数类似于大部分SAS过程中使用的对照变量

　　使用fit.lm <- lm(response ~ trt, data=cholesterol, contrasts=”contr.SAS”)的contrasts参数来改变创建方法，默认的创建方法为contr.treatment)
　　以上是方差分析。
　　
　　
　　
　　