1. Overview and Motivation
A hierarchical model introduces multiple layers of uncertainty.
We have parameters at two levels:
- Lower level (group level):
θ₁, θ₂, …, θⱼ — parameters describing each individual experiment or group. - Upper level (hyperparameters):
φ = (α, β) — parameters describing the population distribution of the θⱼ’s.
We want to make posterior inference for all parameters jointly:
$p(\theta, \phi \mid y) = p(\theta_1, \dots, \theta_J, \alpha, \beta \mid y)$
In practice, this is difficult because of high dimensionality — many θⱼ’s plus the hyperparameters φ.
But in some models, such as the Beta–Binomial, the conjugacy between likelihood and prior makes computation possible with a mix of analytical derivation and numerical integration.
The section uses the rat tumor example, where 71 experiments are performed with binomial data for tumor presence in rats.
2. Hierarchical Model Setup
For each experiment $j = 1, \dots, J$:
$y_j \mid \theta_j \sim \text{Binomial}(n_j, \theta_j)$
$\theta_j \mid \alpha, \beta \sim \text{Beta}(\alpha, \beta)$
Here:
- $y_j$: number of rats with tumors
- $n_j$: total rats in experiment $j$
- $\theta_j$: probability of tumor for that experiment
- $(\alpha, \beta)$: shape parameters of the Beta population distribution
The hyperparameters $(\alpha, \beta)$ themselves have a prior:
$(\alpha, \beta) \sim p(\alpha, \beta)$
Our goal is to find:
$p(\theta, \alpha, \beta \mid y)$
3. Step 1 — Joint Posterior Density (Unnormalized)
By Bayes’ theorem:
$p(\theta, \alpha, \beta \mid y) \propto p(\alpha, \beta)\, p(\theta \mid \alpha, \beta)\, p(y \mid \theta)$
Breakdown of each component:
- Hyperprior:
- $p(\alpha, \beta)$ — prior for hyperparameters (to be specified later).
- Population distribution:
- Each θⱼ follows a Beta(α, β):
- $p(\theta \mid \alpha, \beta) = \prod_{j=1}^J \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \theta_j^{\alpha-1} (1-\theta_j)^{\beta-1}$
- Each θⱼ follows a Beta(α, β):
- Likelihood:
- Each yⱼ given θⱼ is binomial:
- $p(y \mid \theta) = \prod_{j=1}^J \binom{n_j}{y_j} \theta_j^{y_j} (1-\theta_j)^{n_j – y_j}$
- Each yⱼ given θⱼ is binomial:
Ignoring constants independent of θ, α, β (the binomial coefficients), we combine them:
$p(\theta, \alpha, \beta \mid y) \propto p(\alpha, \beta) \prod_{j=1}^J \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\, \theta_j^{\alpha + y_j – 1}\,(1 – \theta_j)^{\beta + n_j – y_j – 1}$
4. Step 2 — Conditional Posterior of θ Given (α, β)
Given α and β, each θⱼ has an independent posterior distribution that’s again Beta (because of conjugacy):
$p(\theta_j \mid \alpha, \beta, y_j) \propto \theta_j^{\alpha + y_j – 1} (1 – \theta_j)^{\beta + n_j – y_j – 1}$
Thus:
$ \theta_j \mid \alpha, \beta, y_j \sim \text{Beta}(\alpha + y_j,\, \beta + n_j – y_j)$
and collectively:
$p(\theta \mid \alpha, \beta, y) = \prod_{j=1}^J \text{Beta}(\theta_j \mid \alpha + y_j,\, \beta + n_j – y_j)$
5. Step 3 — Marginal Posterior of the Hyperparameters (α, β)
We can integrate out all θⱼ’s to obtain $p(\alpha, \beta \mid y)$:
$p(\alpha, \beta \mid y) = \int p(\theta, \alpha, \beta \mid y) \, d\theta$
Substituting the forms above yields:
$p(\alpha, \beta \mid y) \propto p(\alpha, \beta) \prod_{j=1}^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)}$
It can’t be simplified analytically, but is straightforward to compute numerically for given (α, β).
6. Step 4 — Choosing a Hyperprior
We have little prior information about α and β.
However, the Beta distribution’s parameters are not intuitive, so we reparameterize to meaningful quantities:
$m = \frac{\alpha}{\alpha + \beta} \quad \text{(population mean tumor rate)}$
$t = \alpha + \beta \quad \text{(effective “sample size” of the Beta prior)}$
Transform to log scale for numerical stability:
$\text{logit}(m) = \log\frac{m}{1-m} = \log\frac{\alpha}{\beta}, \quad \log(t) = \log(\alpha+\beta)$
We would like a diffuse prior over (logit m, log t).
However, setting a uniform prior in this space causes an improper posterior (diverges as $t \to \infty$).
Hence, they propose a weakly informative but integrable hyperprior:
$p(\alpha, \beta) \propto (\alpha + \beta)^{-5/2}$
This penalizes very large (α + β), preventing infinite precision in the Beta population distribution.
7. Step 5 — Compute the Marginal Posterior p(α, β | y)
We now evaluate:
$p(\alpha, \beta \mid y) \propto (\alpha + \beta)^{-5/2} \prod_{j=1}^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)}$
Practical computation:
- Define a grid in transformed coordinates: $(\log(\alpha / \beta), \log(\alpha + \beta))$ Example:
- initial grid [−2.5, −1] × [1.5, 3], later extended to [−2.3, −1.3] × [1, 5].
- Compute the log posterior at each grid point (for numerical stability):
- Compute log of prior + log of gamma terms.
- Subtract the maximum log density to avoid overflow.
- Exponentiate to get the unnormalized posterior density.
- Normalize the grid so total probability sums to 1.
The resulting contour plot shows a roughly symmetric posterior centered near:
$(\log(\alpha / \beta), \log(\alpha + \beta)) \approx (-1.75, 2.8)$
which corresponds to approximately:
$(\alpha, \beta) \approx (2.4, 14.0)$
or a population mean tumor rate $\alpha / (\alpha + \beta) \approx 0.146$.
8. Step 6 — Posterior Moments of Hyperparameters
Using the grid as a discrete approximation to the posterior, compute expectations like:
$E[\alpha \mid y] = \sum_{\text{grid}} \alpha \, p(\log(\alpha / \beta), \log(\alpha + \beta) \mid y)$
$E[\beta \mid y] = \sum_{\text{grid}} \beta \, p(\log(\alpha / \beta), \log(\alpha + \beta) \mid y)$
Results:
$E[\alpha \mid y] \approx 2.4, \quad E[\beta \mid y] \approx 14.3$
very close to the modal values from the contour plot, confirming approximate symmetry.
This step captures posterior uncertainty in (α, β), unlike a single point estimate.
9. Step 7 — Sampling from the Full Joint Posterior
We now simulate from:
$p(\alpha, \beta, \theta_1, \dots, \theta_J \mid y)$
Algorithm:
Repeat L times (e.g., L = 1000):
- Draw hyperparameters
- Draw (α, β) from the discrete posterior grid $p(\alpha, \beta \mid y)$.
- Draw group parameters
- For each experiment j:
- $\theta_j \mid \alpha, \beta, y_j \sim \text{Beta}(\alpha + y_j,\, \beta + n_j – y_j)$
- For each experiment j:
- (Optional) Draw predictive data $\tilde{y}_j \sim \text{Binomial}(n_j, \theta_j)$.
This yields samples:
$(\theta_1^{(l)}, \dots, \theta_J^{(l)}, \alpha^{(l)}, \beta^{(l)}) \quad \text{for } l = 1, \dots, L$
These represent the full posterior distribution.
10. Step 8 — Results and Interpretation
Posterior summaries:
- For each experiment j, compute:
- Posterior median of $\theta_j$
- 95% posterior interval of $\theta_j$
- Plot against the observed tumor proportion $y_j / n_j$.
Interpretation of the plot:
- The 45° line shows no pooling (raw data proportion).
- Each posterior median lies between its raw proportion and the overall mean ≈ 0.14 — indicating shrinkage.
- Smaller studies (low nⱼ):
- Greater shrinkage toward population mean
- Wider uncertainty intervals
- Larger studies:
- Less shrinkage
- Narrower intervals
- The result is similar to empirical Bayes (plug-in α, β), but slightly wider intervals, reflecting extra uncertainty in α, β.
11. Key Takeaways
- Conjugacy makes analytic simplification possible.
The Beta–Binomial hierarchy allows closed-form conditional posteriors. - Numerical integration over hyperparameters (2D grid) gives a practical way to obtain $p(\alpha, \beta \mid y)$.
- Full Bayesian approach propagates uncertainty in hyperparameters, unlike plug-in (empirical Bayes) methods.
- Shrinkage arises naturally from the hierarchical structure:
- $\theta_j \text{ posterior mean } = E[\theta_j \mid y] \text{ lies between group estimate and population mean.}$
- The posterior plot visually shows the balance between individual information and population sharing — small nⱼ get “pulled” more.
12. Summary of Equations
| Concept | Formula | Comment |
|---|---|---|
| Likelihood | $y_j \mid \theta_j \sim \text{Binomial}(n_j, \theta_j)$ | Observation model |
| Population (prior) | $\theta_j \mid \alpha, \beta \sim \text{Beta}(\alpha, \beta)$ | Group-level model |
| Joint posterior | $p(\theta, \alpha, \beta \mid y) \propto p(\alpha, \beta)\prod_j p(\theta_j \mid \alpha, \beta) p(y_j \mid \theta_j)$ | |
| Conditional posterior | $\theta_j \mid \alpha, \beta, y_j \sim \text{Beta}(\alpha + y_j, \beta + n_j – y_j)$ | |
| Marginal posterior | $p(\alpha, \beta \mid y) \propto p(\alpha, \beta)\prod_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)}$ | |
| Hyperprior | $p(\alpha, \beta) \propto (\alpha + \beta)^{-5/2}$ | |
| Posterior grid space | $(\log(\alpha / \beta), \log(\alpha + \beta))$ | Transform for numerical evaluation |
| Expected hyperparameters | $E[\alpha \mid y] = 2.4,\; E[\beta \mid y] = 14.3$ | From posterior averaging |
In essence:
The Beta–Binomial hierarchical model is the simplest full example of Bayesian partial pooling.
It shows how conjugacy, grid-based marginalization, and simulation can combine to yield complete posterior inference — including shrinkage, uncertainty propagation, and predictive capability — even with dozens of groups.