SAEMIX, an R version of the SAEM algorithm
Emmanuelle Comets (1), Audrey Lavenu (2), Marc Lavielle (3)
(1) INSERM UMR738, Paris, France; University Denis Diderot, Paris, France; (2) University Rennes-I, Rennes, France; INSERM CIC 0203, Rennes, France; (3) INRIA, Saclay, France
Objectives: The Stochastic Approximation Expectation Maximization (SAEM) algorithm has proven very efficient, quickly converging to the maximum likelihood estimators [1] and performing better than linearisation-based algorithms [2]. It has been implemented in the Monolix software [3] which has enjoyed increasingly widespread use over the last few years, more recently in the Statistics toolbox of Matlab (nlmefitsa.m), and is also available in NONMEM version 7 [4]. The objective of the present package was to implement SAEM in the R statistical software [5].
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. The log-likelihood for nonlinear mixed effect models is analytically intractable since it requires integration over the unknown individual parameters. The SAEM algorithm uses an EM algorithm [6], where the unknown individual parameters are treated as missing data, and replaces the usual E-step with a stochastic approximation step [7,8]. The missing parameters are simulated at each iteration via a MCMC procedure, which can be used after the algorithm has converged to obtain the conditional modes, the conditional means and the conditional standard deviations of the individual parameters.
We implemented the SAEM algorithm in R, in the package SAEMIX. It provides estimates of the population parameters and their standard errors for nonlinear mixed effect models expressed in analytical form.
Results: The algorithm was applied to several example datasets and models, and showed good performance. The package provides summaries of the results, individual parameter estimates, standard errors (obtained using a linearised computation of the Fisher information matrix) Wald tests for fixed effects, and a number of diagnostic plots, including VPC plots and npde [9]. The log-likelihood can be computed by three methods: a linearisation of the model, an importance sampling procedure, or a Gaussian quadrature. The diagnostic graphs can be tailored to the user's individual preferences by setting a number of options, and are easily exported to a file.
Conclusion: The SAEMIX package provides the SAEM algorithm for R users, as an alternative to linearisation-based algorithms, implemented in nlme [10], or quadrature methods, implemented in glmmML or lme4 [11,12]. The current version handles models in analytical form, with continuous or binary covariates.
References:
[1] Delyon, B., Lavielle, M., and Moulines, E. Convergence of a stochastic approximation version of the EM algorithm. Annals of Statistics 27 (1999), 94-128.
[2] Girard, P., and Mentré, F. A comparison of estimation methods in nonlinear mixed effects models using a blind analysis (oral presentation). PAGE, Pamplona (2005).
[3] Lavielle, M. MONOLIX (MOdèles NOn LInéaires à effets miXtes) User Guide. MONOLIX group, Orsay, France, 2010. URL: http://software.monolix.org/
[4] Beal, S., Sheiner, L.B., Boeckmann, A., & Bauer, R.J., NONMEM User's Guides. (1989-2009), Icon Development Solutions, Ellicott City, MD, USA, 2009.
[5] R Development Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2006.
[6] Dempster, A. P., Laird, N. M., and Rubin, D. B. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society Series B 39, 1 (1977), 1-38. With discussion.
[7] Kuhn, E., and Lavielle, M. Coupling a stochastic approximation version of EM with a MCMC procedure. ESAIM P&S 8 (2004), 115-31.
[8] Kuhn, E., and Lavielle, M. Maximum likelihood estimation in nonlinear mixed effects models. Computational Statistics and Data Analysis 49 (2005), 1020-38.
[9] Brendel, K., Emmanuelle Comets, Céline Laffont, Christian Laveille, and France Mentré. Metrics for external model evaluation with an application to the population pharmacokinetics of gliclazide. Pharmaceutical Research, 23:203649, 2006.
[10] Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D., and the R Core team. nlme: Linear and Nonlinear Mixed Effects Models, 2009. R package version 3.1-96.
[11] Broström, G. (2009). glmmML: Generalized linear models with clustering. R package version 0.81-6.
[12] Bates, D., Maechler, M. (2010). lme4: Linear mixed-effects models using S4 classes. R package version 0.999375-34.
