Why Refine the Quadrature?
A quadrature-based NLME fit replaces an integral over random effects with a finite weighted sum. The resulting likelihood is therefore an approximation, and the fitted parameters can depend on the numerical integration settings. A fit that converges numerically at one node order is not, by that fact alone, evidence that the integral is sufficiently accurate.
Node refinement provides a direct diagnostic: repeat the calculation with more quadrature nodes (or a finer quadrature level), then compare the estimated parameters and maximized log likelihood. The aim is not to prove that a particular fit is exact. It is to check whether further refinement materially changes the quantities used for inference.
For example, if subject \(i\) has observations \(y_i\), random effects \(b_i\), and conditional density \(p(y_i\mid b_i,\theta)\), the marginal likelihood contribution integrates over the random-effects distribution. Under the usual independent-subject formulation, \(L(\theta)\) is the product of these contributions. Quadrature approximates the integrals; refinement changes that approximation, not the statistical model.
What a Refinement Comparison Can Establish
A useful comparison asks two related questions. First, do the population parameter estimates change enough to affect interpretation or decisions? Second, does the maximized log likelihood stabilize? Both matter: a similar likelihood does not guarantee similar parameter estimates when the likelihood surface is flat, and similar parameters do not guarantee that the likelihood itself is well approximated.
Compare parameters on meaningful scales. A change of 0.05 may be negligible for a parameter near 100 but important for one near 0.1. Relative change can help for nonzero parameters, while an absolute change is more appropriate near zero or for parameters whose scale has a direct interpretation. For a positive parameter \(\theta_j\), one possible relative comparison is \(\lvert\hat{\theta}_{Q_2,j}-\hat{\theta}_{Q_1,j}\rvert/\lvert\hat{\theta}_{Q_2,j}\rvert\). Report the units and retain enough digits to see whether rounding is hiding a change.
For the objective, compare maximized log likelihoods, not likelihoods rounded to a few decimal places or values evaluated at different parameter vectors without explanation. A small change in log likelihood can still accompany an appreciable parameter shift if the surface is shallow. Conversely, a stable estimate may conceal a small but consequential likelihood difference. There is no universal threshold that establishes adequate accuracy for every model; tolerances should reflect the use of the fit and the precision required.
Why Refinement Is a Diagnostic, Not a Proof
For Gaussian random-effects integration, Gauss-Hermite quadrature is exact for polynomial integrands up to a degree determined by its order. In a nonlinear mixed-effects model, however, the integrand is generally not a polynomial. Higher order can improve the approximation for a sufficiently smooth, well-behaved integrand, but the sequence of estimates need not change monotonically. Orders may use different, non-nested nodes, and small differences between two orders do not themselves bound the remaining error.
The effect of integration error on a fitted parameter also depends on the shape of the objective. Suppose the exact log likelihood \(\ell(\theta)\) has a local maximum at \(\theta^\ast\), and write the quadrature approximation as \(\ell_Q(\theta)=\ell(\theta)+e_Q(\theta)\). At an interior maximum, the gradient of the exact log likelihood is zero. A first-order expansion of the approximate score around \(\theta^\ast\) gives
where \(H\) is the Hessian matrix of \(\ell\) at \(\theta^\ast\). This approximation explains an important practical point: even a small integration-induced score error can produce a noticeable parameter shift when the likelihood has weak curvature in some direction. The argument is local; it assumes a smooth objective, a nearby interior optimum, and an invertible Hessian. It is not an error bound for an arbitrary fit.
Worked Example: Refining a One-Dimensional Quadrature Estimate
Consider the hypothetical integral \(I=E[\exp(-Z^2/2)]\), where \(Z\sim N(0,1)\). This simple example isolates quadrature error; it is not an NLME study. Completing the normal-density calculation gives the exact value
For one-point Gauss-Hermite quadrature under the standard-normal convention, the node is 0 with weight 1, so \(I_1=1\), with absolute error \(\lvert1-0.707107\rvert=0.292893\).
For three points, the nodes are \(0,\pm\sqrt{3}\), with weights \(2/3,1/6,1/6\). Substitution gives
Its absolute error is approximately \(\lvert0.741043-0.707107\rvert=0.033937\). For five points, use the standard-normal nodes \(0,\pm1.3556,\pm2.8570\) and corresponding weights \(0.53333,0.22208,0.01126\) for the center, each inner node, and each outer node. Then
Using the displayed weights and nodes, \(I_5\approx0.7109\), and its absolute error is \(\lvert0.7109245-0.7071068\rvert\approx0.003818\). The successive estimates move closer in this example, but that pattern is not a general guarantee for nonlinear NLME integrals. The comparison is informative because the exact value is known here; in a fitted model it is usually not.
A Practical Refinement Workflow
A clean comparison requires more than changing an order setting. Keep the data, parameterization, covariance structure, optimization tolerances, and other model choices fixed. For each order, verify that the optimizer converged and record the maximized log likelihood and all parameters relevant to the intended use. If the optimizer reaches a different local solution, the comparison cannot be attributed to quadrature order alone.
Fit at an initial order and at one or more larger orders that are computationally feasible. Record the actual order or sparse-grid level, since a level and a node count are not interchangeable across schemes.
Confirm convergence, inspect gradients or other available optimization diagnostics, and rule out obvious boundary or failed evaluations before interpreting differences.
Calculate absolute or relative parameter changes on appropriate scales, and compare maximized log likelihoods. Preserve adequate precision and state units.
If changes remain material, increase the order or level again, if feasible. If computation prevents further refinement, report that limitation rather than treating the last result as confirmed.
Worked Example: Comparing Hypothetical NLME Fits
Suppose a hypothetical model is fitted to the same 40 subjects using Gauss-Hermite orders \(Q=5,9,15\). All three optimizations satisfy the stated convergence criteria and use the same model, starting strategy, and tolerances. The estimates are:
| Order \(Q\) | Clearance estimate (L/h) | Volume estimate (L) | Maximized log likelihood |
|---|---|---|---|
| 5 | 4.82 | 31.6 | -128.40 |
| 9 | 4.91 | 31.9 | -127.92 |
| 15 | 4.93 | 31.9 | -127.91 |
From \(Q=5\) to \(Q=9\), clearance changes by \(4.91-4.82=0.09\) L/h, or \(0.09/4.91\approx1.83\%\) relative to the \(Q=9\) estimate. Volume changes by \(0.3\) L, and the maximized log likelihood changes by \(-127.92-(-128.40)=0.48\). These changes suggest that order 5 may be too coarse for a conclusion requiring close agreement.
From \(Q=9\) to \(Q=15\), clearance changes by \(0.02\) L/h, or \(0.02/4.93\approx0.406\%\); volume is unchanged to the displayed precision, and the log likelihood changes by \(0.01\). These results are substantially more stable. Whether that stability is sufficient depends on the application and the required precision; the table alone does not prove the exact integration error is below any specified bound.
Diagnostics, Alternatives, and Limitations
Node refinement should be interpreted alongside basic numerical checks. Parameter changes concentrated in a weakly identified direction may reflect a flat likelihood as well as integration sensitivity. Changes in random-effect variances or correlations can be especially difficult to summarize by a single relative percentage, particularly near a boundary. Inspect the parameter values on their natural scales, the likelihood, and whether the optimizer consistently finds the same solution.
The main alternatives address different computational trade-offs. Adaptive quadrature places nodes using information about the individual-effect integrand; it can be efficient when that integrand is well represented near the chosen center and scale, but its accuracy still needs assessment. Sparse-grid quadrature can reduce the cost of a full tensor product in higher dimensions, but refinement comparisons across sparse-grid levels remain diagnostics rather than guarantees. Importance sampling estimates integrals using weighted draws and therefore has Monte Carlo variability; its error assessment is different from deterministic node refinement. A Laplace approximation is cheaper in many settings but uses a local approximation to the integrand and can be less accurate when that approximation is inadequate. No method is uniformly best for every dimension and likelihood shape.
A common mistake is to reason, “The parameter estimate barely changed between two orders, so the numerical integration must be exact.” The correction is that two close estimates are evidence of local stability across those settings, not proof of accuracy: both orders can share a bias, and a flat likelihood can make estimates appear stable despite approximation error. Another mistake is to treat an increase in the maximized log likelihood at higher order as a universal requirement. Approximate objectives need not change monotonically with order, and optimizer differences can complicate the comparison. Use several refinements when feasible, verify optimization, and describe residual uncertainty honestly.
The diagnostic also has practical limits. High orders can be prohibitively expensive as the number of random effects grows, and numerical underflow, poor scaling, or unstable optimization can obscure the comparison. A stable sequence does not diagnose model misspecification, and it does not replace examination of the statistical model. When results remain sensitive and higher-order computation is infeasible, the sensitivity itself is an important result to report.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- A parameter estimate changes from 2.00 to 2.02 as quadrature order increases. What is the absolute change, and what additional information is needed to judge whether it is material?
- Why is it unsafe to conclude that an integral is accurate merely because two successive quadrature orders give similar estimates?
- A fit's maximized log likelihood changes by 0.30 between orders, while a reported parameter is nearly unchanged. Should the parameter stability alone end the diagnostic? Explain.
- State two conditions that should be held fixed or checked when comparing fits at different node orders.
- Why might a small score error produce a comparatively large parameter change?
Solutions. (1) The absolute change is \(2.02-2.00=0.02\). Its importance depends on the parameter's units, scale, intended use, and changes at further refinements. (2) The approximations may share error; two orders do not provide a general error bound, and the sequence need not be monotone. (3) No. The likelihood change should also be assessed, and its practical importance depends on context and required precision. (4) Keep the data, model, parameterization, and optimization tolerances fixed; check that each optimizer converged and found a comparable solution. (5) The local relation \(\theta_Q-\theta^\ast\approx-H^{-1}\nabla e_Q(\theta^\ast)\) shows that weak curvature, corresponding to a poorly conditioned or nearly singular Hessian, can amplify a small score perturbation.