[R] Within-Subjects ANOVA & comparisons of individual means
Prof Brian Ripley
ripley at stats.ox.ac.uk
Wed Jan 25 14:27:23 CET 2006
On Tue, 24 Jan 2006, Spencer Graves wrote:
> 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.
It is a lot more complicated than that. You can probably adjust TukeyHSD
to work on a single stratum of an "aovlist" object, but it is not so
obvious that you can retrieve the information needed from the aovlist
object so you may need to refit. However, my reason for not pursuing that
is more philosophical: these post hoc tests are intended to allow for
multiple testing, and that would only make sense if you allow for all the
tests done in all the strata. But here there is only stratum, and that
suggests the wrong analysis is being done.
As I understand this design, it is a randomized block design with subject
as block and interval as treatment. The classic analysis is a
fixed-effect one (subject + interval), and TukeyHSD can be applied to
that. E.g.
> fm <- aov(dv ~ subject + interval)
> summary(fm)
Df Sum Sq Mean Sq F value Pr(>F)
subject 6 2.4853 0.4142 1.3067 0.2845
interval 5 13.8174 2.7635 8.7178 3.417e-05 ***
Residuals 30 9.5098 0.3170
(should might look familiar)
> TukeyHSD(fm, "interval")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = dv ~ subject + interval)
$interval
diff lwr upr p adj
2-1 -0.7477000 -1.663060 0.16766002 0.1608612
3-1 -1.0704286 -1.985789 -0.15506855 0.0145883
4-1 -1.4761143 -2.391474 -0.56075426 0.0004035
5-1 -1.5005143 -2.415874 -0.58515426 0.0003225
6-1 -1.6827143 -2.598074 -0.76735426 0.0000601
3-2 -0.3227286 -1.238089 0.59263145 0.8884095
4-2 -0.7284143 -1.643774 0.18694574 0.1814407
5-2 -0.7528143 -1.668174 0.16254574 0.1557224
6-2 -0.9350143 -1.850374 -0.01965426 0.0430664
4-3 -0.4056857 -1.321046 0.50967431 0.7563642
5-3 -0.4300857 -1.345446 0.48527431 0.7095947
6-3 -0.6122857 -1.527646 0.30307431 0.3475861
5-4 -0.0244000 -0.939760 0.89096002 0.9999994
6-4 -0.2066000 -1.121960 0.70876002 0.9821152
6-5 -0.1822000 -1.097560 0.73316002 0.9898235
However, I find it suspicious that your treatment effects are in fact
ordered, and suspect there is something more to the squeezed out here.
You appear to be defining difference contrasts without using contr.sdif
(in MASS). Once you have a single-stratum fit, a lot more tools become
available (se.contrast is one). But for starters
> fm2 <- lm(dv ~ subject+interval, contrasts=list(interval="contr.sdif"))
> summary(fm2)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
...
interval2-1 -0.74770 0.30095 -2.484 0.0188 *
interval3-2 -0.32273 0.30095 -1.072 0.2921
interval4-3 -0.40569 0.30095 -1.348 0.1877
interval5-4 -0.02440 0.30095 -0.081 0.9359
interval6-5 -0.18220 0.30095 -0.605 0.5495
I hope that is enough food for thought (it is not an invitation for
further free consultancy).
> 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
>
> ______________________________________________
> 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
>
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
More information about the R-help
mailing list