###################################################
### 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>")
source('latticedefaults.R')


###################################################
### chunk number 2: 1a
###################################################
library(MASS)
library(lattice)
library(gee)
library(nlme)
set.seed(123)
url <- "http://www.statistik.lmu.de/institut/lehrstuhl/semwiso/longitudinal_ss08/"
ex <- read.table(paste(url,"download/amsterdam.txt", sep = ""), header = T)
ex$smoker <- factor(ex$smoker)
print(xyplot(tc ~ t | smoker*gender,groups=id ,data=ex,type="b", col=1))


###################################################
### chunk number 3: 1a2
###################################################
print(bwplot(tc ~ as.factor(t)| gender,groups=id ,data=ex,type="b",horizontal=FALSE))


###################################################
### chunk number 4: 1b
###################################################
t <- unique(ex$t)
E <- with(ex,tapply(tc,list(t,gender),mean))
SD <- with(ex,tapply(tc,list(t,gender),sd))
res <- cbind(E[,1],SD[,1],E[,2],SD[,2])
colnames(res) <- c("E-FEMALE","SD-FEMALE","E-MALE","SD-MALE")
res


###################################################
### chunk number 5: 1c
###################################################
#reshape to wide format
ex.wide <- reshape(ex,timevar="t",idvar="id",direction="wide",v.names=c("tc","bfatness","smoker"))

#within subject correlation
cor <- cor(ex.wide[,paste("tc.",t,sep="")])
print(cor,digits=2)
round(range(cor[row(cor) != col(cor)]),2)


###################################################
### chunk number 6: 1d
###################################################
lcorstr <- list("independence","stat_M_dep","exchangeable","unstructured")
gees <- lapply(lcorstr,function(corstr) {
  print(corstr)
  gee( tc ~ 1 + t + fitness + bfatness + gender + smoker, corstr=corstr,Mv=ifelse(corstr=="stat_M_dep",5,0),data=ex)
})
estimate <- sapply(1:length(gees),function(i) coef(gees[[i]]))
sd <- sapply(1:length(gees),function(i) sqrt(diag(gees[[i]]$robust.variance)))
colnames(estimate) <- lcorstr
colnames(sd) <- lcorstr

print(estimate,digits=2)
print(sd,digits=2)


###################################################
### chunk number 7: 1e
###################################################
m.slope.lme <- lme(tc ~ 1 + t, random= ~ 1+t|id, data=ex)
summary(m.slope.lme)


###################################################
### chunk number 8: hatbeta
###################################################
print(fixef(m.slope.lme),digits=4)


###################################################
### chunk number 9: 1e2
###################################################
getVarCov(m.slope.lme)


###################################################
### chunk number 10: 1f
###################################################
m.noslope.lme <- update(m.slope.lme, random=~1|id)
anova(m.noslope.lme,m.slope.lme)


###################################################
### chunk number 11: 1g eval=FALSE
###################################################
## cor.pos <- matrix(c(1,0.9,0.9,1),2,2)
## cor.neg <- matrix(c(1,-0.9,-0.9,1),2,2)
## library(mvtnorm)
## 
## do.plot <- function(corM,main) {
##   res <- mvrnorm(n=10,mu=c(5,2),Sigma=corM)
##   plot(NA,xlim=c(0,5),ylim=c(3,10),type="n",xlab="t",ylab="",main=main)
##   abline(coef=c(5,2),lty=1,lwd=4,col='darkgrey')
##   for (i in 1:nrow(res)) abline(coef=res[i,])
## }
## 
## set.seed(123)
## par(mfcol=c(1,2))
## do.plot(cor.pos,"Positive Korrelation")
## do.plot(cor.neg,"Negative Korrelation")


###################################################
### chunk number 12: 1g2
###################################################
cor.pos <- matrix(c(1,0.9,0.9,1),2,2)
cor.neg <- matrix(c(1,-0.9,-0.9,1),2,2)
library(mvtnorm)

do.plot <- function(corM,main) {
  res <- mvrnorm(n=10,mu=c(5,2),Sigma=corM)
  plot(NA,xlim=c(0,5),ylim=c(3,10),type="n",xlab="t",ylab="",main=main)
  abline(coef=c(5,2),lty=1,lwd=4,col='darkgrey')
  for (i in 1:nrow(res)) abline(coef=res[i,])
}
set.seed(123)
par(mfcol=c(1,2))
do.plot(cor.pos,"Positive Korrelation")
do.plot(cor.neg,"Negative Korrelation")
par(mfcol=c(1,1))


###################################################
### chunk number 13: 1h
###################################################
m.lme <- lme( tc ~ 1 + t, random=~1 | id,data=ex,method="ML")
m.gee <- gee( tc ~ 1 + t, id=id, corstr="exchangeable",data=ex)
fixed.effects(m.lme)
vcov(m.lme)
m.gee$robust.variance
m.gee$naive.variance


