# Load lme4 library.
library(lme4)

# Fit Model 4.1 (see chapter4_R_final.R for code to load the data).
model4.1.fit.lmer <- lmer(mathgain ~ 1 + (1|schoolid) + (1|classid),
                            class, REML = T)

# View results and 95% confidence intervals.
summary(model4.1.fit.lmer)
confint(model4.1.fit.lmer)

# Display the random effects (EBLUPs) from the model.
ranef(model4.1.fit.lmer)

# Plot the predicted random effects along with measures of uncertainty.
library(merTools)
REsim(model4.1.fit.lmer)
plotREsim(REsim(model4.1.fit.lmer))

# Model 4.1A.
model4.1A.fit.lmer <- lmer(mathgain ~ 1 + (1|schoolid),
                               class, REML = T)

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

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

summary(model4.2.fit.lmer)

# Model 4.1: ML estimation with lmer().
model4.1.lmer.ml.fit <- lmer(mathgain ~ 1 + (1|schoolid) + (1|classid),
                                 class, REML = F)

# Model 4.2: ML estimation with lmer().
model4.2.lmer.ml.fit <- lmer(mathgain ~ mathkind + sex + minority + ses
                                 + (1|schoolid) + (1|classid),
                                 class, REML = F)

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

# Model 4.3.
model4.3.fit.lmer <- lmer(mathgain ~ mathkind + sex + minority + ses
                              + yearstea + mathprep + mathknow
                              + (1|schoolid) + (1|classid),
                              class, na.action = "na.omit", REML = T)

summary(model4.3.fit.lmer)

# Model 4.4.
model4.4.fit.lmer <- lmer(mathgain ~ mathkind + sex + minority + ses
                              + housepov + (1|schoolid) + (1|classid),
                              class, na.action = "na.omit", REML = T)

summary(model4.4.fit.lmer)