[R] Translating lme code into lmer was: Mixed effect model in R

Douglas Bates bates at stat.wisc.edu
Fri Oct 20 21:05:29 CEST 2006 

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.
>

More information about the R-help mailing list