# ---------------------------------------------------------------
 # Problem 12.2
 # ---------------------------------------------------------------


#  first read the data directly from Dr. Oehlerts website

Prb12.2 <- read.table("http://www.stat.umn.edu/~gary/book/fcdae.data/pr12.2",header=TRUE)

Prb12.2$alloy   <- as.factor(Prb12.2$alloy)
Prb12.2$method  <- as.factor(Prb12.2$method)
Prb12.2$dentist  <- as.factor(Prb12.2$dentist)

 # get the ANOVA table

Prb12.2.lm1 <- lm(hardness ~ (alloy +method +dentist)^2, data = Prb12.2)

anova(Prb12.2.lm1)

 # examine the interaction

interaction.plot(Prb12.2$method,Prb12.2$dentist,
                 Prb12.2$hardness,type="b")


 # check model assumptions (assuming fixed effects, but still useful)

par(mfrow=c(2,2))  
plot(Prb12.2.lm1)
par(mfrow=c(1,1)) 


 # The afex library calls lmer and other functions to give results fairly 
 # similar to those in Proc MIXED

library(afex)

library('emmeans')   # needed for more recent versions of afex

Prb12.2.afex1 <- mixed(hardness ~ alloy +method +alloy:method +(1|dentist) 
                                  +(1|dentist:alloy) +(1|dentist:method), data = Prb12.2)

summary(Prb12.2.afex1)

nice(Prb12.2.afex1)

lsmeans(Prb12.2.afex1, ~method)

 # examine potential interactions

lsmip(Prb12.2.afex1, method ~ alloy)





 # ---------------------------------------------------------------
 # Exercise 13.2
 # ---------------------------------------------------------------

#  first read the data directly from Dr. Oehlerts website

Ex13.2 <- read.table("http://www.stat.umn.edu/~gary/book/fcdae.data/ex13.2",header=TRUE)

Ex13.2

Ex13.2$batch  <- as.factor(Ex13.2$batch)
Ex13.2$temp  <- as.factor(Ex13.2$temp)
Ex13.2$pH    <- as.factor(Ex13.2$pH)
Ex13.2$sugar  <- as.factor(Ex13.2$sugar)

 # perform RCB analysis with a cRF treatment structure

Ex13.2.lm1  <- lm(y ~ batch +temp*pH*sugar, data=Ex13.2)

anova(Ex13.2.lm1)

library('emmeans')

lsmeans(Ex13.2.lm1, ~pH)

lsmeans(Ex13.2.lm1, ~sugar)

 # check model assumptions

par(mfrow=c(2,2))  
plot(Ex13.2.lm1)
par(mfrow=c(1,1))


 # ---------------------------------------------------------------
 # Exercise 13.3
 # ---------------------------------------------------------------

#  first read the data directly from Dr. Oehlerts website

Ex13.3 <- read.table("http://www.stat.umn.edu/~gary/book/fcdae.data/ex13.3",header=TRUE)

Ex13.3

Ex13.3  <-  within(Ex13.3, 
                {
                   square                        <- NA
                   square[col < 5]               <- 1
                   square[col > 4]               <- 2
                })

Ex13.3$row  <- as.factor(Ex13.3$row)
Ex13.3$col  <- as.factor(Ex13.3$col)
Ex13.3$trt    <- as.factor(Ex13.3$trt)

 # perform Latin Square analysis 

Ex13.3.lm1  <- lm(y ~ row +col +trt, data=Ex13.3)

anova(Ex13.3.lm1)

 # alternatively

Ex13.3.lm2  <- lm(y ~ square +square:col +row +trt, data=Ex13.3)

anova(Ex13.3.lm2)

 # check model assumptions

par(mfrow=c(2,2))  
plot(Ex13.3.lm1)
par(mfrow=c(1,1))

 # Tukey tests

library(multcomp)

 # Tukey's multiple comparison method

trt.Tukey  <- glht(Ex13.3.lm1,linfct=mcp(trt="Tukey"))

summary(trt.Tukey)
plot(trt.Tukey,sub="Tukeys method")


 # ---------------------------------------------------------------
 # Exercise 14.4
 # ---------------------------------------------------------------

#  first read the data directly from Dr. Oehlerts website

Ex14.4 <- read.table("http://www.stat.umn.edu/~gary/book/fcdae.data/ex14.4",header=TRUE)

Ex14.4

Ex14.4$computer  <- as.factor(Ex14.4$computer)
Ex14.4$treatment    <- as.factor(Ex14.4$treatment)


 # perform BIBD Square analysis 

Ex14.4.lm1  <- lm(response ~ computer +treatment, data=Ex14.4)

anova(Ex14.4.lm1)

 # check model assumptions

par(mfrow=c(2,2))  
plot(Ex14.4.lm1)
par(mfrow=c(1,1))

 # Tukey tests

library(multcomp)

 # Tukey's multiple comparison method

trt.Tukey14.4  <- glht(Ex14.4.lm1,linfct=mcp(treatment="Tukey"))

summary(trt.Tukey14.4)
plot(trt.Tukey14.4,sub="Tukeys method")


 # ---------------------------------------------------------------
 # Problem 16.8
 # ---------------------------------------------------------------

#  first read the data directly from Dr. Oehlerts website

Pr16.8 <- read.table("http://www.stat.umn.edu/~gary/book/fcdae.data/pr16.8",header=TRUE)

Pr16.8

Pr16.8$sign  <- as.factor(Pr16.8$sign)
Pr16.8$timing  <- as.factor(Pr16.8$timing)
Pr16.8$interchange  <- as.factor(Pr16.8$interchange)


 # view residuals with model assuming fixed effects

Pr16.8.lm1  <- lm(speed ~ sign +sign:interchange +timing +sign:timing, data = Pr16.8)

anova(Pr16.8.lm1)

par(mfrow=c(2,2))  
plot(Pr16.8.lm1)
par(mfrow=c(1,1)) 


 # The afex library calls lmer and other functions to give results fairly 
 # similar to those in Proc MIXED

library(afex)

library('emmeans')   # needed for more recent versions of afex

Pr16.8.afex1  <- mixed(speed ~ sign +(1|sign:interchange) +timing +sign:timing,
                       data = Pr16.8)

summary(Pr16.8.afex1)

nice(Pr16.8.afex1)

lsmeans(Pr16.8.afex1, ~sign)

lsmeans(Pr16.8.afex1, ~timing)

# multiple comparisons with Tukeys method

emmeans(Pr16.8.afex1, list(pairwise ~ timing), adjust = "tukey")


 # examine potential interactions

lsmip(Pr16.8.afex1, sign ~ timing)