Hierarchical Models

Hierarchical Model

A hierarchical model (or multilevel model) is a Bayesian model in which unit-level parameters ${\theta }_1,\dots ,{\theta }_J$ share a common distribution indexed by a hyperparameter $\mathbf{\varphi }$:
$\mathbf{\varphi }\sim p(\mathbf{\varphi }),{\theta }_j\mid \mathbf{\varphi }\stackrel{\text{iid}}{\sim }p(\theta \mid \mathbf{\varphi }),y_j\mid {\theta }_j\sim p(y\mid {\theta }_j),j=1,\dots ,J.$
The top-level prior $p(\mathbf{\varphi })$ is the hyperprior; the distribution $p(\theta \mid \mathbf{\varphi })$ is the population distribution; $p(y\mid {\theta }_j)$ is the unit-level likelihood. The graphical structure is a tree: hyperparameter at the root, unit parameters as children, observations as leaves. Each observation depends on its own ${\theta }_j$, but all ${\theta }_j$‘s are coupled through their shared ancestor $\mathbf{\varphi }$.

Posterior structure. The joint posterior factorizes as
$p(\mathbf{\varphi },\mathbf{\theta }\mid \mathbf{y})\propto p(\mathbf{\varphi }){\prod }_{j=1}^{J}p({\theta }_j\mid \mathbf{\varphi })p(y_j\mid {\theta }_j).$
Conditional on $\mathbf{\varphi }$, the ${\theta }_j$‘s are independent with individual posteriors obtained from Bayes’ rule within each unit. Marginalizing out $\mathbf{\varphi }$ couples the units, which is why hierarchical models need numerical techniques (Gibbs sampling, Hamiltonian Monte Carlo, or empirical Bayes) even when the within-unit updates are conjugate.

Headline result — partial pooling. The posterior mean $\mathbb{E}[{\theta }_j\mid \mathbf{y}]$ is pulled (shrunk) from the unit-level MLE ${\hat{\theta }}_j$ toward the population-level mean, by an amount increasing in the population prior’s concentration and decreasing in the unit’s own sample size. Units with little data borrow strength from the pool; units with abundant data are barely affected. This (Partial Pooling) property is the defining advantage of hierarchical models: it beats both the no-pooling fit (separate MLE per unit, which overfits small units) and the complete-pooling fit (ignore unit identity, which underfits heterogeneity).

Two approaches to inference. Full Bayes puts a hyperprior $p(\mathbf{\varphi })$ and integrates by MCMC. Empirical Bayes estimates $\hat{\mathbf{\varphi }}$ by maximizing the marginal likelihood $p(\mathbf{y}\mid \mathbf{\varphi })={\prod }_j\int p(y_j\mid \theta )p(\theta \mid \mathbf{\varphi })d\theta$ and plugs in $\hat{\mathbf{\varphi }}$ in the unit-level posteriors. Empirical Bayes underestimates uncertainty about $\mathbf{\varphi }$ but is computationally cheap and gives nearly identical posterior means.

James-Stein phenomenon. In the Normal-Normal hierarchy, the empirical-Bayes posterior mean for each ${\theta }_j$ is the James-Stein estimator — an explicit shrinkage estimator that strictly dominates the MLE under squared-error loss whenever $J\ge 3$. This remains the clearest example of a frequentist estimator outperformed by a Bayesian one uniformly in the parameter.

Prerequisites

• 3.1 Bayesian Framework — prior, likelihood, posterior.
• 3.5.1 Beta-Binomial Model — conjugate building block for the binomial hierarchy.
• 3.5.2 Normal-Normal Model — conjugate building block for the Gaussian hierarchy.
• 3.7 Marginal Likelihood — hyperparameter estimation via evidence maximization.
• 1.2.1 Maximum Likelihood Estimation — the “no-pooling” fit that hierarchical models improve on.

Intuition

Suppose we observe cancer mortality rates in each of the $J$ US counties. Some counties are huge (Los Angeles, $n_j\sim 10^7$) and some are tiny ($n_j\sim 10^3$). If we compute a separate maximum-likelihood rate ${\hat{\theta }}_j=y_j/n_j$ for each county, the smallest counties will dominate both the top and the bottom of any ranked list — they have the highest sampling variance, so extreme values are most likely to appear there by chance. Conversely, if we pool all counties together and estimate a single $\hat{\theta }$, we lose all geographic variation and obscure real disparities. Neither extreme is correct.

