Quadrature Accuracy And Apparent Convergence

Why Convergence Can Be Misleading

A joint-model fitting algorithm can stop successfully while the estimates still depend materially on how its random-effects integrals were approximated. These are two different numerical questions. Optimization convergence asks whether the algorithm has located a stationary point of the objective it was given. Quadrature accuracy asks whether that objective is a sufficiently accurate approximation to the marginal likelihood. A small gradient answers the first question, not the second.

This distinction matters because likelihood contributions in shared-random-effects models integrate over subject-specific latent effects. If an integration rule is too coarse, the optimizer may precisely maximize an inaccurate likelihood. Conversely, estimates that change between runs at the same quadrature order may indicate optimizer difficulties rather than integration error. The practical task is to vary the integration accuracy in a controlled way and assess whether the substantive results persist.

Marginal Likelihood and Quadrature

Let \(y_{ij}\) denote the longitudinal outcome for subject \(i\) at visit \(j\), \(T_i\) the observed event or censoring time, and \(\delta_i\) the event indicator. Let \(b_i\) be the subject-specific random-effects vector, with \(b_i\sim N_q(0,D)\). Write \(\theta\) for all model parameters. Conditional on \(b_i\), let \(L_i(b_i;\theta)\) be the subject's joint contribution from the longitudinal observations and survival outcome. For example, the survival part has the form \(\lambda_i(T_i\mid b_i;\theta)^{\delta_i}\exp\{-\Lambda_i(T_i\mid b_i;\theta)\}\), where \(\lambda_i\) is the conditional hazard and \(\Lambda_i\) its cumulative hazard. This factorization presumes the model's conditional likelihood assumptions and an appropriate treatment of censoring.

Definition: The marginal likelihood contribution for subject \(i\) is \(L_i(\theta)=\int L_i(b_i;\theta)\phi_q(b_i;0,D)\,db_i\), where \(\phi_q(\cdot;0,D)\) is the \(q\)-dimensional normal density. The marginal log likelihood is \(\ell(\theta)=\sum_i\log L_i(\theta)\). It integrates out the random effects rather than treating their unobserved values as known.

To see the role of Gaussian quadrature, write \(D=CC^\mathsf{T}\) and set \(b_i=Cz_i\), where \(z_i\sim N_q(0,I)\). Then the subject contribution is an expectation of \(g_i(z;\theta)=L_i(Cz;\theta)\) under a standard multivariate normal distribution. In one dimension, an \(m\)-point standard-normal Gauss-Hermite rule approximates an expectation by a weighted sum:

\int g(z)\phi(z)\,dz \ \approx\ Q_m(g)=\sum_{k=1}^{m}w_k g(x_k), \qquad \sum_{k=1}^{m}w_k=1.

In \(q\) dimensions, a tensor-product rule sums over combinations of nodes in all dimensions; with \(m\) nodes per dimension it uses \(m^q\) combinations for each subject. Increasing \(m\) makes the rule more capable of representing a curved integrand, but the computational cost can rise rapidly. The rule is most effective when the transformed integrand is adequately smooth and the nodes and weights match the specified normal expectation. Its finite-order approximation is not the exact marginal likelihood.

Worked Example: A Five-Point Normal-Expectation Approximation

Consider the standard-normal expectation \(E[\exp(-Z^2/2)]\), where \(Z\sim N(0,1)\). The exact value is \(1/\sqrt{2}\approx0.707107\). Use a five-point standard-normal Gauss-Hermite rule to approximate it.

The nodes are \(0,\ \pm\sqrt{5-\sqrt{10}},\ \pm\sqrt{5+\sqrt{10}}\). The center weight is \(8/15\); the weight for each inner node is \((7+2\sqrt{10})/60\); and the weight for each outer node is \((7-2\sqrt{10})/60\). Since the function is even, the positive and negative node contributions are equal. Thus:

Q_5=\frac{8}{15} +2\left(\frac{7+2\sqrt{10}}{60}\right)e^{-(5-\sqrt{10})/2} +2\left(\frac{7-2\sqrt{10}}{60}\right)e^{-(5+\sqrt{10})/2} \approx 0.710918.

