# Edpsy/Psych/Stat 587
# Spring 2021
# c.j.anderson
#
# For results in lecture notes on infernce for marginal model using the hsb data
#
##########################################################
# Packages that we (may) need                            #
##########################################################
library(lme4)
library(lmerTest)
library(texreg)
library(optimx)
library(xtable)     #for Latex output of data frame (for modifying lectures notes using R)

##########################################################
# Functions that I wrote that compute useful statistics  # 
#  (easier than mannually running through R)             #
##########################################################

#source("C:\\Users\\cja\\Dropbox\\edps587\\lectures\\6 inference1\\R_robust_standard_errors_v7_use_this_one.txt")
#source("C:\\Users\\cja\\Dropbox\\edps587\\lectures\\6 inference1\\hlmRsq_v6.txt")
#source("C:\\Users\\cja\\Dropbox\\edps587\\lectures\\6 inference1\\contrast_v2.txt")
#source("C:\\Users\\cja\\Dropbox\\edps587\\lectures\\6 inference1\\bic_hlm.txt")

source("D:\\Dropbox\\edps587\\lectures\\6 inference1\\R_robust_standard_errors_v7_use_this_one.txt")
source("D:\\Dropbox\\edps587\\lectures\\6 inference1\\hlmRsq_v6.txt")
source("D:\\Dropbox\\edps587\\lectures\\6 inference1\\contrast_v2.txt")
source("D:\\Dropbox\\edps587\\lectures\\6 inference1\\bic_hlm.txt")

##########################################################
# Data steps                                             #
##########################################################
hsb1 <- read.table("D:\\Dropbox\\edps587\\lectures\\6 inference1\\HSB1.txt", header=TRUE)
hsb2 <- read.table("D:\\Dropbox\\edps587\\lectures\\6 inference1\\HSB2.txt", header=TRUE)

#hsb1 <- read.table("D:\\Dropbox\\edps587\\lectures\\6 inference1\\HSB1.txt", header=TRUE)
#hsb2 <- read.table("D:\\Dropbox\\edps587\\lectures\\6 inference1\\HSB2.txt", header=TRUE)
hsb <- merge(hsb1,hsb2, by=c('id'))

# Center ses
hsb$cses <- hsb$ses - hsb$meanses

# This is probably better to use to avoid "warning different scale"
hsb$zsize <- (hsb$size)/sd(hsb$size)

###########################################################
# Fitting 5 models                                        #
###########################################################
# These all use the Sattherwaite df, that's what lmerTest does:

summary(model1 <- lmer(mathach ~ 1 + cses + female + minority + meanses + zsize + sector + (1 |id), data=hsb, REML=FALSE))

summary(model2 <- lmer(mathach ~ 1 + cses + female + minority + meanses + zsize + sector + cses*zsize + (1 |id), data=hsb, REML=FALSE))

summary(model3 <- lmer(mathach ~ 1 + cses + female + minority + meanses + zsize + sector 
                           + cses*zsize + cses*sector + female*meanses + female*zsize + female*sector 
                           + minority*meanses + minority*zsize + minority*sector 
                           + (1 + cses|id), data=hsb, REML=FALSE))

# This work fine; however, for illustration of what you could do or try with a 
#   model that is close to convergence but not quite.
summary(model4 <- lmer(mathach ~ 1  + female + minority + cses + meanses + zsize + sector 
                           + cses*zsize + cses*sector + female*meanses + female*zsize + female*sector 
                           + minority*meanses + minority*zsize + minority*sector 
                           + (1 + female + minority + cses |id), data=hsb, REML=FALSE,
			control= lmerControl( 
			check.conv.singular= .makeCC(action = "message", tol = 1e-4),
			check.conv.hess    = .makeCC(action = "warning", tol = 1e-6),					
			check.conv.grad    = .makeCC("warning", tol = .002, relTol = NULL))))
				
# We could relax convergenece criteria on grad to make it look converged; for 
#   example set change tol = .002 to tol=.007	
	
