> #nonparametric alternative, Kruskal Wallis (not in SW) Verzani 11.1.4

clergy.dat is unstacked  # outcome, knowledge mental health issues

> clergy  = read.table(file="D:\\stat141\\clergy.dat", header = T)
> print(clergy)
   Methodist Catholic Pentecostal
1         32       32          28
2         30       32          21
3         30       26          15
4         29       26          15
5         26       22          14
6         23       20          14
7         20       19          14
8         19       16          11
9         18       14           9
10        12       14           8

> summary(clergy)
   Methodist        Catholic      Pentecostal   
 Min.   :12.00   Min.   :14.00   Min.   : 8.00  
 1st Qu.:19.25   1st Qu.:16.75   1st Qu.:11.75  
 Median :24.50   Median :21.00   Median :14.00  
 Mean   :23.90   Mean   :22.10   Mean   :14.90  
 3rd Qu.:29.75   3rd Qu.:26.00   3rd Qu.:15.00  
 Max.   :32.00   Max.   :32.00   Max.   :28.00  
> attach(clergy)

> cstack = stack(list(M = Methodist, C = Catholic, P = Pentecostal))
> names(cstack)
[1] "values" "ind"   
> attach(cstack)
> detach(clergy)
> kruskal.test(values~ind)

        Kruskal-Wallis rank sum test

data:  values by ind 
Kruskal-Wallis chi-squared = 8.4706, df = 2, p-value = 0.01448


Look at rank transform of mental health data
> crank = rank(values)
> tapply(crank, ind, list)  # breakdown by religious group
$C
 [1] 29.0 29.0 22.0 22.0 19.0 16.5 14.5 12.0  7.0  7.0
$M
 [1] 29.0 26.5 26.5 25.0 22.0 20.0 16.5 14.5 13.0  4.0
$P
 [1] 24.0 18.0 10.5 10.5  7.0  7.0  7.0  3.0  2.0  1.0

> tapply(crank, ind, median)
    C     M     P 
17.75 21.00  7.00 
> tapply(crank, ind, quantile)
$C                                      $M                                     
   0%    25%    50%    75%   100%         0%    25%    50%    75%   100%  
7.000 12.625 17.750 22.000 29.000      4.000 15.000 21.000 26.125 29.000  
 $P                  
 0% 25% 50%  75% 100%
1.0 4.0 7.0 10.5 24.0