###################################################
### chunk number 1: prelim
###################################################
setwd('W:/lehre/longitudinal_sommer08/Übung')
options(width = 70, prompt = "R>")


###################################################
### chunk number 2: 1a
###################################################
rats <- read.table( 'http://www.statistik.lmu.de/institut/lehrstuhl/semwiso/longitudinal_ss08/download/rats.txt', header=T,na.strings='.')
rats$GROUP <- factor(rats$GROUP,labels=c('low','high','control'))
rats$SUBJECT <- factor(rats$SUBJECT)


###################################################
### chunk number 3: 1b
###################################################
rats.wide <- reshape(data=rats, v.names='RESPONSE', idvar='SUBJECT', timevar='TIME',direction='wide')

#fehlende Werte pro Tier
rats.wide$NAs <- rowSums(is.na(rats.wide))

#Tabelle
(na.structure <- xtabs(~NAs+GROUP,rats.wide))
barplot(na.structure, legend.text = rownames(na.structure), main = 'Anzahl fehlender Beobachtungen')


###################################################
### chunk number 4: 1c
###################################################
library(nlme)
keep <- with(rats.wide, SUBJECT[NAs<=4])
rats2 <- groupedData(RESPONSE ~ TIME|GROUP/SUBJECT,
                     data=rats[rats$SUBJECT %in% keep,],
                     labels = list(x = "Age", y = "Distance"), units = list(x = "[Days]", y = "[Pixels]"))
print(plot(rats2, display = 'GROUP', inner = ~SUBJECT))


###################################################
### chunk number 5: 1c2
###################################################
rats$SUBJECT <- factor(rats$SUBJECT)
rats.complete <- rats[complete.cases(rats),]
lm1 <- lm(RESPONSE ~ TIME, data=rats.complete)
xyplot(resid(lm1)~ rats.complete$TIME|rats.complete$SUBJECT,
    panel=function(x,y){panel.points(x,y);panel.lines(x,y);panel.lines(x=x,y=0,col='black')})


###################################################
### chunk number 6: 1c2fig
###################################################
print(xyplot(resid(lm1)~ rats.complete$TIME|rats.complete$SUBJECT,
    panel=function(x,y){panel.points(x,y);panel.lines(x,y);panel.lines(x=x,y=0,col='black')}))



###################################################
### chunk number 7: 1dlin1
###################################################
lmlist1<-lmList(RESPONSE ~I(TIME-45), rats2, na.action=na.omit)

plot(fitted(lmlist1),resid(lmlist1))
abline(h=0,col='grey')
o <- order(fitted(lmlist1))
lines(lowess(fitted(lmlist1)[o],resid(lmlist1)[o]),col=2)

complete <- complete.cases(rats2)
xyplot(resid(lmlist1)~ rats2$TIME[complete]|rats2$SUBJECT[complete],
    panel=function(x,y){panel.points(x,y);panel.lines(x,y);panel.lines(x=x,y=0,col='black')})


###################################################
### chunk number 8: 1dlin2
###################################################
plot(fitted(lmlist1),resid(lmlist1))
abline(h=0,col='grey')
o <- order(fitted(lmlist1))
lines(lowess(fitted(lmlist1)[o],resid(lmlist1)[o]),col=2)


###################################################
### chunk number 9: 1dlin3
###################################################
complete <- complete.cases(rats2)
print(xyplot(resid(lmlist1)~ rats2$TIME[complete]|rats2$SUBJECT[complete],
    panel=function(x,y){panel.points(x,y);panel.lines(x,y);panel.lines(x=x,y=0,col='black')}))


###################################################
### chunk number 10: 1e1 eval=FALSE
###################################################
## #Visualisierung Parameter
## plot(intervals(lmList-Objekt))


###################################################
### chunk number 11: 1d1fig
###################################################
rats2$logT= log(1+(rats2$TIME-45)/10)
lmlist2<-lmList(RESPONSE ~ logT, rats2, na.action=na.omit)
plot(fitted(lmlist2),resid(lmlist2))
abline(h=0,col='grey')
o <- order(fitted(lmlist2))
lines(lowess(fitted(lmlist2)[o],resid(lmlist2)[o]),col=2)


###################################################
### chunk number 12: 1d1fig2
###################################################
print(plot(intervals(lmlist2)))


###################################################
### chunk number 13: 1e
###################################################
get.R2adj<-function(x)
{
    summary(x)$adj.r.squared
}

R2adj.lin <- unlist(lapply(lmlist2,get.R2adj))

lmlist3<-lmList(RESPONSE ~ logT + I(logT^2), rats2, na.action=na.omit)
plot(fitted(lmlist3),resid(lmlist3))
abline(h=0,col='grey')
o <- order(fitted(lmlist3))
lines(lowess(fitted(lmlist3)[o],resid(lmlist3)[o]),col=2)


###################################################
### chunk number 14: 1e21
###################################################
print(plot(intervals(lmlist3)))


###################################################
### chunk number 15: 1e22
###################################################
R2adj.qu <- unlist(lapply(lmlist3,get.R2adj))
plot(R2adj.lin,R2adj.qu)
abline(c(0,1))