[R] Result depends on order of factors in unbalanced designs (lme, anova)?

[Karl Knoblick]

Dear R-Community!

For example I have a study with 4 treatment groups (10 subjects per group) and 4 visits. Additionally, the gender is taken into account. I think - and hope this is a goog idea (!) - this data can be analysed using lme as below.

In a balanced design everything is fine, but in an unbalanced design there are differences depending on fitting y~visit*treat*gender or y~gender*visit*treat - at least with anova (see example). Does this make sense? Which ordering might be the correct one?

Here the example script:
library(nlme)
set.seed(123)
# Random generation of data:
NSubj<-40 # No. of subjects
set.seed(1234)
id<-factor(rep(c(1:NSubj),4)) # ID of subjects
treat<-factor(rep(rep(1:4,each=5),4)) # Treatment 4 Levels
gender<-factor(rep(rep(1:2, each=20),4))
visit<-factor(rep(1:4, each=NSubj))
y<-runif(4*NSubj) # Results
# Add effects
y<-y+0.01*as.integer(visit)
y<-y+0.02*as.integer(gender)
y<-y+0.024*as.integer(treat)
df<-data.frame(id, treat, gender, visit, y)
# groupedData object for lme
gdat<-groupedData(y ~ visit|id, data=df)
# fits - different ordering of factors
fit1<-lme(y ~ visit*treat*gender, data=gdat, random = ~visit|id)
anova(fit1)
fit2<-lme(y ~ gender*treat*visit, data=gdat, random = ~visit|id)
anova(fit2)
# Result: identical (balanced design so far), ok
# Now change gender of subject 1
gdat$gender[c(1,41,81,121)]<-2
# onece more fits with different ordering of factors
fit1<-lme(y ~ visit*treat*gender, data=gdat, random = ~visit|id)
anova(fit1)
fit2<-lme(y ~ gender*treat*visit, data=gdat, random = ~visit|id)
anova(fit2)
# Result: There are differences!!

Hope anybody can help or give me advice how to interpret these results correctly or how to avoid this problem! Is there a better possibility to analyse these data than lme?

Thanks!
Karl