Diagnosing Numerical Integration Failure in NLME Likelihoods

Why Numerical Integration Needs Its Own Diagnostics

In an NLME model, the likelihood for one subject is obtained by integrating over that subject's unobserved random effects. An optimizer can return a converged parameter estimate even when the likelihood values used along the way are inaccurate. The resulting estimates, standard errors, or comparisons may then reflect integration error rather than the data and model. Diagnosis therefore asks not only whether a computation completed, but whether the integral is being represented reliably.

This tutorial focuses on failures in importance sampling and quadrature: inadequate proposal coverage, multiple regions of substantial integrand mass, ill-conditioned curvature, insufficient quadrature order, and numerical overflow or underflow. The key is to examine the integration itself at fixed model parameters before interpreting changes in the fitted model.

Definition: For subject $i$, let $y_i$ denote the observations, $\eta_i$ the random effects, $f_i(y_i\mid\eta_i;\theta)$ the conditional data density, and $g(\eta_i;\Omega)$ the random-effects density. The subject-specific marginal likelihood is $$ L_i(\theta,\Omega)=\int f_i(y_i\mid\eta_i;\theta)g(\eta_i;\Omega)\,d\eta_i. $$ If subjects are independent conditional on the population parameters, the log likelihood is $\ell(\theta,\Omega)=\sum_i\log L_i(\theta,\Omega)$. Numerical integration failure occurs when the computed approximation to one or more $L_i$ is materially inaccurate or unstable for the intended calculation.

Importance sampling and quadrature approximate the same integral in different ways. For importance sampling, draws from a proposal density $q_i(\eta)$ are weighted by $f_i(y_i\mid\eta;\theta)g(\eta;\Omega)/q_i(\eta)$. Its validity requires $q_i(\eta)>0$ wherever the integrand is nonzero, and the integral must be finite. Quadrature instead evaluates the integrand at a prescribed set of nodes and combines the evaluations with weights. Its accuracy depends on how well those nodes and weights resolve the integrand; a finite quadrature rule is not automatically exact for a nonlinear mixed-effects likelihood.

Diagnose Importance-Sampling Coverage

A proposal can be mathematically valid yet practically poor. If it assigns very little probability to a region that contributes substantially to $L_i$, most sampled weights will be small and a few may be extreme. This can cause high Monte Carlo variability and make different runs yield noticeably different likelihood estimates. The effective sample size (ESS) of nonnegative weights $w_{ir}$, $r=1,\ldots,R$, is

\operatorname{ESS}_i=\frac{\left(\sum_{r=1}^R w_{ir}\right)^2}{\sum_{r=1}^R w_{ir}^2}.

ESS is at most $R$; it is close to $R$ when weights are similar and close to 1 when one weight dominates. It is a diagnostic of weight concentration, not a direct measure of integration accuracy. In particular, a high ESS cannot prove that the proposal covers a region it never sampled.

Worked Example: Detecting a Dominant Importance Weight

Suppose four hypothetical importance weights for one subject are $0.2,\ 0.2,\ 0.2,\ 19.4$. The estimated likelihood is the mean weight: $$ \widehat L_i=\frac{0.2+0.2+0.2+19.4}{4}=\frac{20}{4}=5.00. $$

The ESS is $$ \operatorname{ESS}_i=\frac{20^2}{0.2^2+0.2^2+0.2^2+19.4^2} =\frac{400}{0.12+376.36} =\frac{400}{376.48}\approx1.063. $$

Although four draws were used, the weight concentration is comparable to having only about 1.06 equally weighted draws. The estimate is therefore fragile: a new set of draws could produce a materially different result. The ESS alone does not tell us how far $5.00$ is from the true likelihood; replicate runs or a better-covered proposal are needed to assess that.

A useful justification clarifies why repeating the calculation matters. If $\eta_{ir}$ are independent draws from $q_i$, then

\widehat L_i=\frac{1}{R}\sum_{r=1}^R \frac{f_i(y_i\mid\eta_{ir};\theta)g(\eta_{ir};\Omega)} {q_i(\eta_{ir})}.

