#!/usr/local/bin/R

# Xavier Fernández i Marín
# December 2004
# http://xavier-fim.net/

# reply the analysis done with HLM

# load the example datasets from HLM/AppendixA/Examples
# Download them directly. Requires Internet connection
hsb1 <- read.table(file="http://xavier-fim.net/R/scripts/hsb1.dat")
hsb2 <- read.table(file="http://xavier-fim.net/R/scripts/hsb2.dat")

# get the variables and basic information of each dataset
names(hsb1)
length(hsb1$ID)
summary(as.data.frame(hsb1))

names(hsb2)
length(hsb2$ID)
summary(as.data.frame(hsb2))

# load the package to do Linear Mixed-Effects Models
library(nlme)

# create a new dataset with the merged data
ds <- merge(as.data.frame(hsb1), as.data.frame(hsb2))

### model specification

# 1
mod.1 <- lme(fixed = MATHACH ~ 1,
             random = ~ 1 | ID,
             data=ds, method="ML"
             )

# 2
mod.2 <- lme(fixed = MATHACH ~ SES,
             random = ~ 1 | ID,
             data=ds, method="ML"
             )

# 3
mod.3 <- lme(fixed = MATHACH ~ SES + SECTOR,
             random = ~ 1 | ID,
             data=ds, method="ML"
             )

# 4
mod.4 <- lme(fixed = MATHACH ~ SES + SECTOR,
             random = ~ 1 + SES | ID,
             data=ds, method="ML"
             )

# 5
mod.5 <- lme(fixed = MATHACH ~ SES * SECTOR,
             random = ~ 1 + SES | ID,
             data=ds, method="ML"
             )


# Extract information about the models
print(mod.1)      # basic information
summary(mod.1)    # complete information
summary(mod.2)
summary(mod.3)
summary(mod.4)
summary(mod.5)

# attributes of an 'lmeObject'
mod.1$dims     # dimensions 
mod.1$coef     # coefficients
mod.1$call     # the formula used to produce the object
mod.1$residuals     # residuals of the model
mod.1$varFix   # matrix with the fixed effects
mod.1$sigma    # estimated within-group error standard deviation

anova(mod.1)


# Plots
plot(mod.1)
     # by default, the plot of an lme object is 
     # form=resid(., type = "p") ~ fitted(.)
     # that corresponds to the standardized residuals vs predicted values

plot(mod.1, resid(.) ~ fitted(.) | SECTOR)
     # idem, but comparing the sectors 

qqnorm(mod.1)     # normality of the residuals




interaction.plot(ds$SES, ds$SECTOR, ds$MATHACH, col=c("blue", "red"))
     # show one of the plotting abilities of R

# standardized residuals vs predicted values, by sector
plot(mod.3, resid(., type = "p") ~ fitted(.) | SECTOR, abline = 0)

# boxplots with the residuals by school 
#  (it's just an example, there are too many objects in level 2)
plot(mod.4, ID ~ resid(.))