井出草平の研究ノート

ベイジアンロジスティック回帰分析・その2

HardinとHilbeのStata本にあるベイジアンロジスティック回帰分析の続き。前回のエントリはこちら

Generalized Linear Models and Extensions: Fourth Edition

Generalized Linear Models and Extensions: Fourth Edition

メトロポリスヘイスティングス

メトロポリスヘイスティングス法のやり方である。

前回も同じことを行っているが、データの読み込みと変数加工(センタリング)もう一度行っておこう。

. use http://www.stata-press.com/data/hh4/rwm1984
. quietly summarize docvis, meanonly
. generate cdoc = docvis - r(mean)
. quietly summarize age, meanonly
. generate cage = age - r(mean)

医者の往診の変数docvisと年齢ageをセンタリングしている。

MH法では変数の設定が必要である。priorは従属変数を指定し、その後にコンマをうち、分布とパラメーターを指定する構文が必要だ。以下のコマンドの2行目のところである。

. bayesmh outwork female kids cdoc cage, likelihood(logit)
> prior({outwork:}, normal (0,10000)) rseed (2017)
Model summary
------------------------------------------------------------------------------
Likelihood:
  outwork ~ logit(xb_outwork)

Prior:
  {outwork:female kids cdoc cage _cons} ~ normal(0,10000)                  (1)
------------------------------------------------------------------------------
(1) Parameters are elements of the linear form xb_outwork.

Bayesian logistic regression                     MCMC iterations  =     12,500
Random-walk Metropolis-Hastings sampling         Burn-in          =      2,500
                                                 MCMC sample size =     10,000
                                                 Number of obs    =      3,874
                                                 Acceptance rate  =      .2363
                                                 Efficiency:  min =     .02665
                                                              avg =     .05244
Log marginal-likelihood = -2001.0895                          max =       .101

------------------------------------------------------------------------------
             |                                                Equal-tailed
     outwork |      Mean   Std. Dev.     MCSE     Median  [95% Cred. Interval]
-------------+----------------------------------------------------------------
      female |  2.260566   .0803699   .002529    2.25819   2.103102   2.421029
        kids |  .3562706   .0893508   .004221   .3583318   .1783027   .5281715
        cdoc |   .024483   .0061758   .000329   .0244374   .0130205   .0369231
        cage |  .0545722   .0040688   .000174   .0547603   .0464829   .0622289
       _cons |  -2.01027    .082349   .005044   -2.00927  -2.172552  -1.852391
------------------------------------------------------------------------------

プロット

プロットは非収束の徴候が確認できる。

従属変数と各独立変数でプロットを作成する。
従属変数であるoutworkcdocでグラフを作成した。

Stataでは引数なしでグラフ診断を走らせると、4タイプのグラフが描画される。

. bayesgraph diagnostics {outwork: cage}

f:id:iDES:20191111183628p:plain

トレースプロットは各反復におけるパラメータ推定値をグラフ化する。ラグ項との自己相関である。ヒストグラムは事前分布と一致しているか確認する。outworkcdocでのグラフ診断は正規分布に近いが、世正規分布が必須というわけではないようだ。

outworkとセンタリングした年齢の変数であるcageのグラフ診断も走らせてみる。

f:id:iDES:20191111183616p:plain

cdocほどは良いグラフになっていないので、サンプリングに反復を追加し、値の選択を3から5減らすと自己相関を減らすことができるとのこと。

ESSコマンド

ESSコマンドでは、モデル内の指定された各パラメータについて、有効なサンプルサイズ、相関時間、および効率を表示できる。

. bayesstats ess

Efficiency summaries    MCMC sample size =    10,000
                        Efficiency:  min =    .02665
                                     avg =    .05244
                                     max =      .101

----------------------------------------------------
    outwork  |        ESS   Corr. time    Efficiency
-------------+--------------------------------------
      female |    1009.74         9.90        0.1010
        kids |     448.19        22.31        0.0448
        cdoc |     352.37        28.38        0.0352
        cage |     545.17        18.34        0.0545
       _cons |     266.54        37.52        0.0267
----------------------------------------------------

ベイズモデルでは、モデルパラメータの事後平均は、パラメータの事後分布のMCMCサンプリングの平均によって決定される。推定されたサンプル平均の標準誤差はに比例する。

ESS値が低いと注意が必要である。理想的には、ESS値はMCMCサンプルサイズに近い値をとる。効率は、ESS値にまたはを乗算することによって決定される。女性のESS値は1009.74であり、効率は0.101である。

Gibbs法によるサンプリングを行わない限り、平均の効率値0.01が許容できるモデルの最小効率である。0.1以上の値はランダムウォークメトロポリスアルゴリズム(random-walk Metropolis algorithm)にとって効率的であるとされているとのこと。

相関時間(correlation time)は効率の逆数である。自己相関が低い推定遅延と考えられ、値は小さい方が良い。