Comparing Laplace and Adaptive Quadrature in a Binary Example

Why Compare the Two Approximations?

A GLMM marginal likelihood integrates the conditional response likelihood over the random-effects distribution. For binary data, the integrand combines a logistic likelihood with a Gaussian density and is generally not itself Gaussian. The Laplace approximation uses the integrand’s value and curvature at its mode; adaptive Gauss–Hermite quadrature uses a local Gaussian approximation to place its nodes, then evaluates the integrand at those nodes. Comparing both methods on a small, fully specified example makes the role of the mode, curvature, and quadrature weights explicit.

Throughout, consider one cluster with binary responses $Y_j$, a random intercept $b\sim N(0,\sigma^2)$, and fixed intercept $\beta$. Conditional on $b$, responses are independent Bernoulli variables with success probability $p(b)=\operatorname{logit}^{-1}(\beta+b)$. The examples specify $\beta$, $\sigma$, and the observed responses; these are numerical integration exercises, not parameter-estimation exercises.

Definition: For observed binary responses $y_1,\ldots,y_m$, the cluster marginal likelihood is $$ L=\int_{-\infty}^{\infty}\left\{\prod_{j=1}^{m}p(b)^{y_j}[1-p(b)]^{1-y_j}\right\} \frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{b^2}{2\sigma^2}\right)\,db. $$ The log-integrand is $h(b)$, the natural logarithm of the entire integrand, including the normal-density constant. The Laplace approximation at its mode $\hat b$, when $h''(\hat b)<0$, is $$ L_{\mathrm{Lap}}=\exp\{h(\hat b)\}\sqrt{\frac{2\pi}{-h''(\hat b)}}. $$

The assumptions here are a correctly specified conditional Bernoulli-logit model, conditional independence within the cluster given $b$, and the stated normal distribution for $b$. The integrand is positive and smooth. For adaptive quadrature, we also need an identifiable local maximum with positive curvature quantity $A=-h''(\hat b)$. These conditions make the calculations well-defined; they do not guarantee that a low-order quadrature rule is accurate.

Worked Example: Two Successes in One Cluster

Worked Example: Laplace and Five-Node Adaptive Quadrature

Suppose one cluster has $y_1=y_2=1$, $\beta=0$, and $\sigma=1$. Thus $p(b)=1/(1+e^{-b})$, and the likelihood conditional on $b$ is $p(b)^2$. The log-integrand, including the normalizing constant, is $$ h(b)=2\log p(b)-\frac{b^2}{2}-\frac12\log(2\pi). $$

Find the mode. Differentiation gives $h'(b)=2\{1-p(b)\}-b$. Solving $h'(b)=0$ numerically gives $\hat b\approx0.6748$, where $p(\hat b)\approx0.6626$. The second derivative is $h''(b)=-2p(b)\{1-p(b)\}-1$, so $$ A=-h''(\hat b)=1+2(0.6626)(0.3374)\approx1.447. $$ The curvature is positive, confirming a local maximum. Since $h'(b)$ decreases strictly—the second derivative is negative for every $b$—this is also the unique mode.

Laplace approximation. Substitution, including the normal constant, gives $$ h(\hat b)=2\log(0.6626)-\frac{0.6748^2}{2}-\frac12\log(2\pi)\approx-1.970. $$ Therefore $$ L_{\mathrm{Lap}}=\exp(-1.970)\sqrt{\frac{2\pi}{1.447}}\approx0.2906. $$ This approximates the probability of observing two successes in this cluster after averaging over its random intercept.

Five-node adaptive quadrature. For standard-normal-weight Gauss–Hermite nodes $x_i$ and weights $w_i$, the adaptive rule is $$ L_{\mathrm{AGHQ}}\approx \sqrt{\frac{2}{A}}\exp\{h(\hat b)\} \sum_i w_i\exp\left\{x_i^2+h\left(\hat b+\sqrt{\frac{2}{A}}x_i\right)-h(\hat b)\right\}. $$ The five nodes are $0,\ \pm0.95857,\ \pm2.02018$, with corresponding weights $0.945309,\ 0.393619,\ 0.0199532$ (paired symmetrically). Here $\sqrt{2/A}\approx1.1755$. The correction factors $\exp\{x_i^2+h(b_i)-h(\hat b)\}$, from the negative outer node to the positive outer node, are approximately $0.952,\ 0.979,\ 1.000,\ 1.042,\ 1.475$. In particular, the positive outer-node factor is $1.475$.

The weighted sum is $$ 0.0199532(0.952+1.475)+0.393619(0.979+1.042)+0.945309(1.000)\approx1.789. $$ Thus $L_{\mathrm{AGHQ}}\approx \exp(-1.970)\sqrt{2/1.447}(1.789)\approx0.2934$. It is about $0.0028$ higher than the Laplace value. Agreement at this level is a result for this example and rule, not a general accuracy guarantee.

What the Adaptive Rule Changes

The quadrature formula follows by locally matching the integrand to a Gaussian. Around the mode, write $b=\hat b+\sqrt{2/A}\,x$. The quadratic Taylor approximation is $h(b)\approx h(\hat b)-x^2$. Consequently, the integrand can be rewritten as $\exp\{h(\hat b)\}\exp(-x^2)\exp\{x^2+h(b)-h(\hat b)\}$. Integrating in $b$ contributes the Jacobian $\sqrt{2/A}$; applying Gauss–Hermite quadrature to the Gaussian-weighted integral gives the displayed adaptive formula.

If the correction factor were exactly one at every node, the rule would reduce to the Gaussian integral, and the Laplace and adaptive calculations would coincide. In a non-Gaussian integrand, departures of the correction factors from one record how the true log-integrand differs from its quadratic approximation at the selected nodes. This is why adaptive quadrature can improve on Laplace, especially when the integrand is asymmetric or its tails matter.

