Asymptotic Error Order of the Laplace Approximation

Why the Error Has an Asymptotic Order

The Laplace approximation replaces an integrand near its dominant mode by a Gaussian kernel. Its accuracy depends not just on whether that local quadratic approximation is plausible, but on how concentrated the integral becomes. In a GLMM, that concentration is often associated with increasing information from more observations within a cluster. When information grows, the region that contributes most to the integral shrinks, and departures from a quadratic shape matter progressively less.

The rate is not a guarantee for every GLMM integral. It is an asymptotic statement under conditions such as a unique interior mode, sufficiently smooth log integrand, and concentration around that mode. Under these conditions, the leading relative error is typically of order \(N^{-1}\), where \(N\) indexes the increasing information. The coefficient of that error depends on local derivatives beyond the curvature.

Definition: Let \(I_N=\int_{\mathbb{R}} g(x)\exp\{Nf(x)\}\,dx\), where \(N\) is an increasing information index. Suppose \(f\) has a unique interior maximum at \(x_0\), \(f''(x_0)<0\), and the integral is dominated by a neighbourhood of \(x_0\). The leading Laplace approximation is \(L_N=g(x_0)\exp\{Nf(x_0)\}\sqrt{2\pi/[N A]}\), where \(A=-f''(x_0)>0\). A relative error of order \(N^{-1}\) means \((I_N-L_N)/L_N=O(N^{-1})\).

This setup includes log-integrals arising in GLMMs. A cluster’s likelihood contribution and its random-effects density are combined in the log integrand. If the likelihood information increases with the number of observations in that cluster, the log integrand becomes increasingly concentrated. A prior or random-effects density that does not scale with the cluster size may still affect the mode and curvature; the large-information argument does not mean it can be discarded.

For the standard order statement, assume that the mode is isolated, the log integrand is smooth enough for a Taylor expansion through at least fourth order, and the tails contribute negligibly relative to the neighbourhood of the mode. A stronger expansion with a remainder of order \(N^{-2}\) requires further smoothness and suitable control of the tails. The dimension is held fixed when extending the rate to a multidimensional integral.

Deriving the One-Dimensional Correction

The leading Laplace approximation keeps only the quadratic term in the Taylor expansion. To find the first omitted contribution, set \(x=x_0+u/\sqrt{N}\). The \(N^{-1/2}\) scale is appropriate because the quadratic exponent is \(N f''(x_0)(x-x_0)^2/2\). Write \(f_k=f^{(k)}(x_0)\) and \(g_k=g^{(k)}(x_0)\). After removing the factor \(\exp\{Nf(x_0)\}/\sqrt{N}\), the exponent near the mode is

N\{f(x_0+u/\sqrt{N})-f(x_0)\} = -\frac{A u^2}{2} +\frac{f_3u^3}{6\sqrt{N}} +\frac{f_4u^4}{24N} +\cdots .

The amplitude expands as \(g(x_0+u/\sqrt{N})=g_0+g_1u/\sqrt{N}+g_2u^2/(2N)+\cdots\). Integrating against the leading Gaussian kernel, terms with an odd power of \(u\) vanish. In particular, the apparent \(N^{-1/2}\) contributions integrate to zero. At order \(N^{-1}\), the terms come from the quadratic term in the amplitude, the product of the linear amplitude term and the cubic exponent term, the quartic exponent term, and half the square of the cubic exponent term.

For a centred Gaussian kernel proportional to \(\exp(-Au^2/2)\), the relevant moments are \(E(u^2)=1/A\), \(E(u^4)=3/A^2\), and \(E(u^6)=15/A^3\). Substituting those moments gives the first correction:

I_N = L_N \left[ 1+\frac{1}{N} \left\{ \frac{g_2}{2g_0A} +\frac{g_1f_3}{2g_0A^2} +\frac{f_4}{8A^2} +\frac{5f_3^2}{24A^3} \right\} +O(N^{-2}) \right].

This expression assumes the stated smoothness and tail conditions are strong enough for the displayed remainder. If the integrand is simply \(\exp\{Nf(x)\}\), then \(g(x)=1\), so the first two terms in braces are zero. The correction then depends on the fourth derivative and on the square of the third derivative. The dimensions of the terms are consistent: for example, \(f_4/A^2\) and \(f_3^2/A^3\) are dimensionless.

Worked Example: A Non-Gaussian Integrand with a Computable Leading Error

Consider the hypothetical integral \(I_N=\int_{-1}^{1}\exp\{Nf(x)\}\,dx\), with \(f(x)=-x^2/2+0.6(\sin x-x)\). Its derivative is \(f'(x)=-x+0.6(\cos x-1)\). For \(0<x\leq1\), this is negative. For \(-1\leq x<0\), the inequality \(1-\cos x\leq x^2/2\) gives \(f'(x)\geq -x-0.3x^2>0\). Thus \(x_0=0\) is the unique maximum on the integration interval, and the endpoints have smaller values.

