今回は、ベイズファクターを使った線形回帰分析におけるモデル選択について理論を整理します。 ベイズファクターとはなんぞや？という方はこちらを、Savage-Dickey法についてはこちらを参照のこと。

主な参考文献はこちらです。

線形回帰モデル

まずは線形回帰モデルの導出と、周辺尤度について整理します。

$\boldsymbol{Y} = \left( y_1, \ldots, y_N \right)^T$を応答変数、$\boldsymbol{X}_ 1 = \left( x_ {11} - \bar{x}_1, \ldots, x_{1N}- \bar{x}_1 \right)^T, \ldots, \boldsymbol{X}_ p = \left( x_ {p1} - \bar{x}_p, \ldots, x_{pN} - \bar{x}_p \right)^T$を説明変数とします($\bar{x}_i = \mathbb{E}[\boldsymbol{X}_ i], i=1,\ldots,p$ )。

線形回帰モデルでは、応答変数$\boldsymbol{Y}$が平均ベクトル$\boldsymbol{\mu} = \left( \mu_1, \ldots, \mu_N \right)^T$をパラメータにもつ正規分布に従うと仮定します。

$$ p(\boldsymbol{Y}) = \mathrm{Normal}\left(\boldsymbol{\mu}, \sigma^2 \boldsymbol{I}_N\right) = \cfrac{1}{(2\pi)^{n/2}\sigma^N}\exp\left( -\cfrac{1}{2\sigma^{2}}\left( \boldsymbol{Y} - \boldsymbol{\mu} \right)^T\left(\boldsymbol{Y} - \boldsymbol{\mu}\right) \right) \tag{1} $$

ここで、$\boldsymbol{I}_N$は$N×N$の単位行列です。

平均ベクトル$\boldsymbol{\mu}$は、$\boldsymbol{1}_ N$(1で構成されるN行の列ベクトル)、$\boldsymbol{X}_ 1,\ldots,\boldsymbol{X}_ N$のいずれかもしくはすべての線形結合で説明できるものと仮定します。ここで常に考えなければならない問題は、$\boldsymbol{\mu}$の説明に用いる説明変数の組をどうするかです。これが線形回帰分析においての所謂モデル選択の問題と呼ばれます。

本記事では参考文献と同様、説明変数の選び方を一般化して説明するため、$p$次元ベクトル$\boldsymbol{\gamma} = \left( \gamma_1,\ldots, \gamma_p \right)$をおき、$\gamma_i = 1$のとき、$\boldsymbol{X}_ i$を用いる説明変数に含め、$\gamma_i = 0$のときは逆に含めないこととすることで切片と説明変数によって張られるモデル空間を示すことにし、$\gamma_i = 1$の$\boldsymbol{X}_ i$で構成される$N$×$p_{\boldsymbol{\gamma}}$行列を$\boldsymbol{X}_ {\boldsymbol{\gamma}}$と示します。さらに$\boldsymbol{\gamma}$を用いたモデルを$\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}}$とおき、$\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}}$での係数パラメータを下付き文字で同様に示します。

$$ \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}}　: ~~~ \boldsymbol{\mu} = \boldsymbol{1}_ N \alpha + \boldsymbol{X}_ {\boldsymbol{\gamma}} \boldsymbol{\beta}_{\boldsymbol{\gamma}} \tag{2} $$

ここで、$\alpha$はスカラー、$\boldsymbol{\beta}_　{\boldsymbol{\gamma}}$は$p_　{\boldsymbol{\gamma}}$次元ベクトルです。$\boldsymbol{X}_ {\boldsymbol{\gamma}}$は各変数で中心化されているので、変数のばらつきのみが係数パラメータ$\boldsymbol{\beta}_　{\boldsymbol{\gamma}}$に反映され、変数のスケールはすべて$\alpha$で表現されることになります。

また、どの説明変数も用いないモデルをヌルモデルと呼び、$\mathcal{\boldsymbol{M}}_N$と表現します（下付き文字の$N$はサンプル数ではなく、Null-modelの頭文字）。

$$ \mathcal{\boldsymbol{M}}_ N　: ~~~ \boldsymbol{\mu} = \boldsymbol{1}_ N \alpha \tag{3} $$

線形回帰モデルの評価

ベイズ統計におけるモデル評価方法を整理します。ここで述べることは線形回帰に限定されないより一般的な議論です。

モデル確率の更新

