Quadrature Order and Numerical Convergence

Why Quadrature Order Needs to Be Checked

In a generalized linear mixed model, the marginal likelihood involves integrating over random effects. When the integral has no convenient exact form, Gauss–Hermite quadrature replaces it with a finite weighted sum. The number of nodes used in that sum is the quadrature order. A low order is computationally cheaper, but its integral approximation may be inaccurate; a higher order usually costs more and is not automatically a guarantee of accuracy. The practical question is whether increasing the order changes the likelihood or the fitted parameters enough to matter for the analysis.

There are two related but distinct checks. First, hold the model parameters fixed and see whether the approximated marginal log likelihood stabilizes as the order increases. Second, refit the model at each order and see whether the parameter estimates stabilize. The first isolates integration error at a particular parameter value. The second checks whether the integration approximation is affecting the resulting fit. A stable likelihood at one fixed parameter value does not by itself establish that refitted estimates are stable.

For example, if cluster \(i\) has random effect \(b_i\), its exact marginal contribution can be written as $$ L_i(\theta)=\int \exp\{h_i(b_i;\theta)\}\,db_i, $$ where \(h_i\) includes the conditional data log likelihood and the random-effects density on the integration scale. Adaptive Gauss–Hermite quadrature evaluates a weighted sum at nodes selected using a local approximation around the cluster’s conditional mode. In the notation above, the weights include the transformation and scaling needed for that sum. The purpose here is not to derive the nodes again, but to assess how the resulting numerical approximation behaves as \(Q\) changes.

What Numerical Convergence Means

A sequence of quadrature calculations is useful only if the calculations being compared differ in the intended way. Use the same observations, model specification, link, random-effects structure, and parameterization at every order. For a fixed-parameter check, evaluate \(\widehat{\ell}_Q(\theta)\) at exactly the same \(\theta\) for each \(Q\). For a refitting check, optimize the order-specific objective at each \(Q\), using consistent convergence settings. Otherwise, a change in the objective could reflect changed data or an optimizer stopping at a different point rather than quadrature order.

A simple diagnostic is the successive log-likelihood difference:

If \(|\Delta_Q|\) is small relative to a tolerance chosen for the analysis, the log likelihood has stabilized across that comparison. There is no universal tolerance appropriate for every model and purpose. The number of observations, the scale of the log likelihood, and the precision needed for a later comparison all matter. Report the orders and differences rather than declaring convergence solely because software used a particular default.

The same idea applies to estimates. For parameter component \(j\), compare \(\widehat{\theta}_{Q_1,j}\) with \(\widehat{\theta}_{Q_2,j}\), and consider both the absolute difference and its practical size in context. A change of 0.01 in a log-odds coefficient may be consequential in one application and negligible in another. For a standard deviation estimate near zero, a small absolute change can be a large relative change. It is therefore helpful to report the parameter’s units or interpretation as well as the raw numerical difference.

Why Increasing the Order Should Help—and Why It May Not Look Smooth

Gauss–Hermite rules approximate integrals by weighted sums. Under suitable regularity conditions, increasing the number of nodes makes the approximation converge for integrands that are sufficiently smooth and well behaved under the Gaussian weighting. The adaptive transformation aims to make the integrand more nearly Gaussian near its important region, which can make a moderate order effective. But the marginal likelihood integrand need not be well approximated by its local shape everywhere. Skewness, tails, or multiple important regions can require more nodes.

The basic justification for checking successive orders is continuity of the log function and convergence of the integral approximation. Suppose, for a fixed \(\theta\), that \(\widehat{L}_{i,Q}(\theta)\) converges to a finite positive \(L_i(\theta)\) for every cluster. Then:

With a finite number \(m\) of clusters, summing these convergent terms gives \(\widehat{\ell}_Q(\theta)\to\ell(\theta)\). This explains why the log likelihood should stabilize when the individual integral approximations do. It does not imply that successive values must move monotonically toward the limit. Gauss–Hermite rules at different orders are generally not nested: a higher-order rule can produce a value on either side of a lower-order value.

