Hierarchical Bayesian Regression with PyMC: When Groups Share Strength

A multi-line insurer writes auto, home, commercial property, and a dozen other policy types under one roof. Some lines see thousands of claims a year; others might see 50. Every actuary faces the same dilemma: train a separate pricing model for each line and the small ones are pure noise, or pool everything together and pretend a warehouse fire looks like a fender bender. Either way, you lose.

Hierarchical Bayesian regression offers a third way. Each group gets its own parameters, but those parameters are drawn from a shared population distribution. Groups with plenty of data stay close to their own estimates. Groups with little data get "pulled" toward the population average, borrowing statistical strength from the larger groups. This effect is called shrinkage, and it's one of the most elegant ideas in statistics.

By the end of this post, you'll build a hierarchical Bayesian regression model in PyMC, compare it against pooled and unpooled alternatives, and see shrinkage in action on synthetic insurance data.

Let's Build It

First, let's see the hierarchical model in action. Click the badge below to open the full interactive notebook:

We'll generate synthetic insurance claim data for three policy types with deliberately unbalanced sample sizes, then fit a hierarchical model.

import numpy as np
import pandas as pd
import pymc as pm
import arviz as az
import matplotlib.pyplot as plt

np.random.seed(42)

# Three policy types: lots of Auto data, moderate Home, very little Commercial
groups = {
    'Auto':       {'n': 500, 'intercept': 7.5, 'slope': 0.30},
    'Home':       {'n': 300, 'intercept': 8.2, 'slope': 0.50},
    'Commercial': {'n':  50, 'intercept': 9.0, 'slope': 0.70},
}

records = []
for i, (name, p) in enumerate(groups.items()):
    x = np.random.normal(12, 1.5, p['n'])  # log property value (~$160k median)
    noise = np.random.normal(0, 0.8, p['n'])
    y = p['intercept'] + p['slope'] * (x - 12) + noise
    for j in range(p['n']):
        records.append({
            'policy_type': name, 'group_idx': i,
            'log_property_value': x[j], 'log_claim_severity': y[j],
        })

df = pd.DataFrame(records)

Each policy type has a different intercept and slope, but Commercial has just 50 data points. Now let's fit the hierarchical model:

n_types = len(groups)
idx = df['group_idx'].values
x_centered = df['log_property_value'].values - 12  # center the predictor

with pm.Model() as hierarchical_model:
    # Hyperpriors: the "population" distribution that groups are drawn from
    mu_alpha = pm.Normal('mu_alpha', mu=8, sigma=2)
    sigma_alpha = pm.HalfNormal('sigma_alpha', sigma=2)
    mu_beta = pm.Normal('mu_beta', mu=0, sigma=2)
    sigma_beta = pm.HalfNormal('sigma_beta', sigma=2)

    # Group-level parameters, drawn from the population
    alpha = pm.Normal('alpha', mu=mu_alpha, sigma=sigma_alpha, shape=n_types)
    beta = pm.Normal('beta', mu=mu_beta, sigma=sigma_beta, shape=n_types)

    # Observation noise
    sigma = pm.HalfNormal('sigma', sigma=2)

    # Linear model
    mu = alpha[idx] + beta[idx] * x_centered
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma,
                      observed=df['log_claim_severity'].values)

    # Sample the posterior
    hierarchical_trace = pm.sample(
        1000, tune=1500, cores=4, chains=4,
        random_seed=42, target_accept=0.9,
    )

# Summarise the results
print(az.summary(hierarchical_trace, var_names=['alpha', 'beta', 'sigma']))

You just estimated three group-specific regression lines (one per policy type) while letting them share statistical strength through a common population distribution. The Commercial group, despite having only 50 claims, gets a stable estimate because it borrows information from Auto and Home.

What Just Happened?

The Three Pooling Strategies

To understand why the hierarchical model is special, let's compare it against the two extreme alternatives.

Complete pooling ignores group differences entirely. One intercept, one slope for all 850 data points:

