贝叶斯数据分析（9）——贝叶斯层级模型

本笔记基于Bayesian Data Analysis Third Edition， Andrew Gelman et. al. 学习编写，由于是英文教材，可能在学习过程中一些翻译或者内容有误，如有问题或错误，可以发送至邮箱[email protected]反馈。在学习本书中，需要有一定数学基础或者数理统计的基础，并且存在很多需要计算的场景，对于每一条定理，从头开始的证明会让你更明白每一步是如何实现的。

我们之前提及的贝叶斯模型里，我们基本不考虑数据来源的异质性，但实际上，很多贝叶斯模型同时包含一组相互关联的参数。例如，多个医院各有一个生存率，多个学校各有一个教学干预效应，多个临床试验各有一个治疗效应。分别估计这些参数会丢掉组间能够共享的信息；强行令它们相等又会抹去真实差异。层级模型在这两个极端之间建立概率模型，让数据决定信息共享的程度。

在贝叶斯层级模型里，共同分布需要可交换性作为依据，未知超参数会让各组在边缘上相互依赖，这种依赖最终表现为部分汇聚和更完整的不确定性。

层级模型的概率结构

设第 $j$ 个组的数据为 $y_j$，对应的组参数为 $\omega_j$，全部组参数记为 $\omega=(\omega_1,\ldots,\omega_J)$。组参数来自由超参数 $\phi$ 控制的共同分布。最基本的层级结构为

在给定 $\phi$ 时，各组参数通常设为条件独立；给定各自的 $\omega_j$ 时，各组数据也通常条件独立。因此联合分布可以写成

数据通过似然更新组参数，全部组的数据又通过组参数共同更新 $\phi$。于是第 $j$ 组的推断会受到其他组数据的影响。这种影响不是额外添加的校正，而是联合后验分布本身的结果：

三种汇聚方式

当我们面对多个中心的数据时，常常有三个处理方式：不汇聚，即我们对单独中心进行分析；完全汇聚，即我们不考虑组间的异质性，认为其完全相同；部分汇聚，即考虑组件之间的异质性，这也正是我们层级模型所做的，同时也会更加复杂。

需要注意的是，部分汇聚不是固定比例的平均。样本信息精确的组通常保留更多自身特征，样本少或标准误大的组会更明显地向总体中心收缩。组间方差越大，模型允许的组间差异越大，收缩也越弱。

层级模型把组估计收缩到总体中心。收缩强度由组内信息、组间变异和超参数不确定性共同决定；总体中心本身也由数据估计。所有组由此在同一个联合模型中相互借用信息。

可交换性

定义

如果对任意指标置换 $pi$，联合分布均满足

则称 $(\omega_1,\ldots,\omega_J)$ 是可交换的。这个定义表达的是先验联合分布对标签置换不变。它不表示各组在现实中完全相同，也不要求观测值相等。

可交换性的建模依据是：在进入模型的信息中，没有变量足以预先区分各组参数的分布。实验发生在不同地点或时间并不会自动否定可交换性。只有当地点、时间、实验设计或其他变量能够系统地区分组参数时，模型才应将这些信息编码为分组结构或协变量。

可交换性与独立性的区别

独立性要求联合分布可以分解：

可交换性要求联合分布在置换后不变。两者约束的是不同性质。

• 独立但不同分布的变量一般不可交换。
• 独立同分布变量一定可交换。
• 可交换变量不一定独立。
• 层级模型中的组参数通常在给定超参数后独立，在边缘分布下却相互依赖。

层级模型最常用的可交换表示为

第一行是条件独立，第二行对未知的 $\phi$ 做了平均。所有 $\omega_j$ 共享同一个 $\phi$，所以知道某个组参数会改变对 $\phi$ 的判断，并进一步改变其他组参数的分布。

设在给定 $\phi$ 时，$\omega_j$ 与 $\omega_k$ 条件独立且条件均值均为 $m(\phi)$。边缘协方差由全协方差公式得到

$$\begin{aligned} \operatorname{Cov}(\omega_j,\omega_k) &=\mathbb{E}\!\left[\operatorname{Cov}(\omega_j,\omega_k\mid\phi)\right]+\operatorname{Cov}\!\left(\mathbb{E}[\omega_j\mid\phi],\mathbb{E}[\omega_k\mid\phi]\right)\\ &=0+\operatorname{Cov}(m(\phi),m(\phi))\\ &=\operatorname{Var}(m(\phi)). \end{aligned}$$

