# Load the Rat Pup data.
ratpup <- read.table("E:\\LMM Data\\rat_pup.dat", h = T)
attach(ratpup)

ratpup$sex1[sex == "Female"] <- 1
ratpup$sex1[sex == "Male"] <- 0

# Load the nlme library for use of lme().
library(nlme)

# Fit Model 3.1.
model3.1.fit <- lme(weight ~ treatment + sex1 + litsize +
                          treatment:sex1, random = ~ 1 | litter,
                        data = ratpup, method = "REML")

# Optional changing of reference category.
treatment <- relevel(treatment, ref = "High")

# Summary of model fit, along with F-tests.
summary(model3.1.fit)
anova(model3.1.fit)

# Display the random effects (EBLUPs) from the model.
random.effects(model3.1.fit)

# Fit Model 3.1A.
model3.la.fit <- gls(weight ~ treatment + sex1 + litsize +
                           treatment:sex1, data = ratpup)
anova(model3.1.fit, model3.la.fit) # Test Hypothesis 3.1.

# Fit Model 3.2A.
model3.2a.fit <- lme(weight ~ treatment + sex1 + litsize
                         + treatment:sex1, random = ~1 | litter, ratpup, method = "REML",
                         weights  = varIdent(form = ~1 | treatment))

summary(model3.2a.fit)

# Test Hypothesis 3.2.
anova(model3.1.fit, model3.2a.fit)

ratpup$trtgrp[treatment == "Control"] <- 1
ratpup$trtgrp[treatment == "Low" | treatment == "High"] <- 2

model3.2b.fit <- lme(weight ~ treatment + sex1 + litsize + treatment:sex1,
                       random  = ~ 1 | litter, ratpup, method = "REML",
                       weights = varIdent(form = ~1 | trtgrp))

# Test Hypothesis 3.3.
anova(model3.2a.fit, model3.2b.fit)

# Test Hypothesis 3.4.
anova(model3.1.fit, model3.2b.fit)

summary(model3.2b.fit)

# Test Hypothesis 3.5.
anova(model3.2b.fit)

model3.3.ml.fit <- lme(weight ~ treatment + sex1 + litsize,
                         random = ~1 | litter, ratpup, method = "ML",
                         weights = varIdent(form = ~1 | trtgrp))

model3.3a.ml.fit <- lme(weight ~ sex1 + litsize,
                          random = ~1 | litter, ratpup, method = "ML",
                          weights = varIdent(form = ~1 | trtgrp))

# Test Hypothesis 3.6.
anova(model3.3.ml.fit, model3.3a.ml.fit)

# Model 3.3: Final Model.
model3.3.reml.fit <- lme(weight ~ sex1 + litsize + treatment,
                             random  = ~1 | litter, ratpup, method = "REML",
                             weights = varIdent(form = ~1 | trtgrp))
summary(model3.3.reml.fit)
anova(model3.3.reml.fit)