Quadrature Accuracy for Rare Binary Responses

Why Rare Responses Challenge Local Approximation

For a binary generalized linear mixed model, integrating out a random intercept means averaging a cluster’s conditional likelihood over the random-effect distribution. When responses are rare, many clusters contain no events, and the conditional likelihood can change sharply in one direction while flattening in another. The resulting integrand may be asymmetric even though the random-effect distribution is Gaussian.

Laplace approximation replaces the log-integrand near its mode by a quadratic function. Adaptive Gauss-Hermite quadrature also centers and scales its nodes using the mode and local curvature, but evaluates the original integrand at those nodes. The distinction matters: a good local Gaussian fit supports both methods, whereas substantial asymmetry can make a single quadratic approximation inaccurate.

Model and Log-Integrand

Consider one cluster with $m$ binary responses $Y_i$, a random intercept $b$, and a fixed intercept $\beta_0$. Conditional on $b$, suppose the responses are independent Bernoulli variables with probability $p(b)$. Assume $b\sim N(0,\sigma_b^2)$, with $\sigma_b^2>0$, and use the logit link.

Definition: The binary random-intercept log-integrand, with data-only constants omitted, is $$ h(b)=\sum_{i=1}^{m}\left\{y_i(\beta_0+b)-\log\left[1+\exp(\beta_0+b)\right]\right\} -\frac{b^2}{2\sigma_b^2}, \qquad p(b)=\frac{\exp(\beta_0+b)}{1+\exp(\beta_0+b)}. $$ The marginal likelihood contribution is proportional to the integral of $\exp\{h(b)\}$ over $b$; with the Gaussian density normalized, it is $$ I=\frac{1}{\sqrt{2\pi}\sigma_b}\int_{-\infty}^{\infty}\exp\{h(b)\}\,db. $$

Here the data are treated as conditionally independent given the random intercept, and the normal distribution for $b$ is assumed to be the specified mixing distribution. These assumptions define the integration problem; accurate numerical integration cannot correct a misspecified response or random-effects model. Write $y=\sum_i y_i$ for the number of events. Differentiation gives

h'(b)=y-mp(b)-\frac{b}{\sigma_b^2}, \qquad h''(b)=-mp(b)\{1-p(b)\}-\frac{1}{\sigma_b^2}, \qquad h'''(b)=-mp(b)\{1-p(b)\}\{1-2p(b)\}.

Since $h''(b)<0$ everywhere, $h$ is strictly concave and has at most one mode. Its tails are dominated by the negative quadratic normal penalty, so a finite mode exists. For rare probabilities, $p(b)$ is small around relevant values, and $h'''(b)$ is typically negative. A negative third derivative means that, locally, the log-integrand falls more on the positive side of the mode than on the equally distant negative side. This is a warning about asymmetry, not by itself a numerical error measure.

How Laplace and Adaptive Quadrature Use the Integrand

Let $\hat b$ denote the mode and define $H=-h''(\hat b)>0$, the negative curvature there. A second-order Taylor expansion around the mode is

h(b)\approx h(\hat b)-\frac{H}{2}(b-\hat b)^2.

Integrating this Gaussian approximation, including the normalizing constant of the random-effect density, yields the one-dimensional Laplace approximation:

I_{\mathrm{Lap}}\approx \frac{\exp\{h(\hat b)\}}{\sigma_b\sqrt{H}}.

This approximation uses the height and curvature at one point; it does not use the rest of the integrand’s shape. Adaptive Gauss-Hermite quadrature instead places nodes around $\hat b$, scaled by $H^{-1/2}$, and evaluates the actual log-integrand at them. For physicists’ Gauss-Hermite nodes $x_k$ and weights $w_k$, a $K$-node form is

I_{\mathrm{AGH}}\approx \frac{\exp\{h(\hat b)\}}{\sigma_b\sqrt{\pi H}} \sum_{k=1}^{K} w_k \exp\left\{h\left(\hat b+\sqrt{\frac{2}{H}}x_k\right)-h(\hat b)+x_k^2\right\}.

If the log-integrand were exactly quadratic, the exponential correction inside the sum would equal one, and the formula would reduce to the Laplace result because the weights sum to $\sqrt{\pi}$. For a nonquadratic integrand, adaptive quadrature retains information that Laplace discards. It is not automatically exact at a low node count; accuracy must be checked by increasing the number of nodes and comparing the resulting integrals.

Worked Example: Ten Rare-Outcome Non-events

Suppose a hypothetical cluster has $m=10$, $y=0$, $\beta_0=-4$, and $\sigma_b=1$. Its conditional event probability at $b=0$ is $\operatorname{logit}^{-1}(-4)\approx0.01799$. For this all-zero cluster, the score equation is $h'(b)=-10p(b)-b=0$.

Solving numerically gives $\hat b\approx-0.15451$, with $p(\hat b)\approx0.01545$. Substituting into the score gives $h'(\hat b)\approx-10(0.01545)-(-0.15451)\approx0.00001$, which is approximately zero at the displayed precision. The conditional mode is negative: shifting the random intercept down makes an already rare event less likely, consistent with observing no events.

The negative curvature is $H=10p(\hat b)\{1-p(\hat b)\}+1\approx10(0.01545)(0.98455)+1\approx1.152$. The log-integrand at the mode is $h(\hat b)=10\log\{1-p(\hat b)\}-\hat b^2/2\approx10\log(0.98455)-(-0.15451)^2/2\approx-0.1676$. Thus $I_{\mathrm{Lap}}\approx\exp(-0.1676)/\sqrt{1.152}\approx0.788$. This is an approximation to the cluster’s marginal likelihood contribution, with the Bernoulli sequence’s data-only constants omitted.

The score and curvature locate and scale the local Gaussian approximation, but they do not reveal how well a quadratic represents the shoulders and tails. Adaptive quadrature can check those regions by evaluating the full integrand. The example’s $h'''(\hat b)$ is approximately $-10(0.01545)(0.98455)(0.9691)\approx-0.1475$, so the log-integrand is not exactly symmetric around its mode.

