1
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

BUGS Example "Beetle" を解く

Last updated at Posted at 2015-08-25

はじめに

BUGS (Beyesian inference using Gibbs Sampling) のWebサイトにある例題のうち,今回"Beetle"をPythonで解いてみる. ”Beetle"の概要だが,Carbon Disulphide(二硫化炭素)環境下にカブトムシを入れ,ガスの濃度に応じてカブトムシの死亡率(5時間後に何匹死んだか)についてデータ分析するというものである.与えられるのはカウントデータなので,二項分布で離散的な統計モデルを構成する問題となる.

ツールとしては,まず初めにstatsmodelsを用いて分析(GLMの計算)を行い,次にPyMC3で計算を行った.

与えられているBliss(1935)のデータセットを参照する.

Table. Bliss's data

Concentration(xi) Num. of beetles(ni) Num. of killed (ri)
1.6907 59 6
1.7242 60 13
1.7552 62 18
1.7842 56 28
1.8113 63 52
1.8369 59 52
1.8610 62 61
1.8838 60 60

これと一緒に掲載されていたBUGS scriptは,以下の通りである.

model
    {
       for( i in 1 : N ) {
          r[i] ~ dbin(p[i],n[i])
          logit(p[i]) <- alpha.star + beta * (x[i] - mean(x[]))
          rhat[i] <- n[i] * p[i]
         culmative.r[i] <- culmative(r[i], r[i])
       }
       alpha <- alpha.star - beta * mean(x[])
       beta ~ dnorm(0.0,0.001)
       alpha.star ~ dnorm(0.0,0.001)   
    }

やってしまいがちな「ダメ」方法

今回のようなロジスティック回帰の問題では,求めたいモデルが縦軸を 0<p<1 とした確率の関数であるのに対し,与えられるのが「総個体数」と「死亡数」というカウントデータである.よくやってしまうのが,最初にデータの前加工として比率を算出し,それにfitするモデルを求めるというやり方である.これをcodeにすると以下のようになる.

n_betl = np.array([59, 60, 62, 56, 63, 59, 62, 60])
y_betl = np.array([6, 13, 18, 28, 52, 53, 61, 60])
y1 = y_betl *1.0 / n_betl         # calculated the proportion of y : n

x1 = np.array([1.6907, 1.7242, 1.7552, 1.7842,
    1.8113, 1.8369, 1.8610, 1.8839])
x1a = sm.add_constant(x1)

model_g1 = sm.GLM(y1, x1a, family=sm.families.Binomial(sm.families.links.logit))
res1 = model_g1.fit()
print res1.summary()

このようなcodeを書いて出来たと思っていたが,「データ解析のための統計モデリング入門」を開くと久保先生に怒られてしまうことが分かる.

「割り算値(観測データ/観測データ)による回帰や相関分析は現在でも見かける方法です.(中略)このように変形した観測地を統計モデルの応答関数にするのは,不必要であるばかりでなく,場合によっては間違った結果を導きかねないものです.」
(ごもっとです.この本からはいつもお叱り,注意をいただいています.)

割り算値によらない,カウントデータをモデル化する手法として以下の2つが考えられる.

  • Binomial distrubutionとパラメータをきちんと引数としてGLM(=Generalized Linear Model) functionに渡して回帰計算を行う.
  • 集計されているデータを,一つ一つが(Beetleが死亡していない or 死亡を示す)ラベルを持つデータに伸長し,そのデータを回帰分析にかける.

明らかに,前者のやり方がスマートである.因みに"R"のcodeは次の通りである.

# R's model expression
   beetle.glm <- glm(cbind(y, n-y) ~ dose(=x1), family=binomial(link=logit))

R codeで使っている "cbind(y,n-y) ~ dose" のような式表記をうまくPython用に置き換えることができなかった.statsmodels.formula で前記のような統計モデル式を取り扱うことは可能のようだが,情報不足で「組み合わせ」(combination... nCy)を表現する(statsmodels.formula用の)式を作成できず.(質問サイトなどいろいろ調べましたが,結局断念となりました..)仕方がないので,後者のやり方をとることにした.

statsmodelsによるLogistic Regression

書き直したcodeは,以下の通りである.

import statsmodels.api as sm
import statsmodels.formula.api as smf