The hierarchical answer is to acknowledge that each county has its own true rate ${\theta }_j$, but that these rates are drawn from a shared population distribution. Large counties estimate their own rate almost independently; small counties borrow strength from the population. The tool that decides how much each county borrows is Bayes’ rule applied within the unit, using the population distribution as the prior. The population itself has to be estimated from the data too — that is what the hyperprior (full Bayes) or the marginal likelihood (empirical Bayes) provides.

The central tension is between heterogeneity and pooling. Unit parameters are assumed exchangeable at some level — we do not know which county will be high and which low before seeing data — so a shared prior on the ${\theta }_j$‘s is reasonable. But the units are genuinely different; rural and urban counties have different exposures. Pooling too much erases the signal; pooling too little erases the safety. The hyperparameter controls the rate of pooling, and the model learns that rate from data by trading off within-group variance against between-group variance. This is the Bayesian shadow of a random-effects ANOVA.

Properties

(Partial Pooling Beats Both Extremes) Hierarchical posterior means strictly outperform no-pooling and complete-pooling under any reasonable risk measure when heterogeneity is moderate.

Under no pooling (separate MLE per group), small units have huge variance. Under complete pooling (single rate), bias dominates because unit-level variation is ignored. The hierarchical posterior mean interpolates: it shrinks each ${\hat{\theta }}_j$ toward the grand mean $\hat{\mu }$ by a factor that depends on the group’s sample size and the population variance. Under squared-error loss and the Normal-Normal hierarchy, the optimal shrinkage factor is $B_j={\sigma }_j^2/({\sigma }_j^2+{\tau }^2)$, where ${\sigma }_j^2$ is the sampling variance for group $j$ and ${\tau }^2$ is the population variance; small groups have $B_j$ near 1 (shrink fully) and large groups have $B_j$ near 0 (essentially MLE).

(Beta-Binomial Hierarchy) The standard binomial hierarchy has Beta population distribution and county-level Beta-Binomial updating.

Let $y_j\mid {\theta }_j\sim \text{Binomial}(n_j,{\theta }_j)$ and ${\theta }_j\mid \alpha ,\beta \sim \text{Beta}(\alpha ,\beta )$. Conditional on $(\alpha ,\beta )$, the unit-level posterior is ${\theta }_j\mid y_j,\alpha ,\beta \sim \text{Beta}(\alpha +y_j,\beta +n_j-y_j)$ with posterior mean
$\mathbb{E}[{\theta }_j\mid y_j,\alpha ,\beta ]=\frac{\alpha +y_j}{\alpha +\beta +n_j}=\underset{\text{data weight}}{\underbrace{\frac{n_j}{\alpha +\beta +n_j}}}\cdot {\hat{\theta }}_j+\underset{\text{shrinkage weight}}{\underbrace{\frac{\alpha +\beta }{\alpha +\beta +n_j}}}\cdot \frac{\alpha }{\alpha +\beta }.$
The hyperparameters $(\alpha ,\beta )$ are themselves estimated by marginal likelihood: $p(\mathbf{y}\mid \alpha ,\beta )={\prod }_j\text{BetaBin}(y_j;n_j,\alpha ,\beta )$, which is maximized numerically. Maximum-likelihood $\hat{\alpha },\hat{\beta }$ give the empirical-Bayes posterior for each ${\theta }_j$.

(Normal-Normal Hierarchy and James-Stein) The Gaussian hierarchy’s empirical-Bayes posterior mean equals the James-Stein shrinkage estimator.

Let $y_j\mid {\theta }_j\sim N({\theta }_j,{\sigma }_j^2)$ with known sampling variances ${\sigma }_j^2$, and ${\theta }_j\mid \mu ,{\tau }^2\stackrel{\text{iid}}{\sim }N(\mu ,{\tau }^2)$. The conditional posterior is
${\theta }_j\mid y_j,\mu ,{\tau }^2\sim N\left({\tau }^2({\tau }^2+{\sigma }_j^2{)}^{-1}{y}_j+{\sigma }_j^2({\tau }^2+{\sigma }_j^2{)}^{-1}\mu ,{\tau }^2{\sigma }_j^2/({\tau }^2+{\sigma }_j^2)\right).$
With known ${\sigma }_j=\sigma$ and empirical-Bayes estimates $\hat{\mu }=\bar{y}$ and ${\hat{\tau }}^2=\max (0,s_y^2-{\sigma }^2)$ (where $s_y^2$ is the sample variance of the $y_j$‘s), the posterior mean of ${\theta }_j$ is
${\hat{\theta }}_j^{\text{EB}}=\bar{y}+\left(1-\frac{(J-2){\sigma }^2}{{\sum }_j(y_j-\bar{y}{)}^2}\right)(y_j-\bar{y}),$
the classical James-Stein estimator. When $J\ge 3$, this dominates the MLE ${\hat{\theta }}_j=y_j$ under total squared-error risk ${\sum }_j\mathbb{E}({\hat{\theta }}_j-{\theta }_j{)}^2$ for every true $\mathbf{\theta }\in {\mathbb{R}}^J$ — a “strict dominance” result that stunned the frequentist world when Stein proved it in 1956.