只要超参数仍有不确定性，右侧通常大于零。条件独立没有推出边缘独立。数据更新 $\phi$ 后，各组的后验分布也通过 $\phi$ 联系起来。

层级混合产生的相关通常为正，但可交换性本身并不要求正相关。无放回抽样得到的序列可以可交换且负相关。有限维可交换分布还可能受到总和固定等约束，不能写成独立同分布混合。

de Finetti 表示的边界

对无限可交换序列，在适当正则条件下，de Finetti 定理允许将其表示为条件独立同分布模型的混合：

这里的无限性条件不能省略。一个有限可交换分布未必能够扩展为无限可交换序列。实际建模中，条件独立同分布的混合仍然是表达可交换性的主要工具，但它是一种模型选择，不是有限维可交换性的同义定义。

部分可交换与条件可交换

当组具有已知协变量 $x_j$ 时，直接假设 $\omega_j$ 完全可交换可能过于粗糙。可以改为

此时相同或相近协变量条件下的组具有可比较的分布。实验室、地区和时间等离散结构可以形成更高一层的部分可交换模型；连续信息则可以进入组层回归。

对可交换性的有效质疑应落到可观测结构上。若某个变量能够系统解释组间差异，就把它写入模型。仅仅指出每个组都不同，并不足以否定共同分布，因为层级模型本来就允许 $\omega_j$ 彼此不同。

层级模型中的贝叶斯计算

层级模型的参数维数通常较高，但概率结构给出了清晰的计算分解。简单共轭模型可以先解析边缘化超参数，再条件抽取组参数；一般的非共轭模型则通常使用 HMC 对联合后验进行抽样。

一个基本模拟顺序为

第二步在许多模型中可以按组并行，因为

解析边缘化的价值不仅是降低维数。它还揭示了组参数共享信息的具体通道，并能用于检查通用采样器的结果。

实际计算中的 HMC

当超参数的边缘后验不能解析求出，或者模型包含非共轭似然、大量连续组参数、组层回归及相关结构时，实际分析通常直接对联合后验 $p(\omega,\phi\mid y)$ 使用 HMC。相比随机游走式 Metropolis，HMC 在高维且参数相关的后验中能够进行更充分的移动，样本自相关通常更低，单位计算时间内得到的有效样本量也往往更高。HMC 适用于连续且对数密度可微的参数；离散潜变量通常需要先积分掉或采用其他处理方式。

Stan 和 brms 默认使用 HMC 的 NUTS 变体，因此后文的 pybrms 层级模型实际由 HMC/NUTS 完成抽样。HMC 的算法、NUTS 和调参已经在上一节介绍，本章只强调它在层级模型中的使用。层级漏斗仍可能造成发散和低效抽样，因此实际拟合时需要结合非中心化参数化，并检查发散、最大树深、能量诊断、$\widehat R$ 与有效样本量。

在层级模型中使用比例符号时，需要确认被省略的项对当前变量确实是常数。例如，从 $p(\omega,\phi\mid y)$ 推导 $p(\phi\mid y)$ 时，$p(\omega\mid\phi,y)$ 的归一化常数依赖 $\phi$，不能随意删去。最稳妥的做法是显式积分，或者使用完整的规范化条件分布进行相除。

Beta-Binomial 层级模型

我们引入一个例子进行阐述，对于多个组别中，我们使用一种药物诱导若干只大鼠使其产生肿瘤，我们需要知道这个药物成功的概率是多少。因此在这里面，我们的组别即是我们层级模型里的异质性考虑因素。

大鼠肿瘤数据的模型

设第 $j$ 个实验中有 $n_j$ 只大鼠，其中 $y_j$ 只出现肿瘤，肿瘤概率为 $\omega_j$。模型为

数据层描述每个实验内部的二项抽样误差，参数层则让不同实验的肿瘤率来自同一个 Beta 总体分布。$\alpha$ 与 $\beta$ 保留在概率模型中，但直接用它们解释总体位置与组间差异并不方便，因此先改写为总体均值 $m$ 和集中度 $\kappa$。

用 $m$ 与 $\kappa$ 表示总体分布

定义

当 $\alpha,\beta>0$ 时，$(\alpha,\beta)$ 与 $(m,\kappa)$ 一一对应，其中 $0<m<1$、$\kappa>0$。这两个参数的作用可以直接从 Beta 分布的矩看出：

