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

[Jonathan Baron]

Robert Espesser (Robert.Espesser at lpl.univ-aix.fr) answered the
following question, asking me to post to the list if the answer
worked.  It did.  But ... [see below]

My original message:
----------
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.

---------------

The answer that works:
aov2 <- aov( hl ~ task +Error(subject/task) )
summary(aov2)

yields the correct F test for task.  Presumably this is because
task is a within-subject variable, so task is "nested" (I guess
that is the term) within subject.

But ...

1. It does not work with anova(aov2), just with summary(aov2).

2. I still can't figure out how to get get this to give the
   interaction effect.

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._