Fisher の exact test Last modified: Feb 24, 2006
目的
Fisher の正確確率検定を行う。R には fisher.test がある。この fisher 関数は，フィッシャー流の P 値の決め方(fisher.test と同じ）の他に，ピアソン流の P 値の決め方も選べる。また，計算量が多量になり実用的な時間内に計算が終了できないような場合に対応するために，モンテカルロ法に基づく近似検定も選択できる（fisher.test の hybrid オプションと同じ）。ネット上で計算すればもう少し計算時間は短い。
フィッシャー流の P 値の決め方とは，周辺和を固定した全ての分割表の生起確率をもとめ，実際に観察された分割表の生起確率より小さいか等しい分割表の生起確率を合計したものが P 値であるとするものである。
ピアソン流の P 値の決め方とは，周辺和を固定した全ての分割表においてピアソンのカイ二乗統計量と生起確率を求め，実際に観察された分割表のピアソンのカイ二乗統計より小さいか等しい分割表の生起確率を合計したものが P 値であるとするものである。
本来，Fisher の exact test の目的は，独立性の検定にあると思われ，独立性からのずれを評価する統計指標は様々ある。しかるに，フィッシャー流の P 値の決め方は，基盤とする統計量を用いていない。単に生起確率が観察された分割表の生起確率より小さいことのに基づくフィッシャー流の P 値の決め方は不適切であると考えるものである。
使用法
fisher(x, y=NULL, exact=TRUE, method=c("Fisher", "Pearson"), hybrid=FALSE, loop=10000)
引数
x 分割表（合計を含まない） もしくは factor ベクトル
y x が factor ベクトルのときは，factor ベクトル
x が分割表の時には無視される
exact 正確な P 値を求める場合，またはシミュレーションにより近似的な P 値を求めるときには TRUE（デフォルト）
FALSE の場合にはカイ二乗分布による漸近近似のみを行う
method P 値の決め方。Fisher か Pearson を指定する。省略した場合には Fisher が仮定される
hybrid TRUE を指定すれば，シミュレーションにより近似的な P 値を計算する
loop シミュレーションの回数
ソース
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/fisher.R", encoding="euc-jp")
# Fisher の exact test
fisher <- function( x, # 分割表（合計を含まない） もしくは factor ベクトル
y=NULL, # x が factor ベクトルのときは，factor ベクトル
exact=TRUE, # 正確検定を行うかどうか
method=c("Fisher", "Pearson"), # フィッシャーによるか，ピアソンによるか
hybrid=FALSE, # TRUE にすれば，シミュレーションによる
loop=10000) # シミュレーションの回数
{
found <- function() # 周辺度数が同じ分割表が見つかった
{
if (method == "Fisher") { # フィッシャーの基準による
nprod <- temp-sum(perm_table[um+1])
if (nprod <= criterion+EPSILON) {
nntrue <<- nntrue+exp(nprod-nntrue2*log_expmax)
while (nntrue >= EXPMAX) {
nntrue <<- nntrue/EXPMAX
nntrue2 <<- nntrue2+1
}
}
}
else { # ピアソンの基準による
hh <- sum((um-ex)^2/ex) # chisq.test(um)
if (hh >= stat_val || abs(hh-stat_val) <= EPSILON) {
nprod <- temp-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("%s の方法による，正確な P 値 = %.10g\n", method, nntrue/denom*EXPMAX^(nntrue2-denom2))
printf("査察した分割表の個数は %s 個\n", ntables)
}
chisq.test <- function(t) # カイ二乗近似検定
{
return(sum((t-ex)^2/ex))
}
prob <- function(t)
{
return(temp-sum(perm_table[t+1]))
}
monte.carlo <- function() # モンテカルロ検定
{
if (method == "Fisher") { # フィッシャーの基準による
count <- sum(sapply(r2dtable(loop, rt, ct), prob) <= criterion+EPSILON)
}
else { # ピアソンの基準による
count <- sum(sapply(r2dtable(loop, rt, ct), chisq.test) >= stat_val)
}
printf("%i 回のシミュレーション（%s の方法）による P 値 = %g\n", loop, method, count/loop)
}
asymptotic <- function()
{
printf("カイ二乗値 = %g, 自由度 = %i, P 値 = %g\n", stat_val, df, pchisq(stat_val, df, lower.tail=FALSE))
}
if (is.matrix(x)) { # 分割表が与えられたとき
t <- x
}
else { # 2 変数が与えられたとき
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) # 総和
ex <- outer(rt, ct)/n # 期待値
stat_val <- chisq.test(t) # 観察された分割表のカイ二乗値
df <- (nr-1)*(nc-1) # 自由度
asymptotic() # 検定結果を出力
if (exact) { # 正確検定を行うなら，
method <- match.arg(method) # ピアソンの基準かフィッシャーの基準か
perm_table <- cumsum(c(0, log(1:(n+1)))) # 対数を取った階乗の表
temp <- sum(perm_table[rt+1])
criterion <- temp-sum(perm_table[t+1])
if (hybrid) { # モンテカルロ法による検定
monte.carlo()
}
else { # 正確な検定
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)
fisher(x)
fisher(x, method="Pearson")
fisher(x, hybrid=TRUE)
fisher(x, method="Pearson", hybrid=TRUE)
2 つの factor ベクトルを与える場合
x <- c(1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 3, 2, 2, 3, 1)
y <- rep(c("A", "B", "C"), each=5)
table(y, x)
fisher(x, y)
fisher.test(x)
出力結果例
分割表を与える場合
> 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
> fisher(x)
カイ二乗値 = 3.39631, 自由度 = 6, P 値 = 0.757711
Fisher の方法による，正確な P 値 = 0.8091124268
査察した分割表の個数は 24871 個
> fisher(x, method="Pearson")
カイ二乗値 = 3.39631, 自由度 = 6, P 値 = 0.757711
Pearson の方法による，正確な P 値 = 0.7878188077
査察した分割表の個数は 24871 個
> fisher(x, hybrid=TRUE)
カイ二乗値 = 3.39631, 自由度 = 6, P 値 = 0.757711
10000 回のシミュレーション（Fisher の方法）による P 値 = 0.8057
> fisher(x, method="Pearson", hybrid=TRUE)
カイ二乗値 = 3.39631, 自由度 = 6, P 値 = 0.757711
10000 回のシミュレーション（Pearson の方法）による P 値 = 0.7976
> fisher.test(x) # R に用意されている関数
計数データにおけるフィッシャーの正確確率検定
データ: x
P値 = 0.8091
対立仮説: 等しくない
2 つの factor ベクトルを与える場合
> x <- c(1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 3, 2, 2, 3, 1)
> y <- rep(c("A", "B", "C"), each=5)
> table(y, x)
x
y 1 2 3
A 2 2 1
B 1 3 1
C 1 2 2
> fisher(x, y)
カイ二乗値 = 1.28571, 自由度 = 4, P 値 = 0.863795
Fisher の方法による，正確な P 値 = 1
査察した分割表の個数は 180 個
解説ページ(1)
解説ページ(2)
直前のページへ戻る E-mail to Shigenobu AOKI