with pm.Model() as pooled_model:
    alpha = pm.Normal('alpha', mu=8, sigma=5)
    beta = pm.Normal('beta', mu=0, sigma=5)
    sigma = pm.HalfNormal('sigma', sigma=2)

    mu = alpha + beta * x_centered
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma,
                      observed=df['log_claim_severity'].values)

    pooled_trace = pm.sample(1000, tune=1000, cores=4, chains=4, random_seed=42)

No pooling treats each group as completely independent. Three separate intercepts, three separate slopes, with no shared information:

with pm.Model() as unpooled_model:
    alpha = pm.Normal('alpha', mu=8, sigma=5, shape=n_types)
    beta = pm.Normal('beta', mu=0, sigma=5, shape=n_types)
    sigma = pm.HalfNormal('sigma', sigma=2)

    mu = alpha[idx] + beta[idx] * x_centered
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma,
                      observed=df['log_claim_severity'].values)

    unpooled_trace = pm.sample(1000, tune=1000, cores=4, chains=4, random_seed=42)

The comparison reveals the key insight. Complete pooling gives a single line dominated by Auto and Home (which together make up 94% of the data), systematically underestimating Commercial's higher intercept and steeper slope. No pooling gives each group its own line, but Commercial's estimate is noisy because it only has 50 points. Partial pooling (the hierarchical model) sits between the two: each group gets its own line, but the lines are gently pulled toward the population average. Groups with little data get pulled more.

How Hyperpriors Create Partial Pooling

The magic ingredient is the hyperpriors: mu_alpha, sigma_alpha, mu_beta, sigma_beta. These define a "population distribution" from which group-level parameters are drawn.

Think of $\mu_\alpha$ as the average intercept across all policy types, and $\sigma_\alpha$ as how much the types are allowed to differ. If the data supports large differences, $\sigma_\alpha$ will be large and each group behaves almost independently (like no pooling). If the groups are similar, $\sigma_\alpha$ shrinks and the group estimates collapse toward the population mean (like complete pooling).

The sampler learns $\sigma_\alpha$ from the data itself. You don't have to choose between pooling and no pooling; the model figures out the right amount of sharing automatically.

Shrinkage: The Key Insight

Shrinkage is the defining feature of hierarchical models. Compare each group's raw sample mean (what you'd get from no pooling) to its hierarchical posterior mean:

Commercial's intercept gets pulled the most toward the population mean, because it has the least data and therefore the most uncertainty. Auto barely moves, because 500 data points leave little room for the prior to override the evidence. This is exactly the Bayesian compromise between prior and data that we explored in From MLE to Bayesian Inference.

MCMC Diagnostics with ArviZ

Before trusting the results, we need to verify the sampler converged. ArviZ provides the standard toolkit:

az.plot_trace(hierarchical_trace, var_names=['alpha', 'beta', 'sigma'])
plt.tight_layout()

Three things to check:

Trace mixing: The chains should look like "hairy caterpillars", bouncing randomly around a stable mean. If a chain gets stuck or drifts, something is wrong.
R-hat (the Gelman-Rubin statistic): Should be below 1.01 for every parameter. Values above 1.1 indicate the chains haven't converged to the same distribution.
Effective sample size (ESS): Should be at least 400 per chain. Low ESS means the samples are highly autocorrelated and the posterior estimates are unreliable.

summary = az.summary(hierarchical_trace, var_names=['alpha', 'beta', 'sigma'])
print(summary[['mean', 'sd', 'hdi_3%', 'hdi_97%', 'r_hat', 'ess_bulk']])

If you've worked through our MCMC Metropolis-Hastings tutorial, you'll recognise the core idea: the sampler explores the posterior by proposing moves and accepting or rejecting them. PyMC uses the NUTS sampler (No U-Turn Sampler), a sophisticated variant of Hamiltonian Monte Carlo that automatically tunes step sizes and trajectory lengths.

Going Deeper

Why Not a Normal Likelihood?

