R - 协方差分析

R - 现代协方差分析 (ANCOVA)

协方差分析 (ANCOVA) 是一种统计方法，用于评估一个连续预测变量对连续响应变量的影响，同时控制一个或多个分类变量（称为协变量）的影响。它本质上是方差分析 (ANOVA) 和线性回归的结合。

例如，您可能想模拟汽车马力 (hp) 如何影响其油耗里程 (mpg)。然而，您怀疑变速箱类型（自动或手动）也扮演着重要角色。ANCOVA 允许您构建一个模型，其中包含 hp 作为连续预测变量和 am（变速箱类型）作为分类协变量，甚至可以测试 hp 和 mpg 之间的关系对于自动挡和手动挡汽车是否“不同”（这称为交互效应）。

现代设置：Tidyverse

虽然基础 R 具备完整功能，但现代 R 工作流通常会利用 tidyverse，这是一个为数据科学设计的包集合，它们共享相同的设计理念。我们将使用 dplyr 进行数据操作，使用 broom 清理模型输出。

# Install necessary packages if you haven't already
# install.packages(c("tidyverse", "broom"))

library(tidyverse)
library(broom)

步骤 1：数据准备

我们将使用内置的 mtcars 数据集。我们的第一步是选择感兴趣的变量，并确保我们的分类变量 (am) 正确地格式化为因子 (factor)。这会为其提供有意义的标签，并确保统计函数正确处理它。

# Load and prepare the data using dplyr pipes
input_data <- mtcars %>%
  select(mpg, hp, am) %>%
  mutate(am = factor(am, labels = c("Automatic", "Manual")))

# Display the first few rows of our prepared data
head(input_data)

执行上述代码后，将产生以下结果：

                   mpg  hp        am
Mazda RX4         21.0 110    Manual
Mazda RX4 Wag     21.0 110    Manual
Datsun 710        22.8  93    Manual
Hornet 4 Drive    21.4 110 Automatic
Hornet Sportabout 18.7 175 Automatic
Valiant           18.1 105 Automatic

步骤 2：ANCOVA 建模

我们将构建两个模型来检验我们的假设：

交互模型 (mpg ~ hp * am): 此模型测试 hp 和 mpg 之间的关系对于每种变速箱类型是否“不同”。* 符号自动包含 hp、am 及其交互项 (hp:am)。
平行斜率模型 (mpg ~ hp + am): 此模型假设 hp 对 mpg 的影响对于两种变速箱类型是“相同”的，但允许它们有不同的基线 mpg。

模型 1：包含交互项

让我们构建并检查包含交互项的模型。

# Create the ANCOVA model with an interaction term
model_interaction <- aov(mpg ~ hp * am, data = input_data)

# Use broom::tidy() for a clean summary table
summary_interaction <- tidy(model_interaction)
print(summary_interaction)

这将产生一个整洁的、基于 tibble 的输出：

# A tibble: 4 × 6
  term        df    sumsq   meansq statistic      p.value
  <chr>    <dbl>    <dbl>    <dbl>     <dbl>        <dbl>
1 hp           1   678.    678.       77.4   0.00000000150
2 am           1   202.    202.       23.1   0.0000475
3 hp:am        1     0.005   0.005     0.001  0.981
4 Residuals   28   245.      8.77      NA           NA

解释： hp 和 am 的 p.value 都非常小 (远小于 0.05)，表明它们是 mpg 的显著预测变量。然而，交互项 hp:am 的 p 值非常高 (0.981)，这表明马力对油耗里程的影响在自动挡和手动挡汽车之间没有显著差异。

模型 2：不包含交互项（平行斜率）

鉴于不显著的交互项，一个没有交互项的更简单模型可能更好。

# Create the model without an interaction term
model_parallel <- aov(mpg ~ hp + am, data = input_data)

# Get the tidy summary
summary_parallel <- tidy(model_parallel)
print(summary_parallel)

这提供了更简单模型的摘要：

# A tibble: 3 × 6
  term        df   sumsq  meansq statistic      p.value
  <chr>    <dbl>   <dbl>   <dbl>     <dbl>        <dbl>
1 hp           1   678.   678.       80.2   0.000000000763
2 am           1   202.   202.       23.9   0.0000346
3 Residuals   29   245.     8.46      NA           NA

解释： hp 和 am 都仍然是高度显著的预测变量。

步骤 3：比较模型

为了正式决定哪个模型更好，我们可以使用 anova() 函数来比较它们。这里非显著的 p 值意味着更简单的模型（没有交互项）是足够的。

# Compare the two models
anova(model_interaction, model_parallel)

输出证实了我们的猜测：

Analysis of Variance Table

Model 1: mpg ~ hp * am
Model 2: mpg ~ hp + am
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     28 245.43
2     29 245.44 -1 -0.0052515 0.0006 0.9806

Pr(>F) 为 0.9806，我们证实添加交互项并没有显著改善模型。我们应该选择更简单的 model_parallel。

最佳实践和后续步骤

检查假设： 对于任何线性模型，您都必须检查残差的正态性和方差齐性等假设。您可以使用 plot(model_parallel) 或更优雅地使用 ggplot2 来完成。
实际解释： 最终模型 mpg ~ hp + am 表明，在给定马力的情况下，手动挡汽车的平均 mpg 与自动挡汽车不同，并且这种差异在所有马力值范围内都是一致的。
更多资源： 对于更高级的建模，请探索 tidymodels 框架；对于事后分析，emmeans 包是行业标准。