##############################################################
#### 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(nlme)
lme.null<-lme(data=HSB, mathach~1, random = ~ 1 | school)
summary(lme.null)
coef(lme.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(lme.null)~school, HSB, horizontal=T, 
    main="Homogeneity of Variance")
  qqnorm(resid(lme.null), main="QQplot for Null Model")

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

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

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

#### Run the Model adding a Random Component to the Slopes
lme.two<-lme(data=HSB, mathach~cses, random = ~ cses | school)
summary(lme.two)
coef(lme.two)
cbind(coef(lme.one), coef(lme.two))

anova(lme.one, lme.two)

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

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

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

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

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

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

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