Fisherの正確確率
クロス集計表において周辺度数を固定すると、そこから生じうるすべてのセルの値の組合せから任意のクロス集計表の生起確率が計算できる。それをFisherの正確確率という。
正確確率と言うだけあって、Peasonの漸近有意確率よりも「正確」である。ならばこちらを使えばいいやん、というのは全く正しい。ただしこの計算が中々に厄介である。なにしろすべての組合せを数え上げなければならない。組合せ計算には階乗が登場するが、これが時にとてつもなく大きな値となる。そしてコンピュータにも扱える数字には限界がある(例えばExcelは171の階乗は計算できない)。
サンプルサイズがさほど大きくない2*2のクロス集計表ならExcelでも(頑張れば)たぶん計算可能である(Excelでは組合せ関数はCONBIN(総数,抜き取り数))。
Rでは組み合わせはchoose関数を用いる。
実測値 | No | Yes | 総計 |
---|---|---|---|
female | 149 | 107 | 256 |
male | 120 | 124 | 244 |
総計 | 269 | 231 | 500 |
組合せ | 度数 | |
---|---|---|
female | choose(256,149) | 1.836192e+74 |
male | choose(244,120) | 1.396167e+72 |
合計 | choose(500,269) | 2.759088e+148 |
確率 | 0.009291591←女度数*男度数/合計度数 |
女Noを減らす。
可能性1 | No | Yes | 総計 |
---|---|---|---|
female | 148 | 108 | 256 |
male | 121 | 123 | 244 |
総計 | 269 | 231 | 500 |
組合せ | 度数 | |
---|---|---|
female | choose(256,148) | 2.533265e+74 |
male | choose(244,121) | 1.430783e+72 |
合計 | choose(500,269) | 2.759088e+148 |
確率 | 0.01313678←実測値より多いので採択して継続する |
女Noを減らす
可能性2 | No | Yes | 総計 |
---|---|---|---|
female | 147 | 109 | 256 |
male | 122 | 122 | 244 |
総計 | 269 | 231 | 500 |
組合せ | 度数 | |
---|---|---|
female | choose(256,147) | 3.439663e+74 |
male | choose(244,122) | 1.442511e+72 |
合計 | choose(500,269) | 2.759088e+148 |
確率 | 0.0179833←実測値より多いので採択して継続する |
女Noを減らす
可能性3 | No | Yes | 総計 |
---|---|---|---|
female | 146 | 110 | 256 |
male | 123 | 121 | 244 |
総計 | 269 | 231 | 500 |
組合せ | 度数 | |
---|---|---|
female | choose(256,146) | 4.596641e+74 |
male | choose(244,123) | 1.430783e+72 |
合計 | choose(500,269) | 2.759088e+148 |
確率 | 0.02383685←実測値より多いので採択して継続する |
…女Noを減らす
可能性23 | No | Yes | 総計 |
---|---|---|---|
female | 126 | 130 | 256 |
male | 143 | 101 | 244 |
総計 | 269 | 231 | 500 |
組合せ | 度数 | |
---|---|---|
female | choose(256,126) | 5.59185e+75 |
male | choose(244,143) | 3.871852e+70 |
合計 | choose(500,269) | 2.759088e+148 |
確率 | 0.007847089←実測値より少ないので破棄、終了 |
女Noを増やす
可能性24 | No | Yes | 総計 |
---|---|---|---|
female | 150 | 106 | 256 |
male | 119 | 125 | 244 |
総計 | 269 | 231 | 500 |
組合せ | 度数 | |
---|---|---|
female | choose(256,150) | 1.309817e+74 |
male | choose(244,119) | 1.340321e+72 |
合計 | choose(500,269) | 2.759088e+148 |
確率 | 0.006362882←実測値より少ないので破棄、終了 |
この確率の総和が実測値より生起確率が高い組合せ全体の生起確率となる。
p値はその逆(実測値以下の部分の面積)だから1からこの確率を引いた値がp値(両側)となる。なお正確確率においては片側は両側を2でわった値ではない。片側p値は観測値以下(または以上)の生起確率をすべて足したものである。
p.value ← 1 - sum(実測値の生起確率よりも大きい確率)
2*2(自由度1)のクロス集計表で正確確率を求める関数サンプル
fisher.test.sample <- function(crossTable,alternative=c("two.sided","greater","less")){ alternative <- match.arg(alternative) N <- sum(crossTable) if((nrow(crossTable) - 1)*(ncol(crossTable) - 1)==1){ #観測値の生起確率 prob <- choose(apply(crossTable,1,sum)[[1]],crossTable[1])*choose(apply(crossTable,1,sum)[[2]],crossTable[2])/choose(N,apply(crossTable,2,sum)[[1]]) if(is.nan(prob))return(NA) #度数が大きすぎて計算できない時は終了 prob.sum <- 0 #この変数に合計確率を計算していく prob.tmp <- 1 #試行ごとの生起確率、これを観測値の生起確率と比較する x <- 0 if(alternative != "less"){ #1行1列目の観測値を一つずつ減らしていき、生起確率を計算する while(TRUE){ x <- x + 1 prob.tmp <- choose(apply(crossTable,1,sum)[[1]],crossTable[1] - x)*choose(apply(crossTable,1,sum)[[2]],crossTable[2] + x)/choose(N,apply(crossTable,2,sum)[[1]]) if(crossTable[1] - x < 0 || prob.tmp == 0) break if(alternative != "greater"){ #実測値の確率を下回った時点で終了 if(prob.tmp < prob) break } #実測値を上回る生起確率の時は合計確率に加えていく prob.sum <- prob.sum + prob.tmp } } x <- 0 if(alternative != "greater"){ #1行1列目の観測値を一つずつ増やしていき、生起確率を計算する while(TRUE){ x <- x + 1 prob.tmp <- choose(apply(crossTable,1,sum)[[1]],crossTable[1] + x)*choose(apply(crossTable,1,sum)[[2]],crossTable[2] - x)/choose(N,apply(crossTable,2,sum)[[1]]) if(crossTable[2] - x < 0 || prob.tmp == 0) break if(alternative != "less"){ if(prob.tmp < prob) break } prob.sum <- prob.sum + prob.tmp } } p.value <- 1 - prob.sum #1から合計確率を引いた値がp値 }else return(NA) return (p.value) }
- 2*2クロス表限定
- 度数が大きくなると計算できない。
- 計算中p値の小数を丸めてしまうので正確な計算が出来ない(正確確率じゃない!)
あくまでFishet正確確率計算の原理を理解するためのサンプルプログラムなので、実際の運用にはRに初期から入っているfisher.test関数を用いること。ただしこの関数でも度数が大きい時は計算不能。
使用例
ここで公開しているfisher.test.sample関数は途中経過を表示する。
#fisher.test.sample関数の読み込み source("http://kyoto-edu.sakura.ne.jp/weblesson/statistics/fisher.test.R", encoding="UTF-8") #クロス集計をする元データの読み込み data <- read.csv("http://kyoto-edu.sakura.ne.jp/weblesson/statistics/data/gender-opinion.csv", fileEncoding = "utf-8") gender <- data$性別 opinion <- data$Aという意見に #クロス集計表作成 crossTable <- table(gender,opinion) #作成したクロス表のFisherの正確確率を計算する fisher.test.sample(crossTable) #両側 fisher.test.sample(crossTable,alternative="less") #片側(<) fisher.test.sample(crossTable,alternative="greater") #片側(>)