[R-sig-ME] Fwd: Timings with SAMM, lme4, nlme

Douglas Bates bates at stat.wisc.edu
Fri Feb 2 14:45:28 CET 2007


I also meant to cc: the list on this response.  Obviously I need to
work on my email skills.

---------- Forwarded message ----------
From: Douglas Bates <bates at stat.wisc.edu>
Date: Feb 1, 2007 12:39 PM
Subject: Re: [R-sig-ME] Timings with SAMM, lme4, nlme
To: Kevin Wright <kw.statr at gmail.com>


On 2/1/07, Kevin Wright <kw.statr at gmail.com> wrote:
> About a year and a half ago I did some comparison of model fitting with
> SAMM, lme4, and nlme.  Since Doug Bates put out a request for some recent
> timings, I am repeating/extending my comparisons.  In the interim, I have
> switched computers, so the timings of the samm function facilitate comparing
> the speed across the two computers.
>
> First, the results (view with fixed-width font)
>
>            Setup   T30   D620    T30   D620     D620    D620
> Model \ Function  samm   samm   lmer   lmer    lmer2     lme
> ----------------  ----   ----   ----   ----    -----   -----
> yrs|id + yrs|sch     f      f   30.0   kill     65
> yrs|id +   1|sch  10.5    7.8   13.7   kill     23
>   1|id + yrs|sch   9.4    9.5   11.5   kill      7.3
>   1|id +   1|sch   6.5    5.7    5.6   kill      5.4      60
>          yrs|sch   4.8    4.0      f   kill      0.5
>            1|sch   3.4    2.2    0.4    0.7      0.3       4
> yrs|id             5.7    5.7    8.2   kill     15       150
>   1|id             4.1    3.4    2.4   kill      4.2      14
>
> In the table above, "Setup" refers to the following two computer
> configurations:
>
> T30:  IBM Thinkpad T30, 1 GB ram, 1.8 Ghz processor, Windows 2000
>       SAMM version 1.1, lme4 & nlme current on 5.24.2005
>
> D620: Dell Latitude D620, 2 GB ram, 1.8 Ghz duo core processor, WinXP
>       SAMM version 1.1 lme4 & nlme current as of 1.30.2007
>
> I timed most model fits only once.  I did a quick inspection of the results
> from the different modelling functions to persuade myself that I was fitting
> the same model (i.e. that estimates were similar).

Thanks for sending the timings, Kevin.  As Harold Doran mentioned in
his reply to you, it is best to use the control options niterEM = 0
and gradient = FALSE when fitting lmer models to large data sets. The
theory of how to calculate the gradient and how to create the ECME
step for a linear mixed model is some of the most beautiful
mathematics that I have ever done.  Unfortunately it is not very
practical mathematics because one of those calculations is the
equivalent of several hundred function evaluations on models with
large data sets.  It is much better to just let the optimizer work out
an approximation to the gradient and to skip the EM or ECME iterations
altogether.

I repeated your timings on R-forge.R-project.org, an Opteron-based
server running Debian Linux.  I enclose the results.  Once you get rid
of the ECME iterations and the gradient calculations lmer and lmer2
are in the same ballpark with a slight advantage to lmer on most of
the models.  Both functions took about a minute or less to fit models
with up to 20000 random effects to about 25000 observations, which is
not bad at all.

As you know, not all these models can be fit by lme and, as the
timings show, lme is much slower on the models that it can fit.  This
is not surprising, there has been close to a decade's worth of
development of the computational methods since lme was initially
designed.



> Here are the full models I used:
>
> lmer(math~ gr + sx + eth + cltype + (1+yrs|id) + (1+yrs|sch), data=star)
>
> samm(math ~ gr + sx + eth + cltype, random=~ us(link(~1+yrs)):id +
>      us(link(~1+yrs)):sch, data=star, na.method.X="omit")
>
> Starting with these full models, I tried reduced models with simpler random
> effects structures.  These are the different rows in the table above.
>
> Observations:
> (1) The performance of the common version of SAMM on the two computers suggests
> the Dell is slightly faster than the IBM.
>
> (2) On 24 May 2005, lmer and samm has roughly similar timings.  On 31 Jan
> 2007, lmer is nearly unusable for this data (I killed the job after 5-10
> minutes of waiting).
>
> (3) The current version of lmer2 is the only function that appears to fit
> all models.
>
> (4) The current version of lme is slower than sammm/lmer2 for
>     those models I tried to fit.
>
> FYI.  SAMM is officially available at
> http://www.vsni.co.uk/products/samm/(but seems not to be there now,
> perhaps in preparation for release of a new
> version). Unofficially it is available here: ftp://ftp.dpi.nsw.gov.au/

