Divergences and Funnel Geometry

Why Divergences Matter

Hamiltonian Monte Carlo (HMC) uses gradients to propose movement through a posterior distribution. When the posterior’s geometry changes sharply across regions, the numerical approximation used to follow that movement can become unreliable. A divergent transition is a warning that this approximation encountered a region of problematic curvature and that the proposed trajectory could not be trusted. It is not simply another name for high autocorrelation, nor is it proof that every posterior summary is wrong. But divergences can indicate that parts of the posterior are being explored inadequately, so they warrant investigation before interpreting estimates.

Hierarchical joint models are especially susceptible because subject-specific random effects are estimated together with population-level scales and with information from both longitudinal and event-time data. A scale close to zero can constrain many random effects to a narrow region, while a larger scale permits a much wider region. This changing width creates the geometry commonly called a funnel.

Hierarchical Posterior Geometry

Consider a simplified joint model with \(n\) subjects. Let \(b_i\) be a subject-specific random effect and \(\tau\) its population standard deviation. The longitudinal measurements and event outcome for subject \(i\) contribute a likelihood \(L_i(b_i,\theta)\), where \(\theta\) collects other model parameters. A generic posterior kernel is

p(\theta,\tau,b_1,\ldots,b_n\mid \text{data}) \propto p(\theta)p(\tau)\prod_{i=1}^{n} L_i(b_i,\theta)\,p(b_i\mid\tau).

For a normal random-effects distribution, \(b_i\mid\tau\sim N(0,\tau^2)\). The likelihood \(L_i\) may combine a longitudinal contribution and a survival contribution, including their dependence through shared random effects. The funnel phenomenon is already visible in the hierarchical prior, though the data and other model components can change its shape.

Definition: In a centered parameterization, the random effects \(b_i\) are sampled directly, with their distribution written conditionally on the scale, such as \(b_i\mid\tau\sim N(0,\tau^2)\). In a noncentered parameterization, one instead writes \(b_i=\tau z_i\), where \(z_i\sim N(0,1)\), and samples \(z_i\) and \(\tau\).

In the centered form, the conditional density contains a factor \(\tau^{-n}\). When \(\tau\) is small, the random effects must all lie close to zero; when \(\tau\) is larger, they can occupy a much wider range. Thus the narrow end of the funnel corresponds to small \(\tau\), and its broad end to larger \(\tau\). In a joint model, the longitudinal and survival likelihoods can further bend or narrow this geometry. The precise location of difficult regions is model- and data-dependent.

Worked Example: How the Random-Effects Width Changes

Suppose a hypothetical joint model has a random intercept \(b_i\mid\tau\sim N(0,\tau^2)\). Compare the conditional spread of \(b_i\) at two plausible values of \(\tau\), \(0.05\) and \(0.80\).

At \(\tau=0.05\), the conditional standard deviation of \(b_i\) is \(0.05\), and its variance is \(0.05^2=0.0025\). At \(\tau=0.80\), the standard deviation is \(0.80\), and the variance is \(0.80^2=0.64\). The ratio of the standard deviations is \(0.80/0.05=16\); the ratio of the variances is \(0.64/0.0025=256\).

These values do not establish that either scale is likely under the posterior. They show why a single numerical step size may have difficulty traversing the full range: the random-effect distribution can be far narrower at one scale than another. With many subjects, all their random effects are subject to this scale-dependent change.

Why a Funnel Can Cause Divergences

HMC introduces momentum and follows a Hamiltonian trajectory whose potential energy is the negative log posterior. In practice, the trajectory is approximated using discrete numerical steps, commonly leapfrog steps. A divergent transition occurs when the integrator’s energy error becomes too large under the sampler’s divergence criterion. The precise numerical threshold is implementation-dependent; the diagnostic is a warning about numerical integration in a challenging region, not a direct test of a statistical hypothesis.

The narrow end of a funnel requires fine resolution in the random-effect directions. A step size that works in a broad region may be too large there. If the numerical trajectory makes an inaccurate move, the resulting energy error can trigger a divergence. The problem is not that the posterior density is mathematically invalid: it is that the sampler’s numerical path may fail to represent that density accurately in some regions.

A short calculation shows how the funnel arises in the centered prior. Ignoring the likelihood for this geometric argument, the joint density of the random effects and scale is proportional to

p(\tau)\prod_{i=1}^{n}p(b_i\mid\tau) = p(\tau)(2\pi)^{-n/2}\tau^{-n} \exp\!\left(-\frac{\sum_{i=1}^{n}b_i^2}{2\tau^2}\right), \qquad \tau>0.

For small \(\tau\), appreciable density requires \(\sum_i b_i^2\) to be small as well. For larger \(\tau\), the exponent permits a wider range of \(b_i\). This coupled change in width is the funnel. The likelihood multiplies this expression, so conclusions about posterior geometry must be checked in the fitted model rather than inferred from the prior alone.

Worked Example: Interpreting an Energy-Error Warning

Suppose a hypothetical HMC implementation flags a transition when the absolute energy error exceeds 10. For one proposed trajectory, let the change in Hamiltonian energy be \(\Delta H=12\). The absolute error is \(|\Delta H|=12\), which exceeds the stated threshold, so this transition is flagged as divergent.

For an ideal Metropolis correction, the acceptance probability for a proposal with positive energy error 12 would be \(\min(1,\exp(-12))\), or approximately \(6.144\times10^{-6}\). That calculation illustrates how severe a large positive energy error can be; it does not define the divergence rule for every software package. A software-specific diagnostic can use its own threshold and implementation details.

The appropriate interpretation is that the numerical trajectory encountered a region where its approximation was judged unreliable. It is not correct to interpret the flag as a posterior probability, or to conclude from one flag alone that a particular coefficient estimate is biased. Its location and pattern, together with other chain diagnostics, determine how concerning it is.

Centered and Noncentered Parameterizations

The noncentered representation writes \(b_i=\tau z_i\), with \(z_i\sim N(0,1)\). In the prior, this removes the scale-dependent width from the distribution of the standardized effects. To verify the transformation, the centered conditional density contributes a factor \(\tau^{-n}\); changing variables from \(b_1,\ldots,b_n\) to \(z_1,\ldots,z_n\) has Jacobian determinant \(\tau^n\). The factors cancel:

\left\{\prod_{i=1}^{n}p(b_i\mid\tau)\right\} \left|\frac{\partial(b_1,\ldots,b_n)} {\partial(z_1,\ldots,z_n)}\right| = \left\{\prod_{i=1}^{n} \frac{1}{\sqrt{2\pi}\tau} \exp\!\left(-\frac{(\tau z_i)^2}{2\tau^2}\right)\right\}\tau^n = \prod_{i=1}^{n}\frac{1}{\sqrt{2\pi}}\exp(-z_i^2/2).

The transformed likelihood still depends on \(\tau z_i\), so this change does not guarantee a simple or independent posterior. It changes how the hierarchy is represented to the sampler. When individual random effects are weakly informed by the data, noncentering often reduces the tight dependence between the effects and their shared scale. When the data strongly inform each random effect, the centered representation can be more efficient because the likelihood directly constrains \(b_i\). Joint longitudinal and survival information both contribute to how strongly an effect is informed.

Worked Example: Comparing Conditional Scales Under Two Parameterizations

Consider one hypothetical subject-specific effect with the simplified observation model \(y\mid b\sim N(b,0.2^2)\), prior \(b\mid\tau\sim N(0,\tau^2)\), observed value \(y=1\), and fixed \(\tau=0.1\). This is an illustrative calculation of conditional information, not a full joint-model fit.

In the centered parameterization, normal-normal updating gives conditional posterior variance \(V_b=(1/\tau^2+1/0.2^2)^{-1}\). Substituting the values, \(V_b=(1/0.01+1/0.04)^{-1}=(100+25)^{-1}=0.008\), so the conditional posterior standard deviation is \(\sqrt{0.008}\approx0.08944\).

Under the noncentered parameterization \(b=0.1z\), the observation model is \(y\mid z\sim N(0.1z,0.04)\), with prior \(z\sim N(0,1)\). Its posterior precision is \(1+(0.1)^2/0.04=1.25\), so its posterior variance is \(1/1.25=0.8\), and its standard deviation is \(\sqrt{0.8}\approx0.8944\). Mapping back to \(b=0.1z\) gives standard deviation \(0.1(0.8944)=0.08944\), matching the centered calculation.

The two parameterizations describe the same distribution for \(b\). Their computational geometries differ: here, the observation is relatively informative about \(b\) compared with its prior scale, so the centered effect is directly constrained by the likelihood. This example does not prove that centering will be faster in a full model; it illustrates why the best representation depends on the information in the data.

Diagnosing and Responding to Divergences

