[R] degrees of freedom in lme

[Jean-Baptiste Ferdy]

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