Split-Plot ANOVA: model comparison tests in R

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=?).

• Hey @caracal, I think the way you are using "Split-Plot Design" is not the way as Casella, George defines it in his book, Statistical Design. Split Plot definitely talks about nesting, but is a special way of imposing the correlation structure. And most of the time you will end up using lme4 package to fit the model AND NOT lm. But this may be a very specific book-based view. I will let other's comment on it. I can give an example based on how I interpret it which is different from yours. – suncoolsu Feb 4 '11 at 1:31
• @suncoolsu The terminology in the social sciences might be different, but both Kirk (1995, p512) and Maxwell & Delaney (2004, p592) call models with one between- and one within-factor "split-plot". The between-factor provides the "plots" (analogous to the agricultural origin). – caracal Feb 4 '11 at 10:46
• I have a lot of things on my plate at the moment. I will expand my answer to be more specific to your question. I see you have invested a lot of effort in framing your question. Thanks for that. – suncoolsu Feb 8 '11 at 22:01

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

• Thank you so much, that is exactly what I was looking for! I was still interested in anova() because of the problem with Anova() described here. But your last suggestion works just as well and is simpler. (Minor thing: I think the last 2 lines are each missing 1 closing parenthesis, and it should read Error(id/(IVw1*IVw2))) – caracal Aug 13 '11 at 12:37

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

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

• I appreciate the intro to analyzing split-splot designs with mixed-effects models and further background info! It certainly is the preferred way to carry out the analysis. I've updated my question to emphasize that I'd still like to know how to do this "the old way". – caracal Feb 4 '11 at 10:51