Bayesian analysis of conjugate hierarchical models#

Part 1 · Stage 5 · 🏛️ Hierarchical Models · Lesson 035 of 144 · beginner

◀ Previous · Exchangeability and hierarchical models · Next · Normal model with exchangeable parameters ▶ · ↑ Section

Important

✨ AI-generated content. This page was written with the assistance of an AI language model and is provided as a learning aid. Despite careful review, it may still contain mistakes, omissions, or out-of-date information. Whether you are new to the topic, a team lead, or a senior practitioner, treat it as a starting point rather than an authoritative reference: read it critically and independently verify anything you act on (code, commands, figures, and factual claims) against official documentation and primary sources before relying on it.

The rat tumours#

The canonical worked example: \(J = 71\) historical laboratory experiments, each reporting \(y_j\) rats developing tumours out of \(n_j\). Some experiments are tiny. Estimating each rate separately gives wild answers (an experiment with 0 of 5 suggests a rate of exactly zero); pooling them all into one rate denies that experiments differ. Exchangeability says: model the rates as drawn from a common population distribution.

The three-level model#

Binomial likelihood, Beta population, hyperprior on the Beta’s parameters:

\[y_j \mid \theta_j \sim \mathrm{Binomial}(n_j, \theta_j), \qquad \theta_j \mid \alpha, \beta \sim \mathrm{Beta}(\alpha, \beta), \qquad (\alpha, \beta) \sim p(\alpha, \beta) .\]

Because the Beta is conjugate to the binomial, the conditional posterior of each group is immediate:

\[\theta_j \mid \alpha, \beta, y \;\sim\; \mathrm{Beta}(\alpha + y_j,\; \beta + n_j - y_j),\]

which makes \(\mathrm{E}[\theta_j \mid \alpha, \beta, y]\) the familiar weighted average of the group’s own rate \(y_j/n_j\) and the population mean \(\alpha/(\alpha+\beta)\) — with the population now acting as a prior estimated from all 71 experiments. Small experiments are pulled hard toward the population; large ones barely move.

Two levels, two computations#

The hyperparameters are handled by marginalising the group parameters analytically (conjugacy again), leaving a two-dimensional marginal posterior \(p(\alpha, \beta \mid y)\) that can be evaluated on a grid — exactly the bioassay trick. Then draw \((\alpha, \beta)\), and draw each \(\theta_j\) from its conditional Beta. Modern practice simply hands the whole thing to a sampler:

import pymc as pm
with pm.Model():
    # hyperprior on population mean and "prior sample size"
    mu  = pm.Beta("mu", 1, 1)                    # population mean rate
    kap = pm.HalfNormal("kappa", 50)             # concentration = alpha + beta
    theta = pm.Beta("theta", mu * kap, (1 - mu) * kap, shape=J)
    pm.Binomial("y", n=n, p=theta, observed=y)
    idata = pm.sample()

Choosing the hyperprior#

One trap deserves naming. A flat prior on \((\alpha, \beta)\) is improper and yields an improper posterior — the concentration \(\alpha + \beta\) runs off to infinity. Gelman reparameterises to the population mean \(\alpha/(\alpha+\beta)\) and a transformed concentration, placing a proper prior there. The lesson from Stage 2 returns with teeth: for hierarchical variance and concentration parameters, check propriety, and prefer weakly informative hyperpriors.

Hint