# ------------------------------------------------------------------------------------------------ #
# Multilevel workshop: Stata examples
# Benjamin Rosche - benrosche.com
# Spring 2022
# 
# Examples adapted from Germán Rodríguez (Princeton)
# https://data.princeton.edu/pop510/
# 
# The dataset is freely available from:
# https://www.stats.ox.ac.uk/~snijders/mlbook.htm
# ------------------------------------------------------------------------------------------------ #

library(foreign)
library(lme4)
library(dplyr)
library(ggplot2)
library(ggthemes)

rm(list = ls())

# ------------------------------------------------------------------------------------------------ #
# Load data ----
# ------------------------------------------------------------------------------------------------ #

dat <-
  read.dta("snijders.dta") %>% 
  mutate(iqvc = iq_verb - mean(iq_verb))

# ------------------------------------------------------------------------------------------------ #
# Varying intercept model
# ------------------------------------------------------------------------------------------------ #

m1 <- lmer(gpa ~ 1 + ses + schoolses + (1 | schoolnr), REML=F, dat=dat)
summary(m1)

betas.m1 <- fixef(m1)

# ------------------------------------------------------------------------------------------------ #
# Plotting the varying intercept model
# ------------------------------------------------------------------------------------------------ #

dat.pred <- mutate(dat, yhat = predict(m1))

ggplot(dat.pred) + 
  geom_line(aes(x=ses, y=yhat, group = factor(schoolnr))) + 
  geom_abline(intercept=betas.m1["(Intercept)"],slope=betas.m1["ses"], color="red", size=1.3) +
  theme_minimal()

ggsave(filename = "R-varying-intercept.png", width = 7, height = 5, device='png', dpi=300)

# ------------------------------------------------------------------------------------------------ #
# Varying slope model
# ------------------------------------------------------------------------------------------------ #

m2 <- lmer(gpa ~ 1 + ses + schoolses + ses*schoolses + (1 + ses | schoolnr), REML=F, dat=dat)
summary(m2)

# ------------------------------------------------------------------------------------------------ #
# Plotting the varying slope model
# ------------------------------------------------------------------------------------------------ #

betas.m2 <- fixef(m2)

dat.pred <- mutate(dat, yhat = predict(m2))

ggplot(dat.pred) + 
  geom_line(aes(x=ses, y=yhat, group = factor(schoolnr))) + 
  geom_abline(intercept=betas.m2["(Intercept)"],slope=betas.m2["ses"], color="red", size=1.3) +
  theme_minimal()

ggsave(filename = "R-varying-slopes.png", width = 8, height = 4, device='png', dpi=300)


# ------------------------------------------------------------------------------------------------ #
# Varying intercept + slope model
# ------------------------------------------------------------------------------------------------ #

m2 <- lmer(gpa ~ 1 + ses + schoolses + (1 + ses | schoolnr), REML=F, dat=dat)
summary(m2)

# ------------------------------------------------------------------------------------------------ #
# Explaining varying intercept + slope model
# ------------------------------------------------------------------------------------------------ #

m3 <- lmer(gpa ~ 1 + ses + schoolses + ses*schoolses + (1 + ses | schoolnr), REML=F, dat=dat)
summary(m3)

# eof