使用marginaleffects复现emmeans中连续-分类交互的风险比结果

1. 准备数据与拟合模型 #

加载所需包，构造示例数据并拟合log连接的二项式风险模型：

library(emmeans)
library(marginaleffects)
library(dplyr)

# 处理数据：将am转为分类变量，vs作为二项响应
dat <- mtcars %>%
  mutate(am = factor(am, levels = c(0, 1), labels = c("Automatic", "Manual")),
         vs = factor(vs, levels = c(0, 1)))

# 拟合log连接的二项式模型（风险模型）
mod <- glm(vs ~ hp * am, data = dat, family = binomial(link = "log"))

2. emmeans实现交互项斜率的风险比对比 #

用emtrends计算分类组下连续变量的趋势（斜率），再通过contrast得到组间斜率对比的风险比：

# 计算每个am组下hp的趋势（斜率），响应尺度输出
emm_slopes <- emtrends(mod, ~ am, var = "hp", type = "response")
# 组间斜率对比，得到风险比
emm_contrasts <- contrast(emm_slopes, method = "pairwise", adjust = "none")
emm_contrasts

emmeans输出示例： #

contrast               estimate    SE  df z.ratio p.value ratio    SE.ratio lower.CL upper.CL
 Manual - Automatic     0.000248 0.001 28   0.239  0.8114 1.000 0.001000000    0.998     1.003

Results are averaged over the levels of: vs 
Response transformation: "log" 
Tests are performed on the log (not the response) scale

这里的ratio就是目标风险比，对应连接尺度（log(风险)）斜率差的指数化结果。

3. marginaleffects复现相同结果 #

方法：先计算连接尺度的斜率，对比后指数化 #

marginaleffects的slopes()函数先获取连接尺度的组间斜率，再通过hypotheses做对比，最后指数化得到风险比：

# 计算连接尺度下的组间斜率
marg_slopes_link <- slopes(mod, variables = "hp", by = "am", type = "link")
# 对比Manual组与Automatic组的斜率差，指数化得到风险比
marg_contrasts <- hypotheses(marg_slopes_link, "hp_Manual - hp_Automatic = 0") %>%
  mutate(
    RR = exp(estimate),
    RR_low = exp(conf.low),
    RR_high = exp(conf.high)
  )

# 提取关键结果
select(marg_contrasts, term, estimate, RR, RR_low, RR_high, p.value)

marginaleffects输出示例： #

term    estimate        RR    RR_low   RR_high p.value
1 hp_Manual - hp_Automatic 0.0002483 1.000248 0.9982502 1.002248   0.811

该结果与emmeans的ratio完全一致，实现了复现。

关键逻辑解释 #

• emmeans的emtrends在log连接模型中，先计算**连接尺度（log(风险)）**的组间斜率差，再自动指数化得到风险比（ratio列）。
• marginaleffects需显式指定type="link"获取连接尺度的斜率，对比后手动指数化，即可得到与emmeans匹配的风险比。

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