方差分析(ANOVA)

2020-08-06 09:52:30 浏览数 (1)

方差分析主要用于多个样本均数比较的假设检验,因为当我们使用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()。接下来我将和大家讲解非参数假设检验,咱们下期再见!

0 人点赞