What does the anova() command do with a lmer model object?

Question

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?

Answer 1

Use the Source, Luke. We can peek inside the ANOVA function by doing getAnywhere(anova.Mermod). The first part of that function is for comparing two different models. The anova on the fixed effects comes in the big else block in the second half:

 dc <- getME(object, "devcomp")
        X <- getME(object, "X")
        asgn <- attr(X, "assign")
        stopifnot(length(asgn) == (p <- dc$dims[["p"]]))
            ss <- as.vector(object@pp$RX() %*% object@beta)^2
        names(ss) <- colnames(X)
        terms <- terms(object)
        nmeffects <- attr(terms, "term.labels")[unique(asgn)]
        if ("(Intercept)" %in% names(ss)) 
            nmeffects <- c("(Intercept)", nmeffects)
        ss <- unlist(lapply(split(ss, asgn), sum))
        stopifnot(length(ss) == length(nmeffects))
        df <- vapply(split(asgn, asgn), length, 1L)
        ms <- ss/df
        f <- ms/(sigma(object)^2)
        table <- data.frame(df, ss, ms, f)
        dimnames(table) <- list(nmeffects, c("Df", "Sum Sq", 
            "Mean Sq", "F value"))
        if ("(Intercept)" %in% nmeffects) 
            table <- table[-match("(Intercept)", nmeffects), 
                ]
        attr(table, "heading") <- "Analysis of Variance Table"
        class(table) <- c("anova", "data.frame")
        table

object is the lmer output. We start to calculate the sum of squares in line 5: ss <- as.vector ... The code multiplies the fixed parameters (in beta) by an upper triangular matrix; then squares each term. Here is that upper triangular matrix for the irrigation example. Each row corresponds to one of the five fixed effects parameters (intercept, 3 degrees of freedom for irrigation, 1 df for variety).

zapsmall(irrigation.lmer@pp$RX(), digits = 3)
      [,1]  [,2]   [,3]   [,4]  [,5]
[1,] 0.813 0.203  0.203  0.203 0.407
[2,] 0.000 0.352 -0.117 -0.117 0.000
[3,] 0.000 0.000  0.332 -0.166 0.000
[4,] 0.000 0.000  0.000  0.287 0.000
[5,] 0.000 0.000  0.000  0.000 2.000

The first row gives you the sum of squares for the intercept and the last gives you the SS for the within-fields variety effect. Rows 2-4 only involve the 3 parameters for the irrigation levels, so pre-multiplication gives you three pieces of the SS for irrigation.

These pieces are not interesting in themselves since they come from the default treatment contrast in R, but in line ss <- unlist(lapply(split .... Bates scoops up bits of sums of squares according to the number of levels and which factors they refer to. There's a lot of book-keeping going on here. We get the degrees of freedom as well (which are 3 for irrigation). Then, he gets the mean squares which show up on the printout of anova. Finally, he divides all his means squares by the within-groups residual variance, sigma(object)^2.

So what's going on? The philosophy of lmer has nothing to do with the method of moments approach used by aov. The idea in lmer is to maximize a marginal likelihood obtained by integrating out the unseen random effects. In this case, the random fertility level of each field. Chapter 2 of Pinheiro and Bates describes the ugliness of this process. The RX matrix used to get the sums of squares is their matrix $R_{00}$ from equation 2.17, page 70 of the text. This matrix is obtained from a bunch of QR decompositions upon (amongst other things) the design matrix of the random effects and $\sqrt{\sigma^2/\sigma_f^2}$, where $\sigma_f^2$ is the variance of the field effect. This is that missing factor you were asking about, but it does not enter into the solution in a transparent or simple way.

Asymptotically, the estimates of the fixed effects have distribution:

$$ \hat{\beta} \sim \mathcal{N}(\beta, \sigma^2[R_{00}^{-1}R_{00}^{-T}]) $$

which means that the components of $R_{00}\hat{\beta}$ are independently distributed. If $\beta=0$, the squares of those terms (divided by $\sigma^2$) follow a chi-squared distribution. Dividing through by the estimate of $\sigma^2$ (another chi-squared when divided by $\sigma^2$) gives an F statistic. We don't divide by the mean square error for the between-groups analysis because that has nothing to do with what's going on here. We need to compare the sums of squares obtained through $R_{00}$ by comparing them to an estimate of $\sigma^2$.

Note that you would not have got the same F statistics had the data been unbalanced. Nor would you have obtained the same F statistics if you had used ML instead of REML.

The idea behind aov is that the expected mean squares for irrigation is a function of $\sigma^2$, $\sigma_f^2$ and the irrigation effects. The expected mean square for field residuals is a function of $\sigma^2$ and $\sigma_f^2$. When the irrigation effects are 0, both of these quantities estimate the same thing and their ratio follows an F distribution.

Interestingly, Bates and Pinheiro recommend using the ANOVA over fitting two models and doing a likelihood ratio test. The latter tends to be anti-conservative.

Of course, if the data are unbalanced, it is no longer clear exactly what hypothesis the ANOVA is testing. I removed the first observation from the irrigation data and refit. Here is the $R_{00}$ matrix again:

zapsmall(fit2@pp$RX(), digits = 3)
      [,1]  [,2]   [,3]   [,4]   [,5]
[1,] 0.816 0.205  0.205  0.205  0.457
[2,] 0.000 0.354 -0.119 -0.119 -0.029
[3,] 0.000 0.000  0.334 -0.168 -0.040
[4,] 0.000 0.000  0.000  0.288 -0.071
[5,] 0.000 0.000  0.000  0.000  1.874

As you can see, the sums of squares for the irrigation parameters now contain some of the variety effect as well.