Stata与R中贝叶斯逻辑回归结果不一致的排查求助

问题描述 #

我在Stata中运行了针对两个二分变量（exposure暴露因素、outcome结局）的贝叶斯逻辑回归，代码如下：

bayes, nchains(5) ///
        prior({excellent:tx}, normal(0,2))  ///
        rseed(123456) : ///
        logistic outcome exposure   

使用同一数据集，尝试用R的brms包复现分析，代码如下：

priors <- c(
  set_prior("normal(0, 2)", class = "b", coef = "tx"),
  set_prior("normal(0, 10000)", class = "Intercept")
)
bayes_model <- brm(
  formula = outcome ~ exposure,
  family = bernoulli(link = "logit"),
  data = data,
  prior = priors,
  chains = 5,
  iter = 12500,  
  warmup = 2500,
  seed = 123456
)
summary(bayes_model)
plot(bayes_model)

# 输出指数化的比值比
fixed_effects <- fixef(bayes_model)
odds_ratios <- exp(fixed_effects)
odds_ratios

但两者输出的比值比（OR）估计值差异极大，已确认使用同一数据集并尝试匹配Stata默认设置，需排查原因。

核心排查要点 #

• 先验变量名完全不匹配：Stata代码里先验指定的是{excellent:tx}，但R代码中模型公式的自变量是exposure，却给coef = "tx"设置先验——这会导致R中exposure参数实际使用brms的默认先验，而非你设置的normal(0,2)，和Stata的先验配置完全不一致，这是最可能的原因。修正方式是把R里的先验改为coef = "exposure"。
• 迭代与预热设置不匹配：Stata的bayes命令默认每个链迭代2000次，其中500次为预热；你R中设置的是iter=12500、warmup=2500，需确认是否和Stata的实际运行设置一致。若Stata未指定迭代参数，两者的样本量和收敛过程会有差异。
• 截距先验不一致：需确认Stata中bayes logistic的默认截距先验是什么，你R中设置的normal(0,10000)是否和Stata默认一致。比如Stata默认可能采用平坦先验，这会影响截距估计，进而影响整体结果。
• 变量编码反转：检查两个软件中exposure和outcome的0/1编码是否对应同一类别，若出现编码反转（如Stata中暴露为1，R中为0），会导致OR值变为原结果的倒数，出现极大差异。
• 收敛状态差异：分别检查模型收敛情况：R中查看summary(bayes_model)里的rhat值（需接近1）；Stata中运行bayesstats diagnostics。若任一模型未收敛，结果无可靠性，也会出现差异。

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