Hamiltonian Monte Carlo for Bayesian GLMMs

Why Hamiltonian Monte Carlo?

A Bayesian GLMM can have many correlated continuous unknowns: fixed effects, group-level effects, and variance or correlation parameters. A random-walk Metropolis proposal changes these parameters without using the shape of the posterior. In contrast, Hamiltonian Monte Carlo (HMC) uses the posterior gradient to construct a trajectory that tends to move through regions of relatively high posterior density. The trajectory is only a proposal mechanism: an acceptance step corrects for its numerical approximation so the desired posterior remains the stationary distribution.

HMC is designed for a differentiable density on continuous parameters. The observed responses may be binary or counts; what must be continuous and differentiable is the parameter vector being sampled and its log posterior. HMC does not directly update discrete latent variables. The earlier material on Pólya-Gamma sampling describes a different strategy for logistic random-effects models, using augmentation to obtain convenient conditional updates.

The Posterior Target and Its Gradient

For a GLMM, let \(y_{ij}\) be the response for observation \(i\) in group \(j\), \(x_{ij}\) and \(z_{ij}\) be fixed- and random-effect design vectors, \(\beta\) be the fixed-effect vector, and \(b_j\) be the group-level effect. Conditional on \(b_j\), a standard GLMM specifies a response distribution from an exponential family and a linear predictor such as \(\eta_{ij}=x_{ij}^{\mathsf T}\beta+z_{ij}^{\mathsf T}b_j\). The random effects have a distribution governed by covariance parameters, and prior distributions are assigned to the unknowns.

Write \(\theta\) for the continuous unknowns to be sampled, including any random effects and covariance parameters represented on an unconstrained scale. The posterior density, up to a normalizing constant, is the likelihood times the prior density. HMC uses its negative log as a potential energy.

Definition: For a differentiable posterior density \(\pi(\theta\mid y)\), define the potential energy \(U(\theta)=-\log\pi(\theta\mid y)+C\), where \(C\) is any constant independent of \(\theta\). The HMC force is the negative gradient \(-\nabla U(\theta)=\nabla\log\pi(\theta\mid y)\). The constant can be omitted because it does not change the gradient or acceptance probability.

A gradient combines the contributions from the likelihood and prior. For a parameter \(\theta_k\), \(\partial\log\pi(\theta\mid y)/\partial\theta_k\) is the sum of the derivatives of the log likelihood and log prior (and any transformation Jacobian, if a constrained parameter has been transformed). This is the calculation that makes HMC specific to the model: gradients must be derived for the actual GLMM target.

Worked Example: Logistic Random-Intercept Posterior Gradient

Consider a hypothetical model with two groups and three binary observations per group. In group 1, the responses are \(1,0,1\); in group 2, they are \(0,0,1\). Let \(p_j=\operatorname{logit}^{-1}(\beta+b_j)\), with \(\beta\sim N(0,4)\) and independent \(b_j\sim N(0,1)\). Evaluate the log-posterior gradient at \(\beta=0.2,\ b_1=0.3,\ b_2=-0.1\). The group effects are included as unknowns, and the variance parameters are treated as fixed for this calculation.

For a Bernoulli-logit observation, the derivative of its log likelihood with respect to its linear predictor is \(y_{ij}-p_j\). Thus the likelihood contribution for group \(j\) is the number of successes minus \(3p_j\). The normal-prior contributions to the gradients are \(-\beta/4\) for \(\beta\) and \(-b_j\) for \(b_j\).

The fitted probabilities at the stated point are \(p_1=\operatorname{logit}^{-1}(0.2+0.3)=0.622459\) and \(p_2=\operatorname{logit}^{-1}(0.2-0.1)=0.524979\). There are three successes overall, two in group 1, and one in group 2.

Therefore,

\frac{\partial\log\pi}{\partial\beta} =(3-3p_1-3p_2)-\frac{0.2}{4} =3-1.867377-1.574937-0.05 =-0.492314.

For each group effect, the fixed effect \(\beta\) does not enter its prior penalty, so

\frac{\partial\log\pi}{\partial b_1} =2-3(0.622459)-0.3=-0.167377, \qquad \frac{\partial\log\pi}{\partial b_2} =1-3(0.524979)-(-0.1)=-0.474937.

The gradient of the log posterior is approximately \((-0.4923,-0.1674,-0.4749)^{\mathsf T}\), so the gradient of \(U\) has the opposite signs. These values describe the local slope at one parameter point; they are not themselves posterior estimates or standard errors.

Momentum and Simulated Hamiltonian Trajectories

HMC adds an auxiliary momentum vector \(p\) to the parameter vector \(\theta\). Choose a symmetric positive-definite mass matrix \(M\), and draw \(p\sim N(0,M)\). The joint target for \((\theta,p)\) is proportional to \(\exp[-H(\theta,p)]\), where the Hamiltonian is the sum of potential and kinetic energy:

H(\theta,p)=U(\theta)+K(p), \qquad K(p)=\frac{1}{2}p^{\mathsf T}M^{-1}p.

The momentum is auxiliary, not a GLMM parameter of scientific interest. Its role is to provide direction and persistence to the movement through parameter space. In exact Hamiltonian dynamics, position and momentum evolve according to

\frac{d\theta}{dt}=M^{-1}p, \qquad \frac{dp}{dt}=-\nabla U(\theta).

The gradient of the posterior therefore affects how momentum changes along a trajectory. For a differentiable \(U\), these equations conserve energy: by the chain rule,

\frac{dH}{dt} =\nabla U(\theta)^{\mathsf T}M^{-1}p +p^{\mathsf T}M^{-1}[-\nabla U(\theta)] =0.

Exact trajectories are generally unavailable for GLMM posteriors, so HMC approximates them with the leapfrog integrator. Given step size \(\epsilon\), one leapfrog step updates momentum by a half-step, position by a full step, and momentum by another half-step:

p_{t+\epsilon/2}=p_t-\frac{\epsilon}{2}\nabla U(\theta_t), \qquad \theta_{t+\epsilon}=\theta_t+\epsilon M^{-1}p_{t+\epsilon/2}, \qquad p_{t+\epsilon}=p_{t+\epsilon/2}-\frac{\epsilon}{2}\nabla U(\theta_{t+\epsilon}).

Repeating this update \(L\) times produces a proposed endpoint. Leapfrog is useful because it is reversible (with momentum reversal) and volume-preserving; however, its finite-step path generally has a small energy error. The Metropolis step accounts for that error.

Worked Example: A Two-Step Leapfrog Trajectory

To see the arithmetic without GLMM notation obscuring it, consider a one-dimensional target with \(U(\theta)=\theta^2/2\), mass \(M=1\), starting point \(\theta_0=1\), and sampled momentum \(p_0=0.5\). Use step size \(\epsilon=0.1\) and \(L=2\) steps. Here \(\nabla U(\theta)=\theta\).

For the first step, the initial half momentum update is \(p_{1/2}=0.5-(0.1/2)(1)=0.45\). The position update gives \(\theta_1=1+0.1(0.45)=1.045\). The final half update gives \(p_1=0.45-(0.1/2)(1.045)=0.39775\).

For the second step, \(p_{3/2}=0.39775-(0.05)(1.045)=0.34550\), and \(\theta_2=1.045+0.1(0.34550)=1.07955\). The last half update gives \(p_2=0.34550-(0.05)(1.07955)=0.2915225\).

The initial energy is \(H_0=1^2/2+0.5^2/2=0.625\). At the endpoint, \(U(\theta_2)=1.07955^2/2=0.582714\), and \(K(p_2)=0.2915225^2/2=0.0424927\), so \(H_2\approx0.625207\). The energy error is \(\Delta H=H_2-H_0\approx0.000207\). The acceptance probability is

a=\min\{1,\exp(-\Delta H)\} =\min\{1,\exp(-0.000207)\} \approx0.999793.

Thus this proposed endpoint is accepted with probability about \(0.9998\). An actual accept-or-reject decision also requires a uniform random draw; the high probability does not mean acceptance is guaranteed.

The Acceptance Step and Constrained Parameters

After simulating the trajectory, HMC proposes its endpoint and uses the Hamiltonian difference to correct for discretization error. If the starting and proposed states are \((\theta,p)\) and \((\theta^\ast,p^\ast)\), respectively, accept the proposal with probability

a=\min\{1,\exp[-H(\theta^\ast,p^\ast)+H(\theta,p)]\}.

If the proposal is rejected, retain the current parameter value. Momentum is normally reversed in the proposal construction so the transition is reversible. Together, the reversible, volume-preserving leapfrog proposal and this Metropolis correction give the intended joint target, whose marginal distribution for \(\theta\) is the posterior. A trajectory is therefore not accepted merely because it follows the gradient; the correction is essential.

Some GLMM parameters have constrained support, such as positive standard deviations or covariance matrices. HMC commonly operates on unconstrained coordinates instead. For example, if \(\alpha=\log\sigma\), then \(\sigma=\exp(\alpha)\). A density specified for \(\sigma\) must be transformed to a density for \(\alpha\):

\pi_\alpha(\alpha\mid y) =\pi_\sigma(\exp\alpha\mid y)\left|\frac{d\sigma}{d\alpha}\right| =\pi_\sigma(\exp\alpha\mid y)\exp(\alpha).

