A good share of the confusion with this, and the first answer to it, is that different models
have been fitted with lmer (or aov) than the one used for the expected mean squares. So it is
like comparing apples and oranges. The model needed for comparison has random intercepts
for Fac2 and Fac1:Fac2:
> ward.lmer = lmer(Y ~ Fac1 + (1|Fac2) + (1|Fac1:Fac2), data = ward)
Here is the summary:
> ward.lmer
Linear mixed model fit by REML ['lmerMod']
Formula: Y ~ Fac1 + (1 | Fac2) + (1 | Fac1:Fac2)
   Data: ward
REML criterion at convergence: 67.0954
Random effects:
 Groups    Name        Std.Dev.
 Fac1:Fac2 (Intercept) 0.000
 Fac2      (Intercept) 7.681
 Residual              1.584
Number of obs: 18, groups: Fac1:Fac2, 6; Fac2, 2
Fixed Effects:
(Intercept)        Fac1B        Fac1C
  4.179e-14    2.667e+00   -1.333e+00
..  and the ANOVA:
> anova(ward.lmer)
Analysis of Variance Table
     Df Sum Sq Mean Sq F value
Fac1  2 49.778  24.889  9.9241
We see from this that the denominator of the F ratio is 24.89 / 9.9241 = 2.51.
That is not the same as the F ratio from the EMS approach BUT there is a reason:
The variance component estimate for Fac1:Fac2 is zero (because its MS is less
than the MSE). This is tantamount to deleting the interaction from the model.
Observe that if we pool interaction into error, we get
MSE = (1.78 + 33.33) / (2 + 12) = 2.51, i.e. the denominator actually used for F.

To see how this would work out if we don't have a zero variance estimate, I'll
create a new dataset by shifting the data in the first cell so that the interaction is greater:
> ward2 = ward
> ward2$Y[1:3] = ward2$Y[1:3] - 2
> ward2.lm = lm(Y ~ Fac1*Fac2, data = ward2)
> anova(ward2.lm)
Analysis of Variance Table

Response: Y
          Df Sum Sq Mean Sq F value    Pr(>F)
Fac1       2  59.11   29.56   10.64  0.002198
Fac2       1 470.22  470.22  169.28 1.954e-08
Fac1:Fac2  2  11.11    5.56    2.00  0.177979
Residuals 12  33.33    2.78
Thus, using the EMS approach, the correct F ratio for Fac1 is
F = 10.64 / 2.00 = 5.32. The MS for interaction is now larger than MSE.
That will make all the difference.

Here is the model fitted using lmer:
> ward2.lmer = lmer(Y ~ Fac1 + (1|Fac2) + (1|Fac1:Fac2), data = ward2)
> ward2.lmer
Linear mixed model fit by REML ['lmerMod']
Formula: Y ~ Fac1 + (1 | Fac2) + (1 | Fac1:Fac2)
   Data: ward2
REML criterion at convergence: 69.7861
Random effects:
 Groups    Name        Std.Dev.
 Fac1:Fac2 (Intercept) 0.9622
 Fac2      (Intercept) 7.1854
 Residual              1.6667
Number of obs: 18, groups: Fac1:Fac2, 6; Fac2, 2
Fixed Effects:
(Intercept)        Fac1B        Fac1C
    -1.0000       3.6667      -0.3333

> anova(ward2.lmer)
Analysis of Variance Table
     Df Sum Sq Mean Sq F value
Fac1  2 29.556  14.778    5.32
So now we obtain exactly the F ratio as the EMS method. In general, IF NO VARIANCES ARE
ESTIMATED TO BE ZERO in a balanced experiment, the F ratios for fixed effects will be identical
to those based on expected mean squares for the "unrestricted model." The reason for this is
that in this case, the REML estimates of the variance components are the same as the
method-of-moments estimates obtained from the expected mean squares.

Hope this helps.

Russell V. Lenth  -  Professor Emeritus
Department of Statistics and Actuarial Science
The University of Iowa  -  Iowa City, IA 52242  USA
Voice (319)335-0712 (Dept. office)  -  FAX (319)335-3017


	[[alternative HTML version deleted]]

