Hierarchical Bayesian Models

Definition and motivation

A hierarchical Bayesian model (also called a multilevel model) is a probabilistic model in which parameters themselves have probability distributions whose parameters (hyperparameters) are also given distributions. Hierarchical structure encodes the idea of exchangeability and partial pooling: groups are similar but not identical. Rather than estimating each group's parameter independently (no pooling) or assuming all groups share a single parameter (complete pooling), hierarchical models let the data determine how much information to borrow across groups.

Why care: partial pooling reduces estimation variance while controlling bias. Consider stores estimating conversion rates: a tiny store with 2 purchases in 20 visits would have an unreliable raw estimate 0.10. A hierarchical model borrows strength from the chain-level distribution to pull (shrink) that estimate toward the chain mean, often giving dramatically better predictive performance.

Core formalism (two-level example)

The canonical two-level hierarchical model for group-specific parameters $\theta_j$ and observations $y_{ji}$ is:

Here $j=1,\dots,J$ indexes groups (stores), $i=1,\dots,n_j$ indexes observations in group $j$. Exchangeability across groups means that, before seeing data, group labels don't matter and we treat $\theta_1,\dots,\theta_J$ as iid from $p(\theta\mid\phi)$.

Concrete numeric toy: Beta-Binomial intuition

If each store has conversion counts $y_j\sim\text{Binomial}(n_j, \theta_j)$ and we place $\theta_j\sim\text{Beta}(\alpha,\beta)$ with hyperparameters $(\alpha,\beta)$, then posterior inference for each $\theta_j$ is conjugate:

Numeric example: a small store had $y=2$ successes in $n=20$ visits. If chain-level hyperparameters are $\alpha=5,\beta=95$ (prior mean 0.05, variance small), the posterior mean is

which is shrunk toward 0.05 from the raw rate $2/20=0.10$.

Contrast with complete pooling: estimating a single $\theta$ for all stores would ignore between-store variability; no pooling would give $\hat{\theta}_j=y_j/n_j$ with high variance for small $n_j$. The hierarchical model, by introducing hyperparameters and a hyperprior, adapts the amount of pooling to the data.

Philosophy and exchangeability

In [Conjugate Priors] we learned how conjugacy gives closed-form posteriors for single-level models like Beta-Binomial and Normal-Normal. Hierarchical models extend this: conjugacy may allow analytic conditional posteriors for lower-level parameters given hyperparameters. But typically we integrate over hyperparameters (or sample them in MCMC) to get full posterior uncertainty.

Exchangeability justifies the iid assumption on $\theta_j\mid\phi$: if groups are exchangeable, the joint prior over $\theta_{1:J}$ factorizes into iid draws from the same mixture. Partial pooling is the Bayesian implementation of the bias-variance tradeoff: the posterior for each $\theta_j$ is shrunk toward the common prior mean, with the extent of shrinkage determined by the relative sizes of within-group noise and between-group variability.

Key quantities to watch: the hyperparameters (e.g. variance of $\theta_j$ across groups) control shrinkage; the posterior predictive distribution integrates over both $\theta_j$ and $\phi$ and is what you use for predicting new observations or new groups.

Summary sentence

Hierarchical Bayesian models let you model group heterogeneity while borrowing strength across groups via hyperpriors; this produces principled shrinkage estimators that often dominate non-hierarchical alternatives in prediction and mean squared error.

Normal-Normal hierarchical model (canonical analytic case)

A central worked analytic case is the Normal-Normal hierarchical model with known observation variance. This is where the familiar shrinkage formula arises in closed form and links directly to [Conjugate Priors]. Model:

Let $\bar y_j=\frac{1}{n_j}\sum_{i=1}^{n_j} y_{j,i}$ be the group mean and $\sigma^2/n_j$ its sampling variance. Conditional on hyperparameters $(\mu,\tau^2)$ the posterior for $\theta_j$ is Gaussian (by conjugacy):

where the shrinkage weight is

(We prefer the intuitively clearer form below.) Solving algebra shows the posterior mean can be written as:

This $\kappa_j\in(0,1)$ is the effective weight on the data vs the group mean. When $n_j\to\infty$, $\kappa_j\to1$ and we rely on the data; when $\tau^2\to0$ (no between-group variability) $\kappa_j\to0$ and every group collapses to the common mean $\mu$ (complete pooling).

Concrete numeric example (normal-normal)

Suppose $\sigma^2=4$ (observation SD $=2$), $\tau^2=1$ (between-group SD $=1$), $n_j=5$ observations in group $j$, and group sample mean $\bar y_j=10$, hyper-mean $\mu=8$. Then

so posterior mean is

Compare: raw mean 10, pooled mean 8, shrunk estimate 9.11.

Derivation sketch

Conditional posterior precision = prior precision + data precision. Prior precision $=1/\tau^2$, data precision $=n_j/\sigma^2$. Thus posterior mean is a precision-weighted average:

This is the canonical "shrinkage" formula. In [Conjugate Priors] you saw the Normal-Normal conjugacy; here it extends across groups with a hyperprior on $\mu$ and $\tau^2$ if desired.