The model above uses a Normal likelihood, which assumes claim amounts are symmetric around the mean. In practice, insurance claims are heavy-tailed: most claims are small, but a few are enormous. The original code I adapted for this tutorial used a Laplace likelihood to handle this:

The Laplace distribution has heavier tails than the Normal and is more robust to outliers. In PyMC, swapping the likelihood is a single line change:

# Replace:  pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)
# With:     pm.Laplace('y_obs', mu=mu, b=b, observed=y)

Modelling the Spread Too: Heteroscedastic Regression

The original code goes further. It models both the location $\mu$ and the scale $b$ of the Laplace distribution as functions of the covariates. This is heteroscedastic regression: the amount of noise varies across observations.

The $\exp$ ensures both $\mu$ and $b$ are positive (claim severity can't be negative). Each $\beta$ and $\gamma$ coefficient gets its own hierarchical structure:

with pm.Model() as full_model:
    n_groups = 3  # policy types

    # Hyperpriors for intercept
    beta0_mu = pm.Normal('beta0_mu', mu=0, sigma=1)
    beta0_sig = pm.InverseGamma('beta0_sig', alpha=2, beta=5)

    # Group-level intercepts
    beta0 = pm.Normal('beta0', mu=beta0_mu, sigma=beta0_sig, shape=n_groups)

    # ... repeat for each coefficient and for gamma (scale) parameters ...

    mu = pm.math.exp(beta0[group_idx] + beta1[group_idx] * pm.math.log(X[:, 0]) + ...)
    b  = pm.math.exp(gamma0[group_idx] + gamma1[group_idx] * pm.math.log(X[:, 0]) + ...)

    y_obs = pm.Laplace('y_obs', mu=mu, b=b, observed=y)

Notice the pm.InverseGamma hyperprior for the variance parameters. The InverseGamma is the conjugate prior for Normal variance, making it a natural choice. With alpha=2, beta=5, it places mass on moderate variance values while allowing large ones.

The Three-Tier Model

The code also contains a three-tier hierarchy. Instead of just grouping by policy type, it nests policy type within region:

Population → Policy Type → (Region × Policy Type)

At the top level, hyper-hyperpriors define the global population. At the middle level, each policy type gets its own parameters drawn from the population. At the bottom level, each (region, policy type) combination gets parameters drawn from its policy type's distribution. The group-level parameters become 2D arrays with shape (n_regions, n_policy_types):

# Hyper-hyperpriors (population level)
beta0_mu_mu = pm.Normal('beta0_mu_mu', mu=0, sigma=1)
beta0_mu_sig = pm.InverseGamma('beta0_mu_sig', alpha=2, beta=5)

# Hyperpriors (policy type level)
beta0_mu = pm.Normal('beta0_mu', mu=beta0_mu_mu, sigma=beta0_mu_sig, shape=n_types)
beta0_sig = pm.InverseGamma('beta0_sig', alpha=2, beta=5, shape=n_types)

# Priors (region × policy type level)
beta0 = pm.Normal('beta0', mu=beta0_mu, sigma=beta0_sig, shape=(n_regions, n_types))

This allows a commercial policy in an urban area to differ from one in a suburban area, while both borrow strength from the overall commercial distribution, which itself borrows from the global population.

When Not to Use Hierarchical Models

Hierarchical models aren't always necessary. If every group has plenty of data (thousands of observations), no pooling gives nearly identical results to partial pooling because the data overwhelms the prior. The hierarchy adds complexity and sampling time for little benefit.

They can also struggle with very few groups. With only 2 groups, the hyperprior variance $\sigma_\alpha$ is estimated from just 2 data points (the two group-level parameters), making it unreliable. Most practitioners suggest hierarchical models shine with 5 or more groups, though the exact threshold depends on within-group sample sizes.

Where This Comes From

Lindley and Smith (1972)

The mathematical foundation was laid by Dennis Lindley and Adrian Smith in their 1972 paper "Bayes Estimates for the Linear Model." They formalised the multi-level Normal model:

The key result: the posterior mean of $\boldsymbol{\theta}$ is a matrix-weighted average of the group-specific MLE and the prior mean. Groups with more data (higher precision in $C^{-1}$ ) weight their own MLE more heavily; groups with less data lean more on the prior. This is the formal statement of shrinkage.

Efron and Morris (1977): The James-Stein Connection

The frequentist justification for shrinkage came from an unexpected direction. In 1977, Brad Efron and Carl Morris showed that the James-Stein estimator (which shrinks group means toward the grand mean) dominates the usual sample means in terms of total squared error, for three or more groups simultaneously. This was a shocking result: even if the groups have nothing in common, shrinking toward their average reduces total estimation error.

"The James-Stein estimator achieves a smaller total mean squared error than the individual sample means, for any configuration of the true means, provided there are three or more groups."

The hierarchical Bayesian model produces estimates that are closely related to the James-Stein estimator. The Bayesian framework provides a natural explanation: when data is scarce, it's rational to hedge toward the population average rather than fully committing to a noisy local estimate.

Gelman and Hill (2006)

The practical handbook for hierarchical models is Andrew Gelman and Jennifer Hill's Data Analysis Using Regression and Multilevel/Hierarchical Models. Chapter 12 presents the exact three-model comparison we built above (complete pooling, no pooling, partial pooling) using radon measurements across US counties. Their formulation uses the non-centred parameterisation:

This reparameterisation often improves MCMC sampling efficiency because the sampler explores a standard Normal geometry rather than a funnel-shaped one. PyMC can apply this transformation automatically, but it's worth knowing when your model has divergences.

Gelman et al.'s Bayesian Data Analysis (3rd edition, 2013) provides the full mathematical treatment in Chapter 5, including the relationship between hierarchical Bayes, empirical Bayes, and the James-Stein estimator.

Interactive Tools

Bayes' Theorem Calculator — Explore Bayesian updating interactively before diving into hierarchical models
A/B Test Calculator — See Bayesian hypothesis testing in action, a common application of hierarchical models

From MLE to Bayesian Inference: The conceptual foundation for priors, posteriors, and why Bayesian estimates beat point estimates.
MCMC Island Hopping: Understanding Metropolis-Hastings: How the sampler that powers PyMC actually explores the posterior distribution.
Linear Regression: Five Ways: The non-hierarchical regression baseline that this post extends with group structure.

Frequently Asked Questions

What is hierarchical Bayesian regression?

Hierarchical (or multilevel) regression models data that is naturally grouped (students within schools, patients within hospitals) by allowing parameters to vary across groups while sharing a common prior distribution. This "partial pooling" approach borrows strength across groups, producing better estimates for small groups than fitting each group independently.

What is the difference between complete pooling, no pooling, and partial pooling?

Complete pooling ignores group differences entirely (one model for all). No pooling fits a separate model per group (no information sharing). Partial pooling (hierarchical) sits in between: each group gets its own parameters, but they are pulled towards a shared distribution. This is especially valuable when some groups have very few observations.

Why use PyMC for hierarchical models?

PyMC uses MCMC sampling to handle the complex posterior distributions that hierarchical models produce. It naturally propagates uncertainty through all levels of the hierarchy. Frequentist alternatives (like lme4 in R) can fit similar models but do not provide the same rich uncertainty quantification or flexibility for custom model structures.

How do I diagnose convergence in PyMC?

Check the trace plots for good mixing (no trends, no stuck chains), verify that R-hat values are close to 1.0 (below 1.01), and ensure effective sample sizes are sufficiently large (at least 400 per chain). Divergent transitions indicate the sampler is struggling with the posterior geometry and may require reparameterisation.

When should I use a hierarchical model instead of a standard regression?

Use hierarchical models whenever your data has a natural grouping structure and you want to make inferences about individual groups. They are especially valuable when group sizes are unequal: small groups benefit from borrowing strength, and large groups are barely affected by the pooling. If all groups have abundant data, the results will be similar to fitting separate models.