R 语言中的方差分析：单因素和双因素检验示例

什么是方差分析？

方差分析 方差分析 (ANOVA) 是一种用于比较两个或多个组均值的统计技术。该检验的工作原理是将测量值的总变异分解为由组别成员关系解释的部分和剩余的随机噪声部分。因此，R 语言中的 ANOVA 可以告诉你至少有一个组的均值与其他组存在差异，而不是具体是哪个组。它是……的直接扩展。 t检验 适用于因子变量具有两个以上水平的情况。

在进行测试之前，了解 ANOVA 家族中的哪个成员适合你的设计会很有帮助。

R语言中的ANOVA检验类型

“方差分析”是一系列检验方法的统称，而非单一的分析程序。选择合适的检验方法取决于有多少个因子以及数据的收集方式。

本教程涵盖前三种情况。其余变体使用相同的 aov() 接口，因此一旦您能够读取一个输出表，就可以读取所有输出表。

R语言中的方差分析与t检验：主要区别

这两个检验都是比较均值，因此有必要明确一个检验在什么情况下可以替代另一个检验。

面对三组数据，人们很容易会想到进行三次独立的t检验。但请克制这种做法。每次检验都有5%的误差率，因此三次比较会使假阳性率上升到约14%。方差分析（ANOVA）只需一次检验即可解答同样的问题，而Tukey HSD事后检验则可以在控制误差率的情况下处理两两比较的细节。对于两组数据的情况，请参阅相关内容。 t检验教程.

单向方差分析

很多时候都需要比较多个组之间的平均值，比如市场部想知道三个团队的销售业绩是否相同。

• 团队：3 个级别因素：A、B、C
• 销售：绩效衡量标准

方差分析检验可以判断三组的表现是否相似。

为了明确数据是否来自同一人群，你可以执行 单因素方差分析 （以下简称单因素方差分析）。与其他统计检验一样，它提供的证据可以判断原假设（H0）是否可以被拒绝。需要注意的是，未能拒绝原假设并不等同于证明原假设为真。

单向方差分析检验中的假设

• H0：组间均值相同
• H1：至少有一组的均值不同

换句话说，未能拒绝原假设 H0 意味着没有足够的证据得出任何一组均值与其他组均值不同的结论。

此检验与 t 检验类似，但当组数超过两组时，方差分析 (ANOVA) 才是正确的选择。当组数恰好为两组时，这两种检验是等价的，F 统计量等于 t 统计量的平方。

假设

单因素方差分析基于三个条件：观测值是随机抽样且彼此独立；每组内的残差近似服从正态分布；以及各组的方差相同（方差齐性）。下文“检验假设”部分将介绍如何在 R 语言中检验每个假设。

解释方差分析检验

F 统计量用于检验数据是否来自显著不同的总体，即不同的样本均值。

要计算 F 统计量，您需要除以 组间变异性 在 组内变异性.

此 组间 变异性反映了每个组均值与总体均值之间的距离。比较下面的两个图表即可理解。

左图显示三组之间的差异很小，因此三组的平均值都非常接近。 最划算 意思。

右图绘制了三个相距甚远且没有重叠的分布，因此总体均值与各组均值之间的差距很大。

此 组内 变异性衡量的是单个观测值与其所属组均值之间的偏差程度。有些点与组平均值相差甚远，组内项恰好反映了这种偏差，即抽样误差。

为了直观地了解组内变异性的概念，请看下图。

左侧部分绘制了三个不同组的分布情况。你增加了每个样本的离散程度，很明显个体方差很大。F 统计量下降，因此你无法拒绝原假设。

右侧部分显示的是均值相同但离散程度小得多的样本。这提高了F统计量，并支持备择假设。

您可以使用这两种度量来构建 F 统计量。理解 F 统计量非常直观。如果分子增加，则意味着组间变异性很高，样本中的组很可能来自完全不同的分布。

换句话说，较低的 F 统计量表明各组平均值之间几乎没有或根本没有有意义的差异。

单因素方差分析检验示例

您将使用毒药数据集来实现单向方差分析测试。该数据集包含 48 行和 3 个变量：

