如何在R中为交互项添加参考水平（ggstats场景）

我之前在Stack Overflow提问「如何在gtsummary和sjPlot中为交互项添加参考水平」，但未得到满意的解决方案。现提供示例数据与模型：

set.seed(1000)
my_data <- rbind(
  data.frame(time = "Pre", treatment = "Control", response = rnorm(100, mean=1)),
  data.frame(time = "Pre", treatment = "Treatment", response = rnorm(100, mean=2)),
  data.frame(time = "Post", treatment = "Control", response = rnorm(100, mean=1)),
  data.frame(time = "Post", treatment = "Treatment", response = rnorm(100, mean=2))
) %>% 
  mutate(time = factor(time, levels = c("Pre", "Post"))) %>%
  mutate(treatment = factor(treatment, levels = c("Control", "Treatment")))

model3 <- lm(response ~ time * treatment, data = my_data)

我尝试使用ggstats::ggcoef_table(model3)，能显示因子变量的参考水平，但希望交互项也能显示参考水平（即timePre:treatmentControl这一组），让它出现在图表及旁侧的系数表中。请问能用ggstats或其他工具实现吗？

方法1：手动扩展模型结果后用ggstats绘图 #

ggstats本身没有直接添加交互项参考水平的参数，但可以先整理模型结果，手动插入参考组行，再传入绘图函数：

library(tidyverse)
library(ggstats)

# 整理模型结果，添加交互项参考水平行
tidy_model <- broom::tidy(model3, conf.int = TRUE) %>%
  # 插入参考组：timePre:treatmentControl，系数为0，置信区间固定为0
  add_row(
    term = "timePre:treatmentControl",
    estimate = 0,
    std.error = NA,
    statistic = NA,
    p.value = NA,
    conf.low = 0,
    conf.high = 0,
    .before = 1
  ) %>%
  # 重命名term列，优化显示文本
  mutate(
    term = case_when(
      term == "(Intercept)" ~ "timePre:treatmentControl (参考组)",
      term == "timePost" ~ "timePost:treatmentControl",
      term == "treatmentTreatment" ~ "timePre:treatmentTreatment",
      term == "timePost:treatmentTreatment" ~ "timePost:treatmentTreatment",
      TRUE ~ term
    )
  )

# 用整理后的数据集绘图
ggcoef_table(tidy_model, show_p_values = TRUE) +
  labs(title = "包含交互项参考水平的系数表")

方法2：用emmeans计算边际均值，可视化所有交互组 #

如果核心需求是展示所有交互组的差异对比，用emmeans计算每个交互组合的边际均值，再可视化会更直观：

library(emmeans)
library(ggplot2)

# 计算所有交互组的边际均值及置信区间
emm <- emmeans(model3, ~ time * treatment)
emm_tidy <- broom::tidy(emm, conf.int = TRUE)

# 绘制带置信区间的均值对比图
ggplot(emm_tidy, aes(x = interaction(time, treatment, sep = ":"), y = estimate)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  # 添加参考组的基准线
  geom_hline(yintercept = emm_tidy$estimate[emm_tidy$time == "Pre" & emm_tidy$treatment == "Control"], 
             linetype = "dashed", color = "red") +
  labs(x = "交互组", y = "响应均值", title = "所有交互组的响应均值对比") +
  theme_bw()

# 生成对应的统计表格
knitr::kable(emm_tidy, caption = "交互组边际均值及置信区间")

这种方式直接展示每个交互组的实际均值，参考组的均值作为基准线，更契合交互效应的展示需求。

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