如何在ggplot中重现sjPlot的Gamma模型三向交互图并优化样式

核心步骤说明 #

要重现sjPlot的交互图，关键是构建包含连续变量关键水平的预测数据集，再结合ggplot的美学调整实现自定义样式。以下是完整代码流程：

1. 准备预测数据集 #

首先将连续自变量转换为「均值、-1标准差、+1标准差」三个水平，同时控制协变量为参考水平（连续协变量取均值，分类协变量取参考类别）：

library(dplyr)
library(tidyr)
library(ggplot2)

# 假设数据集为df，训练好的Gamma GLM模型为model_gamma
# 提取连续自变量的均值和标准差
cont_vars_stats <- df %>%
  summarise(across(c(iv1, iv2), list(mean = ~mean(.x, na.rm = TRUE), sd = ~sd(.x, na.rm = TRUE))))

# 构建预测网格：
# - 分类自变量iv3保留所有水平
# - 协变量cov1/cov2设为均值（若为分类变量则替换为参考类别）
# - iv1取全范围（用于x轴绘制连续曲线），iv2取三个关键水平
pred_df <- df %>%
  select(iv3, cov1, cov2) %>%
  distinct() %>%
  mutate(across(c(cov1, cov2), ~mean(.x, na.rm = TRUE))) %>% # 控制协变量影响
  expand_grid(
    iv1 = seq(min(df$iv1, na.rm = TRUE), max(df$iv1, na.rm = TRUE), length.out = 100),
    iv2 = c(
      cont_vars_stats$iv2_mean - cont_vars_stats$iv2_sd,
      cont_vars_stats$iv2_mean,
      cont_vars_stats$iv2_mean + cont_vars_stats$iv2_sd
    )
  ) %>%
  # 给iv2的三个水平添加可读标签
  mutate(iv2_group = case_when(
    iv2 == cont_vars_stats$iv2_mean - cont_vars_stats$iv2_sd ~ "-1 SD",
    iv2 == cont_vars_stats$iv2_mean ~ "Mean",
    iv2 == cont_vars_stats$iv2_mean + cont_vars_stats$iv2_sd ~ "+1 SD"
  ))

2. 计算模型预测值与置信区间 #

针对Gamma模型，需使用type="response"将预测值转换回原始因变量尺度：

pred_df <- pred_df %>%
  mutate(
    pred_val = predict(model_gamma, newdata = ., type = "response"),
    # 计算标准误与95%置信区间
    se = predict(model_gamma, newdata = ., type = "response", se.fit = TRUE)$se.fit,
    lcl = pred_val - 1.96 * se,
    ucl = pred_val + 1.96 * se
  )

3. 绘制ggplot图并调整美学 #

实现分面、自定义背景/网格，同时匹配sjPlot的核心展示逻辑：

ggplot(pred_df, aes(x = iv1, y = pred_val, color = iv2_group)) +
  # 绘制预测曲线
  geom_line(linewidth = 1) +
  # 绘制置信区间带
  geom_ribbon(aes(ymin = lcl, ymax = ucl, fill = iv2_group), alpha = 0.2, color = NA) +
  # 按分类变量分面
  facet_wrap(~iv3) +
  # 设置白色背景、移除网格
  theme_bw() +
  theme(
    panel.grid = element_blank(), # 移除所有网格线
    plot.title = element_text(hjust = 0.5, size = 14),
    axis.title = element_text(size = 12)
  ) +
  # 自定义轴标签、标题与图例
  labs(
    x = "自变量1 (IV1)",
    y = "预测因变量 (DV)",
    color = "自变量2 (IV2) 水平",
    fill = "自变量2 (IV2) 水平",
    title = "Gamma GLM三向交互效应预测图"
  ) +
  # 按需调整轴范围
  xlim(c(min(df$iv1, na.rm = TRUE), max(df$iv1, na.rm = TRUE)))

关键细节提示 #

• 如果协变量是分类变量，将mutate(across(c(cov1, cov2), ~mean(.x, na.rm = TRUE)))替换为mutate(cov1 = factor("参考类别"))即可。
• 若需要区分曲线样式，可在geom_line()中添加linetype = iv2_group参数。
• 若需要更精准的置信区间，可使用boot包实现Bootstrap抽样计算。

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