John Maindonald
2004-Aug-11 22:28 UTC
Fwd: [R] Enduring LME confusion… or Psychologists and Mixed-Effects
In my undertstanding of the problem, the model lme1 <- lme(resp~fact1*fact2, random=~1|subj) should be ok, providing that variances are homogenous both between & within subjects. The function will sort out which factors & interactions are to be compared within subjects, & which between subjects. The problem with df's arises (for lme() in nlme, but not in lme4), when random effects are crossed, I believe. It is difficult to give a general rule on the effect of imbalance; it depends on the relative contributions of the two variances and the nature of the imbalance. There should be a rule that people who ask these sorts of questions are required to make their data available either (if the data set is small) as part of their message or (if data are extensive) on a web site. Once the results of the analysis are on display, it is often possible to make an informed guess on the likely impact. Use of simulate.lme() seems like a good idea. John Maindonald. On 11 Aug 2004, at 8:05 PM, r-help-request at stat.math.ethz.ch wrote:> From: Spencer Graves <spencer.graves at pdf.com> > Date: 10 August 2004 8:44:20 PM > To: Gijs Plomp <gplomp at brain.riken.jp> > Cc: r-help at stat.math.ethz.ch > Subject: Re: [R] Enduring LME confusion.... or Psychologists and > Mixed-Effects > > > Have you considered trying a Monte Carlo? The significance > probabilities for unbalanced anovas use approximations. Package nlme > provides "simulate.lme" to facilitate this. I believe this function > is also mentioned in Pinheiro and Bates (2000). > hope this helps. spencer graves > p.s. You could try the same thing in both library(nlme) and > library(lme4). Package lme4 is newer and, at least for most cases, > better. > Gijs Plomp wrote: > >> Dear ExpeRts, >> >> Suppose I have a typical psychological experiment that is a >> within-subjects design with multiple crossed variables and a >> continuous response variable. Subjects are considered a random >> effect. So I could model >> > aov1 <- aov(resp~fact1*fact2+Error(subj/(fact1*fact2)) >> >> However, this only holds for orthogonal designs with equal numbers of >> observation and no missing values. These assumptions are easily >> violated so I seek refuge in fitting a mixed-effects model with the >> nlme library. >> > lme1 <- lme(resp~fact1*fact2, random=~1|subj) >> >> When testing the 'significance? of the effects of my factors, with >> anova(lme1), the degrees of freedom that lme uses in the denominator >> spans all observations and is identical for all factors and their >> interaction. I read in a previous post on the list ("[R] Help with >> lme basics") that this is inherent to lme. I studied the instructive >> book of Pinheiro & Bates and I understand why the degrees of freedom >> are assigned as they are, but think it may not be appropriate in this >> case. Used in this way it seems that lme is more prone to type 1 >> errors than aov. >> >> To get more conservative degrees of freedom one could model >> > lme2 <- lme(resp~fact1*fact2, random=~1|subj/fact1/fact2) >> >> But this is not a correct model because it assumes the factors to be >> hierarchically ordered, which they are not. >> >> Another alternative is to model the random effect using a matrix, as >> seen in "[R] lme and mixed effects" on this list. >> > lme3 <- (resp~fact1*fact2, random=list(subj=pdIdent(form=~fact1-1), >> subj=~1, fact2=~1) >> >> This provides 'correct? degrees of freedom for fact1, but not for the >> other effects and I must confess that I don't understand this use of >> matrices, I?m not a statistician. >> >> My questions thus come down to this: >> >> 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? >> >> I hope that my confusion is all due to an ignorance of statistics and >> that someone on this list will kindly point that out to me. I do >> realize that this type of question has been asked before, but think >> that an illuminating answer can help R spread into the psychological >> community. >>John Maindonald email: john.maindonald at anu.edu.au phone : +61 2 (6125)3473 fax : +61 2(6125)5549 Centre for Bioinformation Science, Room 1194, John Dedman Mathematical Sciences Building (Building 27) Australian National University, Canberra ACT 0200.
Gijs Plomp
2004-Aug-12 08:13 UTC
Fwd: [R] Enduring LME confusion… or Psychologists and Mixed-Effects
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 wouldnt 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 nave skeptic be convinced of
this?
The suggestion to use simulate.lme() to find this out seems good, but I
cant 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 dont change, theyre 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 aovs 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