lvzhaopeng的数字花园

R语言那些事 统计系列（十三） rstatix ANCOVA 去除协变量分析

12分钟阅读

缘由

统计是生信的基石之一

碎碎念

要爱一个具体的人，不要爱一个抽象的人； 要爱生活，不要爱生活的意义。 ——《卡拉马佐夫兄弟》

Ref

ANCOVA in R: The Ultimate Practical Guide - Datanovia
当我们想比较不同处理组在某个结局变量上的差异时，如果存在一个或多个与结局相关但非感兴趣的变量（协变量，covariate），我们可以用 ANCOVA 来在控制协变量影响后检验组效应。相比单纯 ANOVA，ANCOVA 能增加检验功效、减少误差变异并提供更精确的组间比较。

本文通过两个示例演示：

  • 单因素 ANCOVA（一个处理组变量 + 一个协变量）
  • 双向 ANCOVA（两个处理因子 + 一个协变量，并检验交互）

二、示例一：单因素 ANCOVA（anxiety 数据）

library(tidyverse)
library(ggpubr)
library(rstatix)
library(broom)
 
 
data("anxiety", package = "datarium")
anxiety <- anxiety %>%
	select(id, group, t1, t3) %>%
	rename(pretest = t1, posttest = t3)
anxiety[14, "posttest"] <- 19
# Inspect the data by showing one random row by groups
set.seed(123)
anxiety %>% sample_n_by(group, size = 1)
 
 
 
ggscatter(
	anxiety, x = "pretest", y = "posttest",
	color = "group", add = "reg.line"
)+
	stat_regline_equation(
		aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = group)
	)
 
#检查斜率的均匀性 是否有互作; 由于交互作用项不具有统计学意义，因此回归斜率存在同质性；   若回归线斜率差异较大，可能存在 `group * pretest` 的交互（即斜率不齐）。这是 ANCOVA 的关键假设之一：**协变量与因子之间的斜率应当同质**（homogeneity of regression slopes）。当交互项显著（p < .05）时，不能用简单 ANCOVA（因为协变量在不同组的影响不同），需要改用分组回归或报告含交互项的模型。
anxiety %>% anova_test(posttest ~ group*pretest)
 
 
#残差的正态性
model <- lm(posttest ~ pretest + group, data = anxiety)
# Inspect the model diagnostic metrics
model.metrics <- augment(model) %>%
	select(-.hat, -.sigma, -.fitted) # Remove details
head(model.metrics, 3)
shapiro_test(model.metrics$.resid)
 
#方差的同质性
 
model.metrics %>% levene_test(.resid ~ group)
 
#异常值
model.metrics %>%
	filter(abs(.std.resid) > 3) %>%
	as.data.frame()
 
 
##在输入主变量之前　需要消除协变量的影响 重点看 `group` 的 F 和 p 值：如果显著，说明在控制 pretest 后，`group` 仍对 posttest 有显著影响（即处理有效果）。
 
res.aov <- anxiety %>% anova_test(posttest ~ pretest + group)
get_anova_table(res.aov)
 
 
#事后检验  在控制协变量后计算组间均值比较（估计的边际均值，adjusted means），并做多重比较校正（此处用 Bonferroni）
library(emmeans)
pwc <- anxiety %>%
	emmeans_test(
		posttest ~ group, covariate = pretest,
		p.adjust.method = "bonferroni"
	)
pwc
 
 
get_emmeans(pwc) #返回每组的调整后均值。
 
 
pwc <- pwc %>% add_xy_position(x = "group", fun = "mean_se")
ggline(get_emmeans(pwc), x = "group", y = "emmean") +
	geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
	stat_pvalue_manual(pwc, hide.ns = TRUE, tip.length = FALSE) +
	labs(
		subtitle = get_test_label(res.aov, detailed = TRUE),
		caption = get_pwc_label(pwc)
	)

  1. 事后检验（Pairwise comparisons）

    • emmeans_test()：。
    • get_emmeans(pwc) 图上用 ggline() 绘制均值与置信区间，并用 stat_pvalue_manual() 添加显著性标注。

三、双向 ANCOVA（stress 数据）

