# Import Autism data.
autism <- read.csv("E:\\LMM Data\\autism.csv", h = T)
attach(autism)

# Alternative: read data directly from the web.
file <- "http://www-personal.umich.edu/~bwest/autism.csv"
autism <- read.csv(file, h = T)
attach(autism)

sicdegp.f <- factor(sicdegp)
age.f <- factor(age)

# Add the new variables to the data frame object.
autism.updated <- data.frame(autism, sicdegp.f, age.f)

# Number of Observations at each level of AGE.
summary(age.f)

# Number of Observations at each level of AGE within each group
# defined by the SICDEGP factor.
table(sicdegp.f, age.f)

# Overall summary for VSAE.
summary(vsae)

# VSAE means at each AGE.
tapply(vsae, age.f, mean, na.rm=TRUE)

# VSAE minimum values at each AGE.
tapply(vsae, age.f, min, na.rm=TRUE)

# VSAE maximum values at each AGE.
tapply(vsae, age.f, max, na.rm=TRUE)

# Plots of trends in VSAE as a function of AGE.
library(lattice) # Load the library for trellis graphics.
trellis.device(color=F) # Color is turned off.

# Load the nlme library, which is required for the
# plots as well as for subsequent models.
library(nlme)

# Create grouped data object.
autism.g1 <- groupedData(vsae ~ age | childid,
                           outer = ~ sicdegp.f, data = autism.updated)

# Generate individual profiles in Figure 6.1.
plot(autism.g1, display = "childid", outer = TRUE, aspect = 2,
         key = F, xlab = "Age (Years)", ylab = "VSAE",
         main = "Individual Data by SICD Group")

# Create second grouped data object.
autism.g2 <- groupedData(vsae ~ age | sicdegp.f,
                           order.groups = F, data = autism.updated)

# Generate mean profiles in Figure 6.2.
plot(autism.g2, display = "sicdegp", aspect = 2, key = F,
         xlab = "Age (Years)", ylab = "VSAE",
         main = "Mean Profiles by SICD Group")

# Compute age.2 (AGE minus 2) and age.2sq (AGE2 squared).
age.2 <- age - 2
age.2sq <- age.2*age.2

# Recode the SICDEGP factor for model fitting.
sicdegp2 <- sicdegp
sicdegp2[sicdegp == 3] <- 0
sicdegp2[sicdegp == 2] <- 2
sicdegp2[sicdegp == 1] <- 1
sicdegp2.f <- factor(sicdegp2)

# Omit two records with VSAE = NA, and add the recoded
# variables to the new data frame object.
autism.updated <- subset(data.frame(autism, sicdegp2.f, age.2),
                             !is.na(vsae))

# Alternative.
sicdegp2 <- cut(sicdegp, breaks = 0:3, labels= FALSE)

# MODEL FITTING.

# re-create grouped data object.
autism.grouped <- groupedData(vsae ~ age.2 | childid,
                              data = autism.updated, order.groups = F)

# Fit Model 6.1.
model6.1.fit <- lme(vsae ~ age.2 + I(age.2^2) + sicdegp2.f +
                        age.2:sicdegp2.f + I(age.2^2):sicdegp2.f,
                      random = ~ age.2 + I(age.2^2), method= "REML",
                      data = autism.grouped)

# Fit Model 6.2.
model6.2.fit <- lme(vsae ~ age.2 + I(age.2^2) + sicdegp2.f +
                        age.2:sicdegp2.f + I(age.2^2):sicdegp2.f,
                      random = ~ age.2 + I(age.2^2) - 1, method= "REML",
                      data = autism.grouped)
summary(model6.2.fit)

# Fit Model 6.2a.
model6.2a.fit <- update(model6.2.fit, random = ~ age.2 - 1)

# Test Hypothesis 6.1.
h6.1.pvalue <- 0.5*(1-pchisq(83.9,1)) + 0.5*(1-pchisq(83.9,2))
h6.1.pvalue

