Peter Dalgaard <[EMAIL PROTECTED]> 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. > anova(fit,fit2,fit3,test="Spherical",X=~f1+f3,idata=idata) 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 > anova(fit,fit2,fit3,test="Spherical",M=~f1+f3,X=~f3,idata=idata) Analysis of Variance Table Model 1: Y ~ f2a Model 2: Y ~ 1 Model 3: Y ~ 0 Contrasts orthogonal to ~f3 Contrasts spanned by ~f1 + f3 Greenhouse-Geisser epsilon: 0.5598 Huynh-Feldt epsilon: 1.0284 Res.Df Df Gen.var. F num Df den Df Pr(>F) G-G Pr H-F Pr 1 49 315.74 2 55 6 344.98 0.9883 234 1911 0.54 0.52 0.54 3 56 1 347.15 0.8393 39 1911 0.75 0.68 0.75 ## Same thing with f3 > anova(fit,fit2,fit3,test="Spherical",M=~f1+f3,X=~f1,idata=idata) 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: > anova(fit,fit2,fit3,test="Spherical",M=~f3,X=~1,idata=idata) 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. > anova(fit,fit2,fit3,test="Spherical",M=~1,X=~0,idata=idata) Analysis of Variance Table Model 1: Y ~ f2a Model 2: Y ~ 1 Model 3: Y ~ 0 Contrasts orthogonal to ~0 Contrasts spanned by ~1 Greenhouse-Geisser epsilon: 1 Huynh-Feldt epsilon: 1 Res.Df Df Gen.var. F num Df den Df Pr(>F) G-G Pr H-F Pr 1 49 11 2 55 6 12 1.5010e+00 6 49 1.976e-01 1.976e-01 1.976e-01 3 56 1 249956 1.2699e+06 1 49 8.346e-110 8.346e-110 8.346e-110 ## Finally aov() results for comparison: > system.time(aa <- aov(dv~f1*f2*f3+Error(subject/(f1+f3)),data=d)) [1] 526.50 9.27 561.29 0.00 0.00 > summary(aa) Error: subject Df Sum Sq Mean Sq F value Pr(>F) f2 6 15883 2647 1.501 0.1976 Residuals 49 86411 1763 Error: subject:f1 Df Sum Sq Mean Sq F value Pr(>F) f1 39 81906 2100 0.8393 0.7487 f1:f2 234 578666 2473 0.9883 0.5376 Residuals 1911 4781665 2502 Error: subject:f3 Df Sum Sq Mean Sq F value Pr(>F) f3 3 6899 2300 1.0039 0.3930 f2:f3 18 37153 2064 0.9010 0.5782 Residuals 147 336732 2291 Error: Within Df Sum Sq Mean Sq F value Pr(>F) f1:f3 117 308658 2638 1.0460 0.3500 f1:f2:f3 702 1599743 2279 0.9036 0.9603 Residuals 5733 14458856 2522 > -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html