Why Numerical Tolerances Matter
A joint model’s likelihood integrates over subject-specific random effects. In general, this integral cannot be evaluated exactly, so the fitted log likelihood and parameter estimates depend partly on numerical approximations. A computation can report convergence while still being sensitive to the quadrature order or integration accuracy used to calculate its objective. The practical question is therefore not simply whether an algorithm stopped, but whether the reported fit is stable enough for its intended use.
The earlier material on Gauss-Hermite quadrature, adaptive quadrature, and Laplace approximation establishes the main integration approaches. Here the emphasis is on controlling and checking their numerical accuracy. These checks do not establish that the joint model is scientifically appropriate; they address whether the chosen numerical procedure has calculated the fit reliably.
Three Different Numerical Choices
Let \(Y_i\) denote subject \(i\)’s longitudinal data, and let \(T_i\) and \(\delta_i\) denote the observed event or censoring time and event indicator. Write \(b_i\) for that subject’s random-effects vector and \(\theta\) for all model parameters. A general subject-level likelihood contribution is
The conditional contribution and random-effects density must define a finite marginal likelihood for the parameter values considered. The integration method must also be applied consistently across subjects and parameter evaluations. Under these conditions, numerical methods approximate \(L_i(\theta)\), and hence \(\ell(\theta)\).
For fixed-order Gauss-Hermite quadrature, the order determines the nodes and weights; there may be no separate user-specified integration tolerance. An adaptive integrator may instead refine its numerical work until an absolute, relative, or combined error criterion is met. An optimizer tolerance concerns changes in the objective, parameters, or gradient. The exact meaning of each setting depends on the software, so its documentation matters.
A small optimizer stopping criterion cannot compensate for an inaccurate likelihood approximation. If the objective changes irregularly as parameters change because integration is too coarse, the optimizer may stop at a point that appears converged for the approximate objective but is not stable under more accurate integration.
How Integration Error Affects the Log Likelihood
Because likelihood contributions are integrated and then logged, absolute integration error alone is not enough to judge its effect: the same absolute error can be negligible for a large contribution and important for a small one. A useful way to see the relationship is to suppose that the numerical approximation \(Q_i\) satisfies a relative error bound
Then \(Q_i/L_i\) lies between \(1-\varepsilon_i\) and \(1+\varepsilon_i\). Taking logarithms gives
For \(0\leq\varepsilon_i<1\), the magnitude of either endpoint is bounded above by \(-\log(1-\varepsilon_i)\). Summing over subjects therefore yields
For small \(\varepsilon_i\), \(-\log(1-\varepsilon_i)\) is approximately \(\varepsilon_i\). This bound clarifies why many small subject-level errors can accumulate in the sample log likelihood. In practice, software tolerances are often stopping criteria rather than certified error bounds, so this calculation explains the scale of the issue but does not certify a particular fit.
Worked Example: Interpreting Relative Integration Error
Suppose a hypothetical subject has marginal contribution \(L_i=0.0200\), and a numerical approximation is \(Q_i=0.02001\). The absolute error is \(0.00001\), so the relative error is
\[ \frac{|Q_i-L_i|}{L_i}=\frac{0.00001}{0.0200}=0.000500. \]
The change in the subject’s log contribution is
If, hypothetically, 800 subjects each had a relative error of this size in the same direction, the sample log likelihood could differ by about \(800(0.000500)=0.400\). This is an illustration, not a general prediction: actual subject errors need not be equal or have the same sign. It does show why a small error for one contribution does not by itself guarantee a negligible error in the total log likelihood.
Refinement Checks for a Fitted Model
A useful check changes the numerical accuracy while holding the model specification and data fixed. With fixed-order quadrature, increase the order and refit. With an adaptive integration method, tighten the integration criterion and refit. In either case, compare both the maximized log likelihood and the parameter estimates. Use the same parameter constraints and take care that the fits are comparable; different local optima or failed runs can otherwise be mistaken for integration effects.
The following hypothetical fit illustrates a quadrature-order check. The model, data, and optimization settings are held fixed; only the quadrature order changes.
| Nodes per dimension | Maximized log likelihood | Association estimate |
|---|---|---|
| 7 | -842.613 | 0.684 |
| 11 | -842.487 | 0.701 |
| 15 | -842.481 | 0.702 |
| 21 | -842.481 | 0.702 |
Worked Example: Assessing Quadrature-Order Stability
For the hypothetical results above, compare the final two orders. The log-likelihood change from 15 to 21 nodes is
The association estimate changes by \(0.702-0.702=0.000\) at the displayed precision. From 7 to 11 nodes, however, the log likelihood changes by \(-842.487-(-842.613)=0.126\), and the estimate changes by \(0.701-0.684=0.017\). The lower-order results are therefore visibly less stable than the last two.
The table supports the conclusion that these particular reported quantities are stable to increasing the order from 15 to 21, at the displayed precision. It does not prove exactness, establish stability for every parameter, or show that the model is correct. A more demanding analysis would also examine important variance and survival parameters, and would compare changes with the precision needed for the scientific interpretation.
Choosing and Reporting Tolerances
There is no universal quadrature order or tolerance that is adequate for every joint model. Integration difficulty depends on factors such as random-effects dimension, posterior concentration, and the shape of the integrand. A practical choice is justified by refinement: results should not change materially when numerical accuracy is increased beyond the setting used for the reported fit.
A defensible workflow separates integration checks from optimization checks:
Record the quadrature rule and order, integration settings, optimizer stopping criteria, starting values, and convergence information.
Raise the quadrature order or tighten the integration tolerance while keeping the statistical model unchanged. Refit rather than merely reevaluating at an old estimate when assessing the effect on parameter estimates.
Check the maximized log likelihood and the estimates that matter to the analysis, including uncertainty-related parameters when relevant. Compare changes with the reporting precision and inferential purpose.
At the chosen integration accuracy, verify that the optimizer’s convergence criteria are met. If tightening them changes the fit, investigate whether the earlier optimization had stopped too soon.
For an adaptive routine, suppose an absolute error criterion is applied to an integral near \(0.0200\). A target absolute error of \(10^{-5}\) corresponds to a relative scale of about \(10^{-5}/0.0200=0.000500\), or 0.05%. A relative criterion would scale its target with the integral’s magnitude. This conversion helps interpret settings, but the routine’s reported criterion is not necessarily a guaranteed bound on the final sample log likelihood.
Worked Example: Distinguishing Integration and Optimizer Tolerances
In a hypothetical fit, the baseline setting uses 9 quadrature nodes and an optimizer tolerance of \(10^{-6}\). The reported association estimate is 0.72 and the log likelihood is \(-510.24\). A second run keeps the optimizer tolerance unchanged but uses 15 nodes; it gives an estimate of 0.76 and a log likelihood of \(-510.10\). A third run keeps 15 nodes but tightens the optimizer tolerance to \(10^{-8}\); it gives an estimate of 0.760 and a log likelihood of \(-510.100\).
Comparing the first and second runs isolates a substantial quadrature sensitivity: the estimate shifts by \(0.76-0.72=0.04\), and the log likelihood increases by \(-510.10-(-510.24)=0.14\). Comparing the second and third runs changes only the optimizer tolerance. The estimate changes by \(0.760-0.76=0.000\) at the displayed precision, and the log likelihood changes by \(-510.100-(-510.10)=0.000\) at that precision.
These results suggest that tightening the optimizer criterion did not materially change the reported fit, while increasing quadrature order did. The baseline optimizer’s small tolerance did not make its low-order integration accurate. Further quadrature refinement would still be needed to assess whether the 15-node fit is stable.
Diagnostics, Alternatives, and Limitations
Warning signs include estimates that move appreciably as quadrature order increases, a log likelihood that continues to change at the highest tested order, repeated optimizer failures under tighter integration settings, or materially different results from different starting values. Also inspect whether estimated variances approach boundaries or whether the integrand is sharply concentrated in regions poorly represented by the chosen nodes. Such signs do not identify a single cause, but they mean a numerical stability claim needs further investigation.
A fixed quadrature rule is straightforward to reproduce and refine by increasing its order. Adaptive quadrature can concentrate numerical effort around an important region of the integrand, but its accuracy depends on how well its adaptation captures that region. Sparse-grid methods can reduce the cost of high-dimensional integration relative to a tensor product, though their accuracy should also be assessed by refinement. The Laplace approximation is a different approximation to the integral, not merely a tighter quadrature tolerance. Comparing it with quadrature can be informative, but agreement or disagreement alone does not certify either method.
A numerical tolerance cannot repair a misspecified random-effects distribution, an incorrect association structure, or inadequate data information. Nor does stability of the log likelihood alone guarantee stable inference: estimates of a particular parameter may remain sensitive even when changes in the total log likelihood look small. Finally, a reported convergence code is evidence about the algorithm’s stopping rule, not proof that numerical integration error is negligible.
Worked Example: Using a Log-Likelihood Change as a Warning
Suppose a hypothetical fit has maximized log likelihood \(-1{,}204.82\) at 11 nodes and \(-1{,}204.79\) at 17 nodes. The change is
At the same time, the estimated random-effects standard deviation changes from 0.31 to 0.37. The log-likelihood change is numerically small relative to its total magnitude, but the standard deviation changes by \(0.37-0.31=0.06\), about \(0.06/0.31\approx 19.4\%\) of the 11-node estimate. That parameter may matter directly to interpretation or uncertainty quantification. The appropriate conclusion is not that the fit is stable because the total log likelihood barely changed; rather, the variance estimate needs additional refinement checks and investigation.
Common Mistake
Wrong reasoning: “The optimizer converged to six decimal places, so the likelihood and estimates must be accurate.” This confuses the optimizer’s stopping tolerance with the accuracy of the integrated likelihood it was optimizing. An optimizer can precisely maximize an inaccurate approximation.
Correction: Check integration accuracy by changing quadrature order or tightening the integration criterion, and check optimizer convergence separately. Report the numerical settings and assess stability in the log likelihood and the parameters important to the analysis. Use a precision standard appropriate to the inferential task rather than relying on a convergence flag alone.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- In a quadrature fit, what does increasing the quadrature order change, and what does tightening the optimizer tolerance change?
- A subject-level contribution is \(0.0500\), and its approximation is \(0.0501\). Calculate the relative error and the change in its log contribution to a useful approximation.
- Why is a stable total log likelihood not sufficient evidence that every parameter estimate is numerically stable?
- A fit changes little when the optimizer tolerance is tightened but changes noticeably when quadrature order is increased. Which numerical issue is implicated?
Solutions. (1) Increasing quadrature order adds or changes integration nodes and weights, improving the integration approximation when the rule converges appropriately. Tightening optimizer tolerance requires stricter stopping conditions for the parameter search; it does not itself improve the integral approximation. (2) The relative error is \((0.0501-0.0500)/0.0500=0.00200\), or 0.200%. The log change is \(\log(0.0501/0.0500)=\log(1.002)\approx 0.001998\). (3) The total log likelihood aggregates contributions and can be insensitive to a parameter change that remains important for interpretation; estimates must be checked directly. (4) The results point to sensitivity to integration accuracy, particularly quadrature order, rather than a material effect of the optimizer stopping criterion.