###############################################################################
# Two-factor Mixed Model, A fixed, B random: Sunscreen data, Ott Ex 17.3
###############################################################################
library(nlme) # needed for lme function 
### Create a grouped data object for Trellis plotting
###  - we see a lot of variability among tech's and less for sun w/in tech.
Sunscreen=groupedData(burn ~ sun | tech/sun, data=read.csv("Data/Ott5thEd-Data/Chap17/ch17-Example17.3.csv"))
pdf(file="Plots/sun1.pdf",pointsize=9,width=6,height=4.5)
#sun.plots=plot.design(Sunscreen)
plot(Sunscreen, aspect=3)
dev.off()
### Read data for modeling
sunscreen <- read.csv("Data/Ott5thEd-Data/Chap17/ch17-Example17.3.csv")
# first convert numbers to factor variables
sunscreen$sun <- as.factor(sunscreen$sun)
sunscreen$tech <- as.factor(sunscreen$tech)
# model with both effects fixed
sun.lm <- lm(burn~sun+tech+sun*tech, data=sunscreen)
# fit mixed model with: sun (fixed), tech (random), sun*tech (random)
sun.lme <- lme(burn ~ sun, data=sunscreen, random=~1 | tech/sun, method="REML")
# mixed model with no interaction: sun (fixed), tech (random):
#lme(burn ~ sun, data=sunscreen, random=~1 | tech, method="REML")
summary(sun.lme)
# test which model is better fit (NB: must use method="ML" in sun.lme, default is REML)
#anova(sun.lme,sun.lm) # sun.lm is better fit
anova(sun.lme)
# get 95% CI's for the variance components: 
#ranef(sun.lme)
intervals(sun.lme, which="var-cov")
### Resids & fitted diagnostics
# check normality for the 3 variance components
#pdf(file="Plots/sun4.pdf",pointsize=3,width=6,height=5)
qqnorm(sun.lme, ~resid(.))          # sigma_epsilon
qqnorm(sun.lme, ~ranef(., level=1), main="Tech") # sigma_tech
qqnorm(sun.lme, ~ranef(., level=2), main="Sun*Tech") # sigma_(tech*sun)
# check constant variance of sigma_epsilon
plot(sun.lme, resid(.)~fitted(.), abline=0)
# ideally use a pairs plot to assess homoscedasticity of ranef 1 & 2, but it's not working...
pairs(sun.lme, ~ranef(.))
#plot(sun.lme, ranef(., level=1)~fitted(.), abline=0)
#plot(sun.lme, ranef(., level=2)~fitted(.), abline=0)
#dev.off()



###############################################################################
# Two-factor Random Effects Model, B nested in A: Drug data, Ott Ex 17.10
# Note: this ex identical to Section 4.2.3 of Pinheiro & Bates (Thickness of Oxide)
###############################################################################
library(nlme) # needed for lme function 
content <- read.csv("Data/Ott5thEd-Data/Chap17/ch17-Example17.10.csv")
# first convert numbers to factor variables
content$Site <- as.factor(content$Site)
content$Batch <- as.factor(content$Batch)
# fit random effects model with Batch nested in Site: Batch(Site)
drug.lme <- lme(Drug~1, data=content, random=~1 | Site/Batch, method="ML")
# get estimates of variance components
summary(drug.lme)
# get 95% CI's for the variance components: fails here because of
#   non-positive definite approximate variance-covariance matrix...
intervals(drug.lme, which="var-cov")
# test if Batch random effect is sig. (p=9e-04 => sig. Batch effect) 
drug1.lme <- lme(Drug~1, data=content, random=~1 | Site/Batch, method="ML")
drug2.lme <- lme(Drug~1, data=content, random=~1 | Site, method="ML")
anova(drug1.lme, drug2.lme)
# get mean Y values for sites (level=1), use level=2 for batches w/in sites
coef(drug.lme, level=1)
#
qqnorm(drug.lme, ~resid(.))
qqnorm(drug.lme, ~ranef(., level=1))


###############################################################################
# CRD split-plot design: Soybean yield, Ott 5th Ed Ex 17.11 (but ignore blocks)
###############################################################################
library(nlme) # needed for lme function 
soy <- read.csv("Data/Ott5thEd-Data/Chap17/ch17-Example17.11.csv")
### First convert numbers to factor variables
soy$WPlot <- as.factor(soy$WPlot)
soy$Fertilizer <- as.factor(soy$Fertilizer)
soy$Variety <- as.factor(soy$Variety)
### Split-plot model with WPlot nested in Fertilizer (cannot get test for Fert...) 
#soy.lme <- lme(Yield~Fertilizer*Variety, data=soy, random=~1 | Fertilizer/WPlot)
### It suffices to specify groups as WPlot (will know nested in Fertilizer) 
soy.lme <- lme(Yield~Fertilizer*Variety, data=soy, random=~1 | WPlot)
summary(soy.lme)
anova(soy.lme)
# the same thing using aov 
soy.lm <- aov(Yield~Fertilizer*Variety+Error(WPlot), data=soy)
summary(soy.lm)
### get 95% CI's for the variance components: 
intervals(soy.lme, which="var-cov")
### table of LS means
model.tables(soy.lm, type="means")


###############################################################################
# Block split-plot design: Soybean yield, Ott 5th Ed Ex 17.11 (with blocks)
###############################################################################
library(nlme) # needed for lme function 
pdf(file="Plots/soy.pdf",pointsize=9,width=6,height=4.5)
SoyData=groupedData(Yield ~ Variety | Block, data=read.table("Data/Ott5thEd-Data/Chap17/ch17-Example17.11.txt", header=TRUE))
plot(SoyData, inner=~Fertilizer, aspect=1, main="Trellis plot of the Soy Data") 
dev.off()
#
soy <- read.table("Data/Ott5thEd-Data/Chap17/ch17-Example17.11.txt", header=TRUE)
### First convert numbers to factor variables
soy$Block <- as.factor(soy$Block)
soy$WPlot <- as.factor(soy$WPlot)
soy$Fertilizer <- as.factor(soy$Fertilizer)
soy$Variety <- as.factor(soy$Variety)
# Neither Variety nor Fertilizer are sig. under this model
#soy.lme <- lme(Yield~Fertilizer*Variety+Block+Fertilizer:Block, random=~1 | Block, data=soy)
soy.lme <- lme(Yield~Fertilizer*Variety, random=~1 | Block/Fertilizer, data=soy)
anova(soy.lme)
summary(soy.lme)
# Could refine model by omitting least sig factor: Variety, etc.
soy2.lme <- lme(Yield~Variety, random=~1 | Block/Fertilizer, data=soy)
anova(soy2.lme)
#plot( augPred(SoyData) )


###############################################################################
# Repeated Measures Design: root growth under fert/no, Crawley's R Book, p.641
###############################################################################
require(nlme)     # needed for lme function 
require(lattice)  # for trellis plotting
require(multcomp) # for Tukey comparisons
grow <- read.table("Data/Crawley-R-Book/fertilizer.txt", header=TRUE)
# Trellis panel plots of root length of each plant over weeks
grow.group <- groupedData(root~week|plant, outer=~fertilizer, data=grow)
plot(grow.group)         
# Now plot plants over weeks, but aggregate over the treatment (fert added/control)
plot(grow.group, outer=T)
# Pairwise scatterplots of plant growths over weeks to glean within plant corr structure
# Create custom function panel.cor to put corrleation coeffts on upper panel
panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...)
     {
         usr <- par("usr"); on.exit(par(usr))
         par(usr = c(0, 1, 0, 1))
         r <- cor(x, y)
         txt <- format(c(r, 0.123456789), digits=digits)[1]
         txt <- paste(prefix, txt, sep="")
         if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
         text(0.5, 0.5, txt, cex = cex.cor * r)
     }
pdf(file="Plots/root1.pdf",pointsize=9,width=6,height=4.5)
pairs(~ root[week==2] + root[week==4] + root[week==6] + root[week==8] + root[week==10], data=grow, main="Pairwise scatter & correlation of root growths over weeks", lower.panel=panel.cor)
dev.off()
#
# Model with plant (random) nested in fertilizer, fertilizer & week fixed.
#   Assumes IID within plant corr structure by default
#   First convert week to factor (categorical) variable
grow$week <- as.factor(grow$week)
grow.lme <- lme(root~fertilizer*week, data=grow, random=~week | plant)
# Update grow.lme by allowing for compound symmetry corr w/in plants
#   Don't need to do this, will just increase AIC by 2 because of extra param (rho),
#   the corr structure is already captured by the implied structure of grow.lme.
#grow.lme.2 <- lme(root~fertilizer*week, data=grow, random=~1 | plant,
#                correlation=corCompSymm())
# Update grow.lme by allowing for ar(1) corr w/in plants 
grow.lme.3 <- lme(root~fertilizer*week, data=grow, random=~1 | plant,
                correlation=corAR1())
#
summary(grow.lme)
anova(grow.lme)
summary(grow.lme.2)
summary(grow.lme.3)
# Compare two models: p=0.8912 means that simpler (implied compound symmetry) lme1
#  suffices (also has lower AIC). Note: must fit via ML (REML is default).
grow.lme1 <- lme(root~fertilizer*week, data=grow, random=~1 | plant, method="ML")
grow.lme2 <- lme(root~fertilizer*week, data=grow, random=~1 | plant, method="ML",
                 correlation=corAR1())
anova(grow.lme1,grow.lme2)
# Now fit this 2-way anova via AOV just to extract the LS means
grow.lm <- aov(root~fertilizer*week+Error(plant), data=grow)
#summary(grow.lm)
model.tables(grow.lm, type="means")
# Tukey multiple comparisons (cannot get fertilizer*week comps to work...),
#   this will compare weeks by averaging over fertilizer (inappropriate)
summary(glht(grow.lme, linfct=mcp(week="Tukey"))) 
# Profile plot over weeks (fertilizer added vs. control)
week=c(2,     4,     6,     8,    10)
added.means=exp(c(1.667, 3.683, 5.972, 7.450, 9.617))
control.means=exp(c(1.250, 2.250, 4.100, 5.917, 8.333))
plot(week, added.means, type="n", ylab="group means", xlab="week", main="Profile plot for root length")
points(week, added.means, pch=19)
lines(week, added.means, lty=2)
points(week, control.means, pch=1)
lines(week, control.means, lty=3)
legend(x="topleft",legend=c("control","fertilizer added"), pch=c(1,19))
# Finally, diagnostic plots:
#  resids & fitted from "plant" correspond to random effect for plants (beta_j(i))
#  resids & fitted from "fixed" correspond to residual random effect (epslion_ijk)
#  Allow us to assess plausibility of assumed normality & constant variance.
fit.eps=grow.lme.3$fitted[,1]  # fitted for fixed (epsilon)
res.eps=grow.lme.3$resid[,1]   # resids for fixed (epsilon)
fit.bet=grow.lme.3$fitted[,2]  # fitted for plant (beta)
res.bet=grow.lme.3$resid[,2]   # resids for plant (beta)
#  Uncomment the "pdf" & "dev.off" commands to save to a pdf file 
#pdf(file="a.pdf",pointsize=9,width=10,height=8)
par(mfrow=c(2,2))
plot(fit.eps,res.eps, xlab="fitted", ylab="resids", main="Epsilon Residuals"); abline(0,0); qqnorm(res.eps);
plot(fit.bet,res.bet, xlab="fitted", ylab="resids", main="Beta Residuals (plant random effect)"); abline(0,0); qqnorm(res.bet);
#dev.off()