# Load lme4 library.
library(lme4)

# Fit Model 7.1 (see chapter7_R_final.R for code to load the data).

# Need to create a new version of the tooth variable for nesting.
veneer$tooth2 <- as.numeric(paste(factor(veneer$patient),
                                  factor(veneer$tooth),sep=""))

model7.1.fit.lmer <- lmer(gcf ~ time + base_gcf + cda + age +
                            time*base_gcf + time*cda + time*age +
                            (time | patient) + (1 | tooth2),
                          data = veneer, REML = T)

summary(model7.1.fit.lmer)
confint(model7.1.fit.lmer)

# EBLUPs
ranef(model7.1.fit.lmer)$patient
ranef(model7.1.fit.lmer)$tooth2

# Fit Model 7.1a.
model7.1A.fit.lmer <- lmer(gcf ~ time + base_gcf + cda + age +
                             time*base_gcf + time*cda + time*age +
                             (time | patient),
                           data = veneer, REML = T)

# Test Hypothesis 7.1.
anova(model7.1.fit.lmer, model7.1A.fit.lmer)

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

# Test Hypothesis 7.3.
model7.1.ml.fit.lmer <- lmer(gcf ~ time + base_gcf + cda + age +
                               time*base_gcf + time*cda + time*age +
                               (time | patient) + (1 | tooth2),
                             data = veneer, REML = F)

model7.3.ml.fit.lmer <- lmer(gcf ~ time + base_gcf + cda + age +
                               (time | patient) + (1 | tooth2),
                             data = veneer, REML = F)

anova(model7.1.ml.fit.lmer, model7.3.ml.fit.lmer)

# Fit Model 7.3 using REML.
model7.3.fit.lmer <- lmer(gcf ~ time + base_gcf + cda + age +
                            (time | patient) + (1 | tooth2),
                          data = veneer, REML = T)

summary(model7.3.fit.lmer)