$m$ 是总体分布的均值，表示各实验肿瘤率共同围绕的中心；$\kappa$ 是集中度，决定不同实验的 $\omega_j$ 围绕 $m$ 的紧密程度。在 $m$ 固定时，$\kappa$ 越大，组间差异越小；$\kappa$ 越小，Beta 分布越分散，组间异质性越强。

$\kappa$ 有时被称为 Beta 先验的有效样本量，因为 $\alpha$ 与 $\beta$ 可以分别看作成功和失败方向上的先验伪计数，总量为 $\alpha+\beta=\kappa$。它并不是真实观测数，只是在共轭更新中以与样本量相同的代数方式决定先验权重。因此，在理论上，将 $\kappa$ 称为集中度更准确。

联合后验与组参数的条件后验

忽略只依赖观测数据的二项系数，联合后验为

因此给定 $(\alpha,\beta)$ 后，各组条件后验相互独立（即上面等式的后面一部分）：

条件后验均值可以写成部分汇聚的形式：

其中 m&=\frac{\alpha}{\alpha+\beta},&\kappa&=\alpha+\beta,\\

数据比例 $y_j/n_j$ 的权重为 $n_j/(n_j+\kappa)$，总体中心 $m$ 的权重为 $\kappa/(n_j+\kappa)$。小样本实验向总体中心收缩得更多。换句话说，我们可以认为小样本的实验它对总的中心的权重或者贡献较小。

从组参数积分到超参数边缘后验

超参数 $(\alpha,\beta)$ 控制所有组的总体分布。为了只研究这两个超参数，需要将每个组的潜在概率 $\omega_j$ 从联合模型中积分掉。单组积分给出 Beta-Binomial 边缘分布：

利用 Beta 函数 $B(a,b)=\Gamma(a)\Gamma(b)/\Gamma(a+b)$，单组边缘分布为

$$\begin{aligned} p(y_j\mid\alpha,\beta) &=\int_0^1 p(y_j\mid\omega_j)p(\omega_j\mid\alpha,\beta)\,d\omega_j\\ &=\binom{n_j}{y_j}\frac{1}{B(\alpha,\beta)} \int_0^1\omega_j^{\alpha+y_j-1}(1-\omega_j)^{\beta+n_j-y_j-1}\,d\omega_j\\ &=\binom{n_j}{y_j}\frac{B(\alpha+y_j,\beta+n_j-y_j)}{B(\alpha,\beta)}\\ &=\binom{n_j}{y_j} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \frac{\Gamma(\alpha+y_j)\Gamma(\beta+n_j-y_j)}{\Gamma(\alpha+\beta+n_j)}. \end{aligned}$$

各组在给定 $(\alpha,\beta)$ 后条件独立，因此完整的超参数边缘似然是各组边缘分布的乘积。再乘以超先验，才得到超参数的边缘后验：

其中二项系数只依赖观测数据，在写后验核时已经省去。积分后的似然仍然同时包含总体中心和组间异质性的信息：$m$ 决定肿瘤率的总体位置，$\kappa$ 决定各组能否被一个很窄的 Beta 分布共同解释。现在还有一个问题，在于我们的 $p(\alpha,\beta)$, 之前我们可以用均匀分布进行一个非信息先验的假设，但这里，我们不使用这个简单的非信息先验：

超先验、参数变换与 Jacobian

这一部分会涉及三套坐标。$(\alpha,\beta)$ 是模型坐标，$(m,\kappa)$ 是解释坐标，还有我们后面计算时再使用的

其中 $\eta,\lambda\in\mathbb{R}$，适合建立数值积分网格。坐标改变后，同一小块区域所含的概率必须不变，但单位面积对应的密度会改变，Jacobian 就是这两种面积尺度之间的换算因子。

第一步：在可解释参数上定义超先验

为了计算方便，令

其中，$\kappa= \alpha+\beta$, 并在 $(m,s)$ 上采用均匀超先验，即 $p(m,s)\propto1$。这里的 $s$ 是总体分布集中度的反向尺度：$\kappa$ 越大，$s$ 越小。先从 $(m,s)$ 变到 $(m,\kappa)$，再变到 $(\alpha,\beta)$，两次 Jacobian 为

