>假设有三个按时间顺序进行的二项式实验。对于每个实验，我都知道#of 试验以及#of 成功。为了将前两个较旧的实验用于第三个实验，我想"在两个较旧的实验上拟合一个贝叶斯分层模型，并在第三个实验中使用先验形式"。
鉴于我的可用数据(如下)，我的问题是：我下面的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?
TL;DR：如果你直接预测成功的概率，那么模型将是带有参数theta(成功概率)的伯努利似然，其值可能介于0到1之间。在这种情况下，您可以在 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的先验，我们需要一个可以取任何实值的概率分布。如果我们对成功的几率一无所知，我们可能会使用一个非常弱的信息先验，比如正态分布，比如说，平均值=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相同。
回答(现已删除)评论：如果您只知道试验和成功的次数，并且您认为二项式概率是数据如何生成的合理模型，那么如何将数据拆分为"先验"和"可能性"或是否打乱数据的顺序并不重要。生成的模型拟合将相同。