# Edps 590BAY
# Spring 2018
# C.J.Anderson
#
#   MODELS MAY NOT BE CONVERGING
#
# Look at NELS with 23 schools using brms package  which requires Stan.
#     I used "bare-bones" model with brms (i.e., there is a lot more that
#     can be done with the package.  
#
# For more information on the brms package see:
#   Durkner, P.C. (2017).  brms: An R package for Bayesian multilevel models Using Stan
#                      Journal of Statistical Software, 80, DOI:  10:18637/jss.v080.i01
#   https://cran.r-project.org/web/packages/brms/brms.pdf
#
# It install Stan see https://github.com/stan-dev/rstan/wiki/Installing-RStan-on-Windows
#   (there are links to installation for MACs here)
#

#  What the variables are  (data from Kreft & DeLeeuw's text):
#    SCHOOL    = SCHOOL ID
#    STUDENT   = STUDENT ID
#    SEX       = STUDENT SEX
#    RACE      = STUDENT RACE
#    HOMEW     = TIME ON MATH HOMEWORK
#    SCHTYPE   = SCHOOL TYPE
#    SES       = SOCIOECONOMIC STATUS
#    PARED     = PARENTAL EDUCATION
#    MATH      = MATH SCORE
#    CLASSSTR  = CLASS STRUCTURE
#    SCHSIZE   = SCHOOL SIZE
#    URBAN     = URBANICITY
#    GEO       = GEOGRAPHIC REGION
#    MINORITY  = PERCENT MINORITY
#    RATIO     = STUDENT-TEACHER RATIO
#    PUBLIC    = PUBLIC SCHOOL BINARY
#    WHITE     = WHITE RACE BINARY
#
#
library(brm)
library(rstan)
library(lme4)
library(lmerTest)
library(coda)
library(rjags)
library(runjags)
library(mvtnorm)
library(lattice)


#setwd("C:\\Users\\cja\\Dropbox\\edps 590BAY\\Lectures\\10 Multilevel models")
setwd("C:\\Users\\cja\\Documents\\Dropbox\\edps 590BAY\\Lectures\\10 Multilevel models")

# Read in data:
nels <- read.table("school23_data.txt",header=TRUE)
head(nels)
tail(nels)

# sort by schools:
nels <- nels[order(nels$school,nels$student) , ]

# Get information on number schools & create consecutive integer school.id
school <- unique(nels$school)                        ##  
N <- length(school)                                  ## need N and school.id
school.int <- as.data.frame(cbind(school,seq(1:N)))  ##  for dataList & model  
names(school.int) <- c("school", "school.id")        ##
nels <- merge(nels,school.int,by="school")           ##

# total number of students
n <- length(nels$math)                               ## need n for dataList & model

# Double-check set-up
head(nels,n=50)


# You don't need this but might want it
nj <- matrix(999,nrow=N,ncol=1)
for (j in 1:N) {
  nj[j] <- nrow(subset(nels,nels$school.id==j))
}  

#   Lattice plot with regressions
xyplot(math ~ homew | school.id, data=nels, 
       col.line='black', 
	   type=c('p','r'), 
	   main='Variability in math ~ homew| school',
	   xlab="Time Spent Doing Homework",
	   ylab="Math Scores")

# Graph to illustrate need for random intercept and slope
plot(nels$homew,nels$math,
	type="n", 
	col="blue",
	lwd=2,
	main="NELS: Linear Regression by School \n(for where there is data)",
	xlab="Time Spent Doing Homework",
	ylab="Math Scores"
	)
for (j in 1:N) {
  sub <- subset(nels,nels$school.id==j)
  fitted <- fitted(lm(math~homew,sub))
  lines(sub$homew,fitted,col=j)
}  


###################################################################
# Random intercept     math~ 1 + homew + (1 |school.id)           #
###################################################################
modelx.lmer <- lmer(math~ 1 + homew + (1 |school.id), data=nels, REML=TRUE)
summary(modelx.lmer)

modelx.brm <- brm(math ~ homew + (1 | school.id),data=nels,save_all_pars=TRUE,cores=4)
summary(modelx,waic=TRUE)
plot(modelx)



##########################################################################
# Random intercept & slope    math~ 1 + homew + (1 + homew|school.id)    #
##########################################################################
model1.lmer <- lmer(math~ 1 + homew + (1 + homew|school.id), data=nels, REML=TRUE)
summary(model1.lmer)


model1.brm <- brm(math ~ homew + (1 + homew | school.id),data=nels,
   save_all_pars=TRUE,
   iter=10000,
   warmup=2000,
   cores=4,
   thin=5
)

summary(model1.brm,waic=TRUE)           # summary of model information
plot(model1.brm)                        # trace and density plots




##########################################################################
# Random intercept & slope with level 1 and level 2 predictors           #
##########################################################################
nels$public <- ifelse(nels$schtype==1,1,0)
nels$public.fac <- as.factor(nels$public)

model2.lmer <- lmer(math~ 1 + homew + ses + public.fac  + homew*public.fac + (1 + homew|school.id), data=nels, REML=TRUE)
summary(model2.lmer)

model2.brm <- brm(math~ 1 + homew + ses + public.fac  + homew*public.fac + (1 + homew|school.id), 
                data=nels, 
		cores=4,                save_all_pars=TRUE)  # fit the model
summary(model2.brm)                  # summary of model information
plot(model2.brm)                     # trace and density plots
bayes_R2(model2.brm)                 # Bayes version of R2
bayes_factor(model2.brm,model1.brm)  # Computes bayes factor     



# Summary from last model:

#Compiling the C++ model
#Start sampling
#> summary(model2.brm)
# Family: gaussian 
#  Links: mu = identity; sigma = identity 
#Formula: math ~ 1 + homew + ses + public.fac + homew * public.fac + (1 + homew | school.id) 
#   Data: nels (Number of observations: 519) 
#Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; 
#         total post-warmup samples = 4000
#    ICs: LOO = NA; WAIC = NA; R2 = NA
# 
#Group-Level Effects: 
#~school.id (Number of levels: 23) 
#                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
#sd(Intercept)            7.46      1.34     5.33    10.50       1016 1.01
#sd(homew)                4.01      0.74     2.77     5.72        736 1.01
#cor(Intercept,homew)    -0.83      0.08    -0.94    -0.64       1126 1.01
#
#Population-Level Effects: 
#                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
#Intercept            47.57      2.91    41.97    53.64        664 1.01
#homew                 2.33      1.53    -0.78     5.27        665 1.01
#ses                   2.67      0.51     1.65     3.67       2944 1.00
#public.fac1          -0.89      3.65    -8.29     6.05        666 1.01
#homew:public.fac1    -0.72      1.90    -4.36     3.22        708 1.01
#
#Family Specific Parameters: 
#      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
#sigma     7.19      0.24     6.74     7.68       4000 1.00
#
#Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
#is a crude measure of effective sample size, and Rhat is the potential 
#scale reduction factor on split chains (at convergence, Rhat = 1).