モデル選択やモデルの不確実性の評価にあたっては、不明なパラメータ$\boldsymbol{\theta}_ \boldsymbol{\gamma} = \left( \alpha, \boldsymbol{\beta}_{\boldsymbol{\gamma}}, \sigma^2 \right) ∈ \boldsymbol{\Theta}_\boldsymbol{\gamma}$の事前分布を決め、各モデル$\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}}$の事前確率$p\left(\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}}\right)$を新しく得られたデータに基づきアップデートすることが一つのアプローチとして考えられます。

$$ p(\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}} | \boldsymbol{Y}) = \cfrac{p(\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}})p( \boldsymbol{Y} | \mathcal{\boldsymbol{M}}_\boldsymbol{\gamma} )}{\sum_ {\boldsymbol{\gamma}} p(\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}})p( \boldsymbol{Y} | \mathcal{\boldsymbol{M}}_ \boldsymbol{\gamma} ) } \tag{4} $$

$(4)$式はモデル確率の更新方法を示したもので、ベイズの定理そのままですが、右辺に現れた$p(\boldsymbol{Y} | \mathcal{\boldsymbol{M}}_\boldsymbol{\gamma})$が本記事の鍵となる値、周辺尤度です。$(5)$式に示す通り、周辺尤度は尤度をパラメータの事前分布で重みづけ（積分）したもので、データ$\boldsymbol{Y}$に対するモデルの平均的な説明力を示します。

$$ p(\boldsymbol{Y} |\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}}) = \int_{\boldsymbol{\Theta}_{\boldsymbol{\gamma}}} p(\boldsymbol{Y}|\boldsymbol{\theta}_{{\boldsymbol{\gamma}}},\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}})p(\boldsymbol{\theta}_{{\boldsymbol{\gamma}}} | \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}} ) d\boldsymbol{\theta} _ {\boldsymbol{\gamma}} \tag{5} $$

ベイズファクター

ベイズファクターは$(6)$式に示す通り2つのモデル間の周辺尤度比として定義され、データ$\boldsymbol{Y}$に対する2つのモデルの説明力の比と捉えることができます。

$$ BF\lbrack \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}^{'}} : \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}} \rbrack = \cfrac{p(\boldsymbol{Y} |\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}^{'}})}{p(\boldsymbol{Y} |\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}})} = \cfrac{\int_{\boldsymbol{\Theta}_{\boldsymbol{\gamma}^{'}}} p(\boldsymbol{Y}|\boldsymbol{\theta}_{{\boldsymbol{\gamma}^{'}}},\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}^{'}})p(\boldsymbol{\theta}_{{\boldsymbol{\gamma}^{'}}} | \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}^{'}} ) d\boldsymbol{\theta} _ {\boldsymbol{\gamma}^{'}}}{\int_{\boldsymbol{\Theta}_{\boldsymbol{\gamma}}} p(\boldsymbol{Y}|\boldsymbol{\theta}_{{\boldsymbol{\gamma}}},\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}})p(\boldsymbol{\theta}_{{\boldsymbol{\gamma}}} | \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}} ) d\boldsymbol{\theta} _ {\boldsymbol{\gamma}}} \tag{6} $$

ベイズファクターを使うと、$(4)$式は以下のように変形できます。

$$ p(\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}} | \boldsymbol{Y}) = \cfrac{p(\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}})BF\lbrack \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}} : \mathcal{\boldsymbol{M}}_ {b}\rbrack}{\sum_ {\boldsymbol{\gamma}^{'}} p(\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}^{'}})BF\lbrack \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}^{'}} : \mathcal{\boldsymbol{M}}_b \rbrack} \tag{7} $$

$\mathcal{\boldsymbol{M}}_b$はベイズファクターを算出する際のベースとなるモデルで任意のモデルを選択できますが、通常はどの選択をとってもモデル同士がネストされた関係となるよう、ヌルモデルもしくはフルモデル（すべての説明変数を用いたモデル）が用いられます。

Null-Based Model

$$ BF\lbrack \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}} : \mathcal{\boldsymbol{M}}_ {N} \rbrack = \cfrac{p(\boldsymbol{Y} |\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}})}{p(\boldsymbol{Y} |\mathcal{\boldsymbol{M}}_ {N})} \tag{8} $$

また、任意のモデル間を比較するベイズファクターについても、すべてのモデルに対して包含(Encompassing)関係にあるヌルモデルを基準に計算することができます。

Encompassing approach

