Hopefully this is a question that someone here can answer for me on the nature of decomposing sums of squares from a mixed-effects model fit with
lmer (from the lme4 R package).
First off I should say that I am aware of the controversy with using this approach, and in practise I would be more likely to use a bootstrapped LRT to compare models (as suggested by Faraway, 2006). However, I am puzzled at how to replicate the results, and so for my own sanity I thought I would ask here.
Basically, I am getting to grips with using mixed-effects models fit by the
lme4 package. I know that you can use the
anova() command to give a summary of sequentially testing the fixed-effects in the model. As far as I know this is what Faraway (2006) refers to as the 'Expected mean squares' approach. What I want to know is how are the sums of squares calculated?
I know that I could take the estimated values from a particular model (using
coef()), assume that they are fixed, and then make tests using the sums of squares of model residuals with and without the factors of interest. This is fine for a model containing a single within-subject factor. However, when implementing a split-plot design the sums of squares value I get is equivalent to the value produced by R using
aov() with an appropriate
Error() designation. However, this is not the same as the sums of squares produced by the
anova() command on the model object, despite the fact that the F-ratios are the same.
Of course this makes complete sense as there is no need for the
Error() strata in a mixed-model. However, this must mean that the sums of squares are penalised somehow in a mixed-model in order to provide appropriate F-ratios. How is this achieved? And how does the model somehow correct the between-plot sum of squares but not correct the within-plot sum of squares. Evidently this is something that is necessary for a classical split-plot ANOVA that was achieved by designating different error values for the different effects, so how does a mixed-effect model allow for this?
Basically, I want to be able to replicate the results from the
anova() command applied to a lmer model object myself to verify the results and my understanding, however, at present I can achieve this for a normal within-subject design but not for the split-plot design and I can't seem to find out why this is the case.
As an example:
library(faraway) library(lme4) data(irrigation) anova(lmer(yield ~ irrigation + variety + (1|field), data = irrigation)) Analysis of Variance Table Df Sum Sq Mean Sq F value irrigation 3 1.6605 0.5535 0.3882 variety 1 2.2500 2.2500 1.5782 summary(aov(yield ~ irrigation + variety + Error(field/irrigation), data = irrigation)) Error: field Df Sum Sq Mean Sq F value Pr(>F) irrigation 3 40.19 13.40 0.388 0.769 Residuals 4 138.03 34.51 Error: Within Df Sum Sq Mean Sq F value Pr(>F) variety 1 2.25 2.250 1.578 0.249 Residuals 7 9.98 1.426
As can be seen above all the F-ratios agree. The sums of squares for variety also agree. However, the sums of squares for irrigation do not agree, however it appears the lmer output is scaled. So what does the anova() command actually do?