Split-Plot ANOVA: model comparison tests in R

Question

How can I test effects in a Split-Plot ANOVA using suitable model comparisons for use with the X and M arguments of anova.mlm() in R? I'm familiar with ?anova.mlm and Dalgaard (2007)[1]. Unfortunately it only brushes Split-Plot Designs. Doing this in a fully randomized design with two within-subjects factors:

N  <- 20  # 20 subjects total
P  <- 3   # levels within-factor 1
Q  <- 3   # levels within-factor 2
DV <- matrix(rnorm(N* P*Q), ncol=P*Q)           # random data in wide format
id <- expand.grid(IVw1=gl(P, 1), IVw2=gl(Q, 1)) # intra-subjects layout of data matrix

library(car)        # for Anova()
fitA <- lm(DV ~ 1)  # between-subjects design: here no between factor
resA <- Anova(fitA, idata=id, idesign=~IVw1*IVw2)
summary(resA, multivariate=FALSE, univariate=TRUE)  # all tests ...

The following model comparisons lead to the same results. The restricted model doesn't include the effect in question but all other effects of the same order or lower, the full model adds the effect in question.

anova(fitA, idata=id, M=~IVw1 + IVw2, X=~IVw2, test="Spherical") # IVw1
anova(fitA, idata=id, M=~IVw1 + IVw2, X=~IVw1, test="Spherical") # IVw2
anova(fitA, idata=id, M=~IVw1 + IVw2 + IVw1:IVw2,
                      X=~IVw1 + IVw2, test="Spherical")          # IVw1:IVw2

A Split-Splot design with one within and one between-subjects factor:

idB  <- subset(id, IVw2==1, select="IVw1")          # use only first within factor
IVb  <- gl(2, 10, labels=c("A", "B"))               # between-subjects factor
fitB <- lm(DV[ , 1:P] ~ IVb)                        # between-subjects design
resB <- Anova(fitB, idata=idB, idesign=~IVw1)
summary(resB, multivariate=FALSE, univariate=TRUE)  # all tests ...

These are the anova() commands to replicate the tests, but I don't know why they work. Why do the tests of the following model comparisons lead to the same results?

anova(fitB, idata=idB, X=~1, test="Spherical") # IVw1, IVw1:IVb
anova(fitB, idata=idB, M=~1, test="Spherical") # IVb

Two within-subjects factors and one between-subjects factor:

fitC <- lm(DV ~ IVb)  # between-subjects design
resC <- Anova(fitC, idata=id, idesign=~IVw1*IVw2)
summary(resC, multivariate=FALSE, univariate=TRUE)  # all tests ...

How do I replicate the results given above with the corresponding model comparisons for use with the X and M arguments of anova.mlm()? What is the logic behind these model comparisons?

EDIT: suncoolsu pointed out that for all practical purposes, data from these designs should be analyzed using mixed models. However, I'd still like to understand how to replicate the results of summary(Anova()) with anova.mlm(..., X=?, M=?).

[1]: Dalgaard, P. 2007. New Functions for Multivariate Analysis. R News, 7(2), 2-7.

Answer 1

The X and M basically specify the two models you want to compare, but only in terms of the within-subject effects; it then shows results for the interaction of all between-subject effects (including the intercept) with the within-subject effects that changed between X and M.

Your examples on fitB are easier to understand if we add the defaults for X and M:

anova(fitB, idata=idB, M=~1, X=~0, test="Spherical") # IVb
anova(fitB, idata=idB, M=diag(3), X=~1, test="Spherical") # IVw1, IVw1:IVb

The first model is the change from no within subject effects (all have the same mean) to a different mean for each, so we've added the id random effect, which is the right thing to test the overall intercept and the overall between subject effect on.

The second model ads the id:IVw1 interaction, which is the right thing to test IVw1 and the IVw1:IVb terms against. Since there is only one within-subject effect (with three levels) the default of diag(3) in the second model will account for it; it would be equivalent to run

anova(fitB, idata=idB, M=~IVw1, X=~1, test="Spherical") # IVw1, IVw1:IVb

For your fitC, I believe these commands will recreate the Anova summary.

anova(fitC, idata=id, M=~1, X=~0, test="Spherical") #IVb
anova(fitC, idata=id, M=~IVw1 + IVw2, X=~IVw2, test="Spherical") # IVw1
anova(fitC, idata=id, M=~IVw1 + IVw2, X=~IVw1, test="Spherical") # IVw2
anova(fitC, idata=id, M=~IVw1 + IVw2 + IVw1:IVw2,
                  X=~IVw1 + IVw2, test="Spherical")          # IVw1:IVw2

