マン・ホイットニー検定（exact test） Last modified: Apr 13, 2006
目的
正確な P 値を計算するマン・ホイットニー検定（正確確率検定）である。データ（分割表）によっては計算量が多くなり実用的な時間内で計算が終了できないこともあるので，そのような場合にはモンテカルロ法による近似計算もできる。ネット上で計算すればもう少し計算時間は短い。
周辺和を固定した全ての分割表においてマン・ホイットニーの U 統計量と生起確率を求め，実際に観察された分割表の U 統計より小さいか等しい分割表の生起確率を合計したものが P 値であるとするものである。
注：exactRankTests パッケージには wilcox.exact 関数があり，同値のある場合にも正確な P 値を計算できる。
使用法
exact.mw(x, y=NULL, exact=TRUE, hybrid=FALSE, loop=10000)
引数
x 分割表（合計を含まない）または，第 1 群のデータベクトル
y x がデータベクトルの場合は，第 2 群のデータベクトル
x が分割表の場合には無視される
exact 正確な P 値を求める場合，またはシミュレーションにより近似的な P 値を求めるときには TRUE（デフォルト）
FALSE の場合にはカイ二乗分布による漸近近似のみを行う
hybrid TRUE を指定すれば，シミュレーションにより近似的な P 値を計算する
loop シミュレーションの回数
ソース
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/exact-mw.R", encoding="euc-jp")
# マン・ホイットニー検定（exact test）
exact.mw <- function( x, # 分割表（合計を含まない） もしくはデータベクトル
y=NULL, # x がデータベクトルのときは，データベクトル
exact=TRUE, # 正確検定を行うかどうか
hybrid=FALSE, # TRUE にすれば，シミュレーションによる
loop=10000) # シミュレーションの回数
{
found <- function() # 周辺度数が同じ分割表が見つかった
{
hh <-abs(temp2-sum(um[1,]*r)) # u.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)
}
u.test <- function(t) # マン・ホイットニーの U 検定
{
return(abs(temp2-sum(t[1,]*r)))
}
monte.carlo <- function() # モンテカルロ検定
{
printf("%i 回のシミュレーションによる P 値 = %g\n", loop, sum(sapply(r2dtable(loop, rt, ct), u.test) >= stat_val)/loop)
}
asymptotic <- function()
{
z <- stat_val/sqrt(n1n2/(n*(n-1))*(n^3-n-sum(ct^3-ct))/12)
printf("U = %g, Z = %g, P 値 = %g\n", n1n2/2-stat_val, z, pnorm(z, lower.tail=FALSE)*2)
}
if (is.matrix(x)) { # 分割表が与えられたとき
t <- x
}
else { # 2 変数が与えられたとき
nx <- length(x)
ny <- length(y)
t <- table(rep(1:2, c(nx, ny)), c(x, y))
}
EPSILON <- 1e-10
EXPMAX <- 1e100
log_expmax <- log(EXPMAX)
nr <- nrow(t) # 分割表の行数
stopifnot(nr==2) # 2×k 分割表でないといけない
nc <- ncol(t) # 分割表の列数
rt <- rowSums(t) # 行和
ct <- colSums(t) # 列和
n1 <- rt[1] # 第 1 群のケース数
n2 <- rt[2] # 第 2 群のケース数
n1n2 <- n1*n2
n <- n1+n2 # 全ケース数
r <- cumsum(c(0,ct[-nc]))+(ct+1)/2
temp <- n1n2+n1*(n1+1)/2
temp2 <- temp-n1n2/2
stat_val <- abs(temp2-sum(t[1,]*r)) # 観察された分割表に対する統計量
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, 2,3,1,2), byrow=TRUE, ncol=4)
exact.mw(x)
exact.mw(x, hybrid=TRUE)
2 つのデータベクトルが与えられる場合
x <- c(78, 64, 75, 45, 82)
y <- c(110, 70, 53, 51)
exact.mw(x, y)
wilcox.test(x, y, correct=FALSE)
出力結果例
分割表が与えられる場合
> x <- matrix(c(5,3,2,1, 2,3,1,2), byrow=TRUE, ncol=4)
> x
[,1] [,2] [,3] [,4]
[1,] 5 3 2 1
[2,] 2 3 1 2
> exact.mw(x)
U = 33.5, Z = 0.907299, P 値 = 0.364248 # 正規近似による P 値
正確な P 値 = 0.3837421608
査察した分割表の個数は 91 個
> exact.mw(x, hybrid=TRUE)
U = 33.5, Z = 0.907299, P 値 = 0.364248 # 正規近似による P 値
10000 回のシミュレーションによる P 値 = 0.3828
2 つのデータベクトルが与えられる場合
> x <- c(78, 64, 75, 45, 82)
> y <- c(110, 70, 53, 51)
> exact.mw(x, y)
U = 9, Z = 0.244949, P 値 = 0.806496 # 正規近似による P 値
正確な P 値 = 0.9047619048
査察した分割表の個数は 126 個
> wilcox.test(x, y, correct=FALSE) # R に用意されている関数を使う場合
Wilcoxon rank sum test
data: x and y
W = 11, p-value = 0.9048 # 正確な P 値
alternative hypothesis: true mu is not equal to 0
注：wilcox.test は，データ数の合計が 50 未満で，同値がない場合にのみ正確な P 値を計算するが，同値がある場合にはデータ数にかかわらず正規近似による P 値しか与えない。coin パッケージにある wilcox_test で，distribution="exact" を指定することにより，正確な P 値を計算できる。
> x <- c(78, 64, 75, 45, 82)
> y <- c(110, 70, 53, 53)
> wilcox.test(x, y, correct=FALSE)
Wilcoxon rank sum test
data: x and y
W = 11, p-value = 0.8057 # 正規近似による P 値である！
alternative hypothesis: true mu is not equal to 0
警告メッセージ:
In wilcox.test.default(x, y, correct = FALSE) :
タイがあるため、正確な p 値を計算することができません
> exact.mw(x, y)
U = 9, Z = 0.245976, P 値 = 0.805701 # 正規近似による P 値
正確な P 値 = 0.8492063492
査察した分割表の個数は 91 個
> z <- c(x, y)
> g <- factor(c(1, 1, 1, 1, 1, 2, 2, 2, 2))
> library(coin)
> wilcox_test(z~g, distribution="exact") # モデル式で指定する
Exact Wilcoxon Mann-Whitney Rank Sum Test
data: z by g (1, 2)
Z = 0.246, p-value = 0.8492 # 表示されている Z は正規近似のための Z 値である。P 値は正確な P 値
alternative hypothesis: true mu is not equal to 0
解説ページ(1)
解説ページ(2)
直前のページへ戻る E-mail to Shigenobu AOKI