$$\begin{aligned} \left|\frac{\partial s}{\partial\kappa}\right| &=\frac12\kappa^{-3/2},\\ p(m,\kappa) &=p(m,s)\left|\frac{\partial s}{\partial\kappa}\right| \propto\kappa^{-3/2},\\ \left|\frac{\partial(\alpha,\beta)}{\partial(m,\kappa)}\right| &=\left| \begin{matrix} \kappa&m\\-\kappa&1-m \end{matrix} \right|=\kappa,\\ \left|\frac{\partial(m,\kappa)}{\partial(\alpha,\beta)}\right| &=\kappa^{-1},\\ p(\alpha,\beta) &=p(m,\kappa) \left|\frac{\partial(m,\kappa)}{\partial(\alpha,\beta)}\right|\\ &\propto\kappa^{-3/2}\kappa^{-1}=(\alpha+\beta)^{-5/2}. \end{aligned}$$

这样得到的超先验是

$p(m,s)\propto1$ 在 $s>0$ 上是non-proper超先验。我们选择 $\kappa^{-5/2}$ 的尾部衰减，正是为了可以得到一个proper的后验分布。相反，若直接在 $(\eta,\lambda)$ 上采用无限范围的均匀先验，后验在 $\kappa\to\infty$ 时不可积。实际分析也可以改用proper、weak informative的超先验，并进行后验分布是否proper的检查。

第二步：把密度变到计算坐标

由 $m=\operatorname{logit}^{-1}(\eta)$、$\kappa=\exp(\lambda)$，从 $(\eta,\lambda)$ 到 $(\alpha,\beta)$ 的 Jacobian 可以沿着 $(\eta,\lambda)\to(m,\kappa)\to(\alpha,\beta)$ 分解：

$$\begin{aligned} \left|\frac{\partial(m,\kappa)}{\partial(\eta,\lambda)}\right| &=m(1-m)\kappa,\\ \left|\frac{\partial(\alpha,\beta)}{\partial(\eta,\lambda)}\right| &=\left|\frac{\partial(\alpha,\beta)}{\partial(m,\kappa)}\right| \left|\frac{\partial(m,\kappa)}{\partial(\eta,\lambda)}\right|\\ &=\kappa\cdot m(1-m)\kappa\\ &=\alpha\beta. \end{aligned}$$

因此，计算坐标上的超先验和边缘后验核分别为

其中

后验模拟与收缩

共轭模型的后验抽样可以分两步完成：

对 $(\alpha,\beta)$ 的后验不确定性进行平均，会使 $\omega_j$ 的边缘后验比固定超参数的经验贝叶斯近似更宽。后者只使用 $(\widehat\alpha,\widehat\beta)$，遗漏了总体分布估计误差。

图中的 $45^\circ$ 线对应不汇聚估计 $y_j/n_j$。后验中位数整体向总体肿瘤率移动，数据较少的实验区间更宽，收缩也更明显。

正态层级模型

正态层级模型是部分汇聚这种方式更直接的例子，也是许多 Meta 分析和随机效应模型的基础。

数据层与参数层

设第 $j$ 个组的均值统计量为 $\bar y_{\cdot j}$，其已知抽样方差为 $\sigma_j^2$。模型为

其中，第三条里，我们设立一个非信息均匀超先验分布给 $p(\mu|\tau)$, 则认为是一个常数

$\mu$ 是组参数总体的中心，$\tau$ 是组间标准差。$\sigma_j$ 描述第 $j$ 组估计的抽样误差，$\tau$ 描述真实组参数之间的异质性，两者不能混为一项。

联合后验核为

给定超参数时的组参数后验

在给定 $(\mu,\tau)$ 后，$\omega_j$ 的后验仍为正态：

其中

这一点和我们在多参数模型里提到的很类似，这一结果可由配方直接得到。只保留与 $\omega_j$ 有关的项：

$$\begin{aligned} \log p(\omega_j\mid\mu,\tau,y_j) &\doteq-\frac{(\bar y_{\cdot j}-\omega_j)^2}{2\sigma_j^2} -\frac{(\omega_j-\mu)^2}{2\tau^2}\\ &=-\frac12\left(\frac1{\sigma_j^2}+\frac1{\tau^2}\right) \left[\omega_j- \frac{\bar y_{\cdot j}/\sigma_j^2+\mu/\tau^2} {1/\sigma_j^2+1/\tau^2}\right]^2+C, \end{aligned}$$

