Bayesian
hierarchical models

01 / From one person to a population

Similar people, different parameters

A hierarchical model learns what individuals have in common, while allowing them to differ.

In the Bayesian and variational inference introduction, we fitted an unknown parameter using a prior and uncertain data. Now imagine fitting the same model to several people. Each person has their own parameters, but those parameters may have a shared structure.

For example, a decaying response might have a different amplitude and decay rate in each person. A neural model might have different connection strengths. We want to estimate those individual differences and the population they come from.

There are three familiar ways to approach this:

Partial pooling is not an instruction to make everyone identical. It is the consequence of a model in which individual parameters are related. How much they move towards a shared mean depends on the evidence and on the inferred amount of real variation.

02 / Write down the hierarchy

A prior can have its own parameters

Start with one noisy estimate per person. Let $y_i$ be person $i$’s observed estimate, $s_i$ its known measurement standard error, and $\theta_i$ the underlying quantity we want to infer:

The next level says that individual parameters come from a population with mean $\mu$ and standard deviation $\tau$:

We can also be uncertain about the population. Put hyperpriors on its mean and spread. For the browser examples we use:

The half-normal is a normal distribution restricted to positive values. Its argument is a scale in the same units as $\theta$, not a variance or a precision. These are illustrative priors for this synthetic problem; meaningful priors depend on the parameter’s units and plausible range.

The joint posterior is proportional to the likelihoods, the population distributions, and the hyperpriors:

The hierarchy is statistical. Its levels could be trials within sessions, sessions within people, or people within sites. This does not, by itself, imply a hierarchy of brain areas.

03 / The familiar Bayesian update

Partial pooling is precision weighting

For a moment, suppose the population mean and spread are known. Each person then has a Gaussian prior, $\mathcal N(\mu,\tau^2)$. We can reuse the exact update from the companion page:

An equivalent expression makes the shrinkage easy to see:

$w_i$ is the fraction of the individual-to-population difference retained. Reliable measurements have a larger $w_i$. A smaller population spread gives stronger pooling. With a very large spread, the population tells us little about any particular person.

At the default values, A’s posterior mean is 13.2, while B’s is 10.8. B contributes less precise evidence, so the shared prior has more influence. Shrinkage also occurs for well-measured people; its amount is just smaller.

This experiment conditions on μ and τ. Its intervals exclude uncertainty about the population. A fully Bayesian hierarchy learns those quantities and propagates their uncertainty.

04 / Learn the shared distribution

Who sets the population mean and spread?

The group informs each person, and each person informs the group.

In the full hierarchy, $\mu$ and $\tau$ are unknown. A reliable individual estimate provides more information about the population. If reliable estimates differ substantially, the model supports more between-person variation. If much of the observed spread can be explained by measurement error, less genuine variation may be needed.

For this small Gaussian problem, we can integrate out each person’s parameter:

This expression separates two reasons observed estimates can differ: real variation, $\tau^2$, and measurement variance, $s_i^2$. The same observed differences can support different population conclusions when their reliability changes.

Try Reliable differences, then Noisy differences. The observed values remain the same, while their standard errors change. Also reduce the half-normal scale: with only six people, the prior on variation can influence the answer. It is part of the model, rather than a tuning parameter with no scientific meaning.

05 / Connect the model to the algorithm

Hierarchical describes the model. Variational describes inference.

A hierarchical model specifies how observations, individual parameters, and population parameters depend on each other. Variational inference is one way to approximate its joint posterior. The same hierarchy could instead be analysed using quadrature, posterior sampling, or another suitable method.

For a larger problem, collect the unknowns into $z=(\theta_{1:N},\mu,\tau,\ldots)$ and choose an approximate distribution $q(z)$. Fit it by maximising the evidence lower bound:

The first term rewards distributions that put mass where the joint model fits well. The second is entropy: the approximation is fitting a distribution, not just choosing a best point. Maximising this bound is equivalent to minimising $\operatorname{KL}(q(z)\,\|\,p(z\mid y))$.

A simple mean-field family might separate individual factors from group and precision factors. It is computationally convenient, but it removes some posterior dependencies. Population and individual estimates are coupled in the true posterior, so this approximation can understate marginal uncertainty.

The same precision-weighted errors, at two levels

Replace the direct observation by a forward model $f_i(\theta_i)$ and a population regression with design row $x_i^\top$:

There is a data residual $e_i=y_i-f_i(\theta_i)$, and a population residual $r_i=\theta_i-B^\top x_i$. Holding the covariance matrices fixed, the quadratic part of the negative log posterior is

Here $\Pi_{y_i}=\Sigma_{y_i}^{-1}$ and $\Pi_b=\Sigma_b^{-1}$. A MAP-style gradient step for an individual becomes

The first term improves that person’s data prediction; the second improves agreement with the population model. Covariance or precision learning also needs the Gaussian log-determinant terms and hyperpriors. Minimising squared errors alone is not enough.

