# Two factor experiment and analysis

rm(list = ls())  # clear the memory

torque = read.table("http://www.stat.ualberta.ca/~wiens/stat568/datasets/bolt_table_3_8.dat", header = T, quote= "")

torque

y = c(as.matrix(torque[,1:3]))
test = as.factor(rep(c(rep("Mandrel",10), rep("Bolt", 10)),3))
plate = as.factor(c(rep("PO",20), rep("CW", 20), rep("HT",20)))
data = data.frame(y,test,plate)

par(mfrow=c(2,2))
boxplot(y ~ test, xlab = "Factor A: test")
boxplot(y ~ plate, xlab = "Factor B: plate") 

interaction.plot(test, plate, y)
interaction.plot(plate, test, y)
# interaction effect shown

fit = aov(y ~ test + plate + test*plate, data)
summary(fit)

###########################
# Check industry standard:

s2 = vector(length=6)
for (j in 1:6) s2[j] = var(y[(1+10*(j-1)):(10*j)]) # variance in each cell
s2
test2 =  test[1+10*c(0:5)]
plate2 = plate[1+10*c(0:5)]
fits2 = aov(log(s2)~ test2 + plate2 + test2*plate2)
summary(fits2)
sigma2 = exp(fits2$fitted)
meanfits = fit$fitted[1+10*c(0:5)]
limits = meanfits + 2*sqrt(sigma2)

mat = cbind(test2, plate2, meanfits, sqrt(sigma2), limits)
mat


# < 45 for all treatments except test = 1 and plate = 2,3
print(cbind(as.character(data[,2]),as.numeric(data[,2])),quote=F) # test = 1 is Bolt
print(cbind(as.character(data[,3]),as.numeric(data[,3])),quote=F) # plate = 2,3 is HT, PO
###########################

# look at the residuals:
windows()
par(mfrow=c(2,1))
resid <- fit$residuals
qqnorm(resid)
qqline(resid)
text(0, .9*max(resid), paste("p = ", round(shapiro.test(resid)[[2]],3)))

fits <- fit$fitted.values
plot(fits,resid)
abline(h=0)

###########################
# test equality of variances in various ways:
bartlett.test(resid,test)
bartlett.test(resid,plate)

library(car)
leveneTest(y, group = test)
leveneTest(y, group = plate)
###########################
# A rank test:
rankfit = aov(rank(y) ~ test + plate + test*plate)
summary(rankfit)

# Try a transformation:

logfit = aov(log(y) ~ test + plate + test*plate)
summary(logfit)

windows()
par(mfrow=c(2,1))
resid.log <- logfit$residuals
qqnorm(resid.log)
qqline(resid.log)
text(0, .9*max(resid.log), paste("p = ", round(shapiro.test(resid.log)[[2]],3)))

fits.log <- logfit$fitted.values
plot(fits.log,resid.log)
abline(h=0)

bartlett.test(resid.log,test)
bartlett.test(resid.log,plate)

windows()
TukeyHSD(logfit, ordered = TRUE)
par(cex=.5)
plot(TukeyHSD(logfit, ordered = TRUE, which = 'test:plate'))