(Empirical Bayes vs Full Bayes) EB underestimates hyperparameter uncertainty but recovers the same posterior means.

Empirical Bayes plugs in $\hat{\mathbf{\varphi }}=\arg {\max }_{\mathbf{\varphi }}p(\mathbf{y}\mid \mathbf{\varphi })$ and computes unit-level posteriors as $p({\theta }_j\mid y_j,\hat{\mathbf{\varphi }})$. Full Bayes integrates $p({\theta }_j\mid \mathbf{y})=\int p({\theta }_j\mid y_j,\mathbf{\varphi })p(\mathbf{\varphi }\mid \mathbf{y})d\mathbf{\varphi }$. Both give nearly identical posterior means, because the integrand is nearly symmetric around $\hat{\mathbf{\varphi }}$. The difference appears in the variance of the unit-level posterior: EB treats $\hat{\mathbf{\varphi }}$ as known, so its credible intervals are too narrow by a factor reflecting hyperparameter uncertainty. Full Bayes (via MCMC) recovers calibrated intervals at the cost of substantially more computation. See 3.7 Marginal Likelihood for hyperparameter estimation details.

(Gibbs-Sampling Friendly) Conjugate population distributions yield tractable Gibbs updates.

In the Beta-Binomial and Normal-Normal hierarchies, the full conditional for each ${\theta }_j$ given $\mathbf{\varphi }$ is a closed-form conjugate posterior, and the conditional for $\mathbf{\varphi }$ given $\{{\theta }_j\}$ is often also closed form (e.g., Inverse-Gamma for ${\tau }^2$ in the Normal hierarchy). This makes hierarchical models a canonical setting for Gibbs samplers; see 3.11.2 Gibbs Sampling. Non-conjugate hierarchies require Metropolis-within-Gibbs steps or Hamiltonian Monte Carlo for the hyperparameter level.

(Variance Stabilization via the Delta Method) Heterogeneous sampling variances can be homogenized by an appropriate transformation before applying the Normal hierarchy.

Raw sample proportions ${\hat{p}}_j=y_j/n_j$ have sampling variance $p_j(1-p_j)/n_j$ that depends on $p_j$ itself. The arcsine transformation ${\hat{\varphi }}_j=2\arcsin \sqrt{{\hat{p}}_j}$, derived by the delta method applied to stabilize binomial variance, satisfies $\mathrm{Var}({\hat{\varphi }}_j)\approx 1/n_j$ independent of $p_j$. Similarly for Poisson counts, ${\hat{\psi }}_j=2\sqrt{y_j}$ has variance $\approx 1$ independent of the rate. After transformation, Normal-Normal hierarchical machinery (including James-Stein) applies directly. This is how the original James-Stein paper’s kidney-cancer and baseball-batting-average analyses were performed.

(Exchangeability Assumption is Load-Bearing) Hierarchical models assume unit-level parameters are exchangeable; structural covariates should be incorporated.

Exchangeability of ${\theta }_1,\dots ,{\theta }_J$ means no covariate is known a priori to separate them. If we know that ${\theta }_j$ depends systematically on a covariate $x_j$ (e.g., urban vs. rural), a simple exchangeable hierarchy will over-shrink those relationships toward the grand mean. The fix is a regression hierarchy: ${\theta }_j=\alpha +\beta x_j+{ϵ}_j$ with ${ϵ}_j\sim N(0,{\tau }^2)$, replacing the shared mean $\mu$ with a covariate-dependent mean. This extends directly to multilevel regression and is how GLMMs (generalized linear mixed models) are built.

(Three-Level and Deeper Hierarchies Follow the Same Pattern) Adding more layers produces nested partial pooling.

