R语言中geeglm()模型结合margins()计算边际比值比的问题

问题核心 #

使用geeglm()构建含交互项的GEE模型后，需计算exposure_1取0、1时的边际比值比（OR），并与exposure_1的总OR在森林图中对比，但margins()函数不支持geeglm模型，需手动实现边际效应计算与绘图整合。

解决方案 #

1. 手动计算边际OR及置信区间 #

对于含exposure_1*exposure_2交互项的模型，边际OR需基于其他协变量的分布（此处用样本均值）计算：

• 构建exposure_1=1与exposure_1=0时的logit预测值差值的线性组合
• 通过模型系数协方差矩阵计算该差值的标准误，再指数化得到边际OR及95%置信区间

2. 完整实现代码 #

library(ggplot2)
library(dplyr)
library(geepack)
library(broom)
library(forcats)

# 生成模拟数据（与原代码一致）
set.seed(13)
n <- 100
exposure_1 <- sample(0:1, n, replace = TRUE)
exposure_2 <- sample(0:1, n, replace = TRUE)
exposure_3 <- sample(0:1, n, replace = TRUE)
primary_outcome <- rbinom(n, 1, 0.3 + 0.2 * exposure_1 + 0.1 * exposure_2 - 0.1 * exposure_3)
enr_community <- sample(letters[1:5], n, replace = TRUE)
df <- data.frame(
  exposure_1 = exposure_1,
  exposure_2 = exposure_2,
  exposure_3 = exposure_3,
  primary_outcome = primary_outcome,
  enr_community = enr_community
)

# 构建GEE模型
logit_model_1 <- geeglm(formula = primary_outcome ~ exposure_1*exposure_2 + exposure_3 , 
                        family = binomial, id = enr_community, corstr = "independence", data = df)

# 提取exposure_1的总OR（主效应）
total_or <- logit_model_1 %>%
  tidy(exponentiate = TRUE, conf.int = TRUE, conf.level = 0.95) %>%
  filter(term == "exposure_1") %>%
  mutate(
    model = "总OR（exposure_1主效应）",
    exposure = "Exposure 1"
  )

# 计算协变量均值
cov_means <- df %>%
  summarise(
    exposure_2_mean = mean(exposure_2),
    exposure_3_mean = mean(exposure_3)
  )

# 获取模型系数与协方差矩阵
coefs <- coef(logit_model_1)
vcov_mat <- vcov(logit_model_1)

# 计算整体边际OR（exposure_1=1 vs 0，基于协变量均值）
linear_comb <- c(0, 1, cov_means$exposure_2_mean, cov_means$exposure_2_mean, cov_means$exposure_3_mean)
marginal_logor <- sum(linear_comb * coefs)
se_marginal <- sqrt(t(linear_comb) %*% vcov_mat %*% linear_comb)
marginal_or_df <- tibble(
  term = "marginal_exposure_1",
  estimate = exp(marginal_logor),
  conf.low = exp(marginal_logor - 1.96 * se_marginal),
  conf.high = exp(marginal_logor + 1.96 * se_marginal),
  p.value = 2 * pnorm(abs(marginal_logor / se_marginal), lower.tail = FALSE),
  model = "边际OR（exposure_1=1 vs 0）",
  exposure = "Exposure 1"
)

# 计算分层边际OR（按exposure_2取值拆分）
# exposure_2=0时的边际OR
linear_comb_0 <- c(0, 1, 0, 0, cov_means$exposure_3_mean)
logor_0 <- sum(linear_comb_0 * coefs)
se_0 <- sqrt(t(linear_comb_0) %*% vcov_mat %*% linear_comb_0)
# exposure_2=1时的边际OR
linear_comb_1 <- c(0, 1, 1, 1, cov_means$exposure_3_mean)
logor_1 <- sum(linear_comb_1 * coefs)
se_1 <- sqrt(t(linear_comb_1) %*% vcov_mat %*% linear_comb_1)

