基于R中brms包构建两组比例差值后验分布的问题咨询

核心问题答复 #

• 现有代码无法正确生成两组率差的后验分布，存在系数含义误解、非线性变换错误、先验尺度错配三类问题，具体错误点见后续拆解。
• 不需要借助后验预测相关的特殊函数，只要正确理解模型参数的尺度含义，对后验样本做正确变换即可得到目标分布；如果怕手动计算出错，用tidybayes包的条件效应计算函数也可以快速得到正确结果。
• 支持负值的效应量先验不要在0-1有界的概率尺度用Beta分布设置：如果在logit尺度建模组间效应，直接用正态/学生t分布即可天然支持正负取值；如果要直接在率差（Risk Difference, RD）尺度设置先验，用截断在[-1,1]区间的正态分布是可解释性最强、最常用的方案。

现有代码的问题与参数含义解释 #

你对as_draws_df输出的参数含义理解存在偏差，这是后续计算错误的核心原因：

1. 参数含义说明
   模型公式event | trials(n) ~ group中，R会默认按字母序对分组因子做编码，由于con首字母在字母表中排在int之前，对照组con会被设为因子参考水平：

   • b_Intercept是对照组logit尺度的事件率，即log(p_con/(1-p_con))
   • b_groupint是干预组与对照组的logit尺度差值，也就是对数优势比（log OR），计算公式为log(p_int/(1-p_int)) - log(p_con/(1-p_con))，取值范围为全体实数，天然可以为负。

2. 核心计算错误
   逆logit变换是非线性变换，不满足inv_logit(a-b) = inv_logit(a) - inv_logit(b)的运算规则。你当前代码中单独对b_groupint做逆logit变换得到的b_groupint_p没有任何实际统计含义——它既不是干预组的事件率，也不是两组率差，基于这个值计算的new_p = b_Intercept_p - b_groupint_p完全不是目标的两组率差后验。

3. 先验设置错误

   • 尺度错配：你通过beta.select选的是概率尺度的分位点，却把logit变换后的数值直接赋值给logit尺度的回归系数，两个尺度的数值没有直接对应关系，最终得到的先验完全不符合你最初的设定预期。
   • 赋值风险：你靠行索引修改priors矩阵的方式兼容性很差，不同brms版本返回的先验表行顺序可能变化，容易把先验错设到其他参数上。
   • 虽然Beta分布做logit变换后理论上可以取到负值，但由于你设置的Beta分位点集中在0-0.1的正区间，变换后的先验完全集中在负域，相当于强制预设干预只会大幅降低事件率，和你原本想设的“率差中位数0.05、90%分位点0.1”的预期完全相反。

正确实现方案 #

关键步骤说明 #

1. 固定因子参考水平：手动把对照组设为因子第一水平，避免默认排序带来的参数含义混乱。
2. 正确设置先验：
   • 对照组logit率的先验可以保留你原来的分位点匹配逻辑，因为率本身在0-1区间，不需要支持负值。
   • 组间效应先验直接在logit尺度设无界的正态分布即可支持负值，比如设置normal(0, 0.7)，对应log OR的95%区间大概在-1.4到1.4之间，对应率差的范围大概在-0.3到0.3，适合弱信息先验场景；如果有明确的先验信息，可以调整正态分布的均值和标准差，之后跑先验预测检验确认先验下的率差分布符合认知即可。
3. 正确计算率差后验：
   先在logit尺度计算每组的线性预测值，再统一做逆logit变换得到每组的后验事件率，最后做差得到率差：
   • 对照组后验率：inv_logit(b_Intercept)
   • 干预组后验率：inv_logit(b_Intercept + b_groupint)
   • 率差后验：inv_logit(b_Intercept + b_groupint) - inv_logit(b_Intercept)

修正后可运行代码 #

# 加载包
library(tidyverse)
library(ProbBayes)
library(brms)
library(tidybayes)

# 采样设置
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
n.chains = 3
n.iter = 4000
n.warmup = 1000 # 原warmup=500偏短，调整为1000保证采样收敛
n.draws = (n.iter - n.warmup) * n.chains

# 试验数据
con.n = 150
con.event = 50
int.n = 150
int.event = 30

# 构造数据集，手动设置对照组为因子参考水平
data = data.frame(
  group = factor(c("int", "con"), levels = c("con", "int")),
  n = c(int.n, con.n),
  event = c(int.event, con.event)
)

# 自定义变换函数
fun_logit = function(x) log(x/(1-x))
fun_invlogit = function(x) exp(x)/(1+exp(x))

# 生成对照组率的先验（匹配原分位点设定：中位数0.5概率取0.4，90%概率小于0.6）
beta0.val = beta.select(list(x = 0.4, p = 0.5), list(x = 0.6, p = 0.9))
p0_prior_sim = rbeta(n.draws, beta0.val[1], beta0.val[2])
theta0_prior_mean = mean(fun_logit(p0_prior_sim))
theta0_prior_sd = sd(fun_logit(p0_prior_sim))

# 设置组间log OR的先验（支持负值，弱信息先验，95%先验区间对应率差约-0.3~0.3）
theta1_prior_mean = 0
theta1_prior_sd = 0.7

# 按参数名指定先验，避免索引赋值错误
priors = c(
  set_prior(paste0("normal(", theta0_prior_mean, ",", theta0_prior_sd, ")"), 
            class = "Intercept"),
  set_prior(paste0("normal(", theta1_prior_mean, ",", theta1_prior_sd, ")"), 
            class = "b", coef = "groupint")
)

# 拟合模型
model = brm(
  family = binomial,
  event | trials(n) ~ group,
  data = data,
  prior = priors,
  iter = n.iter,
  warmup = n.warmup,
  chains = n.chains,
  refresh = 0,
  seed = 123 # 设随机种子保证结果可复现
)

# 提取后验并正确计算率差
posteriors = as_draws_df(model) %>%
  mutate(
    p_con = fun_invlogit(b_Intercept), # 对照组后验率
    p_int = fun_invlogit(b_Intercept + b_groupint), # 干预组后验率
    rd = p_int - p_con # 两组率差（干预-对照），支持负值
  )

# 计算率差的95%可信区间
quantile(posteriors$rd, c(0.025, 0.975))

# 绘制后验分布
posteriors %>%
  select(rd, p_con, p_int) %>%
  pivot_longer(cols = everything(), names_to = "parameter", values_to = "value") %>%
  ggplot(aes(x = value, fill = parameter)) +
  geom_density(alpha = 0.5) +
  theme_light()

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