# pdf("ex_02_01.pdf")

package <- c(rep(1:4,each=3))
logcount <- c(7.66,6.98,7.80,5.26,5.44,5.80,7.41,7.33,7.04,3.51,2.91,3.66)

bacteria <- data.frame(package,logcount)

attach(bacteria)

package <- factor(package)

bacteria

tapply(logcount,package,mean)

tapply(logcount,package,sd)

ex_02_01.mod <- aov(logcount ~ package)

TukeyHSD(ex_02_01.mod,"package")

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

ybar <- mean(y)
print("Overall Mean")
print(ybar)
ybar.trt <- as.vector(tapply(y,trt,mean))
print("Treatment Means")
print(ybar.trt)
var.trt <- as.vector(tapply(y,trt,var))
print("Treatment Variances")
print(var.trt)
rep.trt <- as.vector(tapply(y,trt,length))
print("Treatment Replicates")
print(rep.trt)
num.trt <- length(ybar.trt)
print("Number of Treatments")
print(num.trt)

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


for (i1 in 1:num.trt) {
sse <- sse + (rep.trt[i1]-1)*var.trt[i1]
}
print("Error Sum of Squares")
print(sse)

mse=sse/(sum(rep.trt)-num.trt)
print("Mean Square Error")
print(mse)



print(c("Bonferroni Simultaneous CI's","Alpha=",alpha))
for(i1 in 1:(num.trt-1)) {
for (i2 in (i1+1):num.trt) {
bonf.msd <- qt(1-alpha/(2*num.pairs),(sum(rep.trt)-num.trt))*sqrt(mse*(1/rep.trt[i1]+1/rep.trt[i2]))
print(c("Bonferroni MSD","Alpha=",alpha))
print(c(round(c(i1,i2,bonf.msd),2)))
print(c(round(c(i1,i2),0),round(c((ybar.trt[i1]-ybar.trt[i2])-bonf.msd,(ybar.trt[i1]-ybar.trt[i2])+bonf.msd),2)))
}}


}

bonf(logcount,package,0.05)



# dev.off()