[R] degrees of freedom in lme

[Manuel Morales]

On Mon, 2007-06-25 at 13:15 -0400, Doran, Harold wrote:
> This is such a common question that it has a an "FAQ-like" response from Doug Bates. Google "lmer p-values and all that" to find the response. 

Isn't this a different question, though, since Jean-Baptiste is using
nlme. 

Details on the calculation of DF in nlme can be found in chapter 4 of
the book by Pinheiro and Bates "Mixed Effects Models in S and S-PLUS.
Using the formula provided, I get denDF of 10 for level 1 and 32 for
level 2. I'm not sure why lme is using the denDF estimated at level 2 in
this example ...

> > -----Original Message-----
> > From: r-help-bounces at stat.math.ethz.ch 
> > [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of 
> > Jean-Baptiste Ferdy
> > Sent: Monday, June 25, 2007 12:26 PM
> > To: r-help at stat.math.ethz.ch
> > Subject: [R] degrees of freedom in lme
> > 
> > Dear all,
> > 
> > I am starting to use the lme package (and plan to teach a 
> > course based on it next semester...). To understand what lme 
> > is doing precisely, I used balanced datasets described in 
> > Pinheiro and Bates and tried to compare the lme outputs to 
> > that of aov. Here is what I obtained:
> > 
> > > data(Machines)
> > > summary(aov(score~Machine+Error(Worker/Machine),data=Machines))
> > Error: Worker
> >           Df  Sum Sq Mean Sq F value Pr(>F) Residuals  5 
> > 1241.89  248.38
> > 
> > Error: Worker:Machine
> >           Df  Sum Sq Mean Sq F value    Pr(>F)
> > Machine    2 1755.26  877.63  20.576 0.0002855 ***
> > Residuals 10  426.53   42.65
> > ---
> > Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> > 
> > Error: Within
> >           Df Sum Sq Mean Sq F value Pr(>F)
> > Residuals 36 33.287   0.925
> > > 
> > anova(lme(fixed=score~Machine,random=~1|Worker/Machine,data=Machines))
> >             numDF denDF  F-value p-value
> > (Intercept)     1    36 773.5709  <.0001
> > Machine         2    10  20.5762   3e-04
> >   
> > No problem here: the results are essentially the same, which 
> > is expected. Now I turn to an ANCOVA with a random grouping factor.
> > 
> > > data(Orthodont)
> > > OrthoFem <- Orthodont[Orthodont$Sex=="Female",];
> > > summary(aov(distance~age+Error(Subject/age),data=OrthoFem))
> > Error: Subject
> >           Df  Sum Sq Mean Sq F value Pr(>F) Residuals 10 
> > 177.227  17.723
> > 
> > Error: Subject:age
> >           Df Sum Sq Mean Sq F value    Pr(>F)
> > age        1 50.592  50.592  52.452 2.783e-05 ***
> > Residuals 10  9.645   0.965
> > ---
> > Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> > 
> > Error: Within
> >           Df Sum Sq Mean Sq F value Pr(>F) Residuals 22 9.8250  0.4466
> > > anova(lme(fixed=distance~age,random=~1+age|Subject,data=OrthoFem))
> >             numDF denDF   F-value p-value
> > (Intercept)     1    32 1269.7764  <.0001
> > age             1    32   52.4517  <.0001
> > 
> > This time the F values are (almost) identical, the numerator 
> > degrees of freedom are the same, but the denominator degrees 
> > of freedom are very different (10 for aov vs. 32 for lme). I 
> > understand that there is an issue with the estimation of that 
> > number, but I would naively expect the number given by lme to 
> > be close to that provided by aov is the case of a balanced 
> > dataset. That's obviously not true in the case of an 
> > ANCOVA... But why?? And how should I interpret the F-test 
> > given by anova.lme?
> > 
> > Thanks in advance for your help !
> > --
> > Jean-Baptiste Ferdy
> > Institut des Sciences de l'Évolution de Montpellier CNRS UMR 
> > 5554 Université Montpellier 2
> > 34 095 Montpellier cedex 05
> > tel. +33 (0)4 67 14 42 27
> > fax  +33 (0)4 67 14 36 22
> > 
> > ______________________________________________
> > R-help at 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
> > and provide commented, minimal, self-contained, reproducible code.
> >
> 
> ______________________________________________
> R-help at 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
> and provide commented, minimal, self-contained, reproducible code.
-- 
Manuel A. Morales
http://mutualism.williams.edu