For estimates, the argument requires more. If the approximate log likelihoods converge uniformly over the relevant parameter region and the exact log likelihood has a unique, well-separated maximizer, then their maximizers can converge to that maximizer. If the objective is flat, has multiple maxima, or is poorly optimized, parameter estimates can remain unstable even when the objective values are very similar. Thus, check likelihood and estimates separately.

Refitting at Successively Higher Orders

A fit at order \(Q\) maximizes \(\widehat{\ell}_Q(\theta)\), not the exact marginal log likelihood. Therefore, changing \(Q\) can change both the objective surface and its maximizer. To check this effect, fit the same model separately at several orders, then compare the estimates. Keep optimization tolerances sufficiently tight that optimizer error is not mistaken for quadrature error. If results differ, repeat fits from sensible starting values and inspect whether the optimizer reports convergence.

Diagnosing Apparent Convergence

A sequence can be irregular before it stabilizes. For instance, a fixed-parameter log likelihood might be \(-74.620\), \(-74.310\), \(-74.350\), \(-74.344\), and \(-74.343\) at orders 5, 9, 13, 21, and 31. The change from order 9 to 13 reverses direction, so the sequence is not monotone. The later values are much closer together, which is evidence of stabilization at those higher orders, not a reason to discard the intermediate values. When successive differences remain appreciable or oscillate, increase the order further if computationally feasible and investigate the integrand and fit.

Several checks help distinguish quadrature behavior from other numerical problems:

• Hold parameters fixed when diagnosing integration error in the objective. Then separately compare fully refitted models.
• Check optimizer convergence at each order. A small change in the returned log likelihood can occur even if an optimizer stopped prematurely.
• Compare more than two orders when results are irregular. One small difference alone may be accidental.
• Inspect all important parameters, including variance or standard-deviation parameters. A stable fixed-effect estimate does not imply that every component is stable.
• Use adequate precision in stored results. Rounding output too early can make genuinely different estimates look identical.

If a likelihood is nearly flat in a parameter direction, two fits may have almost identical log likelihoods but noticeably different parameter estimates. This is not automatically evidence of quadrature failure; it can reflect weak information in the data. Conversely, a parameter estimate that appears stable to three decimals could conceal relevant changes in the likelihood or in another parameter. Review the objective and the full parameter vector together.

Choosing Among Numerical Approaches

Quadrature-order checks complement, rather than replace, other integration approaches. The Laplace approximation uses a local Gaussian approximation at the conditional mode and is often less costly than using many quadrature nodes. It can be accurate when the integrand is close to Gaussian in the important region, but can be less reliable when that local approximation misses skewness or other structure. Adaptive Gauss–Hermite quadrature uses multiple weighted evaluations around a local approximation and allows its accuracy to be assessed by increasing \(Q\). A closed-form integral, when available, provides a useful benchmark; it is not available for many GLMM likelihoods.

These approaches answer different practical needs. Laplace approximation may be attractive when speed is essential, while adaptive quadrature can provide a more direct numerical accuracy check at a computational cost. If increasing order changes the estimates or log likelihood materially, a lower-order result should not be treated as numerically settled. If computation becomes prohibitive as the number of random-effect dimensions grows, the appropriate response is to consider a suitable integration strategy, not to label an unstable low-order result converged.

Most importantly, quadrature convergence is a numerical property. It does not diagnose a wrong link function, omitted covariates, dependence misspecification, or inadequate data. Those are questions of model assessment and study design. A highly stabilized approximation can still describe a poor model, while a sound model can be difficult to fit accurately with insufficient quadrature.

Common Mistake and Key Takeaway

Common mistake: “The log likelihood changed by less than 0.01, so the model has converged and the estimates are reliable.” This reasoning confuses one objective-value comparison with both optimizer convergence and parameter stability. A flat objective can yield substantially different parameter estimates, and one pair of nearby orders may not reveal later changes.

Correction: State the tolerance and compare several increasing orders. Check the fixed-parameter likelihood if the goal is to assess integration at a chosen parameter vector, then refit at each order and compare the full set of estimates. Verify optimizer convergence and interpret changes on meaningful parameter scales. Conclude only that the numerical results are stable to the tested orders and tolerances—not that the model is correct or the estimates have no statistical uncertainty.