Peter Dalgaard <p.dalgaard at biostat.ku.dk> writes:
- (not too sure of this, but R 2.1.0 will have some new code for
multivariate lm, with intra-individual designs, and
tests under the sphericity assumptions; it is possible that
your models can be reformulated in that framework. You'd have to
set it up as a model with 160 responses on 56 subjects, though, and
the code wasn't really designed with that sort of intrasubject
dimensionality in mind.)
OK, I've tried this now and it does seem to work, and much faster than
the aov() approach. I had to fix a buglet in the code (eigenvalues
coming up with small imaginary parts due to numerics), so currently
available versions won't quite work, but it should be in the alpha
releases that are supposed to start Monday.
Here's how it works:
### orig setup, edited so as to actually work...
f1<-gl(40,1,8960,ordered=T)
f2<-gl(7,1280,8960,ordered=T)
f3<-gl(4,40,8960,ordered=T)
subject<-gl(56,160,8960,ordered=T)
dv<-rnorm(8960,mean=500,sd=50)
d <- data.frame(f1,f2,f3,subject,dv)
### Regroup to 56x160 matrix response
f2a <- f2[seq(1,8801,160)]
idata <- d[1:160,] # intrasubj. design, actually need only f1,f3 from this
Y <- matrix(dv,56,byrow=T)
### now fit models with effect of f2a, constant, and empty
fit <- lm(Y~f2a)
fit2 <- lm(Y~1)
fit3 <- lm(Y~0)
## The main trick is to do anova on within-subject effects. First look
## at the interactions, alias residuals from an additive model ~f1+f3.
## The reduction Model 1-> Model 2 corresponds to testing the f1:f2:f3
## interaction in aov (tests whether the f1:f3 interaction depends on
## f2a) and Model 2 -> Model 3 is the test for the f1:f3 interaction
## being zero.
## I'm not quite sure what to make of the correction terms when the
## estimated SSD matrix is singular (it has to be, the dimension is
## 117, but the df is only 49). Probably the G-G epsilon is bogus, but
## the H-F one seems to fare rather well.
Analysis of Variance Table
Model 1: Y ~ f2a
Model 2: Y ~ 1
Model 3: Y ~ 0
Contrasts orthogonal to
~f1 + f3
Greenhouse-Geisser epsilon: 0.2903
Huynh-Feldt epsilon: 0.9639
Res.Df Df Gen.var. F num Df den Df Pr(>F) G-G Pr H-F Pr
1 49 0.00000
2 55 6 0.00000 0.9036 702 5733 0.96034 0.82268 0.95748
3 56 1 0.00000 1.0460 117 5733 0.35003 0.39624 0.35206
## Next, we can look at the f1 effects. This is basically the
## difference between two projections, on ~f1+f3 and ~f3
## This gives us the tests for f2:f1 and f1
Analysis of Variance Table
Model 1: Y ~ f2a
Model 2: Y ~ 1
Model 3: Y ~ 0
Contrasts orthogonal to
~f1
Contrasts spanned by
~f1 + f3
Greenhouse-Geisser epsilon: 0.9482
Huynh-Feldt epsilon: 1.0128
Res.Df Df Gen.var. F num Df den Df Pr(>F) G-G Pr H-F Pr
1 49 35.171
2 55 6 34.679 0.9010 18 147 0.578 0.574 0.578
3 56 1 34.658 1.0039 3 147 0.393 0.390 0.393
## Actually, because of the complete design, we can ignore f1 and get
## the same analysis:
Analysis of Variance Table
Model 1: Y ~ f2a
Model 2: Y ~ 1
Model 3: Y ~ 0
Contrasts orthogonal to
~1
Contrasts spanned by
~f3
Greenhouse-Geisser epsilon: 0.9482
Huynh-Feldt epsilon: 1.0128
Res.Df Df Gen.var. F num Df den Df Pr(>F) G-G Pr H-F Pr
1 49 35.171
2 55 6 34.679 0.9010 18 147 0.578 0.574 0.578
3 56 1 34.658 1.0039 3 147 0.393 0.390 0.393
## Finally we get the main effect of f2a as follows. Notice that the
## Model 2 -> Model 3 reduction is the test for zero overall mean,
## which you might (and aov does) want to omit.