##################################################################
# Edps/Psych/Stat 587
# Fall 2020
# c.j.anderson
#
# Data: Cool kid R-code (minus the data)
#       Illistrates how to fit 3-level model 
#       (kids within peer groups within classrooms)
#
# Need package:
#      lme4
#      GEE
#
#  Note shown, but glmer package would be good to look at
#
#  Notes about data:
#    id = peer group
#    classid = classroom id
#    black1 = student 
#    black3 = classroom
#
#names(cool)
# [1] "id"         "classid"    "zagg"       "zpop"       "tough"      "ideal"     
# [7] "aggressive" "average"    "ta"         "ma"         "l"          "md"        
#[13] "dagg"       "dpop"       "black1"     "gagg"       "gnom"       "sex"       
#[19] "axs"        "nxs"        "axn"        "gppro"      "gpprohc"    "site"      
#[25] "urban"      "suburb"     "smallt"     "rural"      "w"          "black3"    
#[31] "class_agg"  "majority"   "n"
##################################################################

library(lme4)
library(lmerTest)
library(gee)       # General Estimating Equations
library(MASS)      # For confidence limits for parameters
library(texreg)

#laptop:
cool <- read.table(file="D:/Dropbox/edps587/lectures/10 logreg/cool_data.txt", header=TRUE)

names(cool)
head(cool,25)
	
# Rename 'model' student to 'ideal' student
names(cool)[names(cool) == "model"] <- "ideal"

###################################################################
# Ignoring clustering                                             #
###################################################################
summary(single <- glm((ideal/n ~ zagg + black1 + sex, data=cool,family = binomial)) 

###################################################################
# Null multilevel logistic regression                             #
###################################################################
# MLE, which is what is in lecture notes

null.mle <- glmer(ideal/n ~ 1 + (1 | classid),data=cool, weights=n, family=binomial,nAGQ=10)
round(confint(null.mle),digits=2)
round(exp(fixef(null.mle)),digits=2)
tau <- 3.011 
(icc <- tau/(tau + pi**2/3))


# LaPlace
null.lp <- glmer(ideal/n ~ 1 + (1 | classid), data=cool, weights=n, family=binomial)
round(confint(null.lp),digits=2)
round(exp(fixef(null.lp)),digits=2)


# MLE, which is what is in lecture notes

summary(model.1 <- glmer(ideal/n ~ zpop + sex + black1 + (1 | classid),data=cool, weights=n, family=binomial,nAGQ=10))
round(confint(model.1),digits=2)
round(exp(fixef(model.1)),digits=2)



# laplace, should pick mle or laplace...for higher level laplace behaves better

summary(model.1b <- glmer(ideal/n ~ zpop + sex + black3 + (1 | classid),data=cool, weights=n, family=binomial))
round(confint(model.1b),digits=2)
round(exp(fixef(model.1b)),digits=2)
round(1/exp(fixef(model.1b)),digits=2)




# MLE, add class-agression

summary(model.2<- glmer(ideal/n ~ zpop + sex + black1 +  class_agg + (1 | classid),data=cool, weights=n, family=binomial,nAGQ=10))
round(confint(model.2),digits=2)
round(exp(fixef(model.2)),digits=2)
round(1/exp(fixef(model.2)),digits=2)

# MLE with random sex --- fails
summary(model.3<- glmer(ideal/n ~ zpop + sex + black1 + class_agg + (1 +sex | classid),data=cool, weights=n, family=binomial,nAGQ=10))


# laplace with random sex --- converges fine
summary(model.3<- glmer(ideal/n ~ zpop + sex + black1 + class_agg+ (1 +sex | classid),data=cool, weights=n, family=binomial))


# laplace with random sex & drop dpop
summary(model.3<- glmer(ideal/n ~ sex + black1 + class_agg + (1 +sex | classid),data=cool, weights=n, family=binomial))


###################################################
# 3- level models 
###################################################
# 3-level empty 
# 3 -level model laplace with only student level fixed effects
summary(model.3a<- glmer(ideal/n ~ 1 + (1 | classid) + (1 |id) ,data=cool, weights=n, family=binomial))

# 3 level student level variables
summary(model.3b<- glmer(ideal/n ~ zpop + zagg  + black1  +(1 | classid) + (1 |id) ,data=cool, weights=n, family=binomial))


# Note sex is actually a level 2 variable
table(cool$id,cool$sex)

# 3 -level model laplace adding level 2 (peer group) fixed effects
summary(model.3c<- glmer(ideal/n ~ zpop + zagg  + black1 + gnom + gagg + sex +(1 | classid) + (1 |id) ,data=cool, weights=n, family=binomial))


# 3 -level model laplace adding classroom variables
# ***** Does not hit convergence criterion ****
summary(model.3d<- glmer(ideal/n ~ zpop + zagg  + black1 + gnom + gagg + sex +
                        + class_agg + majority +(1 | classid) + (1 |id), data=cool, weights=n, family=binomial))

						

# 3 -level model with cross-level interaction 
# ***** Does no hit convergence criterion *****
summary(model.3e<- glmer(ideal/n ~ zpop + zagg + black1 + gnom + gagg + sex +
                        + class_agg + majority + black1*class_agg + (1 | classid) + (1 |id) ,data=cool, weights=n, family=binomial))