stratified_marginal_df <- tibble(
  term = c("marginal_exposure_1_exposure2_0", "marginal_exposure_1_exposure2_1"),
  estimate = c(exp(logor_0), exp(logor_1)),
  conf.low = c(exp(logor_0 - 1.96 * se_0), exp(logor_1 - 1.96 * se_1)),
  conf.high = c(exp(logor_0 + 1.96 * se_0), exp(logor_1 + 1.96 * se_1)),
  p.value = c(2 * pnorm(abs(logor_0 / se_0), lower.tail = FALSE),
              2 * pnorm(abs(logor_1 / se_1), lower.tail = FALSE)),
  model = c("边际OR（exposure_2=0时）", "边际OR（exposure_2=1时）"),
  exposure = "Exposure 1"
)

# 合并所有OR数据
all_or_data <- bind_rows(total_or, marginal_or_df, stratified_marginal_df)

# 格式化数据（匹配原代码的期刊格式要求）
lcftable <- all_or_data %>%
  mutate(across(
    c(estimate, conf.low, conf.high, p.value),
    ~ str_pad(
      round(.x, 2),
      width = 4,
      pad = "0",
      side = "right"
    )
  )) %>%
  mutate(
    estimate = case_when(estimate == "1000" ~ "1.00", TRUE ~ estimate),
    conf.low = case_when(conf.low == "1000" ~ "1.00", TRUE ~ conf.low),
    conf.high = case_when(conf.high == "2000" ~ "2.00", TRUE ~ conf.high)
  ) %>%
  mutate(model = fct_rev(fct_relevel(model, "总OR（exposure_1主效应）", "边际OR（exposure_1=1 vs 0）", "边际OR（exposure_2=0时）", "边际OR（exposure_2=1时）"))) %>%
  mutate(
    estimate = format(as.numeric(estimate), n.small = 3),
    conf.low = format(as.numeric(conf.low), n.small = 3),
    conf.high = format(as.numeric(conf.high), n.small = 3),
    estimate_lab = paste0(estimate, " (", conf.low, " - ", conf.high, ")")
  ) %>%
  mutate(p.value = as.numeric(p.value)) %>%
  mutate(p.value = case_when(
    p.value < .001 ~ "<0.001",
    round(p.value, 2) == .05 ~ as.character(round(p.value, 3)),
    p.value < .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")
  )) %>%
  mutate(
    estimate = as.numeric(estimate),
    conf.low = as.numeric(conf.low),
    conf.high = as.numeric(conf.high)
  )

# 绘制整合后的森林图
lcfplot1 <- lcftable %>%
  ggplot(mapping = 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.2) +
  facet_wrap(exposure~., ncol = 1, scale = "free_y") +
  xlab("比值比（OR）及95%置信区间") +
  ggtitle("Exposure 1与结局的关联：总OR与边际OR对比") +
  geom_text(aes(x = max(conf.high)*1.1, label = estimate_lab), hjust = 0, size = 3) +
  geom_text(aes(x = max(conf.high)*1.3, label = paste0("P=", p.value)), hjust = 0, size = 3) +
  scale_x_continuous(trans = "pseudo_log", breaks = c(0.5, 1.0, 1.5, 2.0, 3.0)) +
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(fill = NA, color = "black"),
    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),
    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, 4, 0.5, 0.5), "inches")
  ) +
  geom_errorbarh(height = 0.2)

# 输出森林图
print(lcfplot1)

代码说明 #

• 边际OR计算：通过线性组合模型系数，结合协变量均值得到边际logOR，利用协方差矩阵推导标准误，最终指数化得到OR及置信区间
• 分层边际OR：针对exposure_2的不同取值单独计算exposure_1的边际效应，直观展示交互项的影响
• 森林图整合：将总OR、整体边际OR、分层边际OR合并后，沿用原代码的格式化逻辑，在同一图中实现多结果对比

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