18/05/06 22:15:43.42 BK1CxH7U.net
# jonckheereテストを書いてみた
jonckheere <- function(L,
alternative = c("two.sided", "increasing", "decreasing"),
cat=TRUE){
# L : list of vectors A1,A2,...,Ak, with assumed tendency
How.Many.Greater.Pairs <- function(A,B){ # How many pairs of A[i] > B[j], count as 0.5 when equal,
n.a = length(A)
n.b = length(B)
how.many.greater.pairs = 0
for(i in 1:n.a){
for(j in 1:n.b){
how.many.greater.pairs = how.many.greater.pairs+ifelse(A[i]==B[j],0.5,A[i]>B[j])
}
}
return(how.many.greater.pairs)
}
Sum.of.Greater.Pairs <- function(L){ #L=list(A1,,,,Ak),A1 < A2 < A3,..,< Ak : vector
k = length(L)
comb = combn(1:k,2) # possible combinaition of pairs to compare
n.comb = ncol(comb) # how many combinations
J = 0 # sum of greater pairs
for(i in 1:n.comb){
J = J + How.Many.Greater.Pairs(L[[comb[1,i]]],L[[comb[2,i]]])
}
return(J)
}
J = Sum.of.Greater.Pairs(L)
n = sapply(L,length)
N = sum(n)
EJ = (N^2-sum(n^2))/4
VJ = (N^2*(2*N+3)-sum(n^2*(2*n+3)))/72
Z = (J-EJ)/sqrt(VJ)
alternative = match.arg(alternative)
p.value = switch(alternative, 'two.sided' = 2 * min(pnorm(Z), pnorm(-Z), 0.5),
'increasing' = pnorm(Z),
'decreasing' = pnorm(-Z))
if(cat){
cat( 'p.value = ', p.value,'\n')
cat('alternative hypothesis: ' ,alternative,'\n')
}
invisible(p.value)
}