Doran, Harold
2006-Oct-20 12:22 UTC
[R] Translating lme code into lmer was: Mixed effect model in R
This question comes up periodically, probably enough to give it a proper thread and maybe point to this thread for reference (similar to the 'conservative anova' thread not too long ago). Moving from lme syntax, which is the function found in the nlme package, to lmer syntax (found in lme4) is not too difficult. It is probably useful to first explain what the differences are between the two functions for some background. The lme function in the nlme package was originally designed to handle models with nested random effects. For example, students nested in a single school or employees nested in a single organization. It was possible to fit models with crossed random effects using lme, but it was inefficient and the algorithms were not optimized for such situations. The lmer function, on the other hand, can handle models with crossed (or partially crossed) random effects as well as models with nested random effects. Doug Bates has also significantly improved the computational algorithms such that the speed by which very large data problems can be resolved is impressively fast. A model fit using lmer is typically of the form lmer(response ~ covariate(s) + (covariates|group), data) Notice that the random terms in lmer are always of the form (a|b) denoting that a is a random effect conditional on the grouping factor b. Also notice that random effects terms must be enclosed within parentheses. In lme the random term is of the form random = ~ covariate | group To illustrate the difference, here is an example using the sleepstudy data where the first is the lmer model and the second is the lme model. Parameter estimates from both models are the same. fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy, method='ML') fm2 <- lme(Reaction ~ Days, random=~Days|Subject, sleepstudy, method='ML') Now, assume you want to extend this to include additional levels of random variation. The egsingle data set in the mlmRev package has repeated measures on students and those students are nested in schools. We can fit a model with random intercepts and slopes for students and students within schools. The first model is fit using lme and the second is fit using lmer. fm3 <- lme(math ~ year, random=~year|schoolid/childid, egsingle) fm4 <- lmer(math ~ year +(year|schoolid:childid) +(year|schoolid), egsingle, control=list(gradient = FALSE, niterEM = 0)) Both result in parameter estimates that are exactly the same. The newest release of lme4 allows for an even easier transition between lme and lmer. Model syntax in lmer can be specified as: fm5 <- lmer(math ~ year +(year|schoolid/childid), egsingle, control=list(gradient = FALSE, niterEM = 0)) Notice how the grouping structure is used in this case, which resembles the lme syntax for random terms. The reason that the control statement is used is because Doug Bates told me to. I hope he pipes in and explains it, but from what I recall, the way he wrote the program has some fancy math. As it turns out, it is better to skip some of that math by turning off the EM iterations. Just for documention, I mentioned the CPU time is much faster using lmer than lme.> system.time(fm3 <- lme(math~year, random=~year|schoolid/childid,egsingle)) [1] 63.08 0.03 63.69 NA NA> system.time(fm4 <- lmer(math ~ year +(year|schoolid:childid)+(year|schoolid), egsingle, control=list(gradient = FALSE, niterEM 0))) [1] 2.22 0.00 2.22 NA NA On my Windows machine, the first model fit using lme took just over a minute whereas the same model using lmer was estimated in 2.22 seconds. There are other examples of how to use the lmer function in the mlmRev package found using vignette("MlmSoftRev") This is anything but comprehenssive, but I hope others can chime in and share their experiences, or correct any errors I may have made above. Best, Harold> -----Original Message----- > From: r-help-bounces at stat.math.ethz.ch > [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Lina Jansen > Sent: Friday, October 20, 2006 4:58 AM > To: R-help at stat.math.ethz.ch > Subject: Re: [R] Mixed effect model in R > > Thanks for the helping links. Now, I worked out that I have > to use the lme4 package (with the lmer function) for my > analysis. But now I do not understand the input to the lmer function. > > In the lme function (of the nlme package) the correct input > would in my case > be: > > lme(fixed=Ac_LC~cond_ind,random=~img_cond|sub_ind/cond_ind) > > but also after reading the help and the R news I do not > understand the formula I have to use for the lmer function. > Could someone help me translating the lme input to a lmer > input? And does someone know of a good explanation of the > kinds of formulas you can input? In the books they only > explain the lme input. > > Thanks, > Lina > > 2006/10/17, Stefan Grosse <singularitaet at gmx.net>: > > > > Please always reply to the list as well as there always might be > > someone faster/better answering (or it could be that I am wrong, so > > someone might correct me) > > > > Indeed Pinheiro/Bates assume gaussian error terms... but I am not > > really sure whether you meant that with "non normally distributed > > respond variable" resp. "with non-normal data" > > > > however: > > "/ Mixed-effects models: / The recommended nlme > > <http://cran.r-project.org/src/contrib/Descriptions/nlme.html> > > package, associated with Pinheiro and Bates, / > Mixed-Effects Models in > > S and S-PLUS / (Springer, 2000), fits linear and nonlinear > > mixed-effects models, commonly used in the social sciences for > > hierarchical and longitudinal data. Generalized linear > mixed-effects > > models may be fit by the glmmPQL function in the MASS > package, and by > > the lmer function in the Matrix > > <http://cran.r-project.org/src/contrib/Descriptions/Matrix.html> > > package (related to the lme4 > > <http://cran.r-project.org/src/contrib/Descriptions/lme4.html> > > package, which largely supersedes nlme > > > <http://cran.r-project.org/src/contrib/Descriptions/nlme.html> for / > > linear / mixed models). Also see the lmeSplines > > > <http://cran.r-project.org/src/contrib/Descriptions/lmeSplines.html> > > and lmm > <http://cran.r-project.org/src/contrib/Descriptions/lmm.html> > > packages." [ > > http://cran.r-project.org/src/contrib/Views/SocialSciences.html ] > > > > Lina Jansen schrieb: > > > > > > > > > 2006/10/17, Stefan Grosse <singularitaet at gmx.net > > > <mailto:singularitaet at gmx.net>>: > > > > > > Interesting packages for you might be the nlme and > lme4 packages > > > and as > > > a book Pinheiro/Bates, "Mixed-Effects Models in S and S-Plus" > > > > > > > > > Thank you for the answer. I am always unsure concerning the > > > non-normality. Can I use the nlme and lme4 with non-normal data? > > > First, I thought they would work like an ANOVA but with > random and > > > fixed effects. > > > > > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >
Douglas Bates
2006-Oct-20 19:05 UTC
[R] Translating lme code into lmer was: Mixed effect model in R
On 10/20/06, Doran, Harold <HDoran at air.org> wrote:> This question comes up periodically, probably enough to give it a proper > thread and maybe point to this thread for reference (similar to the > 'conservative anova' thread not too long ago).> Moving from lme syntax, which is the function found in the nlme package, > to lmer syntax (found in lme4) is not too difficult. It is probably > useful to first explain what the differences are between the two > functions for some background.> The lme function in the nlme package was originally designed to handle > models with nested random effects. For example, students nested in a > single school or employees nested in a single organization. It was > possible to fit models with crossed random effects using lme, but it was > inefficient and the algorithms were not optimized for such situations.> The lmer function, on the other hand, can handle models with crossed (or > partially crossed) random effects as well as models with nested random > effects. Doug Bates has also significantly improved the computational > algorithms such that the speed by which very large data problems can be > resolved is impressively fast.> A model fit using lmer is typically of the form> lmer(response ~ covariate(s) + (covariates|group), data)> Notice that the random terms in lmer are always of the form (a|b) > denoting that a is a random effect conditional on the grouping factor b. > Also notice that random effects terms must be enclosed within > parentheses. In lme the random term is of the form> random = ~ covariate | group> To illustrate the difference, here is an example using the sleepstudy > data where the first is the lmer model and the second is the lme model. > Parameter estimates from both models are the same.> fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy, method='ML') > fm2 <- lme(Reaction ~ Days, random=~Days|Subject, sleepstudy, > method='ML')> Now, assume you want to extend this to include additional levels of > random variation. The egsingle data set in the mlmRev package has > repeated measures on students and those students are nested in schools. > We can fit a model with random intercepts and slopes for students and > students within schools. The first model is fit using lme and the second > is fit using lmer.> fm3 <- lme(math ~ year, random=~year|schoolid/childid, egsingle) > fm4 <- lmer(math ~ year +(year|schoolid:childid) +(year|schoolid), > egsingle, control=list(gradient = FALSE, niterEM = 0))> Both result in parameter estimates that are exactly the same. The newest > release of lme4 allows for an even easier transition between lme and > lmer. Model syntax in lmer can be specified as:> fm5 <- lmer(math ~ year +(year|schoolid/childid), egsingle, > control=list(gradient = FALSE, niterEM = 0))> Notice how the grouping structure is used in this case, which resembles > the lme syntax for random terms.> The reason that the control statement is used is because Doug Bates told > me to. I hope he pipes in and explains it, but from what I recall, the > way he wrote the program has some fancy math. As it turns out, it is > better to skip some of that math by turning off the EM iterations.Yes, that is correct. The lmer function optimizes the deviance of the linear mixed model using a shortcut known as profiling, which reduces the dimensionality of the optimization problem. In the article @Article{bates04:_linear, author = {Douglas M. Bates and Saikat DebRoy}, title = {Linear mixed models and penalized least squares}, journal = {Journal of Multivariate Analysis}, year = 2004, volume = 91, number = 1, pages = {1--17} } we derive expressions for the profiled deviance and for the gradient and the Hessian of the profiled deviance. We also derive expressions for an ECME iteration, a variation on the EM algorithm. The ECME expressions are closely related to the gradient expressions. For small data sets and simple models (e.g. models with a single grouping factor or with nested grouping factors) the ECME iterations can help to stabilize the initial iterations. Also the use of the gradient helps to stabilize the optimization of the profiled deviance using nlminb, which follows the ECME iterations. However, I have found that with large data sets and crossed or partially crossed grouping factors the ECME iterations and the gradient calculation result in much slower optimization because they both involve a calculation that is equivalent to inverting the Cholesky factor of the crossproduct matrix. The techniques that are used in lmer facilitate the calculation of the profiled deviance which only involves evaluating the Cholesky factor and solution of a small number of linear systems using this factor. Inverting the Cholesky factor is much, much slower than forming the factor. One of the items on my "To Do" list is to change the defaults for the gradient and niterEM control options according to whether there are multiple, non-nested grouping factors. Unfortunately there are many items ahead of that. For the time being it is best to remember that on large complex problems you should turn off the ECME iterations and the calculation of the gradient.> Just for documention, I mentioned the CPU time is much faster using lmer > than lme. > > > system.time(fm3 <- lme(math~year, random=~year|schoolid/childid, > egsingle)) > > [1] 63.08 0.03 63.69 NA NA > > > system.time(fm4 <- lmer(math ~ year +(year|schoolid:childid) > +(year|schoolid), egsingle, control=list(gradient = FALSE, niterEM > 0))) > [1] 2.22 0.00 2.22 NA NA > > On my Windows machine, the first model fit using lme took just over a > minute whereas the same model using lmer was estimated in 2.22 seconds. > > There are other examples of how to use the lmer function in the mlmRev > package found using > > vignette("MlmSoftRev") > > This is anything but comprehenssive, but I hope others can chime in and > share their experiences, or correct any errors I may have made above. > > Best, > > Harold > > > > -----Original Message----- > > From: r-help-bounces at stat.math.ethz.ch > > [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Lina Jansen > > Sent: Friday, October 20, 2006 4:58 AM > > To: R-help at stat.math.ethz.ch > > Subject: Re: [R] Mixed effect model in R > > > > Thanks for the helping links. Now, I worked out that I have > > to use the lme4 package (with the lmer function) for my > > analysis. But now I do not understand the input to the lmer function. > > > > In the lme function (of the nlme package) the correct input > > would in my case > > be: > > > > lme(fixed=Ac_LC~cond_ind,random=~img_cond|sub_ind/cond_ind) > > > > but also after reading the help and the R news I do not > > understand the formula I have to use for the lmer function. > > Could someone help me translating the lme input to a lmer > > input? And does someone know of a good explanation of the > > kinds of formulas you can input? In the books they only > > explain the lme input. > > > > Thanks, > > Lina > > > > 2006/10/17, Stefan Grosse <singularitaet at gmx.net>: > > > > > > Please always reply to the list as well as there always might be > > > someone faster/better answering (or it could be that I am wrong, so > > > someone might correct me) > > > > > > Indeed Pinheiro/Bates assume gaussian error terms... but I am not > > > really sure whether you meant that with "non normally distributed > > > respond variable" resp. "with non-normal data" > > > > > > however: > > > "/ Mixed-effects models: / The recommended nlme > > > <http://cran.r-project.org/src/contrib/Descriptions/nlme.html> > > > package, associated with Pinheiro and Bates, / > > Mixed-Effects Models in > > > S and S-PLUS / (Springer, 2000), fits linear and nonlinear > > > mixed-effects models, commonly used in the social sciences for > > > hierarchical and longitudinal data. Generalized linear > > mixed-effects > > > models may be fit by the glmmPQL function in the MASS > > package, and by > > > the lmer function in the Matrix > > > <http://cran.r-project.org/src/contrib/Descriptions/Matrix.html> > > > package (related to the lme4 > > > <http://cran.r-project.org/src/contrib/Descriptions/lme4.html> > > > package, which largely supersedes nlme > > > > > <http://cran.r-project.org/src/contrib/Descriptions/nlme.html> for / > > > linear / mixed models). Also see the lmeSplines > > > > > <http://cran.r-project.org/src/contrib/Descriptions/lmeSplines.html> > > > and lmm > > <http://cran.r-project.org/src/contrib/Descriptions/lmm.html> > > > packages." [ > > > http://cran.r-project.org/src/contrib/Views/SocialSciences.html ] > > > > > > Lina Jansen schrieb: > > > > > > > > > > > > 2006/10/17, Stefan Grosse <singularitaet at gmx.net > > > > <mailto:singularitaet at gmx.net>>: > > > > > > > > Interesting packages for you might be the nlme and > > lme4 packages > > > > and as > > > > a book Pinheiro/Bates, "Mixed-Effects Models in S and S-Plus" > > > > > > > > > > > > Thank you for the answer. I am always unsure concerning the > > > > non-normality. Can I use the nlme and lme4 with non-normal data? > > > > First, I thought they would work like an ANOVA but with > > random and > > > > fixed effects. > > > > > > > > > > [[alternative HTML version deleted]] > > > > ______________________________________________ > > R-help at stat.math.ethz.ch mailing list > > https://stat.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide > > http://www.R-project.org/posting-guide.html > > and provide commented, minimal, self-contained, reproducible code. > > > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >