Checking Estimates Across Quadrature Resolutions

Why Check Quadrature Resolution?

In a shared-random-effects joint model, each subject's likelihood contribution involves integrating over latent random effects. Numerical quadrature approximates that integral using a finite set of nodes and weights. The resulting likelihood is therefore an approximation, and its accuracy can affect the maximum-likelihood estimates as well as the reported log likelihood.

A fit that has converged computationally is not necessarily accurate with respect to integration. An optimizer may locate a precise maximum of a coarse quadrature approximation, while a finer approximation has a different maximum. Checking estimates across increasing quadrature orders is a practical way to look for this sensitivity. It assesses numerical stability; it does not establish that the statistical model is correct or that its parameters are well identified.

The Likelihood Approximation Being Checked

For illustration, suppose subject \(i\) has a shared random intercept \(b_i\sim N(0,\sigma_b^2)\). Let \(y_i\) denote that subject's longitudinal data, and let \((T_i,\delta_i)\) denote their observed event time and event indicator. Write \(f_Y(y_i\mid b_i;\theta)\) for the conditional longitudinal contribution and \(f_S(T_i,\delta_i\mid b_i;\theta)\) for the conditional survival contribution. Under the model's conditional-independence assumption given \(b_i\), the marginal likelihood contribution is

Here \(\theta\) collects all model parameters, including \(\sigma_b\), and \(\phi\) is the standard normal density. The integral is over the random effect's distribution, not over the observed data. A quadrature rule replaces it with a weighted sum. For Gauss-Hermite quadrature with \(K\) nodes \(x_q\) and weights \(w_q\), one common convention gives

This expression uses weights and nodes for the convention \(\int e^{-x^2}g(x)\,dx\approx\sum_qw_qg(x_q)\). Other software may use a different scaling or normalization, so comparisons should use the same stated quadrature convention. The fitted value at order \(K\), written \(\hat{\theta}_K\), maximizes the approximate log likelihood \(\ell_K(\theta)=\sum_i\log L_{i,K}(\theta)\), subject to the optimizer's convergence.

Why Increasing the Order Is Informative

The transformation \(b=\sqrt{2}\sigma_bx\) converts integration against a normal density into an integral with the \(e^{-x^2}\) weight used by Gauss-Hermite quadrature. As the number of nodes increases, this weighted sum can represent a smooth integrand more accurately. Because the integrand depends on the parameters, improving the integral approximation can change the likelihood surface and hence its maximizer.

This argument motivates a stability check, but it is not a guarantee that a particular finite order is accurate. The likelihood integrand may be sharply peaked, skewed, or concentrated in a tail. Also, Gauss-Hermite node sets at different orders are generally not nested: a higher-order rule is not simply the lower-order rule plus a few extra evaluations. Consequently, a sequence need not change monotonically at every step. What matters in practice is whether successive sufficiently fine resolutions give substantively similar answers.

A Fair Comparison Across Orders

A useful comparison changes the integration resolution while holding everything else fixed. In particular, use the same observations, likelihood definition, parameterization, data scaling, and treatment of censoring at every order. Fit each order to adequate optimizer convergence. Whenever feasible, use the previous order's estimate as one starting value, but also try additional starting values if the likelihood is difficult to optimize. A change in starting values alone can reveal that apparent quadrature sensitivity is actually a failure to find the same maximum.

It is often helpful to report both absolute and relative changes. For example, for parameter \(\theta_j\), one can report \(\hat{\theta}_{j,K_2}-\hat{\theta}_{j,K_1}\), and, when the scale makes it meaningful, the relative change \((\hat{\theta}_{j,K_2}-\hat{\theta}_{j,K_1})/\hat{\theta}_{j,K_1}\). Relative changes are not useful when the denominator is near zero; use an absolute, scientifically interpretable scale instead. A single universal numerical tolerance is not appropriate for every parameter or application.

Interpreting Changes and Diagnosing Problems

A small change in the log likelihood is useful evidence, but it is not sufficient by itself. Parameters can move while the log likelihood barely changes, especially when the likelihood is flat in one direction or when parameters are strongly correlated. Conversely, a small change in a parameter estimate does not establish that the likelihood approximation is accurate for every purpose. Check parameters that drive key scientific conclusions, including variance and association parameters, rather than only the fixed effects.

If estimates change substantially, first check that each fit really reached the intended maximum. Compare optimizer convergence messages, gradients or other optimization diagnostics, and results from alternative starting values. If multiple starts at a fixed quadrature order lead to different likelihoods, do not attribute the difference to quadrature alone. If fits are reproducible at each order but shift systematically as the order increases, the integration approximation is a more plausible source of sensitivity. In a model with several random effects, tensor-product quadrature can become expensive rapidly; that computational burden may require a different integration strategy, but does not remove the need to assess numerical accuracy.

The log likelihood is a sum over subjects. Its numerical change can be small relative to the total even when an individual subject's contribution changes more noticeably. For additional diagnosis, it can help to examine subject-level contributions or their distribution across resolutions, especially if a small number of subjects have extreme outcomes or influential random-effect integrands. That check complements, rather than replaces, comparisons of total likelihood and fitted parameters.

How This Check Relates to Other Approaches

Increasing ordinary Gauss-Hermite order is a direct resolution check for the same quadrature approximation. Adaptive quadrature instead shifts or rescales nodes to better cover the region where the integrand is concentrated. It can improve accuracy with fewer nodes, but its settings and convergence should also be checked. Monte Carlo or importance-sampling integration uses a different numerical approximation; repeated estimates or changes in sampling effort can help assess its simulation error, but it is not the same as a deterministic quadrature-order comparison. Laplace approximation is another alternative, based on a local approximation around an integrand's mode. Agreement between methods can be reassuring, but disagreement requires investigation rather than choosing whichever answer is preferred.

The relevant question is not whether one integration method is universally best. It is whether the method used gives sufficiently stable results for the model, data, and scientific quantities being reported, at a computational cost that is feasible. A comparison should state the method and resolutions used so that the numerical check can be understood and reproduced.

Common Pitfall and Limitations

Wrong reasoning: “The log likelihood changed by less than 0.01, so every parameter estimate is accurate.” This does not follow: a flat likelihood can have a similar value over a range of parameter values, and a small total change can conceal a meaningful shift in an important estimate. Correction: compare the log likelihood and each scientifically important parameter, using sensible scales and explicit tolerances.

A stability check is a numerical diagnostic, not a test of model fit, an assessment of sampling uncertainty, or proof of identifiability. Stable estimates can still be biased by model misspecification, and a poorly identified parameter can remain uncertain even when numerical integration is accurate. Numerical stability also depends on the likelihood being evaluated consistently and the optimizer locating the same maximum at each resolution.