I also follow the advice to use mean scores, somewhat reluctantly though. I know it is common practice in psychology, but wouldn’t it be more elegant if one could use all the data points in an analysis?
The data for 5 subjects (myData) are provided at the bottom of this message. It is a crossed within-subject design with three factors and reaction time (RT) as the dependent variable.
An initial repeated-measures model would be: aov1<-aov(RT~fact1*fact2*fact3+Error(sub/(fact1*fact2*fact3)),data=myData)
Aov complains that the effects involving fact3 are unbalanced:
> aov1
…
Stratum 4: sub:fact3
Terms:
fact3 Residuals
Sum of Squares 4.81639e-07 5.08419e-08
Deg. of Freedom 2 8Residual standard error: 7.971972e-05
6 out of 8 effects not estimable Estimated effects may be unbalanced …
Presumably this is because fact3 has three levels and the other ones only two, making this a ‘non-orthogonal’ design.
I then fit an equivalent mixed-effects model: lme1<-lme(RT~fact1*fact2*fact3,data=meanData2,random=~1|sub)
Subsequently I test if my factors had any effect on RT:
> anova(lme1)
numDF denDF F-value p-value(Intercept) 1 44 105294024 <.0001 fact1 1 44 22 <.0001 fact2 1 44 7 0.0090 fact3 2 44 19 <.0001 fact1:fact2 1 44 9 0.0047 fact1:fact3 2 44 1 0.4436 fact2:fact3 2 44 1 0.2458 fact1:fact2:fact3 2 44 0 0.6660
Firstly, why are the F-values in the output whole numbers?
The effects are estimated over the whole range of the dataset and this is so because all three factors are nested under subjects, on the same level. This, however, seems to inflate the F-values compared to the anova(aov1) output, e.g.
> anova(aov1)
…
Error: sub:fact2
Df Sum Sq Mean Sq F value Pr(>F)
fact2 1 9.2289e-08 9.2289e-08 2.2524 0.2078
Residuals 4 1.6390e-07 4.0974e-08
…
I realize that the (probably faulty) aov model may not be directly compared to the lme model, but my concern is if the lme estimation of the effects is right, and if so, how can a naïve skeptic be convinced of this?
The suggestion to use simulate.lme() to find this out seems good, but I can’t seem to get it working (from "[R] lme: simulate.lme in R" it seems it may never work in R).
I have also followed the suggestion to fit the exact same model with lme4. However, format of the anova output does not give me the estimation in the way nlme does. More importantly, the degrees of freedom in the denominator don’t change, they’re still large:
> library(lme4)
> lme4_1<-lme(RT~fact1*fact2*fact3,random=~1|sub,data=myData)
> anova(lme4_1)
Analysis of Variance Table
Df Sum Sq Mean Sq Denom F value Pr(>F) fact1I 1 2.709e-07 2.709e-07 48 21.9205 2.360e-05 ***
fact2I 1 9.229e-08 9.229e-08 48 7.4665 0.008772 **
fact3L 1 4.906e-08 4.906e-08 48 3.9691 0.052047 .
fact3M 1 4.326e-07 4.326e-07 48 34.9972 3.370e-07 ***
fact1I:fact2I 1 1.095e-07 1.095e-07 48 8.8619 0.004552 **
fact1I:fact3L 1 8.988e-10 8.988e-10 48 0.0727 0.788577 fact1I:fact3M 1 1.957e-08 1.957e-08 48 1.5834 0.214351 fact2I:fact3L 1 3.741e-09 3.741e-09 48 0.3027 0.584749 fact2I:fact3M 1 3.207e-08 3.207e-08 48 2.5949 0.113767 fact1I:fact2I:fact3L 1 2.785e-09 2.785e-09 48 0.2253 0.637162 fact1I:fact2I:fact3M 1 7.357e-09 7.357e-09 48 0.5952 0.444206 ---
I hope that by providing a sample of the data someone can help me out on the questions I asked in my previous mail:
>> 1. When aov’s assumptions are violated, can lme provide the right
>> model for within-subjects designs where multiple fixed effects are
>> NOT hierarchically ordered?
>>
>> 2. Are the degrees of freedom in anova(lme1) the right ones to
>> report? If so, how do I convince a reviewer that, despite the large
>> number of degrees of freedom, lme does provide a conservative
>> evaluation of the effects? If not, how does one get the right denDf
>> in a way that can be easily understood?
If anyone thinks he can help me better by looking at the entire data set, I very much welcome them to email me for further discussion.
In great appreciation of your help and work for the R-community,
Gijs Plomp [EMAIL PROTECTED]
>myData sub fact1 fact2 fact3 RT 1 s1 C C G 0.9972709 2 s2 C C G 0.9981664 3 s3 C C G 0.9976909 4 s4 C C G 0.9976047 5 s5 C C G 0.9974346 6 s1 I C G 0.9976206 7 s2 I C G 0.9981980 8 s3 I C G 0.9980503 9 s4 I C G 0.9980620 10 s5 I C G 0.9977682 11 s1 C I G 0.9976633 12 s2 C I G 0.9981558 13 s3 C I G 0.9979286 14 s4 C I G 0.9980474 15 s5 C I G 0.9976030 16 s1 I I G 0.9977088 17 s2 I I G 0.9981506 18 s3 I I G 0.9980494 19 s4 I I G 0.9981183 20 s5 I I G 0.9976804 21 s1 C C L 0.9975495 22 s2 C C L 0.9981248 23 s3 C C L 0.9979146 24 s4 C C L 0.9974583 25 s5 C C L 0.9976865 26 s1 I C L 0.9977107 27 s2 I C L 0.9982071 28 s3 I C L 0.9980966 29 s4 I C L 0.9980372 30 s5 I C L 0.9976303 31 s1 C I L 0.9976152 32 s2 C I L 0.9982363 33 s3 C I L 0.9978750 34 s4 C I L 0.9981402 35 s5 C I L 0.9977018 36 s1 I I L 0.9978076 37 s2 I I L 0.9981699 38 s3 I I L 0.9980628 39 s4 I I L 0.9981092 40 s5 I I L 0.9977054 41 s1 C C M 0.9978842 42 s2 C C M 0.9982752 43 s3 C C M 0.9980277 44 s4 C C M 0.9978250 45 s5 C C M 0.9978353 46 s1 I C M 0.9979674 47 s2 I C M 0.9983277 48 s3 I C M 0.9981954 49 s4 I C M 0.9981703 50 s5 I C M 0.9980047 51 s1 C I M 0.9976940 52 s2 C I M 0.9983019 53 s3 C I M 0.9982484 54 s4 C I M 0.9981784 55 s5 C I M 0.9978177 56 s1 I I M 0.9978636 57 s2 I I M 0.9982188 58 s3 I I M 0.9982024 59 s4 I I M 0.9982358 60 s5 I I M 0.9978581
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
