[R-sig-ME] Timing for lmer2 versus lmer for chocolate cake data (WinXP)

Douglas Bates bates at stat.wisc.edu
Sun Jan 28 20:49:51 CET 2007


On 1/28/07, Andrew Robinson <A.Robinson at ms.unimelb.edu.au> wrote:
> 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.

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.

> # 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)
Loading required package: lme4
Loading required package: Matrix
Loading required package: lattice
[1] TRUE
> load("~/tmp/cake.rda")
> str(cake)
'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 ...
> head(cake)
  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
> 
> system.time(cake.lmer <- lmer(angle ~ recipe * temperature + 
+                   (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 
> system.time(cake.lmer2 <- lmer2(angle ~ recipe * temperature + 
+                   (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 
> c(0.666667, 0.384900)^2
[1] 0.4444449 0.1481480
> print(cake.lmer, corr = FALSE)
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
> print(cake.lmer2, corr = FALSE)
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
> 
> sessionInfo()
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" 
> 
> 
> proc.time()
   user  system elapsed 
  7.524   0.180   7.749 
> 


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