> kenton.aov <- aov(cases ~ package)
> summary(kenton.aov)
            Df Sum Sq Mean Sq F value    Pr(>F)    
package      3 588.22 196.074  18.591 2.585e-05 ***
Residuals   15 158.20  10.547                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> TukeyHSD(kenton.aov,"package")
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = cases ~ package)

$package
    diff        lwr       upr     p adj
2-1 -1.2 -7.1197584  4.719758 0.9352978
3-1  4.9 -1.3788520 11.178852 0.1548895
4-1 12.6  6.6802416 18.519758 0.0001013
3-2  6.1 -0.1788520 12.378852 0.0582866
4-2 13.8  7.8802416 19.719758 0.0000368
4-3  7.7  1.4211480 13.978852 0.0142180

> 
> 
> 
> bonf_schef <- 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)))
+ }}
+ 
+ print(c("Scheffe' Simultaneous CI's","Alpha=",alpha))
+ for(i1 in 1:(num.trt-1)) {
+ for (i2 in (i1+1):num.trt) {
+ schef.msd <- sqrt((num.trt-1)*qf(1-alpha,num.trt-1,sum(rep.trt)-num.trt))*sqrt(mse*(1/rep.trt[i1]+1/rep.trt[i2]))
+ print(c("Scheffe' MSD","Alpha=",alpha))
+ print(c(round(c(i1,i2,schef.msd),2)))
+ print(c(round(c(i1,i2),0),
+ round(c((ybar.trt[i1]-ybar.trt[i2])-schef.msd,(ybar.trt[i1]-ybar.trt[i2])+schef.msd),2)))
+ }}
+ 
+ 
+ }
> 
> bonf_schef(cases,package,0.05)
[1] "Overall Mean"
[1] 18.63158
[1] "Treatment Means"
[1] 14.6 13.4 19.5 27.2
[1] "Treatment Variances"
[1]  5.3 13.3  7.0 15.7
[1] "Treatment Replicates"
[1] 5 5 4 5
[1] "Number of Treatments"
[1] 4
[1] "Error Sum of Squares"
[1] 158.2
[1] "Mean Square Error"
[1] 10.54667
[1] "Bonferroni Simultaneous CI's" "Alpha="                      
[3] "0.05"                        
[1] "Bonferroni MSD" "Alpha="         "0.05"          
[1] 1.00 2.00 6.24
[1]  1.00  2.00 -5.04  7.44
[1] "Bonferroni MSD" "Alpha="         "0.05"          
[1] 1.00 3.00 6.61
[1]   1.00   3.00 -11.51   1.71
[1] "Bonferroni MSD" "Alpha="         "0.05"          
[1] 1.00 4.00 6.24
[1]   1.00   4.00 -18.84  -6.36
[1] "Bonferroni MSD" "Alpha="         "0.05"          
[1] 2.00 3.00 6.61
[1]   2.00   3.00 -12.71   0.51
[1] "Bonferroni MSD" "Alpha="         "0.05"          
[1] 2.00 4.00 6.24
[1]   2.00   4.00 -20.04  -7.56
[1] "Bonferroni MSD" "Alpha="         "0.05"          
[1] 3.00 4.00 6.61
[1]   3.00   4.00 -14.31  -1.09
[1] "Scheffe' Simultaneous CI's" "Alpha="                    
[3] "0.05"                      
[1] "Scheffe' MSD" "Alpha="       "0.05"        
[1] 1.00 2.00 6.45
[1]  1.00  2.00 -5.25  7.65
[1] "Scheffe' MSD" "Alpha="       "0.05"        
[1] 1.00 3.00 6.84
[1]   1.00   3.00 -11.74   1.94
[1] "Scheffe' MSD" "Alpha="       "0.05"        
[1] 1.00 4.00 6.45
[1]   1.00   4.00 -19.05  -6.15
[1] "Scheffe' MSD" "Alpha="       "0.05"        
[1] 2.00 3.00 6.84
[1]   2.00   3.00 -12.94   0.74
[1] "Scheffe' MSD" "Alpha="       "0.05"        
[1] 2.00 4.00 6.45
[1]   2.00   4.00 -20.25  -7.35
[1] "Scheffe' MSD" "Alpha="       "0.05"        
[1] 3.00 4.00 6.84
[1]   3.00   4.00 -14.54  -0.86
> 
> 
>