R help - Aug 2004 - Enduring LME confusion… or Psychologists and Mixed-Effects

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.

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.
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html

Hello,
>>
>> 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.
>
I suppose that you have, for each subject, enough observations to 
compute his/her average response for each combination of factor1 and 
factor2, no?
If this is the case, you can perform the analysis with the above formula 
on the data obtained by 'aggregate(resp,list(subj,fact1,fact2),mean)'.

This is an analysis with only *within-subject* factors and there 
*cannot* be a problem of unequal number of observation when you have 
only within-subject factors (supposing you have at least one 
observations for each subject in each condition).

I believe the problem with unequal number of observations only  occurs 
when you have at least two crossed *between-subject* (group) variables.

Let's imagine you have two binary group factors (A and B) yielding four 
subgroups of subjects, and for some reason, you do have the same number 
of observations in each subgroup,
Then there are several ways of defining the main effects of A and B.

In many cases, the most reasonable definition of the main effect of A is 
to take the average of A in B1 and in B2 (thus ignoring the number of 
observations, or weithting equally the four subgroups).
To test the null hypothesis of no difference in A when all groups are 
equally weighted, one common approach in psychology is to pretend that 
the number of observation is each group is equal to the harmonic mean of 
the number of observations in each subgroups. The sums of square thud 
obtained can be compared with the error sum of square in the standard 
anova to form an F-test.
This is called the "unweighted" approach.

This can easily be done 'by hand' in R, but there is another approach:

You get equivalent statistics as in the unweighted anova when you use so 
called 'type III' sums of square (I read this in Howell, 1987 
'Statistical methods in psychology',
and in John Fox book 'An R and S-plus companion to appied regression, p. 
140).

It is possible to get type III sums of square using John Fox 'car'
library.

library(car)
contrasts(A)=contr.sum
contrasts(B)=contr.sum
Anova(aov(resp~A*B),type='III')



You can compute the equally weighted cell means defining the effect of A 
with, say:

with(aggregate(resp,list(a=a,b=b),mean),tapply(x,a,mean))


I have seen some people advise against using 'type III' sums of square 
but I do not know their rationale. The important thing, it seems to me, 
is to know
which null hypothesis is  tested in a given test. If indeed the  type 
III sums of square test the effect on equally weighted means, they seem 
okay to me
(when this is indeed the hypothesis I want to test). 

Sorry for not answering any of your questions about the use of 'lme' (I 
hope others will do), but I feel that 'lme' is not needed in the context
of unequal cell frequencies.
(I am happy to be corrected if I am wrong). It seems to me that 'lme' is
useful when some assumptions of standard anova are violate (e.g. with 
repeated measurements when the assumption of sphericity is false), or 
when you have several random factors.

Christophe Pallier
http://www.pallier.org

Dear Gijs:
 
lme fits models using restricted maximum likelihood by default. So, I believe
this is why you have a different DF. If you include method="ML" in the
modeling function the DF should be similar to aov.
 
I think your lme code is incorrect given my understanding of your problem. You
can use all data points using a repeated measures design. Each subject is
measured on multiple occasions. So, observations are nested within individual.
You need to add a time variable to indicate the measurement occasion for each
subject (e.g., t={0,...,T}). You might try the following:
 
fm1.lme<-lme(rt~time*fact1*fact2*fact3, data=mydata, random=~time|sub,
method="ML")
 
Use summary() to get the entire output
 
hope this helps
 
Harold

	-----Original Message----- 
	From: r-help-bounces@stat.math.ethz.ch on behalf of Gijs Plomp 
	Sent: Thu 8/12/2004 4:13 AM 
	To: r-help@stat.math.ethz.ch 
	Cc: 
	Subject: 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 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@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
	
	______________________________________________
	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
	


	[[alternative HTML version deleted]]

On Thu, 12 Aug 2004, Doran, Harold wrote:
>  lme fits models using restricted maximum likelihood by default. So, I
> believe this is why you have a different DF. If you include
method="ML"
> in the modeling function the DF should be similar to aov.
It is REML and not ML that generalizes the classical multistratum AoV.
In any case, REML and ML use the same subspaces, just different methods 
for estimating variances within them.

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