Skip to content

Some timings for lmer2 versus lmer

4 messages · Andrew Robinson, Douglas Bates

Original
Douglas Bates
# 
I enclose an R source file to do some comparative timings on lmer2
fits versus lmer fits and the output generated on the machine
R-forge.R-project.org (Opteron 280 dual-core processors, R internal
BLAS, 13 GB of memory).  You can try running the script on your
computer to get an idea of the timings.

On some machines the lmer fit to the "star" data set will converge in
considerably fewer iterations than on this machine.  There is one
point in the optimization where very small differences in the floating
point operation orders cause a much better step to be taken.

If you look at the verbose output you will see that the
parameterization for lmer2 uses the relative standard deviation (as
described in the Implementation vignette from the lme4 package)
whereas lmer used the relative variance.  Generally the relative
standard deviation is more stable for the optimization.

The other big difference in the optimization, shown in the last
example, is that lmer evaluates the relative precision matrix (the
inverse of the relative variance matrix) and therefore cannot allow
variance components to go to zero.  The value of a variance component
is bounded below at 5e-10, which is why that particular number shows
up in the verbose iterations.  As described in the vignette, the
relative variance matrix is used in lmer2 hence the lower bound on the
variance component is at zero.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: lmer2_test.R
Type: application/octet-stream
Size: 2967 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20070127/1ad3cb3d/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: lmer2_test.Rout
Type: application/octet-stream
Size: 45018 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20070127/1ad3cb3d/attachment-0001.obj>
Andrew Robinson
# 
I've switched from FreeBSD to WinXP temporarily :)

I've attached a comparison of lmer and lmer2 upon the analysis of
Cochran and Cox's chocolate cake data.  Here, it seems that lmer2 is
faster (0.08 vs. 0.15) but the AIC of the fitted model for lmer2 is
higher (1643 vs 1635).  The models are quite different in the random
effects.

The data, script, and output are attached.

Cheers,

Andrew
Douglas Bates
#
On 1/28/07, Andrew Robinson <A.Robinson at ms.unimelb.edu.au> wrote:
Thanks for including that example Andrew.  If you turn on the
msVerbose control setting you will see that it is a problem in the
optimizer behavior near the boundary for the new parameterization
(script and output attached).

It is a good thing to have such an example.  I had only observed the
opposite behavior where the optimizer had boundary problems in the
relative variance scale but not in the relative standard deviation
scale.  I'm not quite sure what I am going to do about it though.
-------------- next part --------------
# The chocolate cake breakage data are referred to in section 5.5 of
# \citet{lee+nelder+pawitan-2006} as an example of a normal-normal
# hiearchical generalized linear model. The data are originally from
# \citet{cochran+cox-1957}.

require(lme4)
load("~/tmp/cake.rda")
str(cake)
head(cake)

system.time(cake.lmer <- lmer(angle ~ recipe * temperature + 
                  (1 | replicate/batch),
                  data = cake,
                  control = list(msV = 1, niterEM = 0, grad = FALSE)))
system.time(cake.lmer2 <- lmer2(angle ~ recipe * temperature + 
                  (1 | replicate/batch),
                  data = cake, control = list(msV = 1)))
c(0.666667, 0.384900)^2
print(cake.lmer, corr = FALSE)
print(cake.lmer2, corr = FALSE)

sessionInfo()
-------------- next part --------------

R version 2.5.0 Under development (unstable) (2007-01-28 r40596)
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.
Loading required package: lme4
Loading required package: Matrix
Loading required package: lattice
[1] TRUE
'data.frame':	270 obs. of  5 variables:
 $ replicate  : Factor w/ 15 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ batch      : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 2 2 2 2 ...
 $ recipe     : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 2 2 2 2 ...
 $ temperature: Factor w/ 6 levels "175","185","195",..: 1 2 3 4 5 6 1 2 3 4 ...
 $ angle      : int  42 46 47 39 53 42 39 46 51 49 ...
replicate batch recipe temperature angle
1         1     1      1         175    42
2         1     1      1         185    46
3         1     1      1         195    47
4         1     1      1         205    39
5         1     1      1         215    53
6         1     1      1         225    42
+                   (1 | replicate/batch),
+                   data = cake,
+                   control = list(msV = 1, niterEM = 0, grad = FALSE)))
  0      1634.73: 0.444444 0.148148
  1      1603.11: 0.608821  1.13455
  2      1597.79: 0.320190  1.21276
  3      1596.71: 0.240022  1.23799
  4      1596.28: 0.164550  1.27497
  5      1596.15: 0.173446  1.30522
  6      1595.82: 0.199696  1.42855
  7      1595.57: 0.196701  1.55461
  8      1595.37: 0.182251  1.75435
  9      1595.35: 0.180630  1.83013
 10      1595.35: 0.181575  1.85835
 11      1595.35: 0.181800  1.86184
   user  system elapsed 
  0.084   0.000   0.085
+                   (1 | replicate/batch),
+                   data = cake, control = list(msV = 1)))
  0      1634.73: 0.666667 0.384900
  1      1606.88: 0.943988  1.34568
  2      1603.44:  0.00000  1.32963
  3      1603.43: 2.67531e-05  1.30813
  4      1603.43: 0.000311112  1.30263
   user  system elapsed 
  0.036   0.004   0.040