Now, as you discovered, these commands are really tricky. Thankfully, there's not much reason to use them anymore. If you're willing to assume sphericity, you should just use aov, or for even easier syntax, just use lm and compute the right F-tests yourself. If you're not willing to assume sphericity, using lme is really the way to go as you get a lot more flexibility than you do with the GG and HF corrections.

For example, here's the aov and lm code for your fitA. You do need to have the data in long format first; here's one way to do that:

library(reshape)
d0 <- data.frame(id=1:nrow(DV), DV)
d0$IVb <- IVb
d0 <- melt(d0, id.vars=c(1,11), measure.vars=2:10)
id0 <- id
id0$variable <- factor(levels(d0$variable), levels=levels(d0$variable))
d <- merge(d0, id0)
d$id <- factor(d$id)

And here's lm andaov` code:

anova(lm(value ~ IVw1*IVw2*id, data=d))
summary(aov(value ~ IVw1*IVw2 + Error(id/(IVw1*IVw2)), data=d))

Answer 2

Split-plot designs originated in agriculture, hence the name. But they frequently occur and I would say -- the workhorse of most of the clinical trials. The main plot is treated with a level of one factor while the levels of some other factor are allowed to vary with the subplots. The design arises as a result of restrictions on a full randomization. For example: a field may be divided into four subplots. It may be possible to plant different varieties in subplots, but only one type of irrigation may be used for the whole field. Not the distinction between splits and blocks. Blocks are features of the experimental units which we have the option to take advantage of in the experimental design, because we know they are there. Splits, on the other hand, impose restriction on what assignments of factors are possible. They impose requirements on the design that prevent a complete randomization.

They are used a lot in clinical trials where when one factor is easy to change while another factor takes much more time to change. If the experimenter must do all runs for each level of the hard-to-change factor consecutively, a split plot design results with the hard-to-change factor representing the whole plot factor.

Here is an example: In an agricultural field trial, the objective was to determine the effects of two crop varieties and four different irrigation methods. Eight fields were available, but only one type of irrigation may be applied to each field. The fields may be divided into two parts with a different variety in each part. The whole plot factor is the irrigation, which should be randomly assigned to fields. Within each field, the variety is assigned.

This is how you do this in R:

install.packages("faraway")
data(irrigation)
summary(irrigation)

library(lme4)

R> (lmer(yield ~ irrigation * variety + (1|field), data = irrigation))
Linear mixed model fit by REML 
Formula: yield ~ irrigation * variety + (1 | field) 
   Data: irrigation 
  AIC  BIC logLik deviance REMLdev
 65.4 73.1  -22.7     68.6    45.4
Random effects:
 Groups   Name        Variance Std.Dev.
 field    (Intercept) 16.20    4.02    
 Residual              2.11    1.45    
Number of obs: 16, groups: field, 8

Fixed effects:
                       Estimate Std. Error t value
(Intercept)               38.50       3.02   12.73
irrigationi2               1.20       4.28    0.28
irrigationi3               0.70       4.28    0.16
irrigationi4               3.50       4.28    0.82
varietyv2                  0.60       1.45    0.41
irrigationi2:varietyv2    -0.40       2.05   -0.19
irrigationi3:varietyv2    -0.20       2.05   -0.10
irrigationi4:varietyv2     1.20       2.05    0.58

Correlation of Fixed Effects:
            (Intr) irrgt2 irrgt3 irrgt4 vrtyv2 irr2:2 irr3:2
irrigation2 -0.707                                          
irrigation3 -0.707  0.500                                   
irrigation4 -0.707  0.500  0.500                            
varietyv2   -0.240  0.170  0.170  0.170                     
irrgtn2:vr2  0.170 -0.240 -0.120 -0.120 -0.707              
irrgtn3:vr2  0.170 -0.120 -0.240 -0.120 -0.707  0.500       
irrgtn4:vr2  0.170 -0.120 -0.120 -0.240 -0.707  0.500  0.500

Basically, what this model says is, irrigation and variety are fixed effects and variety is nested within irrigation. The fields are the random effects and pictorially will be something like

I_1 | I_2 | I_3 | I_4

V_1 V_2 | V_1 V_2 | V_1 V_2 | V_1 V_2

But this was a special variant with fixed whole plot effect and subplot effect. There can be variants in which one or more are random. There can be more complicated designs like split-split .. plot designs. Basically, you can go wild and crazy. But given the underlying structure and distribution ( ie fixed or random, nested or crossed, .. ) is clearly understood, a lmer-Ninja will have no troubles in modeling. May be interpretation will be a mess.

Regarding comparisons, say you have lmer1 and lmer2:

anova(lmer1, lmer2)

will give you the appropriate test based on the chi-sq test statistic with degrees of freedom equal to the difference of parameters.

cf: Faraway, J., Extending Linear Models with R.

Casella, G., Statistical Design