---
title: 'Epsy587:  Random Intercept Models in R'
author: "Carolyn J. Anderson"
date: "12/29/2022"
output: 
   html_document:
     toc: yes
---

In this document, analyses will be done using High School and Beyond data set.  These analyses are reported in the lecture notes on *Random Intercept Models*.


# Set Up

## Libraries

The major R package that we will use is *lme4*.  We will use additional packages later in the semester, but lme4 is sufficient for most of what we'll do in the class.

The output will not include degrees of freedom, test statistics, and p-values for tests of fixed effects, we will also use the package *lmerTest*.

The package *texreg* will be used here to obtain a way to compare models.  An alternative to texreg is *stargazer*. 

The package *lattice* is required for some of the graphics.

```{r libs, message=FALSE, warning=FALSE, include=TRUE}
library(lme4)
library(lmerTest)
library(texreg)
library(stargazer)
library(lattice)
```

I will also set my working directory here
```{r directory}
setwd('D:\\Dropbox\\edps587\\lectures\\3 randomintercept')
```

## Read in Data

In HSB data for this example are contained in two files

1.   HSB1data.txt that contains level 1 (tudent level) data.
2.   HSB2data.txt  contains level 2 (school level) data.

These two files will be read in and then merged.

### hsb1

```{r hsb1}
hsb1 <- read.table(file="HSB1data.txt", header=TRUE)

# variable names
names(hsb1)

# dimensions of data frame
dim(hsb1)

# Top & bottom of the file
head(hsb1)
tail(hsb1)
```

### hsb2
```{r hsb2}
hsb2 <- read.table(file="HSB2data.txt", header=TRUE)

names(hsb2)

dim(hsb2)

# Top & bottom of the file 
head(hsb2)
tail(hsb2)
```

### Merging hsb1 and hsb2

We will make a long file.

I know that the data sets are already sorted by the school id, which will be the matching variable.

```{r merge}
hsb <- merge(hsb1,hsb2, by=c('id'))
dim(hsb)
names(hsb)
```

# Descriptive Statistics

The number of schools can be found in a number of ways, for example

```{r nschools}
(nschools <- nrow(hsb2))
#or 
(nschool2 <- length(unique(hsb$id)))
```

Number of students per school,
```{r perschool}
table(hsb$id)
```

And for a summary table of descriptive statistics ingoring clustering,
```{r sum1}
stargazer(hsb, type="html")
```


Summary statistics for math by school,
```{r byschool}
mmean <- aggregate(mathach~id, data=hsb, FUN="mean")
msd <- aggregate(mathach~id, data=hsb, sd)
mmin <- aggregate(mathach~id, data=hsb, min)
mmax <- aggregate(mathach~id, data=hsb, max)
math.stat <- cbind(mmean,msd[,2], mmin[,2], mmax[,2])

names(math.stat) <- c('Mean Math', 'SD math', 'Min', 'Max')
head(math.stat)
```

# Graphing Data

Below are a number of graphs that we looked at in class.

## Plot of mean ses by id --> variability between schools

```{r}
plot(hsb$id, hsb$meanses, type='p', pch=21,
     main='HSB: Group Mean SES by School ID',
     xlab='School ID',
     ylab='Group Mean SES')
```

## Plot of student ses by id --> variability within schools

```{r}
plot(hsb$id, hsb$ses, type='p', pch=21,
     main='HSB: Student SES by School ID',
     xlab='School ID',
     ylab='Group Mean SES')
```

##  Just one school

```{r}
school.1224 <- hsb[ which(hsb$id==1224), ]

plot(school.1224$ses, school.1224$mathach, type='p', pch=21,
     main="HSB: Math by SES for School 1224",
     xlab='Student SES',
     ylab='Math Achievement')
reg.1224 <- lm(school.1224$mathach~school.1224$ses, data=school.1224)
abline(reg.1224)
```


## Panel plot of a random sample of 20


### Will need this for graph 

```{r}
hsb$cses <- hsb$ses - hsb$meanses

# Random sample of 20 schools
groups <- unique(hsb$id)[sample(1:160,20)]

subset <- hsb[hsb$id%in%groups,]

xyplot(mathach ~ ses | id, data=subset, col.line='black', 
                      type=c('p'),
                      main='Varability in Math ~ SES relationship')
```  


### Panel plot of a random sample of 20 with regression lines

```{r}
xyplot(mathach ~ ses | id, data=subset, col.line='black', 
                      type=c('p','r'),
                      main='Varability in Math ~ SES relationship')
```  

