R中的双向方差分析

学习如何在R中进行双因素方差分析。您还将学习其目的、假设、假定和如何解释结果。

介绍

双因素方差分析（analysis of variance）是一种统计方法，允许评估两个分类变量对一个定量连续变量的同时影响。

双因素方差分析是一元方差分析的扩展，因为它允许评估两个分类变量对数字响应的影响，而不是一个。相比于一元方差分析，双因素方差分析的优点在于我们测试了两个变量之间的关系，同时考虑了第三个变量的影响。此外，它还允许包括两个分类变量对响应的可能交互作用。

双因素方差分析相对一元方差分析的优点类似于相关性相对多元线性回归的优点：

相关性用于衡量两个定量变量之间的关系。多元线性回归也用于衡量两个变量之间的关系，但这次考虑了其他协变量的潜在影响。
一元方差分析测试定量变量在组之间是否不同。双因素方差分析也测试定量变量在组之间是否不同，但这次考虑了另一个定性变量的影响。

之前，我们已经讨论了在R中进行一元方差分析。现在，我们将展示何时、为什么以及如何在R中执行双因素方差分析。

在进一步之前，我想提一下并简要描述一些相关的统计方法和测试，以避免任何混淆：

学生t检验用于评估一个分类变量对定量连续变量的影响，当分类变量具有两个水平时：

独立样本的学生t检验，如果观测值是独立的（例如：如果我们比较女性和男性之间的年龄）
成对样本的学生t检验，如果观测值相关，即它们成对出现（当同一被试在两个不同的时间点（例如，在治疗前后）被测量两次时，就是这种情况）

为了评估一个分类变量对定量变量的影响，当分类变量有3个或更多水平时：

一元方差分析（通常简称为ANOVA），如果组独立（例如接受治疗A的患者组、接受治疗B的患者组以及接受无治疗或安慰剂的患者组）
重复测量ANOVA，如果组相关（当同一被试在三个不同的时间点（例如，在治疗前、期间和后）分别被测量三次时，就是这种情况）

双因素方差分析用于评估两个分类变量（及其潜在交互作用）对定量连续变量的影响。这是本文的主题。

线性回归用于评估定量连续依赖变量与一个或多个独立变量之间的关系：

如果只有一个独立变量（可以是定量或定性），则使用简单线性回归
如果有至少两个独立变量（可以是定量、定性或两者的混合），则使用多元线性回归

ANCOVA（协方差分析）用于评估一个分类变量对定量变量的影响，同时控制另一个定量变量的影响（称为协变量）。实际上，ANCOVA是一种特殊情况的多元线性回归，其中有一个定性和一个定量独立变量的混合。

在本文中，我们首先解释何时和为什么使用双因素方差分析，然后进行一些初步的描述性分析，并介绍如何在R中进行双因素方差分析。最后，我们展示如何解释和可视化结果。我们还简要提到并说明如何验证基本假设。

数据

为了说明如何在R中执行双向方差分析，我们使用来自{palmerpenguins}包中的penguins数据集。

我们不需要导入数据集，但需要先加载该包，然后调用数据集：

# install.packages("palmerpenguins")library(palmerpenguins)

dat <- penguins # 重命名数据集str(dat) # 数据集结构

## tibble [344 × 8] (S3: tbl_df/tbl/data.frame)##  $ species          : 因子 w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...##  $ island           : 因子 w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...##  $ bill_length_mm   : 数值 [1:344] 39.1 39.5 40.3 NA 36.7 39.3 38.9 39.2 34.1 42 ...##  $ bill_depth_mm    : 数值 [1:344] 18.7 17.4 18 NA 19.3 20.6 17.8 19.6 18.1 20.2 ...##  $ flipper_length_mm: 整数 [1:344] 181 186 195 NA 193 190 181 195 193 190 ...##  $ body_mass_g      : 整数 [1:344] 3750 3800 3250 NA 3450 3650 3625 4675 3475 4250 ...##  $ sex              : 因子 w/ 2 levels "female","male": 2 1 1 NA 1 2 1 2 NA NA ...##  $ year             : 整数 [1:344] 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...

