在第六节t检验部分的多组间差异分析提到涉及的内容比较多，并未记录，这次专门系统学习一下相关知识。Analysis of Variance，简称ANOVA

0、相关术语

• 方差分析（ANOVA）：多组别差异的分析；
• 单因素组间方差分析（one-way ANOVA）：实验设计为一个因素，两/多水平的组间因子。例如10位病人分为两组，每组5人，各接受一种治疗方案，评价组间区别。又称为单因素组间方差分析；

特点：一个病人只测量一次；实验分为几组就是几水平的。（较常见我认为）

• 单因素组内方差分析：实验设计为两或多水平的组内因子。例如评价一种治疗方案随时间的变化情况，10位病人为一组在所有时间水平下都进行测量。又因为每个病人不止被测了一次，也称为重复测量方差分析。

特点：一大组病人一起测；测量几次就是几水平的。

• 双因素方差分析：上述的治疗方案与时间均考虑时，即分析不同治疗方案的影响和治疗时间的影响，又分析治疗方案与时间的交互作用。前两个称为主效应，最后一个称为交互效应。将会做分别3次F检验：若疗法结果显著，则说明两种方案治疗效果不同；若时间结果显著，则说明时间长短会影响疗效；若交互效应影响显著，则说明两种方案随着时间变化治疗效果变化不同；
• 双因素混合模型方差分析：两因子一个是组内因子，一个是组间因子。
• 混淆因素：实验设计对象本身其它因素造成的对象起始水平就存在差异，从而也能一定程度解释因变量的组间差异的因素。因为对其不感兴趣，也可以称为干扰常数，可以进行调整；
• 协方差分析（ANCOVA）：约定好混淆因素为协变量，在对目标因素进行结果评价前，对任何混淆因素水平的组间差异进行的统计性调整的设计；
• 多元方差分析（MANOVA）：上述的例子实验结果只记录单个因变量y。当因变量不止一个（结果的多元分析）时，即为多元方差分析。若协变量也存在，即称为多元协方差分析。
• aov方差分析函数表达式中效应的顺序要严格注意：越基础性的效应越需要放在前面；一般首先是协变量--主效应--双因素交互项--三因素交互项

1、单因素方差分析

• 一般是指单因素组间方差分析
• 实验示例来自multcomp包的cholesterol的数据，记载了50位病人分为5组各接受一种trt后的胆固醇含量。

library(multcomp)  #第一次使用需要安装
cholesterol

cholesterol

1.1、数据合理性检验

单因素方差分析中，因变量要服从正态分布；并且各组方差相等。

（1）Q-Q图检验正态分布

library(car)
qqPlot(lm(response ~ trt, data=cholesterol),    #注意qqPlot() 函数要求用lm()拟合
       simulate=TRUE, main="Q-Q Plot", labels=FALSE)
#图形结果表示数据均落在95%的置信区间范围内，符合正态性假设

Q-Q Plot.png

（2）bartlett方差齐性检验

bartlett.test(response ~ trt, data=cholesterol)
#p值越大，表明5组的方差没有显著不同

bartlett.test.png

• 值得注意的是方差齐性分析对离群点分析非常敏感，利用car包的outlierTest()函数检测离群点

outlierTest(fit)
#结果表明 No Studentized residuals with Bonferroni

1.2、aov()单因素方差分析

attach(cholesterol)  
table(trt)     #查看水平分组以及各组样本大小（一般设计相同，为均衡设计）
aggregate(response, by=list(trt), FUN=mean)  #计算各组/水平平均值
aggregate(response, by=list(trt), FUN=sd)    #计算各组/水平方差
fit <- aov(response ~ trt)     #进行单因素方差分析 ，注意格式
summary(fit)

summary(fit)

结果表明ANOVA对trt（治疗方式）的F检验非常显著（p<0.0001）,说明五种疗法的效果不同

#将上述结果可视化
library(gplots)    #绘制带有置信区间的组均值图形
png("1.png")   #刚才在线显示有问题，存在本地显示正常
plotmeans(response ~ trt, xlab="Treatment", ylab="Response", 
          main="Mean Plot\nwith 95% CI")
#绘制组均值+置信区间图形
detach(cholesterol)
dev.off()

1.png

1.3、深入分析多重比较

• 在上述分析的基础上，进一步分析两组的差异显著与否，即具体哪种疗法与其它疗法不同。有以下两种方法--

