## Nonparametric tests

############################
# 1. Kolmogorov-Smirnov test
############################

y1<-c(.2,.3,.4,1.1,2,2.1,3.3,3.8,4.5,4.8,4.9,5,5.3,7.5,9.8,
  10.4,10.9,11.3,12.4,16.2,17.6,18.9,20.7,24,25.4,40,42.2,50,60)
y2<-c(.2,.3,.4,.7,1.2,1.5,1.5,1.9,2,2.4,2.5,2.8,3.6,4.8,4.8,
  5.4,5.7,5.8,7.5,8.7,8.8,9.1,10.3,15.6,16.1,16.5,16.7,20,20.7,33)
#table 3.10

library(survival)

delta1<-rep(1,length(y1))
fit1<-survfit(Surv(y1,delta1)~1)
summary(fit1)
help("survfit.object")
x1<-fit1$time
f1<-1-fit1$surv

delta2<-rep(1,length(y2))
fit2<-survfit(Surv(y2,delta2)~1)
summary(fit2)
x2<-fit2$time
f2<-1-fit2$surv

## make table 3.13
x<-unique(sort(c(x1,x2)))
> length(x)
[1] 50
ff1<-rep(0,50)
ff2<-rep(0,50)
for (i in 1:50) {
  ff1[i]<-ifelse(is.na(match(x[i],x1)),
      f1[ sum(x1<= max(x1[x1<=x[i]]) ) ],  f1[match(x[i],x1)])
  ff2[i]<-ifelse(is.na(match(x[i],x2)),
      f2[ sum(x2<= max(x2[x2<=x[i]]) ) ],  f2[match(x[i],x2)])
}

cbind(x,round(ff1,3) ,round(ff2,3))
 [1,]  0.2 0.034 0.033
 [2,]  0.3 0.069 0.067
 [3,]  0.4 0.103 0.100
 [4,]  0.7 0.103 0.133
 [5,]  1.1 0.138 0.133
 [6,]  1.2 0.138 0.167
 [7,]  1.5 0.138 0.233
 [8,]  1.9 0.138 0.267
 [9,]  2.0 0.172 0.300
[10,]  2.1 0.207 0.300
[11,]  2.4 0.207 0.333
[12,]  2.5 0.207 0.367
[13,]  2.8 0.207 0.400
[14,]  3.3 0.241 0.400
[15,]  3.6 0.241 0.433
[16,]  3.8 0.276 0.433
[17,]  4.5 0.310 0.433
[18,]  4.8 0.345 0.500
[19,]  4.9 0.379 0.500
[20,]  5.0 0.414 0.500
[21,]  5.3 0.448 0.500
[22,]  5.4 0.448 0.533
[23,]  5.7 0.448 0.567
[24,]  5.8 0.448 0.600
[25,]  7.5 0.483 0.633
[26,]  8.7 0.483 0.667
[27,]  8.8 0.483 0.700
[28,]  9.1 0.483 0.733
[29,]  9.8 0.517 0.733
[30,] 10.3 0.517 0.767
[31,] 10.4 0.552 0.767
[32,] 10.9 0.586 0.767
[33,] 11.3 0.621 0.767
[34,] 12.4 0.655 0.767
[35,] 15.6 0.655 0.800
[36,] 16.1 0.655 0.833
[37,] 16.2 0.690 0.833
[38,] 16.5 0.690 0.867
[39,] 16.7 0.690 0.900
[40,] 17.6 0.724 0.900
[41,] 18.9 0.759 0.900
[42,] 20.0 0.759 0.933
[43,] 20.7 0.793 0.967
[44,] 24.0 0.828 0.967
[45,] 25.4 0.862 0.967
[46,] 33.0 0.862 1.000
[47,] 40.0 0.897 1.000
[48,] 42.2 0.931 1.000
[49,] 50.0 0.966 1.000
[50,] 60.0 1.000 1.000
> max(abs(ff1-ff2))
[1] 0.2505747

> ks.test(y1,y2)

        Two-sample Kolmogorov-Smirnov test

data:  y1 and y2 
D = 0.2506, p-value = 0.3127
alternative hypothesis: two.sided 

Warning message: 
cannot compute correct p-values with ties in: ks.test(y1, y2) 

###########################
# 2. Wilcoxon rank sum test
###########################
r<-rank(c(y1,y2))
r1<-sum(r[1:29])
> r1
[1] 976
r2<-sum(r[30:59])
> r1+r2
[1] 1770
> 59*60/2
[1] 1770

> help(wilcox.test)
> wilcox.test(y1,y2)

        Wilcoxon rank sum test with continuity correction

data:  y1 and y2 
W = 541, p-value = 0.1096
alternative hypothesis: true mu is not equal to 0 

> 2*(1-pnorm(1.6072))
[1] 0.1080105

### (approximated) exact test
> choose(59,29)
[1] 5.913229e+16  ## a huge number. We are going to take only a subset
                  ## of all the possible samples
n1<-length(y1)
n2<-length(y2)
p<-NULL
nsim<-10000

set.seed(250)

for (i in 1:nsim){
yy1<-sample(n1+n2,n1)
#yy2<-c(1:59)[is.na(match(1:59,yy1))]
p<-c(p,sum(yy1))}

> sum(p>=976)/nsim
[1] 0.0547
> sum(p<=976)/nsim
[1] 0.9482
> 0.0547*2
[1] 0.1094   ##p-value, similar to the large sample test


########################
# 3. Kruskal-Wallis test
########################
## Example 11.2.2.1, Woolson
## response: average number of pecks per unit time by pigeons
## group 1: amphetamine alone
## group 2: amphetamine plus low dose of chloropromazine
## group 3: amphetamine plus middle dose of chloropromazine
## group 4: amphetamine plus high dose of chloropromazine

yy1<-c(14.1,10.75,9.41,10.65,11.49,13.55,11.11,13.49)
yy2<-c(11,9.73,11.31,13.89,12.1,10.31,11.86)
yy3<-c(9.65,10.05,11.75,11.06,10.51)
yy4<-c(10.53,10.71,11.43,8.99,9.02,10.22,13.03,10.49,11.34)

## ranks
rr<-rank(c(yy1,yy2,yy3,yy4))
rr1<-rr[1:8]
rr2<-rr[9:15]
rr3<-rr[16:20]
rr4<-rr[21:29]

> mean(rr1)
[1] 18.625
> mean(rr2)
[1] 17.28571
> mean(rr3)
[1] 11.6
> mean(rr4)
[1] 11.88889

grp<-rep(c(1:4),c(8,7,5,9))
> kruskal.test(c(yy1,yy2,yy3,yy4)~grp)

        Kruskal-Wallis rank sum test

data:  c(yy1, yy2, yy3, yy4) by grp 
Kruskal-Wallis chi-squared = 3.9532, df = 3, p-value = 0.2666