Laplace Approximation for Nonlinear Mixed-Effects Models

Why the Random-Effects Integral Needs Approximation

In a nonlinear mixed-effects model, an individual’s observations depend on population parameters and that individual’s random effects. Conditional on the random effects, the likelihood contribution may be straightforward to evaluate. The population likelihood, however, requires averaging that conditional contribution over the random-effects distribution. For nonlinear models, this averaging generally produces an integral without a convenient closed form.

The Laplace approximation replaces that integral by a calculation based on the log-integrand near its maximum. It uses both the height of the integrand at the maximum and its local curvature. The curvature matters because two integrands can reach the same maximum value yet place different total probability mass around that maximum.

Definition: For individual $i$, let $y_i$ denote the observed data, $\eta_i$ a $q$-dimensional random-effects vector, and $\theta$ the population parameters. Write the conditional likelihood as $p(y_i\mid\eta_i;\theta)$ and the random-effects density as $p(\eta_i;\theta)$. The individual marginal likelihood contribution is $$ L_i(\theta)=\int_{\mathbb{R}^q}p(y_i\mid\eta_i;\theta)p(\eta_i;\theta)\,d\eta_i. $$ Define the log-integrand $h_i(\eta_i;\theta)=\log p(y_i\mid\eta_i;\theta)+\log p(\eta_i;\theta)$. The conditional mode $\hat{\eta}_i$ maximizes $h_i$ over $\eta_i$, for fixed $y_i$ and $\theta$.

In the usual formulation, random effects have a multivariate normal distribution with mean zero and covariance matrix $\Omega$, and observations are conditionally independent given the individual’s random effects and the model parameters. The density of this random-effects distribution is part of the integrand: omitting it would change both the mode and the curvature. The likelihood also includes any normalization constants in the conditional data density. These constants may not affect the mode, but they do affect the likelihood value.

The approximation requires a well-defined, interior conditional mode and a log-integrand that is sufficiently smooth near that mode. Define the negative Hessian at the mode by $H_i=-\nabla_{\eta_i}^2 h_i(\hat{\eta}_i;\theta)$. For the standard local Gaussian approximation, $H_i$ must be positive definite. Thus the log-integrand has negative curvature in every random-effects direction at the mode.

Deriving the Approximation

The derivation is a local approximation, not a claim that the entire log-integrand is quadratic. Expand $h_i$ in a second-order Taylor series around its conditional mode. Its first derivative is zero at an interior maximum, so the linear term vanishes:

h_i(\eta_i;\theta)\approx h_i(\hat{\eta}_i;\theta) -\frac{1}{2}(\eta_i-\hat{\eta}_i)^\mathsf{T}H_i(\eta_i-\hat{\eta}_i).

Exponentiating this expression gives a constant peak height multiplied by a Gaussian-shaped function. Integrating that local Gaussian over $q$ dimensions yields the determinant term below.

L_i(\theta)\approx \exp\!\left\{h_i(\hat{\eta}_i;\theta)\right\} (2\pi)^{q/2}|H_i|^{-1/2}.

Here $|H_i|$ is the determinant of the negative Hessian. Taking logarithms gives the corresponding approximate individual log-likelihood contribution:

\log L_i(\theta)\approx h_i(\hat{\eta}_i;\theta)+\frac{q}{2}\log(2\pi)-\frac{1}{2}\log|H_i|.

Justification: The Taylor expansion makes the local integrand approximately $\exp\{h_i(\hat{\eta}_i;\theta)\}\exp\{-\tfrac12 u^\mathsf{T}H_i u\}$, where $u=\eta_i-\hat{\eta}_i$. For positive-definite $H_i$, the Gaussian integral is $$ \int_{\mathbb{R}^q}\exp\!\left(-\frac{1}{2}u^\mathsf{T}H_i u\right)\,du =(2\pi)^{q/2}|H_i|^{-1/2}. $$ Multiplying by the peak height proves the displayed Laplace formula under the local quadratic approximation.

The term $|H_i|^{-1/2}$ measures the volume of the local peak. Greater curvature means a narrower peak and, all else equal, a smaller integral. In one dimension, the adjustment is proportional to $1/\sqrt{H_i}$. In several dimensions, the determinant combines curvature in all directions, including any dependence between those directions.

