depr=read.table("http://sites.williams.edu/bklingen/files/2012/02/depression.txt", header=TRUE)
head(depr,10)
tail(depr,10)

## Exploring:
tabledat <- table(depr$severity,depr$drug,depr$time,depr$outcome)
tabledat1 <- as.data.frame(tabledat)
names(tabledat1) <- c("severity", "drug",     "time",     "outcome", "count")
tabledat1
n <- tabledat1$count[1:12]+tabledat1$count[13:24]
tabledat2 <- cbind(tabledat1[13:24,], n, tabledat1$count[13:24]/n)
names(tabledat2)[7] <- "samp_prop"
tabledat2
require(ggplot2)  
ggplot(data=tabledat2, aes(x=time,y=samp_prop, color=drug, group=drug)) + geom_point() + geom_line() + facet_grid(~severity) + theme_bw() 

## Modelling:
## Assuming all observation on the same subject are independent:
fit_logistic <- glm(outcome~severity+drug+time+drug*time,family=binomial,data=depr)
summary(fit_logistic)
anova(fit_logistic, test="LRT")

require(gee) #load GEE package
fit_exch <- gee(outcome~severity+drug+time+drug*time,family=binomial,id=case,data=depr,corstr = "exchangeable")
summary(fit_exch)

fit_unstr <- gee(outcome~severity+drug+time+drug*time,family=binomial,id=case,data=depr,corstr = "unstructured")
summary(fit_unstr)

fit_ind <- gee(outcome~severity+drug+time+drug*time,family=binomial,id=case,data=depr,corstr = "independence")
summary(fit_ind)

#Let's go with with exchangeable model
# Plotting fitted Model
# get design matrix
X <- tabledat2[,1:3]
X <- apply(X,2,as.numeric) #changes factors to numerical values
X
X <- cbind(1, X, X[,1]*X[,2])
coef <- as.matrix(fit_exch$coefficients,ncol=1)
eta <- X %*% coef
fit_prob <- exp(eta)/(1+exp(eta))
tabledat2$fit_prob <- fit_prob
tabledat2
ggplot(data=tabledat2, aes(x=time,y=fit_prob, color=drug, group=drug)) + geom_point() + geom_line() + facet_grid(~severity) + theme_bw()