Worked Example: A Balanced Binary Cluster

Worked Example: Curvature and Five-Node Quadrature

Now let $y_1=1,\ y_2=0$, with $\beta=0$ and $\sigma=1$. The conditional likelihood is $p(b)\{1-p(b)\}$, and $$ h(b)=\log p(b)+\log\{1-p(b)\}-\frac{b^2}{2}-\frac12\log(2\pi). $$

At $b=0$, $p(0)=0.5$, and $h'(0)=0$. Also $h''(b)=-2p(b)\{1-p(b)\}-1$, so the mode is unique and $$ A=-h''(0)=2(0.5)(0.5)+1=1.5. $$ The mode value is $h(0)=\log(0.25)-\tfrac12\log(2\pi)$. The Laplace result simplifies to $$ L_{\mathrm{Lap}}=\frac{0.25}{\sqrt{1.5}}\approx0.2041. $$

For five-node adaptive quadrature, the scaling is $\sqrt{2/1.5}\approx1.1547$. The correction factor at node $x$ is $$ \exp\left\{\frac{x^2}{3}\right\}\frac{4p(1.1547x)\{1-p(1.1547x)\}}{1}. $$ At $x=0,\ \pm0.95857,\ \pm2.02018$, the factors are approximately $1.000,\ 1.015,\ 1.257$, respectively, with each nonzero value applying to both signs. The weighted sum is approximately $$ 0.945309(1.000)+0.393619(2.029)+0.0199532(2.514)\approx1.794. $$ It follows that $L_{\mathrm{AGHQ}}\approx0.1152(1.794)\approx0.2066$, where $0.1152=\sqrt{2/1.5}\,(0.25)/\sqrt{2\pi}$. This is slightly above $0.2041$, the Laplace approximation. The balanced outcomes make the mode particularly transparent, but the integrand is still not exactly Gaussian.

Worked Example: Checking Quadrature Order

Worked Example: Three Nodes Versus Five Nodes

Return to the two-success cluster with $\hat b\approx0.6748$, $A\approx1.447$, and $L_{\mathrm{Lap}}\approx0.2906$. A three-node rule has nodes $0,\ \pm1.224745$ and weights $1.181636,\ 0.295409,\ 0.295409$. Transforming the nonzero nodes gives random-effect values $b\approx2.115$ and $b\approx-0.7650$. Evaluating the same adaptive correction factor gives approximately $1.091$ and $0.965$, respectively; the center factor is $1$.

The weighted sum is $1.181636+0.295409(1.091+0.965)\approx1.789$. Hence the three-node estimate is $\exp(-1.970)\sqrt{2/1.447}(1.789)\approx0.2933$. The five-node estimate was about $0.2934$. Their closeness is a useful numerical check for this particular integral. It does not establish that either value is exact, nor that the same small rule will work for another dataset or parameter choice.

How to Assess the Comparison

A sound comparison checks the model and the numerical calculation separately. First verify the conditional likelihood and the random-effect density: omitting the Gaussian normalization constant changes $h(\hat b)$, even though it does not change the mode or curvature. Then verify the sign and scale of the curvature, the transformed nodes, and the quadrature weights. Finally, increase quadrature order and check whether the marginal-likelihood estimate stabilizes. For sums with very small or large terms, log-sum-exp evaluation avoids numerical underflow or overflow.

The methods answer the same integration problem but use different approximations. Laplace requires a mode and its local curvature and is inexpensive; adaptive Gauss–Hermite quadrature also uses the mode and curvature but evaluates the integrand at multiple nodes. Nonadaptive Gauss–Hermite quadrature instead places nodes around zero in standard-normal coordinates, which can be inefficient when the integrand’s important mass is far from zero. Direct numerical integration can be a useful independent check in one dimension, though it is not the computational strategy for every GLMM.

Diagnostics include checking that the mode is a maximum, looking for a broad or skewed integrand, and comparing successive quadrature orders. Multiple peaks, a very flat mode, or substantial changes as the order increases call for caution: a single local Gaussian approximation may not represent the integral well. Even when two approximations agree, they can share error. These calculations also concern a marginal likelihood at fixed parameters; they do not by themselves establish that the model is appropriate or that fitted parameters are accurate.

Common pitfall: A frequent error is to omit $-\tfrac12\log(2\pi)$ from the log-integrand but then use the full Laplace formula as if that constant were included. The resulting integral is scaled incorrectly. Define $h(b)$ as the log of the entire integrand, including all density constants, and use that same $h$ in both Laplace and adaptive calculations.

Check Your Understanding

Use the stated model and calculations to answer the questions. Short solutions follow.

For the two-success example, why does including the normal-density constant change the Laplace approximation but not the mode?
In the balanced example, what is the curvature quantity $A$, and why is it positive?
For the two-success example, what does the five-node weighted sum approximate in the adaptive formula?
Why is agreement between three- and five-node results evidence of numerical stability but not proof of exactness?

Solutions. (1) The constant adds the same value to $h(b)$ for every $b$, so it does not change the maximizing $b$. It multiplies $\exp\{h(\hat b)\}$, and therefore the approximated integral, by the corresponding density-normalization factor. (2) $A=1+2(0.5)(0.5)=1.5>0$; the second derivative at the mode is $-1.5$, confirming a local maximum. (3) It approximates the Gaussian-weighted integral of the correction factor $\exp\{x^2+h(b_i)-h(\hat b)\}$, with the Gauss–Hermite weights. The prefactor $\exp\{h(\hat b)\}\sqrt{2/A}$ then gives the marginal-likelihood estimate. (4) Both rules are finite numerical approximations, and their errors need not be the same or small. Stability should be assessed by increasing order and, where practical, checking against an independent integration method.