The log density in the new coordinate includes the Jacobian term \(\alpha\); omitting it changes the target distribution. The gradient used by HMC must be the derivative of the transformed log posterior, not simply the old derivative with the parameter renamed.

Worked Example: The Jacobian for a Positive Scale

Suppose a hypothetical random-intercept model has two group effects \(b_1=0.5\) and \(b_2=-1\), conditional on a positive scale \(\sigma\), with \(b_j\mid\sigma\sim N(0,\sigma^2)\). Give \(\sigma\) a half-normal prior with scale 2, whose density is proportional to \(\exp(-\sigma^2/8)\) for \(\sigma>0\). Evaluate the derivative of the transformed log density with respect to \(\alpha=\log\sigma\) at \(\sigma=1\), including these prior contributions.

The two conditional normal densities contribute \(-2\log\sigma-\sum_j b_j^2/(2\sigma^2)\) to the log density, apart from constants. Differentiating with respect to \(\alpha\), using \(d\sigma/d\alpha=\sigma\), gives \(-2+\sum_j b_j^2/\sigma^2\). At the stated values, \(\sum_j b_j^2=0.25+1=1.25\), so this contribution is \(-2+1.25=-0.75\).

The half-normal prior contributes \(-\sigma^2/8\) on the original scale. On the \(\alpha\) scale its log density also includes the Jacobian term \(+\alpha\). Its derivative is therefore \(-\sigma^2/4+1\), which at \(\sigma=1\) equals \(0.75\). Adding the two contributions gives a total derivative of \(-0.75+0.75=0\). This is only the derivative from the stated prior and random-effect density; a data likelihood depending on \(\sigma\) would contribute its own derivative.

Practical Checks, Alternatives, and Limitations

HMC requires a correct gradient and a numerically stable trajectory. Before interpreting output, check the model and implementation. In particular, compare analytic or automatically differentiated gradients with finite-difference approximations at several parameter values; verify that transformed parameters include their Jacobians; and confirm that the log posterior and gradient are finite at plausible starting points. If a step size is too large, the numerical trajectory can accumulate substantial energy error, leading to low acceptance or unstable proposals. If it is extremely small, proposals can be computationally expensive and make little progress per gradient evaluation.

A further practical issue is that posterior coordinates can have very different scales or strong dependence. The mass matrix influences how momentum translates into movement in each direction; a poor choice can make trajectories inefficient. HMC also requires repeated gradient evaluations, which may be costly for a large GLMM. It does not remove the need to assess whether the resulting chain explores the posterior adequately. Detailed warmup choices and specialized HMC diagnostics require separate treatment.

HMC differs from common alternatives in what it uses to propose updates. Random-walk Metropolis–Hastings does not require gradients and can handle targets that are not differentiable, but its local proposals may move slowly in a high-dimensional posterior. Gibbs sampling is attractive when full conditional distributions can be sampled directly, but such updates are not generally available for all GLMM parameters. Pólya-Gamma augmentation provides a useful route for certain logistic models, while HMC instead follows gradients of a continuous target. None is universally superior: model structure, dimension, differentiability, computation, and empirical chain behavior all matter.

Common pitfall: “The trajectory follows the posterior gradient, so its endpoint is automatically a valid draw.” This is incorrect. A finite-step leapfrog trajectory only approximates exact energy-conserving dynamics. The Metropolis acceptance probability based on the Hamiltonian difference is what corrects that numerical approximation; omitting the correction generally changes the target distribution.

HMC is most directly applicable when the posterior can be expressed as a smooth density over continuous coordinates and gradients can be evaluated reliably. It is not a universal solution to weakly identified variance components, difficult posterior geometry, or discrete latent-state sampling. In a GLMM, the method’s benefit is a gradient-guided proposal for the continuous parameters, not a guarantee of rapid or reliable exploration.

Check Your Understanding

Answer each question, then compare with the solutions.

Why does HMC use the negative log posterior as potential energy, and what role does its gradient play?
If a leapfrog proposal has \(H_{\mathrm{new}}-H_{\mathrm{old}}=0.4\), what is its acceptance probability to three significant figures?
For \(\alpha=\log\sigma\), what Jacobian term must be added to the log density when transforming a density from positive \(\sigma\) to \(\alpha\)?
Name two conditions or checks that matter for applying HMC to a GLMM.

Solutions. (1) The negative log posterior gives a potential whose exponential recovers the posterior density up to a constant. Its gradient determines the force, \(-\nabla U=\nabla\log\pi\), that changes momentum along the trajectory. (2) \(a=\min(1,\exp(-0.4))\approx0.670\). (3) Since \(d\sigma/d\alpha=\exp(\alpha)\), the log-Jacobian is \(+\alpha\). (4) Examples include having a differentiable continuous target, computing correct gradients, including transformation Jacobians, and checking that log densities and gradients are finite at plausible starting points.