---
title: 'Edpsy 587: Random Intercept and Slope Models using R'
author: "Carolyn J. Anderson"
date: "1/19/2023"
output: 
   html_document:
     toc: yes
---

The document contains R script used for the high school and beyond analyses given in the lecture notes.  

# Set up

## Libraries

I am adding one package that we haven't used before, "optimx".  This allows us to change the algorithm and solves as least one estimation problem that we run into with these models and data.

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

## Read in Data

This is where the data live on my computer:
```{r}
setwd('D:/Dropbox/edps587/lectures/4 randomslope')
```

The reading in of data and merging files is the same as in the "Random Intercepts Models in R" Rmarkdown file.


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)
```



We will want to center SES (eventually) but let's do it now
```{r, centerses}
hsb$cses <- hsb$ses - hsb$meanses
```


## Some exploratory looks at the data

```{r, nschools}
# Number of schools
(nschools <- nrow(hsb2))
```

Get a list of all the variables,
```{r, variables}
names(hsb)
```

We would like some basic descriptive statistics for numeric variables
```{r, descstats}
x <- hsb[, c(4,5,6,8,9,10,11)]
stargazer(x,type="text")
```

For my discrete variables
```{r, descdiscrete}
table(hsb$minority)
table(hsb$female)
table(hsb$sector)
```

# Graphs of the data

Is there variability in group (school) means for SES?
```{r, graph1}
plot(hsb$id, hsb$meanses, type='p', pch=21,
     main='HSB: Group Mean SES by School ID',
     xlab='School ID',
     ylab='Group Mean SES')
```



How much variability within school in terms of raw ses?
```{r, graph2}
plot(hsb$id, hsb$ses, type='p', pch=21,
     main='HSB: Student SES by School ID',
     xlab='School ID',
     ylab='Group Mean SES')
```



What does data look like for just one school and include linear regression line?
```{r, graph3}
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)
```



We might want to look at a bunch of schools but 160 graphs is a lot, so we'll take a random sample of 20 and just look at these. 

and to get the random sample
```{r, randsample}
groups <- unique(hsb$id)[sample(1:160,20)]
subset <- hsb[hsb$id%in%groups,]
```

We can now plot the data for these 20 schools.
```{r, graph4}
xyplot(mathach ~ ses | id, data=subset, col.line='black', 
       type=c('p'),
       main='Varability in Math ~ SES relationship')
```


However, we would like to also put in regression lines for each of these schools, which mean that we have to change what is in type,
```{r, graph5}
xyplot(mathach ~ ses | id, data=subset, col.line='black', 
       type=c('p','r'),
       main='Varability in Math ~ SES relationship')
``` 

What to graphs look like when we used school mean centered SES?
  ```{r, graph6}

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')
```  


It is also useful to look at all the regressions in a single plot; however, I need consequtive integers to index the school (i.e., index that goes from 1 - 160)
```{r, pregraph7}
sch_id <-  seq(1:nschools)
id <- unique(hsb$id)
sch_info <- as.data.frame(cbind(id,sch_id))
hsb <- merge(hsb,sch_info,by=c("id"))
```


To make the graph, we first set up a data frame that has the specifies the title, axis labels, etc.   Note that  type='n' doesn't plot any data

To make the graph, we first set up a data frame that has the specifies the title, axis labels, etc. ane then add in the regression lines  Note that  type='n' doesn't plot any data
```{r, graph7}
plot(hsb$ses,hsb$mathach,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:nschools){
  sub <- hsb[which(hsb$sch_id==i),]
  fitted <- fitted(lm(mathach~ses,data=sub))
  lines(sub$ses,fitted,col=i)
   }
```

The following code prodcues the graph on page 45 of the lecture note on random intercept and slope models. This is the variance function
```{r, graph8}
x <- seq(from=-2, to=2, by=.1)

total<- 8 + 2*(-0.5)*x + (1.5)*(x**2) + 30

plot(x,total,type = 'l', 
     cex.main = 1.5,
     col='blue',
     xlim=c(-2,2),
     xlab = 'Value of x',  
     ylab = "Estimated Variance of x",
     main = "Variance Function: Illustration of Heteroscedasticity"
)
```

# Create More Variables

The data file already has school centered SES (ie., cses).  But we also with create an aggregate variables for minority (which will be proportion of mintority students in a given school)
```{r, meanminority}
meanminority  <- aggregate(minority ~ id , data=hsb, mean)
names(meanminority) <- c('id', 'meanminority')

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

