方差分析主要用于多个样本均数比较的假设检验,因为当我们使用t检验进行多组样本间均数的假设检验时,常常会增加一类错误。方差分析的主要思想是分解变异,即将总变异分解为处理因素引起的变异和随机误差引起的变异,通过对两者进行比较做出处理因素有无作用的统计推断。在后续的内容中我将会和大家详细讲解方差分析的统计原理。
在R语言进行方差分析是一件非常方便的事,我们通常只需要进行5步即可完成较高质量的方差分析,这五步主要是拟合模型、诊断性绘图、评估模型效应、多重比较和结果可视化。
这里使用的是R里内置的“npk”数据集,该数据集由24行和5列数据组成,第一列代表区组(共6个),N、P和K分别代表氮、磷和钾元素的使用情况,yield代表豌豆产量,该数据集主要是用来研究不同肥料对豌豆产量的影响。
1. 拟合模型
在接下来的例子里,我将会以小写字母表示数值型向量,而大写字母表示因子数据。
代码语言:javascript复制# 完全随机设计的单因素方差分析
# fit <- aov(y ~ A, data=mydataframe) #y是数值向量,A是因子
fit <- aov(yield ~ N, data=npk) #只是示例,实际中不是这么处理的
# 随机区组设计(B代表区组)
# fit <- aov(y ~ A B, data=mydataframe) #y是数值向量,A、B是因子
fit <- aov(yield ~ N block, data=npk) #示例,研究氮元素对豌豆生长的影响
summary(fit) #查看拟合的结果
代码语言:javascript复制# 析因设计
# fit <- aov(y ~ A B A:B, data=mydataframe) # y是数值向量,A、B是因子
# fit <- aov(yield ~ A*B, data=mydataframe) # 和上面的代码的作用是
fit <- aov(yield ~ N*P*K block, data=npk) #三因素析因设计,block位区组,考虑N、P和K之间的交互作用
summary(fit)
2. 查看诊断分析图
诊断图主要是用来评估异方差性、正态性和对结果有影响的异常观测值。
代码语言:javascript复制layout(matrix(c(1,2,3,4),2,2)) # 创建2行2列的画布
plot(fit) # 绘制析因设计结果的诊断图
诊断图的横轴是拟合值,纵轴是残差、标准差或标准差的平方根,一般当各点的标准差集种在0处且分布较为均匀时,则说明拟合结果较好。上图显示2,3,5这三个样本的拟合值可能存在较大误差和,需仔细考虑。
3. 评估模型效应
在R中,我们可以使用函数anova(fit1, fit2)去评估不同模型的效应
代码语言:javascript复制fit1 <- aov(yield ~ N block, data=npk)
fit2 <- aov(yield ~ N*P*K block, data=npk)
anova(fit1,fit2)
4. 多重比较
在这里,你可以使用TukeyHSD()函数来进行Tukey HSD检验,它实际上是在方差分析结论有统计学意义之后进行的两两时候比较。
代码语言:javascript复制TukeyHSD(fit)
5. 结果的可视化
通常情况下,我们可以使用箱线图去可视化各组之间的差异,但我们也可以使用专门用于方差分析的基础可视化函数interaction.plot()和来自“gplots”包的plotmeans()函数。
代码语言:javascript复制# 绘制两因素互作图
attach(mtcars) #固定数据集
gear <- factor(gear) #转换为因子
cyl <- factor(cyl) #转换为因子
interaction.plot(cyl, gear, mpg, type="b", col=c(1:3),
leg.bty="o", leg.bg="beige", lwd=2,pch=c(18,24,22),
xlab="Number of Cylinders",
ylab="Mean Miles Per Gallon",
main="Interaction Plot") #cyl作为x轴,gear作为区组,mpg作为y轴
从上图中我们可以清晰看出mpg和cyl在不同gear组的变化关系,总的来看,随着cyl的增加,mpg在减少,当cyl在4~6范围时,不同gear的mpg差异较大,当cyl>6时,这种差异几乎没有了。
代码语言:javascript复制attach(mtcars)
cyl <- factor(cyl)
plotmeans(mpg~cyl,xlab="Number of Cylinders",
ylab="Miles Per Gallon", main="Mean Plotnwith95% CI")
上图的误差条实际上反映各个均值的95%置信区间,这里就不赘述了。
6. 多元方差分析
假如你有多个因变量,这时你可以使用多元方差分析(MANOVA)的方法来处理,这里因变量通常是一个矩阵,而使用的函数是manova()。
关于方差分析的内容就先讲到这儿,注意方差分析的核心函数是aov()。接下来我将和大家讲解非参数假设检验,咱们下期再见!