• 时间：动物的生存时间
• 毒药：所用毒药类型：因子水平：1,2、3 和 XNUMX
• 治疗：所用治疗类型：因子水平：1,2、3 和 XNUMX

在开始计算方差分析检验之前，您需要按如下方式准备数据：

• 步骤 1：导入数据
• 第 2 步：删除不必要的变量
• 步骤 3：将变量毒药转换为有序级别

library(dplyr)
PATH <- "https://raw.githubusercontent.com/guru99-edu/R-Programming/master/poisons.csv"
df <- read.csv(PATH) %>%
select(-X) %>% 
mutate(poison = factor(poison, ordered = TRUE))
glimpse(df)

输出：

## Observations: 48
## Variables: 3
## $ time   <dbl> 0.31, 0.45, 0.46, 0.43, 0.36, 0.29, 0.40, 0.23, 0.22, 0...
## $ poison <ord> 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 1, 1, 1, 1, 2, 2, 2...
## $ treat  <fctr> A, A, A, A, A, A, A, A, A, A, A, A, B, B, B, B, B, B, ...		

我们的目标是检验以下假设：

• H0：各组间平均生存时间无差异
• H1：至少有一个组的平均生存时间不同。

换句话说，你想知道根据给豚鼠喂食的毒药类型，其平均生存时间是否存在统计学差异。

您将按如下方式进行：

• 步骤 1：检查变量toxic的格式
• 第 2 步：打印汇总统计数据：计数、平均值和标准差
• 步骤 3：绘制箱线图
• 步骤 4：计算单向方差分析检验
• 步骤 5：使用 Tukey HSD 进行成对比较

步骤1） 使用以下代码检查中毒程度。您应该看到三个字符值，因为 mutate 动词将该列转换为有序因子。

levels(df$poison)

输出：

## [1] "1" "2" "3"

步骤2） 计算平均值和标准差。

df %>%
	group_by(poison) %>%
	summarise(
		count_poison = n(),
		mean_time = mean(time, na.rm = TRUE),
		sd_time = sd(time, na.rm = TRUE)
	)

输出：

## 
# A tibble: 3 x 4
##   poison count_poison mean_time    sd_time
##    <ord>        <int>     <dbl>      <dbl>
## 1      1           16  0.617500 0.20942779
## 2      2           16  0.544375 0.28936641
## 3      3           16  0.276250 0.06227627

步骤3） 在第三步中，您可以用图形检查分布之间是否存在差异。请注意，您要包括抖动点。

ggplot(df, aes(x = poison, y = time, fill = poison)) +
    geom_boxplot() +
    geom_jitter(shape = 15,
        color = "steelblue",
        position = position_jitter(0.21)) +
    theme_classic()

输出：

步骤4） 您可以使用命令 aov 运行单向方差分析测试。方差分析测试的基本语法是：

aov(formula, data)
Arguments:			
- formula: The equation you want to estimate
- data: The dataset used	

该公式的语法是：

y ~ X1+ X2+...+Xn # X1 +  X2 +... refers to the independent variables
y ~ . # use all the remaining variables as independent variables

现在你可以回答这个问题：考虑到所给豚鼠的毒药类型，它们的存活时间是否有差异？

将模型存储在一个对象中，并将其传递给 summary() 函数，以获得结果的易读打印输出。

anova_one_way <- aov(time~poison, data = df)
summary(anova_one_way)

Code 说明

• aov(time ~toxic, data = df)：使用以下公式运行方差分析测试
• summary(anova_one_way): 打印测试摘要

输出：

##             Df Sum Sq Mean Sq F value   Pr(>F)
## poison       2  1.033  0.5165   11.79 7.66e-05 ***
## Residuals   45  1.972  0.0438                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p 值为 7.66e-05，远低于通常的阈值 0.05，三个星号表示最显著性。您可以拒绝原假设 H0，并得出结论：至少有一个毒物组的平均生存时间与其他毒物组不同。

如何在 R 语言中检验方差分析的假设

只有当前面列出的三个条件都满足时，ANOVA 的 p 值才可信。每个条件都可以在 R 语言中直接验证，而且所有验证都是在拟合模型对象上运行的。

