第十五讲 R-单因素方差分析2

Original 跟我学投必得医学 2022-05-07

收录于合集 #R与生物统计专题 56个

在“R与生物统计专题”中，我们会从介绍R的基本知识展开到生物统计原理及其在R中的实现。以从浅入深，层层递进的形式在投必得医学公众号更新。

上一讲，我们讲了方差分析的原理和R实现（第十四讲 R-单因素方差分析1），但是，在发现至少存在一组与其他组之间存在差异以后，我们如何知道具体是哪两组间存在差异。方差分析需要满足几个条件，他们的假设检验条件的检验方法是什么？如果数据不满足假设条件，要怎么办？

今天，小编就带大家来一一回答这些问题。

1. 两两成对比较

方差分析各组均值之间的两两成对比较

首先，如第十四讲，我们已经完成了如下操作，发现方差分析P值 < 0.05。

# 导入R内自带的PlantGrowth数据集library(datasets)data(PlantGrowth)

# 将数据存储在变量my_data中my_data <- PlantGrowth

# 计算方差分析res.aov <- aov(weight ~ group, data = my_data)

# 总结分析结果summary(res.aov) #P = 0.016

在单因素方差分析中，显著的p值表示某些组均值不同，但我们不知道具体哪些组之间不同。

这时，可以执行多个成对比较，以确定特定组对之间的平均差异是否具有统计显著性。

我们可以计算Tukey HSD（Tukey Honest Significant Differences，R函数：TukeyHSD()可以实现），以在组均值之间进行多次成对比较。

函数TukeyHD()将之前拟合生成的ANOVA函数作为参数：

TukeyHSD(res.aov)

输出结果

Tukey multiple comparisons of means95% family-wise confidence levelFit: aov(formula = weight ~ group, data = my_data)$groupdiff lwr upr p adjtrt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711trt2-ctrl 0.494 -0.1972161 1.1852161 0.1979960trt2-trt1 0.865 0.1737839 1.5562161 0.0120064

diff：
两组平均值之间的差异
lwr，upr：
置信区间的上下端点为95％（默认值）
p adj：
调整后的多个比较的p值。

从输出中可以看出，只有trt2与trt1之间的差异显着，调整后的p值为0.012。

2. 检查方差分析的假设检验的条件

如开头讲到，ANOVA需要满足一定的条件：

1. 观察值之间在所研究的因子水平，是相互独立且随机从总体中获得的。
2. 每个因子水平的数据（每一组）均呈正态分布。
3. 每个组的样本的方差之间没有差异。

我们需要对假设数据是正态分布的，并且各组之间的方差是均匀的进行验证。

我们可以通过如下方法来检验。

2.1 检查方差假设的同质性

（方差齐性）

2.1.1 残差与拟合曲线图来检查方差齐性

在下面的图中，残差和拟合值（每组的平均值）之间没有明显的关系，因此，我们可以假设方差是同质的。

# 方差齐性检验plot(res.aov, 1)

该图会自动标出明显的异常值。途中点17、15、4被检测为异常值，这会严重影响方差的正态性和同质性。如果此时出现方差不齐，会建议删除异常值以满足假设成立。

2.1.2 Levene 检验

也可以使用Bartlett检验或Levene检验来检验方差齐性。

我们建议使用Levene检验，该检验对偏离正态分布不敏感。可以使用函数leveneTest()[在car包装中]：

library(car)leveneTest(weight ~ group, data = my_data)

输出结果

Levene's Test for Homogeneity of Variance (center = median)Df F value Pr(>F)group 2 1.1192 0.341227

从上面的输出中，我们可以看到p值>0.05的显着性水平。这意味着没有证据表明各组之间的差异在统计学上有显着差异。因此，我们可以假设不同治疗组中方差同质。

2.2 放宽方差假设的同质性

经典的单因素方差分析要求所有组的方差均等。

在我们的示例中，方差假设的同质性满足：Levene检验不显着。

在方差齐性被违反的情况下，我们如何继续方差分析呢？

另一个统计学方法（即：Welch one-way test），它不需要满足方差齐性的假说。函数oneway.test()可以实现。

方差不齐情况下的的方差分析

oneway.test(weight ~ group, data = my_data)

方差不齐情况下的成对t检

pairwise.t.test(my_data$weight, my_data$group,p.adjust.method = "BH", pool.sd = FALSE)

2.3 检查正态性假设

2.3.1 残差的正态图

在下面的图中，相对于正态分布的分位数绘制了残差的QQ图，并且还绘制了45度参考线。

残差的正态概率图用于检查残差呈正态分布的假设。如果它满足正态性，应该要大致分布在45度参考线上。

plot(res.aov, 2)

当所有点都大致沿着该参考线落下时，我们可以假定为正态。

2.3.2 Shapiro-Wilk检验

此外，我们可以使用Shapiro-Wilk检验检查残差是否符合正态分布。

# 提取残差aov_residuals <- residuals(object = res.aov )# Shapiro-Wilk 检验shapiro.test(x = aov_residuals )

输出结果

Shapiro-Wilk normality testdata: aov_residualsW = 0.96607, p-value = 0.4379

以上结论显示P > 0.05, 残差符合正态分布（ p = 0.6），没有发现违反正态性的迹象。

3. 非参数替代单因素方差分析

注意，单因素ANOVA的非参数替代方法是Kruskal-Wallis 秩和检验，当数据不满足ANOVA假设条件时，我们可以使用该检验。

3.1 Kruskal-Wallis 秩和检验

kruskal.test(weight ~ group, data = my_data)

输出结果

Kruskal-Wallis rank sum testdata: weight by groupKruskal-Wallis chi-squared = 7.9882, df = 2, p-value = 0.01842

3.2 多组之间的两两成对比较

从Kruskal-Wallis检验的输出中，我们知道组之间存在显着差异，但我们不知道具体哪些组是不同的。

可这时，以使用pairwise.wilcox.test() 函数对组别之间进行成对比较，并进行多重检验的更正。

pairwise.wilcox.test(PlantGrowth$weight, PlantGrowth$group,p.adjust.method = "BH")

输出结果

Pairwise comparisons using Wilcoxon rank sum testdata: PlantGrowth$weight and PlantGrowth$groupctrl trt1trt1 0.199 -trt2 0.095 0.027P value adjustment method: BH

成对比较显示，只有trt1和trt2显着不同（p <0.05）。

好了，本期讲解就先到这里。小伙伴们赶紧试起来吧。

在之后的更新中，我们会进一步为您介绍R的入门，以及常用生物统计方法和R实现。欢迎关注，投必得医学手把手带您走入R和生物统计的世界。

提前打个预告，下一期我们将学习“双向方差分析及其在R语言中的实现”。

第一讲 R-基本介绍及安装

第二讲 R-编程基础-运算、数据类型和向量等基本介绍