#法1
TukeyHSD(fit)  # p adj 值越小，越显著；即该两组均值差异显著
#将上述数据可视化
opar <- par(no.readonly = TRUE)
par(las=2)      # 旋转轴标签（主要是便于纵坐标便签展示）
par(mar=c(5,8,4,2))   #增大左边界的面积（也是针对纵坐标）
plot(TukeyHSD(fit)) 
#凡是置信区间包含0的（靠左边的）说明差异不显著（p>0.5）
par(opar)

TukeyHSD.png

#法2
library(multcomp)
opar <- par(no.readonly = TRUE)
par(mar=c(5,4,6,2))   #扩大顶部边界面积
tuk <- glht(fit, linfct=mcp(trt="Tukey"))   
plot(cld(tuk, level=.05),col="lightgrey") 
#cld()函数中的level选项设置了使用的显著水平（0.05，即本例中的95%的置信区间）
par(opar)   

如下图，图形上方中，有相同字母的组说明均值差异不显著。该图此外也能反映各组的分布情况，因此信息更多。

箱图.png

2、单因素协方差分析（ANCOVA）

• litter为一实验数据框，列名为剂量、体重，怀孕时间（最后一列应该是产崽数量，不关注）；目的为研究怀孕小鼠接受不同剂量药物产崽体重的变化情况。怀孕时间为协变量。

  
  litter

2.1、数据合理性检验

• ANCOVA也要进行正态分布与齐方差检验，步骤同前，不再赘述；

同质性检验

• 此外还要进行协变量的回归率同质性检验，协变量影响与剂量主因素不存在显著的交互作用，即怀孕时间对体重的影响不会因剂量分组的不同而不同。

library(multcomp)
data(litter, package="multcomp") 
fit <- aov(weight ~ gesttime*dose, data=litter)
summary(fit)

注意，上述表达式 y ~ A*B 为双因素方差分析表达式。我认为这里只是为了单纯看下交互项是否显著，而使用这个表达式的。

summary(fit)

• 返回结果表明交互效应不显著，支持了斜率相等的假设。此外也可绘制因变量、协变量与因子之间的关系图。

library(HH)
ancova(weight ~ gesttime + dose, data=litter)
#可以清楚看出回归线相互平行“同质性”，而截距项不同（剂量因素不同导致）

同质性检验.png

2.2、ANCOVA分析

• 这里强调下表达式顺序，之前为单因素分析仅一个因子。当存在协变量时，协变量要放在前面，其次是主效应，再是双因素交互项

attach(litter)
table(dose)  #按实验剂量分为四组，发现每组的数量都不同（为非均衡设计）
aggregate(weight, by=list(dose), FUN=mean)
fit2 <- aov(weight ~ gesttime + dose)   # 协变量+单因素  注意顺序！
summary(fit2)
# F检验结果表明怀孕时间域产崽体重相关；控制怀孕时间（协变量），药物剂量与出生体重有关。

平均值1 summary(fit2)

2.3、衍生计算

（1）获得去除协变量效应的组均值

我的理解是将小鼠的怀孕时间都变换为相同值时后的剂量影响值

library(effects)
effect("dose", fit2)
#可以对比一下之间用aggregate()计算得到的平均值

平均值2

（2）化繁为简实验目的

• 研究用药与否对产崽体重的影响，即比较第一组（剂量为0）与第二~四组（剂量为5、50,500）的差异

library(multcomp)
contrast <- rbind("no drug vs. drug" = c(3, -1, -1, -1))  
#设定为第一组与其它三组的均值比较
summary(glht(fit2, linfct=mcp(dose=contrast)))
#假设检验的t统计量（2.581）在p<0.05水平下显著，可以得出未用药比用药组体重高的结论。

no drug vs. drug

3、双因素方差分析

• 示例实验来自基础数据ToothGrowth，为研究两种试剂（supp）的各自三种剂量水平（dose）对豚鼠牙齿长度的影响（每个豚鼠只测一次，为组间因子）
  ToothGrowth
  注意：双因素方差分析时，要保证两组样本量相同；否则y ~ A:B与y ~ B:A两个模型结果会有差异。

attach(ToothGrowth)
aggregate(len, by=list(supp,dose), FUN=mean)    #共3*2=6组
aggregate(len, by=list(supp,dose), FUN=sd)
dose <- factor(dose)  
#ToothGrowth数据框中supp已经默认为因子了，这样aov() 函数就会将它当做一个分组变量，而不是一个数值型协变量.
fit <- aov(len ~ supp*dose)  #两因素的函数表达格式
summary(fit)  
#结果表明主效应与交互效应都非常显著

summary(fit)

• 可视化一下~

#法1
interaction.plot(dose, supp, len, type="b", 
                 col=c("red","blue"), pch=c(16, 18),
                 main = "Interaction between Dose and Supplement Type")