A useful diagnosis asks where divergences occur and whether the pattern is consistent with a geometry problem. Plotting divergent draws against \(\tau\) and selected random effects can reveal concentration near a small-scale region or along a curved relationship. In a joint model, inspect relevant shared effects and parameters that enter both the longitudinal and event-time components. A cluster of divergences in a narrow region is more informative than the overall count alone.

Also compare chains and inspect trace plots, \(\hat{R}\), and effective sample sizes. These assess aspects of convergence and Monte Carlo quality, but they do not make divergences harmless. A chain may have apparently reassuring summaries while still failing to explore a region adequately. Conversely, an isolated divergence does not, by itself, establish that a reported estimate is materially wrong. The evidence must be considered together.

Locate the warnings.
Identify which parameters and posterior regions are associated with divergent draws; check whether they concentrate near small random-effects scales or another constrained region.

Check model and computation.
Confirm that the likelihood, constraints, and parameter scales are coded as intended. Poor scaling or a misspecified model can compound difficult geometry.

Try a justified reparameterization.
Fit a noncentered version when random effects are weakly informed, or compare with centering when the data strongly inform them. Preserve the same intended model and priors under the transformation.

Reassess the fit.
Check whether divergences have disappeared or moved, and review chain agreement, effective sample sizes, and the location of any remaining warnings before relying on summaries.

Increasing the target acceptance rate can lead an adaptive HMC sampler to use a smaller step size, which may reduce numerical integration error. This is a computational adjustment, not a guarantee that the geometry has been resolved. It can also increase computation time. Raising the maximum tree depth addresses trajectories that require more integration steps; it is not a general cure for divergences. Neither adjustment substitutes for checking the parameterization and the posterior region where warnings occur.

Common Pitfalls and Limitations

Calling divergences “just low effective sample size.” Slow mixing and divergences are distinct diagnostics. Low effective sample size concerns the amount of independent information in the chain; a divergence signals a potentially unreliable numerical trajectory.
Assuming noncentering always fixes funnels. Noncentering removes scale dependence from the random-effects prior, but the likelihood still couples \(z_i\) and \(\tau\). It may help weakly informed effects and be less efficient for strongly informed ones.
Increasing a tuning setting and stopping there. Fewer warnings after a smaller step size are encouraging, but the resulting fit still needs diagnostic review. A remaining warning should be located and understood, not hidden by a single aggregate count.
Assuming every warning proves biased estimates. A divergence is a warning, not a direct measurement of bias. Its location and the stability of relevant posterior summaries matter; absence of divergences is also not proof of correctness.

Reparameterization does not change the target posterior when the transformation and its probability-density implications are handled correctly. It can change the ease with which HMC explores that posterior. In a joint model, the shared random effects connect multiple likelihood contributions, so the most effective parameterization can depend on the amount and balance of longitudinal and survival information. There is no universal rule that one representation will be best for every effect or every dataset.

Key takeaway: A divergence indicates that HMC encountered a region where its numerical trajectory was unreliable. Funnel-shaped random-effects geometry is a common cause in hierarchical joint models. Locate the warnings, understand the scale-dependent geometry, and assess centered and noncentered representations as alternatives rather than assuming that one is universally superior.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

In a centered model with \(b_i\mid\tau\sim N(0,\tau^2)\), what happens to the conditional spread of the random effects as \(\tau\) approaches zero, and why can this produce a funnel?
A software diagnostic reports a divergent transition. Does this establish that the posterior mean of a parameter is biased? Explain what should be checked next.
For \(b_i=\tau z_i\) with \(z_i\sim N(0,1)\), explain why the centered prior’s scale factors cancel after transformation.
When might a centered parameterization be preferable to a noncentered one?
Why is increasing maximum tree depth not, by itself, a general remedy for divergences?

Solutions. (1) The conditional standard deviation is \(\tau\), so it approaches zero and confines the effects to a narrow region. At larger \(\tau\), they can spread out, producing a scale-dependent width. (2) No. A divergence is a numerical warning, not a direct bias estimate. Locate the divergent draws in parameter space and assess chain agreement, effective sample sizes, and relevant summaries. (3) The centered density contributes \(\tau^{-n}\); the change of variables has Jacobian \(\tau^n\), so their product is one. (4) Centering may be more efficient when data strongly inform individual random effects, because the likelihood directly constrains \(b_i\). (5) Maximum tree depth limits trajectory length; divergences arise from unreliable numerical integration in difficult geometry, which a longer trajectory does not necessarily repair.