Evaluation of software tools for Bayesian estimation on population models: an update based on current software versions
Roberta Bartolucci, Silvia Maria Lavezzi, Elena Maria Tosca, Nicola Melillo, Silvia Grandoni, Elisa Borella, Lorenzo Pasotti, Giuseppe De Nicolao, Paolo Magni
Laboratory of Bioinformatics, Mathematical Modelling and Synthetic Biology, Department of Electrical, Computer and Biomedical Engineering, University of Pavia, via Ferrata 5, Pavia, I-27100, Italy
Objectives: Bayesian modelling based on Markov Chain Monte Carlo (MCMC) methods is acknowledged as a useful instrument in pharmacometrics. This work provides, 3 years after, an updated picture of a previous study [1], in which the performances of several software tools performing Bayesian estimation in a population context were compared in terms of efficiency and reliability of estimates, using as case studies an algebraic model and an ordinary differential equation (ODE) model.
Methods: NONMEM 7.4.1, WinBUGS 1.4.3 (with BlackBox Component Builder 1.5 and WBDiff interface), Stan 2.17, and JAGS 4.3 (with R packages Rstan and RJags) were selected for the present evaluation. In NONMEM, both BAYES and the newly added NUTS methods were evaluated, with and without mu referencing model implementation. The first model selected as a case study was a Poisson count model, describing a clinical trial of an anticonvulsant therapy. Data of seizure attacks, covariates and priors were collected from a published study [2]. The second model was a two-compartment PK ODE model [3] for a Phase I study of a monoclonal antibody for epilepsy. Simulated data were generated via Simulx, and priors were defined based on literature data [4]. For each model and tool, the number of iterations in the burn-in and stationary phases was computed based on the Raftery algorithm [5] (raftery.diag function in R coda package), to obtain a number of independent samples able to describe the posterior distribution with sufficient precision. Posterior distributions were inspected and compared. The capability of the tools to obtain uncorrelated samples was evaluated through the K parameter, i.e. the number of consecutive samples of the generated Markov chain that have to be discharged to obtain a new chain of “independent” samples. The Effective Sample Size per execution time unit (ESS/T) was calculated as an efficiency index.
The study was conducted on a Windows 10 ASUS desktop PC, with Intel Core i5 3.30Ghz 4 cores and 8GB RAM.
Results: For the count model, the posterior distributions of all the tools were similar to the expected ones [2]. In terms of execution times, Stan was the fastest (28 s), followed by NONMEM NUTS and BAYES methods with mu referencing (33 and 38 s), and WinBUGS and JAGS (~100 s). Estimation in NONMEM without mu referencing requested more than 4 min. The NUTS method, implemented both in NONMEM and Stan, showed the lowest K values, demonstrating its capability to generate almost independent samples. As for ESS/T, NONMEM NUTS and BAYES methods with mu referencing showed better performance with respect to the other tools; compared to BAYES, NUTS slightly improved both the efficiency and the estimation results.
For the ODE model, all the tools completed the estimation process, except for NONMEM NUTS method without mu referencing (due to convergence issues), and JAGS (it does not include an ODE solver). No tool was able to recover the expected posterior distributions [4] for all model parameters: variances of residual variability terms were always over/under-estimated. NONMEM BAYES and NUTS methods with mu referencing provided the most reliable results, whereas Stan estimated biased and highly skewed posterior distributions. The lowest execution times were obtained with NONMEM NUTS and BAYES methods with mu referencing (8 and 10 min), followed by WinBUGS (19 min), Stan (63 min), and NONMEM BAYES method without mu referencing (4.8 h). Again, the NUTS method displayed the lowest K values, followed by NONMEM BAYES method with mu referencing, WinBUGS, and NONMEM BAYES method without mu referencing. In terms of ESS/T, the best performances were obtained with NONMEM NUTS and BAYES methods with mu referencing for fixed effects, whereas WinBUGS showed higher ESS/T for random effects. NONMEM BAYES method without mu referencing and Stan showed always considerably lower ESS/T.
Conclusions: Based on the tested count and ODE models, according to the computed posterior distributions and ESS/T, NONMEM with mu referencing (in particular with the NUTS method) appears to be the best choice for reliable and efficient Bayesian estimation. It is followed by WinBUGS (or equivalently JAGS but for algebraic models only), which still represents a more flexible choice for Bayesian analysis compared to NONMEM (despite the significant advances of its version 7.4). Finally, note that significant improvements were made for ODE models Bayesian estimation in the latest release of Stan, even if estimates are not as satisfactory as for other tools.
Acknowledgments: The authors thank Dr. Robert Bauer for his technical support on NONMEM.
References:
[1] E. Borella et al. Evaluation of software tools for Bayesian estimation on population models with count and continuous data. Proceedings of the 2015 PAGE meeting, June 2-5 Hersonissos, Greece
[2] http://www.openbugs.net/Examples/Epil.html
[3] F. Strimenopoulou et al. Bayesian non-linear PK modelling applied to dose escalation studies using WinBUGS. Proceedings of the Bayes 2012 meeting, May 9‐11, Basel, Switzerland
[4] R. Lledo‐Garcia et al. Dose escalation studies for mAb: prior distributions selection and software comparison. Proceedings of the PAGE meeting (2012), June 5‐8, Venice, Italy
[5] A.E. Raftery and S.M. Lewis. Implementing MCMC. Markov Chain Monte Carlo in Practice (1996), London: Chapman and Hall, pp. 115-130