### Panel using school centered ses
 
```{r} 
xyplot(mathach ~ cses | id, data=subset, col.line='black', 
                      type=c('p','r'),
                      main='Math Achievement by School Centered SES',
                      xlab='School Centered SES',
                      ylab='Math Achievement Test Scores')
```  

## All regression lines in one plot

```{r all_inone}
hsb$schid <- as.factor(hsb$id)

mathbysch <- lm(formula = mathach ~ 1*schid + ses*schid, data = hsb)
summary(mathbysch)

fschmath <- fitted(mathbysch)
school.list <- unique(hsb$schid)
length(school.list)


coef <- coefficients(mathbysch)
intercepts <- c(0, coef[3:161])
slopes <- c(0, coef[162:320])


# Plot all regressions in one plot
plot(data = hsb, mathach~ses,type = 'n', 
     ylim = c(-10, 30),
     xlim = c(-3, 3),
     cex.main = 1.5,
     xlab = 'SES', 
     ylab = "Math Scores",
     main = "Separate Regression for Each School "
     )

for(i in 1:length(school.list)){
  abline(a = (coef[1] + intercepts[i]), b = (coef[2]+slopes[i]), col = 'blue')
  par<- par(new=F)
}
abline(v=c(0), lwd=2, color='black')
```

# Create some variables


The data file already has mean SES, but need school mean centered ses


We might want to have a variables that indicates the proportion (mean) of students that are minority within each school.  There will be 160 such proportions and we will need to merge them into our main data set.
```{r}
meanminority  <- aggregate(minority ~ id , data=hsb, FUN="mean")
names(meanminority) <- c('id', 'meanminority')

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

head(hsb)
```

The same goes for gender,
```{r}
meanfemale  <- aggregate(female ~ id , data=hsb, FUN="mean")
names(meanfemale) <- c('id', 'meanfemale')

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

head(hsb)
```

 If you want to save data so that you don't got through steps
  to create is again, you can write the results to a file.
  There are other file formats that you can use (e.g., csv)
  
```{r eval=FALSE}
write.table(hsb, 'hsb.txt', row.names=F, na='.')
```

# Model Fitting

We will start simple and work to more complex

## Empty/Null HLM

This model has no predictor variables.   The formula that is entered is very simple to what you would enter for *lm*; however, what is entered here is a linear mixed model.  The random part of the model is indicated by the "(1 | id)".  The value "1" indicates an interecpt so (1 | id) means random intercept; that is $U_{0j}$.

By default a fixed intercept, $\gamma_{00}$, is included; however, I made it explicit below.  

Note that I included "REML=FALSE".  I want to perform maximum likelihood estimation.  Late we will talk about different estimation methods.
```{r empty}
model0 <- lmer(mathach ~ 1 + (1 | id), data=hsb, REML=FALSE)

summary(model0)
```
The values in the table labeled "Random effects" are the variances and standard deviations of the random effects.  The row labeled  "id (intercept)" contains the $\hat{\tau}_{0}^2$ and the row labeled "residual" contains $\hat{\sigma}^2$.   The columns labeled "Standard deviation" contain $\sqrt{\hat{\tau}_{0}^2}$) and $\sqrt{\hat{\sigma}^2}$).
 

The values in the table labeled "Fixed effects" is the $\hat{\gamma}_{00}$ an, it's standard error, The columns, "df", "t value" and "Pr(>|t|)" were generated by *lmerTest*.  Note there is information about this at the top of the output.

## ICC

I know that I've seen a packge that computes this, but it's easy to do so.  There are (at least) three ways to compute this that range from brute force to writing your own function.

### Cut-n-paste

Cut the values of $\hat{\tau}_{0}^2$ and $\hat{\sigma}^2$ from the summary of the model and use these
```{r}
# Cut and paste
icc <- 8.553/(8.553+39.148)
icc
```
### Extract rather than cut and paste
The next method extracts information and puts it into formula for an ICC
```{r}
vars <- as.data.frame(VarCorr(model0))[4]
total <- sum(vars)
tau00 <- vars[1,1]
icc <- tau00/total
icc
```

### Write a function
The last method, estentially take the last method and turns it into a little function
```{r}
icc.lmer <- function(model){
  vars <- as.data.frame(VarCorr(model))[4]
  total <- sum(vars)
  tau00 <- vars[1,1]
  icc <- tau00/total
  return(icc)
  }

# example of how to use this
icc.lmer(model0)
```

Note that the first method will contain some round off error; however, it's good enough.

