基于marginaleffects的边际效应与全模型森林图及tidy()格式转换

1. 正确计算分水平的边际效应 #

你当前的slopes()调用无需指定primary_outcome = 1（这是结果变量，datagrid()会自动处理协变量的典型值：分类变量取众数，连续变量取均值）。我们需要计算exposure_1在exposure_2=0和exposure_2=1时的边际效应，同时直接指数化得到与主模型一致的优势比（OR）：

# 计算exposure_2不同水平下exposure_1的边际效应，指数化得到OR
meffects <- slopes(logit_model_1,
                   variables = "exposure_1",  # 明确指定目标变量
                   newdata = datagrid(exposure_2 = c(0, 1)),
                   exponentiate = TRUE) %>%  # 直接生成OR值，和主模型匹配
  select(term, estimate, conf.low, conf.high, p.value, exposure_2)

# 调整格式，与主模型结果结构对齐
meffects_tidy <- meffects %>%
  mutate(
    model = case_when(
      exposure_2 == 0 ~ "Exposure 2 = 0",
      exposure_2 == 1 ~ "Exposure 2 = 1"
    ),
    exposure = "Exposure 1"
  ) %>%
  select(-exposure_2)

2. 合并主模型与边际效应结果 #

将主模型的tidy结果和边际效应的tidy结果合并，统一后续处理的数据集结构：

# 主模型结果（保留原逻辑，简化冗余步骤）
lcf <- logit_model_1 %>%
  tidy(exponentiate = TRUE, conf.int = TRUE, conf.level = 0.95) %>%
  mutate(model = "Main Result", exposure = "Exposure 1") %>%
  filter(term == "exposure_1") 

# 合并两个数据集
combined_table <- bind_rows(lcf, meffects_tidy)

3. 统一格式化数据（适配森林图） #

将合并后的数据集应用你原有的格式化规则，确保数值展示符合期刊要求：

library(stringr)
library(forcats)

lcftable <- combined_table %>%
  # 四舍五入并格式化数值
  mutate(across(c(estimate, conf.low, conf.high, p.value), ~ round(.x, 2))) %>%
  mutate(
    estimate = str_pad(estimate, width = 4, pad = "0", side = "right"),
    conf.low = str_pad(conf.low, width = 4, pad = "0", side = "right"),
    conf.high = str_pad(conf.high, width = 4, pad = "0", side = "right"),
    # 修正特殊值显示
    estimate = ifelse(estimate == "1.00", "1.00", estimate),
    conf.low = ifelse(conf.low == "1.00", "1.00", conf.low),
    conf.high = ifelse(conf.high == "2.00", "2.00", conf.high)
  ) %>%
  # 调整模型顺序，保证绘图时主结果在最上方
  mutate(model = fct_rev(fct_relevel(model, "Main Result", "Exposure 2 = 0", "Exposure 2 = 1"))) %>%
  # 生成标签和处理p值
  mutate(
    estimate_lab = paste0(estimate, " (", conf.low, " - ", conf.high, ")"),
    p.value = case_when(
      p.value < 0.001 ~ "<0.001",
      round(p.value, 2) == 0.05 ~ as.character(round(p.value, 3)),
      p.value < 0.01 ~ str_pad(as.character(round(p.value, 3)), width = 4, pad = "0", side = "right"),
      p.value > 0.995 ~ "1.00",
      TRUE ~ str_pad(as.character(round(p.value, 2)), width = 4, pad = "0", side = "right")
    ),
    # 转换回数值类型用于绘图
    estimate = as.numeric(estimate),
    conf.low = as.numeric(conf.low),
    conf.high = as.numeric(conf.high)
  )

4. 绘制包含边际效应的森林图 #

调整绘图代码，同时展示主模型结果和分水平的边际效应：

lcfplot1 <- lcftable %>%
  ggplot(aes(x = estimate, xmin = conf.low, xmax = conf.high, y = model)) +
  geom_pointrange() +
  geom_vline(xintercept = 1, linetype = "dashed", color = "gray50") +
  geom_point(size = 1.0) +
  facet_wrap(exposure~., ncol = 1, scale = "free_y") +
  xlab("Incident Odds Ratio, IOR (95% Confidence Interval)") +
  ggtitle("Associations between exposures and outcome") +
  geom_text(aes(x = 3.5, label = estimate_lab), hjust = 0, size = 3) +
  geom_text(aes(x = 3.5, y = model, label = paste0("p = ", p.value)), hjust = 0, size = 3, vjust = 1.5) +
  coord_cartesian(xlim = c(0, 7.0)) +
  scale_x_continuous(trans = "pseudo_log", breaks = c(1.0, 2.0, 5.0)) +
  theme(
    panel.background = element_blank(),
    panel.spacing = unit(1, "lines"),
    panel.border = element_rect(fill = NA, color = "white"),
    strip.background = element_rect(colour="white", fill="white"),
    strip.placement = "outside",
    text = element_text(size = 12),
    strip.text = element_text(face = "bold.italic", size = 11, hjust = 0.5),
    plot.title = element_text(face = "bold", size = 13, hjust = 0.5),
    axis.title.x = element_text(size = 10, vjust = -1),
    axis.title.y = element_blank(),
    panel.grid.major.x = element_line(color = "#D3D3D3", size = 0.3, linetype = 2),
    plot.margin = unit(c(1, 3, 0.5, 0.5), "inches")
  )

print(lcfplot1)

关键说明 #

• slopes()函数通过newdata = datagrid(exposure_2 = c(0,1))指定exposure_2的两个水平，自动控制其他协变量取典型值，无需手动设置结果变量。
• 加入exponentiate = TRUE直接生成优势比，与主模型的OR结果完全匹配，避免额外转换步骤。
• 合并数据集时保证model、exposure等列结构一致，确保森林图可以统一展示主结果和分水平的边际效应。

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