用R的marginaleffects/emmeans包计算风险比比值及置信区间

下面提供两种实现方法，分别使用marginaleffects和emmeans包：

方法一：使用marginaleffects包直接检验假设 #

利用hypotheses参数可直接指定比值关系，自动生成估计值与置信区间：

# 先获取分性别的边际风险比
comp <- marginaleffects::avg_comparisons(model, comparison="ratio",
                                         variables="condition", by="sex")

# 定义假设：男性RR / 女性RR，计算该比值的估计值和置信区间
hyp <- marginaleffects::hypotheses(comp, "b1 / b2")

# 查看结果
hyp

也可直接在avg_comparisons中通过假设语法指定：

marginaleffects::avg_comparisons(
    model,
    comparison = "ratio",
    variables = "condition",
    by = "sex",
    hypotheses = "condition_1 - condition_0 for sex=1 / condition_1 - condition_0 for sex=0"
)

方法二：使用emmeans包计算 #

通过对数转换的线性性质计算比值的置信区间，步骤如下：

library(emmeans)

# 1. 获取各(condition, sex)组合的预测概率（risk尺度）
emm <- emmeans(model, ~ condition * sex, type = "response")

# 2. 计算分性别的风险比（实验组/对照组）
rr <- pairs(emm, by = "sex", type = "response", ratio = TRUE)

# 3. 提取估计值和方差协方差矩阵，计算比值及其置信区间
rr_est <- coef(rr)
vcov_rr <- vcov(rr)

# 对数转换后做线性组合，再转换回原始尺度
log_ratio <- log(rr_est[1]) - log(rr_est[2])
se_log_ratio <- sqrt(vcov_rr[1,1] + vcov_rr[2,2] - 2*vcov_rr[1,2])

ratio_est <- exp(log_ratio)
ci_lower <- exp(log_ratio - 1.96 * se_log_ratio)
ci_upper <- exp(log_ratio + 1.96 * se_log_ratio)

# 输出结果
data.frame(
    ratio = ratio_est,
    ci_lower = ci_lower,
    ci_upper = ci_upper
)

补充说明 #

• 两种方法均基于对数转换的线性特性：风险比比值的对数等于两个风险比对数的差，以此简化方差计算，再转换回原始尺度得到置信区间。
• 若需更精确的似然置信区间，可在emmeans的pairs结果上执行contrast(rr, list(ratio = c(1, -1)), method = "profile")。

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