用于二项式实验的贝叶斯分层建模的 rstanarm

Question

假设按时间顺序进行了三个二项式实验。对于每个实验，我知道 #of trials 以及 #of successes。要将前两个旧实验用作第三个实验的先验，我想 "fit a Bayesian hierarchical model on the two older experiments and use the posterior form that as prior for the third experiment".

鉴于我的可用数据（如下），我的问题是：我的 rstanarm 下面的代码是否捕获了我上面描述的内容？

Study1_trial = 70
Study1_succs = 27
#==================
Study2_trial = 84
Study2_succs = 31
#==================
Study3_trial = 100
Study3_succs = 55

我在包 rstanarm 中尝试过的内容：

library("rstanarm")

data <- data.frame(n = c(70, 84, 100), y = c(27, 31, 55));
mod <- stan_glm(cbind(y, n - y) ~ 1, prior = NULL, data = data, family = binomial(link = 'logit'))   

## can I use a beta(1.2, 1.2) as prior for the first experiment?

Answer 1

TL;DR: 如果您直接预测成功概率，则该模型将是伯努利似然法，其参数 theta（成功概率）可以取值在零和一之间。在这种情况下，您可以对 theta 使用 Beta 先验。但是使用逻辑回归模型，您实际上是在对成功的对数几率进行建模，它可以采用从 -Inf 到 Inf 的任何值，因此具有正态分布的先验（或其他可以采用任何实际值的先验）由可用的先验信息确定的某个范围）会更合适。

对于唯一参数是截距的模型，先验是对数成功几率的概率分布。在数学上，模型是：

log(p/(1-p)) =  a

其中 p 是成功的概率，a，您正在估计的参数，是截距，它可以是任何实数。如果成功的几率是 1:1（即 p = 0.5），那么 a = 0。如果赔率大于 1:1，则 a 为正。如果赔率小于 1:1，则 a 为负。

因为我们想要 a 的先验，所以我们需要一个可以取任何实数值的概率分布。如果我们对成功的几率一无所知，我们可能会使用信息量非常弱的先验信息，例如正态分布，例如 mean=0 和 sd=10（这是 rstanarm 默认值），这意味着一个标准偏差将包含从大约 22000:1 到 1:22000 的成功几率！所以这个先验基本上是平坦的。

如果我们采用您的前两项研究来构建先验，我们可以使用基于这些研究的概率密度，然后将其转换为对数几率：

# Possible outcomes (that is, the possible number of successes)
s = 0:(70+84)

# Probability density over all possible outcomes
dens = dbinom(s, 70+84, (27+31)/(70+84))

假设我们将对先验使用正态分布，我们需要最有可能的成功概率（这将是先验的均值）和均值的标准差。

# Prior parameters
pp = s[which.max(dens)]/(70+84)  # most likely probability
psd = sum(dens * (s/max(s) - pp)^2)^0.5  # standard deviation

# Convert prior to log odds scale
pp_logodds = log(pp/(1-pp))
psd_logodds = log(pp/(1-pp)) - log((pp-psd)/(1 - (pp-psd)))

c(pp_logodds, psd_logodds)

[1] -0.5039052  0.1702006

您可以通过运行 stan_glm 在前两项研究中使用默认（平坦）先验生成基本相同的先验：

prior = stan_glm(cbind(y, n-y) ~ 1, 
                 data = data[1:2,], 
                 family = binomial(link = 'logit'))   

c(coef(prior), se(prior))

[1] -0.5090579   0.1664091

现在让我们使用研究 3 中的数据使用默认先验和我们刚刚生成的先验来拟合模型。我已经切换到标准数据框，因为当数据框只有一行时 stan_glm 似乎失败（如 data = data[3, ]）。

# Default weakly informative prior
mod1 <- stan_glm(y ~ 1, 
                 data = data.frame(y=rep(0:1, c(45,55))), 
                 family = binomial(link = 'logit'))   

# Prior based on studies 1 & 2
mod2 <- stan_glm(y ~ 1, 
                 data = data.frame(y=rep(0:1, c(45,55))), 
                 prior_intercept = normal(location=pp_logodds, scale=psd_logodds), 
                 family = binomial(link = 'logit'))

为了进行比较，我们还生成一个包含所有三项研究和默认平坦先验的模型。我们希望这个模型给出与 mod2 几乎相同的结果：

mod3 <- stan_glm(cbind(y, n - y) ~ 1, 
                 data = data, 
                 family = binomial(link = 'logit'))

现在让我们比较一下这三个模型：

library(tidyverse)

list(`Study 3, Flat Prior`=mod1, 
     `Study 3, Prior from Studies 1 & 2`=mod2, 
     `All Studies, Flat Prior`=mod3) %>% 
  map_df(~data.frame(log_odds=coef(.x),
                     p_success=predict(.x, type="response")[1]), 
         .id="Model")

                              Model   log_odds p_success
1               Study 3, Flat Prior  0.2008133 0.5500353
2 Study 3, Prior from Studies 1 & 2 -0.2115362 0.4473123
3           All Studies, Flat Prior -0.2206890 0.4450506

对于具有平坦先验（第 1 行）的研究 3，预测的成功概率为 0.55，正如预期的那样，因为数据就是这么说的，而先验没有提供额外的信息。

对于基于研究 1 和 2 的先验研究 3，成功的概率为 0.45。较低的成功概率是由于研究 1 和 2 添加额外信息的成功概率较低。事实上，mod2 的成功概率正是您直接根据数据计算得出的值：with(data, sum(y)/sum(n))。 mod3 将所有信息放入可能性中，而不是将其分为先验和可能性，但在其他方面与 mod2.

基本相同

对（现已删除）评论的回答： 如果您所知道的只是试验次数和成功次数，并且您认为二项式概率是一个合理的数据模型生成，那么您如何将数据拆分为 "prior" 和 "likelihood" 或是否打乱数据的顺序都无关紧要。生成的模型拟合将是相同的。

用于二项式实验的贝叶斯分层建模的 rstanarm

rstanarm for Bayesian hierarchical modeling of binomial experiments

r

bayesian

logistic-regression

rstan

rstanarm