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("./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: Ord.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), cake, + control = list(msV = 1, niterEM = 0, + grad = FALSE))) 0 1640.11: 0.444444 0.148148 1 1608.49: 0.608820 1.13455 2 1603.17: 0.320190 1.21276 3 1602.09: 0.240022 1.23799 4 1601.65: 0.164551 1.27497 5 1601.52: 0.173446 1.30522 6 1601.19: 0.199692 1.42855 7 1600.95: 0.196699 1.55461 8 1600.74: 0.182252 1.75435 9 1600.72: 0.180630 1.83013 10 1600.72: 0.181575 1.85835 11 1600.72: 0.181800 1.86184 user system elapsed 0.084 0.004 0.087 > system.time(cake1.lmer2 <- lmer2(angle ~ recipe + temperature + + (1 | replicate/batch), cake, + control = list(msV = 1))) 0 1677.30: 0.666667 0.384900 1 1649.45: 0.947530 1.34465 2 1645.93: 0.00000 1.34124 3 1645.91: 2.61328e-05 1.31434 4 1645.91: 0.000324928 1.30529 5 1645.91: 0.00299488 1.30485 6 1645.78: 0.0298808 1.30575 7 1637.90: 0.460297 1.30683 8 1637.79: 0.428034 1.31761 9 1637.77: 0.424685 1.34472 10 1637.77: 0.426129 1.36283 11 1637.77: 0.426281 1.36429 12 1637.77: 0.426282 1.36433 user system elapsed 0.044 0.000 0.043 > system.time(cake2.lmer2 <- lmer2(angle ~ recipe * temperature + + (1 | replicate/batch), cake, + control = list(msV = 1))) 0 1640.11: 0.666667 0.384900 1 1612.25: 0.943989 1.34568 2 1608.82: 0.00000 1.32963 3 1608.81: 2.68979e-05 1.30813 4 1608.81: 0.000312802 1.30263 5 1608.81: 0.00298422 1.30136 6 1608.72: 0.0244528 1.30249 7 1601.02: 0.368425 1.30372 8 1600.76: 0.430777 1.31440 9 1600.75: 0.405951 1.37259 10 1600.73: 0.437475 1.37000 11 1600.73: 0.426247 1.34043 12 1600.72: 0.424209 1.37199 13 1600.72: 0.426975 1.37046 14 1600.72: 0.425629 1.36760 15 1600.72: 0.427518 1.36506 16 1600.72: 0.426289 1.36406 17 1600.72: 0.426391 1.36457 18 1600.72: 0.426398 1.36452 user system elapsed 0.044 0.000 0.046 > system.time(cake2.lmer2 <- lmer2(angle ~ recipe * temperature + + (1 | replicate/batch), cake, + start = cake1.lmer2@ST, + control = list(msV = 1))) 0 1600.72: 0.426282 1.36433 1 1600.72: 0.426437 1.36438 2 1600.72: 0.426372 1.36453 3 1600.72: 0.426404 1.36452 4 1600.72: 0.426404 1.36452 user system elapsed 0.040 0.000 0.036 > 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 1641 1713 -800.4 1638 1601 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) 33.12222 1.73680 19.071 recipe2 -1.47778 0.97526 -1.515 recipe3 -1.52222 0.97526 -1.561 temperature.L 6.43033 1.16822 5.504 temperature.Q -0.71285 1.16822 -0.610 temperature.C -2.32551 1.16822 -1.991 temperature^4 -3.35128 1.16822 -2.869 temperature^5 -0.15119 1.16822 -0.129 recipe2:temperature.L 0.45419 1.65211 0.275 recipe3:temperature.L 0.08765 1.65211 0.053 recipe2:temperature.Q -0.23277 1.65211 -0.141 recipe3:temperature.Q 1.21475 1.65211 0.735 recipe2:temperature.C 2.69322 1.65211 1.630 recipe3:temperature.C 2.63856 1.65211 1.597 recipe2:temperature^4 3.02372 1.65211 1.830 recipe3:temperature^4 3.13711 1.65211 1.899 recipe2:temperature^5 -0.66354 1.65211 -0.402 recipe3:temperature^5 -1.62525 1.65211 -0.984 > print(cake2.lmer2, corr = FALSE) Linear mixed-effects model fit by REML Formula: angle ~ recipe * temperature + (1 | replicate/batch) Data: cake AIC BIC logLik MLdeviance REMLdeviance 1641 1713 -800.4 1638 1601 Random effects: Groups Name Variance Std.Dev. batch:replicate (Intercept) 3.722 1.9293 replicate (Intercept) 38.115 6.1738 Residual 20.471 4.5245 Number of obs: 270, groups: batch:replicate, 45; replicate, 15 Fixed effects: Estimate Std. Error t value (Intercept) 33.12222 1.73684 19.070 recipe2 -1.47778 0.97528 -1.515 recipe3 -1.52222 0.97528 -1.561 temperature.L 6.43033 1.16821 5.504 temperature.Q -0.71285 1.16821 -0.610 temperature.C -2.32551 1.16821 -1.991 temperature^4 -3.35128 1.16821 -2.869 temperature^5 -0.15119 1.16821 -0.129 recipe2:temperature.L 0.45419 1.65210 0.275 recipe3:temperature.L 0.08765 1.65210 0.053 recipe2:temperature.Q -0.23277 1.65210 -0.141 recipe3:temperature.Q 1.21475 1.65210 0.735 recipe2:temperature.C 2.69322 1.65210 1.630 recipe3:temperature.C 2.63856 1.65210 1.597 recipe2:temperature^4 3.02372 1.65210 1.830 recipe3:temperature^4 3.13711 1.65210 1.899 recipe2:temperature^5 -0.66354 1.65210 -0.402 recipe3:temperature^5 -1.62525 1.65210 -0.984 > anova(cake1.lmer2, cake2.lmer2) Data: cake Models: cake1.lmer2: angle ~ recipe + temperature + (1 | replicate/batch) cake2.lmer2: angle ~ recipe * temperature + (1 | replicate/batch) Df AIC BIC logLik Chisq Chi Df Pr(>Chisq) cake1.lmer2.p 10 1669.02 1705.00 -824.51 cake2.lmer2.p 20 1678.45 1750.41 -819.22 10.570 10 0.3919 > > 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-11" "0.9975-9" "0.14-16" > > > proc.time() user system elapsed 8.504 0.156 8.693 >