Held, L., & Ott, M. (2018) On p-value and bayes factors. Annual Review of Statistics and Its Application, 5, 393-419.
題名の通り、p値とベイズ・ファクターの関係についての解説。仕事の都合で調べたいことがあって読み始めたんだけど、いやあ、難しかった… 面倒くさかった…
1. イントロダクション
[…中略…] 本論文では、効果がないという点帰無仮説の検定で得られるP値をベイズ・ファクター(BF)に変換する方法について扱う。
1.1 ベイズ・ファクター
\( H_0: \theta = \theta_0 \)のとき(典型的には \( H_0: \theta = 0 \)のとき)、以下の対立仮説が考えられる。
- simple alternative. \(H_1: \theta = \theta_1 \neq \theta_0\)
- composite alternative. 単に\(H_1: \theta \neq \theta_0\). このときベイジアン・アプローチでは事前分布\( f(\theta|H_1) \)が必要になる。
- local alternatives. 事前分布として、\(\theta_0\)のまわりで単峰対称な分布を考える。
- non-local alternatives. \(\theta_0\)で密度が0になるような分布を考える。特に\(H_1: \theta = \theta_1 \neq \theta_0\)だったらsimple alternative。
[simpleとcompositeがどう違うのかわからず面食らったのだが、たとえば2群の分布の期待値の差がないというのが\(H_0\)だとして、なんだか知らんが差があるというのがcompositeな\(H_1\), 3ペリカ(なんだその単位は)の差があるはずです！というのがsimpleな\(H_1\)、ということだろうか。えーと、じゃcomposite alternativeで\(H_1\)についての事前信念が区間上の一様分布のときはどう分類するんだろう？ unimodal symmetricではないからlocal alternativeではないよね…]
BFとは事前オッズを事後オッズに変換するもので、証拠の強さを表す。$$ \frac{Pr(H_0|y)}{Pr(H_1|y)} = BF(y) \frac{Pr(H_0)}{Pr(H_1)} $$ ここで $$ BF(y) = \frac{f(y|H_0)}{f(y|H_1)}$$ と書くと、分子は\(H_0\)のもとでの\(y\)の尤度 \(f(y|H_0) = f(y|\theta = \theta_0)\)、分母は\(H_1\)のもとでの\(y\)の尤度 \( f(y|H_1) = \int f(y|\theta) f(\theta|H_1) d\theta\) である。simple alternativeであれば分母は\(f(y|\theta = \theta_1)\)だから、BFはただの尤度比になる。
以下では\(BF \leq 1\)の場合について考える。
上に述べたBFは\(y\)に基づいているのでdata-based BFという。これに対して、P値から求める\( \frac{f(p|H_0)}{f(p|H_1)} \)のことをP-based BFという。[最初にめくった際にはまずここで混乱したんだけど、p値をデータとして捉えるという発想なわけね]
\(p\)が\(H_0\)のもとでの両側P値だとすると、分子\(f(p|H_0)\)はふつう一様である、つまり\( f(p|H_0) = 1 \)であると考えていることになる。なぜならネイマン・ピアソン的仮説検定というのは任意のType I誤差率 \(\alpha (\in (0,1)) \)について\( Pr(p \leq \alpha | H_0) = \alpha \)を維持しようとするわけだから[← 左辺の\(\alpha\)は名目の\(\alpha\)、右辺の\(\alpha\)は実質の\(\alpha\)ということだろうか。ううう… この理屈、反論できないけれど、なんだか騙されているような気持ち悪さがある…]。いっぽう分母の \( f(p|H_1) \) についてはなんらか仮定を置かないといけない。ベータ分布します、とか。
検定統計量から求める\( BF(t) = \frac{f(t|H_0)}{f(t|H_1)} \)をtest-based BFという。P値と検定統計量が一対一であればP-based BFとtest-based BFは等しい。
1.2 最小ベイズ・ファクター
さきほど触れたように、\(f(p|H_1)\)はふつうなんらかのパラメータ\( \eta \)に依存する。そのMLEを\( \hat{\eta} \)とすると、$$ minBF(p) = \frac{f(y|H_0)}{f(y|\hat{\eta}, H_1)} $$はありうるBFのなかで最小のBFとなる。これを最小P-based BFという。[なるほど。分母\(f(y|\hat{\eta}, H_1)\)は尤度\(f(y|\eta, H_1)\)を最大化する\(\eta\)を埋めた最大尤度だからね]
P値は\(H_0\)を支持する証拠ではなくて支持しない証拠を評価することしかできないという非対称性を持っているが、最小BFも同様の非対称性を持っている。
1.3 例
例その1。
条件間に差なしという\(H_0\)について\(\alpha = 0.05, \beta = 0.2\)となるように設計したRCTで、両側P値が0.01になりました。でも研究者はP値が良くないと知ってたので統計家にBFを出せといいました。[← いまいちリアリティが感じられない設定だ、はっはっは]
統計家はこう考えました。標本サイズが十分に大きければ、検定統計量\(t\)は\(H_0\)のもとで\(N(0,1)\)に従う。\(H_1\)のもとで\(t \sim N(\mu, 1)\)と仮定すると、標準正規CDFを\(\Phi(\cdot)\)として $$ \mu = \Phi^{-1}(1-\alpha/2) + \Phi^{-1}(1-\beta) $$だ。計算してみると、\(\mu = 0.280\)となりました。
[話の本筋ではないのだが、ここで統計家はいったいなにをいっておるのかというと…いま考えているのは条件間の差\(\theta\)について\(H_0: \theta = 0, \ H_1: \theta = \theta_1 \neq 0\)という検定である。標本から得た\(\hat{\theta}\)を検定統計量\(t\)に変換したとき、\(H_0: E(t) = 0, \ H_1: E(t) = \mu\)となると仮定し、さらに \(t \sim N(E(t), 1)\)と仮定する。標本設計時には、臨界値\(\pm c\)について、\(H_0\)のもとで運悪く\(c < t\)となる確率、つまり\(N(0,1)\)の\(c\)より右の面積が\(\alpha/2\)となるように設計したのだから\( \Phi(c) = 1-\alpha/2 \)である。また\(H_1\)のもとで運良く\(c < t\)となる確率、つまり\(N(\mu, 1)\)の\(c\)より右の面積が\(1-\beta\)になるように設計したのだから\( \Phi(c-\mu) = \beta\), つまり \( \Phi(\mu-c) = 1-\beta\)である。ここから\( c= \Phi^{-1}(1-\alpha/2), \mu-c = \Phi^{-1}(1-\beta), \mu = \Phi^{-1}(1-\beta) + \Phi^{-1}(1-\alpha/2) \)である。やっと追いついた…]。
統計家はさらにこう考えました。両側P値 \(p = 2(1-\Phi(|t|)) \)は\(t\)の一対一の関数でなくて\(t^{*} = |t|\)の一対一の関数だ。\(t^{*} = \Phi^{-1}(1-p/2)\)としよう。計算すると、\( t^{*} = 2.58\)となりました。
さて、統計家が考えるに、標準正規PDFを\(\phi(\cdot)\)として、$$ BF(t^{*}) = \frac{f(t^{*}|H_0)}{f(t^{*}|H_1)} = \frac{2 \phi(t^{*})}{ \phi(t^{*}+\mu) + \phi(t^{*}-\mu) } $$ であります[←ここの導出はAppendixをみよとのこと]。計算すると、BFは1/13となりました。
いやまてよ、標本サイズ設計が間違っている可能性もあるぞ。真の処理効果は事前の想定と異なっており、検定力は実は\(1 -\beta\)ではないかもしれない。\(t^{*} = 2.58\)はそのままにして、\(BF(t^{*})\)の最小値 $$ minBF(t^{*}) = min_\mu \{BF(t^{*}) \} $$を求めてみよう。するとBFは1/14となりました。
やれやれ… と思っていたら、同僚の統計家がこんなことをいいだしました。もっと簡単にP値ベースの最小BFを求めることができるよ。Sellke et al.(2001)の「-e p log(p)」調整をしらないの？ いわれるとおりに求めてみたら、こんどはBFが1/8になりました。(実はこれはsimple alternativeではなくlocal alternativeの下界(lower bound)なのであります)
別の同僚はこんなことをいいだしました。Goodman(1999b)の「exp(-t^2/2)」はどう？ 求めてみるとBFは1/28、こんどは小さくなってしまいました。(実はこれはたしかにsimple alternativeではあるんだけど、効果の方向についての余計な仮説がはいっているのであります)
統計家はすっかり混乱してしまい、なにをどう報告したらいいのかわからなくなってしまいましたとさ。[← ブラックな落ちだ…]
例その2。
P値を上のように確証的に使うのではなく、探索的に使うこともある。たとえばあるデータで[説明略]、エンドポイントにおける生存を17個の共変量にロジスティック回帰した。予測子の重要性を評価するために、フルモデルにおける17変数のP値を評価し、ついでにそれぞれの最小BFを求めた。こういう場合はlocal alternativeがふさわしい。
2. 歴史的レビュー
2.1 simple alternative
Edwards et al.(1963)はこう述べた。\(y \sim N(\theta, \sigma^2) \)、\(\sigma^2\)は既知とする。\( t= (y-\theta_0)/\sigma \)とする。すべての可能な事前密度 \(f(\theta|H_1) \)に関して、\( BF(y) = f(y|H_0)/f(y|H_1) \)は下界 $$ minBF(y) = \exp(-t^2/2) $$ を持つ。
よくみると、data-based BFなのに、下界は\(t\)の関数である。\(t | H_0 \sim N(0,1), y|H_1 \sim N(\mu, 1) \)の下での\(BF(t)\)の下界\(minBF(t)\)もこれと等しくなる。このように、検定統計量と事前分布を上手く選べば、test-based BFとdata-based BFは等しくなることが多い。
いっぽう、P値を検定統計量にどう変換するかは明確でない。Edwardsらは片側検定に関心を絞っていたので、\( t = \Phi^{-1} (1-p) \)と考えていた(後述するが、このEdwards下界は、local alternativeでg事前分布を仮定したときのBFに対してシャープな下界を提供する)。
いっぽう、Goodman(1999)は両側検定のP値をこの式に放り込んで「可能な限り最小のBF」を得れば良いと考えた。しかしこの考え方には問題がある。この式の\(t\)と両側検定のP値 \(p = 2(1-\Phi(|t|))\)との関係は一対一でない。\(t\)に基づく最小BFは処理効果の方向を考慮にいれているが(\(t\)の符号として)、\(p\)に基づく最小BFはそうでない。Goodmanのアプローチは、Edwards下界を \(p\)に対応する片側P値 \(p/2\)に適用したもので、効果の方向についての情報が含まれていると考えた方が良いだろう。[←頭が混乱してくるけど、要はGoodman下界はおかしいってことね？]
1.3節の例その1でみたいように、\(t^{*} = |t|\)は\(p\)と一対一の関係があり、\(t^{*}\)ベースのBFを\(\mu\)を通じて最小化すれば、simple alternativeのもとでの最小BFが得られる。これは結局、\(f(\theta|H_1)\)が\(\theta_0\)のまわりで対称であることを求めているのと同じである[←???]。両側P値が0.1より小さいとき、次の近似が成り立つ。$$ minBF(t^{*}) \approx 2 \exp(-t^{*2}/2) $$ P値が十分に小さければ、Goodman下界はこの1/2になるわけだ。
2.2 local alternative
こちらもEdwards et al.(1963)にさかのぼる。\(y \sim N(\theta, \sigma^2) \)、\(\sigma^2\)は既知、\(\theta|H_1 \sim N(\theta_0, \tau^2) \)とする。明確な対立仮説がアプリオリにないときに適切な設定である。
このとき、BFを\(\tau^2\)に関して最小化すると、\(\tau^2 = \sigma^2 \max\{t^2-1, 0\}\)となり、$$ minBF(y) = \begin{cases} |t| \exp(-t^2/2) \sqrt{e} & (|t| > 1) \\ 1 & (otherwise) \end{cases} $$ となる(\(e\)はネイピア数)。simple alternativeの場合は\(t\)を使った場合と\(t^{*}\)を使った場合でBFが異なるが、ここでは同じになる。
よくみると、data-based BFなのに、この下界は\(t^{*} = |t|\)の関数である。実際、\(H_0\)のもとで\(t^{*}\)が\(N(0,1)\)を折りたたんだ分布に従い、\(H_1\)の元で\(t^{*}\)が\( N(0,1+\tau^2/\sigma^2) \)を折りたたんだ分布に従うと仮定すると、test-base BFもこうなるということを示せる。simple alternativeのときは\(H_1\)のもとでの\(t^{*}\)の平均は0でないから\(t\)のBFと\(t^{*}\)のBFはちがうのだが、ここでは\(t\)のBFと\(t^{*}\)のBFは同じになる。
なお、事前分布が非正規の場合まで拡張して考えることももできるが、BFはちょっと小さくなるだけである。
2.3 P-based BF
[この小節、ややこしいのでほぼ全訳]
ここまでの最小BFは\(t\)に依存し、よって\(p\)には間接的にしか依存しない。これに対し、Vovk(1993), Sellke et al.(2001) が提案し広く用いられている「-e p log p」calibrationは、\(p\)に直接依存する。$$ minBF(p) = \begin{cases} -e p \log p & (p < 1/e) \\ 1 & (otherwise) \end{cases} $$ この最小BFのかんたんな導出のためには、まず\(H_0\)のもとで正確P値が[0,1]に一様分布するという仮定が必要である。対立仮説のもとでは小さなP値が期待されるので、[\(H_1\)のもとでのP値の分布として]単調減少する密度関数である ベータ事前分布\(Be(\xi, 1), \xi \leq 1\)のクラスを考える。するとこの最小BFが導出される。Sellke et al.(2001), Held & Ott(2016)を参照のこと。なお、導出においてはMLE \(\hat{\xi}_{ML} = min\{-1/\log(p), 1\}\)が用いられている。
Sellke et al.(2001)は別の導出も示している。そこでは、\(H_1\)のもとでP値がベータ分布のクラスに従うという仮定が不要である。
Held(2010)によれば、上の最小BFは、\(\theta\)が2次元で標本サイズが大きいときの、g事前分布のもとでのtest-based BFとしても導出できる。後述する。
この最小BFは、前節で述べたlocal normal alternativeのもとでの最小BFよりも常に小さく、すべてのlocal alternativesというより一般的なクラスにおける下界にほぼ等しい。
ベータ分布\(Be(\xi, 1)\)の事前標本サイズは \(\xi + 1 \leq 2\)であり、常にほぼ無情報である。従って、\(f(p|H_1, \hat{\xi}_{ML})\)は相対的にフラットであり、最小BF\(minBF(p) = 1/f(p|H_1, \hat{\xi}_{ML})\)は相対的に大きい。しかし、密度関数が単調減少するベータ事前分布は\(Be(\xi, 1)\)以外にもある。別の事前分布として\(Be(1, \kappa), \kappa \geq 1\)というクラスも考えられる(これは我々の知る限りこれまで検討されてこなかったクラスである)。このクラスからのベータ分布は、事前標本サイズが\(1 + \kappa \geq 2\)であり、対立仮説のもとでの尤度は\(Be(\xi, 1)\)を事前分布としたときよりも大きな値になりうる。計算してみると、ここでの\(\kappa\)のMLEは\(\hat{\kappa}_{ML} = \max \{ -1/\log(1-p), 1\}\)となり、最小BFは$$ minBF(p) = \begin{cases} -e (1-p) \log (1-p) & (p < 1-1/e) \\ 1 & (otherwise) \end{cases} $$ となる。いわば「-e q log q」calibrationである。\(p\)が十分に小さければ、\(\log(1-p) \approx -ep\) という近似に基づき、\(minBF(p) \approx e p\)と近似できる。
この最小BFは他のどの下界よりもずっと小さい。P値が0.1未満の場合はGoodman下界よりもなお小さい。これは、\(p\)が小さいのに事前標本サイズが大きいからである。いっぽう「-e p log p」calibrationでは事前標本サイズは2を超えない。しかし後述するように、この最小DFは、標本サイズが非常に小さくても、任意の次元\(d\)におけるg事前分布に基づくBFにシャープな下界を与える[←よくわからない…]。しかし、reasnableな大きさの標本サイズでは、「-e q log q」calibrationはあまりに保守的である。
2.4 test-based BF
data-based BFの欠点は、\(H_0, H_1\)の未知パラメータに割り振られた事前分布にクリティカルに依存してしまうという点である。さらに、計算には多次元の積分の評価が必要になる。
そこでJohnson(2005)はtest-based BFを推した。いわく、\(H_0\)のもとでの検定統計量の標本分布を考え、その中心だけがずれた分布を、\(H_1\)のもとでの検定統計量の標本分布だと思えばいい。このアプローチだと事前分布はいらないし、BFを閉形式で書ける。[…後略]
2.5 比較
P値を動かしながら、{simple alternative, Edwards(片側), Goodman(両側), local normal alternative, -e p log p, -e q log q} による最小BFがどうなるかを比べてみよう。
[話がややこしすぎてわからなくなってしまったので、整理しておこう。いまやりたいことは、検定統計量ないしP値から、BFの下界を求めることである。出場選手は以下の6名である。
- simple alternative: \( minBF(t^{*}) = \min_\mu \left( \frac{2 \phi(t^{*})}{ \phi(t^{*}+\mu) + \phi(t^{*}-\mu) } \right) \).
- Edward下界: 片側検定のP値から \(t = \Phi^{-1}(1-p)\)を逆算して\(minBF(t) = \exp(-t^2/2)\).
- Goodman下界: 両側検定のP値から \(t = \Phi^{-1}(1-p)\)を逆算して\(minBF(t) = \exp(-t^2/2)\).
- local normal alternative: \( minBF(t) = \begin{cases} |t| \exp(-t^2/2) \sqrt{e} & (|t| > 1) \\ 1 & (otherwise) \end{cases} \).
- -e p log p: \( minBF(p) = \begin{cases} -e p \log p & (p < 1/e) \\ 1 & (otherwise) \end{cases} \).
- -e q log q: \( minBF(p) = \begin{cases} -e (1-p) \log (1-p) & (p < 1-1/e) \\ 1 & (otherwise) \end{cases} \).
嗚呼… ややこしい….]
まず、すべての最小BFがP値よりも大きい。simple alternative の最小BFは、Goodman下界の2倍となる。Edwards下界はsimple alternativeに近い。
-e p log p はlocal normal alternativeに近い。-e q log q はもっとも小さい。
… ここまでで本文22頁のうち11頁。そしてもはや時間切れである。涙を飲んでここで終了としたい。(うそだ、別に涙を飲んでない。正直、想像してた内容より100倍くらい難しいのでうんざりしたのである)
本文の最後に載っていた本論文の要点をメモしておく。
- P値は点帰無仮説\(H_0\)に反する証拠の間接的指標である。BFは\(H_0\)に反する証拠の直接的な量的要約を提供する。
- P値は最小BFに変換できる。最小BFは、ある対立仮説のクラスのなかで、P値が持っている\(H_0\)に反する最大の証拠を定量化している。
- P値が持つ最大の証拠は、P値をどうやって計算したかに依存する。それ[minBFのことね]は一般に標本サイズ増大とともに減少し、関心あるパラメータの次元の増大とともに増大する。P値を最小BFに変換する際にはこの特性を考慮する必要がある。
- P値が持つ最大の証拠は研究デザインにも依存する。対立仮説がきちんと定義されたsimple alternativeであるような検証的研究から得られたP値なのか、仮説生成が目的であってlocal alternativeがよりふさわしいような探索的分析から得られたP値なのか、が問題である。
- 広く用いられている「-e p log(p)」calibrationは、大標本において、かつパラメータがスカラーのとき、local alternativeに対するBFの下界を表している。小標本の場合や関心あるパラメータの次元数が大きいときには、必ずしもそうならない。
引き続き精進します、ってことで… (うそだ、精進する気などさらさらない。むしろ、もういい！！わかった！！P値からベイズファクターを出すのはやめろ、元データから出せ！！ とぶちぎれている)