其中 $C$ 与 $\omega_j$ 无关。平方项的中心给出后验均值，二次项系数的倒数给出后验方差。

定义数据权重

则

当 $\sigma_j$ 小时，本组估计精确，$\lambda_j$ 接近 1；当 $\tau$ 小时，组参数更接近共同中心，$\lambda_j$ 接近 0。这给出了部分汇聚的精确权重。

超参数后验

因为

所以卷积后

也可以用矩母函数或直接积分验证：

$$\begin{aligned} p(\bar y_{\cdot j}\mid\mu,\tau) &=\int p(\bar y_{\cdot j}\mid\omega_j) p(\omega_j\mid\mu,\tau)\,d\omega_j\\ &=\int\mathcal N(\bar y_{\cdot j}\mid\omega_j,\sigma_j^2) \mathcal N(\omega_j\mid\mu,\tau^2)\,d\omega_j\\ &=\mathcal N(\bar y_{\cdot j}\mid\mu,\sigma_j^2+\tau^2). \end{aligned}$$

$\sigma_j^2$ 与 $\tau^2$ 相加，是因为前者来自观测误差，后者来自真实组参数的离散，两层变异相互独立。

因此联合超参数后验为

在给定 $\tau$ 并令 $p(\mu\mid\tau)\propto1$，即均匀分布后，有

其中

当我们进一步积分掉 $\mu$，得到 $\tau$ 的一维边缘后验：

上式中的 $V_\mu^{1/2}$ 来自对共同均值的高斯积分。把边缘后验写成二次型后，

$$\begin{aligned} &\sum_{j=1}^{J}\frac{(\bar y_{\cdot j}-\mu)^2}{\sigma_j^2+\tau^2}\\ &\quad=\sum_{j=1}^{J}\frac{(\bar y_{\cdot j}-\widehat\mu)^2}{\sigma_j^2+\tau^2} +\frac{(\mu-\widehat\mu)^2}{V_\mu},\\ &\int_{-\infty}^{\infty} \exp\!\left[-\frac{(\mu-\widehat\mu)^2}{2V_\mu}\right]d\mu =\sqrt{2\pi V_\mu}. \end{aligned}$$

与 $\tau$ 无关的 $\sqrt{2\pi}$ 可以略去，而 $V_\mu$ 依赖 $\tau$，必须保留。

后验模拟

根据上面流程的反向推理，我们可以得到一个具体的后验模拟抽样流程，先抽取超参数 $\tau$ 的分布，然后抽取超参数 $\mu$ 的分布，最后得到了联合超参数的分布后，抽取组参数 $\omega$ 的分布：

低维情况下可以在 $\tau$ 网格上计算归一化密度并用逆 CDF 抽样。通用概率编程则直接在联合参数空间中用 HMC 抽样。两种方法应给出相容的后验结果。

Example

八校数据来自八个独立的 SAT 辅导实验。每个学校提供估计效应 $y_j$ 及近似已知的标准误 $\sigma_j$：

不汇聚分析把 A 校效应中心放在 28，但其他学校的数据削弱了这种解释。完全汇聚分析则把八个真实效应强行设为相等，忽略了学校之间可能存在的差异。正态层级模型采用

对固定的 $\tau$，将 $\mu$ 的后验不确定性也积分进去，可得

其中

第二个方差公式来自全方差公式。给定 $(\mu,\tau)$ 的条件方差是 $V_j$，条件均值中 $\mu$ 的系数是 $1-\lambda_j$，因此

$$\begin{aligned} \operatorname{Var}(\omega_j\mid\tau,y) &=\mathbb E[\operatorname{Var}(\omega_j\mid\mu,\tau,y)\mid\tau,y] +\operatorname{Var}[\mathbb E(\omega_j\mid\mu,\tau,y)\mid\tau,y]\\ &=V_j+(1-\lambda_j)^2V_\mu(\tau). \end{aligned}$$

当 $\tau=0$ 时，八个学校完全汇聚；当 $\tau\to\infty$ 时，估计逐渐接近不汇聚结果。数据中 $p(\tau\mid y)$ 在零附近较高，但对正值仍保留明显质量。仅把 $\tau$ 固定在后验众数零会产生过强汇聚，因为它把关于异质性的全部不确定性删除了。