I'm not sure if I have access to SAMM to check the timings.  Is SAMM
proprietary?  My recollection is that it is proprietary and that it is
Windows-only.  Either one of those characteristics is the kiss of
death for my wanting to use it.

Are you sure that SAMM is fitting a model with partially-crossed
random effects?  If you are only considering students and schools as
grouping factors for the random effects there will be little
difference between assuming nested random effects and correcting for
the partial crossing.  This is because there are very few students in
this study who attend more than one school.  If you considered
teachers as well as students and schools as grouping factors then you
could more easily discern the difference in model fits taking into
account the crossing or assuming that each student:teacher:school
combination is unique.

>
> I hope this information is useful.  Thanks for the progress evident in the
> lme4 package.
>
> Kevin Wright
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
-------------- next part --------------
Random.effects        lmer   lmer2   lme
(yrs|id)+(yrs|sch)   51.203 62.640   NA
(yrs|id)+(1|sch)     14.497 22.533   NA
(1|id)+(yrs|sch)     13.449  5.188   NA
(1|id)+(1|sch)        2.540  2.056   NA
(yrs|id)              6.409  9.992 84.037
(1|id)                0.628  1.020  5.637
(yrs|sch)             0.860  0.533  3.544
(1|sch)               0.544  0.472  2.404
-------------- next part --------------
library(lme4)
data(star, package = 'mlmRev')
system.time(lmer(math ~ gr + sx + eth + cltype + (yrs|id) + (yrs|sch),
                 star, control = list(niterEM = 0, gradient = FALSE)))
system.time(lmer(math ~ gr + sx + eth + cltype + (yrs|id) + (1|sch),
                 star, control = list(niterEM = 0, gradient = FALSE)))
system.time(lmer(math ~ gr + sx + eth + cltype + (1|id) + (yrs|sch),
                 star, control = list(niterEM = 0, gradient = FALSE)))
system.time(lmer(math ~ gr + sx + eth + cltype + (1|id) + (1|sch),
                 star, control = list(niterEM = 0, gradient = FALSE)))
system.time(lmer(math ~ gr + sx + eth + cltype + (yrs|sch),
                 star, control = list(niterEM = 0, gradient = FALSE)))
system.time(lmer(math ~ gr + sx + eth + cltype + (1|sch),
                 star, control = list(niterEM = 0, gradient = FALSE)))
system.time(lmer(math ~ gr + sx + eth + cltype + (yrs|id),
                 star, control = list(niterEM = 0, gradient = FALSE)))
system.time(lmer(math ~ gr + sx + eth + cltype + (1|id),
                 star, control = list(niterEM = 0, gradient = FALSE)))
system.time(lmer2(math ~ gr + sx + eth + cltype + (yrs|id) + (yrs|sch),
                 star, control = list(niterEM = 0, gradient = FALSE)))
system.time(lmer2(math ~ gr + sx + eth + cltype + (yrs|id) + (1|sch),
                 star, control = list(niterEM = 0, gradient = FALSE)))
system.time(lmer2(math ~ gr + sx + eth + cltype + (1|id) + (yrs|sch),
                 star, control = list(niterEM = 0, gradient = FALSE)))
system.time(lmer2(math ~ gr + sx + eth + cltype + (1|id) + (1|sch),
                 star, control = list(niterEM = 0, gradient = FALSE)))
system.time(lmer2(math ~ gr + sx + eth + cltype + (yrs|sch),
                 star, control = list(niterEM = 0, gradient = FALSE)))
system.time(lmer2(math ~ gr + sx + eth + cltype + (1|sch),
                 star, control = list(niterEM = 0, gradient = FALSE)))
system.time(lmer2(math ~ gr + sx + eth + cltype + (yrs|id),
                 star, control = list(niterEM = 0, gradient = FALSE)))
