如何用R的ggplot绘制二分类暴露与years变量的交互效应图

问题背景 #

已构建包含exposure*years交互项的线性混合效应模型，需用ggplot2绘制交互效应图，展示二分类暴露变量（0/1）在二分类变量years（0/1）不同水平下，对结局变量zatt_new的影响效应。

原模型代码：

library(ggplot2)
library(lme4)
library(lmerTest)
library(jtools)

set.seed(123)

# 创建数据集
n_participants <- 100
n_visits <- sample(2:5, n_participants, replace = TRUE)
n_rows <- sum(n_visits)
df <- data.frame(
  subnum = rep(1:n_participants, times = n_visits),
  visit = unlist(lapply(n_visits, seq_len)),
  zatt_new = rnorm(n_rows, mean = 0, sd = 1),
  exposure = sample(0:1, n_rows, replace = TRUE),
  years = sample(0:1, n_rows, replace = TRUE),
  BMI = rep(rnorm(n_participants, mean = 25, sd = 2), times = n_visits),
  Age = rep(rnorm(n_participants, mean = 75, sd = 5), times = n_visits)
)

# 构建混合效应模型
m <- lmer(zatt_new ~ exposure*years +BMI+Age+visit+(1+visit|subnum), data=df, na.action = na.omit)
summ(m)

解决方案 #

以下提供三种常用实现方式，可根据需求选择：

方法1：使用emmeans计算边际均值绘图（推荐） #

该方法会控制其他协变量（BMI、Age、visit）的影响，展示调整后的边际效应，更符合统计推断需求。

# 安装并加载emmeans包
if (!require(emmeans)) install.packages("emmeans")
library(emmeans)

# 计算exposure与years组合的边际均值（固定效应部分）
emm <- emmeans(m, ~ exposure * years, 
               at = list(visit = mean(df$visit),  # 控制visit为均值
                         BMI = mean(df$BMI),      # 控制BMI为均值
                         Age = mean(df$Age)))     # 控制Age为均值

# 转为数据框用于绘图
emm_df <- as.data.frame(emm)

# 将0/1转为易懂标签（可选）
emm_df$exposure <- factor(emm_df$exposure, levels = c(0,1), labels = c("未暴露", "暴露"))
emm_df$years <- factor(emm_df$years, levels = c(0,1), labels = c("基线", "随访"))

# 用ggplot2绘制交互图
ggplot(emm_df, aes(x = exposure, y = emmean, color = years, group = years)) +
  geom_point(position = position_dodge(width = 0.2), size = 3) +
  geom_line(position = position_dodge(width = 0.2)) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), 
                width = 0.1, position = position_dodge(width = 0.2)) +
  labs(x = "暴露状态", y = "zatt_new 预测均值", color = "时间点") +
  theme_bw() +
  theme(legend.position = "top")

方法2：使用jtools的plot_model快速绘图 #

如果你已经加载了jtools，可以直接用内置函数快速生成交互图，操作更简便：

# 绘制交互效应图，控制其他协变量为均值
plot_model(m, type = "pred", terms = c("exposure", "years"),
           pred.arg = list(at = list(visit = mean(df$visit),
                                     BMI = mean(df$BMI),
                                     Age = mean(df$Age))),
           axis.labels = c("暴露状态", "zatt_new 预测值"),
           legend.title = "时间点") +
  theme_bw()

方法3：手动构建预测数据集绘图（灵活定制） #

手动创建包含所有exposure和years组合的数据集，固定其他协变量为均值，再用predict获取预测值：

# 创建预测数据集：所有exposure和years组合，其他协变量取均值
pred_df <- expand.grid(
  exposure = c(0,1),
  years = c(0,1),
  BMI = mean(df$BMI),
  Age = mean(df$Age),
  visit = mean(df$visit),
  subnum = 1  # 随机效应取单个个体，不影响固定效应预测
)

# 获取固定效应预测值（忽略随机效应，用re.form=NA）
pred_df$pred_zatt <- predict(m, newdata = pred_df, re.form = NA)

# 获取置信区间（使用lme4的模拟方法）
set.seed(123)
pred_ci <- predict(m, newdata = pred_df, re.form = NA, interval = "confidence")
pred_df <- cbind(pred_df, pred_ci)

# 转换标签
pred_df$exposure <- factor(pred_df$exposure, levels = c(0,1), labels = c("未暴露", "暴露"))
pred_df$years <- factor(pred_df$years, levels = c(0,1), labels = c("基线", "随访"))

# 绘图
ggplot(pred_df, aes(x = exposure, y = fit, color = years, group = years)) +
  geom_point(position = position_dodge(width = 0.2), size = 3) +
  geom_line(position = position_dodge(width = 0.2)) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                width = 0.1, position = position_dodge(width = 0.2)) +
  labs(x = "暴露状态", y = "zatt_new 预测值", color = "时间点") +
  theme_bw()

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