#### 双向ANCOVA
data("stress", package = "datarium")
stress %>% sample_n_by(treatment, exercise) #快速检查变量类型与分布
 
 
#线性假设
ggscatter(
	stress, x = "age", y = "score",
	facet.by  = c("exercise", "treatment"),
	short.panel.labs = FALSE
)+
	stat_smooth(method = "loess", span = 0.9)
 
 
#回归斜率的均匀性
stress %>%
	anova_test(
		score ~ age + treatment + exercise +
			treatment*exercise + age*treatment +
			age*exercise + age*exercise*treatment
	)
 
stress %>%
	unite(col = "group", treatment, exercise) %>%
	anova_test(score ~ group*age)
 
 
#残差的正态性
model <- lm(score ~ age + treatment*exercise, data = stress)
# Inspect the model diagnostic metrics
model.metrics <- augment(model) %>%
	select(-.hat, -.sigma, -.fitted, -.se.fit) # Remove details
head(model.metrics, 3)
shapiro_test(model.metrics$.resid)
#方差的同质性
#levene_test(.resid ~ treatment*exercisemodel.metrics)
 
model.metrics %>%
	filter(abs(.std.resid) > 3) %>%
	as.data.frame()
 
#调整年龄后，治疗和运动对压力评分存在统计学意义的交互作用
res.aov <- stress %>%
	anova_test(score ~ age + treatment*exercise)
get_anova_table(res.aov)
 
#事后测试
#1. 简单主效应分析
stress %>%
	group_by(exercise) %>%
	anova_test(score ~ age + treatment)
 
#2 成对比较 Pairwise comparisons
pwc <- stress %>%
	group_by(exercise) %>%
	emmeans_test(
		score ~ treatment, covariate = age,
		p.adjust.method = "bonferroni"
	)
pwc %>% filter(exercise == "high")
 
 
lp <- ggline(
	get_emmeans(pwc), x = "exercise", y = "emmean",
	color = "treatment", palette = "jco"
) +
	geom_errorbar(
		aes(ymin = conf.low, ymax = conf.high, color = treatment),
		width = 0.1
	)
 
 
# Comparisons between treatment group at each exercise level
pwc <- pwc %>% add_xy_position(x = "exercise", fun = "mean_se", step.increase = 0.2)
pwc.filtered <- pwc %>% filter(exercise == "high")
lp +
	stat_pvalue_manual(
		pwc.filtered, hide.ns = TRUE, tip.length = 0,
		bracket.size = 0
	) +
	labs(
		subtitle = get_test_label(res.aov,  detailed = TRUE),
		caption = get_pwc_label(pwc)
	)
 
pwc2 <- stress %>%
	group_by(treatment) %>%
	emmeans_test(
		score ~ exercise, covariate = age,
		p.adjust.method = "bonferroni"
	) %>%
	select(-df, -statistic, -p) # Remove details
pwc2 %>% filter(treatment == "yes")
# Comparisons between exercises group at each treatment level
pwc2 <- pwc2 %>% add_xy_position(x = "exercise", fun = "mean_se")
pwc2.filtered <- pwc2 %>% filter(treatment == "yes")
lp +
	stat_pvalue_manual(
		pwc2.filtered, hide.ns = TRUE, tip.length = 0,
		step.group.by = "treatment", color = "treatment"
	) +
	labs(
		subtitle = get_test_label(res.aov,  detailed = TRUE),
		caption = get_pwc_label(pwc2)
	)

  1. 事后检验的两种视角

    • 简单主效应分析（例如在每个 exercise 水平内检验 treatment）：用 group_by(exercise) 后做 anova_test(score ~ age + treatment)。
    • 成对比较（pairwise）：对每个 exercise 层次，用 emmeans_test(score ~ treatment, covariate = age) 计算调整后的均值并进行多重比较校正（这里用 Bonferroni）。结果以图 lp 展示，使用 stat_pvalue_manual() 标注显著性。类似地，也可以把比较方向换成“在每个 treatment 水平比较 exercise”。

四、疑问

1. anova_test(score ~ age + treatment*exercise) 是做什么的？**

这一步做的是 整体模型检验（overall F-test）。
它回答的问题是：

“在控制（调整）年龄 age 的影响后，treatment、exercise 以及它们的交互作用是否显著影响 score？”

score ~ age + treatment*exercise
等价于：score ~ age + treatment + exercise + treatment:exercise