We do the same for females,
```{r, meanfemale}
meanfemale  <- aggregate(female ~ id , data=hsb, mean)
names(meanfemale) <- c('id', 'meanfemale')

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

At this point you might want to consider saving teh data to a file, especially if you plan to use it again,
```{r, savedata, eval=FALSE}
write.table(hsb, 'hsb.txt', row.names=F, na='.')
```

# Model Fitting

## Null HLM
```{r, nullhlm}
#
# Random inercept / empty HLM / Null model / RANDOM effects ANOVA
#

summary(model0 <- lmer(mathach ~ 1 + (1 | id), data=hsb, REML=FALSE))
```

## Computing ICC (and function)

Cut and paste
```{r, icc1}
icc <- 8.553/(8.553+39.148)
icc
```

Or extract the information to compute
```{r, icc2}
vars <- as.data.frame(VarCorr(model0))[4]
total <- sum(vars)
tau00 <- vars[1,1]
icc <- tau00/total
icc
```

Best yet, function since we do this over and over, use the following function that I wrote:
```{r, icc3}
icc <- 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 this
icc(model0)
```


## Random intercept and slope for centered ses

By centering here, we will be looking at the impact of where a student is relavtiave to their peers in terms of SES.
```{r, mod1}
summary(model1 <- lmer(mathach ~ 1 + cses + (1 | id), data=hsb, REML=FALSE))

icc(model1)

# for easier comparison
screenreg(list(model0,model1),
          single.row=TRUE,
          custom.model.names=c("Null","RI ses"),
          custom.header=list("Random Intercept Models"=1:2))
```

## Random intercept and slope & meanses for ses

To add a random slope for a level one predictor varaible, say $x_{1ij}$, we add it to the random part of the model specification; that is, "(1 + $x_{1ij}$ | id)".  In this case, we still have the random intercept but have added a random slope for $x_{1ij}$.

Below we fit a model with both a random intercept and random slope for cses.
```{r, mod2}
summary(model2 <- lmer(mathach ~ 1 + cses  + (1 + cses| id), data=hsb, REML=FALSE))
as.data.frame(VarCorr(model2))[4]
deviance(model2)

# for easier comparison with what we have already fit
screenreg(list(model0,model1,model2),
          single.row=TRUE,
          custom.model.names=c("Null","cses", "ces"),
          custom.header=list("Random Intercept"=1:2, "Random Slope"=3))
```


## Graphs around page 20 of Lecture Notes

```{r}
par(mfrow=c(1,1)) 

fixed <- fixef(model2)
mu2 <- fixed[1] + fixed[2]*hsb$cses 

plot(hsb$cses, mu2, type='l',
     col='black',
     xlim=c(-2, 2),
     ylim=c(0, 25),
     xlab='School Mean Centered SES',
     ylab='Fitted Math Scores',
     main='HSB: Population Average Regressions (Model 2)'
)
```

For the graph on or about page 21 of lecture notes, sample of fitted school specific regression lines 
```{r, graph21}
random <- ranef(model2, level=0:2)
U0j <- random$id[,1]
U1j <- random$id[,2]

plot(hsb$mathach ~ hsb$cses, type = 'n', 
     cex.main = 1.5,
     xlab = 'School Mean Centered SES',  
     ylab = "Fitted Math Scores",
     main = "HSB:  Fitted School Regression Lines (private,white)"
     )
	 
# draw in scho0l specific lines:  Private, white
for (i in 1:40){
    abline(a=(fixed[1]+ U0j[i]), b=(fixed[2]+U1j[i]),  col="blue")
  }
```

The graph on or around page 22 of Side by Side Comparison 
```{r, graph22}
par(mfrow=c(1,2))
plot(hsb$cses, mu2, type='l',
     col='black',
     xlim=c(-2, 2),
     ylim=c(0, 25),
     xlab='School Mean Centered SES',
     ylab='Fitted Math Scores',
     main='HSB: Population Average Regressions (Model 2)'
)

plot(hsb$mathach ~ hsb$cses, type = 'n', 
     cex.main = 1.5,
     xlab = 'School Mean Centered SES',  
     ylab = "Fitted Math Scores",
     main = "HSB:  Fitted School Regression Lines"
     )