数据集包含344只企鹅的8个变量，总结如下：

summary(dat)

##       species          island    bill_length_mm  bill_depth_mm  ##  Adelie   :152   Biscoe   :168   Min.   :32.10   Min.   :13.10  ##  Chinstrap: 68   Dream    :124   1st Qu.:39.23   1st Qu.:15.60  ##  Gentoo   :124   Torgersen: 52   Median :44.45   Median :17.30  ##                                  Mean   :43.92   Mean   :17.15  ##                                  3rd Qu.:48.50   3rd Qu.:18.70  ##                                  Max.   :59.60   Max.   :21.50  ##                                  NA's   :2       NA's   :2      ##  flipper_length_mm  body_mass_g       sex           year     ##  Min.   :172.0     Min.   :2700   female:165   Min.   :2007  ##  1st Qu.:190.0     1st Qu.:3550   male  :168   1st Qu.:2007  ##  Median :197.0     Median :4050   NA's  : 11   Median :2008  ##  Mean   :200.9     Mean   :4202                Mean   :2008  ##  3rd Qu.:213.0     3rd Qu.:4750                3rd Qu.:2009  ##  Max.   :231.0     Max.   :6300                Max.   :2009  ##  NA's   :2         NA's   :2

在本文中，我们将重点关注以下三个变量：

species：企鹅的物种（Adelie、Chinstrap或Gentoo）
sex：企鹅的性别（雌性和雄性）
body_mass_g：企鹅的体重（以克为单位）

如果需要，可以通过在R中运行?penguins来找到有关此数据集的更多信息。

body_mass_g是定量连续变量，将是因变量，而species和sex都是定性变量。

这两个变量将是我们的自变量，也称为因子。请确保它们被 R 读取为因子。如果不是这种情况，它们将需要转换为因子。

双因素方差分析的目的和假设

如上所述，双因素方差分析用于同时评估两个分类变量对一个数量连续变量的影响。

它被称为双因素方差分析，因为我们正在比较由两个独立分类变量形成的组。

在这里，我们想知道体重是否取决于物种和/或性别。特别地，我们感兴趣的是：

测量和检验物种与体重之间的关系，
测量和检验性别与体重之间的关系，以及
潜在地检查物种与体重之间的关系对雌性和雄性是否不同（这等同于检查性别与体重之间的关系是否取决于物种）

前两种关系被称为主效应，而第三点则称为交互作用效应。

主效应测试是否至少有一组与另一组不同（同时控制其他自变量）。另一方面，交互作用效应旨在测试两个变量之间的关系是否取决于第三个变量的水平。

在执行双因素方差分析时，测试交互作用效应并非强制性的。但是，如果存在交互作用效应，则忽略交互作用效应可能会导致错误的结论。

如果我们回到我们的例子，我们有以下假设检验：

性别对体重的主效应：

H0：女性和男性的平均体重相等
H1：女性和男性的平均体重不同

物种对体重的主效应：

H0：所有三种物种的平均体重相等
H1：至少有一种物种的平均体重不同

性别和物种之间的交互作用：

H0：性别和物种之间没有交互作用，这意味着物种与体重之间的关系对于雌性和雄性是相同的（同样，性别和体重之间的关系对于所有三种物种也是相同的）
H1：性别和物种之间存在交互作用，这意味着物种与体重之间的关系对于女性和男性是不同的（同样，性别和体重之间的关系取决于物种）

双因素方差分析的假设

大多数统计测试需要一些假设才能使结果有效，双因素方差分析也不例外。

双因素方差分析的假设与单因素方差分析类似。概括起来：