Under the support and integrability conditions above, the expectation of each summand is $$ E_{q_i}\left[\frac{f_i(y_i\mid\eta_i;\theta)g(\eta_i;\Omega)}{q_i(\eta_i)}\right] =\int \frac{f_i(y_i\mid\eta_i;\theta)g(\eta_i;\Omega)}{q_i(\eta_i)} q_i(\eta_i)\,d\eta_i=L_i. $$ Thus the likelihood estimator is unbiased. But this fact does not guarantee low variance for a finite sample. Moreover, taking logs changes the property: since $\log$ is concave, Jensen's inequality gives $E[\log\widehat L_i]\leq\log L_i$ when the expectations exist. An unbiased likelihood estimate is not an unbiased log-likelihood estimate.

Multiple Modes and Ill-Conditioned Curvature

A single Gaussian proposal centered at a conditional mode can miss another separated region that contributes appreciably to the integral. This risk increases when optimization from different starting points finds distinct modes or when simulated weights show large clusters associated with different regions of random-effects space. A single-mode approximation may then appear stable within its own region while failing to represent the full integral.

Curvature creates a related difficulty. Near a mode, the negative Hessian of the log integrand describes local curvature. If it is nearly singular, or has eigenvalues spanning very different scales, a Gaussian proposal based on that curvature can be extremely wide in some directions and narrow in others. Its orientation or scale may then poorly match the integrand. Curvature is local evidence, not proof that the surrounding distribution is approximately Gaussian.

Worked Example: A Proposal That Misses a Second Mode

Consider a hypothetical one-dimensional integrand with two narrow peaks near $\eta=-3$ and $\eta=3$, each contributing about half the total integral. A Gaussian proposal $q(\eta)=N(-3,0.5^2)$ is centered on the left peak. At the right peak, its density relative to its value at $-3$ is $$ \frac{q(3)}{q(-3)} =\exp\left(-\frac{(3-(-3))^2}{2(0.5)^2}\right) =\exp(-72), $$ which is approximately $5.38\times10^{-32}$.

The proposal therefore makes draws near the right-hand peak extraordinarily rare, even though that peak supplies about half of the integral. In a finite simulation, the computed likelihood may be based almost entirely on the left peak. The failure is indicated by sensitivity to starting mode, extreme weights if the second peak is sampled, or disagreement with an integration method that resolves both regions. A practical response is to compare proposals centered at the distinct modes or use a proposal that covers both; simply increasing the number of draws from this narrow proposal may be inefficient.

When diagnosing this problem, compare conditional-mode searches from multiple starting points and inspect whether the resulting modes have materially different objective values or locations. Also examine how estimates change under broader or mixture proposals. These checks are evidence about coverage, not a guarantee that every relevant region has been found. High-dimensional random effects make visual inspection particularly limited.

Quadrature Order and Numerical Stability

Quadrature failure often appears as substantial movement in the integral when the order or node refinement changes. For adaptive quadrature, compare estimates at successively higher orders while holding the model parameters fixed. A sequence that stabilizes is reassuring, but stability at low orders can be misleading if every rule misses the same narrow or separated region. Compare more than one integration strategy when the stakes warrant it, especially when the integrand is sharply peaked or asymmetric.

Worked Example: Checking Quadrature Convergence

At fixed hypothetical parameters, suppose an adaptive quadrature calculation gives the following log-likelihood contributions for one subject: $-12.418$ at order 5, $-12.401$ at order 9, $-12.400$ at order 15, and $-12.4001$ at order 25.

The change from order 15 to order 25 is $$ |-12.4001-(-12.400)|=0.0001. $$ The change from order 9 to order 15 is $0.001$, while the change from order 5 to order 9 is $0.017$. The decreasing differences suggest stabilization for this subject at these parameters. If a prespecified diagnostic tolerance were $0.001$ log-likelihood units, the order-15 and order-25 values would agree within that tolerance. This is a practical convergence check, not proof that the integral is exact or that the same order is adequate for every subject or parameter value.

If the sequence oscillates, drifts, or changes sharply for particular subjects, investigate their integrands, random-effects scales, and modes. Increasing order can help resolve a smooth, concentrated integrand, but it can also be computationally expensive and may not fix a rule poorly placed relative to separated modes.

Numerical overflow and underflow can corrupt an otherwise appropriate integration rule. If log contributions are $a_1,\ldots,a_R$, exponentiating them directly can overflow for large positive $a_r$ or underflow for very negative values. A stable log-sum-exp calculation uses $m=\max_r a_r$:

\log\left(\sum_{r=1}^R e^{a_r}\right) =m+\log\left(\sum_{r=1}^R e^{a_r-m}\right).

