###################################################
### chunk number 1: prelim
###################################################
rm(list=ls())
try(detach(package:nlme))
try(detach(package:lme4))
setwd('W:/lehre/longitudinal_sommer08/uebung')
options(width = 80, prompt = "R>")


###################################################
### chunk number 2: 
###################################################
library('nlme')
url <- 'http://www.statistik.lmu.de/institut/lehrstuhl/semwiso/longitudinal_ss08'
rats <- read.table( paste(url,'/download/rats.txt',sep=''), header=T,na.strings='.')

rats$GROUP <- factor(rats$GROUP,labels=c('low','high','control'))
rats$SUBJECT <- factor(rats$SUBJECT)
rats$logT= log(1+(rats$TIME-45)/10)

rats <- rats[complete.cases(rats),]
obs <- table(rats$SUBJECT)
zuvieleNAs <- levels(rats$SUBJECT)[obs < 3]
rats <- rats[!(rats$SUBJECT %in% zuvieleNAs),]


###################################################
### chunk number 3: lmlistlin
###################################################
rats <- groupedData(RESPONSE ~ logT|SUBJECT,
                     data=rats,
                     labels = list(x = "log(Age)", y = "Distance"), units = list(x = " ", y = "[Pixels]"))
lmlist.lin <- lmList(RESPONSE ~ logT, rats)
print(plot(intervals(lmlist.lin)))


###################################################
### chunk number 4: lmm1
###################################################
lmm.lin1 <- lme(RESPONSE ~ logT, random = ~1|SUBJECT, data=rats)
summary(lmm.lin1)


###################################################
### chunk number 5: lmm1sol1
###################################################
(getVarCov(lmm.lin1))
(d.sq <- getVarCov(lmm.lin1)[1])
(sigma.sq <- (lmm.lin1$sigma)^2)
(cor.withinGroup <- d.sq/(d.sq + sigma.sq))


###################################################
### chunk number 6: lmm2
###################################################
lmm.lin2 <- lme(RESPONSE ~ logT, random = ~ 1 + logT|SUBJECT, data=rats)
summary(lmm.lin2)


###################################################
### chunk number 7: lmm2sol1
###################################################
(D<-getVarCov(lmm.lin2))
cov2cor(D)


###################################################
### chunk number 8: lmm2sol2
###################################################
lmm.lin3 <- lme(RESPONSE ~ logT, random = pdDiag(~ logT), data=rats)
lmm.lin3 <- lme(RESPONSE ~ logT, random = pdDiag(~ 1 + logT), data=rats)
lmm.lin3 <- lme(RESPONSE ~ logT, random = list(SUBJECT=pdDiag(~1+ logT)), data=rats)
lmm.lin3 <- lme(RESPONSE ~ logT, random = list(SUBJECT=pdDiag(~ logT)), data=rats)
(summary(lmm.lin3))
(getVarCov(lmm.lin3))


###################################################
### chunk number 9: lmm2sol3a
###################################################
print(plot(compareFits(coef(lmlist.lin), coef(lmm.lin3))))


###################################################
### chunk number 10: lmm2sol3b
###################################################
print(pairs(compareFits(coef(lmlist.lin), coef(lmm.lin3))))


###################################################
### chunk number 11: lmm2sol4
###################################################
print(plot(comparePred(lmlist.lin, lmm.lin2)))


###################################################
### chunk number 12: lmmqu1
###################################################
lmlist.qu <- lmList(RESPONSE ~ poly(logT,2), rats)
try(lmm.qu1 <- lme(RESPONSE ~ poly(logT,2), random =~ poly(logT,2), data=rats))
detach(package:nlme)
library(lme4)
lmm.qu1 <- lmer(RESPONSE ~ logT + I(logT^2) + (logT + I(logT^2)|SUBJECT), data=rats)
summary(lmm.qu1)
lmm.qu2 <- lmer(RESPONSE ~ logT + I(logT^2) + (1|SUBJECT) + (-1 + logT + I(logT^2)|SUBJECT), data=rats)
summary(lmm.qu2)
(D1 <- VarCorr(lmm.qu1))
cov2cor(D1[[1]])
(D2 <- VarCorr(lmm.qu2))
cov2cor(D2[[2]])


###################################################
### chunk number 13: lmmqu2fig
###################################################
library(nlme)
lmerlmList.comparePred<-function(lmm,lmlist,primary,level=0,length.out=51)
{
 args <- list(object = lmlist,
                 primary = primary,
                 level = level,
                 length.out = length.out)
 val1 <- do.call("augPred", args)
 val1 <-val1[val1$.type=='original',]
 temp1 <- cbind(val1[,1:2],fitted(lmlist),.type=rep(deparse(substitute(lmlist)),nrow(val1)))
 temp2 <- cbind(val1[,1:2],fitted(lmm),.type=rep(deparse(substitute(lmm)),nrow(val1)))
 names(temp1) <- names(temp2) <- names(val1)
 val <- rbind(temp2,temp1,val1)
 #names(val)[2] <- '.groups'
 class(val) <- c("comparePred", "augPred", class(val))
 attr(val, "labels") <- attr(val1, "labels")
 attr(val, "units") <- attr(val1, "units")
 attr(val, "formula") <- attr(val1, "formula")
 val
}
print(plot(lmerlmList.comparePred(lmm.qu2,lmlist.qu,primary='logT')))