[R] aov for mixed model (fixed and random)?

Jonathan Baron baron at cattell.psych.upenn.edu
Sun Dec 23 18:26:28 CET 2001

I'm starting to understand fixed and random effects, but I'm
puzzled a bit.  Here is an example from Hays's textbook (which is
great at explaining fixed vs. random effects, at least to dummies
like me), from the section on mixed models.  You need
library(nlme) in order to run it.
task <- gl(3,2,36) # Three tasks, a fixed effect.
subj <- gl(6,6,36) # Six subjects, a random effect.

h1 <- c(7.8,8.7,  11.1,12.0, 11.7,10.0, # Each S does each task twice.
        8.0,9.2,  11.3,10.6, 9.8,11.9,
        4.0,6.9,   9.8,10.1, 11.7,12.6,
        10.3,9.4, 11.4,10.5, 7.9,8.1,
        9.3,10.6, 13.0,11.7, 8.3,7.9,
        9.5,9.8,  12.2,12.3, 8.6,10.5)

aov1 <- aov(h1~task*subj)
anova(aov1) # See note below.

lme1 <- lme(h1~task,random=~1|subj)
lme2 <- lme(h1~task,random=~1|subj/task)
anova(lme1,lme2) # for interaction
anova(lme2) # for effect of task

The anova gives very close to the correct F values for subj
(1.68, according to Hays) and for the interaction (7.19), but the
wrong F values for task, because it treats task as a random
effect, using the wrong error term.  Specifically, it gives
F=27.27, instead of 3.78 (Hays's figure).  This is because it
uses the within-cell error as the denominator for F instead of
the interaction term.  (The interaction ms is 5.82.)

The lme results give approximately the same results as Hays, but
they are different, presumably because of using a different
method (REML).  (I still don't understand why I can't just say
anova(lme2) and get everything all at once, but that's another
matter, and at least I've finally gotten this far.)

Now my question: Can I reproduce the results in Hays exactly?
Presumably this will involve instructing aov to use the
interaction as the error term for the denominator of F.  I've
tried everything I can think of.  If you can answer this, it
might go a long way toward helping me understand the (still
somewhat mysterious) "Error()" component of aov, if that turns
out to be part of the answer.

Jon Baron
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch

More information about the R-help mailing list