# check.conv.singular will check whether hat(T) is a good covariance matrix, i.e.,
#   "rules for checking for a singular fit, i.e.  one where some parameters are on the
#   boundary of the feasible space (for example, random effects variances equal to 0 or 
#   correlations between random effects equal to +/- 1.0)"
#
#  Other things that you can try:			
# o control= lmerControl(optimizer = "nloptwrap",restart_edge = FALSE,boundary.tol = 0)
# o control= lmerControl(optimizer = "Nelder-Mead")
# o check to make sure model is not misspecified
# o Try package brms and use Bayesian estimation

# To me this model looks fine:
relgrad <- with(model4@optinfo$derivs,solve(Hessian,gradient))
max(abs(relgrad)) 

isSingular(model4)

############################################
#  --get table and                         #
#    add in Wald statistics                #
############################################
s4 <- summary(model4)
s4 <- as.data.frame(s4[10])
names(s4)<- c("Estimate","StdError","df t","t","Pr(>|t|)")
s4$df.Wald <- rep(1,nrow(s4))
s4$Wald4 <- s4$t**2

# Nice output
options(scipen = 999)
print(s4, type = "html", digits=2)
options(scipen = 0)


# For Carolyn's lecture notes #
xtable(s4,digits=3)

####################################
# Model based confidence intervals #
#   (alternatives are below)       #
####################################

s4$upper <- s4$Estimate - qnorm(.025)*s4$StdError
s4$lower <- s4$Estimate - qnorm(.975)*s4$StdError
s4
round(s4[,8:9],digits=2)

# For LaTex
xtable(s4[,8:9],digits=2)

##############################################
# Big model
##############################################
summary(model5 <- lmer(mathach ~ 1  + female + minority + cses + meanses + zsize + sector 
                           + cses*zsize + cses*sector + female*meanses + female*zsize + female*sector 
                           + minority*meanses + minority*zsize + minority*sector 
                           + (1 + female + cses |id), data=hsb, REML=FALSE) )

# Summarize what we've fit so far
screenreg(list(model1,model2,model3,model4,model5),
          custom.model.names=c("Model 1", "Model 2","Model 3","Model 4", "Model 5")	 )

### Easy way to get likelihood ratio test for nested models
anova(model1,model2)

	

######################################################################
# Constrasts using function Carolyn wrote                            #
#                                                                    #
#   order of elements in contrast matrix (columns) is the same as    #
#   the order in the model that you input                            #
######################################################################

cmodel <- lmer(mathach ~  1 + cses + female + minority + meanses + zsize + sector 
                            + cses*zsize + cses*sector + female*meanses + female*zsize
                            + female*sector + minority*meanses + minority*zsize 
							+  minority*sector   
                          + (1 + cses + female | id), data=hsb, REML=FALSE)
summary(cmodel)

## Wald test for each parameter:
# -- Create a (15 x 15) identity matrix
L1 <- matrix(0,nrow=1,ncol=15)
L1[1,15] <- 1
contrast(cmodel,L1)

# -- for better looking output
round(contrast(cmodel,L1), digits=2)


# Test all cross-level interactions:
L <- matrix(0,nrow=5,ncol=15)
L[5,14] <- 1
L[4,13] <- 1
L[3,11] <- 1
L[2,10] <- 1
L[1, 8] <- 1
L

round(contrast(cmodel,L),digits=2)

			   
################################################################
# Bootstrap confidence intervals-- this may take a long time,  #
#   so I did a simpler model                                   #
################################################################

#
boot.ci <- confint(model1, method="boot", seed=1125, nsim=1000, level=0.99)
round(boot.ci,digits=4)

# for LaTex code
xtable(boot.ci,digits=4)

################################################################
# Profile likelihood confidence intervals --- these are        #
#  faster than bootstrap                                       #     
################################################################

profile.ci<- confint(model1,  nsim=1000, level=0.99)
round(profile.ci,digits=4)

###############################################################
# Easy way to get likelihood ratio test for nested models     #
###############################################################

anova(model1,model2)


##############################################################
# Emprical/Robust/Sandwich standard errors using function    #
# written by carolyn                                         #
##############################################################

r2 <- robust(model2, hsb$mathach,hsb$id,"between/within")
round(r2$table, digits=4)

r3 <- robust(model3, hsb$mathach, hsb$id, "between/within")
round(r3$table, digits=4)

