# Read in Classroom data.
class <- read.csv("E:\\LMM Data\\classroom.csv", h = T)

# Load nlme package.
library(nlme)

# Fit Model 4.1.
model4.1.fit <- lme(mathgain ~ 1, random = ~ 1 | schoolid/classid,
                        class, method = "REML")

summary(model4.1.fit)

random.effects(model4.1.fit)

# Fit Model 4.1A.
model4.1A.fit <- lme(mathgain ~ 1, random = ~1 | schoolid,
                         data = class, method = "REML")

# Test Hypothesis 4.1.
anova(model4.1.fit, model4.1A.fit)

# Model 4.2.
model4.2.fit <- lme(mathgain ~ mathkind + sex + minority + ses,
                        random = ~1 | schoolid/classid, class,
                        na.action = "na.omit", method = "REML")

summary(model4.2.fit)

# Model 4.1: ML estimation with lme().
model4.1.ml.fit <- lme(mathgain ~ 1,
                           random = ~1 | schoolid/classid, class, method = "ML")

# Model 4.2: ML estimation with lme().
model4.2.ml.fit <- lme(mathgain ~ mathkind + sex + minority + ses,
                           random = ~1 | schoolid/classid, class,
                           na.action = "na.omit", method = "ML")

anova(model4.1.ml.fit, model4.2.ml.fit)

# Model 4.3.
model4.3.fit <- update(model4.2.fit,
                           fixed = ~ mathkind + sex + minority + ses + yearstea + mathprep + mathknow)

summary(model4.3.fit)

# Model 4.4.
model4.4.fit <- update(model4.2.fit, fixed = ~ mathkind + sex + minority + ses + housepov)

summary(model4.4.fit)