假设按时间顺序进行了三次二项式实验。对于每个实验,我都知道#of试验以及#of'scescesses。为了使用前两个较旧的实验作为第三个实验的先验,我想“在两个较旧的实验上拟合贝叶斯分层模型,并使用与第三个实验相同的后验形式”。
鉴于我的可用数据(如下),我的问题是:我的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
相同。