What is the difference between using aov() and lme() in analyzing a longitudinal dataset?

Question

Can anyone tell me the difference between using aov() and lme() for analyzing longitudinal data and how to interpret results from these two methods?

Below, I analyze the same dataset using aov() and lme() and got 2 different results. With aov(), I got a significant result in the time by treatment interaction, but fitting a linear mixed model, time by treatment interaction is insignificant.

> UOP.kg.aov <- aov(UOP.kg~time*treat+Error(id), raw3.42)
> summary(UOP.kg.aov)

Error: id
          Df  Sum Sq Mean Sq F value Pr(>F)
treat      1   0.142  0.1421  0.0377 0.8471
Residuals 39 147.129  3.7725               

Error: Within
            Df  Sum Sq Mean Sq  F value  Pr(>F)    
time         1 194.087 194.087 534.3542 < 2e-16 ***
time:treat   1   2.077   2.077   5.7197 0.01792 *  
Residuals  162  58.841   0.363                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> UOP.kg.lme <- lme(UOP.kg~time*treat, random=list(id=pdDiag(~time)), 
                    na.action=na.omit, raw3.42)
> summary(UOP.kg.lme)
Linear mixed-effects model fit by REML
 Data: raw3.42 
       AIC      BIC    logLik
  225.7806 248.9037 -105.8903

Random effects:
 Formula: ~time | id
 Structure: Diagonal
        (Intercept)      time  Residual
StdDev:   0.6817425 0.5121545 0.1780466

Fixed effects: UOP.kg ~ time + treat + time:treat 
                 Value Std.Error  DF   t-value p-value
(Intercept)  0.5901420 0.1480515 162  3.986059  0.0001
time         0.8623864 0.1104533 162  7.807701  0.0000
treat       -0.2144487 0.2174843  39 -0.986042  0.3302
time:treat   0.1979578 0.1622534 162  1.220053  0.2242
 Correlation: 
           (Intr) time   treat 
time       -0.023              
treat      -0.681  0.016       
time:treat  0.016 -0.681 -0.023

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-3.198315285 -0.384858426  0.002705899  0.404637305  2.049705655 

Number of Observations: 205
Number of Groups: 41 

Answer 1

Based on your description, it appears that you have a repeated-measures model with a single treatment factor. Since I do not have access to the dataset (raw3.42), I will use the Orthodont data from the nlme package to illustrate what is going on here. The data structure is the same (repeated measurements for two different groups - in this case, males and females).

If you run the following code:

library(nlme)
data(Orthodont)

res <- lme(distance ~ age*Sex, random = ~ 1 | Subject, data = Orthodont)
anova(res)

you will get the following results:

            numDF denDF  F-value p-value
(Intercept)     1    79 4123.156  <.0001
age             1    79  122.450  <.0001
Sex             1    25    9.292  0.0054
age:Sex         1    79    6.303  0.0141

If you run:

res <- aov(distance ~ age*Sex + Error(Subject), data = Orthodont)
summary(res)

you will get:

Error: Subject
          Df Sum Sq Mean Sq F value   Pr(>F)   
Sex        1 140.46 140.465  9.2921 0.005375 **
Residuals 25 377.91  15.117                    

Error: Within
          Df  Sum Sq Mean Sq  F value  Pr(>F)    
age        1 235.356 235.356 122.4502 < 2e-16 ***
age:Sex    1  12.114  12.114   6.3027 0.01410 *  
Residuals 79 151.842   1.922                     

Note that the F-tests are exactly identical.

For lme(), you used list(id=pdDiag(~time)), which not only adds a random intercept to the model, but also a random slope. Moreover, by using pdDiag, you are setting the correlation between the random intercept and slope to zero. This is a different model than what you specified via aov() and hence you get different results.

Answer 2

I would just add that you might wish to install the car package and use Anova() that this package provides instead of anova() because for aov() and lm() objects, the vanilla anova() uses a sequential sum of squares, which gives the wrong result for unequal sample sizes while for lme() it uses either the type-I or the type-III sum of squares depending on the type argument, but the type-III sum of squares violates marginality-- i.e. it treats interactions no differently than main effects.

The R-help list has nothing good to say about type-I and type-III sums of squares, and yet these are the only options! Go figure.

Edit: Actually, it looks like type-II is not valid if there is a significant interaction term, and it seems the best anybody can do is use type-III when there are interactions. I got tipped off to it by an answer to one of my own questions that in turn pointed me to this post.

Answer 3

It appears to me you have multiple measures for each id at each time. You need to aggregate these for the aov because it unfairly increases power in that analysis. I'm not saying doing the aggregate will make the outcomes the same but it should make them more similar.

dat.agg <- aggregate(UOP.kg ~ time + treat + id, raw3.42, mean)

Then run your aov model as before replacing the data with dat.agg.

Also, I believe anova(lme) is more what you want to do to compare the outcomes. The direction and magnitude of an effect is not the same as the ratio of model variance to error.

(BTW, if you do the lme analysis on the aggregate data, which you shouldn't, and check anova(lme) you'll get almost the same results as aov)