变量类型：因变量必须是数量连续的，而两个自变量必须是分类的（至少有两个水平）。
独立性：组间和每个组内的观察结果应该是独立的。
正态性：
对于小样本，数据应近似遵循正态分布
对于大样本（通常每组/样本 n ≥ 30），不需要正态性（由于中心极限定理）
方差齐性：各组方差应相等。
异常值：任何组中都不应有显著的异常值。

关于这些假设的更多细节可以在单因素方差分析的假设中找到。

现在我们已经看到了双因素方差分析的基本假设，在应用测试和解释结果之前，我们会针对我们的数据集具体审查这些假设。

变量类型

因变量体重是数量连续的，而两个自变量性别和物种都是定性变量（至少有2个水平）。

因此，这个假设得到满足。

独立性

独立性通常基于实验设计和数据收集方式进行检查。

简单来说，观测通常是:

独立的，如果每个实验单位（这里是企鹅）仅被测量一次，并且观测是从人口的代表性和随机选择的部分收集的，或
依赖的，如果每个实验单位至少被测量两次（在医学领域中通常是这种情况，例如对同一受试者进行两次测量; 一次在治疗前，一次在治疗后）。

在我们的情况下，每只企鹅的体重只被测量了一次，并且是从人口的代表性和随机样本收集的，因此满足独立性假设。

正态性

我们在所有子组（即两个因素水平的每个组合，称为单元格）中都有大样本：

table(dat$species, dat$sex)

##            ##             female male##   Adelie        73   73##   Chinstrap     34   34##   Gentoo        58   61

因此不需要检查正态性。

为了完整起见，我们仍然展示如何验证正态性，就好像我们有一个小样本一样。

有几种方法可以测试正态性假设。最常见的方法是：

按组或在残差上进行QQ图，和/或
按组或在残差上进行直方图，和/或
按组或在残差上进行正态性测试（例如Shapiro-Wilk测试）。

最简单/最短的方法是使用残差的QQ图验证正态性。要绘制此图，我们首先需要保存模型：

# save modelmod <- aov(body_mass_g ~ sex * species,  data = dat)

这段代码将进一步解释。

现在我们可以在残差上绘制QQ图。我们展示两种绘制方法，第一种是使用plot()函数，第二种是使用{car}包中的qqPlot()函数：

# method 1plot(mod, which = 2)

# method 2library(car)

