読者です 読者をやめる 読者になる 読者になる

Rを通じて統計学を学ぶ備忘録ブログ

SPSSからRに移行したい私のような人向けのR解説ブログ兼学習用備忘録。

ベイズ統計学その15〜線形回帰モデルとMCMC〜

今回はMCMCを使った線形回帰モデルをRで実践しながら見ていきます。その過程ででくわす基本的な用語も攫っていきます。

ベイズの線形回帰モデルを考える際には、未知パラメタβと誤差項のσ^2を求めることが目標になります。その際に、尤度関数をどう考えるか(定義するか)が大切です。尤度関数は誤差項がどのように振る舞うのか(そのような確率分布に従うのか)によって、定義できますので、そのようなことを念頭に尤度関数を考えます。

前回までに、平均、分散(精度)がわからいない状態な場合、分散が逆ガンマ分布に従い、尤度関数の正規分布を掛け合わせることで、事後分布を得られることがわかりました。

言い換えると、分散σ^2が逆ガンマに従い、分散σ^2が与えられた条件のもと、βが正規分布に従い、分散σ^2とパラメタβが与えられた条件のもと、目的変数yが得られるという感じです。そして、その事後分布を使ってシュミレーションします。 

ベイズ線形回帰モデル~基本用語~

MCMCpackでは色々な分析ができますが、今回は基本的なMCMCregressを使ってベイズ線形回帰モデルを見ていきます。MCMCregress()はMarkov Chain Monte Carlo for Gaussian Linear Regressionであり、MCMCを使った基本的な回帰モデルです。

引数は以下のとおりです。 

MCMCregress(formula, data = NULL, burnin = xxxx, mcmc = xxxx, thin = xxxxx, verbose = xxxxx, b0 = xxxx, B0 = xxxx, c0 = xxxx, d0 = xxxx)

初期値の設定はどうすればよいのでしょうか。事後分布を得るための関係式さえ、適切に定義しておけば、たくさんシュミレーションするなら、初期値は適当で問題ありません。

従い、正規分布の平均(b0 = 0)、分散(B0 = 0)、逆ガンマ分布のシェイプパラメタ、スケールパラメタ( c0 = 0.001, d0 = 0.001)のように設定します。事後分布のパラメタについては正規分布に、分散(精度)については逆ガンマ分布になりますが、当然、事後分布のパラメタの平均、分散(β、σ^2)は事前情報や観測データに依存します。

事前情報に一様分布が当てられている場合や、事前情報を与えなければ(初期値を非常に小さくする)、ベイズ推定と古典統計の回帰分析とほぼ同じ値を得られます。無情報事前分布は積分して1にならないので、非正則事前分布と呼ばれます。古典統計学の線形回帰分析(最小二乗法解)の普遍推定量は、ベイズ推定において、無情報事前分布を仮定した場合とほぼ同じと考えることができます。 

□トレースと自己相関

トレース図はMCMCのサンプルを時系列にプロットしたものです。バーンイン期間やMCMCの質の良さを確認したり、MCMCは時系列に沿ってサンプルを抽出していくので、自己相関の有無も確認します。基本的にマルコフ連鎖は1つ前の状態にのみ依存して、次のステップが決まるので、自己相関がサンプル抽出の数(時間)と共にさがっていくほうが質が良いMCMCと言えます。以下、計算式です。時刻tにおける値をX_tと考え、X_tの平均μ、分散σ^2としたとき、以下のように表せます。

R(t,s) = \frac{E[(X_t - \mu)(X_s - \mu)]}{\sigma^2}\, ,

「平均μに対してある時点の変数xを引いたもの」と「平均μに対してある時点からk時点足した変数xを引いたもの」を分散で割った数です。範囲は[-1~1]となり、1は完全相関、-1は完全反相関となります。自己相関が0に近ずくことが望ましく、以下は良い自己相関の例です。

f:id:teruaki-sugiura:20160206223320p:plain

□バーンイン

MCMCでシュミレーションすると、最初の部分は初期値に依存してしまいます。そのため、この部分を利用すると、事後分布に影響が出るため、不必要と判断して捨てる部分のことです。厳密な表現ではありませんが、分布が定常状態になったあたりから、サンプルとして利用します。以下の例では、-3から0に上がっている部分は初期値の影響を受けているため、捨ててしまいます。 

f:id:teruaki-sugiura:20160206223336p:plain

□シニングインターバル

シニングインターバルはMCMCのサンプルの独立性を担保するためにサンプルを間引いて抽出することで独立性(自己相関を下げる)を維持するために行います。MCMCのサンプルは互いに独立していることが前提条件なので、抽出した部分のみを用いれば、独立しているだろうと考えることができます。

□事後平均・確信区間

事後平均はその名の通り、パラメタの事後分布の平均値のことです。事後分布の確信区間とは、パラメタが振る舞う幅のことで、パラメタが確率的に収まる幅のことです。推測統計学の信頼区間とは異なる意味合いです。

つまり、ベイズ統計では、95%確信区間というのは、パラメータθが確信区間の中に含まれる確率は95%である、と解釈できる一方、推測統計学の95%信頼区間とは、パラメタθを100回推定し、100回の信頼区間を考えると、そのうち95回は信頼区間の中にパラメタθが含まれる、と解釈します。 以下の図の場合、事後平均は0.003593となっており、95%確信区間は−0.41156~0.41213となります。

f:id:teruaki-sugiura:20160206223521p:plain

f:id:teruaki-sugiura:20160206223541p:plain

