井出草平の研究ノート

トップ > Stata > ベイジアンロジスティック回帰分析[Stata]

ベイジアンロジスティック回帰分析[Stata]

Stata ベイズ統計学

HardinとHilbeが書いたglmのStata本をちらちらと読んでいて少し気になったので、自分でも分析をしてみた。

Generalized Linear Models and Extensions: Fourth Edition

Generalized Linear Models and Extensions: Fourth Edition

この分析の出典はこの本の20章2節(20.2)である。

データ

rwm1984.dtaというデータを用いる。1984年のドイツの医療改革データである。下記の変数が含まれている。

変数 ラベル
docvis MD visits/year
hospvis Days in hospital/year
edlevel Level of education
age Ages 25-64
outwork 1=not working; 0=working
female 1=femaie; 0=male
married Married=1; Single=0
kids 1=have children; 0=no children
hhninc Household income (Marks,OECD wgt)
educ Highest years of education
self 1=Self employed; 0=other
edlevel1 edlevel==Not HS grad
edlevel2 edlevel==HS grad
edlevel3 edlevel==Coll/Univ
edlevel4 edlevel==Grad School

ロジスティック回帰分析

まずは普通のロジスティック回帰分析を行う。つまり頻度主義でP値が出力されるものである。

. 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)
. glm outwork female kids cdoc cage, family(binomial) link(logit) nolog

Generalized linear models                         Number of obs   =      3,874
Optimization     : ML                             Residual df     =      3,869
                                                  Scale parameter =          1
Deviance         =   3918.21514                   (1/df) Deviance =    1.01272
Pearson          =  4205.000001                   (1/df) Pearson  =   1.086844

Variance function: V(u) = u*(1-u)                 [Bernoulli]
Link function    : g(u) = ln(u/(1-u))             [Logit]

                                                  AIC             =   1.013995
Log likelihood   =  -1959.10757                   BIC             =  -28047.63

------------------------------------------------------------------------------
             |                 OIM
     outwork |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      female |   2.256804   .0827624    27.27   0.000     2.094593    2.419016
        kids |   .3579762   .0899639     3.98   0.000     .1816503    .5343022
        cdoc |   .0244323   .0062628     3.90   0.000     .0121574    .0367071
        cage |   .0543791   .0041591    13.07   0.000     .0462274    .0625308
       _cons |  -2.010276   .0810901   -24.79   0.000     -2.16921   -1.851342
------------------------------------------------------------------------------

glmの前に変数を加工している行が４行入っている。quietlyは出力抑止のコマンドで、- r(mean)はセンタリングをするコマンドである。centeringのcがついたcdocとcageがglmの独立変数になっている。この処理はベイジアンのために行っており、通常のglmには不要である。結果の比較のために通常のglmでも変数を利用している。

ベイジアンロジスティック回帰分析

次にベイジアンでの分析である。変数はさきほど作ったものと同一である。bayesを接頭語にした場合、nologをつけるとエラーがでるので削除するようだ。デフォルトでは無情報事前分布、一様分布と一様分布と分散normal（0、10000）に設定されている。デフォルトで特に問題がないので、そのまま実行する。シードは2017というおそらく執筆時点の年数を入れていて、別の値でも構わない。

. bayes, rseed (2017): glm outwork female kids cdoc cage, family(binomial)


Model summary
------------------------------------------------------------------------------
Likelihood:
  outwork ~ glm(xb_outwork)

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

Bayesian generalized linear models               MCMC iterations  =     12,500
Random-walk Metropolis-Hastings sampling         Burn-in          =      2,500
                                                 MCMC sample size =     10,000
Family : Bernoulli                               Number of obs    =      3,874
Link   : logit                                   Scale parameter  =          1
                                                 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
------------------------------------------------------------------------------
Note: Default priors are used for model parameters.

結果は基本的には変わらない。

逸脱度情報量規準(DIC)と限界対数尤度値(marginal log-likelihood value)

他のモデルと比較する際はicコマンドが使える。

. bayesstats ic

Bayesian information criteria

----------------------------------------------
             |       DIC    log(ML)    log(BF)
-------------+--------------------------------
      active |  3928.038   -2001.09          .
----------------------------------------------
Note: Marginal likelihood (ML) is computed
      using Laplace-Metropolis approximation.

DICは3928.038、 log(ML)が限界対数尤度値であり、-2001.09である。DICはAICと基本的に同じようなものであるが、事後分布がおおむね多変量正規分布である時に有効な指標となるので、少し注意が必要かもしれない。

余談

このエントリのタイトルはベイジアンロジスティック回帰分析だが、これはBayesian Logistic Regressionを翻訳したものだ。日本語の定訳を探してみたが見つからなかった。Bayesianは翻訳がしにくい用語である。そのままベイジアンであったり(例えばベイジアン・ネットワーク)、ベイズ的であったり(ベイズ的アプローチ)、ベイズ式であったり(ベイズ式分類法)、訳がころころと変わる。

このデータを用いた分析は著者のHilbeの別の本で言及があるそうだ。

Modeling Count Data (English Edition)

Modeling Count Data (English Edition)