Because each $a_r-m\leq0$, the exponentials are bounded by 1. The same scaling principle applies to a weighted sum, by including the log weights in the $a_r$. Stable arithmetic prevents avoidable computational failure, but it cannot correct poor proposal coverage or inadequate quadrature.

Worked Example: Stabilizing a Log-Domain Sum

Suppose two hypothetical log contributions are $a_1=1000$ and $a_2=999$. Directly evaluating $e^{1000}+e^{999}$ may overflow in ordinary floating-point arithmetic. Here $m=1000$, so $$ \log(e^{1000}+e^{999}) =1000+\log(e^0+e^{-1}) =1000+\log(1+0.3679) \approx1000+0.3133 =1000.313. $$

The stable calculation returns a finite result without changing the mathematical sum. In an actual likelihood calculation, density normalizing constants and all integration weights must also be included; omitting them is a model or implementation error, not a numerical stabilization.

A Practical Diagnostic Sequence

Freeze the parameters.
Assess integration at fixed parameter values so changes in the result can be attributed to the integration settings rather than to refitting.

Check numerical validity.
Look for non-finite log contributions, failed factorizations, invalid covariance matrices, and direct exponentiation of extreme log values. Use stable log-domain calculations where appropriate.

Check coverage and repeatability.
For importance sampling, inspect ESS and weight concentration, compare independent runs, and assess sensitivity to proposal location, scale, or mixture coverage.

Refine and cross-check.
For quadrature, increase order and compare successive estimates. Where feasible, compare with a differently configured integration method and investigate subjects with the largest discrepancies.

Reassess the fit.
After diagnosing the integral, check whether parameter estimates and log likelihoods remain materially stable under the improved integration settings.

The alternatives have different failure signatures. Importance sampling provides a direct stochastic approximation and supports proposal adaptation, but it can have high variance and run-to-run fluctuation. Quadrature is deterministic for fixed nodes, but can require many evaluations and miss features not resolved by its nodes. Sparse-grid or other structured quadrature may reduce cost in some dimensions, but it too needs refinement checks. No method is universally reliable: agreement across carefully checked methods is stronger evidence than a single converged software run.

Common Mistakes and Limits

Treating optimizer convergence as integration validation. An optimizer can converge to a maximum of an inaccurate numerical objective. Check integration accuracy at and around the reported solution.
Reading ESS as proof of accuracy. ESS measures observed weight concentration. It cannot reveal a region that the proposal never sampled.
Increasing quadrature order without checking the integrand. Refinement is useful, but separated modes, extreme scaling, or implementation errors can remain.
Calling every discrepancy a numerical failure. Differences between model fits can also arise from parameter changes, coding choices, or model assumptions. Hold parameters fixed when testing the integral itself.

These diagnostics have limits. Monte Carlo variability can persist even with a well-designed proposal, while deterministic rules can agree accidentally. A chosen tolerance is a practical standard, not a universal theorem. Finally, numerical integration checks address computation under a specified model; they do not establish that the model is scientifically appropriate.

Key takeaway: Diagnose an NLME likelihood integral by checking numerical stability, proposal coverage or quadrature refinement, and sensitivity across runs or methods. Then verify that the fitted results are stable under the settings supported by those checks.

Check Your Understanding

Answer each question, then compare with the solutions.

With $R=100$ equally weighted importance draws, what is the ESS? What does that value establish, and what does it not establish?
Why can an importance-sampling estimate of $L_i$ be unbiased while the log of that estimate is not unbiased for $\log L_i$?
A quadrature log likelihood changes from $-8.20$ to $-8.19$ when the order increases. Calculate the absolute change and state what further check would be useful.
What is the main coverage risk of a narrow proposal centered at one of two separated integrand modes?

Solutions. 1. The ESS is $100^2/(100\times1^2)=100$. It establishes that the observed weights are equal; it does not establish that the proposal covered all important regions. 2. The likelihood estimate can have expectation $L_i$ under the support and integrability conditions, but concavity of the logarithm gives $E[\log\widehat L_i]\leq\log E[\widehat L_i]=\log L_i$. 3. The absolute change is $|-8.19-(-8.20)|=0.01$ log-likelihood units. Refine the order further and, if feasible, compare another integration configuration while holding parameters fixed. 4. It may almost never sample the other mode, despite that mode contributing substantially to the integral; assess alternate mode-centered or broader proposals and compare results.