[1] 0.4444449 0.1481480
Linear mixed-effects model fit by REML 
Formula: angle ~ recipe * temperature + (1 | replicate/batch) 
   Data: cake 
  AIC  BIC logLik MLdeviance REMLdeviance
 1635 1707 -797.7       1638         1595
Random effects:
 Groups          Name        Variance Std.Dev.
 batch:replicate (Intercept)  3.7216  1.9292  
 replicate       (Intercept) 38.1137  6.1736  
 Residual                    20.4710  4.5245  
number of obs: 270, groups: batch:replicate, 45; replicate, 15

Fixed effects:
                       Estimate Std. Error t value
(Intercept)             29.1333     2.0381  14.295
recipe2                 -2.2667     1.7960  -1.262
recipe3                 -1.2000     1.7960  -0.668
temperature185           2.4000     1.6521   1.453
temperature195           1.6667     1.6521   1.009
temperature205           4.4000     1.6521   2.663
temperature215           9.5333     1.6521   5.770
temperature225           5.9333     1.6521   3.591
recipe2:temperature185   0.1333     2.3364   0.057
recipe3:temperature185  -1.4000     2.3364  -0.599
recipe2:temperature195   3.2000     2.3364   1.370
recipe3:temperature195   2.1333     2.3364   0.913
recipe2:temperature205   0.8667     2.3364   0.371
recipe3:temperature205  -1.4667     2.3364  -0.628
recipe2:temperature215  -1.9333     2.3364  -0.827
recipe3:temperature215  -3.0667     2.3364  -1.313
recipe2:temperature225   2.4667     2.3364   1.056
recipe3:temperature225   1.8667     2.3364   0.799
Linear mixed-effects model fit by REML 
Formula: angle ~ recipe * temperature + (1 | replicate/batch) 
   Data: cake 
  AIC  BIC logLik MLdeviance REMLdeviance
 1643 1715 -801.7       1647         1603
Random effects:
 Groups          Name        Variance   Std.Dev. 
 batch:replicate (Intercept) 2.2357e-06 0.0014952
 replicate       (Intercept) 3.9195e+01 6.2605706
 Residual                    2.3099e+01 4.8061039
Number of obs: 270, groups: batch:replicate, 45; replicate, 15

Fixed effects:
                       Estimate Std. Error t value
(Intercept)             29.1333     2.0379  14.296
recipe2                 -2.2667     1.7549  -1.292
recipe3                 -1.2000     1.7549  -0.684
temperature185           2.4000     1.7549   1.368
temperature195           1.6667     1.7549   0.950
temperature205           4.4000     1.7549   2.507
temperature215           9.5333     1.7549   5.432
temperature225           5.9333     1.7549   3.381
recipe2:temperature185   0.1333     2.4819   0.054
recipe3:temperature185  -1.4000     2.4819  -0.564
recipe2:temperature195   3.2000     2.4819   1.289
recipe3:temperature195   2.1333     2.4819   0.860
recipe2:temperature205   0.8667     2.4819   0.349
recipe3:temperature205  -1.4667     2.4819  -0.591
recipe2:temperature215  -1.9333     2.4819  -0.779
recipe3:temperature215  -3.0667     2.4819  -1.236
recipe2:temperature225   2.4667     2.4819   0.994
recipe3:temperature225   1.8667     2.4819   0.752
R version 2.5.0 Under development (unstable) (2007-01-28 r40596) 
x86_64-unknown-linux-gnu 

locale:
LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=en_US.UTF-8;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C

attached base packages:
[1] "stats"     "graphics"  "grDevices" "utils"     "datasets"  "methods"  
[7] "base"     

other attached packages:
       lme4      Matrix     lattice 
"0.9975-10"  "0.9975-9"   "0.14-16"
user  system elapsed 
  7.524   0.180   7.749
Douglas Bates
#
On 1/28/07, Douglas Bates <bates at stat.wisc.edu> wrote:
I have just committed a couple of changes to the SVN archive for the
lme4 package

https://svn.r-project.org/R-packages/trunk/lme4

that allow lmer2 to fit this model to these data and obtain the same
estimates as lmer did.

One approach is to fit a simpler model with additive fixed effects
first and use the fitted variance components from that model as the
starting estimates for the model that allows for interaction of the
fixed effects.  The version of lmer2 on the SVN archive allows you to
specify a start argument that should be in the form of the ST slot for
the model.  If you fit two models with the same random effects
specification then you can use the ST slot from the first as the
starting estimate for the second.

The other thing that I did was to change the call in lmer2 to the
nlminb optimizer so that it uses the default setting of the rel.tol
convergence criterion.  In the currently released version of lme4 the
this criterion is reset so that convergence is declared if the
deviance apparently has been determined to an accuracy of 0.001.  I
made this change because we observed that in many cases a substantial
portion of the iterations of the optimizer were spent at the optimum
making very small changes in parameter values that did not have a
substantial impact on the value of the deviance.

It looks like changing that criterion was too risky.  I would rather
have slower optimization with a higher degree of confidence that the
declared optimum is indeed the optimum.

I enclose R code and output for fitting these models with the modified
(and, as yet, unreleased) version of lmer2.  I also modified the cake
data so that the temperature is an ordered factor.  This results in
slightly different values of the REML criterion but you can still see
the pattern of convergence.  My purpose in using an ordered factor is
to see if the linear contrast in the temperature is dominant, which it
is.
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: cake_R.txt
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20070129/cce70675/attachment.txt>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: cake_Rout.txt
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20070129/cce70675/attachment-0001.txt>