Skewness and Sparse-Cluster Error in Laplace Approximation

Why Local Shape Matters

Laplace approximation replaces an integrand near its dominant point by a quadratic, or equivalently replaces its logarithm by a parabola near its mode. This is exact when the log integrand is quadratic. For binary-response GLMMs, however, the conditional log likelihood is generally not quadratic in the random effect. Its local asymmetry and higher-order curvature can matter, especially when a cluster contains little information about its random effect.

The issue is not simply whether a cluster is “small.” The conditional mode, the curvature there, the outcome pattern, and the prior all influence how well a quadratic captures the integrand. Derivatives at the mode provide a useful diagnostic, but they cannot by themselves certify the approximation across the whole integration region.

Binary Random-Intercept Integrand and Local Shape

Consider a binary-response random-intercept model with a logit link. For cluster \(j\), let \(Y_{ij}\) be the binary response for observation \(i\), and let \(b_j\sim N(0,\sigma_b^2)\). Conditional on \(b_j\), suppose the responses are independent with \(P(Y_{ij}=1\mid b_j)=p_{ij}(b_j)\), where \(\operatorname{logit}\{p_{ij}(b_j)\}=x_{ij}^{T}\beta+b_j\). The cluster’s marginal likelihood contribution is an integral of the conditional likelihood against the normal density of \(b_j\).

Write its integrand as \(\exp\{h_j(b)\}\), including the random-effects density in \(h_j\). Constants in \(h_j\) do not affect its mode or derivatives, although they do matter for the integral’s overall value. Let \(\hat b_j\) be an interior conditional mode and define \(J_j=-h_j''(\hat b_j)>0\). The ordinary one-dimensional Laplace approximation is \(I_{Lj}=\exp\{h_j(\hat b_j)\}\sqrt{2\pi/J_j}\).

Definition: The local skewness diagnostic is the standardized third derivative \(h_j'''(\hat b_j)/J_j^{3/2}\). A useful fourth-derivative diagnostic is \(h_j''''(\hat b_j)/J_j^2\). These dimensionless quantities describe local departures from a quadratic log integrand; they are not, in general, the skewness and kurtosis of a probability distribution.

For a binary random-intercept model, put \(p_{ij}=p_{ij}(\hat b_j)\) and \(q_{ij}=1-p_{ij}\). Differentiating the log integrand gives the following mode curvature and higher derivatives:

J_j=\frac{1}{\sigma_b^2}+\sum_i p_{ij}q_{ij},\qquad h_j'''(\hat b_j)=-\sum_i p_{ij}q_{ij}(1-2p_{ij}),\qquad h_j''''(\hat b_j)=-\sum_i p_{ij}q_{ij}(1-6p_{ij}q_{ij}).

The Gaussian prior contributes precision \(1/\sigma_b^2\) to \(J_j\), but contributes zero to the third and fourth derivatives. Each binary observation contributes at most \(1/4\) to the likelihood part of the curvature, since \(p q\leq 1/4\). Thus a small cluster may supply limited likelihood curvature relative to the prior, while a highly unbalanced response pattern can also make the local shape asymmetric.

Where the First Local Correction Comes From

The next calculation explains how the third and fourth derivatives enter. Suppress the cluster index and write \(u=b-\hat b\). At an interior mode, \(h'(\hat b)=0\), so a local Taylor expansion is

h(\hat b+u)=h(\hat b)-\frac{J u^2}{2} +\frac{h'''(\hat b)u^3}{6} +\frac{h''''(\hat b)u^4}{24}+\cdots.

Under the quadratic approximation, \(u\) has the Gaussian kernel proportional to \(\exp(-Ju^2/2)\). Its normalized moments are \(E(u^4)=3/J^2\) and \(E(u^6)=15/J^3\). Expanding the exponential of the cubic and quartic terms, the quartic term contributes \(h''''(\hat b)E(u^4)/24\), while the square of the cubic term contributes \(h'''(\hat b)^2E(u^6)/72\). Consequently, the first local correction is

C=\frac{h''''(\hat b)}{8J^2} +\frac{5\{h'''(\hat b)\}^2}{24J^3},\qquad \int e^{h(b)}\,db\ \approx\ I_L(1+C).

This is a local expansion, not an error bound. It assumes a smooth log integrand, a well-defined dominant interior mode with \(J>0\), and enough concentration near that mode for the Taylor expansion to represent the important part of the integral. The omitted terms can matter, particularly if the integrand is broad or its shape changes substantially away from the mode. The cubic contribution to \(C\) is nonnegative, but the fourth-derivative contribution may have either sign.

Worked Example: One Success in a Sparse Cluster

Take one hypothetical binary observation \(y=1\) in a random-intercept model with fixed linear predictor zero and \(\sigma_b^2=1\). The log integrand, omitting constants, is \(h(b)=b-\log(1+e^b)-b^2/2\). Its score is \(h'(b)=1-p(b)-b\), so the mode solves \(q(b)=b\). Numerically, \(\hat b\approx0.4011\), \(p\approx0.5989\), and \(q\approx0.4011\). These rounded values satisfy the score equation to the displayed precision.

At the mode, \(p q\approx0.2402\), so the curvature is \(J=1+pq\approx1.2402\). The third and fourth derivatives are approximately \(h'''(\hat b)=0.04753\) and \(h''''(\hat b)=0.1060\). Substitution into the first correction gives

C\approx\frac{0.1060}{8(1.2402)^2} +\frac{5(0.04753)^2}{24(1.2402)^3} \approx0.00861+0.000247=0.00886.

The local expansion therefore predicts an integral about \(0.886\%\) above the ordinary Laplace value, since \(1+C\approx1.00886\). This is a diagnostic estimate, not a claim that the exact relative error is \(0.886\%\). With only one response, the cluster provides very little likelihood information about its random intercept; a small correction at the mode should not be mistaken for proof that the full integrand is well approximated.

Worked Example: Two Successes and Stronger Local Asymmetry

Now suppose a hypothetical cluster has two observations, both equal to one, again with fixed linear predictor zero and \(\sigma_b^2=1\). Here \(h'(b)=2q(b)-b\), and the mode satisfies \(b=2q(b)\). Numerically, \(\hat b\approx0.6748\), with \(p\approx0.6626\) and \(q\approx0.3374\); thus \(2q\approx0.6748\), consistent with the mode equation.

Using \(p q\approx0.22356\), the curvature is \(J=1+2pq\approx1.4471\). The derivative formulas give \(h'''(\hat b)\approx0.1454\) and \(h''''(\hat b)\approx0.1526\). Therefore

C\approx\frac{0.1526}{8(1.4471)^2} +\frac{5(0.1454)^2}{24(1.4471)^3} \approx0.00911+0.00145=0.01056.

The first local correction is about \(1.056\%\), so this expansion predicts an integral approximately \(1.01056\) times the Laplace value. The positive third derivative reflects asymmetry around the mode, and its squared contribution increases the correction. The example illustrates why outcome pattern and local derivatives matter alongside the cluster size; it does not establish that the two-observation approximation is accurate to exactly \(1.056\%\).

Worked Example: A Balanced Cluster Can Still Have Higher-Order Curvature

Consider a hypothetical cluster of 20 observations, with 10 successes and 10 failures, fixed predictor zero, and \(\sigma_b^2=1\). By symmetry, the score is zero at \(b=0\), so \(\hat b=0\) and \(p=q=0.5\). The curvature is \(J=1+20(0.25)=6\). The third derivative is zero because each term contains \(1-2p=0\). But the fourth derivative is not zero:

h''''(0)=-20(0.25)\{1-6(0.25)\}=2.5,\qquad C=\frac{2.5}{8(6)^2}+0\approx0.008681.

Here the local cubic asymmetry vanishes, but the quartic term remains. The first correction is about \(0.8681\%\). A zero third derivative alone is therefore not evidence that the log integrand is exactly quadratic or that Laplace approximation is exact.

Reading the Diagnostics and Choosing an Approach

For practical diagnosis, inspect both the standardized derivative terms and the combined correction. A substantial value of \(h'''(\hat b)/J^{3/2}\) signals local asymmetry; a substantial value of \(h''''(\hat b)/J^2\) signals local quartic departure. The correction \(C\) combines these contributions, but cancellation is possible: a small net \(C\) can conceal individually sizeable terms. It is therefore helpful to examine the terms separately as well as their sum.

Sparse clusters deserve attention because their likelihood curvature is limited and their conditional log integrands may not be sharply concentrated. If \(J\) is small, the Gaussian approximation is broad, so behavior farther from the mode can matter. Conversely, many observations do not guarantee a good approximation if the mode lies in a region of pronounced asymmetry or if higher derivatives remain important. These are warnings to check the integral, not automatic rules that classify every small cluster as a failure.

The usual Laplace approximation is computationally economical and often useful when the integrand is well represented by its local quadratic shape. A higher-order local expansion adds derivative information, but its correction is still local and does not remove the need to assess omitted terms. Numerical integration, including sufficiently accurate quadrature, provides a direct way to check a questionable one-dimensional integral. In routine model fitting, such checks can be made on representative clusters or parameter settings. The Gauss-Hermite and adaptive quadrature methods provide systematic alternatives for this purpose, though their numerical accuracy also needs checking.

Key takeaway: In a binary-response GLMM, skewness and sparse within-cluster information can make the local quadratic approximation questionable. Derivatives at the conditional mode quantify local departures and supply a first correction, but only comparison with the integrand away from the mode can reveal whether that local picture captures the integral adequately.

Common Pitfalls

Calling a small correction a guarantee. The reasoning “\(C\) is small, therefore Laplace is accurate” ignores higher-order terms and behavior away from the mode. Treat \(C\) as a local warning or reassurance, then verify when accuracy matters.
Checking only skewness. A symmetric integrand can have a nonzero fourth derivative, as the balanced-cluster example shows. Inspect both third- and fourth-derivative contributions.
Assuming more clusters fix each cluster’s integral. Adding independent clusters can add information about population parameters, but it does not necessarily make each cluster-specific integrand more Gaussian.
Interpreting the correction with the wrong sign. If \(C>0\), the first corrected expansion is \(I_L(1+C)\), above the ordinary Laplace value. This is a prediction from the truncated expansion, not a guaranteed direction for the exact error.

Check Your Understanding

Answer each question, then compare with the solutions below.

For a binary random-intercept model, which part of \(J_j\) comes from the normal random-effects density, and what is its value?
Why does a zero third derivative at the mode not establish exactness of Laplace approximation?
If \(h''''(\hat b)=0.08\), \(J=2\), and \(h'''(\hat b)=0\), calculate the first correction \(C\).
What does a positive first correction \(C=0.02\) say, and what does it not say?

Solutions. (1) The Gaussian density contributes its precision, \(1/\sigma_b^2\), to \(J_j\). (2) The fourth and higher derivatives can still be nonzero, and the integrand away from the mode may depart from its quadratic approximation. (3) \(C=0.08/\{8(2)^2\}+0=0.08/32=0.0025\), or \(0.25\%\); the first corrected integral is \(I_L(1.0025)\). (4) The truncated local expansion predicts an integral about \(2\%\) above \(I_L\). It does not guarantee the exact integral is 2% higher, because omitted terms and nonlocal shape can affect the result.