Re: [R] Within-Subjects ANOVA & comparisons of individual means

From: Prof Brian Ripley <ripley_at_stats.ox.ac.uk>
Date: Thu 26 Jan 2006 - 00:27:23 EST

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@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@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@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

______________________________________________
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
Received on Thu Jan 26 00:44:13 2006

This archive was generated by hypermail 2.1.8 : Thu 26 Jan 2006 - 02:13:04 EST