anova_test() 给出的表包含：

项含义
age协变量 age 是否显著影响 score
treatment在调整 age 后，treatment 主效应是否显著
exercise同上
treatment:exercise在 age 调整后，treatment 在不同 exercise 下是否有不同效果

2. 为什么看起来像普通 ANOVA？ANCOVA 的“协变量调整”体现在哪里？**

ANCOVA 的核心作用 各组平均值已经被协变量（age）调整过了。

控制年龄差异 → 看 treatment/exercise 是否仍有影响。

举例：

  • 如果“年龄越大压力越大（score 越高）”
  • 那么在 ANCOVA 中，所有组的平均压力都会被“校正”为相同年龄水平的估计值

因此 ANCOVA 的结果不是原始均值，而是：

在相同年龄（通常是总体均值）下，各组的预测得分（adjusted means, 也叫 LS-means）。

这种“校正后比较组差异”的能力，就是 ANCOVA 的核心价值。

普通 ANOVA 比较的是“原始均值”。
ANCOVA 比较的是“age 调整后的均值”（这才是公平比较）。

anova_test(score ~ age + treatment*exercise)

3. 什么时候用 emmeans？什么时候用其他方法做事后检验？**

**1. 用 emmeans 的情况 **只要你用了协变量（ANCOVA）→ 就应该用 emmeans。

emmeans 会自动计算：

  • 调整后的均值（LS-means）
  • 考虑协变量的比较
  • 多重比较校正
  • 交互作用下的简单主效应（simple main effects）

它能做 ANCOVA/ANOVA/线性模型中的所有常见事后检验，而且特别适合：

  • 有协变量（ANCOVA）
  • 有交互作用
  • 需要 adjusted means
  • 需要 pairwise 统计推断

2. 不用 emmeans 的情况（可用 TukeyHSD、pairwise.t.test 等）

以下情形可用传统 post-hoc：

  • 只有 ANOVA（没有协变量）
  • 不需要 adjusted means
  • 不关心交互作用下的 simple effects
  • 只是想比较原始组均值

五 、ANOVA vs ANCOVA 结果如何不同？

1. 生成示例数据：协变量不均衡，但与因变量高度相关

set.seed(123)
 
n <- 80
age <- rnorm(n, 40, 10)
 
# group A 年龄略小， group B 年龄略大 → 典型需要 ANCOVA 的情形
group <- factor(rep(c("A", "B"), each = n/2))
age[group=="A"] <- age[group=="A"] - 5
age[group=="B"] <- age[group=="B"] + 5
 
# true model: score = age effect + small group effect
score <- 0.7*age + ifelse(group=="B", 3, 0) + rnorm(n, 0, 6)
 
dat <- data.frame(age, group, score)
head(dat)

特点：

  • age 强烈影响 score（真实模型）
  • A 和 B 年龄不平衡（模拟真实研究中常见混杂）

2. 做 ANOVA（不调整 age）

anova_model <- aov(score ~ group, data = dat)
summary(anova_model)

解释：

  • 由于 group B 年龄更大，它的 score 偏高 → 造成“看起来”组间差异显著
  • 但这其实是年龄差造成的假效应

3. 做 ANCOVA（加入 covariate age）

ancova_model <- aov(score ~ age + group, data = dat)
summary(ancova_model)
分析方式解释
ANOVA（忽略 age）看似显著（假象！年龄混杂）
ANCOVA（控制 age）完全不显著（真实模型）

把由于协变量 age 引起的错误组间差异“纠正”掉，使我们检验的只有组别本身的真实影响。

为什么 F 值会改变？

F = 组间差异 / 组内误差

  • 当不加入 age（ANOVA）
    ➜ age 的影响被当成“组间差异”或“误差”，导致 inflated F
  • 当加入 age（ANCOVA）
    ➜ age 对 score 的解释被从误差中提取出来
    ➜ 误差项变小
    ➜ 组间的真实差异变得清晰（通常更小）
    因此 F 值和 p 值都会发生显著变化。

1) ANOVA 比较的是“原始均值”，容易被协变量误导。
2) ANCOVA 通过控制协变量比较“调整后均值”。
3) 是否加入协变量会改变 F 值与结论，因为误差项不同。

关系图谱

反向链接