# draw in scho0l specific lines:  Private, white
for (i in 1:40){
    abline(a=(fixed[1]+ U0j[i]), b=(fixed[2]+U1j[i]),  col='blue')
  }
```

The graph on or about Page 23 of estimated variance functions for 2 models
```{r, graph23}
par(mfrow=c(1,1))
ses2 <- seq(from=-2, to= 3, by=.1)
vars.ri <- as.data.frame(VarCorr(model1))[4]
total.ri <- rep(sum(vars.ri), 61)

vars.ris <- as.data.frame(VarCorr(model2))[4]
total.ris<- vars.ris[1,] + 2*vars.ris[3,]*ses2 + vars.ris[2,]*(ses2**2) + vars.ris[4,]

plot(ses2,total.ris,type = 'l', 
     cex.main = 1.5,
     col='blue',
     xlim=c(-2,3),
     ylim=c(45,52),
     xlab = 'School Mean Centered SES',  
     ylab = "Estimated Variance of Math Scores",
     main = "HSB: Variance Function"
     )
# random inercept
     abline(a=sum(vars.ri), b=0, col='red')

legend(-2,52, 
       c('Random Intercept & Slope','Random Intercept'),
       lty=c(1,1),
       lwd=c(2.5,2.5),
       col=c('blue','red'))
```


#  Discrete Level 1 predictors (fixed effects)                           

Note female and minority are dummy (i.e., I have not declared them to be factors in this model).

```{r}
model.d1 <- lmer(mathach ~ 1 + cses  + female + minority+ (1 + cses| id), data=hsb, REML=FALSE)
summary(model.d1)
as.data.frame(VarCorr(model.d1))[4]
n2lnlike <- -2* logLik(model.d1)
n2lnlike  # this is deviance
```

Since we have fit a number of models, let's see where we are at.

```{r}
screenreg(list(model1,model2,model.d1),
          custom.model.names=c("RI","RI and RS","Discrete fixed"))
```


## Plot of population average model

```{r}
fixed <- fixef(model.d1)
fem_mi <- fixed[1] + fixed[2]*hsb$cses + fixed[3] + fixed[4]

plot(hsb$cses, fem_mi, type='l',
     col='black',
     ylim = c(0,25),
     xlim = c(-2,2),
     xlab='School Centered SES',
     ylab='Fitted Math Scores',
     main='Population Average Regressions (discrete fixed)'
)
abline(a=(fixed[1] + fixed[3]), b=fixed[2], col='red')
abline(a=(fixed[1] + fixed[4]), b=fixed[2], col='blue')
abline(a=(fixed[1]), b=fixed[2], col='green')
legend(-2,25,
       c('Female, Minority', 'Female, non-Minority', 'Male, Minority', 'Male, non-Minority'),
       lty=c(1,1,1,1),
       col=c('black', 'red', 'blue', 'green')
    )
```


## Plot of variance function
```{r}
ses2 <- seq(from=-2, to= 3, by=.1)
vars.d1 <- as.data.frame(VarCorr(model.d1))[4]
total.d1<- vars.d1[1,] + 2*vars.d1[3,]*ses2 + vars.d1[2,]*(ses2**2) + vars.d1[4,]

plot(ses2,total.d1,type = 'l', 
     cex.main = 1.5,
     col='blue',
     xlim=c(-2,3),
     xlab = 'School Mean Centered SES',  
     ylab = "Estimated Variance of Math Scores",
     main = "HSB: Variance Function (discrete fixed)"
     )
```


## Alternative Treatment of categorical predictors 

Showing effect of random categorical predictors when they are declared to be factors (rather than putting in say dummy codes)
```{r, page66}  
hsb$gender <- as.factor(hsb$female)
hsb$race   <- as.factor(hsb$minority)
hsb$boy <- as.factor(-1*(hsb$female-1))
hsb$male <- -1*(hsb$female-1)

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

