[R-sig-ME] Examples for speed comparisons in fitting linear mixed-effects models
Douglas Bates
bates at stat.wisc.edu
Thu Sep 24 23:56:59 CEST 2009
In some ways this is responding to the messages from Paul, Kevin,
Harold and others but it involves enough new material that I will
start a new thread.
I would very much appreciate examples of data sets and corresponding
models that can be used to assess the speed and reliability of fitting
linear mixed-effects models. As both Harold and I mentioned in the
other thread it is difficult to get a meaningful timing on a model fit
when it only takes a few seconds. I enclose output (star1_Rout.txt)
from fitting a model to the star data. You will see that I timed the
fitting several times with results between 10.3 seconds and 4.3
seconds. Some of the variability is due to garbage collection and the
amount of memory that is currently allocated but not in use. That is
why later timings tend to be faster. But there are many other sources
of variability in addition to memory allocation.
Besides the difficulty in reproducing fast timings, there is also the
"so what" factor. I really don't care whether the model is fit in 5
seconds or in 10 seconds.
But notice that this is a model with 3 partially crossed (i.e.
non-nested) levels of random effects for a total of about 11000 random
effects, 25000 observations and 12 fixed-effects parameters. That
used to be considered big and difficult but it isn't any more. The
second set of timings (star2_Rout.txt) are from an experimental
version of lme4 (called lme4a for the time being) and two functions,
lmer and lmer2. (The lmer2 name is going to change so don't pay too
much attention to it now.) The lmer2 function produces an environment
with certain extractor functions, one of which, setPars, sets new
values of the parameters and causes re-evaluation of the deviance or
the REML criterion. By default it also optimizes the criterion using
the nlminb optimizer. Kate Mullen and others have an R-forge project
on optimization including a package called minqa that provides an
optimizer bobyqa (love those Fortran, 6-character names) which
outperforms nlminb on my examples. One advantage of the lmer2
structure is that the optimization can be separated from the model
specification so you can switch optimizers. With the bobyqa optimizer
and the lmer2 formulation this model can be fit in under 2 seconds.
I do have timings on a much larger and more difficult model fit but
the data are confidential. This is a model with four, non-nested
levels of random effects (in this case student, school, district and
area) and 40 fixed-effects parameters in total. The overall
dimensions are
> print(fm1, corr = FALSE)
Linear mixed model fit by REML
Formula: Rss ~ Yr * fgr + Sx * Ra + Dis + ELP + EC + (1 | id) + (1 |
sch) + (1 | dist) + (1 | area)
Data: dat
REML
16917518
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 1357.394 36.8428
sch (Intercept) 64.341 8.0213
dist (Intercept) 22.408 4.7337
area (Intercept) 6.220 2.4940
Residual 498.701 22.3316
Number of obs: 1714674, groups: id, 735604; sch, 2409; dist, 445; area, 12
In this case it takes about a minute to set up the model structures
and each evaluation of the REML criterion takes about 8 seconds.
Again, bobyqa does much better than nlminb, about 1600 seconds
compared to 2500 seconds and both are much faster than the currently
released lme4.
I can't release these data. I don't even want to post estimates of
the fixed-effects parameters. I could, for the purposes of timing,
generate a dummy data set from model structure and parameter estimates
but I would prefer not to do that. Also, I want to be able to test
the new formulation of lmer on real models fit to real data by real
useRs.
So if you have data to which you want to fit linear mixed-effects
models and they seem to be taking a long time then please send me your
data and model specification. If the model could be fit by HLM6 or
MLWin or SAS PROC MIXED or Stata (I'm reasonably sure that even the
small example could not because of the partially crossed random
effects but I would be happy to hear otherwise) then we could compare
but right now I am mostly interested in the lme4a package methods and
the choice of optimizers bobyqa or nlminb.
-------------- next part --------------
R version 2.9.1 (2009-06-26)
Copyright (C) 2009 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.
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.
> library(lme4)
Loading required package: Matrix
Loading required package: lattice
Attaching package: 'Matrix'
The following object(s) are masked from package:stats :
contr.SAS,
contr.helmert,
contr.poly,
contr.sum,
contr.treatment,
xtabs
The following object(s) are masked from package:base :
rcond
> data(star, package = "mlmRev")
> system.time(fm1 <- lmer(math ~ gr + cltype + sx + eth + (1|id) + (1|tch) + (1|sch), star))
user system elapsed
10.303 0.050 10.824
> system.time(fm1 <- lmer(math ~ gr + cltype + sx + eth + (1|id) + (1|tch) + (1|sch), star))
user system elapsed
7.752 0.056 8.369
> system.time(fm1 <- lmer(math ~ gr + cltype + sx + eth + (1|id) + (1|tch) + (1|sch), star))
user system elapsed
4.311 0.005 4.777
> system.time(fm1 <- lmer(math ~ gr + cltype + sx + eth + (1|id) + (1|tch) + (1|sch), star))
user system elapsed
4.895 0.001 5.036
> print(fm1, corr = FALSE)
Linear mixed model fit by REML
Formula: math ~ gr + cltype + sx + eth + (1 | id) + (1 | tch) + (1 | sch)
Data: star
AIC BIC logLik deviance REMLdev
239236 239365 -119602 239245 239204
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 1001.59 31.648
tch (Intercept) 295.04 17.177
sch (Intercept) 104.76 10.235
Residual 397.32 19.933
Number of obs: 24578, groups: id, 10732; tch, 1374; sch, 80
Fixed effects:
Estimate Std. Error t value
(Intercept) 560.8869 1.5742 356.3
gr.L 95.6128 1.0080 94.9
gr.Q -4.6783 0.9858 -4.7
gr.C -3.2320 0.9729 -3.3
cltypereg -7.7943 1.3322 -5.9
cltypereg+A -7.0055 1.3265 -5.3
sxF 2.8766 0.6836 4.2
ethB -21.9199 1.2218 -17.9
ethA 3.1941 6.7036 0.5
ethH 4.5410 9.6650 0.5
ethI -28.1621 13.5356 -2.1
ethO 3.4358 8.3131 0.4
>
>
> proc.time()
user system elapsed
45.579 0.322 48.042
-------------- next part --------------
R version 2.9.1 (2009-06-26)
Copyright (C) 2009 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.
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.
> library(lme4a)
Loading required package: Matrix
Loading required package: lattice
Attaching package: 'Matrix'
The following object(s) are masked from package:stats :
contr.SAS,
contr.helmert,
contr.poly,
contr.sum,
contr.treatment,
xtabs
The following object(s) are masked from package:base :
rcond
> sessionInfo()
R version 2.9.1 (2009-06-26)
x86_64-unknown-linux-gnu
locale:
C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4a_0.999375-42 Matrix_0.999375-31 lattice_0.17-25
loaded via a namespace (and not attached):
[1] grid_2.9.1
> data(star, package = "mlmRev")
> system.time(fm1 <- lmer(math ~ gr + cltype + sx + eth + (1|id) + (1|tch) + (1|sch), star))
user system elapsed
14.076 0.070 14.960
Warning message:
In merFinalize(rho) : false convergence (8)
> system.time(fm1 <- lmer(math ~ gr + cltype + sx + eth + (1|id) + (1|tch) + (1|sch), star))
user system elapsed
7.897 0.018 8.297
Warning message:
In merFinalize(rho) : false convergence (8)
> system.time(fm1 <- lmer(math ~ gr + cltype + sx + eth + (1|id) + (1|tch) + (1|sch), star))
user system elapsed
7.852 0.001 8.443
Warning message:
In merFinalize(rho) : false convergence (8)
> system.time(fm1 <- lmer(math ~ gr + cltype + sx + eth + (1|id) + (1|tch) + (1|sch), star))
user system elapsed
8.619 0.007 9.208
Warning message:
In merFinalize(rho) : false convergence (8)
> system.time(fm1 <- lmer2(math ~ gr + cltype + sx + eth + (1|id) + (1|tch) + (1|sch), star))
user system elapsed
2.555 0.002 2.565
> system.time(fm1 <- lmer2(math ~ gr + cltype + sx + eth + (1|id) + (1|tch) + (1|sch), star))
user system elapsed
2.166 0.001 2.308
> system.time(fm1 <- lmer2(math ~ gr + cltype + sx + eth + (1|id) + (1|tch) + (1|sch), star))
user system elapsed
2.207 0.001 2.224
> system.time(fm1 <- lmer2(math ~ gr + cltype + sx + eth + (1|id) + (1|tch) + (1|sch), star))
user system elapsed
2.936 0.000 2.960
> print(fm1, corr = FALSE)
Linear mixed model fit by REML
Formula: math ~ gr + cltype + sx + eth + (1 | id) + (1 | tch) + (1 | sch)
Data: star
REML
239204
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 1001.60 31.648
tch (Intercept) 295.04 17.177
sch (Intercept) 104.76 10.235
Residual 397.31 19.933
Number of obs: 24578, groups: id, 10732; tch, 1374; sch, 80
Fixed effects:
Estimate Std. Error t value
(Intercept) 560.8868 1.5747 356.2
gr.L 95.6126 1.0080 94.9
gr.Q -4.6782 0.9858 -4.7
gr.C -3.2320 0.9729 -3.3
cltypereg -7.7945 1.3322 -5.9
cltypereg+A -7.0056 1.3265 -5.3
sxF 2.8766 0.6836 4.2
ethB -21.9196 1.2219 -17.9
ethA 3.1940 6.7036 0.5
ethH 4.5413 9.6650 0.5
ethI -28.1612 13.5356 -2.1
ethO 3.4357 8.3131 0.4
> ## time the setup phase only
> system.time(fm1 <- lmer2(math ~ gr + cltype + sx + eth + (1|id) + (1|tch) + (1|sch), star, doFit = FALSE))
user system elapsed
0.398 0.000 0.398
> ## time the optimization using nlminb
> system.time(nlminb(c(1,1,1), fm1 at setPars, lower = c(0,0,0), control = list(trace = 1)))
0: 240730.91: 1.00000 1.00000 1.00000
1: 239648.32: 1.98282 0.816539 0.979744
2: 239342.22: 1.75919 1.02726 0.914733
3: 239259.59: 1.49223 0.865024 0.882315
4: 239230.77: 1.65620 0.811277 0.619893
5: 239225.00: 1.58040 0.894458 0.326679
6: 239216.42: 1.55704 0.901629 0.639797
7: 239207.82: 1.59249 0.869466 0.630141
8: 239206.23: 1.58324 0.844715 0.589082
9: 239204.54: 1.59280 0.873056 0.550491
10: 239204.24: 1.59525 0.852259 0.506384
11: 239203.77: 1.58797 0.859716 0.507220
12: 239203.75: 1.58832 0.862046 0.517409
13: 239203.74: 1.58797 0.862118 0.513037
14: 239203.74: 1.58772 0.861717 0.513120
15: 239203.74: 1.58773 0.861747 0.513507
16: 239203.74: 1.58774 0.861740 0.513496
user system elapsed
2.400 0.001 2.582
> system.time(nlminb(c(1,1,1), fm1 at setPars, lower = c(0,0,0), control = list(trace = 1)))
0: 240730.91: 1.00000 1.00000 1.00000
1: 239648.32: 1.98282 0.816539 0.979744
2: 239342.22: 1.75919 1.02726 0.914733
3: 239259.59: 1.49223 0.865024 0.882315
4: 239230.77: 1.65620 0.811277 0.619893
5: 239225.00: 1.58040 0.894458 0.326679
6: 239216.42: 1.55704 0.901629 0.639797
7: 239207.82: 1.59249 0.869466 0.630141
8: 239206.23: 1.58324 0.844715 0.589082
9: 239204.54: 1.59280 0.873056 0.550491
10: 239204.24: 1.59525 0.852259 0.506384
11: 239203.77: 1.58797 0.859716 0.507220
12: 239203.75: 1.58832 0.862046 0.517409
13: 239203.74: 1.58797 0.862118 0.513037
14: 239203.74: 1.58772 0.861717 0.513120
15: 239203.74: 1.58773 0.861747 0.513507
16: 239203.74: 1.58774 0.861740 0.513496
user system elapsed
1.947 0.002 2.108
> system.time(nlminb(c(1,1,1), fm1 at setPars, lower = c(0,0,0), control = list(trace = 1)))
0: 240730.91: 1.00000 1.00000 1.00000
1: 239648.32: 1.98282 0.816539 0.979744
2: 239342.22: 1.75919 1.02726 0.914733
3: 239259.59: 1.49223 0.865024 0.882315
4: 239230.77: 1.65620 0.811277 0.619893
5: 239225.00: 1.58040 0.894458 0.326679
6: 239216.42: 1.55704 0.901629 0.639797
7: 239207.82: 1.59249 0.869466 0.630141
8: 239206.23: 1.58324 0.844715 0.589082
9: 239204.54: 1.59280 0.873056 0.550491
10: 239204.24: 1.59525 0.852259 0.506384
11: 239203.77: 1.58797 0.859716 0.507220
12: 239203.75: 1.58832 0.862046 0.517409
13: 239203.74: 1.58797 0.862118 0.513037
14: 239203.74: 1.58772 0.861717 0.513120
15: 239203.74: 1.58773 0.861747 0.513507
16: 239203.74: 1.58774 0.861740 0.513496
user system elapsed
2.015 0.000 2.016
> system.time(nlminb(c(1,1,1), fm1 at setPars, lower = c(0,0,0), control = list(trace = 1)))
0: 240730.91: 1.00000 1.00000 1.00000
1: 239648.32: 1.98282 0.816539 0.979744
2: 239342.22: 1.75919 1.02726 0.914733
3: 239259.59: 1.49223 0.865024 0.882315
4: 239230.77: 1.65620 0.811277 0.619893
5: 239225.00: 1.58040 0.894458 0.326679
6: 239216.42: 1.55704 0.901629 0.639797
7: 239207.82: 1.59249 0.869466 0.630141
8: 239206.23: 1.58324 0.844715 0.589082
9: 239204.54: 1.59280 0.873056 0.550491
10: 239204.24: 1.59525 0.852259 0.506384
11: 239203.77: 1.58797 0.859716 0.507220
12: 239203.75: 1.58832 0.862046 0.517409
13: 239203.74: 1.58797 0.862118 0.513037
14: 239203.74: 1.58772 0.861717 0.513120
15: 239203.74: 1.58773 0.861747 0.513507
16: 239203.74: 1.58774 0.861740 0.513496
user system elapsed
1.896 0.000 1.897
> ## time the optimization using bobyqa
> library(minqa)
> system.time(bobyqa(c(1,1,1), fm1 at setPars, lower = c(0,0,0), control = list(iprint = 2)))
user system elapsed
1.674 0.000 1.795
> system.time(bobyqa(c(1,1,1), fm1 at setPars, lower = c(0,0,0), control = list(iprint = 2)))
user system elapsed
1.577 0.000 1.578
> system.time(bobyqa(c(1,1,1), fm1 at setPars, lower = c(0,0,0), control = list(iprint = 2)))
user system elapsed
1.488 0.000 1.488
> system.time(bobyqa(c(1,1,1), fm1 at setPars, lower = c(0,0,0), control = list(iprint = 2)))
user system elapsed
1.490 0.000 1.602
>
>
>
> proc.time()
user system elapsed
82.897 0.339 87.485
New RHO = 1.0000D-01 Number of function values = 6
Least value of F = 2.396525830667515D+05 The corresponding X is:
2.000000D+00 1.000000D+00 1.000000D+00
New RHO = 1.0000D-02 Number of function values = 16
Least value of F = 2.392271194838589D+05 The corresponding X is:
1.530990D+00 9.271466D-01 5.263914D-01
New RHO = 1.0000D-03 Number of function values = 22
Least value of F = 2.392039829385394D+05 The corresponding X is:
1.586325D+00 8.555261D-01 5.344145D-01
New RHO = 1.0000D-04 Number of function values = 38
Least value of F = 2.392037396723067D+05 The corresponding X is:
1.587629D+00 8.618595D-01 5.139103D-01
New RHO = 1.0000D-05 Number of function values = 45
Least value of F = 2.392037395221720D+05 The corresponding X is:
1.587714D+00 8.617400D-01 5.134964D-01
New RHO = 1.0000D-06 Number of function values = 51
Least value of F = 2.392037395205490D+05 The corresponding X is:
1.587737D+00 8.617412D-01 5.135009D-01
At the return from BOBYQA Number of function values = 58
Least value of F = 2.392037395205379D+05 The corresponding X is:
1.587738D+00 8.617422D-01 5.134982D-01
New RHO = 1.0000D-01 Number of function values = 6
Least value of F = 2.396525830667515D+05 The corresponding X is:
2.000000D+00 1.000000D+00 1.000000D+00
New RHO = 1.0000D-02 Number of function values = 16
Least value of F = 2.392271194838589D+05 The corresponding X is:
1.530990D+00 9.271466D-01 5.263914D-01
New RHO = 1.0000D-03 Number of function values = 22
Least value of F = 2.392039829385394D+05 The corresponding X is:
1.586325D+00 8.555261D-01 5.344145D-01
New RHO = 1.0000D-04 Number of function values = 38
Least value of F = 2.392037396723067D+05 The corresponding X is:
1.587629D+00 8.618595D-01 5.139103D-01
New RHO = 1.0000D-05 Number of function values = 45
Least value of F = 2.392037395221720D+05 The corresponding X is:
1.587714D+00 8.617400D-01 5.134964D-01
New RHO = 1.0000D-06 Number of function values = 51
Least value of F = 2.392037395205490D+05 The corresponding X is:
1.587737D+00 8.617412D-01 5.135009D-01
At the return from BOBYQA Number of function values = 58
Least value of F = 2.392037395205379D+05 The corresponding X is:
1.587738D+00 8.617422D-01 5.134982D-01
New RHO = 1.0000D-01 Number of function values = 6
Least value of F = 2.396525830667515D+05 The corresponding X is:
2.000000D+00 1.000000D+00 1.000000D+00
New RHO = 1.0000D-02 Number of function values = 16
Least value of F = 2.392271194838589D+05 The corresponding X is:
1.530990D+00 9.271466D-01 5.263914D-01
New RHO = 1.0000D-03 Number of function values = 22
Least value of F = 2.392039829385394D+05 The corresponding X is:
1.586325D+00 8.555261D-01 5.344145D-01
New RHO = 1.0000D-04 Number of function values = 38
Least value of F = 2.392037396723067D+05 The corresponding X is:
1.587629D+00 8.618595D-01 5.139103D-01
New RHO = 1.0000D-05 Number of function values = 45
Least value of F = 2.392037395221720D+05 The corresponding X is:
1.587714D+00 8.617400D-01 5.134964D-01
New RHO = 1.0000D-06 Number of function values = 51
Least value of F = 2.392037395205490D+05 The corresponding X is:
1.587737D+00 8.617412D-01 5.135009D-01
At the return from BOBYQA Number of function values = 58
Least value of F = 2.392037395205379D+05 The corresponding X is:
1.587738D+00 8.617422D-01 5.134982D-01
New RHO = 1.0000D-01 Number of function values = 6
Least value of F = 2.396525830667515D+05 The corresponding X is:
2.000000D+00 1.000000D+00 1.000000D+00
New RHO = 1.0000D-02 Number of function values = 16
Least value of F = 2.392271194838589D+05 The corresponding X is:
1.530990D+00 9.271466D-01 5.263914D-01
New RHO = 1.0000D-03 Number of function values = 22
Least value of F = 2.392039829385394D+05 The corresponding X is:
1.586325D+00 8.555261D-01 5.344145D-01
New RHO = 1.0000D-04 Number of function values = 38
Least value of F = 2.392037396723067D+05 The corresponding X is:
1.587629D+00 8.618595D-01 5.139103D-01
New RHO = 1.0000D-05 Number of function values = 45
Least value of F = 2.392037395221720D+05 The corresponding X is:
1.587714D+00 8.617400D-01 5.134964D-01
New RHO = 1.0000D-06 Number of function values = 51
Least value of F = 2.392037395205490D+05 The corresponding X is:
1.587737D+00 8.617412D-01 5.135009D-01
At the return from BOBYQA Number of function values = 58
Least value of F = 2.392037395205379D+05 The corresponding X is:
1.587738D+00 8.617422D-01 5.134982D-01
More information about the R-sig-mixed-models
mailing list