# Import Dental Veneer data.
veneer <- read.delim("E:\\LMM Data\\veneer.dat", h = T)

library(nlme)

# Fit Model 7.1.
model7.1.fit <- lme(gcf ~ time + base_gcf + cda + age +
                      time:base_gcf + time:cda + time:age,
                    random = list(patient = ~time, tooth = ~1),
                    data = veneer, method = "REML")

summary(model7.1.fit)
intervals(model7.1.fit)
random.effects(model7.1.fit)

# Fit Model 7.1a.
model7.1A.fit <- lme(gcf ~ time + base_gcf + cda + age +
                       time:base_gcf + time:cda + time:age,
                     random = list(patient = ~ time),
                     data = veneer, method = "REML")

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

# Fit Model 7.2a.
model7.2A.fit <- lme(gcf ~ time + base_gcf + cda + age +
                       time:base_gcf + time:cda + time:age,
                     random = list(patient = ~ time, tooth = ~1),
                     corr = corCompSymm(0.5, form = ~1 | patient/tooth),
                     weights = varIdent(form = ~1 | time),
                     data = veneer, method = "REML")

summary(model7.2A.fit)
intervals(model7.2A.fit) # Note unstable confidence intervals.

# Fit Model 7.2b.
model7.2B.fit <- lme(gcf ~ time + base_gcf + cda + age +
                       time:base_gcf + time:cda + time:age,
                     random = list(patient = ~ time, tooth = ~1),
                     corr = corCompSymm(0.5, form = ~1 | patient/tooth),
                     data = veneer, method = "REML")

summary(model7.2B.fit)
intervals(model7.2B.fit) # Note unstable confidence intervals.

# Fit Model 7.2c.
model7.2C.fit <- lme(gcf ~ time + base_gcf + cda + age +
                       time:base_gcf + time:cda + time:age,
                     random = list(patient = ~ time, tooth = ~1),
                     weights = varIdent(form = ~1 | time),
                     data = veneer, method = "REML")

summary(model7.2C.fit)
intervals(model7.2C.fit)

# Test Hypothesis 7.2.
anova(model7.2C.fit, model7.1.fit)

# Fit Model 7.1 using ML to test fixed effects.
model7.1.ml.fit <- lme(gcf ~ time + base_gcf + cda + age +
                         time:base_gcf + time:cda + time:age,
                       random = list(patient = ~ time, tooth = ~1),
                       data = veneer, method = "ML")

# Fit Model 7.3 using ML.
model7.3.ml.fit <- lme(gcf ~ time + base_gcf + cda + age,
                       random = list(patient = ~ time, tooth = ~1),
                       data = veneer, method = "ML")

# Test Hypothesis 7.3.
anova(model7.1.ml.fit, model7.3.ml.fit)

# Fit Model 7.3 using REML estimation.
model7.3.fit <- lme(gcf ~ time + base_gcf + cda + age,
                    random = list(patient = ~ time, tooth = ~1),
                    data = veneer, method = "REML")

summary(model7.3.fit)