#Example -  Randomised Block Design 
# Testing hypothesis in two-way ANOVA

#Set working directory
setwd("~/Dropbox/ISI SCB/Data sets")

# read in csv file
tl <-read.csv("tools.csv", TRUE)

Time <- tl$CuttingTime
Tool <- tl$Tool
Material<- tl$Material

#Display the data in boxplots
boxplot(Time~Tool)
boxplot(Time~Material)


#Two-way ANOVA
time1.aov<- aov(Time~Tool+Material)
summary(time1.aov)

#Multiple Comparisons
time1.tukey <- TukeyHSD(time1.aov)
time1.tukey
plot(time1.tukey)

#residuals and fitted values of model
res <- time1.aov$residuals
fit <- time1.aov$fitted.values
#residual plot to assess constant variance assumption
plot(fit,res, main="Residual Plot",xlab= "Fitted Values", ylab="Residuals")


#assess normality assumption
qqnorm(res, ylab = "Residuals Sample Quantiles")
shapiro.test(res)



#Test for constant variance (Levene's Test)
#first transform the data and then perform
#one-way ANOVA on transformed data 

#determine the mean for each treatment
yt1 <- tl$CuttingTime[tl$Tool=="Tool1"]
myt1 <- mean(yt1)

yt2 <- tl$CuttingTime[tl$Tool=="Tool2"]
myt2 <- mean(yt2)

yt3 <- tl$CuttingTime[tl$Tool=="Tool3"]
myt3 <- mean(yt3)

yt4 <- tl$CuttingTime[tl$Tool=="Tool4"]
myt4 <- mean(yt4)

#determine the mean for each block
yb1 <- tl$CuttingTime[tl$Material=="Material1"]
myb1 <- mean(yb1)

yb2 <- tl$CuttingTime[tl$Material=="Material2"]
myb2 <- mean(yb2)

yb3 <- tl$CuttingTime[tl$Material=="Material3"]
myb3 <- mean(yb3)

yb4 <- tl$CuttingTime[tl$Material=="Material4"]
myb4 <- mean(yb4)

yb5 <- tl$CuttingTime[tl$Material=="Material5"]
myb5 <- mean(yb5)

#create vector of treatment means
ytmean <-as.matrix(rbind(myt1,myt2,myt3,myt4))
ymean <-mean(tl$CuttingTime)

#create a data set with the 4 treatments (tool type)
#for each of the 5 blocks (material type)
tt1 <- matrix(nrow=4,ncol=2)
tt2 <- matrix(nrow=4,ncol=2)
tt3 <- matrix(nrow=4,ncol=2)
tt4 <- matrix(nrow=4,ncol=2)
tt5 <- matrix(nrow=4,ncol=2)

#assign the the tool types and material types in the second column
tt1[,2] <- tl[1:4,2]
tt2[,2] <- tl[1:4,2]
tt3[,2] <- tl[1:4,2]
tt4[,2] <- tl[1:4,2]
tt5[,2] <- tl[1:4,2]

# computer the transformed cutting time 
#(cutting time (i,j) - Tool(i) mean - Material (j) mean + Grand mM
tt1[,1] = abs(tl[1:4,1]-ytmean-myb1+ymean)
tt2[,1] = abs(tl[5:8,1]-ytmean-myb2+ymean)
tt3[,1] = abs(tl[9:12,1]-ytmean-myb3+ymean)
tt4[,1] = abs(tl[13:16,1]-ytmean-myb4+ymean)
tt5[,1] = abs(tl[17:20,1]-ytmean-myb5+ymean)

rtime<-rbind(tt1,tt2,tt3,tt4,tt5)

#assess constant variance assumption

#One-way ANOVA on transformed data set (rtime) to get outcome of
#Levene's test for constant variances
rtime.aov<- aov(rtime[,1]~rtime[,2])
summary(rtime.aov)



#One-way ANOVA for only treatments (Tool type) and ignore
#blocks (material type)
time2.aov<- aov(Time~Tool)
summary(time2.aov)