I am trying to replicate the results of an aov command with lmer, to understand
the syntax, but I can't quite figure it out. I have a dataset from
Montgomery p. 520 with a nested and factorial layout. There are 3 fixtures, 2
layouts (the treatments) in a factorial design, but the operators who perform
the runs are nest in layouts (4 operators per layout=8 different operators).
Thus, each operator tries each fixture twice, but only 1 layout. Everything is
balanced so I know what the results should be. I give a dump of the data at the
bottom of the email.
If I use the aov command, I get exactly what I would
expect:> ass.aov<-aov(Time~Fixture*Layout+Error(Operator/Layout),data=ass)
> summary(ass.aov)
Error: Operator
Df Sum Sq Mean Sq F value Pr(>F)
Layout 1 4.083 4.0833 0.3407 0.5807
Residuals 6 71.917 11.9861
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Fixture 2 82.792 41.396 12.2319 8.842e-05 ***
Fixture:Layout 2 19.042 9.521 2.8133 0.07325 .
Residuals 36 121.833 3.384
If I use the lmer command, the layout SS is NOT what I would expect from the aov
(and my calculations). However, the estimates of variance components (not shown
here) are exactly the same. Note, I have coded Operator as 1-8 (not 1-4), so if
understand lmer correctly I can just add a (1|Operator) term into the model and
the nesting is taken care of.
If I try to add a Operator:Fixture interaction (as done in Montgomery) I also do
not get agreement with my manual calculations (which are the same as in the
Montgomery book). Can I recreate the standard anova table broken into the
correct strata using lmer?
Thanks,
Elizabeth
> ass.lmer2<-lmer(Time~Layout*Fixture+(1|Operator),data=ass)
> anova(ass.lmer2)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
Layout 1 1.153 1.153 0.3406
Fixture 2 82.792 41.396 12.2319
Layout:Fixture 2 19.042 9.521 2.8133> summary(ass.lmer2)
Linear mixed model fit by REML
Formula: Time ~ Layout + (1 | Operator) + Fixture + Layout:Fixture
Data: ass
AIC BIC logLik deviance REMLdev
215 230 -99.5 198.4 199
Random effects:
Groups Name Variance Std.Dev.
Operator (Intercept) 1.4336 1.1973
Residual 3.3843 1.8396
Number of obs: 48, groups: Operator, 8
Fixed effects:
Estimate Std. Error t value
(Intercept) 26.0833 0.4997 52.20
Layout1 -0.2917 0.4997 -0.58
Fixture1 -0.8333 0.3755 -2.22
Fixture2 1.8542 0.3755 4.94
Layout1:Fixture1 -0.2083 0.3755 -0.55
Layout1:Fixture2 0.8542 0.3755 2.27
Correlation of Fixed Effects:
(Intr) Layot1 Fixtr1 Fixtr2 Ly1:F1
Layout1 0.000
Fixture1 0.000 0.000
Fixture2 0.000 0.000 -0.500
Layt1:Fxtr1 0.000 0.000 0.000 0.000
Layt1:Fxtr2 0.000 0.000 0.000 0.000 -0.500
ass <-
structure(list(Time = c(22, 24, 30, 27, 25, 21, 23, 24, 29, 28,
24, 22, 28, 29, 30, 32, 27, 25, 25, 23, 27, 25, 26, 23, 26, 28,
29, 28, 27, 25, 27, 25, 30, 27, 26, 24, 28, 25, 24, 23, 24, 27,
24, 23, 28, 30, 28, 27), Operator = structure(c(1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L,
4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L,
7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L), .Label = c("1",
"2", "3", "4", "5", "6",
"7", "8"), class = "factor"), Layout =
structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label =
c("1",
"2"), class = "factor"), Fixture = structure(c(1L, 1L, 2L,
2L,
3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L,
2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L,
1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L), .Label = c("1",
"2", "3"), class = "factor"), opF =
structure(c(1L, 1L, 9L, 9L,
17L, 17L, 2L, 2L, 10L, 10L, 18L, 18L, 3L, 3L, 11L, 11L, 19L,
19L, 4L, 4L, 12L, 12L, 20L, 20L, 5L, 5L, 13L, 13L, 21L, 21L,
6L, 6L, 14L, 14L, 22L, 22L, 7L, 7L, 15L, 15L, 23L, 23L, 8L, 8L,
16L, 16L, 24L, 24L), .Label = c("1:1", "1:2",
"1:3", "1:4", "1:5",
"1:6", "1:7", "1:8", "2:1",
"2:2", "2:3", "2:4", "2:5",
"2:6",
"2:7", "2:8", "3:1", "3:2",
"3:3", "3:4", "3:5", "3:6",
"3:7",
"3:8"), class = "factor"), LF = structure(c(1L, 1L, 3L, 3L,
5L,
5L, 1L, 1L, 3L, 3L, 5L, 5L, 1L, 1L, 3L, 3L, 5L, 5L, 1L, 1L, 3L,
3L, 5L, 5L, 2L, 2L, 4L, 4L, 6L, 6L, 2L, 2L, 4L, 4L, 6L, 6L, 2L,
2L, 4L, 4L, 6L, 6L, 2L, 2L, 4L, 4L, 6L, 6L), .Label = c("1:1",
"1:2", "2:1", "2:2", "3:1",
"3:2"), class = "factor")), .Names = c("Time",
"Operator", "Layout", "Fixture", "opF",
"LF"), row.names = c(NA,
-48L), class = "data.frame")