Worked Example: One Event Among Ten Rare Responses

Now suppose a hypothetical cluster has the same $m=10$, $\beta_0=-4$, and $\sigma_b=1$, but exactly one response is an event. The score is $h'(b)=1-10p(b)-b$. A numerical root is $\hat b\approx0.65832$, where $p(\hat b)\approx0.03417$. Substitution gives $1-10(0.03417)-0.65832\approx-0.00002$, consistent with a zero score to the displayed precision.

The single event shifts the mode above zero, from the negative mode in the previous example to a positive one. The curvature is $H=1+10(0.03417)(0.96583)\approx1.330$. For a specified sequence with one event, the log-likelihood is $(\beta_0+b)-10\log\{1+\exp(\beta_0+b)\}$. At the mode, this gives $h(\hat b)\approx-3.906$, including the normal penalty $-\hat b^2/2$. Therefore $I_{\mathrm{Lap}}\approx\exp(-3.906)/\sqrt{1.330}\approx0.0174$.

If the observations are represented only by the count $y=1$, the binomial likelihood includes the constant $\binom{10}{1}=10$. It does not change the mode or curvature; it multiplies the integral by 10. Keeping track of that convention is essential when comparing marginal likelihood values.

Here $h'''(\hat b)\approx-10(0.03417)(0.96583)(0.93166)\approx-0.307$. The event changes the mode and the local shape, while the negative third derivative still signals local asymmetry. A Laplace value alone cannot tell whether that asymmetry materially affects the integral; adaptive quadrature or a sufficiently accurate reference calculation is needed for that assessment.

A Rare-Event Example with a Convenient Exact Mode

Worked Example: Checking Curvature and Asymmetry at the Mode

Consider a hypothetical cluster with $m=20$, $y=1$, $\sigma_b=1$, and $\beta_0=\log(0.05/0.95)\approx-2.94444$. At $b=0$, $p(b)=0.05$, so the score is $h'(0)=1-20(0.05)-0=0$. Strict concavity therefore makes $\hat b=0$ the unique mode.

The curvature is $H=1+20(0.05)(0.95)=1.95$, and the third derivative is $h'''(0)=-20(0.05)(0.95)(0.90)=-0.855$. The log-integrand at the mode is $h(0)=\beta_0-20\log(1+\exp(\beta_0))\approx-2.94444-20\{ -\log(0.95)\}\approx-3.97031$. Consequently, $I_{\mathrm{Lap}}\approx\exp(-3.97031)/\sqrt{1.95}\approx0.0135$.