层级模型还能直接计算复合后验量，例如最大效应 $\max_j\omega_j$、学校 A 的效应超过学校 C 的概率，以及任一学校效应超过某个实际阈值的概率。这些量只需在每次联合后验抽样中计算相应函数。

极端观测向总体中心收缩能够缓和选择最大值带来的偏差。这项调整已经包含在可交换层级模型的联合后验中，因注意在结果报告时，无须再把它作为传统多重比较校正附加到结果上。

两类后验预测

层级模型需要区分已有组的新观测与新组的新观测。

已有组的新观测

对已有组 $j$，未来观测共享已经学习过的 $\omega_j$：

在正态模型中，若新观测方差为 $\widetilde\sigma_j^2$，则

新组的新观测

对一个尚未出现的新组，先从总体分布抽取新参数 $\widetilde\omega$，再抽取观测：

正态层级模型中，新组真实效应的后验预测方差为

新组预测比总体均值 $\mu$ 的后验更宽，因为它额外包含组间异质性。若还要预测新组的观测统计量，则再加入该统计量的抽样方差。

报告总体平均效应并不能替代新组效应预测。前者描述总体中心的不确定性，后者描述一个具体但尚未观测的可交换组可能落在何处。应用决策通常更关心后者。

层级模型用于 Meta 分析

设第 $j$ 个临床试验有治疗组与对照组。令 $y_{1j}$、$n_{1j}$ 为治疗组事件数和总人数，$y_{0j}$、$n_{0j}$ 为对照组对应数量。对数优势比的估计为

其大样本抽样方差近似为

于是可以使用与八校模型相同的正态层级结构：

$\mu$ 是可交换试验总体的平均对数优势比，$\tau$ 描述研究间异质性，$\omega_j$ 是第 $j$ 个已观测研究的真实效应。新研究效应 $\widetilde\omega$ 的后验预测比 $\mu$ 的后验更宽。

某个 $2\times2$ 表存在零单元格时，上述对数优势比与方差近似失效。简单连续性校正会影响小样本结果。更稳妥的模型直接使用二项似然，并在 logit 尺度上建立层级结构。

可交换性在 Meta 分析中是一项实质性假设。若研究在剂量、结局定义、患者构成或偏倚风险上存在系统差异，这些变量应进入组层回归，或者先划分为部分可交换的研究集合。研究筛选机制也会影响总体解释，不能由层级模型自动修复。

组间标准差 $\tau$ 的先验

在上面的例子里，我们了解到了 $\tau$ 决定部分汇聚强度。当组数较少或数据支持 $\tau$ 接近零时，即完全汇聚时，先验会明显影响后验。但目前为止，我们还没有说过我们的 $\tau$ 应该选用什么先验，在这里，我们进一步讨论 $\tau$ 这个参数，因为其代表了组间标准差，而其先验的选择很重要。在我们之前提到非信息先验的选择时，我们通常会在对数尺度或者一般尺度上使用均匀分布，但在这里我们需要来看看设立这些先验都分别会有什么问题。

对数均匀先验的问题

在 $\log\tau$ 上均匀等价于

当 $\tau\to0$ 时，正态层级模型的边缘似然趋近有限的非零值。因此后验在零附近的积分包含

后验不是proper的。数据无法完全排除组间方差为零，所以不能在该区域放入无限先验质量。

在 $\tau$ 上均匀

$p(\tau)\propto1$ 在零附近可积。对正态层级模型，当 $\tau\to\infty$ 时，

但需要注意的是，只有 $J>2$ 时右尾积分才有限。

由 $\sigma_j^2+\tau^2\asymp\tau^2$ 可得乘积项为 $\tau^{-J}$；同时

$$\begin{aligned} V_\mu(\tau) &=\left\{\sum_{j=1}^{J}(\sigma_j^2+\tau^2)^{-1}\right\}^{-1} \asymp\frac{\tau^2}{J},\\ V_\mu(\tau)^{1/2}\prod_{j=1}^{J}(\sigma_j^2+\tau^2)^{-1/2} &\asymp\tau\cdot\tau^{-J}=\tau^{-(J-1)}. \end{aligned}$$

指数项在该极限下趋近常数，所以尾部阶数由上式决定。

在 $\tau^2$ 上均匀意味着 $p(\tau)\propto\tau$，右尾更厚，需要 $J>3$ 才能得到正常后验，而且更容易产生过大的异质性估计。

小参数 inverse-gamma 先验