James–Stein connection (brief, illuminating)

The James–Stein estimator is a (frequentist) shrinkage estimator for multivariate normal means that dominates the MLE in mean squared error for dimension $\ge 3$. Hierarchical Bayes with an appropriate hyperprior produces shrinkage estimators with similar behavior. A heuristic connection: empirical Bayes or hierarchical Bayes produces estimators that shrink toward a common mean; in large-dimension limits the Bayes estimator approximates a form of James–Stein shrinkage.

Empirical Bayes (ML-II) approximation

If $\mu$ and $\tau^2$ are unknown, one approach is empirical Bayes: estimate them by maximizing the marginal likelihood of the observed $\bar y_j$'s. Marginally,

So we can compute the (log) marginal likelihood and obtain MLEs $\hat\mu,\hat\tau^2$; plug them into the posterior mean formula above to get an empirical Bayes point estimate for each $\theta_j$. This ML-II step is often fast and connects to [MCMC] because you can use these estimates to initialize or inform MCMC priors.

Numeric example of empirical Bayes (method of moments)

Given $J$ observed $\bar y_j$ with known $\sigma^2$ and $n_j$, method-of-moments estimates give

Then estimate between-group variance via

Concrete numbers: if we have three groups with $\bar y=(9,11,10)$, $n_j=5$ each, $\sigma^2=4$, then sample variance of $\bar y$ is $1$; average sampling variance is $\sigma^2/n_j=0.8$, so method-of-moments gives $\hat\tau^2=\max(0,1-0.8)=0.2$.

Summary

This section gave the analytic backbone: conjugate normal-normal hierarchies yield explicit shrinkage weights $\kappa_j$; empirical Bayes estimates hyperparameters by marginalizing the lower-level parameters, producing fast approximate shrinkage estimators. In practice, for nonconjugate or complex multilevel regressions we use MCMC (see [MCMC]) with either centered or non-centered parameterizations to ensure efficient sampling.

Multilevel regression: varying intercepts and slopes

Hierarchical modeling becomes especially powerful when combined with regression: we let intercepts and slopes vary by group. A simple varying-intercept, varying-slope model for outcome $y_{ij}$ with predictor vector $x_{ij}$ is:

where $\alpha_j$ is group-specific intercept, $\beta_j$ group-specific slope vector, and $\Sigma_{\mathrm{group}}$ is a covariance matrix encoding how intercepts and slopes vary and covary across groups. This fully multilevel regression allows partial pooling in both intercepts and slopes and allows complex cross-group shrinkage: groups with scarce data will have their $\alpha_j$ and $\beta_j$ pulled toward population means $(\mu_\alpha,\mu_\beta)$.

Priors on covariance matrices: separation strategy and LKJ

Priors on $\Sigma_{\mathrm{group}}$ are critical. A common strategy (the "separation" prior) decomposes

where $\tau$ is the vector of standard deviations for each random effect and $R$ a correlation matrix. Priors: put half-Cauchy or half-normal priors on components of $\tau$ (weakly informative), and an LKJ prior on $R$ (Lewandowski–Kurowicka–Joe) with concentration parameter $\eta$ favoring identity when $\eta>1$. Example numeric choices: $\tau_k \sim \mathrm{HalfCauchy}(0,2)$, $R\sim\mathrm{LKJ}(\eta=2)$.

Centered vs non-centered parameterizations (important for MCMC efficiency)

In [MCMC] we learned that parameterization affects convergence and mixing. For hierarchical models, two common parameterizations are:

• •Centered: directly model $\theta_j\sim N(\mu,\tau^2)$ and sample $\theta_j$.
• •Non-centered: reparametrize $\theta_j = \mu + \tau\,\eta_j$, with $\eta_j\sim N(0,1)$.

Why it matters: when group-level variance $\tau$ is small and data are informative, centered parameterization mixes well; when $\tau$ is large and/or data per group are small, non-centered often mixes much better. A practical rule: try both or use adaptive diagnostics. Concrete numeric illustration: if $\tau=0.01$ and data are abundant, centered is fine; if $\tau=5$ and $n_j$ small, non-centered reduces funnel-shaped posterior geometry and avoids slow MCMC.

Posterior predictive checks and partial pooling

The full Bayesian workflow uses posterior predictive distributions to judge fit and to predict. For a given group $j$ the posterior predictive for a new observation $y_{j\,\text{new}}$ is

In practice we approximate this by MCMC draws: for each draw of $\phi$ and $\theta_j$ simulate $y_{\text{new}}$. Posterior predictive checks diagnose over/under-pooling: if pooled model systematically mispredicts heterogeneity, consider richer group-level structure (e.g., allow covariates to explain group differences).

Empirical Bayes vs full Bayes: tradeoffs

In empirical Bayes (ML-II) we estimate hyperparameters $\phi$ by marginal likelihood and then condition on $\hat\phi$. This often gives good point estimates and is fast for large hierarchical problems. However, it underestimates hyperparameter uncertainty: it treats $\hat\phi$ as known, which can understate posterior variance for predictions, especially with few groups. Full Bayes places hyperpriors on $\phi$ and samples it in MCMC, which captures uncertainty but costs more computation.

