仮説検定の手法ごとに結論がどれほど変わるか比べてみた

はじめに

本記事では、次の 3 つの統計指標を使った検定を同一の問題設定に対し適用し、結論がどのように変わるのかを比べてみました。

• p 値
• ベイズ信用区間
• ベイズファクター

本記事は各手法の詳細を説明するというよりは、手法間で結論がどう変わるかをまとめることに重きを置いています。したがって、p 値ベースの仮説検定やベイズ推論についての知識はある程度仮定します。

ベイズファクターについてはこちらの資料が個人的にはわかりやすかったです。特徴や注意点も書かれているのでおすすめです。

問題設定

コインが公正なものかを確かめたい状況を考えます。公正であるというのは、表が出る確率と裏が出る確率がどちらも $1/2$ であることを言います。

頻度論的に言うと 1 群の母比率の検定ですね。このコインを $N$ 回投げて $n$ 回だけ表が出たというデータが得られたとします。また、表の出る確率を $q$ で表すことにします。

p 値による検定

まずは最もベーシックであろう p 値を使った方法からです。なにはともあれ帰無仮説 $H_0$ と対立仮説 $H_1$ を設定しましょう。

• $H_0$: $q = 1/2$
• $H_1$: $q \neq 1/2$

有意水準 5% の両側検定をしていくことにします。表の出た回数 $X$ は二項分布 $\mathrm{Binom}(X \mid N, q)$ に従うとみなせます。

p 値とは、帰無仮説のもとで手元のデータよりも極端な値が得られる確率のことでした。したがって $n \leq N/2$ のとき、p 値は次のように計算できます。$n \geq N/2$ のときは不等号の向きが逆になります。

$$2 \cdot p(X \leq n \mid H_0)$$

では具体的な値に対して計算していきましょう。コインを投げ上げた回数 $N$ は $100$ とします。表の出た回数 $n$ の何通りかについて p 値を計算します。

手計算する場合にはラプラスの定理を使って二項分布を正規分布に近似しますが、今回はコンピューターに計算してもらうので二項分布のままで良いです。以前イェーツ補正の必要性について記事を書いたこともありますが、ここではこの辺の細かい話は考えないことにしましょう。

Julia のコード例はこのようになります。

using Distributions

N = 100
ns = 38:41
d = Binomial(N, 0.5)

for n in ns
    @info n
    @info round(cdf(d, n) * 2, digits = 3)
end

出力を表にまとめるとこのようになります。表の出た回数が 39 回以下なら帰無仮説は棄却されますが、40 回では棄却されません。

$n$	p 値
38	0.021
39	0.035
40	0.057
41	0.089

ベイズ信用区間による検定

$q$ の事前分布をベータ分布 $\mathrm{Beta}(q \mid a, b)$、尤度関数を二項分布 $\mathrm{Binom}(X \mid N, q)$ でモデリングします。

計算の詳細は省略しますが、ベータ分布は二項分布の共役事前分布なので、事後分布もベータ分布になります。第一パラメータに表の出た回数を、第二パラメータに裏の出た回数を足せばいいので、事後分布は $\mathrm{Beta}(q \mid a + n, b + N - n)$ となります。

信用区間には Highest Density Interval のようなものもありますが、ここでは簡単のためにパーセント点によるものを採用します。つまり、ベータ分布の下側 $\alpha$ パーセント点を $\beta_\alpha$ と表すことにすれば、$[\beta_{2.5}, \beta_{97.5}]$ が 95% 信用区間です。この区間が $1/2$ を含むかどうかで、コインが公正かどうかを結論づけることができます。

では具体例を考えていきましょう。事前分布としては $\mathrm{Beta}(1, 1)$ を考えます。コインを投げる回数は先ほどと同様に $100$ とします。コードとその出力は以下のようになります。信用区間による検定では 40 回でも帰無仮説が棄却されるようですね。

using Distributions

N = 100
ns = 39:42

for n in ns
    d = Beta(n + 1, N - n + 1)
    @info n
    @info [round(quantile(d, 0.025), digits = 3), round(quantile(d, 0.975), digits = 3)]
end

$n$	信用区間
39	$[0.300, 0.488]$
40	$[0.309, 0.498]$
41	$[0.319, 0.508]$
42	$[0.328, 0.518]$

ベイズファクターによる検定

ベイズファクターは検定というよりはモデル選択のための指標です。ここでは次の2 つの事前分布のモデリングを比較することにします。

• $M_0$: $q \sim \mathrm{Beta}(1, 1) =$ 一様分布
• $M_1$: $q = 1/2$

