#pdf("C:\\Rmisc\\graphs\\amoeba.pdf") amoeba1 <- read.table("http://www.stat.ufl.edu/~winner/data/entozamoeba.dat", header=F, 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) i1.i2 <- matrix(rep(0,2*num.pairs),ncol=2) ybar.diff <- rep(0,num.pairs) bonf.msd <- rep(0,num.pairs) bonf.LL <- rep(0,num.pairs) bonf.UL <- rep(0,num.pairs) pair <- 0 for(i1 in 1:(num.trt-1)) { for (i2 in (i1+1):num.trt) { pair <- pair + 1 i1.i2[pair,] <- cbind(i1,i2) ybar.diff[pair] <- ybar.trt[i1] - ybar.trt[i2] bonf.msd[pair] <- qt(1-alpha/(2*num.pairs),(sum(rep.trt)-num.trt))*sqrt(mse*(1/rep.trt[i1]+1/rep.trt[i2])) bonf.LL[pair] <- (ybar.trt[i1]-ybar.trt[i2])-bonf.msd[pair] bonf.UL[pair] <- (ybar.trt[i1]-ybar.trt[i2])+bonf.msd[pair] }} bonf.out <- cbind(i1.i2,ybar.diff,bonf.msd,bonf.LL,bonf.UL) colnames(bonf.out) <- c("Trt i", "Trt j", "Ybar diff", "MSD", "LL", "UL") bonf.out } bonf(Yield,fTrt,0.05) #dev.off()