使用R语言进行单（双）因素方差分析

9.1 先了解一些术语

方差分析： ANOVA，协方差分析：ANCOVA (analysis of covariance )，多元方差分析： muti-variate ANOVA，多元协方差分析：MANCOVA

9.2 ANOVA 模型拟合

R语言中使用 aov() 函数拟合 ANOVA模型，结果和回归模型中 lm() 函数一样

aov()函数

R 表达式中常见的特殊符号：

R 常见研究设计的表达式：

单因素ANOVA： y ~ A
协变量的单因素ANOVA： y ~ x + A
双因素： y ~ A + B

注意：因子的顺序很重要

9.3 one-way ANOVA

library(multcomp)
attach(cholesterol) # 加载数据
table(trt)
cholesterol
aggregate(response, by = list(trt), FUN = mean) # 求各组均值
aggregate(response, by = list(trt), FUN = sd)# 求标准差fit <- aov(response ~ trt) # 检验组间差异
summary(fit)
install.packages('gplots')
library(gplots)
plotmeans(response ~ trt, xlab = 'treatmeat', ylab = 'response', main = 'mean plot with 95% CI')
detach(cholesterol) # detach 数据

9.3.1 多重比较 – 就是组内再比较

使用TukeyHSD() 函数可以实现多重比较

TukeyHSD(fit) #多重比较
par(las = 2) # 旋转轴标签
par(mar = c(5,8,4,2)) # mar设置下，左，上，右的边界
plot(TukeyHSD(fit))

说明：图中置信区间包含0的代表不显著

multcomp包中的glht()函数更好进行多重比较

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

9.3.2 评价检验的假设条件 – 使用car包中的QQ图

library(car)
qqPlot(lm(response ~ trt, data = cholesterol),simulate = T, main = 'qq plot', labels = T)
# 注意： 这个qqPlot的P大写

点落在95% CI说明符合正态分布

9.4 单因素协方差ANCOAVA分析 –

**使用multcomp包中的litter数据 – **dose是自变量， weight是因变量，gesttime是协变量

data(litter, package = 'multcomp')
attach(litter)
table(dose)
aggregate(weight, by = list(dose), FUN = mean)# 按dose求weight的mean
fit <- aov(weight ~ gesttime + dose)## 注意协变量放在前面
summary(fit)

同ANOVA一样，随后进行多重比较

library(multcomp)
contrast <- rbind('no drug vs. drug' = c(3,-1,-1,-1))
summary(glht(fit, linfct = mcp(dose = contrast)))
#contrast 定义 dose的第一组和其他三组比较
# 结果显示 第一组与其余三组的mean有显著差异

9.4.1 评价检验的假设条件

因为ANOVA和ANCOVA都需要数据服从正态分布，所以要检验假设条件

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

可以看到，gesttime :dose 不显著，表明时间和幼崽的出生体重的关系不依赖于药物剂量的水平

9.4.2 visualize results

使用HH包中的ancova可视化

install.packages('HH')
library(HH)
ancova(weight ~ gesttime + dose, data = litter)

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

前者斜率相同，后者斜率不同

9.5 双因素方差分析

很重要我的实验就是这个

ToothGrowth# R中base安装的数据：两种饮食下各三种处理， 2x3的实验设计

attach(ToothGrowth)
The following object is masked from litter:dose> table(supp, dose)dose
supp 0.5  1  2OJ  10 10 10VC  10 10 10
> aggregate(len, by = list(supp, dose), FUN = mean)Group.1 Group.2     x
1      OJ     0.5 13.23
2      VC     0.5  7.98
3      OJ     1.0 22.70
4      VC     1.0 16.77
5      OJ     2.0 26.06
6      VC     2.0 26.14
> aggregate(len, by = list(supp, dose), FUN = sd)Group.1 Group.2        x
1      OJ     0.5 4.459709
2      VC     0.5 2.746634
3      OJ     1.0 3.910953
4      VC     1.0 2.515309
5      OJ     2.0 2.655058
6      VC     2.0 4.797731
> dose <- factor(ToothGrowth$dose)
> fit <- aov(len ~ supp * dose)
> summary(fit)Df Sum Sq Mean Sq F value   Pr(>F)
supp         1  205.4   205.4  15.572 0.000231 ***
dose         2 2426.4  1213.2  92.000  < 2e-16 ***
supp:dose    2  108.3    54.2   4.107 0.021860 *
Residuals   54  712.1    13.2
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

从结果可以看出，主效应和交互效应都非常显著

下一步就是可视化了 – 有三种可视化方法

使用interaction.plot()
gplots包中的plotmeans()函数
HH包中的interaction2wt()函数

使用interaction.plot()

interaction.plot(dose, supp, len, type = 'b',col = c('red', 'blue'), pch = c(12, 18),main = 'interaction between supp and dose')

gplots包中的plotmeans()函数

library(gplots)
plotmeans(len ~ interaction(supp, dose, sep = ' '),connect = list(c(1,3,5),c(2,4,6)),col = c('red', 'green'),main = 'interaction plot with 95% CI',xlab = 'treatment and dose combination')
## connect()参数  定义 X轴的位置

HH包中的interaction2wt()

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

**可以看出HH包中的interaction2wt函数画出的图最好用，包含的信息更多 **