---
title: 'Edps/Psych/Stat 587:  Statistical Inference for Marginal Model'
author: "Carolyn J. Anderson"
date: 'Feb 3, 2023'
output: 
    html_document:
      toc: yes
---
# Overview
This document reproduces examples in lecture and illustrates how to use various functions that I wrote.  These functions do the following:


*     icc:  Computes interclass correlation for 2-level models
*     bic.hlm(model,cluster.id): Computes BIC statistic that balences between and within
*     robust(model.fit, response, idvar, dftype): Computes Sandwich standard errors and reports model based and sandwich bases statistics.
*     contrast(model,L) where L is c x p matrix of contrast.  This gives Wald and F-tests results 
*     hlmRsq(dataset, model, id): Computes R_1^2 and R_2^2

All of these functions are in a single file called "All_functions.txt". This can be set as a source file.   For more information on these functions see the lectures notes and notes in "All_functions.txt".

For some of these functions, I implicitly assume that the data were such that rows from same group or cluster are grouped together in a block of the data frame. (For large data sets this might reduce computation time needed to fit models to data). 

# Set up

## Set the source file for useful functions 
```{r source_file}
setwd("D:\\Dropbox\\edps587\\lectures\\6 inference1")
source("All_functions.txt")
```
## Package we'll use
```{r packages}
library(lme4)
library(lmerTest)
library(texreg)
library(optimx)
library(xtable)    # for Latex output of data frame
```
## Data steps

```{r data_steps}
# Read in data files
hsb1 <- read.table("HSB1.txt", header=TRUE)
hsb2 <- read.table("HSB2.txt", header=TRUE)

# Merge into one file
hsb.tmp <- merge(hsb1,hsb2, by=c('id'))

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

# Make sure data are ordered/sorted by id.
hsb <- hsb.tmp[order(hsb.tmp$id),]
hsb$id[1:50]
```

One last step for this example: to avoid "warning different re-scale"
```{r scale_size}
hsb$zsize <- (hsb$size)/sd(hsb$size)
```
# Four Models 

We will work with the following set of models
```{r first_three_models}
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))

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))
```	
## Trouble with Convergence
We could relax convergence criteria on the gradient to make it look converged; for example, change tol = .002 to tol=.007.	
	
Other things we can do is to Check.conv.singular; that is, check whether hat(T) (estimated matrix of $\tau$s) is a "good" covariance matrix.  Bad covrariance matrices are non-singular and those 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


```{r one.more.4}
summary(model4.exta <- 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))))
```	

To me this model looks fine, but changing criteria to make things look good if not good practice.
```{r cheating}
relgrad <- with(model4@optinfo$derivs, solve(Hessian, gradient))
max(abs(relgrad)) 

isSingular(model4)
```

# Wald Statistics

## Hypothesis test
The wald statistic for a parameter, for example $H_{o}:  \gamma_{kl}=\gamma_{kl}^{*}$ where $\gamma_{kl}^{*}$ is a fix values (usually 0) equals

$$ z = \frac{\gamma_{ik} - \gamma_{kl}^{*}}{\hat{se}} \sim N(0,1)$$
But usually, 
$$ X^2 = z^2 \sim \chi^2_1 $$
which is what is computed below.
```{r wald.for.each}
# Put summary into object
s4 <- summary(model4)

# Get the cofficients from the summary
s4 <- as.data.frame(s4[10])

# Give them names
names(s4)<- c("Estimate","StdError","df t","t","Pr(>|t|)")

# Compute wald statistics for each
s4$df.Wald <- rep(1,nrow(s4))
s4$Wald4 <- s4$t**2

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


# For Carolyn's lecture notes #
# xtable(s4,digits=3)
```
## Model Based Confidence Limits

```{r ci.wald}
s4$upper <- s4$Estimate - qnorm(.025)*s4$StdError
s4$lower <- s4$Estimate - qnorm(.975)*s4$StdError

options(scipen = 999)  # turn of scientic notation
round(s4[,8:9],digits=2)
options(scipen = 0)    # turn it back one


# For LaTex
# xtable(s4[,8:9],digits=2)
```
# A Big Model (model 5)

```{r model5}
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
```{r all.five}
screenreg(list(model1,model2,model3,model4,model5),
          custom.model.names=c("Model 1", "Model 2","Model 3","Model 4", "Model 5")	 )
