---
title: "Edps/Psych/Stat 587: Fitting Basic HLM using R"
author: "Carolyn J. Anderson"
date: "1/6/2023"
output: 
    html_document:
      toc: yes
---
R code needed to fit HLM using R is discussed and illustrated here.  The code is essentially the same as in the txt file "R_hsb_code_in_notes.txt". 

The major packages that we use are lme4 and lmerTest.  We will use other later in the semester, but these two are used for almost all models fit this semester.  

The data used here is the High School and beyond data.  

# Setup
First we'll go ahead and set our working directory (i.e., this is where the data live on your computer)

```{r}
setwd("D:/Dropbox/edps587/lectures/SAS_MIXED_R")
```

We can also specify the packages that we'll use.  The lme4 package contains command lmer that we'll use to fit the models to data.  The output does not include tests of the fixed effects, but usng lmerTest will provide these.   

```{r setup, message=FALSE, warning=FALSE, include=TRUE}
library(lme4)
library(lmerTest)
library(texreg)
```
## Read in and set up data
 
The data are in two files, a level 1 and a level 2.  The following command reads in the level one data file.

### Read in student level data
```{r}
hsb1 <- read.table(file="HSB1data.txt", header=TRUE)

# Check file
head(hsb1)
```
### Read in the level 2 data.
```{r}
#
# Read in school level data
#
hsb2 <- read.table(file="HSB2data.txt", header=TRUE)

# Check file
head(hsb2)
```

### Merge Level 1 and Level 2
We need to merge the level 1 (student) and level 2 (school) data frames.  We need to match schools and students using the school id.  The following does this.
```{r}
hsb <- merge(hsb1,hsb2, by=c('id'))
names(hsb)
head(hsb)
```

## Creating some new variables

There are some variables that we would like to use, school means for ses and female (mean=proportion), as well as use these to create centered variables
```{r}
###########################################################   
# Create and add mean meanminority & meanfemale to data   #
##########################################################
meanfemale  <- aggregate(female ~ id , data=hsb,  
                         FUN="mean")

names(meanfemale) <- c('id', 'meanfemale')

hsb <- merge(hsb,meanfemale, by =c('id'))

# Since meanses is already in data file
hsb$ses.centered <- hsb$ses - hsb$meanses

head(hsb)
```

If you will reuse these data and don't want to go through all the data setup, you can save the data using 

```{r eval=FALSE}
write.table(hsb, 'hsb.txt', row.names=F, na='.')
```
This should show up in the working directory that you set above.  There are other formats that you can save data, but I like plain text.

# Fitting some models

## Empty/Null/random effects ANOVA
We can start with a simple random intercept model; that is our empty HLM, which is the null model which is equivalent to a random effects ANOVA.

```{r}
model.null <- lmer(mathach ~ 1 + (1 | id), data=hsb, REML=FALSE)

summary(model.null)
```
We will talk though what's in the out put in lecture.

## Computing ICC (and residual ICC)
For models with random intercepts we can compute an interclass correlation (or residual ICC) which interpretable as both (i) the correlation between observations within a cluster or group and (ii) the proportion of variance due to differences between groups.

There are at least 3 ways to do this
1. Brute force
```{r}
# Cut and paste from model summary
(icc <- 8.553/(8.553+39.148))
```
2. Extract the information to compute ICC
```{r}
vars <- as.data.frame(VarCorr(model.null))[4]
total <- sum(vars)
tau00 <- vars[1,1]
(icc <- tau00/total)
```

3. Use a little function that I wrote
```{r}
icc.lmer <- function(modl){
  vars <- as.data.frame(VarCorr(modl))[4]
  total <- sum(vars)
  tau00 <- vars[1,1]
  icc <- tau00/total
  return(icc)
  }

# example of how to use function
(icc.lmer(model.null))
```

Note that they all give the same result.   

You don't really need all those decimal places, so let's go with 
```{r}
round(icc,digits=3)
```
## Simple model with one fixed effect

```{r}
require(lmerTest)

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

round(icc.lmer(model.simple),digits=3)
```
## More complex model
```{r}
summary(model.complex <- lmer(mathach ~ 1 + ses.centered + minority + female +
    meanses + himinty + pracad + disclim +  sector + size + (1 | id), data=hsb, REML=FALSE))

round(icc.lmer(model.complex),digits=3)
```

# Compare Models Fit To Data

## Table
```{r}
screenreg(list(model.null,model.simple, model.complex),
          custom.model.names=c("Null", "One Fixed", "More Complex"),
          single.row=TRUE,
          custom.header=list("Our First HLM models in R"=1:3)
          )
```
## Likelihood Ratio Tests Between Models

The models are nested; that is, the null is a special case of the simple model, and the simple model is a special case of the more complex
```{r}
anova(model.null,model.simple, model.complex)
```

If you only include one model, you will get F tests for each coefficient in the model.

# Final Remarks

Warnings from lmer can sometimes be ignored such as  "Warning: Some predictor variables are on very different scales: rescaling".  It is generally good to follow this advise.

Things you should not ignore are errors message such as 

*     boundary (singular) fit: see help('isSingular')
*     Warning: Model failed to converge with 1 negative eigenvalue: -1.1e+01



