Re: [R] Tukey HSD (or other post hoc tests) following repeated measures ANOVA

From: Mark Difford <mark_difford_at_yahoo.co.uk>
Date: Fri, 30 May 2008 06:48:20 -0700 (PDT)

Hi Ullrich,

>> # what does '~1 | Subj/Cond' mean?

It is equivalent to your aov() error structure [ ... +Error(Subj/Cond) ]. It gives you a set of random intercepts, one for each level of your nesting structure.

## To get some idea of what's being done it helps to have a continuous covariate in your model.
## So look at this example:
mod1 <- lme(distance ~ age, Orthodont, random = ~ 1 | Subject) ## random intercepts, parallel slopes
mod2 <- lme(distance ~ age, Orthodont, random = ~ age | Subject) ## random intercepts, separate slopes

plot(augPred(mod1, primary=~age)) ## parallel slopes plot(augPred(mod2, primary=~age)) ## separate slopes

In my first post I used the example I did because it's used in V&R's MASS (: 283--286), where the equivalence of the two methods is illustrated. So there really can be no argument about whether lme() can do what aov() does. But lme() [and the newer, still somewhat "hard-hat" version, lmer()] can do a whole lot more than aov() can do. So use lme() [or lmer() if you can't get lme() to converge].

Also, you would be much better off using the multcomp package to do TukeyHSD. Amongst its other advantages are that you can use not only the methods it offers for adjusting p-values, but everything offered in p.adjust().
?p.adjust.

HTH, Mark.

Ullrich Ecker wrote:
>
> Great, thanks a lot, guys!
>
> Only thing is I now have two working versions, that yield *slightly*
> different results.
>
>
> ACCaov <- aov(Acc ~ Cond + Error(Subj/Cond), WMU3C)
> ACClme <- lme(Acc ~ Cond, random = ~1 | Subj/Cond, WMU3C) # what does
> '~1 | Subj/Cond' mean?
> summary(glht(ACClme, linfct=mcp(Cond="Tukey")))
>
> yielding
>
> Linear Hypotheses:
> Estimate Std. Error z value p value
> 2 - 1 == 0 0.392560 0.027210 14.427 <1e-05 ***
> 3 - 1 == 0 0.400372 0.027210 14.714 <1e-05 ***
> 3 - 2 == 0 0.007812 0.025442 0.307 0.95
>
>
> and
>
>
> ACCaov <- aov(Acc ~ Cond + Error(Subj/Cond), WMU3C)
> ACCaov2 <- aov(terms(Acc ~ Subj + Cond, WMU3C)) # gives same result
> as 1st aov, but yields aov not aovlist
> ACCtukey <- TukeyHSD(ACCaov2, "Cond"); ACCtukey
>
> yielding
>
> $Cond
> diff lwr upr p adj
> 2-1 0.3637756 0.29889404 0.42865721 0.0000000
> 3-1 0.3715881 0.30670654 0.43646971 0.0000000
> 3-2 0.0078125 -0.05329192 0.06891692 0.9504659
>
>
> I am trying to work my way through this (so I'm one of the good guys,
> at least I'm trying to understand my stats... ;-)), but any hint to
> what solution may be more appropriate to my problem would be much
> appreciated.
>
> Cheers,
>
> Ullrich
>
>
>
> At 03:16 PM 30/05/2008, you wrote:
>

>>Hi Ullrich,
>>
>> >> The model is
>> >> RT.aov <- aov(RT~Cond + Error(Subj/Cond), WMU3C)
>> >> I understand that TukeyHSD only works with an aov object, but that
>> >> RT.aov is an aovlist object.
>>
>>You want to use lme() in package nlme, then glht() in the multcomp
package.
>>This will give you multiplicity adjusted p-values and confidence intervals.
>>
>>## Example
>>require(MASS)         ## for oats data set
>>require(nlme)         ## for lme()
>>require(multcomp)  ## for multiple comparison stuff
>>
>>Aov.mod <- aov(Y ~ N + V + Error(B/V), data = oats)
>>Lme.mod <- lme(Y ~ N + V, random = ~1 | B/V, data = oats)
>>
>>summary(Aov.mod)
>>anova(Lme.mod)
>>
>>summary(Lme.mod)
>>summary(glht(Lme.mod, linfct=mcp(V="Tukey")))
>>
>>HTH, Mark.
>>
>>
>>Ullrich Ecker wrote:
>> >
>> > Hi everyone,
>> >
>> > I am fairly new to R, and I am aware that others have had this
>> > problem before, but I have failed to solve the problem from previous
>> > replies I found in the archives.
>> >
>> > As this is such a standard procedure in psychological science, there
>> > must be an elegant solution to this...I think.
>> >
>> > I would much appreciate a solution that even I could understand... ;-)
>> >
>> >
>> > Now, I want to calculate a post-hoc test following up a within-subjects
>> > ANOVA.
>> >
>> > The dv is reaction time (RT), there is a 3-level Condition factor
>> > (Cond; within-subjects), a number of subjects (Subj), and the
>> > dataframe is called WMU3C.
>> >
>> > The model is
>> >
>> >  > RT.aov <- aov(RT~Cond + Error(Subj/Cond), WMU3C)
>> >
>> > I understand that TukeyHSD only works with an aov object, but that
>> > RT.aov is an aovlist object.
>> >
>> >  > class(RT.aov)
>> > [1] "aovlist" "listof"
>> >
>> > I've tried to work around it using the "maiz" example in the MMC
>> > documentation of the HH package (a solution previously recommended),
>> > but I couldn't get it to work: My best shot was to calculate another
>> > aov avoiding the error term (I don't see how this could be a feasible
>> > approach, but that's how I understood the MMC example) and a contrast
>> > vector (contrasting conditions 2 and 3):
>> >
>> > I have to admit that I don't quite understand what I'm doing here
>> > (not that you couldn't tell)
>> >
>> >  > RT2.aov <- aov(terms(RT~Subj*Cond, WMU3C))
>> >  > Cond.lmat <- c(0,1,-1)
>> >  > Tukey <- glht.mmc(RT2.aov, focus = "Cond", focus.lmat = Cond.lmat)
>> >
>> > yielding
>> >
>> > Error in mvt(lower = carg$lower, upper = carg$upper, df = df, corr =
>> > carg$corr,  :
>> >    NA/NaN/Inf in foreign function call (arg 6)
>> > In addition: Warning message:
>> > In cov2cor(covm) : diagonal has non-finite entries
>> >
>> >  > Tukey
>> >        height
>> >
>> >
>> >
>> > Thank you very much for your help!
>> >
>> > Ullrich
>> >
>> >
>> > Dr Ullrich Ecker
>> > Postdoctoral Research Associate
>> > Cognitive Science Laboratories
>> > School of Psychology (Mailbag M304)
>> > Room 211 Sanders Building
>> > University of Western Australia
>> > 35 Stirling Hwy
>> > Crawley WA 6009
>> > Australia
>> > Office: +61 8 6488 3266
>> > Mobile: +61 4 5822 0072
>> > Fax: +61 8 6488 1006
>> > E-mail: ullrich.ecker_at_uwa.edu.au
>> > Web: www.cogsciwa.com
>> >       [[alternative HTML version deleted]]
>> >
>> > ______________________________________________
>> > 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.
>> >
>> >
>>
>>--
>>View this message in context: 
>>http://www.nabble.com/Tukey-HSD-%28or-other-post-hoc-tests%29-following-repeated-measures-ANOVA-tp17508294p17553029.html
>>Sent from the R help mailing list archive at Nabble.com.
>>
>>______________________________________________
>>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.

>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> 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.
>
>
-- 
View this message in context: http://www.nabble.com/Tukey-HSD-%28or-other-post-hoc-tests%29-following-repeated-measures-ANOVA-tp17508294p17559307.html
Sent from the R help mailing list archive at Nabble.com.

______________________________________________
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.
Received on Mon 02 Jun 2008 - 04:24:12 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 02 Jun 2008 - 04:30:40 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive