# ---------------------------------------------------------- ;
# Exercise 6.3                                               ;
# ---------------------------------------------------------- ;


 input fabric temp shrink @@ ;
 cards ;

data <- scan(what = list(fabric=0,temp=0,shrink=0), multi.line=F)
 1 210 1.8  1 210 2.1   1 215 2.0  1 215 2.1 
 1 220 4.6  1 220 5.0   1 225 7.5  1 225 7.9 
 2 210 2.2  2 210 2.4   2 215 4.2  2 215 4.0 
 2 220 5.4  2 220 5.6   2 225 9.8  2 225 9.2 
 3 210 2.8  3 210 3.2   3 215 4.4  3 215 4.8 
 3 220 8.7  3 220 8.4   3 225 13.2  3 225 13.0 
 4 210 3.2  4 210 3.6   4 215 3.3  4 215 3.5 
 4 220 5.7  4 220 5.8   4 225 10.9  4 225 11.1 

data

hw2_2 <- data.frame(data)
rm(data)

hw2_2$fabric <- as.factor(hw2_2$fabric)
hw2_2$temp <- as.factor(hw2_2$temp)

hw2_2_lm  <- lm(shrink ~ fabric +temp +fabric:temp, data=hw2_2)
anova(hw2_2_lm)

# examine the interaction

interaction.plot(hw2_2$fabric,hw2_2$temp,
                 hw2_2$shrink,type="b")

# exploring polynomial contrasts of temp and their interaction with fabric

library(phia)

lineartemp.contrast <- list(temp = c(-3,-1,1,3))
quadtemp.contrast <- list(temp = c(1,-1,-1,1))
cubictemp.contrast <- list(temp = c(-1,3,-3,1))

testInteractions(hw2_2_lm,custom=lineartemp.contrast)
testInteractions(hw2_2_lm,custom=quadtemp.contrast)
testInteractions(hw2_2_lm,custom=cubictemp.contrast)

testInteractions(hw2_2_lm,custom=lineartemp.contrast,across="fabric")
testInteractions(hw2_2_lm,custom=quadtemp.contrast,across="fabric")
testInteractions(hw2_2_lm,custom=cubictemp.contrast,across="fabric")



# ---------------------------------------------------------- ;
# Exercise 6.11                                              ;
# ---------------------------------------------------------- ;

data <- scan(what = list(base=0,alcohol=0,prod=0), multi.line=F)
1 1 90.7  1 1 91.4  1 2 89.3  1 2 88.1  1 2 90.4  1 3 89.5  1 3 87.6  1 3 88.3  1 3 90.3 
2 1 87.3  2 1 88.3  2 1 91.5  2 2 94.7  2 3 93.1  2 3 90.7  2 3 91.5 

data
 
hw2_3 <- data.frame(data)
rm(data)


hw2_3$base  <- as.factor(hw2_3$base)
hw2_3$alcohol  <- as.factor(hw2_3$alcohol)

hw2_3_lm  <- lm(prod ~ base +alcohol +base:alcohol, data = hw2_3)
anova(hw2_3_lm)

library(phia)

# interactionMeans seems to give least squares means
interactionMeans(hw2_3_lm)
interactionMeans(hw2_3_lm,factors="base")
interactionMeans(hw2_3_lm,factors="alcohol")

# ordinary means
tapply(hw2_3$prod,hw2_3$base,mean)
tapply(hw2_3$prod,hw2_3$alcohol,mean)

# the testInteractions function can produce output like that produced 
# by the SAS ESTIMATE statements

testInteractions(hw2_3_lm,pairwise="base",fixed="alcohol")



# ---------------------------------------------------------- ;
# Table 6.13 data                                            ;
# ---------------------------------------------------------- ;

data <- scan(what =list(t=0, d=0, s=0, weightgain=0), multi.line=F)
  25 80 10 86  25 80 10 52  25 80 10 73
  25 80 25 544 25 80 25 371 25 80 25 482
  25 80 40 390 25 80 40 290 25 80 40 397
  25 160 10 53 25 160 10 73 25 160 10 86
  25 160 25 393 25 160 25 398 25 160 25 208
  25 160 40 249 25 160 40 265 25 160 40 243
  35 80 10 439  35 80 10 436  35 80 10 349
  35 80 25 249  35 80 25 245  35 80 25 330
  35 80 40 247  35 80 40 277  35 80 40 205
  35 160 10 324 35 160 10 305 35 160 10 364
  35 160 25 352 35 160 25 267 35 160 25 316
  35 160 40 188 35 160 40 223 35 160 40 281

data

hw2_4 <- data.frame(data)
rm(data)

hw2_4$t  <- as.factor(hw2_4$t)
hw2_4$d  <- as.factor(hw2_4$d)
hw2_4$s  <- as.factor(hw2_4$s)

hw2_4_lm  <- lm(weightgain ~ t*d*s, data = hw2_4)
anova(hw2_4_lm)


# use the afex library to get expected values of mean squares

library(afex)

ems(r ~ A*B*C, random="ABC")


# For a mixed model, the afex package can be used to obtain an F test
# for any fixed effect. However in this case we have a random effects model and 
# want to test for a variance component. The mixed function in 
# afex does not seem to like fitting models without any fixed effects.
# An alternative that can be used in this case is to fit two lmer models 
# from the lme4 package and then use the anova function to conduct a 
# likelihood-ratio test for the desired variance component:

library(lme4)


hw2_4_lmer1a <- lmer(weightgain ~ (1|t) +(1|d) +(1|s) +(1|t:d) +(1|t:s) +(1|d:s) +(1|t:d:s), 
      data = hw2_4)


hw2_4_lmer1b <- lmer(weightgain ~ (1|d) +(1|s) +(1|t:d) +(1|t:s) +(1|d:s) +(1|t:d:s), 
      data = hw2_4)


anova(hw2_4_lmer1a, hw2_4_lmer1b, test="Chisq")


# ---------------------------------------------------------- ;
# Exercise 7.3                                               ;
# ---------------------------------------------------------- ;

data <- scan(what =list(alloy="", casting=0, strength=0), multi.line=F)
 A 1 13.2  A 1 15.5  A 2 15.2  A 2 15.0  A 3 14.8  A 3 14.2  A 4 14.6  A 4 15.1 
 B 1 17.1  B 1 16.7  B 2 16.5  B 2 17.3  B 3 16.1  B 3 15.4  B 4 17.4  B 4 16.8 
 C 1 14.1  C 1 14.8  C 2 13.2  C 2 13.9  C 3 14.5  C 3 14.7  C 4 13.8  C 4 13.5 

data

hw2_5 <- data.frame(data)
rm(data)

hw2_5$alloy  <- as.factor(hw2_5$alloy)
hw2_5$casting  <- as.factor(hw2_5$casting)

hw2_5_lm  <- lm(strength ~ alloy + casting %in% alloy, data = hw2_5)
# anova here calculates the correct SS and MS terms, but not the right F ratios
anova(hw2_5_lm)

 # gives results much like Proc MIXED

library(afex)

hw2_5.afex1 <- mixed(strength ~ alloy + (1|casting:alloy), data = hw2_5)

summary(hw2_5.afex1)

nice(hw2_5.afex1)

lsmeans(hw2_5.afex1, ~alloy)



# ---------------------------------------------------------- ;
# Exercise 8.1 a,b,c,e,f                                     ;
# ---------------------------------------------------------- ;

data <- scan(what =list(method="", block=0, weight=0), multi.line=F)
  Tr 1 450  Tr 2 469  Tr 3 249  Tr 4 125  Tr 5 280  Tr 6 352  Tr 7 221  Tr 8 251 
  Ba 1 358  Ba 2 512  Ba 3 281  Ba 4  58  Ba 5 352  Ba 6 293  Ba 7 283  Ba 8 186 
  Spy 1 331 Spy 2 402 Spy 3 183 Spy 4 70  Spy 5 258 Spy 6 281 Spy 7 219 Spy 8 46 
  Spk 1 317 Spk 2 423 Spk 3 379 Spk 4 63  Spk 5 289 Spk 6 239 Spk 7 269 Spk 8 357 
  SSp 1 479 SSp 2 341 SSp 3 404 SSp 4 115 SSp 5 182 SSp 6 349 SSp 7 276 SSp 8 182 
  Fld 1 245 Fld 2 380 Fld 3 263 Fld 4 62  Fld 5 336 Fld 6 282 Fld 7 171 Fld 8 98 

data

hw2_6 <- data.frame(data)
rm(data)

hw2_6$method  <- as.factor(hw2_6$method)
hw2_6$method <- relevel(hw2_6$method,ref="Fld")     # setting the reference level for Dunnett

hw2_6$block  <- as.factor(hw2_6$block)

hw2_6_lm  <- lm(weight ~ block +method, data = hw2_6)
anova(hw2_6_lm)

 # Dunnett's multiple comparison method

library(multcomp)

hw2_6.Dunnett <- glht(hw2_6_lm,linfct=mcp(method="Dunnett"))

summary(hw2_6.Dunnett)
plot(hw2_6.Dunnett,sub="Problem 8.1")

 # check model assumptions

par(mfrow=c(2,2))  
plot(hw2_6_lm)
par(mfrow=c(1,1))  



# ---------------------------------------------------------- ;
# Exercise 8.5 a,b,d                                         ;
# ---------------------------------------------------------- ;

data <- scan(what =list(tech=0, timepd=0, construction="", time=0), multi.line=F)
  1 1 C 90  1 2 B 90  1 3 A 89  1 4 D 104 
  2 1 D 96  2 2 C 91  2 3 B 97  2 4 A 100 
  3 1 A 84  3 2 D 96  3 3 C 98  3 4 B 104 
  4 1 B 88  4 2 A 88  4 3 D 98  4 4 C 106 
 
data

hw2_7 <- data.frame(data)
rm(data)

hw2_7$tech  <- as.factor(hw2_7$tech)
hw2_7$timepd  <- as.factor(hw2_7$timepd)
hw2_7$construction  <- as.factor(hw2_7$construction)

hw2_7_lm  <- lm(time ~ tech +timepd +construction, data = hw2_7)
anova(hw2_7_lm)

 # Tukey's multiple comparison method

library(multcomp)

hw2_7.Tukey <- glht(hw2_7_lm,linfct=mcp(construction="Tukey"))

summary(hw2_7.Tukey)
plot(hw2_7.Tukey,sub="Problem 8.5")


 # check model assumptions

par(mfrow=c(2,2))  
plot(hw2_7_lm)
par(mfrow=c(1,1))