lme versus aov repeated measures variance specification

Question

I am performing a repeated measures ANOVA in R, where the experimental design is composed of a between-subjects factor (Prodotto, three levels), the within-subjects factor is time (Tempo, four levels&mdahs;an ordered factor) and each Product * Time is repeated on eleven patients (idPaziente).

The aov() syntax would be:

myAov<-aov(valore ~ Tempo*Prodotto + Error(idPaziente/Tempo),data=mydata)

whilst the lme() (from the nlme package) syntax would be

myModel <- lme(valore ~  Tempo*Prodotto,data=mydata, random= ~ 1 | idPaziente/Tempo)

In Finch (2014) Multilevel Modeling with R, I find a specification of a longitudinal model that would be written as:

myModel <- lme(valore ~  Tempo*Prodotto,data=mydata, random= ~ 1 | idPaziente)

I guess the equivalent aov version would have Error(idPaziente) instead of Error(idPaziente/Tempo).

I would like to ask:

1. What is the correct specification for my repeated measures model, both for aov() and lme()?

2. Why are (1 | idPaziente/Tempo) and Error(idPaziente/Tempo) more appropriate than (1 | idPaziente) and Error(idPaziente)?

3. How can I perform post hoc analysis taking on the Prodotto term?

Answer 1

This depends on your nesting structure (example dataset):

set.seed(123)
df <- data.frame(
  idPaziente = rep(c(1:12), each = 11), 
  Tempo = rep(c(1:4),33),
  Prodotto = rep(letters[1:3], 11, each = 4),
  valore = rnorm(132, 6, 1)  
)  

Using random = ~ 1 | idPaziente/Tempo:

fit <- lme(valore ~  Tempo * Prodotto, data = df, random = ~ 1 | idPaziente/Tempo)

To get the ANOVA table you could do:

library(car)
Anova(fit)

If your idPaziente and Tempo is explicitly nested,

df$newID <- with(df, paste(idPaziente, Tempo, sep = "-"))

you only need to specify random = ~ 1 | newID

fit2 <- lme(valore ~  Tempo * Prodotto, data = df, random = ~ 1 | newID)

Anova(fit2)

then:

summary(fit) and summary(fit2) will give you the same results. Also have a look at the Number of Groups: statement. This is a good check to see whether your nesting worked out correctly and reflects your data structure.

If Prodotto is significant, you could follow up with:

library(lsmeans)
lsmeans(fit, pairwise ~ Prodotto)
lsmeans(fit2, pairwise ~ Prodotto)

And the same should hold for the aov() function, too.