If individuals contribute independently conditional on the population parameters, the marginal likelihood is the product of their contributions. Applying the approximation to each individual gives an approximate total log likelihood by summing the individual expressions. Each individual has a separate conditional mode and Hessian at the current population-parameter values.

Worked Example: A One-Dimensional Integral and the Curvature Factor

Consider the illustrative integral $L=\int_{-\infty}^{\infty}\exp\{-2-(\eta-1)^2\}\,d\eta$. This is a simple log-integrand chosen to make the mode and curvature transparent; it is not a real study. Find the Laplace approximation and check the arithmetic.

The log-integrand is $h(\eta)=-2-(\eta-1)^2$. Its derivative is $h'(\eta)=-2(\eta-1)$, which is zero at $\hat{\eta}=1$. Its second derivative is $h''(\eta)=-2$, so the negative Hessian is $H=2$, which is positive.

For $q=1$, the formula gives $L\approx\exp\{h(\hat{\eta})\}\sqrt{2\pi}/\sqrt{H}$. Substituting $h(1)=-2$ and $H=2$, $L\approx e^{-2}\sqrt{2\pi}/\sqrt{2}=e^{-2}\sqrt{\pi}$. Using $e^{-2}\approx0.1353$ and $\sqrt{\pi}\approx1.7725$, the result is $L\approx0.2399$.

Because this log-integrand is exactly quadratic, the Gaussian integral calculation is exact here. The result also shows why the peak height alone is not enough: the curvature $H=2$ determines the width multiplier.

Worked Example: A Gaussian Random-Intercept Model

Suppose a hypothetical single observation has the model $Y\mid\eta\sim N(\eta,1)$, with $\eta\sim N(0,1)$, and the observed value is $y=2$. Calculate the Laplace approximation to the marginal density of this observation. The conditional normal model and normal random-effects distribution satisfy the smoothness and interior-mode conditions, and the resulting negative Hessian will be positive.

Including both normalizing constants, the log-integrand is $h(\eta)=-\tfrac12(2-\eta)^2-\tfrac12\eta^2-\log(2\pi)$. Its derivative is $h'(\eta)=2-2\eta$, so the conditional mode is $\hat{\eta}=1$. Its second derivative is $-2$, hence $H=2>0$.

At the mode, the two squared residual terms are each 1, giving $h(1)=-\tfrac12-\tfrac12-\log(2\pi)=-1-\log(2\pi)\approx-2.8379$. Therefore, $L\approx\exp(-1-\log(2\pi))\sqrt{2\pi}/\sqrt{2} =e^{-1}/\sqrt{4\pi}\approx0.1038$.

The log-integrand is quadratic, so the Laplace result is exact. As a check, the marginal distribution is $Y\sim N(0,2)$; its density at $y=2$ is $(1/\sqrt{4\pi})\exp\{-2^2/(2\cdot2)\}=e^{-1}/\sqrt{4\pi}\approx0.1038$. Agreement verifies both the normalizing constants and the curvature factor.

Worked Example: Same Peak, Different Curvature in Two Dimensions

Compare two hypothetical log-integrands with the same peak height and mode: $h_A(\eta)=-3-\tfrac12(\eta_1^2+4\eta_2^2)$ and $h_B(\eta)=-3-\tfrac12(\eta_1^2+\eta_2^2)$. Both have mode $(0,0)$, peak height $\exp(-3)$, and smooth, strictly concave log-integrands with positive-definite negative Hessians.

For A, $H_A=\begin{pmatrix}1&0\\0&4\end{pmatrix}$, so $|H_A|=4$. With $q=2$, $L_A\approx e^{-3}(2\pi)/\sqrt{4}=e^{-3}\pi\approx0.1564$. For B, $H_B=\begin{pmatrix}1&0\\0&1\end{pmatrix}$, so $|H_B|=1$, and $L_B\approx e^{-3}(2\pi)/\sqrt{1}=2\pi e^{-3}\approx0.3128$.

The peak heights are identical, but A is more sharply curved along the second coordinate and has half the approximated integral. For these quadratic examples, both calculations are exact. In a nonlinear mixed-effects model, the same determinant adjustment reflects local curvature, while the accuracy of the approximation also depends on how well the quadratic expansion represents the log-integrand around the mode.

Assumptions, Diagnostics, and Failure Modes

