[R] replicating aov results with lmer

[Elizabeth Purdom]

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")