# ------------ problem 17.2 -----------------------------


data <- scan(what = list(diet=0,y=0,x=0), multi.line=F)
 1 48 35  2 65 40  3 79 51  4 59 53  
 1 67 44  2 49 45  3 52 41  4 50 52 
 1 78 44  2 37 37  3 63 47  4 59 52 
 1 69 51  2 73 53  3 65 47  4 42 51 
 1 53 47  2 63 42  3 67 48  4 34 43 

data

prob17.2 <- data.frame(data)
rm(data)

prob17.2$diet <- as.factor(prob17.2$diet)

 # look at the data

plot(prob17.2$x,prob17.2$y,pch=as.numeric(prob17.2$diet),
      col=as.numeric(prob17.2$diet),ylab="y",xlab="x")
legend("bottomright",levels(prob17.2$diet),pch=1:4,col=1:4)

 # standard ANCOVA assuming equal slopes

prob17.2.lm1 <- lm(y ~ x +diet , data=prob17.2)
anova(prob17.2.lm1)
summary(prob17.2.lm1)

 # add lines to plot

abline(-9.29, 1.64,col="black")
abline(-9.29 -4.29, 1.64,col="red")
abline(-9.29 -2.05, 1.64,col="green")
abline(-9.29 -24.01, 1.64,col="blue")


 # calculate adjusted treatment means and se's

library(effects)

effect("diet",prob17.2.lm1) 
effect("diet",prob17.2.lm1)$se


 # check assumptions

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


 # to compute efficiency and average s.e. of differences

Txx <- anova(lm(x ~ diet , data=prob17.2))$"Sum Sq"[1]

Exx <- anova(lm(x ~ diet , data=prob17.2))$"Sum Sq"[2]

MSE <- anova(lm(y ~ x +diet , data=prob17.2))$"Mean Sq"[3]

MSEr <- anova(lm(y ~ diet , data=prob17.2))$"Mean Sq"[2]

t <- 4
r <- 5

c(Txx,Exx,MSE,MSEr)

Efficiency <- MSEr/(MSE*(1 + Txx/((t-1)*Exx) ) )

AvgSEDiff <- sqrt(2*(MSE/r)*(1 + Txx/((t-1)*Exx) ) )

c(Efficiency,AvgSEDiff)



 # check equal slopes assumption

prob17.2.lm2 <- lm(y ~ x +diet +x:diet, data=prob17.2)
anova(prob17.2.lm2)





 # ------------ problem 17.1 -----------------------------


data <- scan(what = list(alloy=0,y=0,x=0), multi.line=F)
1 37.5 12.5  1 40.5 14.0  1 49.0 16.0  
1 51.0 15.0  1 61.5 18.0  1 63.0 19.5 
2 57.5 16.5  2 69.5 17.5  2 87.0 19.0 
2 92.0 19.5  2 107.0 24.0  2 119.5 22.5 
3 38.0 15.5  3 44.5 16.0  3 53.0 19.0 
3 55.0 18.0  3 58.5 19.0  3 60.0 20.5 

data

prob17.1 <- data.frame(data)
rm(data)

prob17.1$alloy <- as.factor(prob17.1$alloy)


 # look at the data

plot(prob17.1$x,prob17.1$y,pch=as.numeric(prob17.1$alloy),
      col=as.numeric(prob17.1$alloy),ylab="y",xlab="x")
legend("bottomright",levels(prob17.1$alloy),pch=1:3,col=1:3)

 # standard ANCOVA assuming equal slopes

prob17.1.lm1 <- lm(y ~ x +alloy , data=prob17.1)
anova(prob17.1.lm1)
summary(prob17.1.lm1)


 # add lines to plot

abline(-36.751,5.505,col="black")
abline(-36.751 +16.312,5.505,col="red")
abline(-36.751 -10.845,5.505,col="green")

 # calculate adjusted treatment means and se's

library(effects)

effect("alloy",prob17.1.lm1) 
effect("alloy",prob17.1.lm1)$se

 # check assumptions

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

 # to compute efficiency and average s.e. of differences

Txx <- anova(lm(x ~ alloy , data=prob17.1))$"Sum Sq"[1]

Exx <- anova(lm(x ~ alloy , data=prob17.1))$"Sum Sq"[2]

MSE <- anova(lm(y ~ x +alloy , data=prob17.1))$"Mean Sq"[3]

MSEr <- anova(lm(y ~ alloy , data=prob17.1))$"Mean Sq"[2]

t <- 3
r <- 6

c(Txx,Exx,MSE,MSEr)

Efficiency <- MSEr/(MSE*(1 + Txx/((t-1)*Exx) ) )

AvgSEDiff <- sqrt(2*(MSE/r)*(1 + Txx/((t-1)*Exx) ) )


 # check equal slopes assumption

prob17.1.lm2 <- lm(y ~ x +alloy +x:alloy, data=prob17.1)
anova(prob17.1.lm2)

 # differences at three weld diameter quartiles

summary(prob17.1$x)


library(multcomp)

 # at x = 16.

contrdiff16 <- rbind("one - two at 16."  <- c(0,0,-1,0,-16.,0),
                    "two - three at 16."     <- c(0,0,1,-1,16.,-16.),
                    "one - three at 16."      <- c(0,0,0,-1,0,-16.))
   
summary(glht(prob17.1.lm2, linfct= contrdiff16))

 # at x = 18.

contrdiff18 <- rbind("one - two at 18."  <- c(0,0,-1,0,-18.,0),
                    "two - three at 18."     <- c(0,0,1,-1,18.,-18.),
                    "one - three at 18."      <- c(0,0,0,-1,0,-18.))
   
summary(glht(prob17.1.lm2, linfct= contrdiff18))


# at x = 19.4:

contrdiff19.4 <- rbind("one - two at 19.4"  <- c(0,0,-1,0,-19.4,0),
                    "two - three at 19.4"     <- c(0,0,1,-1,19.4,-19.4),
                    "one - three at 19.4"      <- c(0,0,0,-1,0,-19.4))
   
summary(glht(prob17.1.lm2, linfct= contrdiff19.4))