MCMCの収束判定方法(詳しくは後日)

・Gelman & Rubinの収束判定方法

Potential scale reduction factor (PSRF)値を計算し、PSRF値が「1.05以下(または1近辺)」の場合、収束したと考えます。言い換えると、十分に長いシュミレーション期間だと判断することになります。 

・Gewekeの診断方法

バーンインを除いた2つのMCMCのサンプリング期間(最初の10%、最後の50%)について、Z 統計量を計算し、平均値の差を検定します。例えば、1100回サンプル抽出した場合、バーンイン期間の100回を除き、101~200回のサンプルと501~1000回のサンプルの平均値を比較。検定で「分布が等しいという帰無仮説を受容」した場合、分布が同じと考えて収束していると判断。トレース図で考えるとわかりやすく、パラメタがあっちに行ったり、こっちに行ったりと上下している感じだと収束していないと考えれそうです。

・Raftery & Lewisの収束判定方法

Dependence factorが「5」より大きい場合、MCMCが初期値に依存しているため、収束していないと考えます。そのため、5以下であればMCMCは収束していると考えます。

□使用データセット(swiss)

データセットの説明は、Rドキュメントより。このデータセットは、1888年頃のスイスの47のフランス語圏地域の標準化された出生率と社会経済指標で構成されています。6つの変量に対する47の観測値からなるデータフレームで、各々はパーセント値である、 つまり [0,100] 中の数値。

Fertility:標準化された通常出生率

Agriculture:農業従事男性

Examination:軍隊の検査で最高点を得た被徴兵者

Education:被徴兵者の小学校以後の学歴

Catholic:カソリック (プロテスタントと対比した)

Infant.Mortality:一歳未満死亡幼児

「Fertility」を除く全ての変数は人口に対する割合を与える。

f:id:teruaki-sugiura:20160206223724p:plain

ベイズ線形回帰モデル~R実践~

まずは基本的な回帰分析から。

> lm001 <- lm(Fertility ~ Agriculture + Examination + Education, data = swiss)
> summary(lm001)

f:id:teruaki-sugiura:20160206224057p:plain

次にギブスサンプリングでの回帰分析。今回は無情報事前分布を仮定しますので、回帰係数の数値はほぼ同じになります。

> library(MCMCpack)
> lm002 <- MCMCregress(Fertility ~ Agriculture + Examination + Education, data=swiss, burnin = 1000, mcmc = 10000, thin = 2, b0 = 0, B0 = 0, c0 = 0.001, do = 0.001)
> summary(lm002)
> plot(lm002)

f:id:teruaki-sugiura:20160206224123p:plain

f:id:teruaki-sugiura:20160207000227p:plain

f:id:teruaki-sugiura:20160207000234p:plain

通常の回帰分析の回帰係数とほぼ一緒になりました。係数が有意かどうかの判定のようなものは、パラメタの確信区間を見て判断します。基本的には「0」を跨いでいないものが「有意」というような扱いになります。あくまで、検定をしているわけではないので、有意とは言いません。収束診断と自己相関を見ていきましょう。

> library(coda)
> autocorr.plot(lm002)

f:id:teruaki-sugiura:20160206225803p:plain

自己相関については、Lagが進むごとに0に近づいていることがわかります。

> geweke.diag(lm002)

> raftery.diag(lm002)

f:id:teruaki-sugiura:20160206230605p:plain

geweke.diagは、|Z| =<1.96 のとき定常状態と判断します。または、|Z|の値が1よりも十分に小さいとき、収束と判断します。

raftery.diag(lm002)

f:id:teruaki-sugiura:20160206231936p:plain

Dependence Factorが「5」以下であれば収束と判断されるので、収束していると判断できます。

> lm003 <- MCMCregress(Fertility ~ Agriculture + Examination + Education, data=swiss, burnin = 1000, mcmc = 10000, thin = 2, b0 = 0.001, B0 = 0.001, c0 = 0.0001, do = 0.0001)
> lm002.003 <- mcmc.list(lm002, lm003)
> gelman.diag(lm002.003)

f:id:teruaki-sugiura:20160206225734p:plain

gelman.diagをみていきます。Point Setが1.05以下の場合、収束と判断するので、Educationとsigm2のみ収束していると判断できます。

最後にベイズファクターを用いてモデルの診断を行います。比較のために、変数を1つづつ増やしたモデルを作成します。

> model1 <- MCMCregress(Fertility ~ Agriculture, data=swiss, burnin = 10000, mcmc = 100000,b0 = 0, B0 = 0.001, c0 = 0.001, d0 = 0.001, marginal.likelihood="Chib95")
> model2 <- MCMCregress(Fertility ~ Agriculture + Examination, data=swiss, burnin = 10000, mcmc = 100000,b0 = 0, B0 = 0.001, c0 = 0.001, d0 = 0.001, marginal.likelihood="Chib95")
> model3 <- MCMCregress(Fertility ~ Agriculture + Examination + Education, data=swiss, burnin = 10000, mcmc = 100000,b0 = 0, B0 = 0.001, c0 = 0.001, d0 = 0.001, marginal.likelihood="Chib95")
> BayesFactor(model1, model2, model3)

f:id:teruaki-sugiura:20160206235834p:plain

対数ベイズファクターのところを見ると(対数でなくてもよいが。。。)Model2が最もよいモデルと判断できそうですが、Model3とはどっこいどっこいというかんじですね。

以上で今回はおしまい。長かった。

 

 

広告を非表示にする