Numerical Integration Error in ODE Likelihoods

Why integration accuracy matters in a likelihood

In an ODE-based pharmacometric model, predicted concentrations or amounts are obtained by solving differential equations. A likelihood calculation then compares those predictions with observations. If the numerical solution changes with solver settings, the likelihood can change too. An optimizer may respond to that change by moving parameters, even when the underlying scientific model and data have not changed.

Integration error is therefore a possible source of unstable likelihoods, gradients and estimates. It is distinct from residual error: residual error represents variation described by the observation model, whereas integration error is introduced while numerically approximating the model-predicted state. Making the observation-error model more flexible does not, by itself, correct an inaccurate ODE solution.

The basic diagnostic question is not simply whether a solver reports success. It is whether reasonable changes in solver settings or numerical method leave the quantities relevant to inference effectively unchanged. Those quantities include the state trajectory, observation predictions, objective-function values, gradients and parameter estimates.

Likelihoods computed from numerical predictions

Suppose, for illustration, that observations \(y_i\) at times \(t_i\) have independent Gaussian errors with known standard deviation \(\sigma\), and the model prediction is \(m_i(\theta)=g(x(t_i;\theta),\theta)\). The log likelihood, up to any terms that do not depend on \(\theta\), is

In computation, the exact state is generally unavailable, so the solver supplies \(\tilde{m}_i(\theta)=g(\tilde{x}(t_i;\theta),\theta)\). The calculated likelihood is then an approximation to the likelihood defined by the continuous model. A small state error does not automatically imply a negligible likelihood error: the effect depends on the observation scale, the residuals and how sensitive the predictions are to the parameters.

The usual numerical interpretation assumes that the ODE has a well-defined solution over the interval of interest and that the chosen solver approximates that solution adequately at the parameters being evaluated. For comparisons across parameters, the solver must also behave consistently enough that the computed objective represents the same model throughout the optimization. Solver tolerances are numerical controls, not statistical assumptions and not substitutes for a scientifically justified observation-error model.

Why small trajectory errors can alter gradients

Optimization uses changes in the computed objective to decide which direction to move parameters. If the numerical prediction has an error that varies with \(\theta\), its derivative contributes to the computed likelihood gradient. Near a well-fitting parameter, even a small prediction error can matter because the true likelihood gradient may be close to zero.

For the single-observation Gaussian model above, the score, or derivative of the log likelihood with respect to \(k\), is

This follows by differentiating \(-\{y-m(k)\}^2/(2\sigma^2)\) with respect to \(k\). In a numerical calculation, both the residual and the prediction derivative may differ from their continuous-model values. Consequently, a solver can produce a small-looking discrepancy in concentration while still noticeably perturbing the score when \(\sigma\) is small.

Check refinement, not just solver status

A solver’s successful return indicates that it completed according to its numerical rules. It does not establish that the solution is accurate enough for a particular likelihood or parameter estimate. A practical check is to tighten tolerances, refine the numerical resolution where appropriate, or compare a second suitable solver, then inspect quantities at identical parameter values as well as after refitting.

For adaptive solvers, relative tolerance controls error relative to state magnitude, while absolute tolerance provides a scale for error when a state is small. Their effects depend on the units and magnitudes of the states. A tolerance appropriate for a large state may be inadequate for a small one. Set numerical controls with state scaling in mind, and verify the result empirically by refinement rather than assuming a chosen tolerance guarantees inferential accuracy.

Why refinement works in this example

The direction of the Euler bias above can be justified directly. For \(0<k<N\), take logarithms of the Euler prediction:

All terms after \(-k\) are negative on this branch. Therefore the Euler endpoint is less than \(e^{-k}\) for \(k>0\). To match a fixed endpoint \(y=0.80\), the Euler fit must use a smaller \(k\) than the exact estimate in this example. Also,

This limit explains the convergence shown in the table. The argument applies to this decay equation and the nonnegative Euler branch; it should not be treated as a general claim about the direction of error for other models or numerical methods.

Diagnostics, alternatives and limitations

A useful numerical audit separates changes in the forward solution from changes caused by refitting. First evaluate the model at a fixed parameter vector under more than one suitably accurate numerical setting. Compare state trajectories at observation times and, where relevant, between observations. Then compare predicted observations, the likelihood or objective components, and gradients. Finally refit under the refined settings and compare estimates, uncertainty summaries and scientific conclusions.

Warning signs include trajectories that change visibly under refinement, large changes in predictions relative to the observation-error scale, likelihood or gradient changes that affect optimizer decisions, and fitted parameters that move materially as tolerances tighten. Inspect intermediate states as well as endpoints: the even-step Euler example shows how an implausible trajectory can return to a plausible endpoint. Save solver settings and report the refinement checks so that the numerical calculation can be reproduced.

There are several complementary ways to check an ODE likelihood:

• Use an analytic solution when available. It provides a direct reference for numerical predictions and can reveal error even when the numerical solver reports success.
• Refine a suitable numerical solver. Tighten tolerances or reduce step sizes and check convergence of the quantities used for inference.
• Compare an independent numerical method. Agreement between methods with different numerical behavior strengthens a diagnosis, although agreement alone is not proof of accuracy.

These checks have limits. Tighter tolerances increase computational cost and do not guarantee accuracy if the model is poorly scaled or the solver is being used outside its effective range. An independent method can share a common error. A close match in predictions at observation times may conceal differences in the trajectory between them. Moreover, if estimates remain unstable after integration has been shown to be adequate, other sources of instability should be investigated rather than continuing to tighten tolerances without a diagnostic purpose.