Validating Quadrature Accuracy

Why Quadrature Needs Validation

A selection-model likelihood often contains integrals over a latent variable. Numerical quadrature replaces each integral with a finite weighted sum, so the likelihood being optimized is an approximation to the intended likelihood. A fit can appear to converge while still depending materially on that approximation. Validation asks whether increasing the numerical accuracy changes the likelihood, the estimates, or the inferential conclusions enough to matter.

The earlier material on Numerical Integration in Selection-Model Likelihoods introduced quadrature as a way to approximate these integrals. Here the emphasis is on checking the approximation. A single node count, successful optimizer termination, or a plausible-looking estimate is not evidence by itself that numerical integration is adequate.

For illustration, suppose the contribution of subject \(i\) to an observed-data likelihood involves a latent standard-normal variable \(Z\). Write that contribution as \(I_i(\theta)=\int g_i(z;\theta)\phi(z)\,dz\), where \(\theta\) is the model parameter vector, \(g_i\) is the relevant observed-data contribution conditional on \(Z=z\), and \(\phi\) is the standard normal density. The log likelihood is \(\ell(\theta)=\sum_{i=1}^n \log I_i(\theta)\). A quadrature rule with \(K\) nodes approximates each integral and therefore defines an approximate log likelihood, \(\ell_K(\theta)\). The fitted estimate is the maximizer of that approximate function.

For a standard-normal expectation, the nodes and weights are scaled so that the weighted sum approximates integration against \(\phi\), with the weights summing to one. If the latent variable has a different normal mean or variance, the transformation to the quadrature scale must reflect that distribution. A validation comparison is meaningful only when the likelihood, data, parameter constraints, and other model settings are held fixed while the integration calculation is changed.

What Should Be Compared?

The most direct check is to refit the model using a sequence of increasing node counts. Record more than the final log likelihood: compare parameter estimates, relevant standard errors or profile-likelihood results, and the quantities used in the scientific conclusion. Quadrature approximations need not change monotonically with \(K\), so a small difference at one step is not enough to establish stability. Several successive comparisons are more informative.

A second check uses a different integration rule, such as adaptive integration, at the same fitted parameter values or in a separate refit. Agreement between methods that have different numerical error patterns is useful corroboration. It is not a proof of exactness: two calculations can share a misspecification or a common implementation error. In either comparison, use unrounded output for calculations and report enough digits to make the comparison interpretable.

The relevant tolerance depends on the analysis. A tiny change in the log likelihood can be important when a profile is flat, a constraint is near a decision boundary, or an interval endpoint moves appreciably. Conversely, a change in the last displayed digit may have no practical consequence if estimates and inferential conclusions are unchanged. Numerical tolerances should be judged against the intended use, not selected as universal constants.

Why Small Integration Errors Can Matter

A short approximation explains why checking only the integral values can be misleading. Suppose a quadrature calculation gives \(I_i+\delta_i\) instead of \(I_i\), with \(|\delta_i/I_i|\) small and \(I_i>0\). A first-order expansion of the logarithm gives \(\log(I_i+\delta_i)=\log I_i+\delta_i/I_i+O((\delta_i/I_i)^2)\). Thus the approximate log likelihood can accumulate small relative integration errors across subjects. Moreover, if the likelihood is flat in some parameter direction, a small numerical perturbation can move the maximizer noticeably.

The effect on an estimate can also be seen from the score equation. Let \(\widehat{\theta}\) maximize the exact log likelihood locally, and suppose the numerical score differs from the exact score there by \(e\). If the Hessian \(H\) is nonsingular and the perturbation is small, a first-order expansion of the numerical score around \(\widehat{\theta}\) yields \(0\approx H\Delta\theta+e\), so \(\Delta\theta\approx-H^{-1}e\). When curvature is weak, \(H^{-1}\) can be large. This is why likelihood ridges and weak identification make careful numerical validation especially important.

Check the Estimates and the Inference

Numerical stability is an inferential question as well as a likelihood question. For each node count, inspect the parameters that matter to the scientific conclusion, including any selection parameter that is weakly identified. Compare standard errors, profile-likelihood curves or endpoints, and test statistics as appropriate. If the conclusion relies on whether an interval includes a particular value or a statistic crosses a prespecified cutoff, check whether the numerical changes could reverse that conclusion.

The numerical calculation does not determine which inferential method or calibration is appropriate. Use the intended method and assumptions consistently across quadrature comparisons. In particular, do not assume that a familiar reference distribution applies merely because the likelihood values are stable; selection-model constraints and weak identification can require careful calibration.

Compare with an Independent Integration Rule

Increasing node counts tests convergence within a chosen family of quadrature rules. A second rule provides a different diagnostic. For example, a fixed-node Gaussian rule and an adaptive rule may respond differently to a sharply concentrated or tail-heavy integrand. Compare them first at the same parameter values to isolate integration differences; then, if needed, compare separately optimized fits to assess their effect on estimates.

Diagnostics, Alternatives, and Limitations

Evidence that quadrature may be inadequate includes estimates or interval endpoints that continue to drift as nodes increase, unstable score or Hessian calculations, and changes in a test or interval conclusion across reasonable numerical settings. Large changes between two node counts are a clear warning. Small changes do not automatically clear the calculation, especially when a parameter is weakly identified or the decision is close to a cutoff.

Separate integration error from optimization error. Refit with suitable convergence checks at each node count; otherwise, an apparent quadrature effect may be an optimizer stopping effect. Conversely, optimizer convergence at each setting does not show that the integrated likelihood is accurate. Keep parameter constraints and starting strategies comparable, and check that the fitted values are not merely changing because one run reached a different local maximum.

Alternative integration rules can help, but they have limits. Adaptive methods may be more effective when the integrand is concentrated in a region that a fixed rule represents poorly, yet they also have tolerances and implementation choices that require checking. Simulation-based integration can be useful in some settings but may introduce Monte Carlo noise into the likelihood and its derivatives. No method should be treated as a numerical gold standard without diagnostics.

Quadrature validation cannot resolve model misspecification, weak identification, or uncertainty due to the observed sample. It assesses numerical approximation, not whether the selection model is scientifically adequate. Also, increasing nodes in multiple latent dimensions can be computationally expensive. The practical goal is not an unattainable declaration of exactness; it is a documented level of numerical stability adequate for the intended estimates and conclusions.