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

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

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.