基于lmerMod分层模型的ggplot2交互图：添加组间比较显著性标识

核心思路 #

通过emmeans高效提取模型边际均值（对应系数估计），结合ggsignif添加分组内的显著性对比标注，解决大数据集运算效率和可视化标注的问题。

步骤1：拟合模型并提取关键统计量 #

使用lmerTest拟合模型，借助emmeans快速获取交互项的边际均值及分组对比的显著性结果（比ggeffects更适配大数据集）：

# 加载依赖包
library(lmerTest)
library(ggplot2)
library(emmeans)
library(ggsignif)
library(dplyr)

# 生成模拟数据（复用示例代码）
set.seed(123)
N <- 1000
G1 <- 20
G2 <- 10

df <- data.frame(
  outcome = rnorm(N),
  cat1 = factor(rep(1:2, each = N / 2)),
  cat2 = factor(rep(c("Control", "T1", "T2", "T3"), each = N / 4)),
  G1 = factor(rep(1:G1, each = N / G1)),
  G2 = factor(rep(1:G2, each = N / G2))
)

# 拟合分层模型
fit <- lmerTest::lmer(outcome ~ cat1 * cat2 + (1|G1) + (1|G2), data = df)

# 提取cat1×cat2组合的边际均值（含置信区间）
emm <- emmeans(fit, ~ cat2 | cat1)
emm_df <- as.data.frame(emm)

# 在每个cat1水平内，做cat2的两两对比（Tukey校正p值）
pairs_emm <- pairs(emm, adjust = "tukey")
pairs_df <- as.data.frame(pairs_emm) %>%
  # 添加显著性星号规则
  mutate(
    signif = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01 ~ "**",
      p.value < 0.05 ~ "*",
      TRUE ~ ""
    )
  ) %>%
  # 过滤无显著性的对比，减少冗余标注
  filter(signif != "") %>%
  # 拆分对比组名称，方便后续定位
  separate(contrast, into = c("group1", "group2"), sep = " - ")

步骤2：绘制基础交互图（pointrange） #

用geom_pointrange实现要求的可视化样式，通过position_dodge确保同组内的cat2水平分开显示：

dodge_width <- 0.5 # 统一dodge宽度，保证图形对齐

base_plot <- ggplot(emm_df, aes(x = cat1, y = emmean, color = cat2)) +
  geom_pointrange(
    aes(ymin = lower.CL, ymax = upper.CL),
    position = position_dodge(width = dodge_width),
    size = 1
  ) +
  scale_color_brewer(palette = "Set1") +
  labs(x = "cat1", y = "模型系数估计值", color = "cat2") +
  theme_bw()

步骤3：添加显著性标注（括号+星号） #

计算标注的位置坐标，用ggsignif::geom_signif手动添加分组内的对比括号和星号：

# 为每个cat2水平分配x轴偏移量（匹配dodge宽度）
cat2_levels <- levels(emm_df$cat2)
x_offsets <- seq(-dodge_width/2, dodge_width/2, length.out = length(cat2_levels))
names(x_offsets) <- cat2_levels

# 完善标注的位置信息
pairs_df <- pairs_df %>%
  mutate(
    x = as.numeric(cat1), # 转换cat1为数值型，对应x轴位置
    x_start = x + x_offsets[group1], # 对比起始点x坐标
    x_end = x + x_offsets[group2], # 对比终点x坐标
    # 标注y轴位置：在置信区间最大值上方预留空间
    y_position = max(emm_df$upper.CL) + 0.1 * diff(range(emm_df$upper.CL, emm_df$lower.CL))
  )

# 合并标注到基础图
final_plot <- base_plot +
  geom_signif(
    data = pairs_df,
    aes(xmin = x_start, xmax = x_end, annotations = signif, y_position = y_position),
    textsize = 5, vjust = -0.2,
    manual = TRUE # 手动指定标注位置
  )

print(final_plot)

大数据集优化说明 #

• 替换ggeffects::hypothesis_test()为emmeans::pairs()：emmeans直接基于模型固定效应矩阵计算对比，避免冗余运算，在N=46000的数据集上效率提升明显。
• 若仍有速度瓶颈，可简化对比逻辑（如仅对比Control与其他处理组），减少运算量。

完整可运行代码 #

library(lmerTest)
library(ggplot2)
library(emmeans)
library(ggsignif)
library(dplyr)

# 生成模拟数据
set.seed(123)
N <- 1000
G1 <- 20
G2 <- 10

df <- data.frame(
  outcome = rnorm(N),
  cat1 = factor(rep(1:2, each = N / 2)),
  cat2 = factor(rep(c("Control", "T1", "T2", "T3"), each = N / 4)),
  G1 = factor(rep(1:G1, each = N / G1)),
  G2 = factor(rep(1:G2, each = N / G2))
)

# 拟合分层模型
fit <- lmerTest::lmer(outcome ~ cat1 * cat2 + (1|G1) + (1|G2), data = df)

# 提取边际均值及置信区间
emm <- emmeans(fit, ~ cat2 | cat1)
emm_df <- as.data.frame(emm)

# 计算分组内的两两对比及显著性
pairs_emm <- pairs(emm, adjust = "tukey")
pairs_df <- as.data.frame(pairs_emm) %>%
  mutate(
    signif = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01 ~ "**",
      p.value < 0.05 ~ "*",
      TRUE ~ ""
    )
  ) %>%
  filter(signif != "") %>%
  separate(contrast, into = c("group1", "group2"), sep = " - ")

# 计算标注位置
dodge_width <- 0.5
cat2_levels <- levels(emm_df$cat2)
x_offsets <- seq(-dodge_width/2, dodge_width/2, length.out = length(cat2_levels))
names(x_offsets) <- cat2_levels

pairs_df <- pairs_df %>%
  mutate(
    x = as.numeric(cat1),
    x_start = x + x_offsets[group1],
    x_end = x + x_offsets[group2],
    y_position = max(emm_df$upper.CL) + 0.1 * diff(range(emm_df$upper.CL, emm_df$lower.CL))
  )

# 绘制最终图形
ggplot(emm_df, aes(x = cat1, y = emmean, color = cat2)) +
  geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL), 
                  position = position_dodge(width = dodge_width), size = 1) +
  scale_color_brewer(palette = "Set1") +
  labs(x = "cat1", y = "模型系数估计值", color = "cat2") +
  theme_bw() +
  geom_signif(
    data = pairs_df,
    aes(xmin = x_start, xmax = x_end, annotations = signif, y_position = y_position),
    textsize = 5, vjust = -0.2,
    manual = TRUE
  )

内容的提问来源于stack exchange，提问作者CompSocialSciR