Fisherの正確確率検定

僕が作った架空のデータで練習しましょう。結果についてはご意見があると思いますが、統計学の練習のためのサンプルですのでご了承ください。

例題)新たなトレーニング方法(Z法)が開発されました。果たしてこのトレーニングは運動障害が重度群と軽度群で効果に差があるのか?

データの準備と要約

Googleでのヒット数(2024.1.24調べ)

Fisherの正確確率検定-> 9200件
Fisherの正確検定->5220件
Fisherの直接確率検定-> 4230件
Fisherの直接法-> 3220件

Fisherの正確検定、Fisherの正確確率検定、Fisherの直接確率検定、Fisherの直接法などと呼ばれています(統一してくれれば良いのにといつも思います)。Fisherの正確確率検定が多かったので、このタイトルに使ってます。

定義

\(\displaystyle オッズ=(\frac{効果あり}{効果なし})\)

  • ここではオッズを効果率と定義します
  • このようにオッズやオッズ比の定義を確認しておくことが重要です

サンプル(架空のサンプルです)

効果あり効果なし
重度群12315
軽度群6814
合計181129

重度群の効果率は4、軽度群の効果率は0.75となります

この効果率の比が「1」か「1ではない」かということを調べるために、下記のオッズ比を利用します

\(\displaystyle (オッズ比 =)(\dfrac{重度群の効果率}{軽度群の効果率})(= 5.33)\)

5.33なのでオッズ比は1より大きいようですが・・・

オッズ比が高くなると重度群の効果率が高くなることを意味しており、オッズ比が1より大きくなると軽度群の効果率より重度群の効果率が高いことを意味します

この検定の仮説は以下のようになります

