发病率的分层贝叶斯模型

机器算法验证 贝叶斯 分层贝叶斯 pymc 概率规划
2022-03-12 14:22:09

Kevin Murphy 的讨论了一个经典的分层贝叶斯问题(最初在 中讨论Johnson and Albert, 1999, p24):

假设我们试图估计癌症发病率N城市。在每个城市,我们抽样了一些人Ni并测量患有癌症的人数xiBin(Ni,θi), 在哪里θi是城市的真实癌症发病率。

我们想估计θi同时允许数据贫乏的城市从数据丰富的城市借用统计实力。

为此,他建模θiBeta(a,b)使所有城市共享相同的先验,因此最终模型如下所示:

p(D,θ,η|N)=p(η)i=1NBin(xi|Ni,θi)Beta(θi|η)

在哪里η=(a,b).

这个模型的关键部分当然是(我引用),“我们推断η=(a,b)从数据中,因为如果我们只是将其固定为一个常数,则θi将是有条件独立的,它们之间不会有信息流动”


我正在尝试在PyMC中对此进行建模,但据我所知,我需要先验ab(我相信这是p(η)多于)。

这个模型有什么好的先验?

如果它有帮助,我现在拥有的代码是:

bins = dict()
ps   = dict()
for i in range(N_cities):
    ps[i]   = pm.Beta("p_{}".format(i), alpha=a, beta=b)
    bins[i] = pm.Binomial('bin_{}'.format(i), p=ps[i],n=N_trials[i],  value=N_yes[i], observed=True)

mcmc = pm.MCMC([bins, ps])
1个回答

在 Gelman,贝叶斯数据分析,(第 2 版,第 128 页;第 3 版,第 110 页)中讨论了类似的问题。格尔曼建议先p(a,b)(a+b)5/2,这有效地限制了“先验样本量”a+b,因此 beta 超先验本身不太可能提供大量信息。(作为数量a+b增长,β分布的方差缩小;在这种情况下,较小的先验方差会限制后验数据的“权重”。)此外,此先验不设置是否a>b,或相反,所以适当的分布对(a,b)是从所有数据中推断出来的,正如您在这个问题中所希望的那样。

Gelman 还建议根据均值的 logit 重新参数化模型θ和先验的“样本量”。所以而不是做推断(a,b)直接,问题是关于转换量的推断logit(aa+b)log(a+b). 这允许在实平面上转换的先验值,而不是必须严格为正的未转换的先验值。此外,这实现了绘制时更加分散的后验密度。这使随附的图表更加清晰,我觉得这很有帮助。