## Random intercept with Fixed Slope for SES

Just add "ses" to the fixed part of the model,
```{r} 
model1 <- lmer(mathach ~ 1 + ses + (1 | id), data=hsb, REML=FALSE)
summary(model1)
icc.lmer(model1)
```

Since we now have two models, it would be nice to have an easy way to compare them
```{r}
screenreg(list(model0,model1), single.row=TRUE,
          custom.model.names=c("Null","+SES"),
          custom.header=list("Random Intercept Models"=1:2)
)
```

## Random intercept with level 1 and level 2 predictors

This is where being able to write the linear mixed model is good thing.  The level 2 variable, meanses, is just added to the fixed part of the model.

```{r}
model2 <- lmer(mathach ~ 1 + ses + meanses + (1 | id), data=hsb, REML=FALSE)

summary(model2)
```

## Use School centered SES

as a level 1 predictor,
```{r}
hsb$cses <- hsb$ses - hsb$meanses

model3 <- lmer(mathach ~ 1 + cses  + (1 | id), data=hsb, REML=FALSE)

summary(model3)
icc.lmer(model3)
```

## Some Graphics:  

We will plot both the population average models as well as school specific ("conditional") regression line(s).

We need to get fixed coefficients from the output,
```{r}
intercept <- fixef(model3)[1]
slope <- fixef(model3)[2]
```

For the conditional/school specific we will need estimates of the random effects, $\hat{U}_{0j}$.

```{r}
intercept.j <- ranef(model3)
student.int <- intercept.j$id[1]
```

We will start by setting up a graphics window (i.e., type='n').  In the first graph, we only plot the population average regression line.  In the 2nd figure, we graph both the popuation average line and the school specific regressions.

### Popuation Average Regression Line
```{r, fitted, warnings=FALSE, message="FALSE"}
plot(mathach~cses,type = 'n', data = hsb, 
     ylim = c(-10, 30),
     xlim = c(-3, 3),
     cex.main = 1.5,
     xlab = 'School Centered SES',  
     ylab = "Fitted Math Achievement Scores",
     main = "Population Average Regression "
     )
     
# Population average line
for(i in min(hsb$cses):max(hsb$cses)){
   abline(a=intercept, b=slope, color='black',lwd=4)
}
```

### Population and School Specific

```{r, warning=FALSE, comment=FALSE, warning=FALSE}
plot(mathach~cses,type = 'n', data = hsb, 
     ylim = c(-10, 30),
     xlim = c(-3, 3),
     cex.main = 1.5,
     xlab = 'School Centered SES',  
     ylab = "Fitted Math Achievement Scores",
     main = "Population & School Specific Regressions "
     )
# draw in school specific line
for (i in 1:length(school.list)){
   abline(a=(intercept+student.int[i,]), b=slope, type='l', col=i)
}

# Population average line
for(i in min(hsb$cses):max(hsb$cses)){
   abline(a=intercept, b=slope, col='black',lwd=4)
}
```
The black line is the population average.
 
# Model 4: center ses & add back in mean ses

```{r}
model4 <- lmer(mathach ~ 1 + cses + meanses + (1 | id), data=hsb, REML=FALSE)

summary(model4)
icc.lmer(model4)
```

# Model 5:  Lots of level 1 variables

```{r}
model5 <- lmer(mathach ~ 1 + cses + female + minority + meanses + (1 | id), data=hsb, REML=FALSE)

summary(model5)
icc.lmer(model5)
```

# Model 6: Lots of level 1  & 2 variables

```{r}
model6 <- lmer(mathach ~ 1 + cses + female + minority + meanses + himinty  + pracad 
                         + disclim + sector + size +(1 | id), data=hsb, REML=FALSE)

summary(model6)
icc.lmer(model6)
```


# Model 7: Kitchen Sink

We add even more variables.  Note this is not how we should go about modeling the data; however, this is for showing you how to fit complex models.

We will ignore the warning messages for now...
```{r}
model7 <- lmer(mathach ~ 1 + cses + female + minority + meanses + meanminority + meanfemale   + himinty  + pracad + disclim + sector + size +(1 | id), data=hsb, REML=FALSE)

summary(model7)
icc.lmer(model7)
```
# Quick comparison

```{r}
screenreg(list(model0,model5,model6,model7),
     single.row=TRUE,
     custom.model.names=c("Empty/Null", "Add Level 1", "Level 1 & 2","Kitchen Sink"),
     custom.head=list("Simple"=1, "Complex"=2:4)
          )
```