$$ BF\lbrack \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}^{'}} : \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}} \rbrack = \cfrac{BF\lbrack \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}^{'}} : \mathcal{\boldsymbol{M}}_ {N} \rbrack}{BF\lbrack \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}} : \mathcal{\boldsymbol{M}}_ {N} \rbrack} \tag{9} $$

事前分布

これまでの議論から、ベイズファクターは事前分布の影響を強く受けることは明白です。客観ベイズの立場から見れば、客観的で研究者間で合意のとれた事前分布の設定方法を確立することで、ベイズファクターやモデル確率の不確定性を克服することが重要です。以下では、望ましい事前分布の設定について整理します。

事前分布に求められる性質

事前分布に求められる性質を列挙します。

・ Location and Scale Invariance

気温や距離など、説明変数の単位や値の大きさに影響しないこと。

・ Consistency

サンプル数が極大に近づくと、ベイズファクターが適切な値に収束すること。つまり$\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}}$と$\mathcal{\boldsymbol{M}}_ N$を比較するとき、「真の」モデルが$\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}}$であるならば、$BF\lbrack \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}} : \mathcal{\boldsymbol{M}}_ {N} \rbrack \rightarrow \infty$に、「真の」モデルが$\mathcal{\boldsymbol{M}}_ N$ならば$BF\lbrack \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}} : \mathcal{\boldsymbol{M}}_ {N} \rbrack \rightarrow 0$にそれぞれ収束すること。

・ Consistent in Information

データが、モデルの比較に用いられる値（線形回帰分析の場合、決定係数）を介してのみベイズファクターに影響すること。線形回帰分析の場合、$R^2 = 1$のとき$BF\lbrack \mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}} : \mathcal{\boldsymbol{M}}_ {N} \rbrack \rightarrow \infty$となること。

・ Computationally Convenience

ベイズファクターが解析的に計算可能であること。ただし近年ではMCMC法やブリッジサンプリング等の推定手法が進歩しているので、必ずしもこの条件にとらわれる必要はないかもしれません。

Zellner-Siow’s Priors

前節で挙げた事前分布の望ましい性質を満たす線形回帰分析におけるパラメータの事前分布として、Zellner-Siow’s Priorsが提案されています。

$(10)$、$(11)$式はZellner(1986)が提案した事前分布で、Zellner’s g-Priorsと呼ばれます。この事前分布は次元のペナルティとしての役割を持つパラメータ$g$について特定の値を与える必要があり、それが前述のConsistencyを満たせない要因となっていました。一方、ZellnerとSiow(1980)が提案したZellner-Siow’s Priorsは、パラメータgについて事前分布を与えることで、事前分布の望ましい性質を満たすことに成功しています。

また、下記$(11)$式に示すように、$(10)$、$(11)$式は$\boldsymbol{\beta}_\boldsymbol{{\gamma}}$の事前分布としてコーシー分布を採用することと等価です。これは$g$について積分することで確認できます（筆者未導出）。

$$ \begin{split} \int_{0}^{∞} \boldsymbol{\beta}_\boldsymbol{{\gamma}} | \sigma^2,g ~~ dg &=& \int_{0}^{∞} \mathrm{Normal}\left(\boldsymbol{0}, g\sigma^2\left( \cfrac{\boldsymbol{X}_ \boldsymbol{{\gamma}}^{T} \boldsymbol{X}_ \boldsymbol{{\gamma}}}{N} \right)^{-1}\right) dg \\ &=& \mathrm{Multivariate-Cauchy}\left( \boldsymbol{\beta} _\boldsymbol{\gamma} | \boldsymbol{0}, N r \sigma^2 \left( \boldsymbol{X}_ \boldsymbol{{\gamma}}^{T} \boldsymbol{X}_ \boldsymbol{{\gamma}} \right) ^{-1} \right) \\
&=& \cfrac{\Gamma\left( \cfrac{1 + p _\boldsymbol{\gamma}}{2} \right)}{\Gamma\left(\cfrac{1}{2}\right) \pi^{\cfrac{p _\boldsymbol{\gamma}}{2}} \left| N r \sigma^2 \left( \boldsymbol{X}_ \boldsymbol{{\gamma}}^{T} \boldsymbol{X}_ \boldsymbol{{\gamma}} \right) ^{-1} \right|^{\cfrac{1}{2}} \left( 1 + \cfrac{1}{N r \sigma^2} \boldsymbol{\beta} _\boldsymbol{\gamma}^{T} \boldsymbol{X}_ \boldsymbol{{\gamma}}^{T} \boldsymbol{X}_ \boldsymbol{{\gamma}} \boldsymbol{\beta} _\boldsymbol{\gamma} \right)^{\cfrac{1 + p _\boldsymbol{\gamma}}{2}} } \end{split} \tag{13} $$