A three-level hierarchy has ${\theta }_{jk}\sim p(\theta \mid {\varphi }_j)$, ${\varphi }_j\sim p(\varphi \mid \psi )$, $\psi \sim p(\psi )$. A student test-score model might have students nested in classrooms nested in schools. The pooling pattern recurses: each ${\theta }_{jk}$ is pulled toward its classroom mean ${\varphi }_j$, which is pulled toward the school mean $\psi$, which is pulled toward nothing (top-level prior). The graphical structure stays a tree. Full Bayes via MCMC is the standard inference method.

Examples

Example 1: Kidney Cancer Mortality (Beta-Binomial Hierarchy)

The US has roughly 3{,}100 counties. Kidney cancer mortality rates ${\hat{p}}_j=y_j/n_j$ vary enormously across counties, and the counties with the highest rates are also the counties with the lowest rates — small, rural counties with low population dominate both tails purely by sampling noise. (See Gelman et al. BDA3 Figure 2.7.)

Model: $y_j\mid {\theta }_j\sim \text{Binomial}(n_j,{\theta }_j)$, ${\theta }_j\mid \alpha ,\beta \stackrel{\text{iid}}{\sim }\text{Beta}(\alpha ,\beta )$.
Empirical-Bayes hyperparameters: maximizing $p(\mathbf{y}\mid \alpha ,\beta )={\prod }_j\text{BetaBin}(y_j;n_j,\alpha ,\beta )$ over $(\alpha ,\beta )$ yields $\hat{\alpha }\approx 20$, $\hat{\beta }\approx 400,000$ (numbers illustrative). The prior mean is $20/400,020\approx 5\times 10^{-5}$, the national baseline rate.
Unit-level posterior mean:
$\mathbb{E}[{\theta }_j\mid y_j,\hat{\alpha },\hat{\beta }]=\frac{\hat{\alpha }+y_j}{\hat{\alpha }+\hat{\beta }+n_j}.$
A county with $n_j=1000$ and $y_j=1$ has MLE $10^{-3}$ but posterior mean $\approx 21/401,001\approx 5.2\times 10^{-5}$ — essentially the national rate, because the data evidence is barely audible against the hierarchical prior. A county with $n_j=10^6$ and $y_j=40$ has MLE $4\times 10^{-5}$ and posterior mean $\approx 60/1,400,020\approx 4.3\times 10^{-5}$ — mostly its own data. The hierarchical model reshuffles the ranking of counties by pooling out the noise in small ones while leaving large ones nearly alone.

Example 2: Eight Schools (Normal-Normal Hierarchy)

The classic “eight schools” dataset (Rubin 1981) has $J=8$ schools each reporting the effect of an SAT coaching program, with treatment-effect estimate $y_j$ and known standard error ${\sigma }_j$:

$j$	$y_j$	${\sigma }_j$
1	28	15
2	8	10
3	$-3$	16
4	7	11
5	$-1$	9
6	1	11
7	18	10
8	12	18

Model: $y_j\mid {\theta }_j\sim N({\theta }_j,{\sigma }_j^2)$, ${\theta }_j\mid \mu ,\tau \stackrel{\text{iid}}{\sim }N(\mu ,{\tau }^2)$, with $\mu$ and $\tau$ unknown.
Full Bayesian fit (MCMC): posterior mean of $\tau \approx 6$, indicating only modest between-school variance; posterior means of the ${\theta }_j$‘s shrink substantially toward the grand mean $\bar{y}\approx 8$. School 1, with MLE 28, has posterior mean $\approx 11$ — the model disbelieves the extreme value.
Interpretation: the raw MLEs over-dramatize between-school differences. After shrinkage, the coaching program’s estimated effect is broadly positive but modest; extreme schools are attributed to sampling noise.

Example 3: James-Stein Beats the MLE

Suppose $J=10$ independent units with $y_j\sim N({\theta }_j,1)$, and the true ${\theta }_j$ are drawn from $N(0,{\tau }^2)$. The MLE ${\hat{\theta }}_j^{\text{MLE}}=y_j$ has total expected squared error ${\sum }_j\mathbb{E}({\hat{\theta }}_j-{\theta }_j{)}^2=J=10$, regardless of $\tau$.

The James-Stein estimator
${\hat{\theta }}_j^{\text{JS}}=\left(1-\frac{J-2}{{\sum }_k{y}_k^2}\right)y_j$
shrinks all $y_j$‘s toward zero by a data-dependent factor. Its expected total squared error is strictly less than $J$ for every $\mathbf{\theta }$, improving most when the true $\mathbf{\theta }$ is near zero and vanishingly when it is far. A simple simulation with $\tau =1$ gives expected error $\approx 5$ for JS and $=10$ for MLE — a factor-of-two improvement from pure shrinkage logic, with no prior expertise about $\mathbf{\theta }$ needed.