1. 观测结果的独立性。 这是研究设计本身的特性，而非数据的特性，因此任何检验都无法弥补。每只豚鼠都必须测量一次，并随机分配到相应的组别。如果同一只豚鼠出现在多行数据中，则需要使用重复测量模型。

2. 残差的正态性。 方差分析假设残差（而非原始数据）近似服从正态分布。请检查 QQ 图，并使用 Shapiro-Wilk 检验进行确认：

par(mfrow = c(2, 2))
plot(anova_one_way)          # four diagnostic plots

shapiro.test(residuals(anova_one_way))

正态QQ图中靠近对角线的点表明残差服从正态分布。Shapiro-Wilk检验的p值大于0.05意味着不能拒绝正态性假设。

3. 方差齐性。 每个组的分布应该相似。残差与拟合值的关系图应该呈平坦的带状，而不是漏斗状。使用 Levene 检验来验证这一点，该检验比 Bartlett 检验对非正态性更稳健：

library(car)
leveneTest(time ~ poison, data = df)

bartlett.test(time ~ poison, data = df)

p 值大于 0.05 支持方差相等。

当假设不成立时该怎么办？ 如果方差不相等，则运行 `oneway.test(time ~ poison, data = df, var.equal = FALSE)`，进行 Welch 校正。如果残差明显不符合正态分布且样本量较小，则切换到 Kruskal-Wallis 秩和检验，即 `kruskal.test(time ~ poison, data = df)`。对于大样本且样本平衡的情况，ANOVA 对偏离正态分布的程度相当稳健，因此 Shapiro-Wilk 检验结果接近临界值通常不会造成致命影响。

成对比较

显著性F检验表明各组均值并非完全相等，但无法指出具体是哪两组之间存在差异。Tukey 诚实显著性差异检验通过比较每一对均值，同时控制族错误率，解决了这个问题。

TukeyHSD(anova_one_way)

输出：

逐行读取输出结果。 差异 这一列保存的是两组均值之间的差值， 轻水反应堆 和 上升 确定该差异的 95% 置信区间，并且 p 调整 p 值是针对多重比较进行调整后的 p 值。当 p 值与 p 值的比较区间不包含零时，即 p 值小于 0.05 时，两组数据差异显著。在本数据集中，涉及毒物 3 的比较结果具有统计学意义，这与箱线图的结果一致：第 3 组的平均生存时间明显低于第 1 组和第 2 组，而第 1 组和第 2 组之间无统计学差异。

两向方差分析

双因素方差分析在公式中增加了一个因素。它的工作原理与单因素方差分析完全相同，只是公式有所不同：

y ~ x1 + x2

这里 y 是定量响应变量，而 x1 和 x2 都是分类因子。

双向方差分析检验中的假设

• H0：两个因子变量的组均值相等
• H1：至少有一个组的均值存在差异，且至少在两个因素中的一个方面存在差异。

您将处理变量添加到模型中。该变量记录豚鼠接受的处理方式。下面的加性公式检验在考虑其他因素后，每个因素是否单独影响生存时间。

修改代码，在第一个自变量旁边添加 treat。

anova_two_way <- aov(time~poison + treat, data = df)
summary(anova_two_way)

输出：

##             Df Sum Sq Mean Sq F value  Pr(>F)    
## poison       2 1.0330  0.5165   20.64 5.7e-07 ***
## treat        3 0.9212  0.3071   12.27 6.7e-06 ***
## Residuals   42 1.0509  0.0250                    
## ---

两个 p 值（毒药为 5.7e-07，治疗为 6.7e-06）都远小于 0.05，因此拒绝了两个因素的 H0，并得出结论：改变毒药或治疗都会影响生存时间。

添加交互项

上述加性模型假设毒物的影响与处理方式无关。为了检验这一假设，请将加号替换为星号，这样既可以拟合主效应，也可以拟合交互效应：

anova_interaction <- aov(time~poison * treat, data = df)
summary(anova_interaction)

如果中毒：治疗行不显著，则加性模型是更好的选择，因为它消耗的自由度更少。

R 语言中的方差分析：快速测试参考

下表列出了上面使用的每个测试、运行该测试的 R 调用以及它所评估的假设：

aov(y ~ X, data = df)

TukeyHSD(ANOVA summary)

aov(y ~ X1 + X2, data = df)

常见问题