在接下来的例子里,我将会以小写字母表示数值型向量,而大写字母表示因子数据。
# 完全随机设计的单因素方差分析
# 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) #查看拟合的结果
# 析因设计
# 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)
诊断图主要是用来评估异方差性、正态性和对结果有影响的异常观测值。
layout(matrix(c(1,2,3,4),2,2)) # 创建2行2列的画布
plot(fit) # 绘制析因设计结果的诊断图
诊断图的横轴是拟合值,纵轴是残差、标准差或标准差的平方根,一般当各点的标准差集种在0处且分布较为均匀时,则说明拟合结果较好。上图显示2,3,5这三个样本的拟合值可能存在较大误差和,需仔细考虑。
在R中,我们可以使用函数anova(fit1, fit2)去评估不同模型的效应
fit1 <- aov(yield ~ N + block, data=npk)
fit2 <- aov(yield ~ N*P*K + block, data=npk)
anova(fit1,fit2)
在这里,你可以使用TukeyHSD()函数来进行Tukey HSD检验,它实际上是在方差分析结论有统计学意义之后进行的两两时候比较。
TukeyHSD(fit)
通常情况下,我们可以使用箱线图去可视化各组之间的差异,但我们也可以使用专门用于方差分析的基础可视化函数interaction.plot()和来自“gplots”包的plotmeans()函数。
# 绘制两因素互作图
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时,这种差异几乎没有了。
attach(mtcars)
cyl <- factor(cyl)
plotmeans(mpg~cyl,xlab="Number of Cylinders",
ylab="Miles Per Gallon", main="Mean Plot\nwith95% CI")
上图的误差条实际上反映各个均值的95%置信区间,这里就不赘述了。
假如你有多个因变量,这时你可以使用多元方差分析(MANOVA)的方法来处理,这里因变量通常是一个矩阵,而使用的函数是manova()。
关于方差分析的内容就先讲到这儿,注意方差分析的核心函数是aov()。接下来我将和大家讲解非参数假设检验,咱们下期再见!