pdf("juice_ls.pdf")

juice_ls <- read.fwf("C:\\data\\juice_ls.txt", width=c(8,8,8,8),
col.names=c("month","container","store","sales"))

attach(juice_ls)

month <- factor(month)
container <- factor(container)
store <- factor(store)

model <- aov(sales ~ container + month + store)

summary(model)

TukeyHSD(model,"container")

bonf <- function(y,trt,rowblk,colblk,alpha) {

ybar <- mean(y)
print("Overall Mean")
print(ybar)


ybar.trt <- as.vector(tapply(y,trt,mean))
print("Treatment Means")
print(ybar.trt)
num.trt <- length(ybar.trt)
print("Number of Treatments")
print(num.trt)

ybar.rblk <- as.vector(tapply(y,rowblk,mean))
print("Row Block Means")
print(ybar.rblk)
num.rblk <- length(ybar.rblk)
print("Number of Row Blocks")
print(num.rblk)

ybar.cblk <- as.vector(tapply(y,colblk,mean))
print("Column Block Means")
print(ybar.cblk)
num.cblk <- length(ybar.cblk)
print("Number of Column Blocks")
print(num.cblk)

ss.total <- sum((y-ybar)^2)
ss.trt <- num.rblk*sum((ybar.trt-ybar)^2)
ss.cblk <- num.rblk*sum((ybar.cblk-ybar)^2)
ss.rblk <- num.cblk*sum((ybar.rblk-ybar)^2)
sse <- ss.total-ss.trt-ss.cblk-ss.rblk  

print("Error Sum of Squares")
print(sse)

mse=sse/((num.trt-1)*(num.trt-2))
print("Mean Square Error")
print(mse)

num.pairs <- num.trt*(num.trt-1)/2

bonf.msd <- qt(1-alpha/(2*num.pairs),(num.trt-1)*(num.rblk-2))*sqrt(mse*(2/num.rblk))
print(c("Bonferroni MSD","Alpha=",alpha))
print(bonf.msd)

print(c("Bonferroni Simultaneous CIs","Alpha=",alpha))
for(i1 in 1:(num.trt-1)) {
for (i2 in (i1+1):num.trt) {
print(c(round(c(i1,i2,(ybar.trt[i1]-ybar.trt[i2])-bonf.msd,(ybar.trt[i1]-ybar.trt[i2])+bonf.msd),2)))
}}


}

bonf(sales,container,month,store,0.05)



dev.off()