クラスカル・ウォリス検定（exact test） Last modified: Feb 24, 2006
目的
正確な P 値を与えるクラスカル・ウォリス検定（正確確率検定）である。データ（およびそれから導かれる分割表）によっては計算量が多くなり実用的な時間内で計算が終了できないこともあるので，そのような場合にはモンテカルロ法による近似計算を用いる必要があるかもしれない。ネット上で計算すればもう少し計算時間は短い。
周辺和を固定した全ての分割表においてクラスカル・ウォリス検定の統計量と生起確率を求め，実際に観察された分割表の 統計量より小さいか等しい分割表の生起確率を合計したものが P 値であるとするものである。
使用法
exact.kw(x，y=NULL, exact=TRUE, hybrid=FALSE, loop=10000)
引数
x 分割表（合計を含まない） もしくは データベクトルもしくはリスト（使用例参照）
y x がデータベクトルのときは，そのデータがどの群のものかを表す factor ベクトル
x が分割表やリストの時には無視される
exact 正確な P 値を求める場合，またはシミュレーションにより近似的な P 値を求めるときには TRUE（デフォルト）
FALSE の場合にはカイ二乗分布による漸近近似のみを行う
hybrid TRUE を指定すれば，シミュレーションにより近似的な P 値を計算する
loop シミュレーションの回数
ソース
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/exact-kw.R", encoding="euc-jp")
# クラスカル・ウォリス検定（exact test）
exact.kw <- function( x, # 分割表（合計を含まない） もしくはデータベクトル
y=NULL, # x がデータベクトルのときは，factor ベクトル
exact=TRUE, # 正確検定を行うかどうか
hybrid=FALSE, # TRUE にすれば，シミュレーションによる
loop=10000) # シミュレーションの回数
{
found <- function() # 周辺度数が同じ分割表が見つかった
{
hh <- sum((um%*%q)^2/rt) # kw_test(um)
if (hh >= stat_val || abs(hh-stat_val) <= EPSILON) {
nprod <- sum(perm_table[rt+1])-sum(perm_table[um+1])
nntrue <<- nntrue+exp(nprod-nntrue2*log_expmax)
while (nntrue >= EXPMAX) {
nntrue <<- nntrue/EXPMAX
nntrue2 <<- nntrue2+1
}
}
ntables <<- ntables+1
}
search <- function(x, y) # 分割表の探索
{
if (y == 1) { # 見つかった
found()
}
else if (x == 1) {
if ((t <- um[1, 1]-um[y, 1]) >= 0) {
um[1, 1] <<- t
search(nc, y-1)
um[1, 1] <<- um[1, 1]+um[y, 1]
}
}
else {
search(x-1, y)
while (um[y, 1] && um[1, x]) {
um[y, x] <<- um[y, x]+1
um[y, 1] <<- um[y, 1]-1
um[1, x] <<- um[1, x]-1
search(x-1, y)
}
um[y, 1] <<- um[y, 1]+um[y, x]
um[1, x] <<- um[1, x]+um[y, x]
um[y, x] <<- 0
}
}
exact.test <- function() # 正確検定
{
denom2 <- 0
denom <- perm_table[n+1]-sum(perm_table[ct+1])
while (denom > log_expmax) {
denom <- denom-log_expmax
denom2 <- denom2+1
}
denom <- exp(denom)
um[,1] <<- rt
um[1,] <<- ct
search(nc, nr)
printf("正確な P 値 = %.10g\n", nntrue/denom*EXPMAX^(nntrue2-denom2))
printf("査察した分割表の個数は %s 個\n", ntables)
}
kw.test <- function(u) # クラスカル・ウォリス検定
{
return(sum((u%*%q)^2/rt))
}
monte.carlo <- function() # モンテカルロ検定
{
printf("%i 回のシミュレーションによる P 値 = %g\n", loop, sum(sapply(r2dtable(loop, rt, ct), kw.test) >= stat_val)/loop)
}
asymptotic <- function()
{
chisq <- (stat_val*12/(n*(n+1))-3*(n+1))/(1-sum(ct^3-ct)/(n^3-n))
printf("カイ二乗値 = %g, 自由度 = %i, P 値 = %g\n", chisq, nr-1, pchisq(chisq, nr-1, lower.tail=FALSE))
}
if (is.list(x)) {
y <- factor(rep(1:length(x), sapply(x, length)))
t <- table(y, unlist(x))
}
else if (is.matrix(x)) {
t <- x
}
else {
t <- table(y, x)
}
EPSILON <- 1e-10
EXPMAX <- 1e100
log_expmax <- log(EXPMAX)
nr <- nrow(t) # 分割表の行数
nc <- ncol(t) # 分割表の列数
rt <- rowSums(t) # 分割表の行和
ct <- colSums(t) # 分割表の列和
n <- sum(t) # 総和
q <- cumsum(c(0,ct[-nc]))+(ct+1)*0.5
half <- (n+1)*0.5
stat_val <- kw.test(t) # クラスカル・ウォリス検定統計量
asymptotic() # 検定結果を出力
if (exact) {
if (hybrid) { # モンテカルロ法による検定
monte.carlo()
}
else { # 正確な検定
perm_table <- cumsum(c(0, log(1:(n+1))))
ntables <- nntrue <- nntrue2 <- 0
um <- matrix(0, nr, nc)
exact.test()
}
}
}
使用例
分割表を与える場合
x <- matrix(c(5,3,2,1, 4,3,5,2, 2,3,1,2), byrow=TRUE, ncol=4)
exact.kw(x)
exact.kw(x, hybrid=TRUE)
データベクトルと factor ベクトルを与える場合
data <- c(
3.42, 3.84, 3.96, 3.76,
3.17, 3.63, 3.47, 3.44, 3.39,
3.64, 3.72, 3.91
)
group <- rep(1:3, c(4, 5, 3))
exact.kw(data, group)
リストで与える場合
exact.kw(list(c(3.42, 3.84, 3.96, 3.76), c(3.17, 3.63, 3.47, 3.44, 3.39), c(3.64, 3.72, 3.91)))
出力結果例
分割表を与える場合
> x <- matrix(c(5,3,2,1, 4,3,5,2, 2,3,1,2), byrow=TRUE, ncol=4)
> x
[,1] [,2] [,3] [,4]
[1,] 5 3 2 1
[2,] 4 3 5 2
[3,] 2 3 1 2
> exact.kw(x)
カイ二乗値 = 1.32485, 自由度 = 2, P 値 = 0.5156 # カイ自乗近似による P 値
正確な P 値 = 0.5268191237
査察した分割表の個数は 24871 個
> exact.kw(x, hybrid=TRUE)
カイ二乗値 = 1.32485, 自由度 = 2, P 値 = 0.5156 # カイ自乗近似による P 値
10000 回のシミュレーションによる P 値 = 0.525
データベクトルと factor ベクトルを与える場合
> data <- c(
+ 3.42, 3.84, 3.96, 3.76,
+ 3.17, 3.63, 3.47, 3.44, 3.39,
+ 3.64, 3.72, 3.91
+ )
> group <- rep(1:3, c(4, 5, 3))
> exact.kw(data, group)
カイ二乗値 = 5.54872, 自由度 = 2, P 値 = 0.0623895 # カイ自乗近似による P 値
正確な P 値 = 0.0538961039
査察した分割表の個数は 27720 個
リストで与える場合
> exact.kw(list(c(3.42, 3.84, 3.96, 3.76), c(3.17, 3.63, 3.47, 3.44, 3.39), c(3.64, 3.72, 3.91)))
カイ二乗値 = 5.54872, 自由度 = 2, P 値 = 0.0623895 # カイ自乗近似による P 値
正確な P 値 = 0.0538961039
査察した分割表の個数は 27720 個
解説ページ(1)
解説ページ(2)
直前のページへ戻る E-mail to Shigenobu AOKI
