Pythonで実装しながら緑本を学ぶ (第10章階層ベイズモデル)

データ解析のための統計モデリング入門(通称、緑本)を読み進めています。
述べられている理論を整理しつつ、Rでの実装をPythonに置き換えた際のポイントなども深掘りしていきます。

今回は第10章です。PyMC3を使って階層ベイズモデルを表現します。実装は以下で公開しています。

introduction_to_machine_learning_with_python/chapter10.ipynb at master · ohke/introduction_to_machine_learning_with_python · GitHub

10 GLMのベイズモデル化と事後分布の推定 -GLMMのベイズモデル化-

階層事前分布を使って一般化線形混合モデル(GLMM)を階層ベイズモデルとして扱うことが、本章のテーマです。

10.1 例題：個体差と生存種子数(個体差あり)

1個体あたりの種子数8個を持つ架空植物100個体について、生存種子数を計測した。個体 $i$ における生存種子数を $y_i$ とする。

8個の内 $y_i$ 個というデータですので、二項分布が使えそうです。
しかし個体に由来するランダム効果によって過分散が生じています。平均4.03に対して分散9.93となっており、二項分布が期待する分散2.00(=8×(4.03/8)×(1-(4.03/8)))よりも、かなり大きいです。

f:id:ohke:20180310095005p:plain

10.2 GLMMの階層ベイズモデル化

第7章(こちらの投稿)では、GLMにこうした原因不明の個体差を組み込んだGLMMを導入し、パラメータを最尤推定で求めた。
パラメータが増えると最尤推定では困難になるため、MCMCサンプリングによる事後分布推定を行いたい。そこで、GLMMのベイズモデル化を行う。

まずリンク関数と線形予測子を決めます。二項分布ですので、リンク関数にはロジットを使います(第6章も参照)。 $r_i$ が個体差です。

$$ logit(q_i)=\beta+r_i $$

8章・9章で追ってきたとおり、事後分布は $p({\bf Y} \mid \beta, \{r_i\}) \times \mbox{事前分布}$ に比例します。

データ ${\bf Y }$ が得られる条件付き確率(尤度)は全個体の二項分布の積です。

$$ p({\bf Y} \mid \beta, {r_i})=\prod_i \binom{8}{y_i} q_i^{y_i}(1-q_i)^{8-y_i} $$

$\beta$ の事前分布は、特に制約がないので平均0・標準偏差100の正規分布に従うこととする(無情報事前分布)。

$$ p(\beta)=\frac{1}{\sqrt{2\pi \times 100^2}} \exp(\frac{-\beta^2}{2 \times 100^2}) $$

個体差 $r_i$ は、個体差のばらつきを表す変数 $s$ を導入し、平均0・標準偏差 $s$ の正規分布に従うことにする。

$$ p(r_i \mid s)=\frac{1}{\sqrt{2\pi s^2}} \exp(\frac{-r_i^2}{2s^2}) $$

$s$ は正であれば良いので、0から10000の連続一様分布に従う無情報事前分布とする。

$$ p(s)=U(0, 10^4) $$

このように、別の事前分布 $p(s)$ によって定まる事前分布 $p(r_i \mid s)$ を階層事前分布と呼ばれます。これに対して、 $s$ を超パラメータ、 $p(s)$ を超事前分布)と言います。

階層事前分布を使うベイズモデルが、階層ベイズです。

10.3 階層ベイズモデルの推定・予測

PyMC3を使い、10.2のモデルを構築します。超事前分布 $p(s)$ 、階層事前分布 $p(r \mid s)$ の順に作っていることを確認してください。

Normalで正規分布、Uniformで連続一様分布、Binomialで二項分布を設定できます
- rは個体数分必要ですので、shape=len(data.y)としています
- invlogitでロジット関数の逆関数(つまりロジスティック関数)を二項分布のpに設定しています
NUTSで1600点(最初の100点は捨てる)をサンプリングしたサンプリング列を3つ作る
- 合計4500点をサンプリング

import pymc3 as pm

# モデルの作成
with pm.Model() as model:
    # βの事前分布をN(0, 100)の正規分布で設定(無情報事前分布)
    beta = pm.Normal('beta', mu=0, sd=100)
    
    # 超パラメータsの(超)事前分布をU(0, 10000)の連続一様分布で設定(無情報事前分布)
    s = pm.Uniform('s', lower=0, upper=10000)
    
    # パラメータrの事前分布をN(0, s)の正規分布で設定(階層事前分布)
    r = pm.Normal('r', mu=0, sd=s, shape=len(data.y))
    
    # ロジットリンク関数を設定し、二項分布で推定する
    y = pm.Binomial('y', n=8, p=pm.invlogit(beta + r), observed=data.y.values)

# サンプリング
with model:
    # NUTSで101個目からサンプルを取得するチェインを3つ作る
    trace = pm.sample(1500, step=pm.NUTS(), tune=100, njobs=3, random_seed=0)
    