system.time(lmer2(math ~ gr + sx + eth + cltype + (1|id),
                 star, control = list(niterEM = 0, gradient = FALSE)))
-------------- next part --------------
library(nlme)
data(star, package = 'mlmRev')
system.time(lme(math ~ gr + sx + eth + cltype, star, random = ~yrs|sch, na.action = na.omit))
system.time(lme(math ~ gr + sx + eth + cltype, star, random = ~1|sch, na.action = na.omit))
system.time(lme(math ~ gr + sx + eth + cltype, star, random = ~yrs|id, na.action = na.omit))
system.time(lme(math ~ gr + sx + eth + cltype, star, random = ~1|id, na.action = na.omit))
-------------- next part --------------

R version 2.5.0 Under development (unstable) (2007-01-31 r40628)
Copyright (C) 2007 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

[Previously saved workspace restored]

> library(nlme)
> data(star, package = 'mlmRev')
> system.time(lme(math ~ gr + sx + eth + cltype, star, random = ~yrs|sch, na.action = na.omit))
   user  system elapsed 
  3.544   0.168   3.713 
> system.time(lme(math ~ gr + sx + eth + cltype, star, random = ~1|sch, na.action = na.omit))
   user  system elapsed 
  2.404   0.128   2.532 
> system.time(lme(math ~ gr + sx + eth + cltype, star, random = ~yrs|id, na.action = na.omit))
   user  system elapsed 
 84.037   1.588  85.672 
> system.time(lme(math ~ gr + sx + eth + cltype, star, random = ~1|id, na.action = na.omit))
   user  system elapsed 
  5.637   0.124   5.759 
> 
> proc.time()
   user  system elapsed 
 97.006   2.064  99.109 
-------------- next part --------------

R version 2.5.0 Under development (unstable) (2007-01-31 r40628)
Copyright (C) 2007 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

[Previously saved workspace restored]

> library(lme4)
Loading required package: Matrix
Loading required package: lattice
> data(star, package = 'mlmRev')
> system.time(lmer(math ~ gr + sx + eth + cltype + (yrs|id) + (yrs|sch),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
 51.203   6.944  58.218 
> system.time(lmer(math ~ gr + sx + eth + cltype + (yrs|id) + (1|sch),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
 14.497   0.172  14.669 
> system.time(lmer(math ~ gr + sx + eth + cltype + (1|id) + (yrs|sch),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
 13.449   0.140  13.624 
> system.time(lmer(math ~ gr + sx + eth + cltype + (1|id) + (1|sch),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
  2.540   0.080   2.623 
> system.time(lmer(math ~ gr + sx + eth + cltype + (yrs|sch),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
  0.860   0.064   0.925 
> system.time(lmer(math ~ gr + sx + eth + cltype + (1|sch),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
  0.544   0.060   0.606 
> system.time(lmer(math ~ gr + sx + eth + cltype + (yrs|id),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
  6.409   0.392   6.802 
> system.time(lmer(math ~ gr + sx + eth + cltype + (1|id),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
  0.628   0.084   0.711 
> system.time(lmer2(math ~ gr + sx + eth + cltype + (yrs|id) + (yrs|sch),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
 62.640   4.204  66.843 
> system.time(lmer2(math ~ gr + sx + eth + cltype + (yrs|id) + (1|sch),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
 22.533   0.944  23.492 
> system.time(lmer2(math ~ gr + sx + eth + cltype + (1|id) + (yrs|sch),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
  5.188   0.216   5.406 
> system.time(lmer2(math ~ gr + sx + eth + cltype + (1|id) + (1|sch),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
  2.056   0.032   2.089 
> system.time(lmer2(math ~ gr + sx + eth + cltype + (yrs|sch),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
  0.533   0.032   0.565 
> system.time(lmer2(math ~ gr + sx + eth + cltype + (1|sch),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
  0.472   0.016   0.486 
> system.time(lmer2(math ~ gr + sx + eth + cltype + (yrs|id),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
  9.992   0.348  10.364 
> system.time(lmer2(math ~ gr + sx + eth + cltype + (1|id),
+                  star, control = list(niterEM = 0, gradient = FALSE)))
   user  system elapsed 
  1.020   0.024   1.046 
> 
> proc.time()
   user  system elapsed 
204.576  13.948 218.748 


More information about the R-sig-mixed-models mailing list