Although the mode is exactly at zero, the nonzero third derivative shows that the log-integrand is not locally symmetric. A centered Gaussian approximation is not made exact simply by centering it correctly. This is a useful distinction: mode location and curvature determine where adaptive nodes go, but the correction terms from evaluating the full integrand determine whether quadrature improves on the local approximation.

Assessing Accuracy and Recognizing Failure Modes

A rare event rate does not imply that Laplace approximation must fail. The relevant question is whether the log-integrand is close to quadratic over the region that contributes appreciably to its integral. Cluster size, the event pattern, $\beta_0$, and $\sigma_b^2$ all affect that shape. In particular, an all-zero cluster and a cluster with one event can have different mode locations even under the same fixed effects and random-effect variance.

A practical accuracy assessment compares the Laplace approximation with adaptive quadrature at increasing node counts, using the same likelihood convention and model parameters. The sequence of quadrature values should stabilize; a large and persistent difference from Laplace indicates that its local Gaussian approximation may be inadequate for that cluster. Also check that the reported mode actually solves $h'(b)=0$ and that $H=-h''(\hat b)$ is positive. In this model strict concavity guarantees a unique mode, so a purported second mode or nonpositive curvature points to a calculation or implementation problem.

The main alternatives answer different needs. Laplace approximation is computationally inexpensive and uses one mode and curvature. Adaptive quadrature costs more but can represent nonquadratic shape around the mode, with accuracy depending on the number of nodes. Nonadaptive Gauss-Hermite quadrature uses fixed nodes rather than centering them on the cluster’s conditional mode; when rare observations shift that mode from the prior center, fixed nodes may represent the important region poorly unless the rule is sufficiently rich. None of these integration methods addresses misspecification of the random-effects distribution or conditional-independence assumption.

Numerical integration error is distinct from sampling uncertainty in estimated model parameters. An integral that changes with quadrature order is a numerical warning; a stable integral does not establish that the model fits the data. For a full model with many clusters, each cluster’s contribution can differ, so agreement for one illustrative cluster should not be treated as a universal guarantee.

Common Pitfalls

Equating a rare event rate with a guaranteed poor Laplace approximation. Rarity alone does not determine approximation error; assess the shape over the region contributing to the integral.
Using the prior center as the conditional mode. The mode solves $h'(b)=0$, which includes the cluster’s observed responses as well as the random-effect penalty.
Assuming adaptive quadrature is exact because it is centered adaptively. Centering and scaling improve node placement, but convergence still depends on node count and integrand shape.
Comparing likelihoods with different constants. A binomial count likelihood may contain a combinatorial factor absent from the likelihood of a specified response sequence. It leaves the mode unchanged but changes the integral’s value.

Key takeaway: Rare binary responses can shift a random-intercept mode and produce asymmetric log-integrands. Laplace approximation summarizes only the local height and curvature; adaptive quadrature evaluates the nonlinear integrand and should be checked for numerical stability when accuracy matters.

Check Your Understanding

Answer each question, then compare with the solutions below.

For $m=8$, $y=0$, and $\sigma_b=1$, write the score equation for the conditional mode.
Why does $h''(b)<0$ guarantee that a stationary point is the unique mode in this model?
For $m=20$, $y=1$, $\sigma_b=1$, and $\beta_0=\log(0.05/0.95)$, calculate the mode and curvature.
What feature of adaptive quadrature can make it more accurate than Laplace approximation for an asymmetric integrand, and what check is still needed?

Solutions. (1) With $y=0$, the score is $h'(b)=-8p(b)-b$, so the mode solves $-8p(b)-b=0$. (2) Strict concavity means the derivative is strictly decreasing; it can cross zero at most once, and that stationary point is the global maximum. (3) At $b=0$, $p=0.05$, so $h'(0)=1-20(0.05)=0$; hence $\hat b=0$. The curvature is $H=20(0.05)(0.95)+1=1.95$. (4) Adaptive quadrature evaluates the original integrand at nodes centered and scaled using the mode and curvature, rather than replacing it everywhere by one quadratic. Increase the node count and verify that the integral stabilizes; adaptivity alone does not guarantee accuracy.