pdf("C:\\Rmisc\\graphs\\amoeba.pdf")

amoeba1 <- read.fwf("C:\\data\\amoeba.txt", width=c(8,8), col.names=c("Trt", "Yield"))

fTrt <- factor(amoeba1$Trt, levels=1:5)
levels(fTrt) <- c("None", "H", "F", "HF", "FH")

amoeba <- data.frame(amoeba1, fTrt)
attach(amoeba)

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(Yield,Trt,0.05)





dev.off()