ハイパーパラメータ$r$は、$\boldsymbol{\beta} _\boldsymbol{\gamma}$の事前分布（$(13)$式の多変量コーシー分布）のscale matrixの大きさを調整する役割を持ちます。Richard D. Morey氏は現在のところ$r = \cfrac{\sqrt{2}}{4}$を標準値として推奨しているようです。

$\boldsymbol{\beta} _\boldsymbol{\gamma}$について説明変数や応答変数のスケール等に配慮した事前分布を与えている$(11)$式の解釈は直感的です。まず$g$は各スケールについて標準化した係数パラメータの分散の意味をもちます。次に$\sigma^2$はその値を応答変数の単位にスケーリングします。

$\left( \cfrac{\boldsymbol{X}_ \boldsymbol{{\gamma}}^{T} \boldsymbol{X}_ \boldsymbol{{\gamma}}}{N} \right)^{-1}$ですが、これは$\boldsymbol{X}_ \boldsymbol{{\gamma}}$の分散共分散行列になるので、$g$を$\boldsymbol{X}_ \boldsymbol{{\gamma}}$の単位にスケーリングする役割をもちます。

ベイズファクターの計算

以降ではZellner-Siow’s Priorsをおいたときのベイズファクターの計算方法を整理します。

1変数積分近似による推定

概要を以下に示します。

$(14)$式の$\pi(g)$以外の部分は、$(10)$$(11)$式で表現されるZellner’s g-Priorsにおける周辺尤度を解析的に求めることで算出されたものです。但し、Zellner’s g-Priorsはパラメータ$g$については固定値を想定していたようで、さらに$g$の事前分布$(10)$式が加わったZellner-Siow’s Priorsについては、$g$についてのみ解析的な計算が不可能となってしまっています。そのため、$(14)$式は$g$についての一変数積分が残っていますが、これはガウス求積法などの手法を使って高精度に推定することが可能です。

ちなみに、本手法を使ったベイズファクターの計算はRichard D. Morey氏らが開発しているBayesFactorパッケージのregressionBF()関数で可能です。

Savage-Dickey法による推定

Savage-Dickey法の概要を以下に示します。

$(11)$式より、 $$ p(\boldsymbol{\beta}_ {\boldsymbol{\gamma}} = \boldsymbol{0} |\mathcal{\boldsymbol{M}}_ {\boldsymbol{\gamma}}) = \cfrac{\Gamma\left(\cfrac{p_\boldsymbol{\gamma}+1}{2}\right)}{ \pi^{\cfrac{p_\boldsymbol{\gamma}+1}{2}} \left| N r \sigma^2 \left( \boldsymbol{X}_ \boldsymbol{\gamma}^T \boldsymbol{X}_ \boldsymbol{\gamma}\right) ^{-1} \right|^{-\cfrac{1}{2}}} $$

であることは明らかですから、あとは、$p(\boldsymbol{\beta} _\boldsymbol{\gamma} = \boldsymbol{0} | \boldsymbol{Y}, \mathcal{\boldsymbol{M}} _\boldsymbol{\gamma})$をどうやって求めるかが問題です。

ここでは、以前の記事でも採用したConditional Marginal Density Estimator(CMDE) を使って$p(\boldsymbol{\beta} _{\boldsymbol{\gamma}} = \boldsymbol{0} | \boldsymbol{Y}, \mathcal{\boldsymbol{M}} _{\boldsymbol{\gamma}})$を推定する方法を提案します。

※上記の証明はストレートかつ冗長なので割愛します。自分で導出したい人はこちらが大いに参考になります。

まとめ

本記事では以下の内容について整理しました。

• ベイズファクターを使ったモデル確率の更新

• ベイズファクターを計算する際に求められる事前分布の性質

• 線形回帰モデルにおいて提案された適切な事前分布

実践編はこちらです。ベイズファクターを実際に算出し、モデルをどのように評価できるのか見ていきます。
さらにこちらの記事ではベイズファクターを推定するための手法を複数実践しているので、読んでみてください。