Why HMC Warnings Matter in Hierarchical Models
Hamiltonian Monte Carlo (HMC) can explore a posterior efficiently when its numerical trajectories follow the posterior’s geometry. In a hierarchical GLMM, however, a group-level scale may be weakly informed by the data, or different parameter combinations may have very different curvature. The resulting geometry can be difficult for a numerical integrator to follow. A sampler may then report divergent transitions or poor energy exploration.
These warnings are not simply measures of whether the chain has run long enough. They can indicate that the sampler is failing to explore some parts of the posterior reliably. Conversely, a warning is not by itself proof that a particular estimate is biased or that the model is unusable. Its meaning depends on where it occurs, whether it persists after appropriate changes, and what other diagnostics show.
Hamiltonian Energy and Divergent Transitions
Consider a hierarchical GLMM with observations \(y_{ij}\) in group \(j\), fixed-effect parameters \(\beta\), and group-level effects \(b_j\). For example, its conditional mean can be related to predictors through \(g(\operatorname{E}[y_{ij}\mid b_j])=x_{ij}^{T}\beta+z_{ij}^{T}b_j\), with a specified sampling distribution and a distribution for the group effects. HMC targets the posterior distribution of the model parameters. In computation, constrained parameters are often transformed to unconstrained coordinates; the transformed log density must include the appropriate change-of-variables adjustment.
Write \(u\) for the unconstrained parameter vector and \(p\) for its auxiliary momentum. The negative log posterior, up to an additive constant, is the potential energy \(U(u)\). Given a positive-definite mass matrix \(M\), the kinetic energy and Hamiltonian are
Exact Hamiltonian dynamics conserve \(H\). HMC instead uses a numerical integrator, commonly the leapfrog integrator, to approximate a trajectory. Let \(H_{\text{start}}\) and \(H_{\text{end}}\) denote the Hamiltonian at the beginning and end of a proposed trajectory, and define \(\Delta H=H_{\text{end}}-H_{\text{start}}\). A large absolute energy error signals that the numerical trajectory did not adequately track the Hamiltonian dynamics. The Metropolis acceptance probability for a proposal with this energy change is
This expression follows from the target joint density of position and momentum, which is proportional to \(\exp[-H(u,p)]\). The ratio of endpoint to starting joint densities is therefore \(\exp[-(H_{\text{end}}-H_{\text{start}})]\). The Metropolis rule uses this ratio, capped at one, to correct for numerical integration error. A divergent transition is a software warning that the integrator encountered an excessive numerical error or otherwise could not safely continue the trajectory. The precise detection rule and threshold are implementation-specific.
In hierarchical models, a common source of difficult geometry is a group-level scale such as a random-effect standard deviation. When that scale is near zero, the effects it scales may be tightly constrained; at larger values they may vary more broadly. In a centered parameterization, these changing scales can create a narrow, curved region often described as a funnel. A trajectory that uses steps suitable for a broad region may be too coarse in a narrow one. This is a geometric explanation, not a claim that every divergence comes from a funnel.
Worked Example: Interpreting a Divergence Rate
Suppose a hypothetical binary random-intercept GLMM is sampled for four chains, with 500 post-warmup transitions per chain. The sampler flags 14 divergent transitions in total. Calculate the pooled divergence rate and interpret it.
There are \(4\times500=2{,}000\) post-warmup transitions. The flagged proportion is
Thus, 0.70% of the recorded transitions were flagged. That percentage is a description of this run, not a guarantee that the remaining 99.30% represent the posterior adequately. Suppose further inspection shows that 11 of the 14 flagged transitions occur when the random-intercept standard deviation is below 0.10 on the model’s linear-predictor scale. This concentration is consistent with a difficult near-zero scale region, but the plot alone does not prove the cause. The next checks should include the parameterization, the model’s information about that scale, and whether the warning persists under an appropriate computational change. A longer run by itself may produce more flagged transitions without resolving the geometry.
Energy-Based Diagnostics and E-BFMI
A chain can have few or no divergences and still move poorly through the posterior’s energy distribution. The energy Bayesian fraction of missing information (E-BFMI) assesses how effectively momentum refreshment and HMC transitions move between energy levels. For a sequence of sampled energy values \(E_1,\ldots,E_N\), a commonly used estimate is
where \(\widehat{\operatorname{Var}}(E)\) is the sample variance of the energy values in that chain. Software may use a slightly different finite-sample convention, so use its reported definition when comparing values. E-BFMI is normally considered chain by chain; a pooled value can conceal a poorly exploring chain. Values below about 0.3 are commonly treated as a warning threshold in Stan workflows, not as a universal mathematical boundary between adequate and inadequate exploration.
The numerator measures squared changes in energy between successive draws. If successive energies change very little relative to the overall spread of energy values, the chain may be moving slowly across the energy distribution. A low value can therefore signal poor exploration, even if no individual trajectory is flagged as divergent.
Worked Example: Calculating E-BFMI
A hypothetical chain produces six consecutive energy values, in arbitrary energy units: 10.0, 10.2, 10.4, 10.6, 10.8, and 11.0. Calculate the sample E-BFMI using the formula above.
The five successive differences are each 0.2, so the mean squared difference is
The energy mean is 10.5. The squared deviations from that mean sum to \(0.5^2+0.3^2+0.1^2+0.1^2+0.3^2+0.5^2=0.70\). The sample variance is therefore \(0.70/(6-1)=0.140\). Hence
This value is slightly below the commonly used 0.3 warning threshold. It suggests that this chain may be exploring energy levels inefficiently, but six values are far too few for a reliable practical assessment. In an actual analysis, inspect the sampler’s chain-specific diagnostic over the full sampling run and consider it alongside divergences, trace plots, and other chain diagnostics.
Reading the Diagnostics Together
Divergences and E-BFMI describe related but distinct aspects of HMC computation. A divergence points to a trajectory with problematic numerical integration at a particular location. E-BFMI summarizes movement across energy levels over a chain. Neither replaces checks of between-chain mixing or Monte Carlo precision. For example, \(\widehat{R}\), effective sample size, and Monte Carlo standard errors address different questions; a favorable value for one does not cancel a warning in another.
A useful investigation connects the warnings to model parameters and posterior geometry. Plot divergent draws against relevant parameters, especially group-level scales and correlations, and compare them with non-divergent draws. Look for concentration in a narrow region, near a boundary, or along a strongly curved relationship. Also inspect chains separately: a pooled summary can hide a chain with low E-BFMI or many divergences. If the same region repeatedly produces divergences, treat posterior summaries that depend on that region cautiously until the computational problem is understood.
Record the number of divergent transitions and check which parameters or parameter combinations are associated with them. Calculate rates only as descriptive summaries, not as proof of validity.
Review E-BFMI for each chain and identify whether a low value is isolated or repeated. Inspect energy behavior alongside trace plots and standard mixing diagnostics.
Ask whether a weakly informed scale, a centered random-effects representation, or a parameter boundary could produce sharp curvature. Use an alternative parameterization when it is scientifically and computationally appropriate.
After a targeted change, rerun the sampler and reassess all chains. A warning count that falls is encouraging only if mixing and energy exploration are also adequate and the model’s posterior target has not been unintentionally changed.
Worked Example: Comparing Centered and Noncentered Fits
A hypothetical study has many observations in some groups but very little information about the random-intercept variation in other groups. Two fits use the same likelihood, priors, data, and HMC settings, but different centered and noncentered representations of the random effects. Across four chains and 2,000 post-warmup transitions, the centered fit has 46 divergences and a minimum chain E-BFMI of 0.18. The noncentered fit has no divergences and a minimum chain E-BFMI of 0.72.
For the centered fit, the pooled divergence rate is \(46/2{,}000=0.023\), or 2.30%. For the noncentered fit it is \(0/2{,}000=0\%\). The minimum E-BFMI changes from 0.18, below the commonly used warning threshold, to 0.72. The pattern is consistent with the noncentered representation making this weakly informed hierarchical geometry easier for HMC to explore.
This comparison does not establish that noncentering is always preferable, nor does it alone prove that the second run is adequate. Centering can work better when group-level effects are strongly informed. Check the new run’s chain mixing, effective sample sizes, and Monte Carlo error, and confirm that both parameterizations represent the same intended posterior. The diagnostic improvement is evidence about computational behavior, not evidence that the scientific model is true.
Remedies, Alternatives, and Limitations
A smaller HMC step size can reduce integration error and may remove divergences. In adaptive samplers, increasing the target acceptance probability generally encourages a smaller step size. This is a reasonable diagnostic experiment, but it is not a universal cure: very small steps can make computation expensive, and persistent warnings may reflect deeper geometry. Increasing the maximum tree depth in a NUTS implementation addresses trajectories that grow too long; it does not directly resolve a trajectory that diverges because the step size cannot follow sharp curvature.
Reparameterization is often more directly relevant when the warning pattern suggests hierarchical funnel geometry. The earlier material on centered and noncentered parameterizations explains their forms and when each may suit the information in the data. Other checks include verifying that predictor scales and priors are intentional, and confirming that the model is sufficiently identified by the data. Changing a prior merely to suppress a warning is not appropriate unless the revised prior is substantively justified.
These diagnostics also have limitations. A chain with no divergences can still mix poorly, while occasional divergences do not quantify the exact size or direction of any resulting inferential error. E-BFMI is a warning summary, not a guarantee when it is above a threshold. Gibbs sampling or specialized samplers such as Pólya–Gamma methods may be useful alternatives for models with suitable conditional structure, but their availability does not make HMC diagnostics irrelevant when HMC is used.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- In a four-chain run with 750 post-warmup transitions per chain, 9 transitions diverge. What is the pooled divergence rate, and what can it establish?
- Why can a hierarchical GLMM have difficult HMC geometry when a random-effect scale is near zero?
- A chain’s E-BFMI is 0.22. What does this suggest, and what does it not prove?
- Why does increasing the maximum tree depth not directly fix a divergence caused by sharp local curvature?
- After changing parameterization, divergences disappear. Name two further checks needed before treating the run as satisfactory.
Solutions. (1) There are \(4\times750=3{,}000\) transitions, so \(9/3{,}000=0.003\), or 0.30%. This describes the run but does not establish that posterior exploration is adequate. (2) The scale changes the range of plausible group effects; in a centered representation, the posterior can narrow and curve as the scale approaches zero, challenging a fixed numerical step size. (3) It is below the common 0.3 warning threshold and suggests inefficient movement across energy levels. It does not prove bias, nonconvergence, or the cause of the problem. (4) Tree depth limits trajectory length; it does not make a numerical step fine enough to resolve sharp curvature. (5) Check chain mixing and effective sample sizes, and assess Monte Carlo error; also verify the intended posterior is unchanged and inspect chain-specific energy diagnostics.