Single-parameter models 单参数模型
模型中只含有一个待估参数,也就是说,估计值 $\theta$ 是一维的,包括二项分布、正态分布、泊松分布和指数分布四类经典的模型。同时介绍了贝叶斯数据分析的重要概念和计算方法。
在简单的二项式模型中,目标是从一系列“伯努利试验”的结果中估计一个未知的人口比例,即数据$y_1,...,y_n$,每个数据都是$\theta$或 1。这个问题为贝叶斯推断的讨论提供了一个相对简单但重要的起点。
二项分布提供了一个自然的数据模型,这些数据来自一系列 $n$ 个可交换的试验,或者来自一个大群体,每个试验产生两种可能的结果之一,通常标记为“成功”和“败”由于可交换性,数据可以总结为 $n$ 个试验的成功总数,我们在这里用 $y$ 表示。通过让参数 $\theta$ 表示总体成功的比例,或者等价地表示每个试验成功的概率,可以很自然地从可交换试验的公式转换为使用独立同分布的随机变量的公式。二项式抽样模型是,
$$
P(y|θ)=Bin(y|n,θ)=\left(\begin{matrix}
n\\
y
\end{matrix}\right)\theta^y(1-\theta)^{n-y}
$$
在左边,我们抑制了对$n$的依赖性,因为它被认为是实验设计的一部分,被认为是固定的; 这个问题讨论的所有概率都假定是以 $n$ 为条件的。
作为二项式模型的一个具体应用,我们考虑人口出生性别比的估计。女性出生的比例长期以来一直是科学和普通公众感兴趣的话题。200 年前,欧洲人口中女性出生的比例低于0.5,而在本世纪,人们的兴趣集中在可能影响性别比例的因素上。目前欧洲人口中女性出生比例的公认值为 0.485。
对于这个例子,我们将参数 $\theta$ 定义为女性出生率的比例,但是解释这个参数的另一种方法是将其作为男女出生率的比例,$\phi=(1-\theta)/\theta$。 设$y$是$n$个出生记录中的女孩数。应用二项模型:
==图 1==: 二项参数 $\theta$ 的未归一化后验密度,基于均匀的先验分布和 $n$ 次试验的成功,显示 $n$ 和 $y$ 的几个值的曲线。
(1) 中,我们假设 $n$ 个出生是条件独立的,给定 $\theta$ ,在所有情况下,女性出生的概率等于 $\theta$ 。当我们没有可能影响婴儿性别的解释性信息(例如,区分同一家庭中的多胞胎或多胞胎)时,可以判断出可交换性,从而推动了这种建模假设。
为了在二项式模型中执行贝叶斯推断,必须指定 $\theta$ 的先验分布。我们将多次讨论与指定先验分布相关的问题,但为了简单起见,在这一点上,我们假设 $\theta$ 的先验分布在区间[0,1]上是一致的。
初步应用贝叶斯规则,给出 $\theta$ 为
$$
p(\theta|y)∝\theta^y(1-\theta)^{n-y}
$$
在 $n$ 和 $y$ 不变的情况下,因子$\left(\begin{matrix}n\y\end{matrix}\right)$不依赖于未知参数 $\theta$ ,因此在计算 $\theta$ 的后验概率分布时,可以将其视为一个常数。作为许多例子的典型,后验密度可以立即以封闭形式写入,直到一个比例常数。
在单参数问题中,这允许立即以图形方式表示后验概率分布。
例如,在图1 中,显示了几个不同实验的未归一化密度,即 $n$ 和 $y$ 的不同值。在本例中,我们可以认为$\beta$ 分布的非标准化形式为 ,
$$
\theta|y\sim Beta(y+1,n-y+1)
$$
贝叶斯推断的过程包括从先验概率 $p(\theta)$ 到后验概率 $p(\theta|y)$ 的过程,人们很自然地期望这两个分布之间可能存在一些一般关系。例如,我们可以预期,由于后验概率分布包含了来自数据的信息,因此它比前验概率分布的可变性更小。这个概念在下面的第二个表达式中被形式化了:
$$
E(\theta)=E(E(\theta|y)) ①
$$
还有
$$
var(\theta)=E(var(\theta|y))+var(E(\theta|y)) ②
$$
方程①所表达的结果并不令人惊讶: $\theta$ 的先验平均值是可能数据分布上所有可能的后验平均值的平均值。方差公式②更有趣,因为它说后验方差平均小于前验方差,其数量取决于后验平均值在可能数据分布上的变化。后者的变化越大,就越有可能减少 $\theta$ 的不确定性。均值和方差关系只描述预期,在特定情况下,后验方差可能类似于或甚至大于先验方差(尽管这可能表明抽样模型和先验分布之间的冲突或不一致)。
在具有均匀先验分布的二项式例子中,先验均值为$\frac{1}{2}$,先验方差为$\frac{1}{12}$。后验平均值$\frac{y+1}{n+2}$是先验平均值与样本比例$\frac{y}{n}$之间的折衷,在后验平均值中,随着样本数量的增加,先验平均值的作用越来越小。这是贝叶斯推断的一个普遍特征: 后验概率分布集中在代表先验信息和数据之间的折衷点上,随着样本量的增加,折衷在很大程度上受到数据的控制。
后验概率分布包含了当前关于参数 $\theta$ 的所有信息。理想情况下,可以报告整个后验概率$p(\theta|y)$; 如我们所见
图 2假设 95% 的中心间隔和 95% 的最高后部密度区域有显著差异: (a)中心后部间隔,(b)最高后部密度区域。
通过模拟实现的贝叶斯方法的一个关键优势是,即使在复杂的转换之后,后验推理也可以被总结的灵活性。然而,出于许多实际目的,对分布的各种数值概括是可取的。常用的位置概括是分布的平均值、中位数和模式; 变异通常用标准差、四分差范围和其他分位数来概括。每个总结都有自己的解释: 例如,平均值是参数的后验期望值,模式可以解释为给定数据(和模型)的单个“最有可能”值。此外,正如我们将看到的,许多实际推理依赖于正态近似的使用,通常通过对 $\theta$ 应用对称变换来改进,在这里,均值和标准差起着关键作用。这种模式在更复杂问题的计算策略中很重要,因为它通常比平均值或中值更容易计算。
当后验概率分布具有封闭形式时,如本例中的 $\beta$ 分布,则后验概率分布的平均值、中值和标准差通常以封闭形式提供。
除了点总结之外,报告后验不确定性几乎总是重要的。我们通常的方法是给出感兴趣的估计值的后验分布的分位数,或者,如果需要一个区间汇总,一个后验概率的中心区间,在 $100(1-\alpha)%$ 的区间中,这个区间对应于正好位于后验概率的 $100(\alpha/2)%$ 以上和以下的值的范围。这样的区间估计被称为后验区间。对于简单的模型,如二项式和正态分布,后验区间可以直接从累积分布函数计算出来,通常使用对标准计算机函数的调用。一般来说,间隔可以通过计算机模拟从后验概率分布计算出来。
对后验不确定性的一个略有不同的概括是最高的后验概率区域: 包含后验概率 $100(1-\alpha)%$ 的一组值,并具有区域内的密度从不低于区域外的密度的特点。如果后验概率分布是单峰和对称的,则这样的区域与中心后验概率相同。在目前的实践中,中央后部间隔是常用的,部分原因是它有一个直接的解释作为后部 $\alpha/2$ 和 $1-\alpha/2$ 分位数,部分原因是它是直接计算使用后部模拟。
图 2显示不同的后验概括看起来差别很大的一种情况: 95% 的中心区域包括分布中心的零概率区域,而 95% 的最高后验密度区域包括两个不相交的区域。在这种情况下,最高后验密度区域比中心间隔更麻烦,但传递更多的信息; 然而,最好不要试图用任何单一间隔来总结这种双峰密度。当后验密度高度偏斜时,中心间隔和最高后验密度区域也可能大不相同。
正态分布是大多数统计建模的基础。中心极限定理有助于证明在许多统计问题中使用正态似然,作为解析上不太方便的实际似然的近似。此外,正如我们将在后面的章节中看到的,即使正态分布本身并不能提供一个很好的模型拟合,它也可以作为涉及 $t$ 或有限混合分布的更复杂模型的一个组成部分。现在,我们只是在假设正态模型是合适的情况下通过贝叶斯结果进行工作。我们首先推导出单个数据点的结果,然后推导出多个数据点的一般情况。
作为最简单的第一种情况,考虑一个由平均 $\theta$ 和方差 $\sigma_2$ 参数化的正态分布的单个标量观测 $y$ ,对于这个初始发展,我们假设 $\sigma_2$ 是已知的。抽样分布是
$$
p(y|\theta)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{1}{2\sigma^2}(y-\theta)^2}
$$
作为 $\theta$ 的函数,可能性是 $\theta$ 中二次型的指数,因此共轭先验密度族看起来像
$$
p(\theta)=e^{A\theta^2+B\theta+C}
$$
我们将这个族参数化为
$$
p(\theta)\propto exp\left(-\frac{1}{2\tau^2_0}(\theta-\mu_o)^2\right)
$$
即 $\theta\sim n(\mu_0,\tau^2)$ ,超参数 $\mu_0$ 和 $\tau_0^2$ 。像通常一样,在这个初步的发展,我们假设超参数是已知的。
共轭先验密度意味着 $\theta$ 的后验概率分布是二次型的指数,因此是正态分布,但是需要一些代数来揭示它的特殊形式。在后验密度中,除了 $\theta$ 以外的所有变量都被视为常数,给出条件密度,
$$
p(\theta|y)\propto exp\left(-\frac{1}{2}\left(\frac{(y-\theta)^2}{\sigma^2}+\frac{(\theta-\mu_0)^2}{\tau_0^2}\right)\right)
$$
扩展指数,合并项,然后在 $\theta$ 中完成平方
$$
p(\theta|y)\propto exp\left(-\frac{1}{2\tau_1^2}(\theta-\mu_1)^2\right)
$$
即$\theta|y\sim N(\mu_1,\tau_1^2)$ ,其中
$$
\mu_1=\frac{\frac{1}{\tau_0^2}\mu_0+\frac{1}{\sigma^2}y}{\frac{1}{\tau_0^2}+\frac{1}{\sigma^2}} 和 \frac{1}{\tau_1^2}=\frac{1}{\tau_0^2}+\frac{1}{\sigma^2}
$$
在操纵正态分布时,方差的逆起着重要作用,称为精度。上面的代数表明,对于正态数据和正态先验分布(每个都具有已知精度) ,后验精度等于先验精度加上数据精度。
有几种不同的方法来解释后验平均值 $\mu_1$ 的形式。后验平均数表示为先验平均数和观测值 $y$ 的加权平均数,加权平均数与精度成正比。或者,我们可以将$\mu_1$表示为向观测 $y$ 调整的先验平均值,
$$
\mu_1=\mu_0+(y-\mu_0)\frac{\tau_0^2}{\sigma^2+\tau_0^2}
$$
或者当数据“收缩”到先前的平均值时,
$$
\mu_1=y-(y-\mu_0)\frac{\sigma^2}{\sigma^2+\tau_0^2}
$$
每个公式都代表后验平均值,作为先验平均值和观测值之间的折衷。
在极端情况下,后验平均值等于先前平均值或观察到的数据:
$$
\mu_1=\mu_0 若 y=\mu_0 或 \tau_0^2=0
$$
$$
\mu_1=y 若 y=\mu_0 或 \sigma^2=0
$$
如果$\tau_0^2=0$,则先验分布比数据精确无限,因此后验分布和先验分布相同且集中在$\mu_0$处。如果 $\sigma^2=0$ ,数据是完全精确的,并且后验概率分布集中在观察值 $y$ 上。如果 $y=\mu_0$ ,则先验平均值和数据平均值一致,后验平均值也必须在此处下降。
标准分布——二项分布、正态分布、泊松分布和指数分布——都是从简单概率模型自然导出的。
二项分布是由计算可交换结果而产生的,正态分布适用于一个随机变量,这个随机变量是许多可交换或独立项的和。我们还将有机会将正态分布应用于所有正数据的对数,这将自然地适用于被建模为许多独立乘法因子的乘积的观测。
泊松分布和指数分布随着事件的计数次数和等待时间的增加而出现,事件模型在所有时间间隔内都是可交换的,也就是说,在时间上是独立的,出现率是恒定的。我们通常会通过这些基本分布的组合来构建更复杂结果的现实概率模型。
这些标准模型中的每一个都有一个相关的共轭先验分布家族,我们依次讨论。
具有已知均值 $\theta$ 和未知方差的正态模型是一个重要的例子,不一定是因为它的直接应用价值,而是作为更复杂,有用的模型。此外,具有已知均值但未知方差的正态分布提供了尺度参数估计的介绍性示例。
对于已知 $\theta$ 且$\sigma^2$未知的$p(y|\theta,\sigma^2)=N(y|\theta,\sigma^2)$,$n$ 个独立同分布观测值的向量$y$ 的可能性为
$$
p(y|\sigma^2)\propto\sigma^{-n}exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\theta)^2\right)=(\sigma^2)^{-n/2}exp\left(-\frac{n}{2\sigma^2}\upsilon\right)
$$
充分的统计数据是
$$
\upsilon=\frac{1}{n}\sum^n_{i=1}(y_i-\theta)^2
$$
对应的共轭先验密度是反伽马,
$$
p(\sigma^2)\propto(\sigma^2)^{-(\alpha+1)}e^{-\beta/\sigma^2}
$$
它有超参数$(\alpha,\beta)$。一种简便的参量化方法是将$\sigma^2$的先验分布假定为反$\chi^2$的分布,刻度为$\sigma_0^2$和$\nu_0$自由度(见附录 A) 其中$\chi$是一个$\chi_{\nu_0}^2$ 的随机变量。我们使用方便但非标准的符号,$\sigma^2\sim Inv-\chi^2(\nu_0,\sigma_0^2)$。
$\sigma^2$ 的后验密度是
$$
p(\sigma^2|y)\propto p(\sigma^2)p(y|\sigma^2)
$$
$$
\propto \left(\frac{\sigma_0^2}{\sigma^2}\right)^{\nu/2+1}exp\left(-\frac{\nu_0\sigma_0^2}{2\sigma^2}\right)\cdot(\sigma^2)^{-n/2}exp\left(-\frac{n}{2}\frac{\nu}{\sigma^2}\right)
$$
$$
\propto(\sigma^2)^{-((n+\nu_0)/2+1)}exp\left(-\frac{1}{2\sigma^2}(\nu_0\sigma_0^2+n\nu)\right)
$$
是一个刻度反$\chi^2$分布,刻度等于先验的自由度加权平均值,数据刻度和自由度等于先验和数据自由度之和。先验分布可以被认为是提供了相当于$\nu_0$观测值的平均平方差 $\sigma_0^2$ 的信息。
泊松分布自然地出现在以计数形式研究数据的过程中; 例如,一个主要的应用领域是研究疾病发病率的流行病学。
如果一个数据点 $y$ 遵循泊松分布,其速率为 $\theta$ ,那么一个观测点 $y$ 的概率分布就是概率分布
$$
p(y|\theta)=\frac{\theta^ye^{-\theta}}{y!},对于y=0,1,2,...
$$
对于独立同分布观测的向量 $y= (y_1,. . y_n)$ ,
可能性是
$$
p(y|\theta)=\prod^n_{i=1}\frac{1}{y_i!}\theta^{y_i}e^{-\theta}\propto \theta^{t(y)}e^{-n\theta}
$$
其中 $t(y)=\sum^n_{i=1}y_i$ 是充分的统计量。我们可以将指数族形式的可能性改写为
$$
p(y|\theta)\propto e^{-n\theta}e^{t(y)log\theta}
$$
显示自然参数是$\phi(\theta)=log\theta$,自然共轭先验分布是
$$
p(\theta)\propto (e^{-\theta})^{\eta}e^{\nu log\theta}
$$
是超参数指标$(\eta,\nu)$。换句话说,这个论点的可能性是$\theta^\alpha e^{-b\theta}$,因此共轭先验密度必须是$p(\theta)\propto\theta^Ae^{-B\theta}$。在一个更常规的参量化中,
$$
p(\theta)\propto e^{-\beta\theta}\theta^{\alpha-1}
$$
这是一个伽马密度,参数为 $\alpha$ 和 $\beta$ ,$Gamma(\alpha,\beta)$。通过比较$p($y$ |\theta)$和$p(\theta)$,发现先验密度在某种意义上等同于 $\beta$ 先验观测中 $\alpha-1$ 的总计数。用这种共轭先验分布,后验概率分布是
$$
\theta|y\sim Gamma(\alpha+n\overline{y},\beta+n)
$$
对于共轭族,先验密度和后验密度的已知形式可以用公式求出边缘分布
$$
p(y)=\frac{p(y|\theta)p(\theta)}{p(\theta|y)}
$$
例如,对于单个观测点的泊松模型,$y$,具有先验预测分布泊松
$$
p(y)=\frac{Poisson(y|\theta)Gamma(\theta|\alpha,\beta)}{Gamma(\theta|\alpha+y,1+\beta)}=\frac{\Gamma(\alpha+y)\beta^\alpha}{\Gamma(\alpha)y!(1+\beta)^{\alpha+y}}
$$
减少到
$$
p(y)=\left(\begin{matrix}{\alpha+y-1}\y\end{matrix}\right)\left(\frac{\beta}{\beta+1}\right)^\alpha\left(\frac{1}{\beta+1}\right)^y
$$
即负二项密度:
$$
y\sim Neg-bin(\alpha,\beta)
$$
上述推导表明,负二项分布是泊松分布和 $\theta$ 的结合, $\theta$ 遵循伽马分布:
$$
Neg-bin(y|\alpha,\beta)=\int Poisson(y|\theta)Gamma(\theta|\alpha,\beta)d\theta
$$
在许多应用中,将数据点$y_1,...,y_n$的泊松模型扩展到形式是很方便的
$$
y_i\sim Poisson(x_i\theta)
$$
其中值 $x_i$ 是解释变量的已知正值,$x$ 和 $\theta$ 是感兴趣的未知参数。在流行病学中,参数 $\theta$ 通常被称为速率,而 $x_i$ 被称为第 $i$ 单位的暴露。这个模型在 $y_i$ 中不可交换,但在 $(x,y)i$对中可交换。扩展泊松模型中 $\theta$ 的可能性是
$$
p(y|\theta)\propto \theta(\sum {i=1}^ny_i)e^{-\sum^n_{i=1}x_i\theta}
$$
(忽略不依赖于 $\theta$ 的因子) ,因此 $\theta$ 的伽马分布是共轭的。先验分布
$$
\theta\sim Gamma(\alpha,\beta)
$$
得到的后验概率分布是
$$
\theta|y\sim Gamma\left(\alpha+\sum_{i=1}^ny_i,\beta+\sum_{i+1}^nx_i\right)
$$
假设美国某个城市的死亡原因被详细审查了一年。据发现,在 20 万人口中有 3 人死于哮喘,粗略估计该市每年每 10 万人中有 1.5 人死于哮喘。泊松试验模型经常用于这种形式的流行病学数据。泊松模型来源于所有小的暴露间隔之间的可交换性假设。在泊松模型下, $y$ 的抽样分布,即一个城市一年内 200,000 人的死亡人数,可以用$Poisson(2.0 \theta)$表示,其中 $\theta$ 代表我们城市真实的长期哮喘死亡率(以每 100000 人每年的病例计算)。
在上面的表示法中,$y=3$是一个单一的观察结果,暴露$x=2.0$(因为 $\theta$ 是以 100000 人为单位定义的) ,未知的速率 $\theta$ 。我们可以利用世界各地关于哮喘死亡率的知识来构造 $\theta$ 的先验概率分布,然后将 $y=3$ 和先验概率分布结合起来得到后验概率分布。
$\theta$ 的合理先验分布是什么?
对世界各地哮喘死亡率的回顾表明,西方国家的死亡率高于每 10 万人 1.5 人的情况很少见,典型的哮喘死亡率约为每 10 万 人 0.6 人。对这个问题的共轭先验家族——伽马分布的性质进行的反复探索表明,如果我们假设这个城市和其他城市之间以及今年和其他年份之间存在可交换性,那么在这个例子中,$Gamma(3.0,5.0)$密度为哮喘死亡率提供了一个合理的先验密度。先验分布的平均值为 0.6(模式为 0.4) ,97.5% 的密度质量在 1.44 以下。在实际应用中,指定先验均值可以设定两个伽马参数的比值,然后通过反复试验改变形状参数,使之与分布尾部的先验知识相匹配。后验概率。结果表明,$Gamma(\alpha,\beta)$先验分布 $\theta$ 的后验概率为 $Gamma(\alpha+$y$,\beta+x)$。根据先验分布和数据描述,$\theta$ 的后验概率为$Gamma(6.0,7.0)$,平均后验概率为 0.86ー向先验分布发生了实质性收缩。1000 的直方图来自后验概率分布。例如,根据 伽马后验概率计算,我市哮喘长期死亡率大于每 10 万人每年 1.0 人的后验概率0.30。
图 3 $\theta$ 的后验密度,每 10 万人每年的哮喘死亡率, 以$Gamma(3.0,5.0)$为先验分布: (a)给定$y=3$万人中的 3 人死亡; (b)给定 $y=30$ 万恒定人口 10 年中的 30 人死亡。直方图显示参差不齐,因为
它们仅从每种情况下的后验概率分布中的 1000 个随机抽取构建。
附加数据的后验概率分布。
为了考虑额外数据的影响,假设在我们的例子中获得了 10 年的城市数据,而不仅仅是 1 年,我们发现死亡率维持在每 10 万人中 1.5 人,我们发现 10年内 $y=30$ 人死亡。假设人口在 200,000 不变,并且假设 10 年的结果与恒定的长期比率是独立的$\theta$,$\theta$ 的后验概率是 $Gamma(33.0,25.0)$ ; 图 3b显示了从这个分布中得出的 1000 个结果。后验概率比以前集中得多,它仍然位于先验分布和数据之间。经过 10 年的数据,$\theta$ 的后验概率为 1.32,$\theta$ 超过 1.0 的后验概率为 0.93。
指数分布通常用于模拟“等待时间”和其他连续的、正的、实值的随机变量,通常是在一个时间尺度上测量的。给定参数 $\theta$ 的结果 $y$ 的抽样分布是
$$
p(y|\theta)=\theta exp(-y\theta),对于y>0,
$$
$\theta=1/E(y|\theta)$ 称为“速率”。在数学上,指数是伽马分布的一种特殊情况,其参数$(\alpha,\beta)=(1,\theta)$。然而,在这种情况下,它被用作结果 $y$ 的抽样分布,而不是参数的先验分布$\theta$,就像泊松的例子。
指数分布具有“无记忆”的特性,使其成为生存或生命周期数据的自然模型; 一个物体在一个额外时间 $t$ 中存活的概率与到达这一点的时间无关: 对于任何 $s$ , $t$ ,$Pr(y>t+s|y>s,\theta)=Pr(y>t|\theta)$。指数参数 $\theta$ 与泊松均值的共轭先验概率分布为$Gamma(\theta|\alpha,\beta)$与相应的后验概率分布$Gamma(\theta|\alpha+1,\beta+y)$。$n$ 个独立指数观测值的抽样分布,$y=(y_1,...,y_n)$,恒定速率 $\theta$ 是
$$
p(y|\theta)=\theta^nexp(-n\overline{y}\theta),对于\overline{y}\geq0,
$$
当看作 $\theta$ 的可能性时,对于固定的 $y$ ,与 $Gamma(n+1,n\overline{y})$ 密度成正比。因此,$\theta$ 的$Gamma(\alpha,\beta)$先验分布可以看作是总等待时间 $\beta$ 的 $\alpha-1$ 指数观测值。
图 4 1980-1989 年美国白人男性肾癌/输尿管癌年龄标准化死亡率最高的 10% 的县。
上一节考虑了在给定数据量的情况下,先验分布对推理的影响。相比之下,在这里,每个推论基于不同的数据,但具有共同的先验分布。除了说明先验分布的作用之外,这个例子还引入了层次建模。
图4显示了 20 世纪 80 年代美国肾癌死亡率最高的县。地图中最显著的模式是,位于美国中部大平原的许多县被遮盖了,但是靠近海岸的县相对较少。
当展示这张地图时,人们提出了许多理论来解释大平原上不适当的阴影: 也许空气或水被污染了,或者人们倾向于不去寻求医疗,因此癌症被发现得太晚而无法治疗,或者也许他们的饮食不健康。这些猜测可能都是正确的,但是它们实际上并不需要解释图4中的模式。要了解这一点,请看图5,它绘制了肾癌死亡率最低的 10% 的县。这些也大部分位于该国的中部。所以现在我们需要解释为什么这些地区的死亡率最低,也是最高的。
问题是样本的大小。考虑一个人口 1000 的县。肾癌是一种罕见的疾病,在任何一个10 年期间,一个 1000 人的县可能没有任何肾癌死亡人数,因此,它将与全国最低的肾癌死亡率并列,并将在图 5 中被遮盖。然而,这个县也有可能在十年内死于一例肾癌。如果是这样的话,那么它的比率将是每 10,000 人中有 1 人每年,这个比率足够高,可以把它放在前 10% 的位置,因此它将在图 4中被阴影处理。大平原有许多人口较少的县,因此在两张地图上都显示过多。从这些地图上没有证据表明那里的癌症发病率特别高。
图5:1980-1989 年美国白人男性肾癌/输尿管癌年龄标准化死亡率最低的 10% 的县。令人惊讶的是,该模式与图 4所示的最高发病率地图有些相似。
贝叶斯推断癌症死亡率
原始利率地图中的误导模式表明,基于模型的方法估计真实的基础利率可能是有帮助的。特别是,使用该模型估计每个县的潜在癌症死亡率是很自然的
$$
y_j\sim Poisson(10n_j\theta_j)
$$
其中$y_j$是 1980-1989 年 j 县的肾癌死亡人数,$n_j$是该县的人口,$\theta_j$是每年人均死亡单位的潜在比率。在这个符号中,图4和5中的地图绘制了原始死亡率$\frac{y_j}{10n_j}$(在这里,我们忽略了年龄标准化,尽管对模型进行概括以考虑到这一点是可能的。)
这个模型不同于 $\theta_j$ 在不同的县之间是不同的,所以上式是美国每个县的独立模型。我们使用上式中的下标 $j$ (而不是 $i$ )来强调这些是独立的参数,每个参数都是从自己的数据中估计出来的。如果我们只对其中一个县进行推断,我们只需要写 $y\sim Poisson(10n\theta)$ 。
为了执行贝叶斯推断,我们需要一个未知比率 $\theta_j$ 的先验分布。为了方便起见,我们使用了与泊松共轭的伽马分布。正如我们将在后面讨论的,参数 $\alpha= 20$ 和 $\beta= 430,000$ 的伽马分布是在此期间美国各县潜在肾癌死亡率的合理先验分布。这个先验分布的平均值为$\frac{\alpha}{\beta}=4.65\times10^{-5}$、标准差 为 $\frac{\sqrt \alpha}{\beta}=1.04\times10^{-5}\theta_j$ 的后验分布为,
$$
\theta|y_j \sim Gamma(20+y_j,430000+10n_j)
$$
其均值与方差为,
$$
E(\theta_j|y_j)=\frac{20+y_i}{430000+10n_j}
$$
$$
var(\theta_j|y_j)=\frac{20+y_i}{(430000+10n_j)^2}
$$
后验平均值可以看作是原始比率, $\frac{y_i}{10n_j}$ ,先验平均值, $\frac{\alpha}{\beta}=4.65\times10^{-5}$ 的加权平均数。
本地数据和先验分布的相对重要性
推断为一个小县。 先前信息和数据的相对权重取决于人口规模 $n_j$ 。例如,考虑一个 $n_j=1000$ 的小县:
对这个县来说,如果 $y_j=0$ ,那么原始死亡率是$\theta$ ,但是后面的平均值是 $\frac{20}{440000}=4.55\times10^{-5}$
如果 $y_i=1$ ,那么原始死亡率是每 10 年每 1000 人中有 1 人死亡,或者每人年 $10^{-4}$ 人死亡(大约是全国平均值的两倍) ,但后面的平均值仅为 $\frac{21}{440000}=4.77\times10^{-5}$
如果 $y_j=2$ ,那么原始死亡率是极高的 $2\times10^{-4}$ 每人-年,但是后验平均值仍然只是 $\frac{22}{440000}=5.00\times10{-5}$
在如此小的人口规模下,数据被先验分布所支配。
但是有多大的可能性?对于这个 $n_j=1000$ 的县来说,先验的 $y_j$ 会等于 0,1,2 等等?这是由预测分布决定的, $y_j$ 的边缘分布,平均于 $\theta_j$ 的先验分布。
如上一节所讨论的,具有伽马先验分布的泊松模型具有负二项预测分布:
$$
y_j\sim Neg-bin(\alpha,\frac{\beta}{10n_j})
$$
直接模拟 yj 的预测分布或许更简单,如下所示: (1)从 Gamma (20,430,000)分布中得出 $\theta_j$ 的 500 个值; (2)从参数为 $10000\theta_j$ 的泊松分布中得出每个值 $y_j$ 。用这种方法生产的 500 个 $y_j$ 模拟物中,319 个为 0,141 个为 1,33 个为 2,5个为 3。
*对一个大县的推论。
*现在考虑一个 $n_j=1000000$ 万的大县。在 10 年的时间里,我们期望看到多少癌症死亡?同样,我们可以使用 $Gamma(20,430000)$ 和 $Poisson (10^7 \theta_j)$ 分布来模拟预测分布的 500 个值 $y_j$ 。这样做,我们发现中位数为 473,间隔为 50% [393,545]。因此,这样一个县的原始死亡率有可能或不可能下降到 $3.94\times10^{-5}$ 和 $5.45\times10^{-5}$ 之间。
贝叶斯估计或“贝叶斯调整”死亡率怎么样?例如,如果 $y_i$ 取低值 393,则原始死亡率为 $3.94\times10^{-5}$ , $\theta_j$ 的后验平均值为 $\frac{20+393}{10^7+430000}=3.96\times10^{-5}$ ,如果 $y_i=545$ ,则原始死亡率为 $5.45\times10^{-5}$ ,后验均值为 $5.41\times10^{-5}$ 。在这个大县,数据较于先验分布占优势。
*比较不同大小的县。
*在泊松模型中, $\frac{y_i}{10n_j}$ 的方差与暴露参数 $n_j$ 成反比,因此可以认为是该县的“样本量”。图 6 显示了原始肾癌死亡率随人口而变化的情况。极高和极低的比率都在低人口县。相比之下,图 7a 显示贝叶斯估计的比率变化很小。最后,图7b 显示了一个县样本的 50% 的区间估计值(之所以选择这个区间估计值是因为很难在一个图中显示所有的 3071)。较小的县提供的信息较少,因此具有较宽的后验间隔。
图 6(a)肾癌死亡率 $y_i/(10n_j)$ 与人口规模 $n_j$ 。(b)按 $\log_{10}$ 人口的比例重新绘制,以更清楚地看到数据。模式来自数据的离散性 $(n_j=0,1,2,...)$ 。
图 7(a) 贝叶斯估计的肾癌后期平均死亡率,$$E(\theta|y_i)=\frac{20+y_i}{430000+10n_j}$$
美国的 3071 个县(b)100 个县(j)的后验中位数和 $\theta$j 的 50% 区间。$y$ 轴上的比例与图 6 中的图不同。
构建先验分布
现在我们退后一步,讨论基础利率的 Gamma (20,430000)先验分布。正如我们在介绍模型时所讨论的,我们选择伽马分布是为了便于数学计算。我们现在解释如何从数据中估计两个参数 $\alpha$ ,$\beta$ 以匹配观察到的癌症死亡率$$\frac{y_j}{10n_j}$$的分布。使用数据来设置先验分布似乎是不恰当的,但是我们认为这是一个有用的近似于我们首选的分层建模方法 ,其中分布参数如 $\alpha$ ,$\beta$ 在这个例子中被视为待估计的未知数。
在该模型下,任何县 $j$ 的观测计数 $y_j$ 来自预测, $p(y_i)=\int p(y_j|\theta_j)p(\theta_j)d\theta_j$ ,在这种情况下 是 $Neg-bin(\alpha,\frac{\beta}{10n_j})$ 。从附录A 中,我们可以
找到这个分布的平均值和方差:
$$
E(y_j)=10n_j\frac{\alpha}{\beta}
$$
$$
var(y_j)=10n_j\frac{\alpha}{\beta}+(10n_j)^2\frac{\alpha}{\beta^2}
$$
将观察到的平均值和方差与其预期匹配,并求解 $\alpha$ 和 $\beta$ 产生先验分布的参数。实际的计算更复杂。
图8美国 3071 个县的年龄调整后肾癌死亡率 $\frac{y_i}{10n_j}$ 的经验分布,以及基础癌症死亡率 $\theta_j$的$Gamma (20,430000)$先验分布。
因为我们必须处理年龄调整的问题,而且利用比率 $\frac{y_i}{10n_j}$ 的平均值和方差也更有效率:
$$
E(\frac{y_j}{10n_j})=\frac{\alpha}{\beta}
$$
$$
var(\frac{y_j}{10n_j})=\frac{1}{10n_j}\frac{\alpha}{\beta}+\frac{\alpha}{\beta^2}
$$
在处理了年龄调整之后,我们把观察到的时刻和理论上的时刻等同起来,设定的值的均值 $\frac{y_j}{10n_j}$ 比 $\frac{\alpha}{\beta}$ 。值的方差为 $\frac{y_j}{10n_j}$ 比 $E(\frac{1}{10n_j})\frac{\alpha}{\beta}+\frac{\alpha}{\beta^2}$ ,使用值的样本平均值 $\frac{1}{10n_j}$ 代替了 $E(\frac{1}{10n_j})$ .
图8显示了原始癌症发病率的经验分布,以及基础癌症发病率 $\theta_j$ 的估计 $Gamma(20,430000)$ 的先验分布。原始比率的分布要广泛得多,这是有意义的,因为它们包括泊松变异性以及各县之间的变化。
在这个例子中,我们的先验分布是合理的,但是这种通过匹配矩来构造它的方法有点草率,一般情况下可能难以应用。这个模型可以改进的一个更重要的方法是包括可以预测癌症发病率变化的县级信息。这将使该模型朝着分层泊松回归方向发展。