# ---------------------------------------------------------- 
# Exercise 9.1 a,c,d,e                                       
# ---------------------------------------------------------- 

data <- scan(what = list(run=0,temp=0,germrate=0),multi.line=F)
 1 25 24.65  1 40 1.34  
 2 30 24.38  2 40 2.24 
 3 25 29.17  3 30 21.25 
 4 35 5.90   4 40 1.83 
 5 25 28.90  5 35 18.27 
 6 30 25.53  6 35 8.42 

data

prob9.1 <- data.frame(data)


prob9.1 <- within(prob9.1, 
   { run  <- factor(run)
     temp <- factor(temp) } )


prob9.1.rcb  <- lm(germrate ~ run +temp, data=prob9.1)
anova(prob9.1.rcb)

 # residual plots

par(mfrow=c(2,2)) 
plot(prob9.1.rcb)
par(mfrow=c(1,1)) 


 # we can use the lsmeans package to get standard errors of means

library(lsmeans)
lsmeans(prob9.1.rcb,"temp")


 # by using the block as a random effect, we can use the lme4 package
 # along with the lsmeans package to obtain standard errors from an interblock analysis

library(lme4)
library(lsmeans)

prob9.1.lmer1  <- lmer(germrate ~ temp +(1|run), data=prob9.1)
lsmeans(prob9.1.lmer1, "temp")


# ---------------------------------------------------------- 
# Exercise 14.6                                              
# ---------------------------------------------------------- 


data <- scan(what = list(container=0,temp=0,seafood=0,logcount=0),multi.line=F)
1 0 1 3.6882
1 0 2 0.3565
2 0 1 1.8275
2 0 2 1.7023
3 0 1 5.2327
3 0 2 4.5780
4 5 1 7.1950
4 5 2 5.0169
5 5 1 9.3224
5 5 2 7.9519
6 5 1 7.4195
6 5 2 6.3861
7 10 1 9.7842
7 10 2 10.1352
8 10 1 6.4703
8 10 2 5.0482
9 10 1 9.4442
9 10 2 11.0329

data

prob14.6 <- data.frame(data)

prob14.6  <- within(prob14.6 , 
   { container  <- factor(container)
     temp <- factor(temp)
     seafood <- factor(seafood)   } )


 # we can check the ANOVA table 

anova(lm(logcount ~ temp +temp:container +seafood +seafood:temp, data = prob14.6))


 # now using lme4 package


library(lme4)
library(lmerTest)

prob14.6.lme1  <- lmer(logcount ~ temp +(1|temp:container) +seafood +seafood:temp, data = prob14.6)

summary(prob14.6.lme1)

anova(prob14.6.lme1)

lsmeans(prob14.6.lme1)

difflsmeans(prob14.6.lme1)

 # residual plots

par(mfrow=c(2,1)) 
plot(fitted(prob14.6.lme1),resid(prob14.6.lme1),xlab="Predicted Values",
     ylab="Residuals",main="Residual by Predicted plot")
qqnorm(resid(prob14.6.lme1))
qqline(resid(prob14.6.lme1))
par(mfrow=c(1,1)) 


 # the mixed model can also be fitted in afex

library(afex)

prob14.6.afex1  <- mixed(logcount ~ temp +(1|temp:container) +seafood +seafood:temp, data = prob14.6)

summary(prob14.6.afex1)

nice(prob14.6.afex1)


# ---------------------------------------------------------- 
# Exercise 15.1                                              
# ---------------------------------------------------------- 

data <- scan(what = list(diet=0,subject=0,time=0,glucose=0), multi.line=F)
 1 1 15 28  1 1 30 34  1 1 45 32 
 1 2 15 15  1 2 30 29  1 2 45 27 
 1 3 15 12  1 3 30 33  1 3 45 28 
 1 4 15 21  1 4 30 44  1 4 45 39 
 2 5 15 22  2 5 30 18  2 5 45 12 
 2 6 15 23  2 6 30 22  2 6 45 10 
 2 7 15 18  2 7 30 16  2 7 45  9 
 2 8 15 25  2 8 30 24  2 8 45 15 
 3 9 15 31  3 9 30 30  3 9 45 39 
 3 10 15 28  3 10 30 27 3 10 45 36 
 3 11 15 24  3 11 30 26 3 11 45 36 
 3 12 15 21  3 12 30 26 3 12 45 32 

data

prob15.1 <- data.frame(data)

 # using the within function to convert variables to factor
prob15.1 <- within(prob15.1, 
   { diet  <- factor(diet)
     subject <- factor(subject)
     time       <- factor(time) } )


 # get an interaction plot (profile plot)

with(prob15.1,interaction.plot(diet,time,glucose))


 # first check the ANOVA table 

anova(lm(glucose ~ diet +subject:diet +time +diet:time, data = prob15.1))



 # -----------------------------------------------------------------------
 # Split Plot Analysis - gives mostly the same information as Proc Mixed 
 # -----------------------------------------------------------------------

library(lme4)
library(lmerTest)

prob15.1.lme1  <-  lmer(glucose ~ diet +(1|subject:diet) +time +diet:time, data = prob15.1)

anova(prob15.1.lme1)

lsmeans(prob15.1.lme1)
difflsmeans(prob15.1.lme1)


 # residual plots

par(mfrow=c(2,1)) 
plot(fitted(prob15.1.lme1),resid(prob15.1.lme1),xlab="Predicted Values",
     ylab="Residuals",main="Residual by Predicted plot")
qqnorm(resid(prob15.1.lme1))
qqline(resid(prob15.1.lme1))
par(mfrow=c(1,1)) 



 # or with afex

library(afex)

prob15.1.afex1  <-  mixed(glucose ~ diet +(1|subject:diet) +time +diet:time, data = prob15.1)


 # -----------------------------------------------------------------------
 # Mauchly tests and adjusted p values 
 # -----------------------------------------------------------------------

 # using the aov_ez function from the afex package

prob15.1.afex2 <- aov_ez(id="subject",dv="glucose",data=prob15.1,between=c("diet"),within=c("time"))

summary(prob15.1.afex2)

nice(prob15.1.afex2)