#########################
# Edps/Psych/Stat 587
# Fall 2020
# c.j.anderson
#
# Depression data from Agresti
#
# glmer is an alternative package to try
#
###########################

library(lme4)
lirbary(lmerTest)
library(gee)

depress <- read.table(file="D:/Dropbox/edps587/lectures/10 logreg/depression_data_long.txt", header=TRUE)

head(depress)

######################################################
# Ignore heirarchical structure
######################################################

# standard logistic regression

model1a <- glm(y ~ severe +  Rx + time + Rx*time, data=depress, family=binomial)
summary(model1a)

# Check: in biomial distirbution var= n*p*(1-p)*phi.  The esimated value of phi=0.99
#                                                     which suggests binomial is OK.

model1b <- glm(y ~ severe +  Rx + time + Rx*time, data=depress, family=quasibinomial)
summary(model1b)

#####################################################
#  GEE 
#####################################################

# exchangable
model.gee <- gee(y ~ severe  + Rx + time + Rx*time, id, data=depress, family=binomial, corstr="exchangeable")
summary(model.gee)

# Unstructured
summary(model.gee <- gee(y ~ severe  + Rx + time + Rx*time, id, data=depress, family=binomial, corstr="unstructured"))

##################################################
#  GLIMM using glmer which is in package lme4
###################################################

# default estimation alogrithm:  laplace

model2 <- glmer(y ~ severe +  Rx + time + Rx*time + (1 |id ), data=depress, family=binomial)
summary(model2)


# estimation alogrithm:  quad   (gold standard)

model3 <- glmer(y ~ severe +  Rx + time + Rx*time + (1 |id ), data=depress, family=binomial, nAGQ=10)
summary(model3)