kw.test <- function(data, cp1, cp2)
{
r <- rank(x <- data[,cp1])
n = sum(ni<- table(g <- data[,cp2]))
Ri <- tapply(r, g, sum)
S <- 12*sum(Ri^2/ni)/(n*(n+1))-3*(n+1)
if (length(unique(x)) != n) {
tie <- table(x)
S <- S/(1-sum(tie^3-tie)/(n^3-n))
}
df <- length(ni)-1
p <- pchisq(S, df, low=F)
a.m <- (n+1)/2
R.m <- Ri/ni
V <- sum((r-a.m)^2)/(n-1)
comp <- outer(R.m, R.m, function(x, y) { (x-y)^2 })/outer(ni, ni, function(n1, n2) { (1/n1+1/n2)*V })
pm <- pchisq(comp, df, lower=F)
result <- c(S, df, p)
names(result) <- c("Statistics(Chi-sq. value)", "d.f.", "P value")
list(result=result, comparison=comp, pvalue=pm)
}
# 使用例（データ）
data <- matrix(c(
3.42, 1,
3.84, 1,
3.96, 1,
3.76, 1,
3.17, 2,
3.63, 2,
3.47, 2,
3.44, 2,
3.39, 2,
3.64, 3,
3.72, 3,
3.91, 3
), byrow=TRUE, ncol=2)
# 使用例
kw.test(data, 1, 2)
第一引数はデータ行列
第二引数は測定値のある列番号
第三引数はグループ識別値のある列番号
$result （検定結果）
Statistics(Chi-sq. value) d.f. P value
5.54871795 2.00000000 0.06238946
$comparison （対比較の統計量）
1 2 3
1 0.000000000 4.104274 0.003663004
2 4.104273504 0.000000 3.702564103
3 0.003663004 3.702564 0.000000000
$pvalue （対比較のP値）
1 2 3
1 1.0000000 0.1284601 0.9981702
2 0.1284601 1.0000000 0.1570357
3 0.9981702 0.1570357 1.0000000