僕が作った架空のデータで練習しましょう。結果についてはご意見があると思いますが、統計学の練習のためのサンプルですのでご了承ください。
例題）新たなトレーニング方法（Z法）が開発されました。果たしてこのトレーニングは運動障害が重度群と軽度群で効果に差があるのか？
データの準備と要約
Googleでのヒット数(2022.5.19調べ)
| Fisherの正確検定 | 68,000 |
| Fisherの正確確率検定 | 7,060 |
| Fisherの直接確率検定 | 4,770 |
| Fisherの直接法 | 4,190 |
Fisherの正確検定、Fisherの正確確率検定、Fisherの直接確率検定、Fisherの直接法などと呼ばれています（統一してくれれば良いのにといつも思います）。Fisherの正確検定が多かったので、このタイトルに使ってます。
定義
\(\displaystyle オッズ=(\frac{効果あり}{効果なし})\)
- ここではオッズを効果率と定義します
- このようにオッズやオッズ比の定義を確認しておくことが重要です
サンプル（架空のサンプルです）
| 効果あり | 効果なし | 計 | |
| 重度群 | 12 | 3 | 15 |
| 軽度群 | 6 | 8 | 14 |
| 合計 | 18 | 11 | 29 |
重度群の効果率は4、軽度群の効果率は0.75となります
この効果率の比が「1」か「1ではない」かということを調べるために、下記のオッズ比を利用します
\(\displaystyle (オッズ比 =)(\dfrac{重度群の効果率}{軽度群の効果率})(= 5.33)\)
5.33なのでオッズ比は1より大きいようですが・・・
オッズ比が高くなると重度群の効果率が高くなることを意味しており、オッズ比が１より大きくなると軽度群の効果率より重度群の効果率が高いことを意味します
この検定の仮説は以下のようになります
帰無仮説H0: オッズ比＝1（重度群の効果率 ＝ 軽度群の効果率）
対立仮説H1: オッズ比≠1（重度群の効果率 ≠ ）軽度群の効果率
Rで分割表の作成
tabという単語の中にサンプルの分割表を挿入してみます
tab <- matrix(c(12, 6, 3, 8), 2, 2)
tab
tabと入力してEnterをクリックして以下のようになれば成功です
これで「tab」という単語（なんでも構いません）の中に分割表が入りました
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つの違いについて
ｘ：重度群かつ効果あり、効果あり=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\)
＊この確率は超幾何分布を使って求めます (後で出てきます)
ｐ値は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以上になるのは稀であることを意味します
つまり、効果率は１以上と考え重症群の効果率が有意に高いと判断します
fisher.test(tab, alternative = “two.sided”)
P=0.06043 なので有意差はありません
このｐ値は後半で説明します
fisher.test(tab, alternative = “greater”) の結果の説明
- p-value（ｐ値）
0.04601384 - alternative hypothesis（対立仮説）
Rは対立仮説を出力します
真のオッズ比は１ではない
すなわち・・・p値<0.05なので「重症群と軽度群には有意差があり、
AM実施群の効果率は有意に高いという結果になります
注意）両側検定の場合Ｐ値が大きくなり有意ではなくなります - 95 percent confidence interval
下記のオッズ比の95%信頼区間 - 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) $$
なかなか大変です・・・
なのでｎが大きくなると昔は計算が大変だった思います・・・
でも今はRですぐに計算できるので感謝です
fisher.test(tab, alternative = “greater”) のｐ値
# 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”) のｐ値
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])
ｐ値と信頼区間の整合性
fisher.test では、p 値は tsmethod=“minlike” 、信頼区間はtsmethod=“central” に基づいて求められており、整合性に問題が生じます
上記の例teb3をfisher.testで検定した場合
fisher.test(tab3, alternative = "two.sided")
p<0.05なので信頼区間が1を跨いでしまっています
ですので上述したtsmethod=”minlike”で実行することでp値と信頼区間の整合性が担保されます
間違った解釈をしないように！
再掲
| 効果あり | 効果なし | 計 | |
| 重度群 | 12 | 3 | 15 |
| 軽度群 | 6 | 8 | 14 |
| 合計 | 18 | 11 | 29 |
AM実施群が有効であることを証明したい場合
帰無仮説H0: AM群の効果率 ＝ PM群の効果率
対立仮説H1: AM群の効果率 ＞ 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“
)
ということになります
もちろんどちらも同じｐ値になります
＊ただし分割表の基準が異なるのでオッズ比も異なります
3群以上の場合
帰無仮説：3群における効果率は等しい
対立仮説：3群における効果率は等しくない（3群のどこかで効果率が異なる）
または
帰無仮説：クラメールの連関係数 = 0
対立仮説：クラメールの連関係数 > 0
要するにｶｲ二乗検定の独立性の検定です・・・
tab3 <- matrix(c(15, 12, 6, 3, 3, 10), 3, 2)
tab3
以下のような分割表を検定してみます
| Z法の効果 | 効果あり | 効果なし | 計 |
| 重度群 | 12 | 3 | 15 |
| 中等度群 | 21 | 28 | 49 |
| 軽度群 | 6 | 8 | 14 |
| 合計 | 29 | 13 | 42 |
fisherの正確検定
fisher.test(tab3)
fisherの正確検定は大変複雑な計算になりますが、Rでは一瞬で算出してくれます。両側検定をした結果は、p=0.048 となり統計学的に効果率には差があるという結果になりました。
fisher.test(tab3, alternative = "less")
fisher.test(tab3, alternative = "greater")
３群以上の場合は、３群の割合の差の二乗の和が0か0じゃないかの検定なので、片側検定の意味はありません。なので上記のように片側検定を行っても同じ答えになります。
ｶｲ二乗検定も同様ですが、３群以上のFisherの正確検定は、3群における効果の有無に差があるかどうかを検定しています。以下、ｶｲ二乗検定の結果です。
ｶｲ二乗検定
chisq.test(tab3)
p値は少し異なりますが、ｶｲ二乗検定もやはり有意として答えてくれません。
fisherの正確検定、ｶｲ二乗検定ともに３群の効果率には違いがあるという結果になりました。計算は非常に面倒なことになるので（2✕2分割表でも大変でした）、他の参考書などをご参照ください。
多重比較
グループ間の割合に差があるということなので、多重比較を行う場合が多いと思われます。でも、個人的には井口先生のご意見に賛成です。以下をご参照ください。
このことを理解した上で、一応多重比較の方法のみ記載しておきます。
使用するパッケージは、RVAideMemoire
インストールしてない場合は
install.packages(“RVAideMemoire”)
library(RVAideMemoire)
以下のようなｐ値調整の方法があります
ホルム
実際にやってみましょう
A=１行目、B＝2行目、C＝３行目
#ホルム
fisher.multcomp(tab3, p.method="holm")
ボンフェローニ
#ボンフェローニ
fisher.multcomp(tab3, p.method="bonferroni")
どちらの方法も１行目と3行目に有意差ありという判定でした。
ダメ出し 間違い、分かりにくい部分などのご意見をお待ちします