Numeric demonstration of marginal likelihood (for normal-normal)

Recall marginal of $\bar y_j$ is $\bar y_j\sim N(\mu, \tau^2 + \sigma^2/n_j)$. For known $\sigma^2$ the log marginal likelihood is

For example, with $J=3$, $\bar y=(9,11,10)$, $n_j=5$, $\sigma^2=4$, we numerically maximize $\ell$ over $\mu,\tau^2$ to get empirical Bayes estimates used in the shrinkage formula.

Practical advice and diagnostics

• •Standardize predictors before hierarchical modeling; it stabilizes priors and sampling.
• •Use weakly informative priors (e.g., half-normal for scales) to regularize and avoid pathologies.
• •Compare centered vs non-centered parameterizations; non-centered is often superior when group-level variance is small relative to observation noise.
• •Check posterior predictive distributions for each group to identify misfit.

Summary

This section extended the core shrinkage idea to multilevel regressions with random effects, explained priors on covariance matrices, and provided practical sampling and parameterization advice. Non-centered parameterization and proper hyperpriors are essential tools for efficient, robust multilevel Bayesian inference.

Store-level vs chain-level parameter sharing (concrete business example)

Imagine an online retailer with many stores (or merchants) belonging to a chain. We want to estimate store conversion rates $\theta_j$ and predict revenue. Options:

• •No pooling: estimate each store separately, $\hat\theta_j = y_j/n_j$, high variance for small stores.
• •Complete pooling: estimate a single chain-level $\theta$, ignores useful heterogeneity.
• •Hierarchical model: $\theta_j\sim p(\theta\mid\phi)$ with chain-level hyperprior on $\phi$. This automatically shares information across stores: large stores dominate inference for themselves; tiny stores are shrunk toward the chain-level mean.

Concrete Beta-Binomial numeric example

Suppose store A: $y_A=2$, $n_A=20$; store B: $y_B=40$, $n_B=400$. With hyperparameters $\alpha=10,\beta=190$ (prior mean 0.05), posterior means are

• •store A: $E[\theta_A]=(10+2)/(200+20)=12/220\approx0.0545$ (shrunk from 0.10 to 0.0545),
• •store B: $E[\theta_B]=(10+40)/(200+400)=50/600\approx0.0833$ (raw 0.10 shrunk only slightly).

This demonstrates adaptive pooling: big-store estimates stay near their data; small-store estimates move toward the chain-level prior.

Empirical Bayes in large-scale problems

For thousands of stores one may use empirical Bayes to estimate hyperparameters quickly via marginal likelihood. A common pipeline:

1. 1)Aggregate store-level sufficient statistics (e.g., counts for Beta-Binomial or sample means and variances for Normal-Normal).
2. 2)Optimize marginal likelihood for hyperparameters (e.g., $\alpha,\beta$ or $\mu,\tau^2$).
3. 3)Compute posterior summaries for each store conditional on the estimated hyperparameters.

This ML-II approach scales because it reduces high-dimensional posterior inference to a moderate-dimensional optimization; it often suffices for point predictions in industry systems. But be mindful: it ignores hyperparameter uncertainty, which can matter for decision-making, particularly when the number of groups $J$ is small.

Real-world use cases

• •A/B testing across multiple sites: hierarchical models estimate site-specific treatment effects while pooling information to improve detection power.
• •Education: estimate school-level effects (test scores) with partial pooling across schools; used in value-added models.
• •Healthcare: hospital-level outcome rates with case-mix adjustment; hierarchical models stabilize hospital estimates and reduce false positives.
• •Small-area estimation: census and survey domains with sparse data benefit from Bayesian small-area models.

Downstream tasks enabled

• •Better prediction: partial pooling often reduces out-of-sample error compared to no pooling.
• •Decision-making with uncertainty: full Bayesian posterior distributions for group-level parameters enable decision rules that account for both within-group and between-group uncertainty.
• •Hierarchical modeling of complex structures: nested hierarchies (patients within clinicians within hospitals), crossed random effects, and multi-membership models are natural extensions.

Limitations and modeling caveats

• •Model misspecification: if groups are not exchangeable (e.g., unmodeled covariates explain between-group differences), naive pooling may be inappropriate. Solution: include group-level covariates or model non-exchangeable structure.
• •Outliers and heavy tails: a Normal prior for $\theta_j$ may be too light-tailed; consider Student-t for robustness.
• •Over-shrinkage: overly-informative hyperpriors can force excessive shrinkage. Use weakly informative priors on hyperparameters and check sensitivity.

Computational practice and modern tooling

Modern probabilistic programming (Stan, PyMC, Turing, etc.) makes hierarchical Bayes accessible. Use non-centered parameterizations, LKJ priors for correlations, and adapt step sizes/warmups in HMC. For very large problems consider variational inference or empirical Bayes for speed; validate approximations against full MCMC on subsamples.

Summary

Hierarchical Bayesian models operationalize the balance between pooling and group-specific estimation. In practice, they are the default approach for multi-group inference when you want improved estimates, calibrated uncertainty, and principled borrowing of strength across related units.