##############################################################
#### High Schools and Beyond Example using R
##############################################################

#### Read in the dataset
HSB<-read.table("http://www.hlm-online.com/software/hsbdataset.txt", header=T)
#### Make "GROUP" a factor
HSB$school<-factor(HSB$school)

#### Run the linear model (OLS; Ordinary Least Squares)
lm.model<-lm(mathach~ses, HSB)
summary(lm.model)
plot(mathach~ses, HSB)
abline(lm.model)

  ## Checks Assumptions
par(mfrow=c(2,2))
plot(lm.model)

#### Run the Null Multilevel Model
  ## This is also called the multilevel ANOVA
library(lme4)
lmer.null<-lmer(mathach~1 + (1 | school), HSB)
summary(lmer.null)
coef(lmer.null)

  ## Compare Estimates of the Empirical Bayes means to actual means
  cbind(means=tapply(HSB$mathach, HSB$school, mean), coef(lme.null))

  ## Check assumptions
  par(mfrow=c(1,2))
  boxplot(resid(lmer.null)~school, HSB, horizontal=T, 
    main="Homogeneity of Variance")
  qqnorm(resid(lmer.null), main="QQplot for Null Model")

#### Run the Model with Fixed Effect for Slopes
lmer.one<-lmer(mathach~cses + (1 | school), HSB)
summary(lmer.one)
coef(lmer.one)

  ## Comparing Models
  anova(lmer.null, lmer.one)

  ## Check assumptions
  library(lattice)
  qqmath(resid(lmer.one), main="QQplot for First Model")

#### Run the Model adding a Random Component to the Slopes
lmer.two<-lmer(mathach~cses + (cses | school), HSB)
summary(lmer.two)
coef(lmer.two)
cbind(coef(lmer.one)$school, coef(lmer.two)$school)

anova(lmer.one, lmer.two)

  ## Check assumptions
  library(lattice)
  qqmath(resid(lmer.two), main="QQplot for Second Model")
  xyplot(resid(lmer.two)~fitted(lmer.two))

#### Adding "meanses" as a predictor of both the intercept and the slope of ses
lmer.three<-lmer(mathach~cses*meanses + (cses | school), HSB)
summary(lmer.three)
coef(lmer.three)

  ## Check assumptions
  qqmath(resid(lmer.three), main="QQplot for Third Model")
  xyplot(resid(lmer.three)~fitted(lmer.three))

anova(lmer.one, lmer.two, lmer.three)

#### Adding "sector" as a predictor of both the intercept and the slope of ses
lmer.four<-lmer(mathach~cses*meanses + cses*sector + (cses | school), HSB)
summary(lmer.four)
coef(lmer.four)

  ## Check assumptions
  qqmath(resid(lmer.four), main="QQplot for Fourth Model")
  xyplot(resid(lmer.four)~fitted(lmer.four))

anova(lmer.one, lmer.two, lmer.three, lmer.four)