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

-
Thanks @Wolfgang; your explanation helps a lot. A follow up question that I have then is this. I am indeed analyzing a repeated-measures model with a single treatment factor. Each subject is randomly assigned to either treatment A or B. Then they are measured at 0 mins, 15 mins, 30 mins, 60 mins, 120 mins and 180 mins. From my understanding, time should be a random factor because it is just samples from time 0 to 180 mins. So, should I do: lme(UOP.kg~time*treat, random=~time|id, raw3.42)? – biostat_newbie Feb 15 '11 at 1:59
Yes, but I would think of it this way: You are essentially allowing the intercept and slope of the regression line (of UOP.kg on time) to differ (randomly) between subjects within the same treatment group. This is what random=~time|id will do. What the model will then tell you is the estimated amount of variability in the intercepts and the slopes. Moreover, the time:treat interaction term indicates whether the average slope is different for A and B. – Wolfgang Feb 15 '11 at 10:13

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.

-

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)

-