ベイズファクター $BF_{01}$ は次のように定義されます。

$$BF_{01} = \frac{p(X \mid M_0)}{p(X \mid M_1)}$$

まずは分子から計算していきます。$\mathrm{Beta}(1, 1)$ が一様分布であることに注意しましょう。また、$B$ はベータ関数です。

$$\begin{aligned} p(X \mid M_0) &= \int_0^1 p(X, q \mid M_0) dq \\[.5em] &= \int_0^1 p(X \mid q, M_0) \ p(q \mid M_0) dq \\[.5em] &= \int_0^1 \mathrm{Binom}(X \mid N, q) \ \mathrm{Beta}(q \mid 1, 1) dq \\ &= \int_0^1 \binom{N}{n} q^n (1-q)^{N-n} dq \\[.5em] &= \binom{N}{n} B(n+1, N-n+1) \\[.5em] &= \binom{N}{n} \frac{n!(N-n)!}{(N+1)!} \end{aligned}$$

続いて分母です。

$$\begin{aligned} p(X \mid M_1) &= p(X, q = 1/2 \mid M_1) \\[.5em] &= \mathrm{Binom}(X \mid q = 1/2, M_1) \\[.5em] &= \binom{N}{n} \left(\frac{1}{2}\right)^n \left(\frac{1}{2}\right)^{N-n} \end{aligned}$$

分母と分子をひっくり返したほうが計算しやすいので、$BF_{01}$ ではなく $BF_{10}$ を考えていきます。

$$\begin{aligned} BF_{10} &= \frac{(N+1)!}{n!(N-n)!} \left(\frac{1}{2}\right)^n \left(\frac{1}{2}\right)^{N-n} \\[.5em] &= (N+1) \frac{N!}{n!(N-n)!} \left(\frac{1}{2}\right)^n \left(\frac{1}{2}\right)^{N-n} \\[.5em] &= (N+1) \cdot \mathrm{Binom}(X = n \mid N, q = 1/2) \end{aligned}$$

ベイズファクターは確率の比であるため、0 以上の任意の実数を取ります。1 より大きければ、分子のモデルを支持する根拠が強いことを意味し、1 より小さければ、分母のモデルを支持する根拠が強いことを意味します。

いくつ以上なら $M_0$ を採用する、のような明確な基準はありませんが、参考となる基準は Jeffreys さんが提唱しています。以下の表はこちらのページから引用させていただきました。

$BF_{10}$	Interpretation
$> 100$	Extreme evidence for $M_1$
$30 \sim 100$	Very strong evidence for $M_1$
$10 \sim 30$	Strong evidence for $M_1$
$3 \sim 10$	Moderate evidence for $M_1$
$1 \sim 3$	Anecdotal evidence for $M_1$
$1$	No evidence
$1 \sim 1/3$	Anecdotal evidence for $M_0$
$1/3 \sim 1/10$	Moderate evidence for $M_0$
$1/10 \sim 1/30$	Strong evidence for $M_0$
$1/30 \sim 1/100$	Very strong evidence for $M_0$
$< 1/100$	Extreme evidence for $M_0$

では具体例を計算していきます。以下のようなコードで $BF_{10}$ が計算できます。

using Distributions

N = 100
ns = 30:40
d = Binomial(N, 0.5)

for n in ns
    @info n
    @info (N+1) * pdf(d, n)
end

表の出た回数が 40 回くらいだと、どちらのモデルも大差ないことになります。35 回くらいになると $M_0$ に対する強い根拠が得られます。

$n$	$BF_{10}$
30	0.00234
31	0.00528
32	0.01138
33	0.02348
34	0.04627
35	0.08725
36	0.15753
37	0.27249
38	0.45176
39	0.71818
40	1.09523

まとめ

アメリカ統計協会 (ASA) が発表した声明の通り、p 値は正しく取り扱うのが難しい概念です。そのため、検定するならベイズ的手法を使うほうがいいのかなと漠然と考えていました。しかし手法間でどれくらい結論が変わるのかというのが気になり、比較記事を書くに至りました。

今回の設定ではベイズファクターによる結論が多少保守的になりました。ベイズファクターは p 値のいくつかの欠点 (再現性の問題など) を解消する指標として注目されてるんじゃないかなと思います。ただ、明確な基準が無いのが難しい点かもしれません (まあ信頼区間や信用区間の 95% という基準も単なる慣例であって明確な理由があるわけではないですし、なんでも有意かそうでないかという二値に持ち込むことに問題があったりするわけですが)。

しかし使わないと使いこなせるようにはならないので、今後 A/B テストなどする機会があればベイズファクターも計算してみようと思います。