Multimodal Integrands and Failure of a Single-Mode Approximation

Why a Single Mode Can Be Misleading

Laplace approximation replaces an integrand near a mode by a Gaussian-shaped curve determined by its height and curvature. This can work well when nearly all the integral is concentrated in one neighborhood. But the random-effects integrand may instead have several separated peaks. A calculation centered on just one then describes only one part of the mass; it can give a poor approximation even if the mode and its local curvature have been computed accurately.

Multimodality is not inevitable. In particular, a useful sufficient condition for a unique mode is strict concavity of the log integrand. The warning is that this condition is model-dependent: checking that one optimizer converged, or that one Hessian has the expected curvature, does not establish that other modes do not exist.

Definition: For a random-effects vector \(b\), write the marginal likelihood contribution for one cluster as \(I=\int \exp\{h(b)\}\,db\), where \(h(b)\) is the cluster's conditional log likelihood plus the log random-effects density, with any terms constant in \(b\) retained or omitted consistently. A mode is a local maximum of \(h(b)\). The integrand is multimodal if \(h(b)\) has more than one local maximum.

When Uniqueness Is Guaranteed

Some common GLMM integrands have a shape that rules out separated modes. For example, for a binomial response with a canonical logit link and an affine predictor, the conditional log likelihood is concave in the random effects. For a Poisson response with a log link and an affine predictor, it is also concave. An affine transformation of the random effects preserves concavity. If the random-effects distribution is proper Gaussian, its log density contributes a negative-definite quadratic term. The sum is then strictly concave, so it has at most one mode.

This reasoning has a boundary: it applies when the conditional log likelihood is concave in the linear predictors. It does not apply to every response distribution and link. In particular, a noncanonical or otherwise nonconcave conditional likelihood may remain nonconcave after the Gaussian penalty is added. A positive-definite Gaussian precision matrix supplies negative curvature, but that curvature need not dominate the likelihood's positive curvature everywhere.

Worked Example: Checking Concavity for a Binomial Logit Cluster

Consider a hypothetical cluster with \(m=8\) binary trials, \(y=3\) successes, and linear predictor \(\eta=x^\mathsf{T}\beta+b\), where \(b\sim N(0,\tau^2)\) and \(\tau^2=4\). Ignoring terms independent of \(b\), the log integrand is

h(b)=y\eta-m\log(1+\exp(\eta))-\frac{b^2}{2\tau^2}.

Let \(p=\exp(\eta)/(1+\exp(\eta))\). Differentiating twice with respect to \(b\), noting \(d\eta/db=1\), gives

h''(b)=-mp(1-p)-\frac{1}{\tau^2}.

At \(b=0\), suppose \(x^\mathsf{T}\beta=0\), so \(p=0.5\). Then \(h''(0)=-8(0.5)(0.5)-1/4=-2.25\). More importantly, for every finite \(b\), \(p(1-p)>0\), so \(h''(b)<0\) everywhere. The log integrand is strictly concave and cannot have multiple separated modes. This check is a model-specific justification, not a general property of all GLMMs.

How Multiple Modes Can Arise

A nonconcave likelihood can create regions of positive curvature. If those regions are not outweighed by the Gaussian penalty, the log integrand can have a local minimum between two peaks. The following hypothetical Gaussian-response example uses a log link and a random slope. It demonstrates why the concavity guarantee above cannot be assumed for every GLMM.

Worked Example: Two Modes in a Random-Slope Integrand

Suppose two independent hypothetical Gaussian responses have variance \(\sigma^2=0.25\), observed values \(y_1=y_2=4\), and conditional means \(\mu_i=\exp(x_i b)\), with \(x_1=1\), \(x_2=-1\). Let the random slope be \(b\sim N(0,100)\). Apart from constants independent of \(b\), the log integrand is

h(b)=-2\left[(4-\exp(b))^2+(4-\exp(-b))^2\right]-\frac{b^2}{200}.

The expression is symmetric: \(h(b)=h(-b)\). At \(b=0\), both fitted means are 1, so

$$ h(0)=-2[(4-1)^2+(4-1)^2]=-36. $$

At \(b=\log(4)\), the fitted means are 4 and \(1/4\). Thus

h(\log 4)=-2\left[0^2+\left(4-\frac14\right)^2\right]-\frac{(\log 4)^2}{200} =-28.125-0.0096090603\ldots=-28.1346090603\ldots\approx-28.1346.