The absolute error is approximately \(|0.710918-0.707107|=0.003811\), about \(0.54\%\) of the exact value. The approximation is close, but it is not exact. This calculation illustrates why a plausible-looking quadrature result should be checked for stability at higher orders rather than assumed accurate from its appearance.

What Increasing Quadrature Order Can Establish

For a sufficiently regular integrand, Gaussian quadrature approximations improve as the number of points increases; under appropriate conditions, they converge to the integral. This is an asymptotic statement, not a guarantee that every finite increase produces a smaller error. In a fitted model, the integrands vary across subjects and depend on \(\theta\). A point count that works well for one subject or parameter value may be inadequate for another, particularly when the integrand is sharply concentrated or has substantial tail contribution.

Let \(\ell_m(\theta)\) denote the log likelihood computed using \(m\) quadrature points per dimension, and let \(\hat\theta_m\) be a maximizer found for that approximate objective. The optimizer may satisfy \(\nabla\ell_m(\hat\theta_m)\approx0\), while \(\hat\theta_m\) is still noticeably different from the estimate obtained under a more accurate integral. A short justification shows why. Write \(\ell_m(\theta)=\ell(\theta)+\Delta_m(\theta)\), where \(\Delta_m\) is the integration-induced error. Near a well-behaved maximum \(\hat\theta\) of \(\ell\), a first-order expansion of the score gives:

0\approx\nabla\ell_m(\hat\theta_m) \approx \nabla\ell(\hat\theta) +\nabla^2\ell(\hat\theta)(\hat\theta_m-\hat\theta) +\nabla\Delta_m(\hat\theta).

Because \(\nabla\ell(\hat\theta)=0\), and if the Hessian is nonsingular, this implies approximately \(\hat\theta_m-\hat\theta\approx-[\nabla^2\ell(\hat\theta)]^{-1}\nabla\Delta_m(\hat\theta)\). Thus even a modest integration error can shift estimates, especially in directions with little likelihood curvature. A converged optimizer only establishes near-stationarity for \(\ell_m\); it does not make \(\Delta_m\) vanish.

Worked Example: Checking Whether Estimates Stabilize

Suppose a hypothetical joint-model analysis estimates an association parameter \(\alpha\), measuring the change in log hazard per one-unit increase in the modeled marker value. The same model is fitted at successively larger quadrature orders. Assume the optimizer has converged at each order. The results are:

Points per dimension, \(m\)	\(\hat\alpha_m\)	Approximate log likelihood	Gradient norm
5	0.742	-812.46	0.0008
9	0.716	-811.92	0.0006
15	0.711	-811.89	0.0007
25	0.710	-811.89	0.0009

The gradient norms are small, which supports convergence of each optimization, but do not by themselves validate the integration. Compare estimates at the larger orders: \(|0.710-0.711|=0.001\), small relative to the estimate. The approximate log likelihoods at 15 and 25 points agree to the displayed precision. These results support numerical stability over the tested refinement, while the shift from 5 to 15 points shows that the lower-order estimate was not as stable.

If the longitudinal marker is measured in its original units, \(\alpha\) has units of inverse marker units; the hazard ratio for a one-unit increase is \(\exp(\alpha)\). For \(\hat\alpha_{25}=0.710\), this is \(\exp(0.710)\approx2.03\). The numerical check concerns stability of the estimate and its implied hazard ratio, not whether the association is causal or whether the model is otherwise adequate.

Separating Integration Error from Optimizer Failure

A useful diagnostic keeps the two numerical components distinct. First, at a fixed quadrature order, check whether repeated starts or appropriate optimization runs reach the same solution with small gradients. Second, assess integration by refining the quadrature while monitoring estimates and the objective. To sharpen the comparison, evaluate the higher-order objective at the lower-order estimate as well as at the higher-order optimum. At any one fixed order, these objective values can be compared directly; values from different orders are evaluations of different approximations and should not be treated as if they were identical likelihood calculations.

Worked Example: A Coarse-Rule Optimum That Does Not Persist

Consider a hypothetical fit with one association parameter \(\alpha\). At nine points, two starts return \(\alpha=0.48\) with approximate log likelihood \(-500.21\) and gradient norm \(0.12\), and \(\alpha=0.90\) with approximate log likelihood \(-500.06\) and gradient norm \(0.0006\). A further start reaches the second solution. The first result is not persuasive evidence of a competing optimum: its gradient norm is much larger, while the repeated solution has a better objective at the same quadrature order.

