M2<-lme(score ~ disease* time* week, random=~1|treatment/Id, method = "REML",
        data = Dat1)
M3<-gls(score ~ disease * time * week, method = "REML", data = Dat1)

df         AIC
M2  21 -2662.715
M3  19 -2612.308

Model   df    AIC        BIC      logLik     Test    L.Ratio    p-value
M2     1  21  -2662.715  -2550.233  1352.358
M3     2  19  -2612.308  -2510.538  1325.154   1 vs 2  54.40752   <.0001