Generalized FOCE with Pumas
Andreas Noack (1), Patrick Kofod Mogensen (1), Joakim Nyberg (2), and Vijay Ivaturi (3)
(1) Julia Computing (2) Pharmetheus, Uppsala, Sweden (3) Center for Translational Medicine, University of Maryland School of Pharmacy, Baltimore, MD, USA
Introduction: Traditionally, the first-order estimation method (FOCE) [1] has been associated with the Gaussian error models. However, when FOCE is interpreted as an approximation of the quadratic coefficient in the Laplace approximation of the marginal likelihood function the method applies to any parametric error model. Specifically, it is possible to benefit from the cheaper FOCE approximation compared to a Laplace model using the actual Hessian when estimating Generalized Linear Models (GLM) type models such as logistic and Poisson regression. Furthermore, for data where the standard error of the dependent variable is proportional to the mean, as in many pharmacokinetic models, it is possible to use the FOCE method to fit Gamma models as proposed in [2] instead of relying on the traditional proportional (Gaussian) error model.
The generalized FOCE method has recently been added to the new JuliaPro package Pumas [3]. To our knowledge, this feature is not available in any other software package for PKPD modeling. Pumas is implemented in the Julia programming language and builds on top of many open source packages from the Julia ecosystem such as Distributions.jl [4], DifferentialEquations.jl [5], Optim.jl [6], and ForwardDiff.jl. By being implemented in Julia, Pumas is able to run at speed comparable to C programs while easily allowing for utilization of automatic differentiation which eliminates the requirement for symbolic derivative computations even for models based on ordinary differential equations.
Objectives: Examine the speed and accuracy of the newly added generalized FOCE support in the Pumas package on a PKPD model for hepatitis C virus (HCV) kinetics defined by five ordinary differential equations.
Methods: The HCV model in [7] was modified such the dependent PK variable was modeled with a Gamma error model. Seven subjects were simulated for four dosing intervals of 24 hours each and 11 samples were taken in each interval. The ODEs were solved with a hybrid integrator, AutoVern7(Rodas5()), from the DifferentialEquations.jl package. The inner optimization for the computation of the empirical Bayes estimates used a Newton trust region solver and the outer optimization used a quasi Newton (BFGS) solver both from Optim.jl. All gradients were computed with automatic differentiation through the ForwardDiff.jl package. The same dataset was then estimated with Pumas's FOCE method and full Hessian based LaplaceI (with interaction) method.
Results: The estimated parameters from the two estimations were similar and the number of iterations required for convergence was 42 and 40 for the FOCE and LaplaceI estimation respectively. The FOCE completed in 322 seconds (including 154 seconds of compilation) and the LaplaceI estimation completed in 3891 seconds (including 215 seconds of compilation).
Conclusions: The FOCE approximation of the marginal likelihood is much simpler than the Hessian based Laplace approximation both in terms of the required number of floating-point operations as well as the compiler overhead associated with automatic differentiation. The difference clearly showed in the timings of the two estimation methods where FOCE was 12 times faster than the Laplace estimation (22 times faster when excluding compile time). While it might be possible to reduce some of this difference with improvements to the Julia compiler, we suspect that FOCE will remain significantly faster than Laplace. Using the FOCE approximation in an adaptive Gauss-Hermite estimation procedure, where the Hessian based Laplace approximation has traditionally been used, is an interesting future direction of work.
References:
[1] Beal SL, Sheiner LB, Boeckmann A, Bauer RJ. NONMEM users guides. NONMEM Project Group, University of California, San Francisco. 1992.
[2] Nelder JA, Wedderburn RW. Generalized linear models. Journal of the Royal Statistical Society: Series A (General). 1972 May;135(3):370-84.
[3] Besançon M, Anthoff D, Arslan A, Byrne S, Lin D, Papamarkou T, Pearson J. Distributions. jl: Definition and modeling of probability distributions in the JuliaStats ecosystem. arXiv preprint arXiv:1907.08611. 2019 Jul 19.
[4] Rackauckas, C. & Nie, Q., (2017). DifferentialEquations.jl – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia. Journal of Open Research Software. 5(1), p.15. DOI: http://doi.org/10.5334/jors.151
[5] Mogensen PK, Riseth AN. Optim: A mathematical optimization package for Julia. Journal of Open Source Software. 2018 Apr 5;3(24).
[6] pumas.ai
[7] Nyberg J. et al. Methods and software tools for design evaluation in population pharmacokinetics–pharmacodynamics studies. Br J Clin Pharmacol. 2015 Jan; 79(1): 6–17.