> # 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) The following object(s) are masked _by_ '.GlobalEnv': logcount, package > > package <- factor(package) > > bacteria package logcount 1 1 7.66 2 1 6.98 3 1 7.80 4 2 5.26 5 2 5.44 6 2 5.80 7 3 7.41 8 3 7.33 9 3 7.04 10 4 3.51 11 4 2.91 12 4 3.66 > > tapply(logcount,package,mean) 1 2 3 4 7.48 5.50 7.26 3.36 > > tapply(logcount,package,sd) 1 2 3 4 0.4386342 0.2749545 0.1946792 0.3968627 > > ex_02_01.mod <- aov(logcount ~ package) > > TukeyHSD(ex_02_01.mod,"package") Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = logcount ~ package) $package diff lwr upr p adj 2-1 -1.98 -2.869962 -1.090038 0.0004549 3-1 -0.22 -1.109962 0.669962 0.8563618 4-1 -4.12 -5.009962 -3.230038 0.0000020 3-2 1.76 0.870038 2.649962 0.0010160 4-2 -2.14 -3.029962 -1.250038 0.0002639 4-3 -3.90 -4.789962 -3.010038 0.0000031 > > 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) [1] "Overall Mean" [1] 5.9 [1] "Treatment Means" [1] 7.48 5.50 7.26 3.36 [1] "Treatment Variances" [1] 0.1924 0.0756 0.0379 0.1575 [1] "Treatment Replicates" [1] 3 3 3 3 [1] "Number of Treatments" [1] 4 [1] "Error Sum of Squares" [1] 0.9268 [1] "Mean Square Error" [1] 0.11585 [1] "Bonferroni Simultaneous CI's" "Alpha=" [3] "0.05" [1] "Bonferroni MSD" "Alpha=" "0.05" [1] 1.00 2.00 0.97 [1] 1.00 2.00 1.01 2.95 [1] "Bonferroni MSD" "Alpha=" "0.05" [1] 1.00 3.00 0.97 [1] 1.00 3.00 -0.75 1.19 [1] "Bonferroni MSD" "Alpha=" "0.05" [1] 1.00 4.00 0.97 [1] 1.00 4.00 3.15 5.09 [1] "Bonferroni MSD" "Alpha=" "0.05" [1] 2.00 3.00 0.97 [1] 2.00 3.00 -2.73 -0.79 [1] "Bonferroni MSD" "Alpha=" "0.05" [1] 2.00 4.00 0.97 [1] 2.00 4.00 1.17 3.11 [1] "Bonferroni MSD" "Alpha=" "0.05" [1] 3.00 4.00 0.97 [1] 3.00 4.00 2.93 4.87 >