The Bayesian interpretation: the JS formula is the empirical-Bayes posterior mean under the (nearly) Normal-Normal hierarchy with $\mu =0$ and ${\tau }^2$ estimated from ${\sum }_k{y}_k^2$. Stein’s dominance theorem says this hierarchical Bayesian improvement over the MLE is not a lucky coincidence — it is structural.

Derivation

Goal: derive the shrinkage weight $B_j={\sigma }_j^2/({\sigma }_j^2+{\tau }^2)$ in the Normal-Normal hierarchy, conditional on known $(\mu ,{\tau }^2)$.

Starting from the unit-level joint for ${\theta }_j$ given $y_j$ and hyperparameters:
$\log p({\theta }_j\mid y_j,\mu ,{\tau }^2)=-\frac{1}{2{\sigma }_j^2}(y_j-{\theta }_j{)}^2-\frac{1}{2{\tau }^2}({\theta }_j-\mu {)}^2+\text{const}.$
Both terms are quadratic in ${\theta }_j$, so by precision addition (3.5.2 Normal-Normal Model), the posterior is Gaussian with precision ${\kappa }_j={\sigma }_j^{-2}+{\tau }^{-2}$ and mean
$\mathbb{E}[{\theta }_j\mid y_j,\mu ,{\tau }^2]=\frac{{\sigma }_j^{-2}{y}_j+{\tau }^{-2}\mu }{{\sigma }_j^{-2}+{\tau }^{-2}}=\frac{{\tau }^2{y}_j+{\sigma }_j^2\mu }{{\sigma }_j^2+{\tau }^2}.$
Rearranging as a convex combination:
$\mathbb{E}[{\theta }_j\mid y_j,\mu ,{\tau }^2]=(1-B_j)y_j+B_j\mu ,B_j=\frac{{\sigma }_j^2}{{\sigma }_j^2+{\tau }^2}.$
Shrinkage weight $B_j\to 1$ (pull fully to $\mu$) when ${\sigma }_j^2\gg {\tau }^2$, i.e. unit noisier than population spread. $B_j\to 0$ (ignore the pool) when ${\sigma }_j^2\ll {\tau }^2$, i.e. unit precise relative to heterogeneity.

James-Stein plug-in: estimate $\mu$ by the weighted mean $\hat{\mu }=({\sum }_j{y}_j/{\sigma }_j^2)/({\sum }_{j}1/{\sigma }_j^2)$ (for equal ${\sigma }_j=\sigma$, this is $\bar{y}$), and estimate ${\tau }^2$ by ${\hat{\tau }}^2=\max (0,s_y^2-{\sigma }^2)$ where $s_y^2$ is the sample variance of the $y_j$‘s (for equal ${\sigma }_j$). With equal ${\sigma }_j=\sigma$ and $J\ge 3$, the (slightly adjusted) estimate ${\hat{\tau }}^2$ substituted into the shrinkage factor recovers exactly $1-(J-2){\sigma }^2/{\sum }_j(y_j-\bar{y}{)}^2$, yielding the classical James-Stein formula.

Empirical Bayes vs Full Bayes coincidence. Both produce the same posterior mean up to $O(1/J)$ corrections because the marginal likelihood $p(\mathbf{y}\mid \mathbf{\varphi })$ is approximately quadratic in $\mathbf{\varphi }$ around its maximum (BvM applied at the hyperparameter level), so integrating against a diffuse hyperprior returns the mode to leading order.

See Also

• 3.1 Bayesian Framework — underlying Bayesian machinery.
• 3.5.1 Beta-Binomial Model — unit-level building block for binomial hierarchies.
• 3.5.2 Normal-Normal Model — unit-level building block for Gaussian hierarchies.
• 3.5.3 Gamma-Poisson Model — unit-level building block for count-data hierarchies (e.g., disease incidence per county).
• 3.7 Marginal Likelihood — hyperparameter estimation via evidence maximization (empirical Bayes).
• 3.11 MCMC — standard inference engine for full Bayes hierarchies.
• 3.11.2 Gibbs Sampling — natural sampler when population distributions are conditionally conjugate.
• 3.8 Bayesian Regression — regression extensions yield multilevel (mixed-effects) models.
• 1.2.1 Maximum Likelihood Estimation — no-pooling baseline that hierarchical models improve on.