csrbd <- read.table("http://www.stat.ufl.edu/~winner/data/chopstick2_rcb.dat", header=F, col.names=c("pinch.Y", "trt.cs", "block.cs")) attach(csrbd) trt.cs <- factor(trt.cs) block.cs <- factor(block.cs) cs.mod1 <- aov(pinch.Y ~ trt.cs + block.cs) summary(cs.mod1) TukeyHSD(cs.mod1, "trt.cs") 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=F,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) i1.i2 <- matrix(rep(0,2*num.pairs),ncol=2) ybar.diff <- rep(0,num.pairs) bonf.msd1 <- 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.msd1[pair] <- bonf.msd bonf.LL[pair] <- (ybar.trt[i1]-ybar.trt[i2])-bonf.msd1[pair] bonf.UL[pair] <- (ybar.trt[i1]-ybar.trt[i2])+bonf.msd1[pair] }} bonf.out <- cbind(i1.i2,ybar.diff,bonf.msd1,bonf.LL,bonf.UL) colnames(bonf.out) <- c("Trt i", "Trt j", "Ybar diff", "MSD", "LL", "UL") bonf.out } bonf(pinch.Y,trt.cs,block.cs,0.05) friedman.test(pinch.Y ~ trt.cs | block.cs)