Bayesian Hierarchical Linear Regression

Model

A supermarket chain wants to know how far its cheese sales fall when it raises the price. It could put that question to each of its 88 stores separately — but no single store changes its price often enough to answer it, and fitted store by store the estimates come back wildly scattered, several of them with the wrong sign. It could instead pool all 88 into a single number, at the cost of forcing a shop in a wealthy suburb and one in a small town to share one answer. The hierarchical model is the middle course, and this page is its plainest form: every store keeps its own coefficients, but those coefficients are treated as draws from a common population whose shape is estimated alongside them. Each store's estimate is then informed by the other 87 — heavily where its own data is thin, barely at all where it is rich.

Unit-level regression

Let $i = 1,\ldots,n_{\text{reg}}$ index units (e.g. stores, individuals). Each unit has $n_i$ observations and its own coefficient vector $\boldsymbol{\beta}_i$:

• $y_i$ — $n_i \times 1$ vector of outcomes for unit $i$
• $\mathbf{X}_i$ — $n_i \times p$ design matrix for unit $i$
• $\boldsymbol{\beta}_i$ — $p \times 1$ unit-specific coefficients
• $\sigma_i^2$ — unit-specific error variance

Population (hierarchy)

The key feature of the hierarchical model is that unit-level coefficients are drawn from a common population distribution governed by $\boldsymbol{\Delta}$ and $\mathbf{V}_\beta$:

where $\mathbf{z}_i$ is a vector of $n_z$ observed unit characteristics (always including a constant). $\boldsymbol{\Delta}$ is an $n_z \times p$ matrix of population regression coefficients, and $\mathbf{V}_\beta$ is the $p \times p$ cross-unit covariance. Units with no covariates ($n_z=1$) reduce to simple partial pooling.

Put in words: $\boldsymbol{\Delta}$ says how a store's observable characteristics — its size, its market — shift that store's entire coefficient vector, while $\mathbf{V}_\beta$ says how much stores still differ from one another after those characteristics have been accounted for. The second is the more interesting of the two. A large $\mathbf{V}_\beta$ means the population is genuinely heterogeneous, each store's own data has to carry the weight, and little should be pooled; a small one means the stores are near-copies of each other and pooling is close to free. The sampler does not have to be told which case it is in — it estimates $\mathbf{V}_\beta$ from the data and shrinks by exactly that much.

Priors

The conjugate prior structure mirrors bayesm::rhierLinearModel:

so $\boldsymbol{\Delta}$ is matrix-normal with row covariance $\mathbf{A}^{-1}$ and column covariance $\mathbf{V}_\beta$ — the natural prior for a matrix of coefficients, being a multivariate Normal on the flattened matrix whose covariance factors into one part across rows and another across columns, so dependence in the two directions can be specified separately instead of through a single unwieldy joint covariance. A scalar $\mathbf{A} = a \mathbf{I}_{n_z}$ with small $a$ is diffuse. The error variance prior $\sigma_i^2 \sim \mathcal{IG}(\nu_e/2,\, \nu_e s_i^2/2)$ is exchangeable across units — meaning the prior treats the stores as interchangeable, encoding no belief about which particular store is which, only that they are draws from a shared population. Exchangeability is the assumption that licenses pooling in the first place.

Gibbs Sampler

Initialise $\{\boldsymbol{\beta}_i\}, \boldsymbol{\Delta}, \mathbf{V}_\beta, \{\sigma_i^2\}$. At each iteration:

Block 1 — Unit coefficients $\boldsymbol{\beta}_i$

Posterior precision combines the population prior and the unit likelihood:

Block 2 — Error variances $\sigma_i^2$

The residual sum of squares $e_i^\top e_i$ updates the IG scale:

Block 3 — Population coefficients $\boldsymbol{\Delta}$

Stacking all $n_{\text{reg}}$ units: let $\mathbf{B}$ be the $n_{\text{reg}} \times p$ matrix of current $\boldsymbol{\beta}_i$ draws. The posterior is matrix-normal:

Block 4 — Cross-unit covariance $\mathbf{V}_\beta$

Residuals from the population mean $\mathbf{R} = \mathbf{B} - \mathbf{Z}\boldsymbol{\Delta}$ update the IW scale:

Notebooks

Two notebooks, matched section by section. The Python notebook builds the Gibbs sampler from scratch (hlm_gibbs.py) and runs the full arc. Synthetic validation first (60 units × 80 observations, $k=3$, $n_z=2$): the population matrix $\Delta$ is recovered on top of its true values and the unit-level $\beta_i$ come back with a correlation of 0.999. Then the IRI cheese scanner panel (bayesm::cheese — 88 retailers, 5,555 store-weeks of supermarket checkout records), modelling $\log(\text{volume})$ on log-price and display intensity. Exploratory plots expose the core problem: each store sees only a narrow range of its own prices, so a store-specific elasticity is barely identified — and indeed the per-store OLS baseline flies apart (log-price elasticity mean $-2.27$, SD $1.77$, range $[-12.65,\ +2.86]$ — several stores come out positively sloped, which demand theory forbids). Model A adds the hierarchy with no unit covariates: partial pooling drags the noisy extremes back toward a pooled elasticity of $-2.15$ and collapses the cross-store spread. Model B then lets the population mean differ by market size (24 of 88 stores are large metros) and delivers the notebook's most interesting result — large markets are not more price-sensitive (the elasticity increment is $-0.001$, 95% CrI $[-0.46,\ 0.45]$, straddling zero) but they do sell more at any price (intercept increment $+0.78$, CrI $[0.22,\ 1.34]$). A closing comparison puts OLS, Model A and Model B side by side. The R notebook is the reference implementation, reproducing every step with bayesm::rhierLinearModel.

