Fisherの正確確率検定

1ページ　2ページ　3ページ

超幾何分布

P(X | x=12, m=15, n=14, k=18) の条件付き確率は、超幾何分布から特定の成功数が抽出される確率を計算します。

超幾何分布

超幾何分布は、有限個の要素から成る母集団より無作為に抽出したサンプルにおいて、特定の属性を持つ要素が何個含まれているかを示す確率変数の分布です。

y2pt.com

2023.01.01

上記の記事より

P(X=12)=$\dfrac{_{15}C_{12}\times_{14}C_{6}}{_{29}C_{18}}$=0.03949341

Rを使用したら、以下のようになります

choose(15, 12)*choose(14, 6)/choose(29, 18)

> choose(15, 12)*choose(14, 6)/choose(29, 18)
[1] 0.03949341

Fisherの正確検定のp値の求め方

Rの関数を使用して、P(X | x=12, m=15, n=14, k=18) を求めます。dhyperは超幾何分布の確率質量を求める関数です（”d”は通常 “density” を意味しますが、ここでは離散分布の確率質量を指します）。

dhyper関数のパラメータ説明
　x : 成功が観測された回数（”severe” and “yes” の数）
　m : 母集団における成功の数（”severe” の数）
　n : 母集団における失敗の数（”mild” の数）
　k : 抽出回数、サンプルサイズ（”yes” の数）
＊重度群に対する効果を検証するので、この例では、成功を”重度群”、また観察を”効果あり”として扱います。少し分かりにくいので、超幾何分布をご参照ください。

x <- 12  # Observed number of successes
m <- 15  # Total number of successes in the population
n <- 14  # Total number of failures in the population
k <- 18  # Size of the sample to be drawn
dhyper(x, m, n, k)

> dhyper(x, m, n, k)
[1] 0.03949341

片側検定（上側），alternative = “greater”　のｐ値

再掲

> fisher.test(tab, alternative = "greater")

        Fisher's Exact Test for Count Data

data:  tab
p-value = 0.04601
alternative hypothesis: true odds ratio is greater than 1
95 percent confidence interval:
 1.029929      Inf
sample estimates:
odds ratio 
  5.003853

sumは集計結果を算出する関数です。

x <- c(12:15)
m <- 15
n <- 14
k <- 18
sum(dhyper(x, m, n, k))

> sum(dhyper(x, m, n, k))
[1] 0.04601384

P(X≧12) = P(X=12) + P(X=13) + P(X=14) + P(X=15)

X＝12, 13, 14, 15 の確率を直接求めることになるので、なかなか大変です。なので、ｎが大きくなると昔は計算が大変だった思います。今はパソコンですぐに計算できます。

両側検定, alternative = “two.sided” のｐ値

再掲

fisher.test(tab, alternative = "two.sided")

> fisher.test(tab, alternative = "two.sided")

        Fisher's Exact Test for Count Data

data:  tab
p-value = 0.06043
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  0.8163182 40.7244582
sample estimates:
odds ratio 
  5.003853 

dhyper関数は指定されたxの値での超幾何分布の確率を計算します。ここではxが15から4の範囲で、各xに対する確率を求めます。phyperは累積分布関数です。qパラメータで指定された各点での累積確率を計算し、その値以下の確率の総和を返します。lower.tail=Fオプションで、qパラメータで指定された値よりも大きな値の確率の総和を求めています。puはq = 14:3で計算されているため、各xの値に対してその次の数以上の成功を得る確率を求めています。また、plは、下側の確率を求めます。

x <- c(15:4)
m <- 15
n <- 14
k <- 18

d  <- dhyper(x = 15:4, m, n, k)
pu  <- phyper(q = 14:3, m, n, k, lower.tail=F)
pl  <- phyper(q = 15:4, m, n, k)
df <- data.frame(x, d, pu, pl)
round(df,6)

観測されたオッズ比の出現確率と、それよりも低い確率を合計します。 周辺合計を固定して、セルの値を再配分して組合せを求めるます。黄色ハイライト部分の合計がfisher.testの両側検定のp値。

0.046014 + 0.014419 = 0.060433

d1 <- d[1:4]
d2 <- d[10:12]
sum(d1, d2)

> d1 <- d[1:4]
> d2 <- d[10:12]
> sum(d1, d2)
[1] 0.06043294

ここでは、両側検定と片側検定の累積確率の位置が同じなっていますが、以下のような例になるとズレが生じるので注意が必要です。

累積確率と異なる値になる場合

tab2 <- matrix(c(7, 2, 8, 10), 2, 2)
tab2

> tab2
     [,1] [,2]
[1,]    7    8
[2,]    2   10

x <- 7  # Observed number of successes
m <- 15  # Total number of successes in the population
n <- 12  # Total number of failures in the population
k <- 9  # Size of the sample to be drawn
dhyper(x, m, n, k)

> dhyper(x, m, n, k)
[1] 0.09061785

fisher.test(tab2, alternative = "two.sided")$p.value

> fisher.test(tab2, alternative = "two.sided")$p.value
[1] 0.2172387

X <- c(9:0)
m <- 15
n <- 12
k <- 9

d  <- dhyper(x = 9:0, m, n, k)
pu  <- phyper(q = 8:-1, m, n, k, lower.tail=F)
pl  <- phyper(q = 9:0, m, n, k)
df2 <- data.frame(X, d, pu, pl)
round(df2,6)

tabの例では確率と累積確率のXが同じ値になるのですが、このように異なる場合もあるので注意が必要です。

> 0.001068+0.016476+0.090618+0.089703+0.017743+0.001584+0.000047
[1] 0.217239
> 0.108162+0.019375
[1] 0.127537

念のためにRで確認します

sum(df2[c(1:3, 7:10), 2])
df2[3, 3] + df2[8, 4]

> sum(df2[c(1:3, 7:10), 2])
[1] 0.2172387
> df2[3, 3] + df2[8, 4]
[1] 0.1275362

黒木玄様（@genkuroki）よりご助言いただき、サンプルも使用させていただきました。ありがとうございました。