用于二项分布实验贝叶斯层次建模的rstanarm方法

Question

假设有三个二项式实验按时间顺序进行。对于每一个实验，我都知道试验的#以及成功的#。为了在第三个实验中使用前两个老实验，我想“在两个老实验上拟合一个贝叶斯层次模型，并使用后验形式作为第三个实验的先验形式”。

考虑到我的可用数据(下面)，我的问题是:下面的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:，如果你直接预测成功的概率，这个模型将是一个贝努利似然，参数θ(成功的概率)，可以取0到1之间的值。在这种情况下，你可以用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

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

对于先前基于研究1和2的研究3，成功的概率为0.45。成功的可能性较低是因为在研究1和2中添加了更多的信息，成功的概率较低。事实上，来自mod2的成功概率正是您从数据中直接计算出来的：with(data, sum(y)/sum(n))。mod3将所有信息放入可能性中，而不是在先验和可能性之间进行分割，但在其他方面与mod2本质上是相同的。

回答(现在删除了)评论：，如果你只知道试验和成功的数量，并且你认为二项式概率是数据生成的合理模型，那么如何将数据划分为“先验”和“可能性”，或者是否调整数据的顺序，都无关紧要。得到的模型拟合将是相同的。