This step updates a point estimate. A Gaussian approximation based on local curvature adds uncertainty, but MAP plus an inverse Hessian is not automatically a complete variational Laplace algorithm. Variational Laplace also uses a variational objective and the appropriate updates for its model.

For linear Gaussian models, conjugate variational updates are often available. For nonlinear neural or behavioural models, local Gaussian approximations can be useful, while multimodality and strong nonlinearities may require richer approximations or sampling. A successful optimiser does not establish that an approximate posterior is well calibrated.

06 / Ask what you are predicting

The average person is not a new person

A precise population mean can coexist with substantial individual variation.

Estimating $\mu$ asks where the population is centred. Predicting a new person’s latent parameter $\theta_{\mathrm{new}}$ also includes real between-person variation. Predicting their measurement $y_{\mathrm{new}}$ adds observation noise:

These identities apply to this scalar model with conditionally independent, zero-mean new-person deviations and known future measurement variance. More people can sharpen the group mean without making the population itself homogeneous.

For someone whose data we have already observed, use their individual posterior instead. Also consider who a “new person” is: a new site, a different clinical population, or a changed task may need a different model.

07 / Neural models and group effects

DCM + PEB is already a hierarchical model

In Dynamic Causal Modelling, each person’s generative model relates neural parameters to measured data. Parametric Empirical Bayes (PEB) then models those parameters at the group level, including group means and effects of covariates. Importantly, it uses uncertainty from the individual models, rather than treating their fitted parameters as perfectly known observations.

The population mean need not be the same for everybody. With an intercept, a group indicator, and a centred age covariate, for example,

The hierarchy pools residual differences around the appropriate predicted mean. It does not require shrinking a patient and a control towards one identical grand mean. The design matrix determines which systematic differences the model can represent.

The SPM PEB documentation explains the DCM-to-group workflow, Bayesian model reduction, and prediction. The browser’s scalar example illustrates a statistical principle; it is not a replacement for SPM’s DCM/PEB implementation.

A hierarchy assumes that units are suitably comparable after accounting for included predictors. Missing covariates, site effects, heavy tails, or unmodelled subgroups can make pooling misleading. Partial pooling improves estimation under useful assumptions; it does not guarantee improvement for every individual or dataset.

08 / Build, check, and extend

A reusable MATLAB framework

Start by separating the forward model from the inference machinery. The forward function produces predictions for one person. The inference function handles individual distributions, the population regression, precision updates, and convergence.

The accompanying framework uses this interface:

fit = hbayes_fit(model, subjects, X, priors, opts);

subjects contains each person’s data and inputs. X is the population design matrix. priors defines the group and precision priors. model.forward can describe a constant, a regression, or a nonlinear response.

Unzip, open the matlab_hierarchical_bayes folder in MATLAB, and run:

addpath(genpath(pwd))
report = run_tests();
examples = run_examples(true);

For a small Gaussian summary-data problem, the model definition is just:

y = [14, 14, 7, 9, 12, 11];
s = [1, 4, 1, 2, 0.5, 3];
subjects = cell(numel(y), 1);
for i = 1:numel(y)
    subjects{i} = struct('y', y(i), 'covariance', s(i)^2);
end

model.forward = @(theta, subject) theta(1) + zeros(size(subject.y));
model.jacobian = @(theta, subject) ones(numel(subject.y), 1);
model.isLinear = true;

X = ones(numel(y), 1);             % Population intercept
priors.betaMean = 10;
priors.betaCov = 25;
priors.betweenShape = 2;           % Gamma prior on PRECISION
priors.betweenRate = 4;            % Shape/rate convention
opts.fixedNoisePrecision = 1;     % Covariances s(i)^2 are known

fit = hbayes_fit(model, subjects, X, priors, opts);
disp(fit.betaMean)
disp(fit.thetaMean)
disp(fit.objectiveKind)
assert(fit.converged, 'Inspect the convergence diagnostics.');

Four problems to try

For the decay model, use $\theta_i=(\log A_i,\log k_i)$ and define the forward function as

This gives positive physical parameters while placing the hierarchy explicitly on log parameters. A difference on the log scale becomes a multiplicative difference on the physical scale.

Check the model, as well as the fit

Before fitting, simulate from the prior: do the generated parameters and responses make sense? After fitting, simulate replicated data from the posterior predictive distribution: does the model reproduce the features that matter? Compare against analytic results in simple cases, and check sensitivity to priors, initialisation, and approximation choices.

For prediction across people, hold out people. For prediction across sites, hold out sites where possible. More observations from the same person do not provide the same validation as genuinely new individuals.

For the inference background, see Blei, Kucukelbir & McAuliffe (2017) on variational inference, and Friston et al. (2007) on variational Laplace. All experiments use synthetic data. Equations and explanations work offline; the interactive calculations run in your browser.