At the mode, \(f(0)=0\), \(f''(0)=-1\), \(f'''(0)=-0.6\), and \(f''''(0)=0\). Therefore \(A=1\), and the leading approximation is \(L_N=\sqrt{2\pi/N}\). Since \(g=1\), the relative correction coefficient is

\frac{f_4}{8A^2}+\frac{5f_3^2}{24A^3} = 0+\frac{5(0.6)^2}{24} =0.075.

At \(N=20\), the leading value is \(L_{20}=\sqrt{2\pi/20}\approx0.5605\). The first corrected value is \(0.5605(1+0.075/20)=0.5605(1.00375)\approx0.5626\). This is an asymptotic correction, not a claim that the corrected value is exact at \(N=20\). It illustrates that a nonzero third derivative contributes at order \(N^{-1}\) even when the fourth derivative at the mode is zero.

What the Error Order Means

If \(g(x_0)>0\) and the leading approximation is nonzero, an \(O(N^{-1})\) relative error corresponds to an absolute error of order \(L_N/N\). Since \(L_N\) itself contains a factor \(N^{-1/2}\) in one dimension, the absolute error in this setting is of order \(\exp\{Nf(x_0)\}N^{-3/2}\), up to fixed multiplicative factors. The logarithm has a related but distinct interpretation: if \(I_N=L_N(1+c/N+O(N^{-2}))\), then \(\log I_N-\log L_N=c/N+O(N^{-2})\). Thus the log-integral error is also of order \(N^{-1}\), not of order \(N^{-3/2}\).

For a fixed-dimensional integral over \(d\) random effects, the leading approximation has a factor \(N^{-d/2}\). Under analogous assumptions, the relative error is still typically \(O(N^{-1})\), so the absolute error is of order the leading integral times \(N^{-1}\). The coefficient involves contractions of higher-order derivative tensors with the inverse curvature matrix. A growing dimension or poorly conditioned curvature can invalidate this fixed-dimensional order statement.

In practice, a GLMM log integrand is often written directly as \(h(b)\), with mode \(\hat b\) and negative Hessian \(H=-\nabla^2h(\hat b)\). When the curvature grows with information and suitably scaled higher derivatives remain controlled, the leading approximation has the usual Gaussian determinant term. The first relative correction is governed by cubic and quartic terms in the local Taylor expansion. The scalar formula above makes their role explicit; the multidimensional version combines corresponding derivatives across coordinates.

Worked Example: Leading Correction for a Binomial Random-Intercept Integral

Suppose, hypothetically, one cluster has \(m=40\) binary observations, with 20 successes and 20 failures. Let its conditional success probability be \(\operatorname{logit}^{-1}(b)\), and let \(b\sim N(0,1)\). Ignoring no parameter-dependent terms, the log integrand for the marginal probability of the observed count is

h(b)=\log {40\choose20}+20b-40\log(1+e^b)-\frac{b^2}{2}-\frac12\log(2\pi).

The binomial count coefficient and prior normalizing constant are included so the integral is a marginal probability. By symmetry, \(h'(0)=20-40(1/2)=0\). The log likelihood has second derivative \(-40p(1-p)\), which at \(b=0\), where \(p=1/2\), equals \(-10\). The normal log density contributes another \(-1\), so \(h''(0)=-11\), and the mode has curvature \(A_h=11\). The likelihood’s third derivative is zero at \(b=0\), and the prior’s third derivative is also zero. The likelihood’s fourth derivative there is \(40/8=5\), while the prior’s fourth derivative is zero.

For a one-dimensional log integrand expanded at its mode, the first relative correction is \(h''''(\hat b)/(8A_h^2)+5\{h'''(\hat b)\}^2/(24A_h^3)\). Here this gives \(5/(8\cdot11^2)=5/968\approx0.005165\). The leading Laplace value is

L=\exp\{h(0)\}\sqrt{\frac{2\pi}{11}} =\frac{{40\choose20}}{2^{40}\sqrt{11}} \approx 0.0378.

Multiplying by the first correction gives \(L(1+5/968)\approx0.0380\). The correction is about \(0.5165\%\) of the leading value. This calculation uses the full log integrand, including the Gaussian random-effects density, when finding the mode and curvature. The result is an asymptotic expansion in the information contributed by the cluster observations; with only one hypothetical cluster size, its numerical accuracy is not guaranteed by the order statement alone.

Worked Example: Fixed-Dimensional Error Order in Two Coordinates

Consider a hypothetical separable integral on \([-1,1]^2\) with \(f(x,y)=-(x^2+y^2)/2+0.6(\sin x-x)+0.3(\sin y-y)\), and integrand \(\exp\{Nf(x,y)\}\). The same derivative argument as in the one-dimensional example shows that the unique maximum is \((0,0)\). The negative Hessian there is the identity matrix, so the leading approximation is \(L_N=2\pi/N\).

Because the function is separable, the integral is a product of two one-dimensional integrals. The first coordinate contributes a relative correction \(0.075/N\). The second has third derivative \(-0.3\), fourth derivative zero, and therefore contributes \(5(0.3)^2/(24N)=0.01875/N\). Multiplying the two expansions gives

\frac{I_N}{L_N} = \left(1+\frac{0.075}{N}+O(N^{-2})\right) \left(1+\frac{0.01875}{N}+O(N^{-2})\right) = 1+\frac{0.09375}{N}+O(N^{-2}).

The product of the two order-\(N^{-1}\) corrections is of order \(N^{-2}\), which is why it does not enter the leading coefficient. This example illustrates the fixed-dimension result and shows why a higher-dimensional coefficient must account for derivatives in all coordinates, rather than inspecting only one marginal direction.

Checking the Approximation and Choosing an Alternative

An asymptotic order is a statement about a sequence of integrals as information increases. It does not provide a reliable error bound for a particular fitted model. Useful checks include verifying that the mode is interior, inspecting the eigenvalues of the negative Hessian, and confirming that the log integrand is smooth and well behaved in the region contributing to the integral. A very small curvature eigenvalue indicates a broad direction in which the local concentration argument may be weak. Large scaled third- or fourth-order derivatives indicate that the quadratic approximation may leave an appreciable correction at the available information level.

When feasible, compare the Laplace value with a more accurate numerical integration calculation for representative clusters or parameter values. Agreement as the numerical integration is refined supports the numerical calculation; it does not prove that the Laplace approximation has met a specified error tolerance. Disagreement is a signal to investigate the local shape, curvature, and integration settings rather than to assume that the asymptotic rate applies at the observed cluster size.

Approach	What it provides	Important limitation
Laplace approximation	A local Gaussian approximation based on the mode and curvature, with low computational cost.	Its usual error order requires increasing information and regularity; the leading approximation alone omits higher-derivative corrections.
Exact Gaussian integration	An exact integral when the log integrand is quadratic over the integration domain.	It does not extend as an exact calculation to a general nonquadratic GLMM integrand.
Numerical quadrature	A direct numerical approximation to the integral that can be used for comparison.	Accuracy and computational cost depend on the numerical method and, in multiple dimensions, the integration problem can become demanding.

The Gaussian case is an important limit: if the log integrand is exactly quadratic and the domain and normalization are handled correctly, there are no higher-order terms and Laplace integration is exact. Away from that case, the leading error is not determined by curvature alone. In GLMMs, information may grow through more observations within a cluster, but simply increasing the number of clusters does not automatically make each cluster-specific random-effects integral more concentrated. The relevant information must grow for the stated within-integral asymptotic argument to apply.

Common Pitfall

A common mistake is to read “error of order \(N^{-1}\)” as “the error is exactly \(1/N\),” or to conclude that any dataset with a large total number of observations must yield an accurate Laplace approximation. The order notation describes how the relative error behaves in a suitable sequence of integrals; it does not give its coefficient or certify a finite-sample tolerance. The coefficient depends on local derivatives, and the information relevant to a random-effects integral is the information concentrating that particular integral. Check the cluster-level curvature and higher derivatives, and use a numerical comparison when accuracy is consequential.

Key takeaway: Under regularity and concentration conditions, the leading Laplace approximation has relative error of order \(N^{-1}\). The first correction reflects third- and fourth-order local shape, while curvature determines the Gaussian scale. The asymptotic rate is informative, but it is not a finite-sample error guarantee.

Check Your Understanding

Answer the questions, then check each result in the solutions.

For \(I_N=\int \exp\{Nf(x)\}\,dx\), why does the local variable scale as \(N^{-1/2}\) when \(f''(x_0)<0\)?
If \(g=1\), \(A=2\), \(f_3=0\), and \(f_4=-4\), what is the coefficient of \(N^{-1}\) in the relative correction?
For a one-dimensional integral, the leading approximation is \(L_N\propto\exp\{Nf(x_0)\}N^{-1/2}\). If its relative error is \(O(N^{-1})\), what is the order of its absolute error?
Why does increasing the number of clusters alone not necessarily justify the large-information error order for each random-effects integral?

Solutions. (1) The quadratic exponent is approximately \(Nf''(x_0)(x-x_0)^2/2\). Keeping this term of constant order requires \(x-x_0=u/\sqrt{N}\). (2) The coefficient is \(f_4/(8A^2)+5f_3^2/(24A^3)=-4/(8\cdot4)+0=-0.125\), so the first relative correction is \(-0.125/N\). (3) The absolute error is the leading scale times \(N^{-1}\), hence \(O(\exp\{Nf(x_0)\}N^{-3/2})\), up to fixed factors. (4) Each integral is associated with a cluster’s random effects; if the information within each cluster remains fixed, that integral need not become increasingly concentrated even when the number of separate clusters grows.