###################################################
### 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: 1data eval=FALSE
###################################################
## url <- "http://www.statistik.lmu.de/institut/lehrstuhl/semwiso/longitudinal_ss08/"
## leprosy <- read.table(paste(url,"download/leprosylong.txt", sep = ""), header = T)
## leprosy$Drug <- factor(leprosy$Drug, levels = c("C", "A", "B"))


###################################################
### chunk number 3: 1data2
###################################################
leprosy <- read.table('leprosy_long.txt')
leprosy$Drug <- factor(leprosy$Drug, levels = c("C", "A", "B"))


###################################################
### chunk number 4: 1a
###################################################
print(xyplot(count ~ time | Drug, groups = id, data = leprosy, type = "b",  layout = c(3, 1)))
(mean.Pre <- with(subset(leprosy, time=='0'),
            tapply(count,Drug,mean)))
(var.Pre <- with(subset(leprosy, time=='0'),
            tapply(count,Drug,var)))
(mean.Post <- with(subset(leprosy, time=='1'),
            tapply(count,Drug,mean)))
(var.Post <- with(subset(leprosy, time=='1'),
            tapply(count,Drug,var)))


###################################################
### chunk number 5: 1b1
###################################################
library(gee)
leprosy.gee.1 <- gee(count ~ time + time:Drug, data = leprosy, family = poisson(link = "log"),  id = id, corstr = "independence")
summary(leprosy.gee.1)


###################################################
### chunk number 6: 1b2
###################################################
summary(glm(formula = count ~ time + time:Drug, family = quasipoisson(link = "log"), data = leprosy))


###################################################
### chunk number 7: 1c1
###################################################
leprosy.gee.2 <- gee(count ~ time + time:Drug, data = leprosy, family = poisson(link = "log"), id = id, corstr = "unstructured")
summary(leprosy.gee.2)


###################################################
### chunk number 8: 1c2
###################################################
summary(leprosy.gee.1)$coefficients
summary(leprosy.gee.2)$coefficients


###################################################
### chunk number 9: 1d1
###################################################
beta <- as.matrix(coef(leprosy.gee.2))
L1 <- matrix(c(0, 0, 1, 0, 0, 0, 0, 1), ncol = 4, byrow = TRUE)
Sigma.star <- leprosy.gee.2$robust.variance
Sigma1 <- L1 %*% Sigma.star %*% t(L1)
(W1 <- t(L1 %*% beta) %*% solve(Sigma1) %*% (L1 %*% beta))
1 - pchisq(W1, 2)


###################################################
### chunk number 10: 1e1
###################################################
L2 <- matrix(c(0, 0, 1, -1), ncol = 4, byrow = TRUE)
Sigma2 <- L2 %*% Sigma.star %*% t(L2)
(W2 <- t(L2 %*% beta) %*% solve(Sigma2) %*% (L2 %*% beta))
1 - pchisq(W2, 1)


###################################################
### chunk number 11: 1f1
###################################################
leprosy$trt <- factor(ifelse(leprosy$Drug == "C", "Placebo", "Treatment"))
leprosy.gee.3 <- gee(count ~ time + time:trt, data = leprosy, family = poisson(link = "log"), id = id, corstr = "unstructured")
summary(leprosy.gee.3)


###################################################
### chunk number 12: 1f2
###################################################
exp(coef(leprosy.gee.3)[3] + c(0,-1, 1) * qnorm(0.975) * sqrt(leprosy.gee.3$robust.variance[3,3]))