用R语言做单因素方差分析及多重比较

SPSS方差分析的应用已经做得非常好了，绝大多数的方差分析问题均可通过SPSS“点菜单”的方式得以解决，R语言在统计和可视化方面有自己的特色，

我们不妨来对比着学习。

选用R语言自带案例数据集PlantGrowth，研究两个处理和一个对照组对植物产量的影响，每组10例共3记录，主要考察处理对提高植物产量有无影响。

数据构成：因变量weight，因子变量group，三个水平依次为ctrl、trt1、trt2。

数据正态性检验

用夏皮罗-威尔克检验3组数据是否服从正态分布，样本容量3~5000均可，这一点比SPSS说得更明确。

语法代码：

data <- PlantGrowth[,1]

shapiro.test(data[1:10])

shapiro.test(data[11:20])

shapiro.test(data[21:30])

ctrl、trt1、trt2三个分组的shapiro.test检验概率p依次为0.7475、0.4519、0.5643，均大于0.05，原假设成立（H0：假设数据服从正态分布），表明3组数据均来自正态分布总体。

方差齐次检验

用bartlett.test检验3个分组数据方差是否一致。

语法代码：

bartlett.test(weight~group,data = PlantGrowth)

直接看p-value = 0.2371>0.05，原假设成立（H0：假设3组数据方差相等），表明3组数据的方差齐次。

单因素方差分析，必须对数据正态性和方差齐次做出判断，如果不满足正态或方差齐次，则需要做出有关相应。（在《SPSS从入门到实践提高》课程中有具体讲述）

图形可视化

用箱图观察一下3组数据的分布情况。

单因素方差分析

用R语言aov函数完成方差分析。

语法代码：

dfc <- aov(weight~group,data = PlantGrowth)

summary(dfc)

R语言是用文本来输出结果，这一点比SPSS逊色多了，读取结果不够直观方便。

方差分析表的p-value=0.0159<0.05，拒绝原假设（H0：3组样本数据均值相等，无差异），说明3组样本均值有差异，不同处理方式对植物产量有显著影响。

多重比较

既然不同处理方式对产量有显著影响，那么我们就有必要深究一下，与控制组ctrl相比，处理1组trt1和处理2组trt2哪一个方案更佳。

选择Dunnett比较法。

语法代码：

library(multcomp)

dcbj <- glht(dfc,linfct = mcp(group = 'Dunnett'), alternative = 'two.side')

summary(dcbj)

与控制组相比，trt1的产量下降了0.371，p=0.323>0.05，原假设成立（H0：trt1和ctrl产量无差异），从绝对值来看trt产量下滑，但统计学上认为二者无差异。

与控制组相比，trt2的产量提高了0.494，p=0.153>0.05，原假设成立（H0：trt2和ctrl产量无差异），从绝对值来看trt产量虽然有所提升，但统计学上认为二者仍无差异。

大家辛苦研究，统计上认为两个处理方式竟然和对照组相比没有啥差异，有点小尴尬啊。不过这种被统计拒绝，未达到显著的情况也比较常见。

多重比较的方法非常多，Dunnett法没有LSD法灵敏，我们换LSD看看啥情况。

语法代码：

install.packages("agricolae")

dfclsd <- LSD.test(dfc,'group',p.adj = 'none')

dfclsd$groups

输出一个可视化效果

plot(dfclsd)

可见，两个处理与控制组相比均无差异。和Dunnett的结论一致。

如果用SPSS来做LSD的话，trt2与ctrl相比，差异检验的p值=0.088，比较接近拒绝域临界值水平0.05，如果我们将临界值调整为0.1的话，可说明处理2与对照组产量差异显著，那么在这项研究中，统计检验的拒绝域临界值水平调整为0.1是否恰当或可行，可由这个领域的研究者们讨论确定。

本文认为，可以确定拒绝域显著性水平a为0.1（有90%的把握），此时可认为处理2组trt2控制组相比差异明显，有统计学意义。

本文完

文/图=数据小兵

更多R统计文章

car包spm函数矩阵散点图
用R语言pie函数做饼图

ggplot2统计图形：常见的4种箱线图

ggplot2统计图形：常见的4种直方图

按列索引按列名称删除指定的列数据
R语言Levene方差齐次检验
谁说菜鸟不会数据分析（R语言篇）
table函数：分类数据的频数与频率统计

jiebaR包中文分词及词云制作

R语言带文字标签的散点图
如何获取R自带数据集与R包数据集说明文档？

如何修改R数据框的列名称？
二元正态分布及双变量相关分析简单案例演示

ggplot2统计图形：常见的4种散点图
R语言scale()函数实现数据标准化

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。