> 
> trt <- factor(trt)
> 
> ranks_y <- rank(y)
> 
> (mean_y <- as.vector(tapply(y,trt,mean)))
[1]  50.370  22.130 121.208
> 
> (std_y <- as.vector(tapply(y,trt,sd)))
[1]  42.29353  33.21828 127.15128
> 
> (n <- as.vector(tapply(y,trt,length)))
[1] 5 5 5
> r <- length(n)
> nT <- sum(n)
> 
> 
> print(cbind(trt,y,ranks_y))
      trt      y ranks_y
 [1,]   1   4.41       2
 [2,]   1 100.65      13
 [3,]   1  14.45       6
 [4,]   1  47.13       9
 [5,]   1  85.21      12
 [6,]   2   8.24       4
 [7,]   2  81.16      11
 [8,]   2   7.35       3
 [9,]   2  12.29       5
[10,]   2   1.61       1
[11,]   3 106.19      14
[12,]   3  33.83       7
[13,]   3  78.88      10
[14,]   3 342.81      15
[15,]   3  44.33       8
> 
> (std_y^2/mean_y)
[1]  35.51206  49.86237 133.38597
> 
> (std_y/mean_y)
[1] 0.839657 1.501052 1.049034
> 
> (std_y/mean_y^2)
[1] 0.016669785 0.067828835 0.008654822
> 
> 
> library(MASS)
> 
> boxcox(y~trt)
> 
> rank.aov <- aov(ranks_y ~ trt)
> summary(rank.aov)
            Df Sum Sq Mean Sq F value  Pr(>F)  
trt          2   91.2  45.600  2.8983 0.09399 .
Residuals   12  188.8  15.733                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> 
> kruskal.test(y~trt)

        Kruskal-Wallis rank sum test

data:  y by trt 
Kruskal-Wallis chi-squared = 4.56, df = 2, p-value = 0.1023

> 
> z_bon <- qnorm(1-(.10/(2*r*(r-1)/2)))
> for (i1 in 1:(r-1)) {
+ for (i2 in (i1+1):r) {
+ diff <- (sum(ranks_y[trt==i1])/n[i1]) - (sum(ranks_y[trt==i2])/n[i2])
+ se_diff <- sqrt((nT*(nT+1)/12)*((1/n[i1])+(1/n[i2])))
+ print(cbind(i1,i2,diff-z_bon*se_diff,diff+z_bon*se_diff))
+ }}
     i1 i2                  
[1,]  1  2 -2.419021 9.61902
     i1 i2                   
[1,]  1  3 -8.419021 3.619021
     i1 i2                     
[1,]  2  3 -12.01902 0.01902086
>