##############################################################
#### Introductory Example of Multilevel Analysis
##############################################################

#### Read in the dataset
example.data<-read.table("http://www.hlm-online.com/training/example1.txt", header=T)
#### Make "GROUP" a factor
example.data$GROUP<-factor(example.data$GROUP)

#### Run the linear model (OLS; Ordinary Least Squares)
lm.model<-lm(SCIENCE~URBAN, example.data)
summary(lm.model)
plot(SCIENCE~URBAN, example.data)
abline(lm.model)

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

#### Run the Null Multilevel Model
  ## This is also called the multilevel ANOVA
library(nlme)
lme.null<-lme(data=example.data, SCIENCE~1, random = ~ 1 | GROUP)
summary(lme.null)
coef(lme.null)

  ## Compare Estimates
  cbind(means=tapply(example.data$SCIENCE, example.data$GROUP, mean), coef(lme.null))

  ## Check assumptions
  par(mfrow=c(1,2))
  boxplot(resid(lme.null)~GROUP, example.data, 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=example.data, SCIENCE~URBAN, random = ~ 1 | GROUP)
summary(lme.one)
coef(lme.one)

  ## Plotting the Fitted Regression Lines for Each Person
  fit1 <- unlist(fitted.values(lme.one))
  interaction.plot(example.data$URBAN, example.data$GROUP, fit1, 
    xlab="Urbanicity", ylab="Science", col=1:16)

  ## 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=example.data, SCIENCE~URBAN, random = ~ URBAN | GROUP)
summary(lme.two)
coef(lme.two)
cbind(coef(lme.one), coef(lme.two))

fit2 <- unlist(fitted.values(lme.two))
interaction.plot(example.data$URBAN, example.data$GROUP, fit2, 
  xlab="Urbanicity", ylab="Science", col=1:16)

anova(lme.one, lme.two)

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