pdf("caffrbd.pdf")

caffrbd <- read.fwf("C:\\data\\caffeine.txt", width=c(1,5,8), 
col.names=c("subject","dose","enduretime"))
attach(caffrbd)


subject <- factor(subject, levels=1:9, 
 labels=c("S1","S2","S3","S4","S5","S6","S7","S8","S9"))
dose <- factor(dose, levels=c(0,5,9,13), 
labels=c("0mg", "5mg", "9mg", "13mg"))
 

bonf <- function(y,trt,blk,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.blk <- as.vector(tapply(y,blk,mean))
print("Block Means")
print(ybar.blk)
num.blk <- length(ybar.blk)
print("Number of Blocks")
print(num.blk)
sse <- 0
num.pairs <- num.trt*(num.trt-1)/2
y.mat <- matrix(y,byrow=T,ncol=num.trt)
#print(y.mat)

for (i1 in 1:num.blk) {
for (i2 in 1:num.trt) {
sse <- sse + (y.mat[i1,i2]-ybar.trt[i2]-ybar.blk[i1]+ybar)^2
}}
print("Error Sum of Squares")
print(sse)

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

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

print(c("Bonferroni Simultaneous Confidence Intervals","Alpha=",alpha))
for(i1 in 1:(num.trt-1)) {
for (i2 in (i1+1):num.trt) {
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(enduretime,dose,subject,0.05)



dev.off()