[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