# prep data arrays by expanding tabulated data
resp_expand = np.array([])
x1_expand = np.array([])
for i in range(len(n_betl)):
ni = n_betl[i] ; yi = y_betl[i]
resp_expand = np.append(resp_expand, np.repeat(0, (ni-yi)))
resp_expand = np.append(resp_expand, np.repeat(1, yi))
x1_expand = np.append(x1_expand, np.repeat(x1[i], ni))

mydf = pd.DataFrame()
mydf['x1_ex'] = x1_expand
mydf['resp'] = resp_expand

myformula = 'resp ~ x1_ex'
model_g1 = smf.logit(myformula, data=mydf)
# model_g1 = smf.glm(myformula, data=mydf, family=sm.families.Binomial()) これでもOK.
result1 = model_g1.fit()
print result1.summary()


下から4行めの通り logit() を使い,うまく計算できた.summary() を示す.

                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                   resp   No. Observations:                  481
Model:                            GLM   Df Residuals:                      479
Model Family:                Binomial   Df Model:                            1
Link Function:                  logit   Scale:                             1.0
Method:                          IRLS   Log-Likelihood:                -186.24
Date:                Tue, 25 Aug 2015   Deviance:                       372.47
Time:                        09:46:41   Pearson chi2:                     436.
No. Iterations:                     8                                         
==============================================================================
                 coef    std err          z      P>|z|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept    -60.7175      5.181    -11.720      0.000       -70.871   -50.563
x1_ex         34.2703      2.912     11.768      0.000        28.563    39.978
==============================================================================

PyMC3 によるLogistic Regression

PyMC3ではベイズ推定の作法に従い,事前分布(prior distribution)と尤度関数(likelihood)を与えてモデル化し,それをMCMC (Marcov chain Monte Carlo)で解くというやり方になるが,ここでは上述のような「割り算値を使ってしまう」こともなく,自然にBinomial distributionを記述してモデル化できてしまう.
(参照したPyMC Tutorial コードがきちんとできているおかげとも言えますが...)

以下のようなcodeを用意し,シミュレーションを実行した.

import pymc3 as pm

n_betl = np.array([59, 60, 62, 56, 63, 59, 62, 60])
y_betl = np.array([6, 13, 18, 28, 52, 53, 61, 60])
x1 = np.array([1.6907, 1.7242, 1.7552, 1.7842,
    1.8113, 1.8369, 1.8610, 1.8839])

mymodel1 = pm.Model()
def invlogit(x):
    return pm.exp(x) / (1 + pm.exp(x))

with mymodel1:  # model definition
    mytheta = pm.Normal('theta', mu=0, sd=32, shape=2)   # 2-dim array
    p = invlogit(mytheta[0] + mytheta[1] * x1)
    y_obs = pm.Binomial('y_obs', n=n_betl, p=p, observed=y_betl)
    
with mymodel1:  # running simulation
    start = pm.find_MAP()
    step = pm.NUTS()
    trace1 = pm.sample(10000, step, start=start)

pm.summary(trace1)
pm.traceplot(trace1)

traceplot()は,以下の通り.(初稿にあったsummary()は削除しました.)

Fig. Traceplot (Logistic Regression)
beetle_traceplot1.png

今回,機械的に10000 step と設定してMCMC計算を行ったが,一応,Geweke法収束を確認する.(因みに計算時間は,どのMCMC計算も 30秒前後でした.)

Fig. Geweke Diagnostic
beetle_geweke01.png

収束に関し,問題なし.

Trace結果を用いて作図したFitting Plotは以下.

Fig. Logistic Regression plot
beetle_model1_fitting.png

計算結果比較

(追記:初稿では,probitモデル,c-loglogモデルの説明を記載していましたが,記事が冗長でしたので省略することにしました.)

Table 計算結果の比較

Method Link func. theta1 (mean) theta1 (SD) theta2 (mean) theta2 (SD)
BUGS (reference) logit -60.79 5.147 34.31 2.893
BUGS (reference) probit -35.04 2.646 19.79 1.488
BUGS (reference) c-loglog -39.77 3.211 22.15 1.788
GLM, statsmodels logit -60.71 - 34.27 -
MCMC, PyMC3 logit -59.56 4.911 33.62 2.761

3者間で(BUGS(logit) vs. statsmodels vs. PyMC3)で概略,整合性がとれているように見える.(今回は,PyMC3よりstatsmodelsの方でだいぶ苦労されられました.)

参考文献 (web site)

1
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
1
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?