# To get variance 
as.data.frame(VarCorr(model.00))[4]
```
This should yield that same fit statistics, interpretation of effects, significance of effects, etc.  

# Random Descrete Predictors 
There are a number of ways to do this that lead to essentially the same result.  Below we look at some.

```{r}
summary(model.b1 <- lmer(mathach ~ 1 + cses  + female + race + (1 + cses + female| id), data=hsb, REML=FALSE))
as.data.frame(VarCorr(model.b1))[4]
```


## Alternative Model b1 
```{r}
summary(model.ab1 <- lmer(mathach ~ 1 + cses  + gender + race + (1 + cses + gender| id), data=hsb, REML=FALSE))
as.data.frame(VarCorr(model.ab1))[4]
```

## Another parameterization

```{r}
summary(model.b2 <- lmer(mathach ~ 1 + cses  + hsb$male + race + (1 + cses + hsb$male| id), data=hsb, REML=FALSE))
as.data.frame(VarCorr(model.b2))[4]
```

## Bad Parameteriation

Note that this one does not work --- does not converge.  Why?
```{r}
# Model c.
summary(model.c <- lmer(mathach ~ 1 + cses  + hsb$male + race + (0 + cses + male + female + minority| id), data=hsb, REML=FALSE))
as.data.frame(VarCorr(model.b2))[4]
```

## This Works.

Why?
```{r}
summary(model.all <- lmer(mathach ~ 1 + cses  + female + minority + (1 + cses + female + minority | id), data=hsb, REML=FALSE))
```

## And this works.


```{r}
summary(model.c <- lmer(mathach ~ 1 + cses  + female + minority + (0 + cses  + female + race| id), data=hsb, REML=FALSE))
```


You're on your own for the figures --- cut, paste & edit code for similar figures  


# Cross-level interactions      

When a level 2 variable, say $z_{j}$ is a predictor of a random effect (or fixed level 1 variable), we add the variable to the fixed/structural part of the formula. If $z_{j}$ is a predictor of the slope for $x_{ij}$, the formula would be 
$$y_{ij} \leftarrow 1 + x_{ij} + z_{j} + z_{j}*x_{ij} + (1 + x_{ij} | \mbox{id}) $$
In other words, we add an interaction between $z_{j}$ and $x_{ij}$.

In the example below, we add meanses as a predictor of the effect of cses:
```{r, cross-level}
summary(model.3 <- lmer(mathach ~ 1 + cses  + meanses + cses*meanses + (1 + cses| id), data=hsb, REML=FALSE))
as.data.frame(VarCorr(model.3))[4]
```
# Models Refinements

## Drop non-signicant Interaction

```{r, restOfThem}
summary(model.4 <- lmer(mathach ~ 1 + cses + meanses  + (1 + cses| id), data=hsb, REML=FALSE))
```

## p127 Big model  
R starts to get slow, can switch to microsoft's OpenR (free), which is about 4 times faster. The speed up is due to the use of a faster BLAS (basic linear algebra system), in particular an Intel MKL libraries.  Other alternative BLAS's are OpenBIAS and ATLAS.


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

We can do a bit better by rescaling school size which has a very different scale
````{r}
hsb$xsize <- scale(hsb$size,center=FALSE,scale=TRUE)
summary(model.5 <- lmer(mathach ~ 1 + cses  + female + minority + meanses  + xsize + sector 
               + pracad + disclim + himinty +(1 + cses + female + minority| id), data=hsb, REML=FALSE))
as.data.frame(VarCorr(model.5))[4]
```
# Other tricks to get models to converge 
ALthough this is a topic for another lecture, but I did need an additional one to get the complex model above to converge.

In this model (see aroun) p127 of lecture notes), we do make some changes to the model by adding sector as cross-level (i.e. predictor at level 2) and dropping random gender.

(This does lead to some differences in estimates of tau relative to SAS.)

The thing that really help is changing the optimizer.  This is done below using  "control" 
```{r}
hsb$xhiminty <- scale(hsb$himinty,center=FALSE,scale=TRUE)
summary(model.6 <- lmer(mathach ~ 1 + cses  + female + minority + meanses  + xsize + sector  + pracad + disclim + xhiminty + cses*sector + minority*sector + (1 + cses + minority| id), data=hsb, REML=FALSE,
                  control = lmerControl(optimizer ="Nelder_Mead") ))
as.data.frame(VarCorr(model.6))[4]
```

Round p127 A bit simpler model  --- only random intercept

```{r}
summary(model.7 <- lmer(mathach ~ 1 + cses  + female + minority + meanses  + xsize + sector 
               + pracad + disclim  + cses*sector + minority*sector 
               + (1| id), data=hsb, REML=FALSE))
as.data.frame(VarCorr(model.7))[4]
```



