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检验
新年口号（st(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)
Df Sum Sq Mean Sq F value Pr(>F)
gesttime 1 134.3 134.30 8.049 0.00597 **
do 3 137.1 45.71 2.739 0.04988 *
Residuals 69 1151.3 16.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("no drug vs drug"=c(3,-1,-1,-1))
> summary(glht(fit,linfct=mcp(do=contrast)))
Multiple Comparisons of Means: Ur-defined Contrasts
Fit: aov(formula = weight ~ gesttime + do)
Linear Hypothes:
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+do,data=litter)
需要⾃⾏加载HH包。图形如下，此项表明，在不同的剂量下，重量均随着怀孕时间的增加⽽增加，且各剂量的线相互平⾏（注意：各直线平⾏是由于公式的设⽽⾮结果为平⾏）。
第六，假设验证，与单因素⽅差分析⼀样，均需要要求因变量的正态性及各⽔平的⽅差齐性，使⽤QQ图及bartlett.Test来验证，另外，协⽅差分析还需要对协变量与因变量的回归斜率相同（⽬的是为了验证协变量对因变量的影响是不受其它因素的变化⽽影响的），此项可以直接使⽤⽅差分析中的交互项来判定，若交互项不显著，可以认为协变量与因变量是回归斜率相同，若交互项显著，则认为因变量与协变量的关系依赖于⾃变量的⽔平：
> qqPlot(lm(weight~gesttime+do,data=litter))
> st(weight~do,data=litter)
Bartlett's K-squared = 9.6497, df = 3, p-value = 0.02179
> summary(aov(weight~gesttime*do,data=litter))
Df Sum Sq Mean Sq F value Pr(>F)
gesttime 1 134.3 134.30 8.289 0.00537 **
do 3 137.1 45.71 2.821 0.04556 *
gesttime:do 3 81.9 27.29 1.684 0.17889
Residuals 66 1069.4 16.20
从结果上看，正态性及⽅差齐性均符合假设，回归斜率检验中，交互项不显著，可以说明回归斜率相同。
双因素⽅差分析
与单因素协⽅差分析类似，双因素⽅差分析的区别就在于它们有交互效应，其实协⽅差分析中后⾯的回归斜率检验就是双因素分析，⼀个协变量的单因素协⽅差分析与双因素⽅差分析差异就在于，交互效应是否显著，当然从本质上不是这样的，不要混淆了。
> fit<-aov(len~supp*do,data=ToothGrowth)
> summary(fit)
Df Sum Sq Mean Sq F value Pr(>F)一百的英文
supp 1 205.4 205.4 12.317 0.000894 ***
do 1 2224.3 2224.3 133.415 < 2e-16 ***
supp:do 1 88.9 88.9 5.333 0.024631 *
Residuals 56 933.6 16.7
由此可以看到，使⽤的药物、剂量都对⽛齿的长度有影响，⽽且，药物和剂量的交互效应也⼗分显著。
箭头怎么打 R语⾔实战对双因素⽅差分析的交互效应绘图给出三个⽅法：
蜂蜜和什么不能一起吃> #interaction.plot()函数，在基础包{stats}中
> interaction.plot(do,supp,len,type="b",col=c("blue","red"),pch=c(16,18),main="Interaction between do and supplement type") > #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+…))代码如下：
