Fwd: [R] Enduring LME confusion… or Psychologists and Mixed-Effects

Gijs Plomp gplomp at brain.riken.jp
Thu Aug 12 10:13:14 CEST 2004


I will follow the suggestion of John Maindonald and present the problem 
by example with some data.

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           8

Residual 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
g.plomp at brain.riken.jp

 

 >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




More information about the R-help mailing list