Multiple Conditional Modes and Laplace Error

Why Multiple Conditional Modes Matter

A Laplace approximation replaces an integral by a Gaussian approximation near a mode. This is effective when one mode accounts for nearly all of the integral and the log integrand is well approximated by a quadratic in its influential region. But an NLME conditional integrand can have two or more well-separated modes. If a calculation finds only one, it may miss an appreciable part of the marginal likelihood, even when the mode it finds and its local curvature are both computed accurately.

For subject \(i\), write the conditional log integrand, at fixed population parameters \(\theta\), as \(h_i(b_i;\theta)=\log p(y_i\mid b_i,\theta)+\log p(b_i\mid\theta)\). Here \(y_i\) is the subject’s data and \(b_i\) is the subject-specific random-effects vector. The subject’s marginal likelihood contribution is the integral of \(\exp\{h_i(b_i;\theta)\}\) over \(b_i\). A mode is a local maximum of this conditional log integrand. The central issue here is that there may be several such maxima, not just one.

This contribution depends on both the height of the log integrand and its curvature. A high, narrow mode can contribute less area than a lower, broader one. Thus, comparing mode heights alone does not establish which mode matters most for the integral.

From Local Curvature to Several Contributions

Near a nondegenerate mode \(\hat b_{ij}\), a second-order Taylor expansion gives \(h_i(b_i;\theta)\approx h_i(\hat b_{ij};\theta)-\tfrac12(b_i-\hat b_{ij})^\mathsf{T}K_{ij}(b_i-\hat b_{ij})\). The linear term vanishes because the gradient is zero at the mode. Exponentiating this quadratic approximation gives a Gaussian-shaped integrand centered at \(\hat b_{ij}\).

Extending the local Gaussian approximation over all of \(\mathbb R^q\) gives \(C_{ij}\). This approximates the neighborhood’s contribution when the Gaussian mass outside that neighborhood is negligible. If there are separated modes and the regions around them account for nearly all the integral, the total integral is approximated by adding their contributions. In contrast, a single-mode calculation uses just one \(C_{ij}\). The multi-mode sum is an approximation too: it can be poor if the local quadratic shapes are inaccurate, if substantial mass lies outside the selected neighborhoods, or if the modes are not sufficiently separated for local contributions to be treated as distinct.

For independent subjects conditional on the population parameters, the marginal likelihood is the product of the subject-specific integrals. Therefore the contribution from each subject’s modes must first be summed on the likelihood scale. Taking logarithms then gives a sum of log subject-specific integrals. Keeping only one mode for a subject changes that subject’s integral before it enters the overall log likelihood.

Height and Width Can Tell Different Stories

A mode’s contribution is proportional to its height on the integrand scale, \(\exp\{h(\hat b)\}\), multiplied by a width factor determined by curvature. In one dimension, the width factor is \(\sqrt{2\pi/K}\): smaller positive curvature means a broader peak and a larger factor. In several dimensions, the determinant summarizes the local volume: a smaller \(|K|\) gives a larger local volume factor, all else equal.

Justification and Conditions for the Sum

The multi-mode calculation follows by partitioning the integration region into neighborhoods around the modes and the remaining region. Within each neighborhood, use the quadratic expansion and integrate it; if the remainder of the region has negligible integral, the whole integral is approximately the sum of those local Gaussian integrals. This reasoning also identifies the assumptions: the modes should be interior, nondegenerate, and sufficiently separated; the Hessian at each should yield a positive-definite \(K_{ij}\); and each local quadratic approximation should represent its neighborhood adequately.

The approximation is most persuasive when each mode is locally close to quadratic and the integrand is concentrated near the modes. With weakly identified random effects, flat regions, boundary modes, or substantial mass between modes, a local Gaussian calculation can be unreliable even if all prominent maxima have been found. Multiple modes and asymmetric local curvature can also occur together; discovering several modes does not remove local approximation error.

Diagnostics, Alternatives, and Common Pitfalls

A practical diagnostic is to repeat conditional optimization from several dispersed starting values and compare the resulting mode locations, log-integrand values, and curvature matrices. Distinct converged locations with near-zero gradients and positive-definite \(K\) indicate candidate modes. Evaluating the log integrand along directions joining the modes can help reveal whether the modes are genuinely separated and whether a saddle or broad connecting region may matter. These checks are evidence, not a guarantee: an optimizer can repeatedly find the same basin while missing another.

The main alternatives answer slightly different computational needs. A single-mode Laplace calculation is economical and can be adequate when one mode dominates. A multi-mode Laplace calculation sums local contributions and is preferable when several separated modes have appreciable mass and local Gaussian shapes are reasonable. Direct numerical integration can serve as a check in low-dimensional examples, but becomes harder as the dimension increases. The key limitation of any mode-based calculation is that it approximates local neighborhoods; it does not automatically account for unlocated modes or non-negligible mass away from them.

• Pitfall: keeping the mode with the highest peak. The incorrect reasoning is that the largest value of \(h\) must supply almost all the integral. The correction is to compare \(C_j\), which includes both height and curvature, and then assess the total of the relevant contributions.
• Pitfall: treating optimizer convergence as proof of uniqueness. Convergence shows that an algorithm found a stationary solution from its starting point; it does not rule out other basins. Use dispersed starts and inspect distinct solutions.
• Pitfall: adding log contributions directly. If local contributions are \(C_1\) and \(C_2\), the integral approximation is \(C_1+C_2\), not \(\log C_1+\log C_2\). The log of the approximate integral is \(\log(C_1+C_2)\).