```

```{r lr.test}
# Easy way to get likelihood ratio test for nested models
anova(model1,model2)
```
This is a test whether the interaction is significant or not.
	
# Contrasts
This will use the function "contrast' from the  "All_functions.txt" source file.

The model we will work with here
```{r cmode.fit}
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)
```
## To Help with Defining L
 You need to create a matrix where the number of columns equals the number of fixed effect.  A fast way to get number of coefficents is 
```{r nfixed}
length(fixef(cmodel)) 
```
And if you want the names either check output summary table of fixed effects, or 

```{r names_fixed}
rownames(summary(cmodel)$coefficients)
```
The number of rows equals the number of tests you want to perform (simultaneously).

## Test for just one $\gamma_{kl}$
We will test the $\gamma_{kl}$ for the interaction between minority and sector.  The contrast matrix L for this
```{r L1}
L1 <- matrix(0,nrow=1,ncol=15)
L1[1,15] <- 1

# This give us
L1

# And the test
contrast(cmodel,L1)

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

## Test all cross-level interactions:
```{r wald.cross.level}
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)
```

# More Ways to do inference for Fixed Effects

## Bootstrap CI
This can take longer depending on method.  We did the wald based above and these were pretty fast. The bootstrap CIs an take a very long time

```{r boot.ci, eval=FALSE}
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 CIs
An alternative are Profile likelihood confidence intervals, which takes less time to compute than bootstrap.            
```{r profile.ci}
profile.ci<- confint(model1,  nsim=1000, level=0.99)
round(profile.ci,digits=4)
```

## Likelihood ratio test for nested models 
This is easy
```{r}
anova(model1,model2)
```

# Sandwich Standard Errors

For these we'll use the function "robust" from the source file. 
```{r Sandwich2}
r2 <- robust(model2, hsb$mathach, hsb$id,"between/within")
round(r2, digits=4)
```
```{r sandwich3}
r3 <- robust(model3, hsb$mathach, hsb$id, "between/within")
round(r3, digits=4)

# LaTex
# r3.data <- r3[1]$table      # turn list into matrix
# xtable(r3.data,digits=2)
```
Compare model and robust using model 3
```{r model.vs.robust.model4}
r4 <- robust(model4,  hsb$mathach,hsb$id,"between/within")
round(r4, digits=4)
```

# BIC For Multilevel

For these we use the bic.html function from the source file
```{r bics}
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)
```

# R squares for multilevel
```{r}
# Null or empty HLM
summary(model.a <- lmer(mathach ~ 1 + (1 | id),data=hsb,REML=FALSE))

# Random intercept with cses are fixed effect
summary(model.b <- lmer(mathach ~ 1 + cses + (1 | id),data=hsb,REML=FALSE))
hlmRsq(hsb, model.b, hsb$id)

# Random intercept with more micro level variables
summary(model.c <- lmer(mathach ~ 1 + cses + minority + (1 | id),data=hsb,REML=FALSE))
hlmRsq(hsb,model.c, hsb$id)

# Random intercept with micro and marco varialbles
summary(model.d <- model.d <- lmer(mathach ~ 1 + cses + minority +  female + sector + meanses + zsize + (1 | id),data=hsb,REML=FALSE))
hlmRsq(hsb,model.d, hsb$id)

# Random intercept with some cross-level (fixed) interactions
summary(model.e <- lmer(mathach ~ 1 + cses + minority +  female + sector + meanses + zsize 
          + cses*sector + female*zsize + minority*sector +(1 | id),data=hsb,REML=FALSE))
hlmRsq(hsb,model.e, hsb$id)

# Random intercept with lots of variables and interactions
summary(model.f <- 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.f, hsb$id)


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

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

# Random intercept and slope for cses
summary(model.bis<- lmer(mathach ~ 1 + cses + meanses + (1 + cses| id),data=hsb,REML=FALSE))
hlmRsq(hsb, model.bis, hsb$id)


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


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, hsb$id)
```

# Tests on Variance Components

Computing p-value for mixture of Chi-squares.  

This one is for $H_o: \tau^2_1 = \tau_{01} = 0$ from page 97 of lecture notes (i.e., testing null model (random intercept) vs random intercept & slope for female):

```{r mixture.chi.squares}
(p <- .5*(pchisq(5.40,1, lower.tail=FALSE) + pchisq(5.40,2, lower.tail=FALSE)))   
```

This is the test reported on page 101 of lecture notes. 
```{r another.variance.test}
(p <- .5*(pchisq(1.55, 3, lower.tail=FALSE) + pchisq(1.55, 4, lower.tail=FALSE)) )  
```