对 $\tau^2$ 使用 $\operatorname{Inv\text{-}Gamma}(\epsilon,\epsilon)$ 看似接近无信息，但当 $\epsilon\to0$ 时不存在稳定的正常后验极限。只要似然在 $\tau=0$ 附近仍有质量，后验就可能对 $\epsilon$ 极为敏感。把 $\epsilon$ 从 $0.01$ 改为 $0.001$ 并不是单纯地减少先验信息，它会改变近零区域的形状。

half-Cauchy 与 half-$t$

half-Cauchy 先验的密度为

它在零附近有限，右尾较厚，并由尺度 $A$ 表达异质性的合理数量级。更一般的 half-$t$ 先验也常用于组层标准差。尺度不应脱离结果变量的单位设置。

八校数据中，在 $\tau$ 上均匀的先验仍能得到可用后验，而 inverse-gamma 先验明显约束了近零区域。只保留三个学校时，均匀先验产生很长的右尾；具有有限尺度的 half-Cauchy 去除了缺乏实际意义的超大异质性，同时保留了高似然区域。

弱信息先验仍然是模型的一部分。拟合后应比较先验与后验、检查不同合理尺度下的敏感性，并把实际不可接受的 $\tau$ 范围在拟合前表达出来。用极宽的有限区间替代不正常先验，可能也只是把问题推到一个更遥远的边界，仍然没有解决问题。

使用 pybrms 拟合层级模型

pybrms 是对 R 包 brms 的 Python 接口。其公开版本停留在 0.0.33，发布于 2020 年，内部依赖 rpy2、brms 和旧版 PyStan。

例子里的正态层级模型

brms 公式中的 se(se) 表示每个效应估计有已知标准误且不另加残差尺度，(1 | school) 表示学校层随机截距。它对应

import pandas as pd
import pybrms
import arviz as az

schools = pd.DataFrame({
    'school': ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H'],
    'effect': [28, 8, -3, 7, -1, 1, 18, 12],
    'se': [15, 10, 16, 11, 9, 11, 10, 18],
})

fit_8schools = pybrms.fit(
    formula='effect | se(se) ~ 1 + (1 | school)',
    data=schools,
    family='gaussian',
    priors=[
        ('normal(0, 50)', 'Intercept'),
        ('student_t(3, 0, 25)', 'sd'),
    ],
    chains=4,
    iter=4000,
    warmup=2000,
    seed=20260812,
)

idata_8schools = az.from_pystan(fit_8schools)
print(az.summary(idata_8schools, round_to=2))
az.plot_trace(idata_8schools)
az.plot_posterior(idata_8schools)

student_t(3, 0, 25) 作用于标准差类参数时由 brms 截断到正半轴，相当于 half-$t$ 先验。尺度 25 继承了八校问题中以 SAT 分数为单位的量级，不应直接搬到其他结果变量。

二项分布层级模型

pybrms 中最直接的多组二项模型使用 logit-normal 组层分布：

import pandas as pd
import pybrms

rat_data = pd.DataFrame({
    'experiment': ['e01', 'e02', 'e03', 'e04', 'e05', 'e06'],
    'tumor': [0, 1, 2, 3, 1, 4],
    'total': [20, 20, 20, 20, 14, 20],
})

fit_rat = pybrms.fit(
    formula='tumor | trials(total) ~ 1 + (1 | experiment)',
    data=rat_data,
    family='binomial',
    priors=[
        ('normal(-2, 1.5)', 'Intercept'),
        ('student_t(3, 0, 1)', 'sd'),
    ],
    chains=4,
    iter=4000,
    warmup=2000,
    seed=20260812,
)

print(fit_rat)

计算结果的检查

在实际应用里，我们构建层级模型后对其进行采样需要检查下面几点：

1. 检查所有总体参数、组层标准差和代表性组参数的 $\widehat R$ 与有效样本量。
2. 检查发散、最大树深和能量诊断；漏斗形后验在小 $\tau$ 时尤其容易暴露参数化问题。
3. 比较中心化与非中心化参数化。数据较弱且组数少时，非中心化通常更稳定。
4. 做先验预测，确认 $\mu$、$\tau$ 和观测尺度上的先验能够生成合理数据。
5. 做后验预测，分别检查组内拟合、组间离散和极端组。
6. 对 $\tau$ 的先验尺度进行敏感性分析，并报告它对收缩和新组预测的影响。