Laplace Approximation from a Quadratic Taylor Expansion

Why Approximate an Integral Locally?

In a GLMM, integrating over random effects can produce an integral that is difficult to evaluate directly. The Laplace approximation replaces the log integrand near its dominant mode with a quadratic function. Exponentiating that quadratic gives a Gaussian-shaped function whose integral is available analytically. The approximation therefore uses two local features: the height of the integrand at the mode and its curvature there.

This tutorial derives that calculation for a general integral over a vector \(b\) of dimension \(q\). Here \(b\) can represent a vector of random effects, but the derivation itself only requires a suitable integrand. The result is a local approximation, not an identity for an arbitrary integrand.

From a Taylor Expansion to a Gaussian Integral

For an interior mode at which \(h\) is differentiable, the gradient satisfies \(\nabla h(b^*)=0\). Let \(\delta=b-b^*\). The second-order Taylor expansion around the mode is

where \(r(\delta)\) is the Taylor remainder, which is small relative to \(\|\delta\|^2\) as \(\delta\) approaches zero when \(h\) is twice differentiable at the mode. The gradient term vanishes, and \(\nabla^2 h(b^*)=-H\). Neglecting the remainder gives the local approximation

If \(H\) is positive definite, this quadratic decreases in every direction away from the mode. Exponentiating the approximation gives a peak of height \(\exp\{h(b^*)\}\), multiplied by a Gaussian-shaped term. Integrating that term over all \(q\) dimensions yields the Laplace approximation.

For completeness, diagonalize the symmetric positive-definite matrix \(H\) as \(H=Q\Lambda Q^\mathsf{T}\), where \(Q\) is orthogonal and \(\Lambda\) has positive eigenvalues. Under the change of variables \(z=\Lambda^{1/2}Q^\mathsf{T}\delta\), the quadratic form becomes \(z^\mathsf{T}z\), and the volume element changes by the factor \((\det H)^{-1/2}\). The integral of \(\exp(-z^\mathsf{T}z/2)\) over \(\mathbb{R}^q\) is \((2\pi)^{q/2}\), since it factors into \(q\) one-dimensional standard Gaussian integrals. Thus the local approximation to \(I=\int_{\mathbb{R}^q}\exp\{h(b)\}\,db\) is

The determinant term is not an arbitrary correction: it measures the volume of the local Gaussian-shaped peak. A sharply curved peak has a larger determinant and a smaller integral for a fixed peak height. A flatter peak has a smaller determinant and a larger integral.

Example in a Binary Random-Intercept Setting

The same steps apply when the integration variable is a random effect. Consider hypothetical binary responses from four observations in one cluster, with two successes and two failures. Suppose the fixed intercept is zero and the cluster random intercept \(u\) has a \(N(0,1)\) density. For this example, define \(h(u)\) as the log of the product of the four conditional Bernoulli probabilities and the random-effect density. The conditional success probability is \(p(u)=e^u/(1+e^u)\), so

This integrand is positive after exponentiation and integrable: the normal density factor decays in the tails, while the conditional probabilities are at most one. To locate its mode, use \(\log p(u)=u-\log(1+e^u)\) and \(\log\{1-p(u)\}=-\log(1+e^u)\). Then

At \(u=0\), \(p(0)=1/2\), so \(h'(0)=2-4(1/2)-0=0\). Also \(h''(u)=-4p(u)\{1-p(u)\}-1<0\) for every finite \(u\). Thus \(h\) is strictly concave, and \(u^*=0\) is its unique mode. At the mode, \(H=-h''(0)=4(1/2)(1/2)+1=2\), and \(h(0)=-4\log 2-\tfrac12\log(2\pi)\approx-3.6915\). The one-dimensional Laplace calculation is

The numerical value is an approximation to the integral of the specified integrand over \(u\). The example checks the relevant conditions directly: the mode is interior and unique, the curvature is positive after changing the Hessian’s sign, and the integrand is integrable.

Conditions, Diagnostics, and Alternatives

The calculation is most defensible when the integral is finite, the relevant mode is interior, the log integrand is sufficiently smooth near it, and \(H\) is positive definite. The mode should also capture the region that contributes most of the integral. These are conditions for interpreting a local approximation sensibly; positive curvature at one point alone does not guarantee accuracy over the entire integration domain.

In computation, check the mode and curvature rather than treating the output as self-validating. A materially nonzero gradient suggests that the proposed point is not a well-located mode. A negative or near-zero eigenvalue of \(H\) indicates that the quadratic does not define a well-behaved local Gaussian peak; it may signal a saddle point, an inadequately located mode, or weak curvature. A very small determinant can make the approximation extremely sensitive to numerical error in the Hessian. In multiple dimensions, the determinant must be calculated from the full Hessian, including off-diagonal entries.

Direct numerical integration can be preferable when the integral is low-dimensional and can be evaluated accurately over the important region. Quadrature methods approximate an integral using weighted evaluations of the integrand rather than relying only on its local second-order shape. The Laplace approximation is attractive when a compact local calculation is useful, particularly in higher-dimensional settings, but neither method should be selected without considering the integrand’s shape and the required accuracy.

A single local quadratic may be unreliable if substantial mass lies far from the mode, the integrand is strongly asymmetric near the mode, or there are multiple important modes. A boundary mode also requires care, because the derivation above integrates a full Gaussian over \(\mathbb{R}^q\), whereas a constrained domain may contain only part of that local Gaussian. These cases call for further investigation or a method that represents more of the integrand than its value and curvature at one point.

Common Mistakes

• Using the Hessian with the wrong sign. At a twice-differentiable local maximum, \(\nabla^2h(b^*)\) is negative semidefinite; for the standard nondegenerate Laplace approximation, it is assumed negative definite, so \(H=-\nabla^2h(b^*)\) is positive definite.
• Forgetting the square root. The denominator is \(\sqrt{\det H}\), not \(\det H\).
• Using only the peak height. The approximation also depends on curvature. Two peaks with equal height but different widths generally have different approximate integrals.
• Assuming a valid local calculation proves global accuracy. The formula describes the integrand near a mode. It does not, by itself, rule out important contributions elsewhere.