Gauss-Newton algorithm for parameter estimation of nonlinear mixed-effects models
Felix Jost (1) and Sebastian Sager (1)
(1) Institute of Mathematical Optimization, Faculty of Mathematics, Otto-von-Guericke University Magdeburg, Germany
Objectives: Solving (nonlinear) least squares problems (LS problems) with a Gauss-Newton (GN) algorithm is a state-of-the-art approach and the algorithm's benefits are mathematically proven, since GN algorithms provide reliable estimates, not affected by residuals compared to Quasi-Newton (QN) and Newton algorithms.
In NONMEM, the individual parameters ETA are estimated by a GN algorithm, whereas the population parameters THETA, OMEGA and SIGMA are estimated by solving an optimization problem using a variable metric method [1] (QN Method [2]).
The resulting Hessian matrix contains residual-dependent terms, which can drive the iterative optimization procedure into solutions with undesired properties, such as solutions which are strongly influenced by large residual terms in the Hessian matrix.
For this reason, we propose a GN algorithm for the first-order (FO) and first-order conditional estimation (FOCE) approximation of parameter estimation problems for nonlinear mixed-effects models. The algorithm is implemented as a prototype in the software CasADi [3]. CasADi is a symbolic framework for automatic differentiation and numerical optimization with interfaces e.g. to ODE (CVODES and IDAS) and NLP (IPOPT, SNOPT, BLOCKSQP) solvers.
Methods: In a first step, we compared Hessian matrices from a Newton and GN algorithm for LS problems in terms of the residuals' influence on the Hessians' positive definiteness. Therefore, we solved 1000 LS problems for a one-compartment pharmacokinetic model with first-order absorption and oral administration of 320 mg theophylline published in [2].
Based on the results from the first step, we propose a GN algorithm for parameter estimation of nonlinear mixed-effects models. We reimplemented NONMEM's FO and FOCE methods [1,2,4] as a prototype in CasADi and modified the standard algorithm towards a GN algorithm. As Hessian matrix, the approximated population Fisher Information Matrix from [5] is implemented and passed to the interior point algorithm IPOPT [6] interfaced in CasADi [3].
We validated our FO and FOCE reimplementation by comparing parameter estimates, objective function values and Hessian matrix values from NONMEM's FO and FOCE method for two examples.
We used the theophylline model [2] consisting of 12 patients with 10 measurements each (interindividual variability implemented on volume of distribution, elimination rate and absorption rate constant, combined error model for residual variability and full OMEGA matrix) and a one-compartment dummy example with single intravenous bolus administration for 10 subjects with 2 measurements each (interindividual variability implemented on eliminate rate constant, additive error model for residual variability) from [4].
Finally, we tested our new GN FO method by comparing parameter estimates, objective function values and the Hessian matrix dependence on residuals for the theophylline model.
Our GN FOCE method was tested by estimating parameters for the dummy example and comparing the objective function values from our GN and QN FOCE algorithm with those obtained from NONMEM.
Results: The results of the 1000 LS problems demonstrate the superiority of a GN algorithm with respect to the Hessian matrices' independence of residual-containing terms.
Our reimplementation of NONMEM's FO and FOCE method was successfully validated in terms of almost identical objective function values with respect of the two published examples.
The benefit achieved from our GN algorithm applied on the FO approximation is shown using the theophylline example, in which the Hessian matrix from the GN algorithm compared to the Newton algorithm is not affected by large residuals within the solution.
For the dummy example our QN and GN algorithms converge to the same solution which slightly deviates, with a 0.36 smaller objective function value, from the reference solution from NONMEM.
NONMEM: sum( OBJ_i ) = sum( 8.56, -1.04, -2.04, 1.78, -2.07, -2.27, -1.69, -1.93, -2.18, -1.37 ) = -4.25
QN/GN: sum( OBJ_i ) = sum( 9.03, -1.13, -1.99, 1.88, -2.24, -2.37, -1.92, -2.02, -2.34, -1.51 ) = -4,61
Conclusions: GN algorithms are a promising alternative to QN and Newton algorithms for parameter estimation of nonlinear mixed-effects models, as the Hessian matrix, and thus the computed solution, is not affected by random residual terms. Therefore, we provide a prototype reimplementation of NONMEM's FO and FOCE method in CasADi, which allows to choose between a Newton, QN and GN algorithm.
References:
[1] Bauer R., NONMEM7 Technical Guide, 2010.
[2] Kim, M., Yim, D. and Bae, K., R-based reproduction of the estimation process hidden behind NONMEM Part 1: first-order approximation method, Translational and Clinical Pharmacology, 2015.
[3] Andersson, J., A General-Purpose Software Framework for Dynamic Optimization. Ph.D. Thesis, Arenberg Doctoral School, KU Leuven, Leuven, Belgium, 2013.
[4] Wang, Y., Derivation of various NONMEM estimation methods, Journal of Pharmacokinetics and pharmacodynamics, 2007.
[5] Nyberg, J., Karlsson, M. O. and Hooker, A. C., Simultaneous optimal experimental design on dose and sample times, Journal of Pharmacokinetics and Pharmacodynamics, 2009.
[6] Wächter, A. and Biegler, L., On the Implementation of an Interior-Point Filter Line-Search Algorithm for Large-Scale Nonlinear Programming. Mathematical Programming, 2006.