第9章 方差分析

注：R语言的再复习之路
 
 

1. ANOVA模型拟合

1.1 aov()函数

语法为aov(formula, data = dataframe)

符号	用法
~	分隔符号，左边为响应变量，右边为解释变量
:	表示变量的交互项，例如y ~ A + B + C
*	表示所有可能交互项，例如y ~ A * B * C等价于y ~ A + B + C + A:B + B:C + A:C + A:B:C
^	表示交互项达到某个次数，例如y ~ (A + B + C)^2等价于y ~ A + B + C + A:B + B:C + A:C
.	表示包含除因变量外的所有变量，例如y ~ .等价于y ~ A + B + C

设计	表达式
单因素ANOVA	y ~ A
含单个协变量的单因素ANCOVA	y ~ x + A
双因素ANOVA	y ~ A * B
含两个协变量的双因素ANCOVA	y ~ x1 + x2 + A * B
单因素组内ANOVA	y ~ A + Error(Subject/A)
含单个组内因子(W)和单个组间因子(B)的重复测量ANOVA	y ~ B * W + Error(Subject/W)

1.2 表达式中各项的顺序

表达式中效应的顺序在两种情况下会造成影响：（a）因子不只一个，并且是非均衡设计；（b）存在协变量。
表达式的抒写规则：首先是协变量，然后是主效应，接着是双因素的交互项，再接着是三因素的交互项，以此类推。
 
 

2. 单因素方差分析

library(multcomp)
attach(cholesterol)

# 基本的统计分析
table(trt)
aggregate(response, by = list(trt), FUN = mean)
aggregate(response, by = list(trt), FUN = sd)

# 单因素方差分析
fit <- aov(response ~ trt)
summary(fit)

# 绘制带有置信区间的组均值图形
library(gplots)
plotmeans(response ~ trt, xlab = 'Treatment', ylab = 'Response', main = 'Mean Plot \n with 95% CI')
detach(cholesterol)

2.1 多重比较

1. TukeyHSD()函数

TukeyHSD(fit)
par(las = 2)
par(mar = c(2, 2))
plot(TukeyHSD(fit))

若图形中的置信区间包含0，则说明对应的两个水平之间差异不显著
2. glht()函数

library(multcomp)
par(mar = c(5, 4, 6, 2))
tuk <- glht(fit, linfct = mcp(trt = 'Tukey'))
plot(cld(tuk, level = .05), col = 'lightgrey')

2.2 评估检验的假设条件

1. 因变量服从正态分布（正态性检验）

library(car)
qqPlot(lm(response ~ trt, data = cholesterol), simulate = TRUE, main = 'Q-Q Plot', labels = FALSE)

2. 各组方差相等（方差齐性检验）

barlett.test(response ~ trt, data = cholesterol)

# 注：方差齐性分析对离群点非常敏感，需要检测离群点
library(car)
outlierTest(fit)

 
 

3. 单因素协方差分析

data(litter, package = 'multcomp')
attach(litter)
fit <- aov(weight ~ gesttime + dose)
summary(fit)

由于使用了协变量，你可能想要获取调整的组均值，即去除协变量效应后的组均值，代码如下：

library(effects)
effect('dose', fit)

3.1 多重比较

library(multcomp)
contrast <- rbind('no drug vs. drug' = c(3, -1, -1, -1))
summary(glht(fit, linfct = mcp(dose = contrast)))

3.2 评估检验的假设条件

1. 正态性检验
   同上
2. 方差齐性检验
   同上
3. 检验回归斜率的同质性
   ANCOVA还假定回归斜率相同，原假设假定各组的回归斜率相等。

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

3.3 结果可视化

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

若图形上的直线平行，则表明回归的斜率相等。
 
 

4. 双因素方差分析

attach(ToothGrowth)
dose <- factor(dose)
fit <- aov(len ~ supp * dose)
summary(fit)
detach(ToothGrowth)

注：这里需要将dose变量转换为因子变量，这样aov()函数才会将它当做一个分组变量，否则函数会把它算作数值型协变量，会发生错误。
注：所以以后用R语言做方差分析，一定要注意将分组变量因子化。

4.1 结果可视化

1. interaction.plot()函数

interaction.plot(dose, supp, len, type = 'b', col = c('red', 'blue'), pch = c(16, 18), main = 'Interaction between Dose and Supplement Type')

注：前三个参数分别为x变量（这里指dose），分组变量（这里指supp），y变量（这里指len）

2. plotmeans()函数

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')

3. interaction2wt()函数

library(HH)
interaction2wt(len ~ supp*dose)

 
 

5. 重复测量方差分析

# 构建模型
CO2$conc <- factor(CO2$conc)
wlbl <- subset(CO2, Treatment == 'chilled')
fit <- aov(uptake ~ conc*Type + Error(Plant/conc), wlbl)
summary(fit)

# 可视化
par(las = 2)
par(mar = c(10, 4, 4, 2))
with(wlbl, interaction.plot(conc, Type, uptake, type = 'b', col = c('red', 'blue'), pch = c(16, 18), main = 'Interaction Plot for Plant Type and Concentration'))
boxplot(uptake ~ Type*conc, data = wlbl, col = c('gold', 'green'), main = 'Chilled Quebec and Mississippi Plants', ylab = 'Carbon dioxide uptake rate (umol/m^2 sec)')