#法2
library(gplots)
plotmeans(len ~ interaction(supp, dose, sep=" "),
          connect=list(c(1, 3, 5),c(2, 4, 6)), 
          col=c("red","darkgreen"),
          main = "Interaction Plot with 95% CIs", 
          xlab="Treatment and Dose Combination")
#观察与上图的区别（显示要问题的话配合 png("2.png")  dev.off()下载到本地）
#法3
library(HH)
interaction2wt(len~supp*dose)
#显示为四幅图的汇总大图；箱线图为控制一种水平下，另一种因素不同水平的评价；折线图为交互效应互相影响变化图。

1.png
2.png
3.png

4、重复测量方差分析（组内因子）

• 常见的实验设计为含一个组内因子，一个组间因子重复测量方差分析。比如本节一开始的论述两种治疗方法随时间长短变化的影响。
• 这里示例实验来自基础函数CO2，此处我们去研究两种植物（Type，为组间因子）在不同的CO2浓度范围（conc，为组内因子，重复测量）下二氧化碳吸收量（uptake）的变化。

CO2$conc <- factor(CO2$conc)
w1b1 <- subset(CO2, Treatment=='chilled')  #抽取目标行数据
fit <- aov(uptake ~ (conc*Type) + Error(Plant/(conc)), w1b1) 
#注意含单个组间因子（Type）、单个组内因子（conc）的重复测量方差分析的表达格式。
#(conc*Type)里两者的顺序可以颠倒；
#Error(Plant/(conc))中Plant是相对于每组conc中的独有的标识变量。即重复测量的单位
summary(fit)
# 结果表明在0.01的水平下，主效应type与conc以及交互效应都非常显著

summary(fit)

• 可视化一下~

par(las=2)
par(mar=c(10,4,4,2))
#法1：点线图
with(w1b1, 
     interaction.plot(conc,Type,uptake, 
                      type="b", col=c("red","blue"), pch=c(16,18),
                      main="Interaction Plot for Plant Type and Concentration"))
#法2：箱线图
boxplot(uptake ~ Type*conc, data=w1b1, col=(c("gold","green")),
        main="Chilled Quebec and Mississippi Plants", 
        ylab="Carbon dioxide uptake rate (umol/m^2 sec)")
par(opar)

1.png
2.png

结合上图及分析结果可得结论：两种植物吸收二氧化碳能力差异显著，且随二氧化碳浓度升高，差异越明显。

5、多元方差分析（MANOVA）

• 因变量（结果变量）不止一个时（y1，y2...），可用多元方差对它们同时进行分析
• 示例实验设计为不同储架层谷物（自变量）里，卡路里calories、脂肪fat以及糖sugars含量（3个因变量）是否存在差异

5.1、多元正态性检验

library(MASS)
attach(UScereal)    #数据包
shelf <- factor(shelf) #
y <- cbind(calories, fat, sugars)  #将3个因变量合并为一个目标矩阵
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="QQ Plot Assessing Multivariate Normality",
                ylab="Mahalanobis D2")
abline(a=0,b=1)  
#绘制y=x 的参考线；若数据服从多元正态分布，则点将落在直线上。
#比如本图中Wheaties Honey Gold点与Wheaties点异常，违反多元正态性。可以删除这两个点再重新分析。
identify(coord$x, coord$y, labels=row.names(UScereal))  #交互式标注重点关注的点

Q-Q Plot

5.2、MANOVA分析

aggregate(y, by=list(shelf), FUN=mean)  #获取各个货架的各个均值
cov(y)  #输出个谷物间的方差和协方差（是什么？存疑）
fit <- manova(y ~ shelf)  #多元方差分析
summary(fit) 
#返回结果说明三个货架的营养成分测量值不同
summary.aov(fit)
#对每一个因变量做单因素方差分析，返回结果表明三个货架中每种营养成分的测量值都是不同的

summary(fit)
summary.aov(fit)

5.3、稳健单因素MANOVA

• 在存在不满足多元正态性、方差-协方差均值假设，以及多元离群点，可以考虑使用本方法--稳健或非参数版本的MANOVA检验

library(rrcov)
Wilks.test(y,shelf, method="mcd")   #计算过程需要等几十秒时间（也是第一次遇到）~ 

稳健MANOVA

• 返回结果表明稳健检验对离群点和违反MANOVA假设的情况不敏感，也验证了存储在不同层货架的营养含量不同。

教材中还提到一种用回归方法做ANOVA的介绍（P220）
参考书籍《R语言实战（第2版）》