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


###################################################
### chunk number 2: cd4
###################################################
url<-'http://www.statistik.lmu.de/institut/lehrstuhl/semwiso/longitudinal_ss08/'
cd4<-read.table(paste(url,'download/cd4_Diggle.txt',sep=''),header=T)
cd4$drug<-factor(cd4$drug,labels=c('yes','no'))
cd4$ID <- factor(cd4$ID)


###################################################
### chunk number 3: cd4a1
###################################################
library('lattice')
xyplot(CD4~Time,data=cd4, group=ID,
    panel=function(x,y){
      panel.superpose(x, y,..., col="grey",type='l',lwd=.1)
      panel.xyplot(x, y, col="black",type='p', pch=19, cex=.3)
      panel.loess(x, y, lwd=2.5, col="blue")
      panel.loess(x, y, lwd=2.5, span=1, degree=2, col="green")
      panel.loess(x, y, lwd=2.5, span=0.1, evaluation=100, col="red")
    })


###################################################
### chunk number 4: cd4afig
###################################################
library('lattice')
print(xyplot(CD4~Time,data=cd4, group=ID,
    panel=function(x,y,...){
      panel.superpose(x, y,..., col="grey",type='l',lwd=.1)
      panel.xyplot(x, y, col="black",type='p', pch=19, cex=.3)
      panel.loess(x, y, lwd=2.5, col="blue")
      panel.loess(x, y, lwd=2.5, span=1, degree=2, col="green")
      panel.loess(x, y, lwd=2.5, span=0.1, evaluation=100, col="red")
    }))


###################################################
### chunk number 5: cd4b1 eval=FALSE
###################################################
## plot(table(round(4*cd4$Time)/4))
## timepoints<-c(floor(min(cd4$Time)),-1.25,-.75,-.25,0,.25,.75,1.25,1.75,2.25,2.75,3.25,ceiling(max(cd4$Time)))
## cd4$Time2 <- cut(cd4$Time, breaks=timepoints)
## (table(cd4$Time2))


###################################################
### chunk number 6: cd4b1fig
###################################################
plot(table(round(4*cd4$Time)/4))
timepoints<-c(floor(min(cd4$Time)),-1.25,-.75,-.25,0,.25,.75,1.25,1.75,2.25,2.75,3.25,ceiling(max(cd4$Time)))
cd4$Time2 <- cut(cd4$Time, breaks=timepoints)
(table(cd4$Time2))


###################################################
### chunk number 7: cd4c1
###################################################
m <- lm(CD4 ~ Time , data=cd4)
cd4$eps <- resid(m)
o <- order(cd4$Time)
cd4.kat <- reshape(cd4[o,],v.names=c('CD4','eps'),timevar='Time2',idvar='ID',direction='wide',drop=c("Age",'Time',"cigpacks","drug","sexpartners","cesd"))
colnames(cd4.kat)[seq(2,24,by=2)] <- levels(cd4$Time2)
colnames(cd4.kat)[seq(3,25,by=2)] <- paste('e.',levels(cd4$Time2),sep='')


###################################################
### chunk number 8: cd4c1s
###################################################
cd4[cd4$ID=='10005',]
cd4.kat[cd4.kat$ID=='10005',]


###################################################
### chunk number 9: cd4c2
###################################################
trellis.par.set(list(add.text = list(cex = .6)))
print(splom(cd4.kat[,seq(2,24,by=2)],pscale=0,cex=.5))
round(cor(cd4.kat[,seq(2,24,by=2)],use='pairwise'),2)


###################################################
### chunk number 10: cd4c3
###################################################
trellis.par.set(list(add.text = list(cex = .5)))
print(splom(cd4.kat[,seq(3,25,by=2)],pscale=0,cex=.5))
round(cor(cd4.kat[,seq(3,25,by=2)],use='pairwise'),2)


###################################################
### chunk number 11: cd4c4
###################################################
int.means <- with(cd4, tapply(Time, Time2, mean) )
boxplot(CD4~Time2, data=cd4, at = int.means, varwidth=T, xlim=c(-2.1,4.2), boxwex=.3, axis.cex=.5, main='CD4')


###################################################
### chunk number 12: cd4c5
###################################################
var.time2<-with(cd4, tapply(CD4, Time2, var))
plot(int.means, sqrt(var.time2), type='l', main='sd(CD4)', xlab='Time' , ylab='sd(CD4)')