Now evaluate the nine-point estimate \(\alpha=0.90\) using 25 points, obtaining \(\ell_{25}(0.90)=-499.72\). Then optimize the 25-point objective, obtaining \(\hat\alpha_{25}=0.76\), \(\ell_{25}(0.76)=-499.61\), and gradient norm \(0.0005\). The comparison at 25 points is direct: \(-499.61>-499.72\), so the refined fit improves that objective. The estimate has shifted by \(|0.76-0.90|=0.14\).

The nine-point optimizer had found a stationary point of its approximate objective, but its estimate did not persist under refined integration. The gradient at the 25-point solution supports optimization convergence there; stability still requires checking another refinement or otherwise establishing that further integration changes are negligible for the intended use. The log likelihoods \(-500.06\) and \(-499.61\) come from different approximations, so their difference alone is not a clean measure of integration error.

Diagnostics, Alternatives, and Limitations

Evidence for quadrature stability is strongest when estimates of primary interest, derived quantities such as hazard ratios, and relevant uncertainty summaries change negligibly across further refinements. A useful check is to track these alongside the approximate log likelihood and optimization diagnostics. If estimates move but gradients remain small at each order, integration accuracy is a leading concern. If repeated fits at the same order disagree, or gradients remain large, optimization is not yet adequately resolved. Both problems can occur together.

Parameter shifts across orders: a material change suggests that the reported result may be sensitive to integration approximation.
Small gradient but unstable estimates: the optimizer may have converged accurately to different approximate objectives.
Unstable gradients or different same-order endpoints: investigate optimization behavior before attributing the difference to quadrature.
Uncertainty changes: standard errors based on curvature inherit numerical error in the likelihood and its derivatives, so stable point estimates alone are not sufficient evidence that uncertainty estimates are stable.

Quadrature is not the only way to approximate random-effects integrals. The Laplace approximation uses a local quadratic approximation near an integrand's mode and can be computationally attractive, but it may be inaccurate when that approximation represents the integrand poorly. Monte Carlo integration uses random draws and can be useful when tensor-product quadrature becomes costly in higher dimensions, at the price of simulation variability. These alternatives have their own numerical diagnostics; neither makes an optimization stopping criterion a substitute for checking the integral approximation.

The main limitation of increasing tensor-product order is cost: with \(q\) random effects, \(m\) points per dimension require \(m^q\) combinations per subject. High dimensionality can therefore make repeated refinement impractical. Quadrature stability also does not establish model correctness or resolve weak identification. When the likelihood is flat in a parameter direction, small numerical perturbations can cause substantial estimate changes; interpretation should take that fragility into account.

Common Mistake

Pitfall: “The algorithm stopped, the gradient is nearly zero, and the software reported convergence; therefore the random-effects integrals must be accurate.” The reasoning confuses optimization of the computed objective with accuracy of that objective. A small gradient means the fitted value is near a stationary point of \(\ell_m\), not necessarily of the exact marginal log likelihood \(\ell\). Refine the quadrature and check whether the estimates and conclusions persist; separately verify optimization at each order.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

A fit has a gradient norm of \(0.0004\) at 7 quadrature points. What does this establish, and what does it not establish?
An association estimate is 0.52 at 7 points and 0.51 at 15 points. Calculate the absolute change and state what this comparison can suggest.
At 20 points, the log likelihood at a candidate estimate is \(-320.4\), and after optimization it is \(-319.9\). Which value is larger, and what can be concluded at this fixed order?
A tensor-product rule uses 9 points in each of four random-effects dimensions. How many node combinations does it use per subject?

Solutions. (1) The small gradient supports near-stationarity of the 7-point approximate objective. It does not establish that the integral approximation is accurate or that estimates are stable at higher orders. (2) The absolute change is \(|0.52-0.51|=0.01\). This supports stability over those tested orders, but does not prove further refinements will have no effect. (3) \(-319.9>-320.4\). At 20 points, the optimized candidate has the larger objective by \(0.5\); this comparison is valid because both values use the same rule, though optimization diagnostics are still needed. (4) The number is \(9^4=6{,}561\) node combinations per subject.