C:\Rmisc>C:\R-2.9.0\bin\Rterm --vanilla 

R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 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)
[1] "Overall Mean"
[1] 55.23667
[1] "Treatment Means"
[1] 46.44000 57.67667 58.68111 58.14889
[1] "Number of Treatments"
[1] 4
[1] "Block Means"
[1] 41.8925 65.4800 67.9925 55.0075 57.3600 71.3975 43.0250 61.7850 33.1900
[1] "Number of Blocks"
[1] 9
[1] "Error Sum of Squares"
[1] 1261.657
[1] "Mean Square Error"
[1] 52.56904
[1] "Bonferroni MSD" "Alpha="         "0.05"          
[1] 9.826772
[1] "Bonferroni Simultaneous Confidence Intervals"
[2] "Alpha="                                      
[3] "0.05"                                        
[1]   1.00   2.00 -21.06  -1.41
[1]   1.00   3.00 -22.07  -2.41
[1]   1.00   4.00 -21.54  -1.88
[1]   2.00   3.00 -10.83   8.82
[1]   2.00   4.00 -10.30   9.35
[1]  3.00  4.00 -9.29 10.36
> 
> 
> 
> dev.off()
null device 
          1 
>