Downloads

Gibbs Sampler — Source Code

Download hlm_gibbs.py

"""
Bayesian Hierarchical Linear Model -- Gibbs sampler.
Replicates bayesm::rhierLinearModel (Rossi, Allenby & McCulloch 2005, Ch. 3).

Model
-----
  y_i  = X_i @ beta_i + eps_i,       eps_i ~ N(0, sigma2_i * I_ni)
  beta_i = Delta.T @ z_i + u_i,      u_i   ~ N(0, Vbeta)

Priors
------
  vec(Delta) | Vbeta ~ N(vec(Deltabar), Vbeta (x) A^{-1})
  Vbeta  ~ IW(nu, V)
  sigma2_i ~ IG(nu_e/2, nu_e*ssq_i/2)   iid, ssq_i may be unit-specific
"""

import numpy as np
from scipy.stats import invwishart


# ── private helpers ───────────────────────────────────────────────────────────

def _chol_draw(Prec, rhs, rng):
    """x ~ N(Prec^{-1} rhs, Prec^{-1}) via Cholesky of the precision matrix."""
    L          = np.linalg.cholesky(Prec)           # L L' = Prec
    L_inv_rhs  = np.linalg.solve(L,   rhs)
    mu_post    = np.linalg.solve(L.T, L_inv_rhs)    # = Prec^{-1} rhs
    return mu_post + np.linalg.solve(L.T, rng.standard_normal(len(rhs)))


def _draw_Delta(betas, Z, Vbeta, Deltabar, A, rng):
    """Delta | Vbeta, betas ~ MN(Delta_tilde, Vbeta, A_tilde^{-1})."""
    A_tilde     = A + Z.T @ Z
    A_tilde_inv = np.linalg.inv(A_tilde)
    Delta_tilde = A_tilde_inv @ (A @ Deltabar + Z.T @ betas)
    L_A = np.linalg.cholesky(A_tilde_inv)
    L_V = np.linalg.cholesky(Vbeta)
    return Delta_tilde + L_A @ rng.standard_normal(Delta_tilde.shape) @ L_V.T


def _draw_Vbeta(betas, Z, Delta, nu, V, rng):
    """Vbeta | betas, Delta ~ IW(nu + nreg, V + resid' resid)."""
    resid = betas - Z @ Delta
    return invwishart.rvs(df=nu + len(betas),
                          scale=V + resid.T @ resid,
                          random_state=rng)


# ── public API ────────────────────────────────────────────────────────────────

