[R] Within-Subjects ANOVA & comparisons of individual means
Spencer Graves
spencer.graves at pdf.com
Tue Jan 24 17:52:21 CET 2006
1. Did you try "summary(aovRes, ...)" rather than
"summary.aov(aovRes, ...)"? From "summary.aov", I got the same error
message you did, but from "summary(aovRes, ...)", I got something that
looked like what you were expecting.
2. To understand this, look at the code for "summary": It consists
essentially of 'UseMethod("summary")'. To trace methods dispatch here,
note that 'class(aovRes)' is a 2-vector: c("aovlist", "listof"). When
I requested 'methods("summary"), I got a long list, the first two were
'summary.aov' and 'summary.aovlist'. Since aovRes was NOT of class
'aov' but instead of class 'aovlist', summary(aovRes, ...) uses
'summary.aovlist'. You got an error message, because summary.aov was
expecting an argument of class 'aov', and 'aovRes' didn't have the
structure it required.
3. To find an attribute "contrasts" in aovRes, I requested
'str(aovRes)' and searched for "contrasts". I found it in two places:
attr(attr(aovRes, "error.qr")$qr, "contrasts")
attr(aovRes, "contrasts")
4. I wondered why you said, "From my understanding, TukeyHSD is not
appropriate in this context." Then I discovered that TukeyHSD is
officially a generic function with a method only for objects of class
"aov". Since aovRes is NOT of class "aov", it can't use that function.
However, it looks to me like the same algorithm could be used for what
you want, e.g., by copying "TukeyHSD.aov" into a script file, changing
the name to "TukeyHSD.aovlist" and then walking through the code line by
line, and changing things as necessary to produce the desired results.
If I wanted to do this, I might first check with the author of TukeyHSD
(Douglas Bates). He might know some reason why it is inappropriate or
some other special concerns of which I'm unaware. Or he might have a
not-quite-fully debugged version someplace that could help you.
Best Wishes,
Spencer Graves
Steffen Katzner wrote:
> I am having problems with comparing individual means in a
> within-subjects ANOVA. From my understanding, TukeyHSD is not
> appropriate in this context. So I am trying to compute contrasts, as
> follows:
>
> seven subjects participated in each of 6 conditions (intervals).
>
>>subject = factor(rep(c(1:7), each = 6))
>>interval = factor(rep(c(1:6), 7))
>
> and here is the dependent variable:
>
>>dv = c(3.3868, 3.1068, 1.7261, 1.5415, 1.7356, 0.7538,
>
> + 2.5957, 1.5666, 1.1984, 1.2761, 1.0022, 0.8597,
> + 3.9819, 3.1506, 1.5824, 1.7400, 1.4248, 0.6519,
> + 2.2521, 1.5248, 1.1209, 1.2193, 1.1994, 2.0910,
> + 2.4661, 1.3863, 1.3591, 0.9163, 1.3976, 1.7471,
> + 3.2486, 1.9492, 2.4228, 1.1276, 1.2836, 0.9814,
> + 1.7148, 1.7278, 2.7433, 1.4924, 1.0992, 0.7821)
>
>
>>d = data.frame(subject, interval, dv)
>
> next I'm defining a contrast matrix:
>
>>con = matrix(c(1, -1, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 1, -1, 0,
>
> 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 1, -1), nrow=6, ncol=5, byrow=F)
>
>>con
> [,1] [,2] [,3] [,4] [,5]
> [1,] 1 0 0 0 0
> [2,] -1 1 0 0 0
> [3,] 0 -1 1 0 0
> [4,] 0 0 -1 1 0
> [5,] 0 0 0 -1 1
> [6,] 0 0 0 0 -1
>
>>contrasts(d$interval)=con
>
> and then I'm doing the ANOVA
>
>>aovRes = aov(dv~interval+Error(subject/interval), data=d)
>
>>summary(aovRes)
>
> Error: subject
> Df Sum Sq Mean Sq F value Pr(>F)
> Residuals 6 2.48531 0.41422
>
> Error: subject:interval
> Df Sum Sq Mean Sq F value Pr(>F)
> interval 5 13.8174 2.7635 8.7178 3.417e-05 ***
> Residuals 30 9.5098 0.3170
> ---
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> but if I want to look at the contrasts, something has gone wrong:
>
> summary.aov(aovRes, split=list(interval = list("i1 vs i2" = 1, "i2 vs
> i3" = 2, "i3 vs i4" = 3, "i4 vs i5" = 4, "i5 vs i6" = 5)))
>
> Error in 1:object$rank : NA/NaN argument
>
>>aovRes$contrasts
>
> NULL
>
> Can anybody help?
> Thank you very much, -Steffen
>
> ______________________________________________
> R-help at 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
More information about the R-help
mailing list