###################################################
### chunk number 1: pg 388
###################################################
driscoll1 <- driscoll
driscoll1[9,4] <- NA


###################################################
### chunk number 2: pg 388
###################################################
driscoll1.aov <- aov(CALLS~Error(BLOCK)+YEAR, data=driscoll1, subset=BLOCK!="newpipe")
summary(driscoll1.aov)


###################################################
### chunk number 3: pg 389
###################################################
#calculate the mean of the newpipe block
BM<-with(driscoll1, tapply(CALLS, BLOCK, mean, na.rm=T))["newpipe"]
#calculate the mean of year 2
YM<-with(driscoll1, tapply(CALLS, YEAR, mean, na.rm=T))["2"]
#calculate the overall mean
M<-mean(driscoll1$CALLS,na.rm=T)
#duplicate the data set and work on the duplicate
driscoll2 <- driscoll1
#substitute the new value into the data frame
driscoll2[9,3]<-YM+BM-M
#fit the linear model
driscoll2.aov <- aov(CALLS~Error(BLOCK)+YEAR, data=driscoll2)
summary(driscoll2.aov)
#then make adjustments to the F-ratio and Pvalue (to reflect a reduction) 
#in residual degrees of freedom by one for each substituted value)
(MSresid <- summary(driscoll2.aov)[[2]][[1]]["Residuals","Sum Sq"]/9)
(Fyear <- summary(driscoll2.aov)[[2]][[1]]["YEAR","Mean Sq"]/MSresid)
(Pvalue <- 1-pf(Fyear, 2,8))


###################################################
### chunk number 4: pg 389
###################################################
#fit full model
driscoll1.aovF <- aov(CALLS~BLOCK+YEAR, data=driscoll1)
#fit reduced model
driscoll1.aovR <- aov(CALLS~BLOCK, data=driscoll1)
anova(driscoll1.aovF,driscoll1.aovR)


###################################################
### chunk number 5: pg 390
###################################################
anova(driscoll1.aovF)


###################################################
### chunk number 6: pg 390
###################################################
library(nlme)
#No structure
driscoll1.lme1 <- lme(CALLS ~ YEAR, random = ~1 | BLOCK, data = driscoll1, subset = !is.na(CALLS))
#Unstructured 
driscoll1.lme2 <- lme(CALLS ~ YEAR, random = ~1 | BLOCK, data = driscoll1, subset = !is.na(CALLS), correlation = corSymm(form = ~1 | BLOCK))
#Compound symmetry
driscoll1.lme3 <- update(driscoll1.lme1, correlation = corCompSymm(form = ~1 | BLOCK))
#First order autoregressive
driscoll1.lme4 <- lme(CALLS ~ YEAR, random = ~1 | BLOCK, data = driscoll1, subset = !is.na(CALLS), correlation = corAR1(form = ~1 | BLOCK))
driscoll1.lme4 <- update(driscoll1.lme1, correlation = corAR1(form = ~1 | BLOCK))
#Compare each to compound symmetry
#technically,only models fitted with ML (not REML)
#should be compared via anova
driscoll1.lmeML1<-update(driscoll1.lme1, method="ML")
driscoll1.lmeML2<-update(driscoll1.lme2, method="ML")
driscoll1.lmeML3<-update(driscoll1.lme3, method="ML")
driscoll1.lmeML4<-update(driscoll1.lme4, method="ML")
anova(driscoll1.lmeML3, driscoll1.lmeML1, driscoll1.lmeML2, driscoll1.lmeML4)
anova(driscoll1.lme3)