## Chapter 11 R code.

require(ggplot2)
require(lme4)
require(AICcmodavg)

##

FittedFE <- function(x) model.matrix(x) %*% fixef(x)

lmer.1 <- lmer(read ~ grade5 + (grade5 | subid), MPLS.LS, REML = FALSE)
head(FittedFE(lmer.1))

plotdata <- data.frame(lmer.1@frame, fitted = FittedFE(lmer.1))
plotdata$grade <- plotdata$grade5 + 5
head(plotdata, n = 11)

theme_set(theme_bw())
myx <- scale_x_continuous(breaks=5:8)
g1 <- ggplot(plotdata, aes(x = grade, y = read)) + geom_point(shape = 19)
g2 <- g1 + stat_summary(fun.y = mean, geom = "point", size = 5, shape = 1)
g3 <- g2 + geom_line(aes(y = fitted), lwd = 1.5) + myx + opts("aspect.ratio" = 1)
print(g3)

lmer.2 <- lmer(read ~ grade5 * risk2 + grade5 * eth2 + (grade5 | subid),
               MPLS.LS, REML = FALSE)
plotdata <- data.frame(lmer.2@frame, fitted = FittedFE(lmer.2))
plotdata$grade <- plotdata$grade5 + 5
head(plotdata)

my.m <- ddply(data.frame(plotdata$fitted), .(plotdata$risk2, plotdata$grade), mean)
colnames(my.m) <- c("risk2", "grade", "mfitted")
my.m

g1 <- ggplot(plotdata, aes(x = grade, y = read, shape = risk2))
g2 <- g1 + stat_summary(fun.y = mean, geom = "point")
g3 <- g2 + stat_summary(fun.y = mean, geom = "line", aes(y = fitted, linetype = risk2))
g4 <- g3 + opts("aspect.ratio" = 1) + myx
print(g4)

levels(MPLS.LS$risk)

lmer.3 <- lmer(read ~ grade5 * risk + (grade5 | subid), MPLS.LS, REML = FALSE)
print(lmer.3, cor = FALSE)

rrisk <- relevel(MPLS.LS$risk, ref = "HHM")
lmer.3A <- lmer(read ~ grade5 * rrisk + (grade5 | subid), MPLS.LS, REML = FALSE)
print(lmer.3A, corr = FALSE)

lmer.3AA <- lmer(read ~ grade5 * C(risk, base = 2) + (grade5 | subid),
                 MPLS.LS, REML = FALSE)
summary(lmer.3AA)@coefs

summary(lmer.3)@AICtab
summary(lmer.3A)@AICtab

plotdata <- data.frame(lmer.3@frame, fitted = FittedFE(lmer.3))
plotdata$grade <- plotdata$grade5 + 5

g1 <- ggplot(plotdata, aes(x = grade, y = read, shape = risk))
g2 <- g1 + stat_summary(fun.y = mean, geom = "point")
g3 <- g2 + stat_summary(fun.y = mean, geom = "line", aes(y = fitted, linetype = risk))
g4 <- g3 + myx + opts("aspect.ratio" = 1) + scale_shape(solid = FALSE)
print(g4)

reduced <- lmer(read ~ grade5 + (grade5 | subid), MPLS.LS, REML = FALSE)
full <- lmer(read ~ grade5 * risk + (grade5 | subid), MPLS.LS, REML = FALSE)
anova(reduced, full)

full.1 <- lmer(read ~ grade5 + risk + (grade5 | subid), MPLS.LS, REML = FALSE)
full.2 <- lmer(read ~ grade5 * risk + (grade5 | subid), MPLS.LS, REML = FALSE)
anova(reduced, full.1, full.2)

require(multcomp)
## Create contrast coefficients.
K <- rbind("Int: POV-ADV" = c(0,0,1,0,0,0),
           "Int: HHM-ADV" = c(0,0,0,1,0,0),
           "Int: HHM-POV" = c(0,0,1,-1,0,0),
           "Slo: POV-ADV" = c(0,0,0,0,1,0),
           "Slo: HHM-ADV" = c(0,0,0,0,0,1),
           "Slo: HHM-POV" = c(0,0,0,0,1,-1))
## Save glht() output object.
glht.3 <- glht(lmer.3, K)

summary(glht.3, test = adjusted(type = "none"))

confint(glht.3, calpha = 1.96)

with(MPLS.LS[MPLS.LS$grade == 5, ], table(risk2, eth2))

lmer.4 <- lmer(read ~ grade5 * risk2 * eth2 + (grade5 | subid),
               MPLS.LS, REML = FALSE)
print(lmer.4, cor = FALSE)

lmer.0 <- lmer(read ~ grade5 * eth2 + grade5 * risk2 + (grade5 | subid),
               MPLS.LS, REML = FALSE)
anova(lmer.0, lmer.4)

plotdata <- data.frame(lmer.4@frame, fitted = FittedFE(lmer.4))
plotdata$group4 <- with(plotdata, interaction(eth2, risk2))
plotdata$grade <- plotdata$grade5 + 5
head(plotdata)
levels(plotdata$group4)

g1 <- ggplot(plotdata, aes(x = grade, y = read))
g2 <- g1 + geom_line(aes(groups = subid), colour = "grey80")
g3 <- g2 + stat_summary(fun.y = mean, geom = "point", aes(groups = 1))
g4 <- g3 + stat_summary(fun.y = mean, geom = "line", aes(y = fitted))
g5 <- g4 + facet_grid(eth2 ~ risk2) + myx
print(g5)

g1 <- ggplot(plotdata, aes(x = grade, y = fitted, shape = group4, linetype = group4))
g2 <- g1 + stat_summary(fun.y = mean, geom = "point", size = 2.5) + myx
g3 <- g2 + stat_summary(fun.y = mean, geom = "line") + opts("aspect.ratio" = 1)
g4 <- g3 + scale_shape(solid = FALSE)
print(g4)

## Estimate models.
lmer.1 <- lmer(read ~ grade5 + (grade5 | subid), MPLS.LS, REML = FALSE)
lmer.2 <- lmer(read ~ grade5 * risk2 + (grade5 | subid), MPLS.LS, REML = FALSE)
lmer.3 <- lmer(read ~ grade5 * risk2 + grade5 * eth2 + (grade5 | subid),
               MPLS.LS, REML = FALSE)
## Compute R-squared.
my.rsq <- c(cor(y = lmer.1@frame$read, x = FittedFE(lmer.1)) ^ 2,
            cor(y = lmer.2@frame$read, x = FittedFE(lmer.2)) ^ 2,
            cor(y = lmer.3@frame$read, x = FittedFE(lmer.3)) ^ 2)
data.frame(model = c("lmer.1", "lmer.2", "lmer.3"), rsq = my.rsq)

## aictab().
myaicc <- as.data.frame(aictab(list(lmer.1, lmer.2, lmer.3), sort = FALSE,
                       c("lmer.1", "lmer.2", "lmer.3"))[ ,-c(5,7)])
## Compute R-squared.
myaicc$rsq <- c(cor(y = lmer.1@frame$read, x = FittedFE(lmer.1)) ^ 2,
                cor(y = lmer.2@frame$read, x = FittedFE(lmer.2)) ^ 2,
                cor(y = lmer.3@frame$read, x = FittedFE(lmer.3)) ^ 2)
myaicc

reduced <- lmer(read ~ grade5 + (1 | subid), MPLS.LS, REML = FALSE)
full <- lmer(read ~ grade5 + risk2 + (1 | subid), MPLS.LS, REML = FALSE)
v.r <- attr(VarCorr(reduced)$subid, "stddev")^2
v.f <- attr(VarCorr(full)$subid, "stddev")^2
v.r; v.f

PRV <- (v.r - v.f) / v.r
PRV

reduced <- lmer(read ~ grade5 + risk2 + (grade5 | subid), MPLS.LS, REML = FALSE)
full <- lmer(read ~ grade5 * risk2 + (grade5 | subid), MPLS.LS, REML = FALSE)
v.r <- attr(VarCorr(reduced)$subid, "stddev")[2]^2
v.f <- attr(VarCorr(full)$subid, "stddev")[2]^2
PRV <- (v.r - v.f) / v.r
PRV

reduced <- lmer(read ~ grade5 + (1 | subid), MPLS.LS, REML = FALSE)
full <- lmer(read ~ grade5 * risk2 + (1 | subid), MPLS.LS, REML = FALSE)
## Determinant.
v.r <- det(as.matrix(VarCorr(reduced)$subid))
v.f <- det(as.matrix(VarCorr(full)$subid))
v.r; v.f
## PRV.
PRV <- (v.r - v.f) / v.r
PRV

mywide <- reshape(MPLS.LS, v.names = "read", timevar = "grade",
                  idvar = "subid", direction = "wide", drop = "grade5")
set.seed(1234)       # Use to reproduce this example.
mywide$X1 <- round(rnorm(nrow(mywide), mean = 100, sd = 15), 0)
mylong <- reshape(mywide, direction = "long")
mylong <- mylong[order(mylong$subid),]
head(mylong)

## Estimate the models.
sim.1 <- lmer(read ~ grade + (grade | subid), MPLS.LS, REML = FALSE)
sim.2 <- lmer(read ~ grade * risk2 + (grade | subid), MPLS.LS, REML = FALSE)
sim.3 <- lmer(read ~ grade * risk2 + grade * eth2 + (grade | subid),
              MPLS.LS, REML = FALSE)
sim.4 <- lmer(read ~ grade * risk2 + grade * eth2 + grade * X1 + (grade | subid),
              mylong, REML = FALSE)
## Extract the variances.
myvar <- data.frame(rbind(attr(VarCorr(sim.1)$subid,"stddev") ^ 2,
                          attr(VarCorr(sim.2)$subid,"stddev") ^ 2,
                          attr(VarCorr(sim.3)$subid,"stddev") ^ 2,
                          attr(VarCorr(sim.4)$subid,"stddev") ^ 2))
rownames(myvar) <- c("(none) ", "risk2", "risk2.eth2", "risk2.eth2.X1")
## Compute the determinant of G.
myvar$det <- c(det(as.matrix(VarCorr(sim.1)$subid)),
               det(as.matrix(VarCorr(sim.2)$subid)),
               det(as.matrix(VarCorr(sim.3)$subid)),
               det(as.matrix(VarCorr(sim.4)$subid)))
## Compute R-squared.
myvar$rsq <- c(cor(x = sim.1@frame$read, y = FittedFE(sim.1))^2,
               cor(x = sim.2@frame$read, y = FittedFE(sim.2))^2,
               cor(x = sim.3@frame$read, y = FittedFE(sim.3))^2,
               cor(x = sim.4@frame$read, y = FittedFE(sim.4))^2)
## Deviance.
myvar$deviance <- c(deviance(sim.1), deviance(sim.2), deviance(sim.3), deviance(sim.4))
## AIC.
myvar$AIC <- c(AIC(sim.1), AIC(sim.2), AIC(sim.3), AIC(sim.4))
## Print the results.
colnames(myvar)[1:2] <- c("Var(int)", "Var(slopes)")
myvar

year <- MPLS.LS$grade5
month <- year * 12
day <- year * 360
head(data.frame(year, month, day), n = 4)

## Estimate models.
year.out <- lmer(read ~ year + (year | subid), MPLS.LS, REML = TRUE)
month.out <- lmer(read ~ month + (month | subid), MPLS.LS, REML = TRUE)
day.out <- lmer(read ~ day + (day | subid), MPLS.LS, REML = TRUE)
## Concatenate and print.
my.f <- rbind(summary(year.out)@coefs[2, ],
              summary(month.out)@coefs[2, ],
              summary(day.out)@coefs[2, ])
row.names(my.f) <- c("year", "month", "day")
my.f

my.v <- data.frame(rbind(attr(VarCorr(year.out)$subid,  "stddev") ^ 2,
                         attr(VarCorr(month.out)$subid, "stddev") ^ 2,
                         attr(VarCorr(day.out)$subid,   "stddev") ^ 2),
                   Cor = rbind(attr(VarCorr(year.out)$subid, "correlation")[1,2],
                         attr(VarCorr(month.out)$subid,      "correlation")[1,2],
                         attr(VarCorr(day.out)$subid,        "correlation")[1,2]))
row.names(my.v) <- c("year", "month", "day")
colnames(my.v)[1:2] <- c("Var(int)", "Var(slopes)")
round(my.v, 6)

day.out0 <- lmer(read ~ day + (1 | subid), MPLS.LS)
anova(day.out0, day.out)

MPLS.LS$read.norm <- scale(MPLS.LS$read, center = 206.7, scale = 15.6)
with(MPLS.LS, head(data.frame(subid, grade, read, read.norm)))

norm.1 <- lmer(read.norm ~ grade5 + (grade5 | subid), MPLS.LS, REML = FALSE)
summary(norm.1)@coefs


## MTWO <- read.table("C:/Mine/mytab.txt", header = TRUE) ## Tailor to your system.
head(MTWO)

MTWO.L <- reshape(MTWO, varying = 1:4, direction = "long")
head(MTWO.L)

desc <- ddply(MTWO.L, .(time = MTWO.L$time), function(x)
              c(Mean.a = mean(x$a), SD.a = sd(x$a),
                Mean.b = mean(x$b), SD.b = sd(x$b)))
round(desc, 2)

# Composite then standardize.
MTWO.L$c.z <- scale(MTWO.L$a + MTWO.L$b)
# Standardize then composite.
MTWO.L$z.c <- scale(MTWO.L$a) + scale(MTWO.L$b)
# Means.
desc <- ddply(MTWO.L[ , 5:6], .(time = MTWO.L$time), mean)
round(desc, 2)
# Row differences.
round(data.frame(c.z.diff = diff(desc[, 2]), z.c.diff = diff(desc[ , 3])), 2)