Comparing Numerical Integration Settings

Why Compare Integration Settings?

A fitted shared-random-effects model can look well behaved while still depending on how its likelihood integral was evaluated. A small change in an integral approximation can shift the likelihood surface, move the reported estimates, or alter the curvature used to calculate standard errors. Therefore, convergence under one numerical setting is not, by itself, evidence that the numerical integration is sufficiently accurate.

The practical question is comparative: if the same model is refitted with more accurate or differently constructed integration settings, do the log likelihood, estimates, and uncertainty measures remain materially unchanged? This is a numerical-stability assessment, not a test of a scientific hypothesis. The earlier material on quadrature accuracy and apparent convergence motivates checking integration error; here the emphasis is on a controlled comparison across complete refits.

Likelihood and Integration Settings

For subject \(i\), let \(y_{ij}\) be the longitudinal response at time \(t_{ij}\), \(T_i\) the observed event or censoring time, and \(\delta_i\) the event indicator. Let \(b_i\) be the subject-specific random-effects vector. A representative shared-random-effects model is

with conditionally independent measurement errors \(\epsilon_{ij}\sim N(0,\sigma^2)\), and a conditional hazard such as

Here \(m_i(t;b_i)\) is the modeled longitudinal value linked to the hazard, and \(\alpha\) is the association parameter. The exact specification may differ, but the comparison principle is the same. Assume subjects are independent; conditional on \(b_i\) and covariates, the longitudinal and event processes have the stated model distributions; and censoring is noninformative under the model. The baseline hazard and all other model components must be held fixed during the numerical comparison.

The subject contribution integrates out the unobserved random effects. Writing \(f_y(y_i\mid b_i)\) for the longitudinal density, and \(h_i\) and \(S_i\) for the conditional hazard and survival function, respectively, gives

The parameter vector \(\theta\) includes the longitudinal, survival, association, and random-effects parameters. Except in special cases, this integral is evaluated numerically. Settings to compare can include the quadrature rule, the number of nodes, whether nodes are adapted to the subject-specific integrand, and any integration tolerance used by the software. A higher node count is not automatically better in every implementation: the useful evidence is whether independently tightened settings give stable results.

A Controlled Comparison

A fair comparison changes integration settings deliberately, not several aspects of the analysis at once. Use the same observations, model formula, parameter constraints, baseline-hazard specification, and likelihood constants. Use sufficiently reliable optimization settings for every fit, and check that each fit reached a comparable optimum. If starting values vary, distinguish variation due to competing optima from variation due to integration. The earlier material on multiple starting values and gradient checks is relevant to those checks, but integration comparisons do not replace them.

For each refit, record at least the maximized log likelihood, the parameters of substantive interest, and their standard errors. It is also useful to record integration rule and settings, optimizer status, and elapsed computation time. Compare estimates both in their own units and relative to their standard errors. A difference that looks small in raw units may matter for a parameter with a very small standard error, while a larger raw difference may be inconsequential for a poorly identified parameter. Neither comparison should replace subject-matter judgment.

If the same likelihood, data, and constants are used, the maximized log likelihoods should approach one another as the integral approximation becomes more accurate. They need not increase monotonically across settings: each approximation defines a slightly different numerical objective, and the optimizer can locate different approximate maxima. A log-likelihood difference should therefore be read alongside parameter and uncertainty changes, not as a standalone pass/fail test.

Why Estimates and Standard Errors Can Behave Differently

The numerical approximation affects both the likelihood value and its shape. Estimates depend on where the approximated likelihood is maximized. Standard errors from the observed information depend on its local curvature at that maximum. Consequently, estimates can appear stable while standard errors remain sensitive to integration settings, or the reverse can occur.

A basic justification for checking curvature is the local quadratic approximation to the log likelihood near an interior maximum \(\hat{\theta}\). If the observed information is \(J(\hat{\theta})=-\partial^2\ell(\theta)/\partial\theta\partial\theta^{\mathsf T}\) evaluated at the maximum, the usual covariance estimate is \(J(\hat{\theta})^{-1}\). A small integration error in \(\ell\) does not guarantee a small error in its second derivatives: differentiation can amplify small numerical fluctuations. Thus, standard errors should be compared directly, and a suspicious Hessian should be recomputed using a sufficiently accurate likelihood evaluation.

Diagnostics and Failure Modes

A comparison is informative only if differences can reasonably be attributed to integration. Check that every fit used the same data and likelihood definition, that optimizer convergence is credible, and that no parameter has moved to a boundary or a different local optimum. If the likelihood changes substantially but parameters barely move, inspect estimates not shown in a short summary, especially variance and association parameters. If estimates change substantially while likelihoods look nearly identical, a flat ridge or weak identification may be present; a similar likelihood does not imply a unique or precise estimate.

Numerical settings can also expose implementation difficulties. Higher-order quadrature may cost much more, particularly as the random-effects dimension grows. A rule that works well for one random-effects distribution or data configuration may perform poorly for another. Large changes in standard errors can indicate that the Hessian is sensitive to numerical noise, that the likelihood is weakly curved, or both. Re-evaluating the likelihood at the reported estimate with tighter settings, and recomputing derivatives consistently, helps separate these possibilities.

Choosing Among Approaches

Product Gauss–Hermite quadrature is a deterministic option with a transparent node-count control, but the number of node combinations grows rapidly with random-effects dimension. Adaptive Gauss–Hermite quadrature shifts and scales nodes to better match the subject-specific integrand; this can improve efficiency, though results still require a settings comparison. Other deterministic integration or Monte Carlo approaches may be useful when the dimension makes a product rule impractical, but they introduce their own controls, such as tolerances or simulation variability. A Laplace approximation can be computationally attractive, but it is an approximation of a different form and should not be treated as numerically interchangeable with a highly resolved quadrature fit.

For a focused sensitivity check, first compare progressively tighter settings within one integration approach. If results do not stabilize or the dimension makes that approach impractical, compare a suitable alternative and report the method used. Agreement across two methods is reassuring, but both can share modeling or optimization problems. Conversely, disagreement does not identify which result is correct; it signals a need for further numerical investigation.

Common Mistake and Practical Limits

A common mistake is to conclude, “The optimizer converged, so the integration must be accurate.” Optimizer convergence means that the algorithm met its stopping criteria for the objective it was given; it does not certify that the objective accurately approximates the intended likelihood. The correction is to assess integration sensitivity separately and to inspect estimates and standard errors, not only the optimizer message or log likelihood.

There is no universal numerical threshold that makes every comparison acceptable. A difference small relative to an estimate's standard error is useful context, but scientific consequences depend on the parameter scale, the intended use, and reporting precision. Nor does this exercise diagnose every source of instability: weak identification, model misspecification, sparse follow-up, and influential observations require their own assessment. Numerical stability supports interpretation of a fitted model; it does not establish that the model is scientifically adequate.