The value at \(\log(4)\) is higher than at zero. Also, \(h'(0)=0\). For a Gaussian observation with mean \(\exp(xb)\), the likelihood's second derivative is \(x^2\exp(xb)(y-2\exp(xb))/\sigma^2\). Therefore, at zero the two observations contribute \(2[1(4-2)/0.25]=16\), while the Gaussian density contributes \(-1/100=-0.01\). Hence \(h''(0)=15.99>0\): zero is a local minimum, not a mode. Since the log integrand tends to \(-\infty\) in both tails and is symmetric, it must have a maximum on each side of zero. These are distinct modes. The data are deliberately hypothetical; the point is the shape of the resulting integrand.

Why the Local Laplace Formula Can Miss Mass

If \(b_j\) is a mode and the negative Hessian \(H_j=-\nabla^2h(b_j)\) is positive definite, the local Laplace contribution in \(d\) dimensions is

I_j\approx \exp\{h(b_j)\}(2\pi)^{d/2}|H_j|^{-1/2}.

The height \(\exp\{h(b_j)\}\) alone does not determine how much integral lies near that mode. A lower, broader peak can contribute more than a higher, narrower one because its curvature determinant is smaller. When several well-separated modes contribute appreciable mass, the corresponding leading approximation is the sum of their local contributions, rather than the contribution from just one.

Justification: Choose disjoint neighborhoods around the modes and split the integral into the integrals over those neighborhoods plus the remaining region. Near each mode, a second-order Taylor expansion gives \(h(b)\approx h(b_j)-\tfrac12(b-b_j)^\mathsf{T}H_j(b-b_j)\). Integrating the exponential of this quadratic over the neighborhood, extended to all of \(\mathbb{R}^d\) when the omitted tails are negligible, gives \(\exp\{h(b_j)\}(2\pi)^{d/2}|H_j|^{-1/2}\). Adding the neighborhood integrals gives the sum. This approximation requires each important mode to be locally well represented by its quadratic expansion and the mass outside the neighborhoods to be negligible. It does not justify discarding an important neighborhood.

Worked Example: A Broader Peak Can Matter More

Suppose a hypothetical one-dimensional random-effects calculation reports two separated local maxima. The first has log height \(h_1=-3\) and curvature \(H_1=100\); the second has log height \(h_2=-3.2\) and curvature \(H_2=4\). These are hypothetical numerical mode summaries, not measurements from a study.

The local contributions are

C_1\approx e^{-3}\sqrt{\frac{2\pi}{100}} =0.04979(0.2507)=0.01248,

C_2\approx e^{-3.2}\sqrt{\frac{2\pi}{4}} =0.040762(1.253314)=0.05109.

Although the second peak is lower by \(0.2\) log units, its local contribution is about \(0.05109/0.01248=4.094\) times the first's because it is much broader. The two-mode approximation is \(C_1+C_2\approx0.06357\). Keeping only the first gives \(0.01248\), about \(80.4\%\) below this two-mode sum. This comparison illustrates the local approximations' relative contributions; it is not a guarantee that the sum equals the exact integral.

Warning Signs and Practical Checks

A mode finder is a local procedure: its convergence confirms a stationary point has been found, not that it is the only one. In a scalar random effect, a plot of \(h(b)\) over a sufficiently broad range can reveal separated peaks and a valley. In multiple dimensions, useful checks include repeated optimization from dispersed starting values, inspecting distinct converged solutions and their log-integrand values, and checking that each candidate's negative Hessian is positive definite. A stationary point with the wrong curvature may be a minimum or saddle point rather than a mode.

Numerical results should be treated cautiously if fitted modes change substantially with starting values, if a Hessian is nearly singular, or if the integrand remains appreciable far from the mode being used. A very low valley between peaks is a particularly clear warning that one local quadratic expansion cannot describe the full shape. Conversely, repeated convergence to the same mode is not proof of uniqueness: starting points may all lie in the same basin of attraction.

Worked Example: Interpreting a Multi-Start Diagnostic

A hypothetical optimizer is run from five starting values for a scalar random effect. Three runs converge to \(b=-1.3\) with \(h(b)=-28.1\), and two converge to \(b=1.3\) with \(h(b)=-28.1\). At both solutions, the computed second derivative is \(-15\).

Both points pass the local maximum check because \(h''(b)<0\) in one dimension. The equal log heights suggest symmetry, but the run counts do not estimate the modes' relative integral contributions: they depend on the chosen starting values and optimizer. A single-mode calculation at \(b=-1.3\) would ignore the distinct maximum at \(1.3\). The appropriate warning is therefore not “the optimizer failed”; rather, the integrand has at least two candidate modes, and a one-mode approximation requires additional evidence that the other mode contributes negligibly. Here the equal heights and curvatures give no such evidence.

Alternatives, Limits, and a Common Mistake

The main alternative to a single-mode calculation is to account for all important modes, either by summing local Laplace contributions or by using a numerical integration method that can resolve separated regions of mass. Adaptive quadrature is one such alternative; its effectiveness depends on whether the integration procedure explores all relevant regions rather than concentrating only near one mode. The introductory material on Laplace approximation explains the local quadratic construction; the key issue here is whether one such local construction covers the important mass.

Multi-mode checks also have limits. Searching from multiple starting values can miss a narrow mode, and a local quadratic approximation can be poor when a mode is skewed or its neighborhood is not close to Gaussian. In high dimensions, locating and integrating around every mode can be computationally difficult. A single-mode approximation may still be adequate if one mode's contribution dominates, but that conclusion should be supported by checking plausible alternatives, not inferred solely from local convergence.

Common mistake: “The optimizer converged and the Hessian at the solution is negative, so the log integrand has one mode.”
Correction: These facts establish only a local maximum. They do not rule out other modes. Compare solutions from dispersed starts, assess the shape over a suitable region, and consider whether other modes could carry appreciable integral mass.

Key takeaway: Laplace approximation is local. A reliable single-mode calculation requires both a well-behaved local peak and evidence that other peaks do not contribute materially to the random-effects integral.

Check Your Understanding

Answer each question, then compare with the brief solutions.

If a stationary point in one dimension has \(h''(b)>0\), can it be used as a local Laplace mode? Why or why not?
Using the two hypothetical contributions \(C_1=0.01248\) and \(C_2=0.05108\), which mode contributes more, and by approximately what factor?
Why does a proper Gaussian random-effects density not guarantee a unique mode for every GLMM?
For a symmetric integrand with two well-separated modes of equal curvature and height, what does a one-mode Laplace approximation omit?

Solutions. (1) No. Positive second derivative indicates a local minimum in one dimension; a Laplace mode requires negative second derivative. (2) The second mode contributes more, by \(0.05108/0.01248\approx4.093\). (3) The Gaussian term adds negative quadratic curvature, but a nonconcave conditional likelihood can have positive curvature that the Gaussian term does not outweigh everywhere. (4) It omits the other mode's local contribution; under equal heights and curvatures the two leading local contributions are equal, so using one captures only one of the two contributions in the summed approximation.