参考资料:R语言实战【第2版】
单因素方差分析中,沃恩感兴趣的是比较分类因子定义的两个或多个组别中的因变量均值。以multcomp包中的cholesterol数据集为例,进行分析:
# 加载multcomp包
library(multcomp)
# 查看数据
head(cholesterol)
attach(cholesterol)
# 查看处理类型
table(trt)
# 查看不同处理的均值
aggregate(response~trt,data=cholesterol,FUN=mean)
# 或者
aggregate(response,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\nwith 95% CI")
detach(cholesterol)
1、多重比较
上面的方差分析告诉我们这5种疗法的效果存在显著差异,但并没有告诉我们那种疗法与其他疗法存在不同。而多重比较可以解决这个问题。例如:TukeyHSD()函数提供了对各组均值差异的成对检验。
# 对方差分析结果进行TukeyHSD多重比较
TukeyHSD(fit)
# 绘图
par(las=2)
par(mar=c(5,8,4,2))
plot(TukeyHSD(fit))
结果显示:1times和2times、2times和4times、drugD和4times之间的差异不显著。其余两两比较均达到显著水平。
multcomp包中的glht()函数提供了多重均值比较更为全面的方法,既适用于线性模型,也适用于广义线性模型。下面的代码可以重现TukeyHSD检验:
# 加载multcomp包
library(multcomp)
# 设置图形参数
par(mar=c(5,4,6,2))
# 进行多重比较
tuk<-glht(fit,linfct = mcp(trt="Tukey"))
# 图形展示
plot(cld(tuk,level=0.05),col="lightgrey")
2、评估检验的假设条件
单因素方差分析中,我们假设因变量服从正态分布,各组方差相等。我们可以使用Q-Q图来检验正态性假设,使用car包中的qqPlot()函数(注意:qqPlot()函数要求lm()拟合):
# 加载car包
library(car)
# 绘制Q-Q图
qqPlot(lm(response~trt,data=cholesterol),
simulate=TRUE,
main="Q-Q Plot",
labels=FALSE)
R语言提供了一些可用来做方差齐性检验的函数。例如Bartlett检验:
bartlett.test(response~trt,data=cholesterol)
Bartlett检验表明这5组数据的方差没有显著差异(p=0.9653)。当然还有其他检验如Fligner-Killeen检验(fligner.test()函数)和Brown-Forsythe检验(HH包中的hov()函数)。
不过方差齐性分析对离群点非常敏感。可以利用car包中outlierTest()函数来检验高离群点:
# 加载car包
library(car)
outlierTest(fit)
没有数据点的Bonferonni p值小于0.05,也就是说数据中没有离群点。因此根据Q-Q图、Bartlett检验和离群点检验,可以认为该数据可以用ANOVA模型拟合得很好。