[R] How to get Greenhouse-Geisser epsilons from anova?
John Fox
jfox at mcmaster.ca
Tue Dec 9 00:07:13 CET 2008
Dear Peter and Nils,
I hesitate to repeat this (though I'm going to do it anyway!), but it's
quite simple to get these tests from Anova() in the car package. Here's an
example from ?Anova of a repeated-measures ANOVA with two within and two
between-subject factors:
-------------- snip ---------------
Anova> ## a multivariate linear model for repeated-measures data
Anova> ## See ?OBrienKaiser for a description of the data set used in this
example.
Anova>
Anova> phase <- factor(rep(c("pretest", "posttest", "followup"), c(5, 5,
5)),
Anova+ levels=c("pretest", "posttest", "followup"))
Anova> hour <- ordered(rep(1:5, 3))
Anova> idata <- data.frame(phase, hour)
Anova> idata
phase hour
1 pretest 1
2 pretest 2
3 pretest 3
4 pretest 4
5 pretest 5
6 posttest 1
7 posttest 2
8 posttest 3
9 posttest 4
10 posttest 5
11 followup 1
12 followup 2
13 followup 3
14 followup 4
15 followup 5
Anova> mod.ok <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5,
Anova+ post.1, post.2, post.3, post.4, post.5,
Anova+ fup.1, fup.2, fup.3, fup.4, fup.5) ~
treatment*gender,
Anova+ data=OBrienKaiser)
Anova> (av.ok <- Anova(mod.ok, idata=idata, idesign=~phase*hour))
Type II Repeated Measures MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df Pr(>F)
treatment 2 0.4809 4.6323 2 10 0.0376868 *
gender 1 0.2036 2.5558 1 10 0.1409735
treatment:gender 2 0.3635 2.8555 2 10 0.1044692
phase 1 0.8505 25.6053 2 9 0.0001930
***
treatment:phase 2 0.6852 2.6056 4 20 0.0667354 .
gender:phase 1 0.0431 0.2029 2 9 0.8199968
treatment:gender:phase 2 0.3106 0.9193 4 20 0.4721498
hour 1 0.9347 25.0401 4 7 0.0003043
***
treatment:hour 2 0.3014 0.3549 8 16 0.9295212
gender:hour 1 0.2927 0.7243 4 7 0.6023742
treatment:gender:hour 2 0.5702 0.7976 8 16 0.6131884
phase:hour 1 0.5496 0.4576 8 3 0.8324517
treatment:phase:hour 2 0.6637 0.2483 16 8 0.9914415
gender:phase:hour 1 0.6950 0.8547 8 3 0.6202076
treatment:gender:phase:hour 2 0.7928 0.3283 16 8 0.9723693
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Anova> summary(av.ok, multivariate=FALSE)
Univariate Type II Repeated-Measures ANOVA Assuming Sphericity
SS num Df Error SS den Df F Pr(>F)
treatment 211.286 2 228.056 10 4.6323 0.037687
*
gender 58.286 1 228.056 10 2.5558 0.140974
treatment:gender 130.241 2 228.056 10 2.8555 0.104469
phase 167.500 2 80.278 20 20.8651 1.274e-05
***
treatment:phase 78.668 4 80.278 20 4.8997 0.006426
**
gender:phase 1.668 2 80.278 20 0.2078 0.814130
treatment:gender:phase 10.221 4 80.278 20 0.6366 0.642369
hour 106.292 4 62.500 40 17.0067 3.191e-08
***
treatment:hour 1.161 8 62.500 40 0.0929 0.999257
gender:hour 2.559 4 62.500 40 0.4094 0.800772
treatment:gender:hour 7.755 8 62.500 40 0.6204 0.755484
phase:hour 11.083 8 96.167 80 1.1525 0.338317
treatment:phase:hour 6.262 16 96.167 80 0.3256 0.992814
gender:phase:hour 6.636 8 96.167 80 0.6900 0.699124
treatment:gender:phase:hour 14.155 16 96.167 80 0.7359 0.749562
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Mauchly Tests for Sphericity
Test statistic p-value
phase 0.74927 0.27282
treatment:phase 0.74927 0.27282
gender:phase 0.74927 0.27282
treatment:gender:phase 0.74927 0.27282
hour 0.06607 0.00760
treatment:hour 0.06607 0.00760
gender:hour 0.06607 0.00760
treatment:gender:hour 0.06607 0.00760
phase:hour 0.00478 0.44939
treatment:phase:hour 0.00478 0.44939
gender:phase:hour 0.00478 0.44939
treatment:gender:phase:hour 0.00478 0.44939
Greenhouse-Geisser and Huynh-Feldt Corrections
for Departure from Sphericity
GG eps Pr(>F[GG])
phase 0.79953 7.323e-05 ***
treatment:phase 0.79953 0.01223 *
gender:phase 0.79953 0.76616
treatment:gender:phase 0.79953 0.61162
hour 0.46028 8.741e-05 ***
treatment:hour 0.46028 0.97879
gender:hour 0.46028 0.65346
treatment:gender:hour 0.46028 0.64136
phase:hour 0.44950 0.34573
treatment:phase:hour 0.44950 0.94019
gender:phase:hour 0.44950 0.58903
treatment:gender:phase:hour 0.44950 0.64634
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
HF eps Pr(>F[HF])
phase 0.92786 2.388e-05 ***
treatment:phase 0.92786 0.00809 **
gender:phase 0.92786 0.79845
treatment:gender:phase 0.92786 0.63200
hour 0.55928 2.014e-05 ***
treatment:hour 0.55928 0.98877
gender:hour 0.55928 0.69115
treatment:gender:hour 0.55928 0.66930
phase:hour 0.73306 0.34405
treatment:phase:hour 0.73306 0.98047
gender:phase:hour 0.73306 0.65524
treatment:gender:phase:hour 0.73306 0.70801
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-------------- snip ---------------
The default is so-called "Type-II" tests, but with proper contrast coding
one can get so-called "Type-III" tests as well.
Regards,
John
------------------------------
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On
> Behalf Of Peter Dalgaard
> Sent: December-08-08 1:11 PM
> To: Skotara
> Cc: r-help at r-project.org
> Subject: Re: [R] How to get Greenhouse-Geisser epsilons from anova?
>
> Skotara wrote:
> > Thank you for your help!
> > Sorry, for bothering you again..
> > I still have trouble combining within and between subject factors.
> > Interactions of within factors and D having only 2 levels work well.
> >
> > How can I get the main effect of D? I have tried anova(mlmfitD, mlmfit).
> > With D having 3 levels I would expect the dfs to be 2 and 33. However,
> > the output states 84,24??
>
> That's not a main effect, it's a simultaneous test of all 28 components.
> What you need is an analysis of the average of all
> within-measurements. In principle that is obtained by X=~0, M=~1 or
> T=matrix(1,1,28) but there's a bug that prevents it from working (fixed
> in R-patched a couple of days ago).
>
> >
> > As long as the between factor has only 2 levels the between/within
> > interactions fit well with SPSS, but if D has 3 levels, the mismatch is
> > immense.
> > If I calculate the within effects with myma having not 12 subjects from
> > one group but for example 24 from 2 groups, the output treats it as if
> > all subjects came from the same group, for example for main effect A the
> > dfs are 1 and 35. SPSS puts out 1 and 33 which is what I would have
> > expected.. ..
>
> Hmm, there's a generic problem in that you can't get some of the
> traditional ANOVA table F tests by comparing two models, and in your
> case, SPSS is de facto using the residuals from a model with A:D
> interaction when testing for A. It might help if you try
>
> anova(mlmfitD, X=~..., M=~...)
>
> Look at the (Intercept) line.
>
> >
> >
> > Peter Dalgaard schrieb:
> >> Nils Skotara wrote:
> >>> Thank you, this helped me a lot!
> >>> All within effects and interactions work well!
> >>>
> >>> Sorry, but I still could not get how to include the between factor..
> >>> If I include D with 2 levels, then myma is 24 by 28. (another 12 by
> >>> 28 for the
> >>> second group of subjects.)
> >>> mlmfitD <- lm(myma~D) is no problem, but whatever I tried afterwards
> >>> did not seem logical to me.
> >>> I am afraid I do not understand how to include the between factor. I
> >>> cannot include ~D into M or X because it has length 24 whereas the
other
> >>> factors have 28...
> >>
> >> Just do the same as before, but comparing mlmfitD to mlmfit:
> >>
> >> anova(mlmfitD, mlmfit, X=~A+B, M=~A+B+C)
> >> # or anova(mlmfitD, mlmfit, X=~1, M=~C), as long as things are balanced
> >>
> >>
> >> gives the D:C interaction test (by testing whether the C contrasts
> >> depend on D). The four-factor interaction is
> >>
> >> anova(mlmfitD, mlmfit, X=~(A+B+C)^2, M=~A*B*C)
> >>
> >>
> >>>
> >>>
> >>> Zitat von Peter Dalgaard <p.dalgaard at biostat.ku.dk>:
> >>>
> >>>> Skotara wrote:
> >>>>> Dear Mr. Daalgard.
> >>>>>
> >>>>> thank you very much for your reply, it helped me to progress a bit.
> >>>>>
> >>>>> The following works fine:
> >>>>> dd <- expand.grid(C = 1:7, B= c("r", "l"), A= c("c", "f"))
> >>>>> myma <- as.matrix(myma) #myma is a 12 by 28 list
> >>>>> mlmfit <- lm(myma~1)
> >>>>> mlmfit0 <- update(mlmfit, ~0)
> >>>>> anova(mlmfit, mlmfit0, X= ~C+B, M = ~A+C+B, idata = dd,
> >>>>> test="Spherical"), which tests the main effect of A.
> >>>>> anova(mlmfit, mlmfit0, X= ~A+C, M = ~A+C+B, idata = dd,
> >>>>> test="Spherical"), which tests the main effect of B.
> >>>>>
> >>>>>
> >>>>> However, I can not figure out how this works for the other effects.
> >>>>> If I try:
> >>>>> anova(mlmfit, mlmfit0, X= ~A+B, M = ~A+C+B, idata = dd,
> >>>>> test="Spherical")
> >>>>>
> >>>>> I get:
> >>>>> Fehler in function (object, ..., test = c("Pillai", "Wilks",
> >>>>> "Hotelling-Lawley", :
> >>>>> residuals have rank 1 < 4
> >>>> dd$C is not a factor with that construction. It works for me after
> >>>>
> >>>> dd$C <- factor(dd$C)
> >>>>
> >>>> (The other message is nasty, though. It's slightly different in
> >>>> R-patched:
> >>>>
> >>>> > anova(mlmfit, mlmfit0, X= ~A+B, M = ~A+C+B, idata = dd,
> >>>> test="Spherical")
> >>>> Error in solve.default(Psi, B) :
> >>>> system is computationally singular: reciprocal condition number =
> >>>> 2.17955e-34
> >>>>
> >>>> but it shouldn't happen...
> >>>> Looks like it is a failure of the internal Thin.row function. Ick!
> >>>> )
> >>>>
> >>>>> I also don't know how I can calculate the various interactions..
> >>>>> My read is I should change the second argument mlmfit0, too, but I
> >>>>> can't
> >>>>> figure out how...
> >>>>
> >>>> The "within" interactions should be straightforward, e.g.
> >>>>
> >>>> M=~A*B*C
> >>>> X=~A*B*C-A:B:C
> >>>>
> >>>> etc.
> >>>>
> >>>> The within/between interactions are otained from the similar tests of
> >>>> the between factor(s)
> >>>>
> >>>> e.g.
> >>>>
> >>>> mlmfitD <- lm(myma~D)
> >>>>
> >>>> and then
> >>>>
> >>>> anova(mlmfitD, mlmfit,....)
> >>>>
> >>>>
> >>>>
> >>
> >>
>
>
> --
> O__ ---- Peter Dalgaard Øster Farimagsgade 5, Entr.B
> c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
> (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
> ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
More information about the R-help
mailing list