Emmanuelle Comets (1), Audrey Lavenu (2), Marc Lavielle (3)
(1) INSERM UMR1137 IAME, Paris, France; University Denis Diderot, Paris, France; (2) University Rennes-I, Rennes, France; INSERM CIC 1414, Rennes, France; (3) INRIA Saclay, Popix University Paris Sud, France
Objectives: The saemix package for R [1] provides an implementation of the SAEM (Stochastic Approximation Expectation Maximization) algorithm. In the present paper, we assess the operational characteristics of saemix in terms of performance and scalability, using simulated data.
Methods: The SAEM algorithm is used to obtain maximum likelihood estimates of the parameters of nonlinear mixed effects models without any linearisation of the model [2,3]. The SAEM algorithm uses an EM algorithm, where the unknown individual parameters are treated as missing data, and replaces the usual E-step with a stochastic approximation step [4,5]. The saemix package makes use of the S4 classes in R to provide a user-friendly functions for estimation, diagnostics and summary.
We applied saemix on simulated data from Plan et al., who created it to compare the performance of various software [6]. A sigmoid Emax model was fitted to dose-response data simulated with relatively large curvature (gamma=3), under a rich and a sparse design (respectively 4 and 2 points per subject). We compared the performance of saemix with results from nlme [7] and nlmer from lme4 [8].
In parallel, we investigated the scalability of the algorithm by showing runtimes across models with varying numbers of parameters and across different designs.
Results: In the rich design, saemix was able to provide unbiased estimates of the population parameters, while both nlme and nlmer had trouble estimating ED50 and its variability. In the sparse design, the three algorithms exhibited bias but saemix showed the best performance. nlmer and to a lesser extent nlme exhibited convergence issues, especially for the sparse design.
As expected due to the stochastic nature of the algorithm, runtimes for saemix increased as a function of the number of random effects in the model and of the number of subjects in the dataset. We also implemented an ODE model using the standard R solver (deSolve) but this proved extremely slow.
Conclusions: The saemix package provides the SAEM algorithm for R users, as an alternative to linearisation-based algorithms, implemented in nlme [7], or quadrature methods, implemented in glmmML or lme4 [8]. Current development focuses on extending the capabilities of the package via new models, such as ODE systems or hidden Markov models [9], extended diagnostics, and automating covariate handling.
References:
[1] R Development Core Team (2006). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
[2] Delyon B, Lavielle M, Moulines E (1999). Convergence of a stochastic approximation version of the EM algorithm. Annals of Statistics 27:94-128.
[3] Kuhn E, Lavielle M (2004). Coupling a stochastic approximation version of EM with a MCMC procedure. ESAIM P&S 8:115-31.
[4] Kuhn E, Lavielle M (2005). Maximum likelihood estimation in nonlinear mixed effects models. Comput Stat Data Anal 49:1020-38.
[5] Lavielle M (2014). Mixed Effects Models for the Population Approach – Models, Tasks, Methods and Tools. Chapman & Hall/CRC Biostatistics Series, CRC Press, Boca Raton, Florida, USA
[6] Plan E, Maloney A, Mentré F, Karlsson M, Bertrand J (2012). Performance Comparison of Various Maximum Likelihood Nonlinear Mixed-Effects Estimation Methods for Dose-Response Models. AAPS J 14:420-32.
[7] Pinheiro J, Bates D, DebRoy S, Sarkar D, R Core team (2009). nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-96.
[8] Bates D, Maechler M (2010). lme4: Linear mixed-effects models using S4 classes. R package version 0.999375-34.
[9] Delattre M (2010). Inference in mixed Hidden Markov models and applications to medical studies. J Soc Franc Stat 151:90-105.
Reference: PAGE 25 (2016) Abstr 5775 [www.page-meeting.org/?abstract=5775]
Poster: Methodology - Estimation Methods