Hamiltonian Monte Carlo for Joint Model Posteriors

Why Use Hamiltonian Monte Carlo?

A joint-model posterior can contain many correlated continuous quantities: longitudinal fixed effects, association parameters, and a latent random effect for each subject. A random-walk proposal may move only a short distance before its acceptance rate becomes poor. Hamiltonian Monte Carlo (HMC) uses the gradient of the log posterior to guide proposals across the posterior, potentially making larger moves while retaining a Metropolis acceptance step to correct for numerical error.

HMC is a computation method, not a different joint model. It targets the same posterior distribution as other valid MCMC updates. Its main requirements are a continuously differentiable log posterior over the parameters being updated, a way to evaluate its gradient, and a proper target distribution. Discrete quantities cannot be included in an ordinary HMC trajectory. Parameters constrained to a positive range or another restricted space require a suitable transformation or a method that handles the constraint.

A Joint-Model Target and Its Gradient

For a concrete derivation, consider a shared-random-intercept model. Subject \(i\) has longitudinal measurements \(y_{ij}\) at times \(t_{ij}\), and an observed event or censoring time \(T_i\) with event indicator \(\delta_i\). The longitudinal mean is \(\eta_i=\beta+b_i\), where \(\beta\) is a fixed intercept and \(b_i\) is a subject-specific random intercept. Suppose \(y_{ij}\mid b_i,\beta\sim N(\eta_i,\sigma^2)\), \(b_i\sim N(0,\tau^2)\), and \(\beta\sim N(0,s_\beta^2)\). For illustration, let the event hazard be constant over time for subject \(i\), conditional on their shared level: \(h_i(t)=\lambda\exp(\alpha\eta_i)\). Here \(\alpha\) is the association parameter, and \(\sigma,\tau,s_\beta,\lambda\) are treated as known positive constants. Take \(\alpha\sim N(0,s_\alpha^2)\). This deliberately simple specification makes the gradient visible without changing the central HMC construction.

The event contribution is \(\delta_i\log h_i-T_i h_i\): an observed event contributes the log hazard as well as the cumulative hazard term, while a censored observation contributes only the survival term. Up to constants not depending on \((\beta,\alpha,b_1,\ldots,b_n)\), the log posterior is

The normalizing constants omitted here do not depend on the HMC parameters, so they do not affect the gradient or acceptance ratio. Differentiating gives

These expressions show how the sources of information combine. For example, each \(b_i\) receives a longitudinal contribution from that subject's residuals, a shrinkage contribution from its normal prior, and an event contribution through the shared level \(\eta_i\). The fixed effect \(\beta\) appears in every subject's level, so its gradient collects contributions across subjects as well as its prior term.

Hamiltonian Dynamics and Leapfrog Updates

Let \(\mathbf q\) denote the vector of parameters being updated, such as \((\beta,\alpha,b_1,\ldots,b_n)\). Write \(\ell(\mathbf q)\) for its log posterior. HMC introduces an auxiliary momentum \(\mathbf p\sim N(\mathbf 0,M)\), where \(M\) is a positive-definite mass matrix chosen by the implementation. The Hamiltonian, a negative log density on the augmented space up to an additive constant, is

The corresponding dynamics are \(d\mathbf q/dt=M^{-1}\mathbf p\) and \(d\mathbf p/dt=\nabla\ell(\mathbf q)\). If these equations could be followed exactly, they would conserve \(H\). In practice, HMC uses a numerical integrator. With step size \(\epsilon\), one leapfrog step is

Repeating these steps for a chosen number of steps produces a proposed endpoint. The Metropolis acceptance probability is \(\min\{1,\exp[H(\mathbf q_0,\mathbf p_0)-H(\mathbf q_1,\mathbf p_1)]\}\). Implementations commonly reverse the endpoint momentum when forming the proposal; this makes the numerical trajectory reversible. Resampling momentum for each trajectory and applying the acceptance correction are essential parts of the update, not optional refinements.

Implementation Choices for Joint Models

For continuous fixed effects and latent random effects, one possible HMC state is the entire vector \((\beta,\alpha,\mathbf b)\), with the joint log-posterior gradient evaluated at every leapfrog position. The subject-level structure in the gradient can make evaluation efficient: longitudinal and event contributions can be accumulated subject by subject. The latent effects are part of the state, not integrated out in this update. Which parameters are updated together is an implementation decision; the gradient formula itself does not prescribe a particular block strategy.

Scaling matters because HMC advances all coordinates using a common step size. If one coordinate has much higher posterior curvature than another, a step size large enough to move effectively in the flatter direction may be too large for the sharply curved direction. The mass matrix changes the momentum-to-position scaling. In the convention above, for a one-dimensional Gaussian target with variance \(s^2\), choosing \(M=1/s^2\) makes the local oscillation frequency \(\sqrt{(1/s^2)/M}=1\). In practice, a diagonal or more general estimated mass matrix can help account for different scales and correlations. Software may describe its tuning quantity using an inverse-mass convention, so verify the convention rather than transferring matrix values by name alone.

If a parameter must be positive, such as a standard deviation or hazard rate, one option is to sample an unconstrained log parameter \(z\) and set \(\theta=\exp(z)\). The target density in \(z\) must include the change-of-variables Jacobian: \(\log|d\theta/dz|=z\). Omitting this term targets the wrong posterior. More generally, gradients must be taken with respect to the coordinates actually used by the sampler.

Diagnostics, Alternatives, and Limitations

A high rejection rate or frequent large changes in Hamiltonian energy can signal that the step size is too large, gradients are inaccurate, or the posterior has difficult curvature. Very small steps may avoid numerical failure but require many gradient evaluations to travel a useful distance. Check analytic gradients against finite-difference checks at representative parameter values, and inspect whether the log posterior and gradient remain finite over the trajectory. Extreme exponentials in an event likelihood, coding errors in the censoring contribution, and an omitted transformation Jacobian are particularly consequential failure modes.

HMC is most useful when the target is continuous and gradients are available and reliable. Random-walk Metropolis requires no gradients and can accommodate a broader range of proposal designs, but often explores strongly correlated continuous targets inefficiently. Gibbs updates can be effective when full conditional distributions are available in a tractable form; HMC is not a substitute for an exact conditional update that is already convenient. In a joint model, HMC can also be less effective when the posterior geometry is poorly scaled, contains sharp or narrow regions, or has nondifferentiable boundaries. Discrete model indicators require other update mechanisms.

A common mistake is to interpret a large HMC proposal as evidence that the chain is mixing well. A long numerical trajectory can instead accumulate large integration error and be rejected, while an accepted trajectory can still explore the posterior poorly if the parameterization or scaling is unfavorable. Acceptance alone does not establish convergence. Evaluate chain behavior with appropriate convergence and Monte Carlo diagnostics, and investigate divergences or unstable energy behavior rather than treating them as cosmetic warnings.