Heck, D. W. (2018) A caveat on the Savage–Dickey density ratio: The case of computing Bayes factors for regression parameters. British J. Mathematical and Statistical Psychology.
仕事の都合で読んだノンパラ検定のベイズ・ファクターについての論文(van Doorn, et al. 2020) で、ネストされたモデルのベイズファクターを点密度だけで簡単に求めるというくだりが出てきて (Savege-Dickey 密度比)、よくわからんかったので探して読んでみた。
論文の主旨は、Savege-Dickey 密度比を使ってはいけない場面があるよ、特に回帰係数の検定のときには気をつけなさい、というもの。
1. イントロダクション
ベイジアン・モデル選択において、ベイズ・ファクター(モデルの周辺尤度の比; 以下BF)は直感的かつ筋の通った指標を提供してくれる。
ベイズ・ファクターを求める方法としてよく使われているのがSavage-Dickey密度比(以下SDDR)である。SDDRが主張するのは、BFは「ある条件」の下で、検定したいパラメータの、ある制約(たとえば\(\theta =0\))のもとでの事前密度と事後密度の比になっている、ということである。つまり積分しなくてもいいわけで、これは助かる。
問題はこの「ある条件」が成り立たない場合である。実は、検定したいパラメータと局外パラメータが独立でない場合には、この「ある条件」は成り立たない。そういう事態は、たとえば、Jeffreys-Zellner-Siow(JZS)事前分布を使った線形重回帰で、いくつかの共変量がネストされたモデルと一般的モデルの両方に含まれているときによく生じる。ベイジアンANOVAも、一般化線形モデルも、潜在過程と行動データを回帰構造でリンクさせる認知過程モデルも、みなこの問題と無縁でいられない。
2. Savage-Dickey密度比
データを\(x\), 尤度関数を\(f(x|\theta, \Psi)\)とする。\(H_0: \theta = \theta_0\), \(H_1: \theta \neq \theta_0\)について検定したい。
BFは、\(H_1\)に対する\(H_0\)の証拠の強さを表す値で、周辺尤度の比として定義される。$$ B_{01} = \frac{\int f(x|\theta = \theta_0, \Psi) \pi_0(\Psi) d\Psi}{\int f(x|\theta, \Psi) \pi_1(\theta, \Psi) d\theta d\Psi} $$ ただし\(\pi_1(\theta, \Psi)\)は\(H_1\)のもとでのすべてのパラメータの同時事前分布, \(\pi_0(\theta)\) は\(H_0\)のもとでの局外パラメータの事前分布である。[分母はメモの誤記ではなくて、二重積分のはずなのにインテグラルはひとつだ。数学からっきしだめなのでわかんないんだけど、こういう略記ってありなんすか]
これは積分があるので計算が難しい。そこでSDDRの登場である。$$ B_{01} = \frac{\pi_1(\theta_0|x)}{\pi_1(\theta_0)} $$ つまり、(包含するほうのモデルのなかにある)検定したいパラメータの、事後の周辺密度と、制約の下での事前の周辺密度の比である。ふつう事前分布はよく知られた分布族に属しているから(t検定とか回帰の場合は単変量コーシー分布を使うことが多い)、分母を求めるのはかんたんである。分子はフルモデルの下で\(\theta\)の周辺事後密度から求めるんだけど、解析的には難しいことが多いので、まずMCMCで事後標本をサンプリングし、カーネル密度推定なり正規近似なりで\(\theta = \theta_0\)の密度を求める。
3. 一般化Savage-Dickey密度比
SDDRの導出においては、$$ \pi_0(\Psi) = \pi_1 (\Psi | \theta = \theta_0)$$ すなわち、\(H_0\)の下での\(\Psi\)の事前分布が、\(H_1\)の下での\(\Psi\)の条件付き事前分布と同一であるという仮定が必要である。従って、等値制約\(\theta = \theta_0\)を置いたときには、事前分布\( \pi_1(\theta, \Psi) \)のフルモデル\(H_1\)は事前分布\(\pi_0(\Psi)\)の帰無仮説\(H_0\)に縮約されなければならない。いいかえると、局外パラメータは\(H_0\)でも\(H_1\)でも全く同じ役割を果たさなければならない。
この仮定は維持されることも多い。ふつう\(\theta\)は効果パラメータで(平均の差とか)、\(\Psi\)は残差分散のような局外パラメータで、後者は\(H_0\)でも\(H_1\)でも同じ役割を持つからだ。一般に、この仮定は$$ \pi_1(\theta, \Psi) = \pi_1(\theta) \pi_1(\Psi)$$のときには自動的に満たされる。
次節では、心理学によくあるシナリオで、かつこの仮定が満たされない場合を示しましょう。
この仮定が満たされない場合にBFを求める方法として、以下の一般化SDDRが提案されている: $$ B_{01} = \frac{\pi_1(\theta_0|x)}{\pi_1(\theta_0)} E\left[ \frac{\pi_0(\Psi)}{\pi_1(\Psi|\theta = \theta_0)} \right]$$ただし\(E[\cdot]\)は\(\pi_1(\Psi|\theta = \theta_0)\)を通じた期待値である。この修正式は元のBFと違って\(H_0\)にも依存するので、2つのモデルからのサンプリングが必要になり、ちょっと遅くなる。
4. 回帰パラメータのベイズ・ファクター
次の重回帰モデルについて考えよう。$$ y_i = \mu + X_i \beta + \epsilon_i$$ \(X_i\)は計画行列\(X\)の\(i\)行目。共変量は\(P\)個ですべて中心化済(\(X\)の列は平均0)としよう。\(\epsilon_i\)は\(N(0,\sigma^2)\)にiidに従うとする。
傾きパラメータのうち一部だけについて、それがゼロでないかどうかに関心があるとしよう。関心がある\(K\)個のパラメータを\(\beta_K\)、残り\(R\)個のパラメータを\(\beta_R\)と書く。仮説を以下とする:
- \(H_P: \beta \neq 0_p\) … すべてのパラメータを自由推定するフルモデル
- \(H_R: \beta_K = 0_K, \beta_R \neq 0_R\) … 関心あるパラメータがゼロだというネストされたモデル
- \(H_0: \beta = 0_P\) … グローバル帰無仮説 Z
ベイジアンモデル選択では、すべてのモデルに含まれているパラメータを共通パラメータと呼ぶことが多い。共通パラメータの事前分布はBFにはあまり影響しないことが知られている。重回帰の場合でいうと、\(\mu, \sigma^2\)が共通パラメータである。事後分布の導出のためには$$ \pi(\mu, \sigma^2, \beta) = \pi(\beta|\mu, \sigma^2) \pi(\mu, \sigma^2)$$ 無情報事前分布としては \(\pi(\mu, \sigma^2) \propto 1/\sigma^2\)が用いられることが多い。ご存じJeffrey事前分布である。
さて、\(\beta_K = 0\)について検定したい場合、\(\beta\)の正則事前分布が必要になる。もし\(H_R\)と\(H_P\)を比べたいだけなら\(\beta_K\)の事前分布さえ正則であればよいが、複数のネストされたモデルを比較したり、\(H_0\)とも比較したい場合には、\(\beta\)のすべてについて正則でなければならない。
こういうときに用いられるのがJZS事前分布 $$ (\beta | \sigma^2) \sim MVC(0_P, \gamma^2 \sigma^2 C^{-1}) $$ である。MVCとは多変量コーシー分布。\(\gamma\)はユーザがアプリオリに決める定数。\(C = X’X / N\)は中心化ずみ計画行列の共分散行列。上のは\(H_P\)の事前分布だけど、\(H_R\)だったら\(R\)次元のJZS事前分布をつかうわけである。
JZS事前分布は中心が0で対称。尺度不変(従属変数や共変量のスケールを変えてもBFが変わらない)。\(\gamma\)で期待sている効果量についての信念を表現できる(BayesFactorパッケージでは中程度の効果量\(\gamma = \sqrt{2}/4\)になっている)。他にもいろいろ美点がある。詳しくはRouder & Morey(2012 Multivariate Behav.Res.), Liang et al.(2008 JASA)をみよ。
さて、ここで\(B_RP\)を求めるとき、検定したいパラメータ\(\beta_K\)と残りのパラメータ\(\beta_R\)は独立でないから、SDDRは使えない。
細かく見てみてみよう。SDDRを使うためには$$ \pi_1(\theta, \Psi) = \pi_1(\theta) \pi_1(\Psi)$$が満たされていればいいわけだ。ネステッド・モデル(左辺)では、$$ (\beta_R | \sigma^2) \sim MVC(0_R, \gamma^2 \sigma^2 C_R^{-1}) $$である。ところがフルモデル(右辺)では、\(\beta_K = 0_K\)と条件づけたときの\(\beta_R\)の条件付き事前分布 \( (\beta_R|H_P, \beta_K = 0_K, \sigma^2) \)は多変量コーシー分布にならない。自由度\(1+K\)の多変量t分布になってしまうのだ。\(K\)が大きくなるほど、多変量コーシー分布に比べてゼロ近辺が盛り上がることになる。
数値例… [面倒なの読み飛ばした]
考察。
[大幅に中略…] このように、JZS事前分布で一部の傾きパラメータについて検定する際にSDDRは使えない。なお、すべての傾きパラメータが0かどうかを検定する際にはSDDRでよい。またg事前分布の場合もSDDRでよい。
云々。
… 勉強にはなったけど、正直なところ偏回帰係数の検定にはいま興味がなくて、ほんとはSavage-Dickey密度比の導出について知りたかったのであった。Wagenmakers, Lodewyckx, Kuriyal, & Grasman (2010 CogPsy)をみるといいらしい。でもCogPsyの論文なんて… いまさら読みたかねえなあ…