_ = pm.traceplot(trace)
pm.summary(trace).loc[['beta', 's']]

各パラメータの事後分布(95%信用区間)を見ると $\beta$ は-0.600〜0.723、 $s$ は0.010〜3.802となっており、いずれも $\hat{R}$ は十分1に近い(収束している)。

f:id:ohke:20180310111301p:plain

$r_i$ の事後分布(下図中段左)を見ると、 $r_i$ 毎に異なる分布となっており、個体差が表現されていることがわかる。

f:id:ohke:20180310111413p:plain

パラメータの事後分布から、生存確率を予測してみます。

3つのサンプル列の平均値を求める( $\beta$ と $s$ のそれぞれ1500点ずつ)
$\beta$ と $s$ を1つずつ取り出し、 $s$ から100個体分の $\{r_i \}$ を正規分布で生成する
$\beta$ と $\{r_i \}$ から二項分布で種子数の生存確率を求める
1500サンプル・100個体分の生存確率の平均値を計算する

こうして得られた生存確率から、100個体を調べた時の生存種子数の頻度をプロットしたのが以下となります。過分散を予測できていることがわかります。

f:id:ohke:20180310115057p:plain

10.4 ベイズモデルで使うさまざまな事前分布

事前分布の決め方は3パターン。

ある分布を信じて、主観的に決める (基本的には使わない)
- 他のデータセットや統計モデルで得られた分布を使う場合など
無情報事前分布
階層事前分布

無情報事前分布は、個体ごとに自由に決めてしまうと、フルモデルと同等になるデメリットがある。そのため10.2のように、個体ごとに異なるパラメータなどは階層事前分布とし、階層事前分布の(超)パラメータのみを無情報事前分布とするのが基本戦略となる。
したがって、パラメータの影響が局所的な場合は階層事前分布、大域的な場合は無情報事前分布と使い分ける。

10.5 個体差＋場所差の階層ベイズモデル

個体差に加えて、場所差を組み込んだモデルを見ていきます。

架空植物100個体の種子数を計測し、施肥処理と種子数の影響について調べたい、というケースを想定します。
架空植物は10個の植木鉢A〜Jで10個体ずつ植えられてます。植木鉢A〜Eの50個体は施肥処理無し('C')、F〜Jの50個体は施肥処理有り('T')で育てられています。

個体ごとの種子数をプロットしてみますと、施肥処理無し(左側50個体)の方が若干種子数が多いように見えます。

f:id:ohke:20180310121717p:plain

しかし、植木鉢ごとの種子数を箱ひげ図でプロットしてみますと、種子数の分布が異なっており、植木鉢による差が少なからずあるものと思われます。

f:id:ohke:20180310121731p:plain

個体差と植木鉢差を考慮したモデルを定義していきます。

種子数は離散値・上限なしのため、ポアソン分布を使います。
個体 $i$ におけるパラメータ $\lambda_i$ は、線形予測子と対数リンク関数でモデルを定義します。

$$ \log \lambda_i = \beta_1 + \beta_2 f_i + r_i + r_{j(i)} $$

は、平均0・標準偏差100のの正規分布
- $f_i$ は個体 $i$ の施肥処理の有無(有りの場合に1となる)
$r_i$ は個体差で、平均0・標準偏差 $s$ の正規分布 (階層事前分布)
は植木鉢差で、平均0・標準偏差の正規分布 (階層事前分布)
- パラメータ数は10個( $j$ は0〜9)
$s, s_p$ は0〜10000の連続一様分布

PyMC3での実装は以下のとおりです。

with pm.Model() as model:
    # β1とβ2の事前分布をN(0, 100)の正規分布で設定(無情報事前分布)
    beta1 = pm.Normal('beta1', mu=0, sd=100)
    beta2 = pm.Normal('beta2', mu=0, sd=100)
    
    # sとspの(超)事前分布をU(0, 10000)の連続一様分布で設定(無情報事前分布)
    s = pm.Uniform('s', lower=0, upper=10000)
    sp = pm.Uniform('sp', lower=0, upper=10000)
    
    # 個体差パラメータrの事前分布をN(0, s)の正規分布で設定(階層事前分布)
    r = pm.Normal('r', mu=0, sd=s, shape=len(data.y))
    
    # 場所差パラメータrpの事前分布をN(0, sp)の正規分布で設定(階層事前分布)
    rp = pm.Normal('rp', mu=0, sd=sp, shape=len(data.POT.unique()))
    
    # ログ関数を設定し、ポアソン分布で推定する
    y = pm.Poisson('y', mu=np.exp(beta1 + beta2 * data.F + r + rp[data.POT]), observed=data.y.values)
    
# HMCでサンプリング
with model:
    # 101個目からサンプルを取得するチェインを3つ作る
    trace = pm.sample(1500, step=pm.HamiltonianMC(), tune=100, njobs=3, random_seed=0)

pm.summary(trace).loc[['beta1', 'beta2', 's', 'sp']]