def rhier_linear_model(regdata, Z, Prior=None, Mcmc=None, seed=42):
    """
    Gibbs sampler for the Hierarchical Linear Model.

    Parameters
    ----------
    regdata : list of dicts  {'y': (n_i,), 'X': (n_i, k)}
    Z       : (nreg, nz) unit-level covariates  (include intercept column)
    Prior   : dict with optional keys
                Deltabar : (nz, k)   prior mean of Delta,   default 0
                A        : (nz, nz) prior precision for Delta, default 0.01*I
                nu       : int       IW df for Vbeta,        default k+3
                V        : (k, k)   IW scale for Vbeta,     default (k+2)*I
                nu_e     : float     IG df for sigma2,       default 3
                ssq      : float or (nreg,)  IG scale,      default 1
    Mcmc    : dict with keys  R, burn (default 0.2*R), keep (default 1),
                nprint (default R//10)

    Returns
    -------
    dict with
      betadraw    : (nreg, k, R_kept)
      Deltadraw   : (R_kept, nz*k)   column-major (matches R vec())
      Vbetadraw   : (R_kept, k*k)    column-major
      sigmasqdraw : (R_kept, nreg)
    """
    if Prior is None:
        Prior = {}
    if Mcmc is None:
        Mcmc = {'R': 10000}

    rng  = np.random.default_rng(seed)
    nreg = len(regdata)
    k    = regdata[0]['X'].shape[1]
    nz   = Z.shape[1]

    Deltabar = np.asarray(Prior.get('Deltabar', np.zeros((nz, k))),      float)
    A        = np.asarray(Prior.get('A',        0.01 * np.eye(nz)),      float)
    nu       = int(Prior.get('nu',   k + 3))
    V        = np.asarray(Prior.get('V',        float(k + 2) * np.eye(k)), float)
    nu_e     = float(Prior.get('nu_e', 3.0))
    ssq_raw  = Prior.get('ssq', 1.0)
    ssq      = (np.full(nreg, ssq_raw, float) if np.isscalar(ssq_raw)
                else np.asarray(ssq_raw, float))

    R      = int(Mcmc['R'])
    burn   = int(Mcmc.get('burn',   0.2 * R))
    keep   = int(Mcmc.get('keep',   1))
    nprint = int(Mcmc.get('nprint', max(R // 10, 1)))
    R_kept = (R - burn) // keep

    # precompute sufficient statistics once
    ys   = [d['y']         for d in regdata]
    Xs   = [d['X']         for d in regdata]
    XtXs = [X.T @ X        for X in Xs]
    Xtys = [X.T @ y        for X, y in zip(Xs, ys)]
    ns   = [len(y)         for y in ys]

    # storage
    betadraw    = np.zeros((nreg, k, R_kept))
    Deltadraw   = np.zeros((R_kept, nz * k))
    Vbetadraw   = np.zeros((R_kept, k * k))
    sigmasqdraw = np.zeros((R_kept, nreg))

    # initialise
    betas  = np.zeros((nreg, k))
    sigma2 = ssq.copy()
    Delta  = Deltabar.copy()
    Vbeta  = V / max(nu - k - 1, 1)

    i_kept = 0
    for r in range(R):
        Vbeta_inv = np.linalg.inv(Vbeta)
        mu_prior  = Z @ Delta              # nreg × k

        # 1. beta_i | Delta, Vbeta, sigma2_i
        for i in range(nreg):
            Prec     = Vbeta_inv + XtXs[i] / sigma2[i]
            rhs      = Vbeta_inv @ mu_prior[i] + Xtys[i] / sigma2[i]
            betas[i] = _chol_draw(Prec, rhs, rng)

        # 2. sigma2_i | beta_i  (IG draw via Gamma)
        for i in range(nreg):
            e         = ys[i] - Xs[i] @ betas[i]
            scale_i   = (nu_e * ssq[i] + float(e @ e)) / 2.0
            sigma2[i] = scale_i / rng.gamma((nu_e + ns[i]) / 2.0)

        # 3. Delta | betas, Vbeta
        Delta = _draw_Delta(betas, Z, Vbeta, Deltabar, A, rng)

        # 4. Vbeta | betas, Delta
        Vbeta = _draw_Vbeta(betas, Z, Delta, nu, V, rng)

        if nprint and (r + 1) % nprint == 0:
            print(f'  Iter {r+1:6d}/{R}')

        if r >= burn and (r - burn) % keep == 0:
            betadraw[:, :, i_kept] = betas
            Deltadraw[i_kept]      = Delta.ravel(order='F')
            Vbetadraw[i_kept]      = Vbeta.ravel(order='F')
            sigmasqdraw[i_kept]    = sigma2
            i_kept += 1

    return {
        'betadraw':    betadraw,
        'Deltadraw':   Deltadraw,
        'Vbetadraw':   Vbetadraw,
        'sigmasqdraw': sigmasqdraw,
    }


def simulate_hlm(nreg, k, Z, true_Delta, true_Vbeta, true_sigma2=1.0,
                 ni=20, seed=0):
    """Simulate data from the HLM with known parameters."""
    rng   = np.random.default_rng(seed)
    L_V   = np.linalg.cholesky(true_Vbeta)
    betas = Z @ true_Delta + rng.standard_normal((nreg, k)) @ L_V.T

    if np.isscalar(ni):
        ni = [ni] * nreg
    sig = (np.full(nreg, true_sigma2, float) if np.isscalar(true_sigma2)
           else np.asarray(true_sigma2, float))

    regdata = []
    for i in range(nreg):
        X = np.column_stack([np.ones(ni[i]),
                             rng.standard_normal((ni[i], k - 1))])
        y = X @ betas[i] + rng.normal(0.0, np.sqrt(sig[i]), ni[i])
        regdata.append({'y': y, 'X': X})

    return regdata, betas


def posterior_summary(draws, param_names=None, true_vals=None):
    """Posterior mean, SD, 95% CrI for each column of a (R, p) draws array."""
    import pandas as pd
    p = draws.shape[1]
    if param_names is None:
        param_names = [f'p{j}' for j in range(p)]
    rows = []
    for j, name in enumerate(param_names):
        d      = draws[:, j]
        lo, hi = np.percentile(d, [2.5, 97.5])
        row    = {'param': name, 'mean': d.mean(), 'sd': d.std(),
                  'q025': lo, 'q975': hi}
        if true_vals is not None:
            row['true']    = float(true_vals[j])
            row['covered'] = bool(lo <= true_vals[j] <= hi)
        rows.append(row)
    return pd.DataFrame(rows).set_index('param')

References

• Lindley, D. V. & Smith, A. F. M. (1972). Bayes estimates for the linear model. Journal of the Royal Statistical Society: Series B 34(1), 1–41. — the classic hierarchical linear model
• Rossi, P. E., Allenby, G. M. & McCulloch, R. (2005). Bayesian Statistics and Marketing. Wiley. — bayesm::rhierLinearModel, and the cheese scanner panel
• Gelman, A. & Hill, J. (2007). Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press. — the standard reference on partial pooling and varying-intercept models