帰無仮説H0: オッズ比=1(重度群の効果率 = 軽度群の効果率
対立仮説H1: オッズ比1(重度群の効果率 ≠ 軽度群の効果率

Rで分割表の作成

tabという単語の中にサンプルの分割表を挿入してみます

tab <- matrix(c(12, 6, 3, 8), 2, 2)
print(tab)

これで「tab」という単語(なんでも構いません)の中に分割表が入りました

ベクトルから作成すると以下のようになります

障害 <- c(rep("重症", 12), rep("軽症", 6), rep("重症", 3), rep("軽症", 8))
障害 <- factor(障害)
障害 <- relevel(障害 , ref=c("重症")) 
効果 <- c(rep("あり", 18), rep("なし", 11))
効果 <- factor(効果)

xtabs(
    ~障害 + 効果
)

Fisherの正確検定では、「効果あり=18」となる12パターンの確率を求めます
つまり「効果あり=18」を条件とした場合の確率です

これでFisherの正確検定の準備は完了、有意水準5%で検定してみましょう

研究計画(何を立証したいのか)に従い以下から適切な検定方法を選択します

Rの関数で検定

tabには上記の分割表が入っています

fisher.test(tab, alternative = "less")

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

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

ごちゃごちゃしてますが、以下のように出力されれば成功です

上記結果の3つの違いについて

x:重度群かつ効果あり、効果あり=18は固定
つまり、もとめる確率は、P(X | X=x, 効果あり=18)という条件付き確率となります 

fisher.test(tab, alternative = “less”)

X≦12となる確率を求めます
X=12, X=11, X=10, ・・・, X=4となる確率
つまり効果率がX≦12になる場合の確率です

P(X | x≦12, 効果あり=18) 

\(P(X=12) + P(X=11) + ・・・ + P(X=4)=0.9935\)

*この確率は超幾何分布を使って求めます (後で出てきます)

p値は0.9935なので重度群と軽度群の効果率の差はありません

fisher.test(tab, alternative = “greater”)

X≧12となる確率を求めます
X=12, x=13, X=14, X=15となる確率
つまり効果率がX≧12になる場合の確率です

P(X | x≧12, 効果あり=18) 

\(P(X=12) + P(X=13) + P(X=14) + P(X=15)=0.04601384\)

p-value = 0.04601384
有意水準が0.5%の場合
X が12以上になるのは稀であることを意味します
つまり、効果率は1以上と考え重症群の効果率が有意に高いと判断します

fisher.test(tab, alternative = “two.sided”)

P=0.06043 なので有意差はありません
このp値は後半で説明します

fisher.test(tab, alternative = “greater”) の結果の説明
  1. p-value(p値)
    0.04601384
  2. alternative hypothesis(対立仮説)
    Rは対立仮説を出力します
    真のオッズ比は1ではない
    すなわち・・・p値<0.05なので「重症群と軽度群には有意差があり、
    AM実施群の効果率は有意に高いという結果になります
    注意)両側検定の場合P値が大きくなり有意ではなくなります
  3. 95 percent confidence interval
    下記のオッズ比の95%信頼区間
  4. odds ratio(オッズ比)5.003853
    注意
    これは単純な割り算のオッズ比(5.33)ではありません
    条件付き最尤推定から算出されています(効果あり=18)
    詳細は下記の奥村晴彦先生のページをご参照ください
    https://oku.edu.mie-u.ac.jp/~okumura/stat/fishertest.html

超幾何分布

超幾何分布から\(P(X)\)を算出します

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

\(P(X=12)=0.03949341\)

ということが分かりました
つまり数あるパターンの中から上の分割表のパターンになる確率のことです

Rでやってみましょう
dhyperは幾何分布の確率密度を求めるための関数です(”d” gives the probability density)

# dhyperは超幾何分布の確率密度関数

a <- 15
b <- 14
n <- 18
x <- 12

dhyper(
    x = x, m = a, n = b, k = n
    )

以下のように表示されていれば成功です

さきほど求めた解と同じ確率が算出されます

X≧12やX≦12となる確率はどうなるでしょうか?

X=0~15の確率を直接求めるのです!
全て足せばよいのでやってみましょう
$$ P(X≧12) = P(X=12) + P(X=13) + P(X=14) + P(X=15) $$
なかなか大変です・・・
なのでnが大きくなると昔は計算が大変だった思います・・・
でも今はRですぐに計算できるので感謝です

fisher.test(tab, alternative = “greater”) のp値
# sumは集計結果を算出する関数です
a <- 15
b <- 14
n <- 18
x <- c(12:15)

sum(
    dhyper(x = x, m = a, n = b, k = n)
    )

以下のように表示されていれば成功です

fisher.test(tab, alternative = “two.sided”) のp値

Rでやってみましょう

tab <- matrix(c(12, 6, 3, 8), 2, 2)
tab
fisher.test(tab, alternative = "two.sided")
a <- 15
b <- 14
n <- 18
X <- c(15:4)

確率  <- dhyper(x = 15:4, m = a, n = b, k = n)
累積1  <- phyper(q = 14:3, m = a, n = b, k = n, lower.tail=F)
累積2  <- phyper(q = 15:4, m = a, n = b, k = n)
df <- data.frame(X, 確率, 累積1, 累積2)
round(df,6)

赤塗り部分の合計がfisher.testのp値

sum(df[c(1:4,10:12),2])

累積確率で考えたら間違いますので注意が必要です

例1:これまで使用してきた例では累積確率と同じ値になる場合

赤の合計と青の合計が同じ値になります

#赤部分=確率の和
sum(df[c(1:4,10:12),2])
#青部分=累積確率の和
df[4, 3]+df[10, 4]

累積確率と異なる値になる場合もあるので注意が必要です

例2:累積確率と異なる値になる

tab2 <- matrix(c(7, 2, 8, 10), 2, 2)
tab2
fisher.test(tab2, alternative = "two.sided")

赤の合計と青の合計が異なる値になります

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

library(exact2x2)

サンプル(tab3)を使用して3つのp値についてまとめます

tab3 <- matrix(c(16, 4, 4, 6), 2, 2)
tab3
a <- 20
b <- 10
n <- 20
X <- c(19:10)
確率  <- dhyper(x = 19:10, m = a, n = b, k = n)
累積1  <- phyper(q = 18:9, m = a, n = b, k = n, lower.tail=F)
累積2  <- phyper(q = 19:10, m = a, n = b, k = n)
df3 <- data.frame(X, 確率, 累積1, 累積2)
round(df3,6)
tsmethod=”minlike”

観測されたクロス表の確率以下の確率を全て合計したもの(Sterne の方法)

fisher.testで算出されるP値です

exact2x2(tab3, tsmethod="minlike")
sum(df3[c(1:4,10), 2])
tsmethod=”central”

片側(小さい方) p 値の 2 倍

exact2x2(tab3, tsmethod="central")

赤塗り部分の合計の2倍(=青塗り部分の2倍)

df3[4, 3]*2
tsmethod=”blaker”

観測された裾の確率に、裾の確率よりも大きくない反対側の裾の最大の裾の確率を加えた合計

exact2x2(tab3, tsmethod="blaker")

赤塗り部分の合計(この例ではminlikeと同じ値になりますが書いておきます)

sum(df3[c(1:4,10), 2])

p値と信頼区間の整合性

fisher.test では、p 値は tsmethod=“minlike” 、信頼区間はtsmethod=“central” に基づいて求められており、整合性に問題が生じます

上記の例teb3をfisher.testで検定した場合

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

p<0.05なので信頼区間が1を跨いでしまっています

ですので上述したtsmethod=”minlike”で実行することでp値と信頼区間の整合性が担保されます

間違った解釈をしないように!

再掲

効果あり効果なし
重度群12315
軽度群6814
合計181129

AM実施群が有効であることを証明したい場合

帰無仮説H0AM群の効果率 = PM群の効果率
対立仮説H1AM群の効果率 > AM群の効果率

上記のような実験結果から試験群の有効率が対照群の有効率より大きいことを証明します.
Fisherの直接法により有意水準5%で検定しましょう.

tab <- matrix(c(12, 6, 3, 8), 2, 2)
fisher.test(tab)

p=0.06となり有意差はありません・・・
という回答は間違っています.

fisher.testのデフォルトを確認しましょう

?fisher.test

alternative = “two.sided”
両側検定が初期設定として機能しています.
しかし今回の検定は 「片側検定を5%水準で行う」と宣言しています.
alternative = “lessまたは“greater”を設定しなければなりません.
さてどちらを選択すればよいでしょうか・・・?

  • less

先に述べましたようにX≦12となる確率算出します.
p-value = 0.9935

Rは分割表の左上を基準にしています

たまたまAM実施群の効果ありが左上にあるので、そこが基準となっています

それでは左列が効果なしなっている分割表に換えてみましょう

# c(2, 1) 2列目と1列目を入れ替えるという指示です
tab2 <- tab[, c(2, 1)] 
tab2

これで列が入れ替わりました

検定してみましょう

fisher.test(
    tab2, alternative = "less"
    )

fisher.test(
tab, alternative = “less”
)


のときはp-value = 0.9935でしたが

fisher.test(
tab2, alternative = “less”
)

ではp-value = 0.04601となっています.

なので…基準を間違うと検定結果も間違うということになります.


今回は「 AM実施が有効 > PM実施が有効」を証明したいので・・・

tab2のまま検定するのであれば、

fisher.test(
tab2, alternative = “less
)

tabで検定するのであれば

fisher.test(
tab, alternative = “greater
)

ということになります

もちろんどちらも同じp値になります
*ただし分割表の基準が異なるのでオッズ比も異なります

3群以上の場合

帰無仮説:3群における効果率は等しい
対立仮説:3群における効果率は等しくない(3群のどこかで効果率が異なる)
または
帰無仮説:クラメールの連関係数 = 0
対立仮説:クラメールの連関係数 > 0

要するにカイ二乗検定の独立性の検定です・・・

tab3 <- matrix(c(15, 12, 6, 3, 3, 10), 3, 2)
tab3

以下のような分割表を検定してみます

Z法の効果効果あり効果なし
重度群12315
中等度群212849
軽度群6814
合計291342
fisherの正確検定
fisher.test(tab3)

fisherの正確検定は大変複雑な計算になりますが、Rでは一瞬で算出してくれます。両側検定をした結果は、p=0.048 となり統計学的に効果率には差があるという結果になりました。

fisher.test(tab3, alternative = "less")
fisher.test(tab3, alternative = "greater")

3群以上の場合は、3群の割合の差の二乗の和が0か0じゃないかの検定なので、片側検定の意味はありません。なので上記のように片側検定を行っても同じ答えになります。

カイ二乗検定も同様ですが、3群以上のFisherの正確検定は、3群における効果の有無に差があるかどうかを検定しています。以下、カイ二乗検定の結果です。

カイ二乗検定
chisq.test(tab3)

p値は少し異なりますが、カイ二乗検定もやはり有意として答えてくれません。

fisherの正確検定、カイ二乗検定ともに3群の効果率には違いがあるという結果になりました。計算は非常に面倒なことになるので(2✕2分割表でも大変でした)、他の参考書などをご参照ください。

多重比較

グループ間の割合に差があるということなので、多重比較を行う場合が多いと思われます。でも、個人的には井口先生のご意見に賛成です。以下をご参照ください。

Fisher 正確検定の後に多重比較するな
...

このことを理解した上で、一応多重比較の方法のみ記載しておきます。
使用するパッケージは、RVAideMemoire

インストールしてない場合は

install.packages(“RVAideMemoire”)

library(RVAideMemoire)

以下のようなp値調整の方法があります

ホルム

実際にやってみましょう
A=1行目、B=2行目、C=3行目

#ホルム
fisher.multcomp(tab3, p.method="holm")
ボンフェローニ
#ボンフェローニ
fisher.multcomp(tab3, p.method="bonferroni")

どちらの方法も1行目と3行目に有意差ありという判定でした。

参考ページ ↑↑↑

ダメ出し 間違い、分かりにくい部分などのご意見をお待ちします

タイトルとURLをコピーしました