The formula is most reliable when most of the integral’s mass lies near one well-behaved interior mode and the log-integrand is close to quadratic over that region. A valid calculation therefore involves more than finding a numerical maximum. The mode must be located accurately, and the local curvature must be estimated reliably.

Check the mode. The gradient with respect to the random effects should be close to zero at the reported mode. A large gradient suggests incomplete optimization or an incorrectly evaluated derivative.
Check curvature. The negative Hessian should be symmetric and positive definite at a strict local maximum. A zero or negative eigenvalue indicates a flat or upward-curving direction; the usual Gaussian integral factor is then not justified.
Check numerical scale. A very small eigenvalue makes $|H_i|^{-1/2}$ large and can signal a weakly identified direction or unstable determinant calculation. A very large eigenvalue indicates a narrow direction that may require careful numerical resolution.
Check local shape. Strong skewness, heavy tails, or multiple separated modes can make a single local quadratic approximation unrepresentative. A positive-definite Hessian at one mode does not establish that no other important mode exists.
Check likelihood stability. Where feasible, compare with a more direct integration calculation for selected individuals or parameter settings. Material disagreement can reveal that local curvature is not adequately capturing the integral.

The inverse Hessian $H_i^{-1}$ describes the covariance matrix of the local Gaussian approximation around the conditional mode. It is not automatically the exact conditional covariance of the random effects. That distinction matters especially when the log-integrand is notably nonquadratic.

How Laplace Differs from Linearization

First-order and FOCE methods approximate the likelihood using linearization of the model prediction with respect to random effects, with FOCE centering the expansion at a conditional estimate. The Laplace construction instead approximates the integral by expanding the log-integrand to second order around its conditional mode and integrating the resulting Gaussian form. It therefore includes a determinant adjustment based on the local curvature of the full log-integrand.

These are different approximations to the marginal likelihood, not interchangeable names for the same calculation. Laplace can account for local curvature in the random-effects integral, but it does not guarantee accuracy when the integral is dominated by broad, skewed, or multiple regions. Nor does a positive Hessian alone guarantee a good approximation: it verifies local curvature at the mode, not global quadratic shape. The choice among estimation methods depends on model and design features and is treated separately in the subsequent comparison of FO, FOCE, and Laplace.

Common Mistakes

Using the conditional likelihood without the random-effects density. The integral is over the product $p(y_i\mid\eta_i;\theta)p(\eta_i;\theta)$. Leaving out the density changes the quantity being integrated.
Using the Hessian with the wrong sign. At a maximum of $h_i$, its Hessian is negative definite. The positive-definite matrix in the formula is $H_i=-\nabla^2 h_i$, not $\nabla^2 h_i$.
Ignoring the determinant because the mode is correct. The mode gives peak height, while the determinant adjusts for the local volume around that peak. A likelihood contribution cannot generally be recovered from the peak alone.
Treating the approximation as exact for every model. It is exact for a quadratic log-integrand over the full integration domain. For a nonquadratic log-integrand, it is a local approximation whose quality depends on the shape and concentration of the integral.

Key takeaway: For each individual, Laplace approximation combines the log-integrand at its conditional mode with a curvature correction, $-\tfrac12\log|H_i|$. Correct modes and positive-definite Hessians are essential, but the approximation is reliable only when the local quadratic shape adequately represents the important mass of the integral.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

In one dimension, a log-integrand has mode height $h(\hat{\eta})=-1$ and negative Hessian $H=8$. Write the Laplace approximation to the integral in terms of these values.
Why must the random-effects density be included in the log-integrand, even if the conditional data likelihood is the primary object of interest?
A two-dimensional negative Hessian has eigenvalues 2 and 5. Is it positive definite, and what is its determinant?
What does it mean if a numerical conditional mode has a negative eigenvalue in its negative Hessian?

Solutions. (1) $L\approx e^{-1}\sqrt{2\pi}/\sqrt{8}=e^{-1}\sqrt{\pi/4}\approx0.3260$. (2) The marginal likelihood integrates the product of the conditional data density and the random-effects density; omitting the latter changes the integrand, its mode, and its curvature. (3) Both eigenvalues are positive, so the Hessian is positive definite; its determinant is their product, $2\cdot5=10$. (4) The log-integrand has an upward-curving direction at that point, so it is not a strict local maximum in all directions and the standard Laplace Gaussian integral is not justified there.