qqPlot(mod$residuals,  id = FALSE # remove point identification)

方法1的代码略短，但它错过了参考线周围的置信区间。

如果点遵循直线（称为亨利线）并落在置信带内，我们可以假定正态性。在这里是这种情况。

如果您更喜欢基于残差的直方图来验证正态性，则可以使用以下代码：

# histogramhist(mod$residuals)

残差的直方图显示高斯分布，与QQ图的结论相符。

尽管QQ图和直方图在很大程度上足以验证正态性，但是如果您想更正式地进行统计测试来验证它，则可以在残差上应用Shapiro-Wilk测试：

# normality testshapiro.test(mod$residuals)

## ##  Shapiro-Wilk normality test## ## data:  mod$residuals## W = 0.99776, p-value = 0.9367

⇒ 我们不拒绝残差遵循正态分布的零假设（p-value = 0.937）。

从 QQ 图，直方图和 Shapiro-Wilk 检验中我们得出结论，不拒绝残差正态性的零假设。

因此，正态性假设得到验证，现在我们可以检查方差的相等性。

方差的同质性

方差的相等性，也称为同质性或同方差性，可以通过 plot() 函数进行可视化检验：

plot(mod, which = 3)

由于残差的散布是恒定的，因此红色光滑线是水平的和平的，因此在这里似乎满足恒定方差的假设。

上述诊断图足够，但如果您愿意，也可以通过 Levene 检验（也来自 {car} 包）进行更正式的测试：

leveneTest(mod)

## 方差同质性的 Levene 检验（中心 = 中位数）##        Df F value Pr(>F)## group   5  1.3908 0.2272##       327

⇒ 我们不拒绝方差相等的零假设（p 值 = 0.227）。

可视化和正式方法均得出相同的结论；我们不拒绝方差同质性的假设。

异常值

检测异常值的最简单和最常见的方法是通过组别的箱线图进行可视化。

对于雌性和雄性：

library(ggplot2)

# 按性别绘制箱线图ggplot(dat) +  aes(x = sex, y = body_mass_g) +  geom_boxplot()

对于三个物种：

# 按物种绘制箱线图ggplot(dat) +  aes(x = species, y = body_mass_g) +  geom_boxplot()

根据四分位距准则，有两个 Chinstrap 物种的异常值。然而，这些点不足以偏离结果。

因此，我们认为没有显著的异常值假设得到满足。

双因素方差分析

我们已经证明了所有假设的满足，因此现在我们可以在 R 中实现双因素方差分析。

这将使我们能够回答以下研究问题：

在控制物种的情况下，两性之间的体重是否存在显着差异？
在控制性别的情况下，至少有一种物种的体重是否存在显着差异？
雌性和雄性企鹅之间的物种和体重之间的关系是否不同？

初步分析

在执行任何统计测试之前，最好进行一些描述性统计，以便对数据有一个初步的概述，并可能了解预期的结果。

这可以通过描述性统计或图表来完成。

描述性统计

如果我们想保持简单，我们可以仅计算每个子组的平均值：

# 按组计算平均值aggregate(body_mass_g ~ species + sex,  data = dat,  FUN = mean)

##     species    sex body_mass_g## 1    Adelie female    3368.836## 2 Chinstrap female    3527.206## 3    Gentoo female    4679.741## 4    Adelie   male    4043.493## 5 Chinstrap   male    3938.971## 6    Gentoo   male    5484.836

或者，使用{dplyr}包计算每个子组的平均值和标准差：

# 按组计算平均值和标准差library(dplyr)

group_by(dat, sex, species) %>%  summarise(    mean = round(mean(body_mass_g, na.rm = TRUE)),    sd = round(sd(body_mass_g, na.rm = TRUE))  )

## # A tibble: 8 × 4## # Groups:   sex [3]##   sex    species    mean    sd##   <fct>  <fct>     <dbl> <dbl>## 1 female Adelie     3369   269## 2 female Chinstrap  3527   285## 3 female Gentoo     4680   282## 4 male   Adelie     4043   347## 5 male   Chinstrap  3939   362## 6 male   Gentoo     5485   313## 7 <NA>   Adelie     3540   477## 8 <NA>   Gentoo     4588   338

绘图

如果您经常阅读博客，您就会知道我喜欢在解释测试结果之前绘制图表来可视化数据。

当我们有一个定量变量和两个定性变量时，最适合的图表是按组绘制的箱线图。这可以使用{ggplot2}包轻松实现：

# 按组绘制箱线图library(ggplot2)

ggplot(dat) +  aes(x = species, y = body_mass_g, fill = sex) +  geom_boxplot()

对于性别存在缺失值的观测值，我们可以删除它们以获得更简洁的图表：

dat %>%  filter(!is.na(sex)) %>%  ggplot() +  aes(x = species, y = body_mass_g, fill = sex) +  geom_boxplot()

请注意，我们还可以制作以下图表：

dat %>%  filter(!is.na(sex)) %>%  ggplot() +  aes(x = sex, y = body_mass_g, fill = species) +  geom_boxplot()

但是，为了使图表更易读，我倾向于将具有最少类别数的变量作为颜色（实际上是在aes()层中的fill参数），将具有最多类别的变量放在x轴上（即，aes()层中的x参数）。

根据子组的平均值和箱线图，我们可以看出，在我们的样本中：

雌性企鹅的体重通常比雄性企鹅低，对于所有考虑的物种而言都是如此，以及
Gentoo企鹅的体重高于其他两个物种。

请记住，这些结论仅在我们的样本内有效！要将这些结论推广到人群中，我们需要执行双因素方差分析，并检查解释变量的显着性。这是下一节的目的。

R中的双因素方差分析

如前所述，在双因素方差分析中包含交互作用效应并非强制性。但是，为了避免错误的结论，建议首先检查交互作用是否显着，然后根据结果是否包含它。

如果交互作用不显著，可以将其从最终模型中删除。相反，如果交互作用显著，应该将其包括在最终模型中，以用于解释结果。

因此，我们从一个包括两个主要效应（即性别和物种）和交互作用的模型开始：

# 带有交互作用的双向方差分析# 保存模型mod <- aov(body_mass_g ~ sex * species,  data = dat)

# 打印结果summary(mod)

##              Df    Sum Sq  Mean Sq F value   Pr(>F)    ## sex           1  38878897 38878897 406.145  < 2e-16 ***## species       2 143401584 71700792 749.016  < 2e-16 ***## sex:species   2   1676557   838278   8.757 0.000197 ***## Residuals   327  31302628    95727                     ## ---## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## 由于缺失值删除了11个观测值

平方和（列Sum Sq）显示物种解释了体重变异的大部分。它是解释这种变异最重要的因素。

上面输出的最后一列显示了p值（Pr(>F)）。从这些p值可以得出以下结论（在5%的显著性水平下）：

在控制物种的情况下，两性之间的体重显著不同
在控制性别的情况下，至少有一种物种的体重显著不同
性别和物种之间的交互作用（在上面输出的sex:species行中显示）是显著的。

因此，从显著的交互作用效应可以看出，体重和物种之间的关系在雄性和雌性之间有所不同。由于这是显著的，我们必须将其保留在模型中，并应该从该模型中解释结果。

相反，如果交互作用不显著（即p值≥0.05），我们将从模型中删除此交互作用效应。为了说明目的，下面是没有交互作用的双向方差分析的代码，称为加性模型：

# 无交互作用的双向方差分析aov(body_mass_g ~ sex + species,  data = dat)

对于习惯于在R中执行线性回归的读者，您将注意到双向方差分析的代码结构实际上是类似的：

公式为因变量 ~ 自变量
使用+符号包括没有交互作用的自变量
使用*符号包括具有交互作用的自变量

与线性回归的相似之处并不令人意外，因为双向方差分析（像所有方差分析一样）实际上是一个线性模型。

请注意，以下代码也有效，并且给出相同的结果：

# 方法2mod2 <- lm(body_mass_g ~ sex * species,  data = dat)

Anova(mod2)

## Anova Table (Type II tests)## ## Response: body_mass_g##                Sum Sq  Df F value    Pr(>F)    ## sex          37090262   1 387.460 < 2.2e-16 ***## species     143401584   2 749.016 < 2.2e-16 ***## sex:species   1676557   2   8.757 0.0001973 ***## Residuals    31302628 327                      ## ---## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

请注意，aov()函数假定平衡设计，这意味着我们在独立分组变量的水平内具有相等的样本大小。

对于非平衡设计，即不同子组中的受试者数量不相等，建议使用以下方法：

当不存在显著的交互作用时，可以使用 Anova(mod, type = "II") 在 R 中进行二型方差分析，得出 5
当存在显著的交互作用时，可以使用 Anova(mod, type = "III") 在 R 中进行三型方差分析。

这已经超出了本篇文章的范围，我们假设这里是一个平衡的设计。对于感兴趣的读者，请参阅有关一型、二型和三型方差分析的详细讨论。

成对比较

通过两个主效应显著，我们得出结论：

在控制物种的情况下，雌性和雄性的体重不同，
在控制性别的情况下，至少有一种物种的体重不同。

如果两性的体重不同，而且正好有两性，那么这一定是因为雌性和雄性的体重显著不同。

如果想知道哪个性别的体重最高，可以从子组的均值和/或箱线图中推断出来。在这里，很明显雄性的体重显著高于雌性。

然而，对于物种来说并不是那么简单。让我解释一下为什么对于性别来说并不像对于物种那么容易。

有三种物种 (Adelie, Chinstrap 和 Gentoo)，所以有 3 对物种：

Adelie 和 Chinstrap
Adelie 和 Gentoo
Chinstrap 和 Gentoo

如果至少有一种物种的体重显著不同，那么可能是：

Adelie 和 Chinstrap 的体重显著不同，但 Adelie 和 Gentoo 以及 Chinstrap 和 Gentoo 的体重没有显著不同，或
Adelie 和 Gentoo 的体重显著不同，但 Adelie 和 Chinstrap 以及 Chinstrap 和 Gentoo 的体重没有显著不同，或
Chinstrap 和 Gentoo 的体重显著不同，但 Adelie 和 Chinstrap 以及 Adelie 和 Gentoo 的体重没有显著不同。

或者，也可能是：

Adelie 和 Chinstrap 的体重显著不同，Adelie 和 Gentoo 的体重显著不同，但 Chinstrap 和 Gentoo 的体重没有显著不同，或
Adelie 和 Chinstrap 的体重显著不同，Chinstrap 和 Gentoo 的体重显著不同，但 Adelie 和 Gentoo 的体重没有显著不同，或
Chinstrap 和 Gentoo 的体重显著不同，Adelie 和 Gentoo 的体重显著不同，但 Adelie 和 Chinstrap 的体重没有显著不同。

最后，所有物种的体重也可能显著不同。

与单向方差分析一样，我们在这个阶段无法确切地知道哪个物种在体重方面与哪个物种不同。为了知道这一点，我们需要通过事后检验 (也称为成对比较) 来逐一比较每两种物种。

有几种事后检验方法，最常见的是 Tukey HSD，它测试所有可能的组合。正如前面提到的，这个测试只需要在物种变量上进行，因为性别只有两个水平。

与单向方差分析一样，Tukey HSD 可以在 R 中进行如下操作：

# 方法 1TukeyHSD(mod,  which = "species")

##   Tukey multiple comparisons of means##     95% family-wise confidence level## ## Fit: aov(formula = body_mass_g ~ sex * species, data = dat)## ## $species##                        diff       lwr       upr     p adj## Chinstrap-Adelie   26.92385  -80.0258  133.8735 0.8241288## Gentoo-Adelie    1377.65816 1287.6926 1467.6237 0.0000000## Gentoo-Chinstrap 1350.73431 1239.9964 1461.4722 0.0000000

或使用{multcomp}包：

# 方法 2library(multcomp)

summary(glht(  aov(body_mass_g ~ sex + species,    data = dat  ),  linfct = mcp(species = "Tukey")))

## ##   广义线性假设的同时检验## ## 多组均值比较：Tukey对比## ## ## 模型：aov(formula = body_mass_g ~ sex + species, data = dat)## ## 线性假设：##                         Estimate Std. Error t value Pr(>|t|)    ## Chinstrap - Adelie == 0    26.92      46.48   0.579     0.83    ## Gentoo - Adelie == 0     1377.86      39.10  35.236   <1e-05 ***## Gentoo - Chinstrap == 0  1350.93      48.13  28.067   <1e-05 ***## ---## 显著性水平:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## (校正后的p值报告 -- 单步法)

或使用您选择的p值调整方法的pairwise.t.test()函数：

# 方法 3pairwise.t.test(dat$body_mass_g, dat$species,  p.adjust.method = "BH")

## ## 使用合并SD的t检验进行成对比较## ## 数据：dat$body_mass_g 和 dat$species## ##           Adelie Chinstrap## Chinstrap 0.63   -        ## Gentoo    <2e-16 <2e-16   ## ## p值校正方法：BH

请注意，当使用第二种方法时，即使交互作用显著，也需要将没有交互作用的模型指定到glht()函数中。此外，不要忘记在我的代码中用您的模型的名称和自变量的名称替换mod和species。

这两种方法都给出了相同的结果，即：

帝企鹅和鹈鹕的体重没有显著差异（校正后的p值为0.83），
小企鹅和帝企鹅的体重显著不同（校正后的p值小于0.001），
小企鹅和鹈鹕的体重显著不同（校正后的p值小于0.001）。

请记住，报告的是校正后的 p值，以防止比较多组群时出现多重测试的问题。

如果您想比较所有组合的群体，可以使用TukeyHSD()函数并在which参数中指定交互作用：

# 所有性别和物种的组合TukeyHSD(mod,  which = "sex:species")

##   Tukey多重均值比较##     95%系列内置信度水平## ## 模型：aov(formula = body_mass_g ~ sex * species, data = dat)## ## $`sex:species`##                                      diff       lwr       upr     p adj## male:Adelie-female:Adelie        674.6575  527.8486  821.4664 0.0000000## female:Chinstrap-female:Adelie   158.3703  -25.7874  342.5279 0.1376213## male:Chinstrap-female:Adelie     570.1350  385.9773  754.2926 0.0000000## female:Gentoo-female:Adelie     1310.9058 1154.8934 1466.9181 0.0000000## male:Gentoo-female:Adelie       2116.0004 1962.1408 2269.8601 0.0000000## female:Chinstrap-male:Adelie    -516.2873 -700.4449 -332.1296 0.0000000## male:Chinstrap-male:Adelie      -104.5226 -288.6802   79.6351 0.5812048## female:Gentoo-male:Adelie        636.2482  480.2359  792.2606 0.0000000## male:Gentoo-male:Adelie         1441.3429 1287.4832 1595.2026 0.0000000## male:Chinstrap-female:Chinstrap  411.7647  196.6479  626.8815 0.0000012## female:Gentoo-female:Chinstrap  1152.5355  960.9603 1344.1107 0.0000000## male:Gentoo-female:Chinstrap    1957.6302 1767.8040 2147.4564 0.0000000## female:Gentoo-male:Chinstrap     740.7708  549.1956  932.3460 0.0000000## male:Gentoo-male:Chinstrap      1545.8655 1356.0392 1735.6917 0.0000000## male:Gentoo-female:Gentoo        805.0947  642.4300  967.7594 0.0000000

或者使用{agricolae}包中的HSD.test()函数，它用相同的字母表示不显着不同的子组：

library(agricolae)

HSD.test(mod,  trt = c("sex", "species"),  console = TRUE # print results)

## ## 研究：mod ~ c("sex", "species")## ## HSD Test for body_mass_g ## ## 均方误差：95726.69 ## ## sex:species,  means## ##                  body_mass_g      std  r  Min  Max## female:Adelie       3368.836 269.3801 73 2850 3900## female:Chinstrap    3527.206 285.3339 34 2700 4150## female:Gentoo       4679.741 281.5783 58 3950 5200## male:Adelie         4043.493 346.8116 73 3325 4775## male:Chinstrap      3938.971 362.1376 34 3250 4800## male:Gentoo         5484.836 313.1586 61 4750 6300## ## Alpha: 0.05 ; DF Error: 327 ## Studentized Range的临界值: 4.054126 ## ## Treatments with the same letter are not significantly different.## ##                  body_mass_g groups## male:Gentoo         5484.836      a## female:Gentoo       4679.741      b## male:Adelie         4043.493      c## male:Chinstrap      3938.971      c## female:Chinstrap    3527.206      d## female:Adelie       3368.836      d

如果您有很多组要比较，绘制它们可能更容易理解：

# 设置轴边距，以便标签不被切断par(mar = c(4.1, 13.5, 4.1, 2.1))

# 为每个比较创建置信区间plot(TukeyHSD(mod, which = "sex:species"),  las = 2 # 旋转x轴刻度)

从上面的输出和图中可以得出结论，除了女性Chinstrap和女性Adelie（p值=0.138）以及男性Chinstrap和男性Adelie（p值=0.581）之间的组合之外，所有性别和物种的组合都显着不同。

这些结果与上面显示的箱线图一致，并将在下面的可视化中得到证实，从而得出了R中的双向ANOVA的结论。

可视化

如果您想以不同于初步分析中已呈现的方式可视化结果，则以下是一些有用绘图的想法。

首先，使用{effects}包中的allEffects()函数显示每个子组的均值和标准误差：

# 方法1library(effects)

plot(allEffects(mod))

或使用{ggpubr}包：

# 方法2library(ggpubr)

ggline(subset(dat, !is.na(sex)), # 删除sex的NA水平  x = "species",  y = "body_mass_g",  color = "sex",  add = c("mean_se") # 添加均值和标准误) +  labs(y = "Mean of body mass (g)")

或者，使用{Rmisc}和{ggplot2}：

library(Rmisc)

# 按子组计算平均值和平均值的标准误差summary_stat <- summarySE(dat, measurevar = "body_mass_g", groupvars = c("species", "sex"))

# 绘制平均值和平均值的标准误差ggplot(subset(summary_stat, !is.na(sex)), aes(x = species, y = body_mass_g, colour = sex)) + geom_errorbar(aes(ymin = body_mass_g - se, ymax = body_mass_g + se), width = 0.1) + geom_point() + labs(y = "平均体重（克）")

其次，如果您只想绘制子组的平均值：

with(dat, interaction.plot(species, sex, body_mass_g))

最后但并非最不重要的，对于那些熟悉GraphPad的人，您可能已经熟悉了如下绘制平均值和误差线的方法：

# 作为条形图绘制平均值和平均值的标准误差ggplot(subset(summary_stat, !is.na(sex)), aes(x = species, y = body_mass_g, fill = sex)) + geom_bar(position = position_dodge(), stat = "identity") + geom_errorbar(aes(ymin = body_mass_g - se, ymax = body_mass_g + se), width = 0.25, position = position_dodge(.9)) + labs(y = "平均体重（克）")

R中的双向方差分析 AI 技术第16张

结论

在本文中，我们首先回顾了比较不同组之间定量变量的不同测试方法。然后，我们重点介绍了双向方差分析，从其目标和假设到在R中的实现，以及一些可视化和解释。我们还简要介绍了其基本假设和一种事后比较所有子组的测试方法。

所有这些都是使用{palmerpenguins}软件包中提供的penguins数据集进行说明的。

感谢您的阅读。

我希望本文能帮助您进行双向方差分析。

与往常一样，如果您有与本文主题相关的问题或建议，请在评论中提出，以便其他读者也能从讨论中受益。

理论上，单向方差分析也可以用于比较2组，而不仅仅是3组或更多组。然而，在实践中，通常会使用学生t检验比较2组，使用单向方差分析比较3组或更多组。使用独立样本的学生t检验和使用2组的单向方差分析得出的结论将是相似的。 ↩︎
请注意，如果不满足正态性假设，则可以应用许多转换来改善它，其中最常见的是对数转换（R中的log()函数）。 ↩︎
请注意，Bartlett的测试也适用于测试等方差性假设。 ↩︎
加法模型假定2个解释变量是独立的，它们不会相互作用。 ↩︎
其中mod是您保存的模型名称。 ↩︎
在这里，我们使用Benjamini和Hochberg（1995）校正，但您可以在几种方法之间选择。有关更多详细信息，请参见?p.adjust。 ↩︎

在R中进行方差分析
在R中和手动进行学生t检验：如何比较不同情况下的两组数据？
在R中进行单样本Wilcoxon检验
Kruskal-Wallis检验，或ANOVA的非参数版本
手动进行假设检验

最初发布于https://statsandr.com于2023年6月19日。