# Test Hypothesis 6.2.
model6.2.ml.fit <- update(model6.2.fit, method = "ML")
model6.3.ml.fit <- update(model6.2.ml.fit,
                          fixed = ~ age.2 + I(age.2^2) + sicdegp2.f + age.2:sicdegp2.f)
anova(model6.2.ml.fit, model6.3.ml.fit)

# Fit Model 6.3.
model6.3.fit <- update(model6.2.fit,
                       fixed = ~ age.2 + I(age.2^2) + sicdegp2.f + age.2:sicdegp2.f)

# Examine results of final model.
summary(model6.3.fit)
intervals(model6.3.fit)
getVarCov(model6.3.fit, individual= "1", type= "marginal")

curve(0.11*x^2 + 6.15*x + 13.46, 0, 11,
        xlab = "AGE minus 2", ylab = "Marginal Predicted VSAE",
        lty = 3, ylim=c(0,100), lwd = 2)

curve(0.11*x^2+ 2.65*x + 9.84, 0, 11, add=T, lty = 2, lwd = 2)

curve(0.11*x^2 + 2.08*x + 8.47, 0, 11, add=T, lty = 1, lwd = 2)

# Add a legend to the plot; R will wait for the user to click
# on the point in the plot where the legend is desired.
  
legend(locator(1),
           c("SICDEGP = 1", "SICDEGP = 2", "SICDEGP = 3"),
    lty = c(1, 2, 3), lwd = c(2, 2, 2))

plot(augPred(model6.3.fit, level = 0:1),
     layout = c(4, 3, 1), xlab = "AGE minus 2", ylab = "Predicted VSAE",
     key = list(
       lines = list(lty = c(1, 2), col = c(1, 1), lwd = c(1, 1) ),
       text = list(c("marginal mean profile", "subject-specific profile")),
       columns = 2))

# MODEL DIAGNOSTICS

# Residual-fitted plots.
plot(model6.3.fit,
     resid(., type = "p") ~ fitted(.) | factor(sicdegp),
     layout=c(3,1), aspect=2, abline=0)

plot(model6.3.fit,
     resid(., type= "p") ~ fitted(.) | factor(sicdegp),
     id = 0.05, layout = c(3,1), aspect = 2, abline = 0)

plot(model6.3.fit, resid(.) ~ age.2, abline = 0)

# Assess assumptions of normality.
qqnorm(model6.3.fit,
       ~resid(.) | factor(sicdegp) ,
       layout = c(3,1), aspect = 2, id = 0.05)

# Check random effects.
qqnorm(model6.3.fit, ~ranef(.) , id = 0.10)

pairs(model6.3.fit,
      ~ranef(.) | factor(sicdegp),
      id = ~childid == 124, layout = c(3, 1), aspect = 2)

# Fitted values against predicted values.
plot(model6.3.fit, vsae ~ fitted(.) | factor(sicdegp),
     id  = 0.05, layout = c(3,1) , aspect = 2)

# See if results are sensitive to outlier(s).
autism.grouped2 <- autism.grouped[(autism.grouped$childid != 124 & autism.grouped$childid != 46),]

model6.3.fit.out <- update(model6.3.fit, data = autism.grouped2)


# ALTERNATIVE MARGINAL MODELING APPROACH

index <- age.2
index[age.2 ==  0] <- 1
index[age.2 ==  1] <- 2
index[age.2 ==  3] <- 3
index[age.2 ==  7] <- 4
index[age.2 == 11] <- 5

autism.updated <- subset(data.frame(
  autism,sicdegp2.f, age.2, index), !is.na(vsae))

marg.model.fit <- gls(
  vsae ~ age.2 + I(age.2^2) + sicdegp2.f + age.2:sicdegp2.f,
  correlation=corSymm(form = ~ index | childid),
  weights = varIdent(form = ~ 1 | age.2) , autism.updated)

summary(marg.model.fit)

getVarCov(marg.model.fit)