###################################################
### chunk number 14: 2
###################################################
t <- unique(ex$t)
q75 <- sapply(t,function(time) { quantile(subset(ex, t==time)$tc,0.75) })
tc.dich <- (ex$tc > rep(q75,length(unique(ex$id))))*1
ex.dich <- cbind(ex,tc.dich)


###################################################
### chunk number 15: 2abc
###################################################
set.seed(123)
MCAR <- MAR <- NMAR <- ex.dich

for(i in unique(ex$id)){
            for(j in unique(ex$t)){
              #missing completely at random with P(missing)=1/3
              if(runif(1)<.3333) MCAR[(ex$id==i & ex$t ==j),'tc'] <- NA

              #missing at random:  subjects drop out AFTER they enter risk group
              if(MAR[(ex$id==i & ex$t == j),'tc.dich']==1) MAR[(ex$id==i & ex$t > j),'tc'] <- NA

              #missing not at random:  subjects drop out with P=.8 BEFORE they enter risk group
              if((NMAR[(ex$id==i & ex$t == j),'tc.dich']==1) & (runif(1) <.8)) NMAR[(ex$id==i & ex$t >= j),'tc'] <- NA
   }
}
MCAR <- groupedData(tc ~ t | id, MCAR)
MAR <- groupedData(tc ~ t | id, MAR)
NMAR <- groupedData(tc ~ t | id, NMAR)

m.noslope.lme
mMCAR <-  update(m.noslope.lme, data=MCAR, na.action=na.exclude)
mMAR <-  update(m.noslope.lme, data=MAR, na.action=na.exclude)
mNMAR <-  update(m.noslope.lme, data=NMAR, na.action=na.exclude)

compareEstimates <-
cbind(NoMiss=c(fixef(m.noslope.lme),as.numeric(VarCorr(m.noslope.lme)[,2])),
      MCAR  =c(fixef(mMCAR),as.numeric(VarCorr(mMCAR)[,2])),
                        MAR  =c(fixef(mMAR),as.numeric(VarCorr(mMAR)[,2])),
                        NMAR  =c(fixef(mNMAR),as.numeric(VarCorr(mNMAR)[,2])))
rownames(compareEstimates) <- c('Intercept','t','sd(RandomIntercept)','sd(Error)')
compareEstimates


###################################################
### chunk number 16: 3a
###################################################
library(glmmML)
#Default values
m.glme.1 <- glmmML(tc.dich ~ 1 + t + smoker + gender, cluster=id, family=binomial,data=ex.dich)
#Number of quadrature points makes a big difference! Laplace approximation might not be very good.
m.glme.2 <- glmmML(tc.dich ~ 1 + t + smoker + gender, cluster=id, data=ex.dich,method="ghq",n.points=25)
round((coef(m.glme.2) - coef(m.glme.1))/coef(m.glme.2),4)
summary(m.glme.1)
summary(m.glme.2)

#Alternative using LME4
detach(package:nlme)
library(lme4)
m.glme.3 <- lmer(tc.dich ~ 1 + t + smoker + gender + (1|id), family=binomial, data=ex.dich,method="Laplace")
m.glme.4 <- lmer(tc.dich ~ 1 + t + smoker + gender + (1|id), family=binomial, data=ex.dich,method="PQL")
summary(m.glme.3)
summary(m.glme.4)
m.glme.1$deviance
m.glme.2$deviance
deviance(m.glme.3)
deviance(m.glme.4)


###################################################
### chunk number 17: 3a2
###################################################
#This should be same
round((fixef(m.glme.3) - coef(m.glme.1))/fixef(m.glme.3),4)
#This sucks
round((fixef(m.glme.3) - coef(m.glme.2))/coef(m.glme.2),4)
#This sucks worse
round((fixef(m.glme.4) - coef(m.glme.2))/coef(m.glme.2),4)


###################################################
### chunk number 18: 
###################################################
coef(m.glme.2)


###################################################
### chunk number 19: 3b2
###################################################
(smokerCI.glmm<-exp(coef(m.glme.2)[3] + c(est=0,lo=-1,hi=1)*qnorm(0.975)* m.glme.2[["coef.sd"]][3]))
(maleCI.glmm <- exp(coef(m.glme.2)[4] + c(est=0,lo=-1,hi=1)*qnorm(0.975)* m.glme.2[["coef.sd"]][4]))


###################################################
### chunk number 20: 3c
###################################################
library(gee)
m.gee <- gee(tc.dich ~ 1 + t + smoker + gender, id = id, data = ex, family = binomial, corstr = "unstructured")
summary(m.gee)
(smokerCI.gee<-exp(coef(m.gee)[3] + c(est=0,lo=-1,hi=1)*qnorm(0.975)* sqrt(m.gee$robust.variance[3,3])))
(maleCI.gee <- exp(coef(m.gee)[4] + c(est=0,lo=-1,hi=1)*qnorm(0.975)* sqrt(m.gee$robust.variance[4,4])))