# LaTex
r3.data <- r3[1]$table            # turned list into matrix
xtable(r3.data,digits=2)

# Compare model and robust using model 3

r4 <- robust(model4,  hsb$mathach,hsb$id,"between/within")
round(r4$table, digits=4)

#
# Alternative method for computing sandwich...output isn't as
#   nice as mine. 
# ....... it is not working anymore ......
#
#library(sandwich)
#library(merDeriv)
#sand <- sandwich(model3, bread. = bread(model3, full = TRUE), meat. = meat(model3, level = 2))
# This yields the full covariance matrix for fixed and random effects)
#var <- diag(sand)
# Need to do a bit of work to get se of fixed effects and even more it want t-test, 
#  p-values, compare with model based, etc.
# Note that the last 4 (in this case) are for variance components
#se <- sqrt(var)
#se <- se[1:(length(se)-4)]


############################################################
# Compute BIC appropriate for HLMs ...and AIC              #
############################################################

bic.hlm(model1,hsb$id)
bic.hlm(model2,hsb$id)
bic.hlm(model3,hsb$id)
bic.hlm(model4,hsb$id)
bic.hlm(model5,hsb$id)


##############################################################
#  R1sq and R2sq using function written by Carolyn           # 
##############################################################
#
# REQUIRES that group/cluster/etc id is actually called "id"
#
id <- hsb$id

summary(model.a <- lmer(mathach ~ 1 + (1 | id),data=hsb,REML=FALSE))

summary(model.b <- lmer(mathach ~ 1 + cses + (1 | id),data=hsb,REML=FALSE))
hlmRsq(hsb,model.b)


summary(model.c <- lmer(mathach ~ 1 + cses + minority + (1 | id),data=hsb,REML=FALSE))
hlmRsq(hsb,model.c)


summary(model.d <- lmer(mathach ~ 1 + cses + minority + female + (1 | id),data=hsb,REML=FALSE))
hlmRsq(hsb,model.d)

summary(model.e <- model.e <- lmer(mathach ~ 1 + cses + minority +  female + sector + meanses + zsize + (1 | id),data=hsb,REML=FALSE))
hlmRsq(hsb,model.e)


summary(model.h <- lmer(mathach ~ 1 + cses + minority +  female + sector + meanses + zsize 
          + cses*sector + female*zsize + minority*sector +(1 | id),data=hsb,REML=FALSE))
hlmRsq(hsb,model.h)


summary(model.i <- lmer(mathach ~ 1 + cses + minority +  female + sector + meanses + zsize 
          + cses*sector + female*zsize + minority*sector 
          + cses*zsize + sector*female + meanses*female + zsize*minority + meanses*minority + (1 | id),data=hsb,REML=FALSE))
hlmRsq(hsb,model.i)


summary(model.bplus <- lmer(mathach ~ 1 + cses + meanses + (1 | id),data=hsb,REML=FALSE))
hlmRsq(hsb,model.bplus)

summary(model.extra <- lmer(mathach ~ 1 + cses + meanses + female + minority + pracad + sector +(1 | id),data=hsb,REML=FALSE))
hlmRsq(hsb,model.extra)


summary(model.bis<- lmer(mathach ~ 1 + cses + meanses + (1 + cses| id),data=hsb,REML=FALSE))
hlmRsq(hsb,model.bis)


summary(model.ris <- lmer(mathach ~ 1 + cses + meanses + female + minority + pracad + sector +(1 +cses| id),data=hsb,REML=FALSE))
hlmRsq(hsb,model.ris)


summary(model.ris2 <- lmer(mathach ~ 1 + cses + meanses + minority + female + minority + pracad + sector +(1+cses + minority | id),data=hsb,REML=FALSE))
hlmRsq(hsb,model.ris2)


######################################################################
#     Computing p-value for Mixture of chi-square for testing random #
######################################################################


(p <- .5*(pchisq(5.40,1, lower.tail=FALSE) + pchisq(5.40,2, lower.tail=FALSE)))   

(p <- .5*(pchisq(1.55, 3, lower.tail=FALSE) + pchisq(1.55, 4, lower.tail=FALSE)) )