Stanで統計モデリングを学ぶ(4): とりあえず階層ベイズモデルを試してみる（基本編）

だいぶ間が空いちゃいましたね（汗）。これまでの記事で大体Stanで何ができるか分かったので、ぼちぼちStanらしいことをやってみようと思います。一応過去記事のリスト出しておきますので、良かったら復習も兼ねてお読みください。

Stanらしいことをやってみようと言っても、美味しいところはとっくの昔に@berobero11さんが持って行ってしまわれているので、僕は延々とチマチマ基本的なところを自分の理解のスピードの範囲でやっていくつもりです。

ということで、今回はGLMMのベイズ化という久保先生の緑本と同じテーマでやってみようかと思います。ついでなので、データも久保先生の緑本のサンプルから完全に拝借してしまいますｗ　今回は第7章のデータに対して、第10章の内容でやっていきます。

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

作者: 久保拓弥
出版社/メーカー: 岩波書店
発売日: 2012/05/19
メディア: 単行本
購入: 16人クリック: 163回
この商品を含むブログ (19件) を見る

もっと言ってしまうと、これって思いっきりITO Hirokiさんのブログで既にやられているネタなので、ぶっちゃけ全く新鮮味はありません（泣）。

R: Stanをためす：Taglibro de H：So-netブログ

ただしITOさんのStanコードはわずかに久保先生のBUGSコードと整合していないので、僕の手元でちょっと試行錯誤してできるだけBUGSコード側に整合させてみたものを紹介しておきます。

階層ベイズって何ぞ

実は結構な数のテキストが「階層ベイズ＝個体差のあるモデル」みたいな説明をしていると思うんですが、これって別にGLMMでもできるのであまり意味のある説明じゃないんですよね。これは多分丹後-Becque本のpp.148-149にある説明を踏まえるのがベストでしょう。

ベイジアン統計解析の実際 (医学統計学シリーズ)

作者: 丹後俊郎,Taeko Becque
出版社/メーカー: 朝倉書店
発売日: 2011/10/14
メディア: 単行本
クリック: 8回
この商品を含むブログを見る

頻度論者の母数モデル(fixed-effects model)では、それぞれの研究結果が共通の平均 $\theta$ をもつ
$y_i \sim N(\theta, s^2_i)$
と仮定する。一方、変量モデル(random-effects model)では、それぞれ異なる平均 $\theta_i$ をもち、その平均が全体平均 $\mu$ の周りに分布する
$y_i \sim N(\theta_i, s^2_i)$
$\theta_i \sim N(\mu, \sigma^2_i)$
と仮定するモデルである。これに対し、ベイジアン階層的モデルでは変量モデルの $\theta_i$ の確率分布のパラメータ $\mu, \sigma^2_\theta$ にさらなる確率分布（事前分布）を仮定する。例えば、
$\mu \sim N(0,100)$
$1/{\sigma^2_i} \sim Gamma(0.001,0,001)$
などと仮定するモデルである。この例のモデルは、二つのレベルの階層構造をもつという。すなわち

第1レベル： $\theta_i$ に確率分布が仮定されている

第2レベル：第1レベルの確率分布のパラメータ $\mu, \sigma^2_i$ に確率分布が仮定されている

という構造を意味するのである。

ということで、端的に言えば頻度主義における未知母数の変量モデルにおけるパラメータに対して、それ自体が何かまた別の事前分布に従うという、「未知母数→変量モデル→そのまたパラメータという順番で確率分布が階層的になっている」モデルのことを指している、と理解する方が良さそうです。

この辺のことが分かってないと、いかなStanコーディングを覚えても全く階層ベイズモデルを構築できないということになりかねないので、上に挙げたような「流れ」を覚えておくことが大事なんじゃないかなぁ、と初学者としては思ってます。

実際に階層ベイズをStanでやってみる

さて、四の五の理屈を垂れるのはこの辺にして、実際にStanで階層ベイズモデルを推定してみようと思います。お題は既に冒頭で触れたように、久保先生の緑本第10章です。

まずStanコードを組みます。以下のコードをテキストエディタなどで書いて、"hbayes.stan"という名前で保存しておきます。ちなみに、ただのロジスティック回帰であればStanではbernoulli()関数で簡単に推定できるんですが、階層モデルにすると何故かぐちゃぐちゃになったりすることが多いようで。。。一旦transformed parametersブロックや、modelブロックのサンプリング手前部分などで手書きでロジット変換してから、改めてnormal()関数でサンプリングした方が安定することが多いようです（僕の経験から言っても）。

data {
	int<lower=0> N; // サンプルサイズ
	int<lower=0> y[N]; // 種子8個当たりの生存数（目的変数）
}
parameters {
	real beta; // 全個体共通のロジスティック偏回帰係数
	real r[N]; // 個体差
	real<lower=0> s; // 個体差のばらつき
}
transformed parameters {
	real q[N];
	for (i in 1:N)
		q[i] <- inv_logit(beta+r[i]); // 生存確率を個体差でロジット変換
}
model {
	for (i in 1:N)
		y[i] ~ binomial(8,q[i]); // 二項分布で生存確率をモデリング
	beta~normal(0,1.0e-4); // ロジスティック偏回帰係数の無情報事前分布
	for (i in 1:N)
		r[i]~normal(0,1/(s*s)); // 個体差の階層事前分布
	s~uniform(0,1.0e+4); // r[i]を表現するための無情報事前分布
}

事前分布の設定の仕方などは基本的に緑本のやり方を踏襲してます。ここで使われているデータは各個体から8個の種子を選び、その中で生存していた種子の数を記録したというもの。そこに個体差があると仮定するのがここでの階層ベイズモデルのポイントです。緑本p.229に分かりやすい図がありますが、要は

二項分布で生存確率をモデリングする
生存確率は「個体差」＋共通偏回帰係数のロジット変換で表される
共通偏回帰係数はよく分からないので標準偏差 $10^2$ の正規分布を無情報事前分布とする
個体差の階層事前分布は適当に正規分布で置く
その個体差の階層事前分布（＝正規分布）の標準偏差は $0 < s < 10^4$ の一様分布たる無情報事前分布で与える

というロジックから成り立っているモデルです。全く分からないところは広い一様分布を無情報事前分布として突っ込み、ある程度正規分布が仮定できるものは無情報事前分布にせよ階層事前分布にせよ正規分布を突っ込む、という感じですね。では、第7章のデータを直接読み込んで、以下のようにStanを回してみましょう。

> d <- read.csv(url("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/hbm/data7a.csv"))
> head(d)
  id y
1  1 0
2  2 2
3  3 7
4  4 8
5  5 1
6  6 7
> dat<-list(N=nrow(d),y=d$y)
> d.fit<-stan(file='hbayes.stan',data=dat,iter=1000,chains=4)

TRANSLATING MODEL 'hbayes' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'hbayes' NOW.
cygwin warning:
  MS-DOS style path detected: C:/PROGRA~1/R/R-30~1.2/etc/i386/Makeconf
  Preferred POSIX equivalent is: /cygdrive/c/PROGRA~1/R/R-30~1.2/etc/i386/Makeconf
  CYGWIN environment variable option "nodosfilewarning" turns off this warning.
  Consult the user's guide for more details about POSIX paths:
    http://cygwin.com/cygwin-ug-net/using.html#using-pathnames
C:/Program Files/R/R-3.0.2/library/rstan/include//stansrc/stan/agrad/rev/var_stack.hpp:49:17: warning: 'void stan::agrad::free_memory()' defined but not used [-Wunused-function]
SAMPLING FOR MODEL 'hbayes' NOW (CHAIN 1).
Iteration: 1000 / 1000 [100%]  (Sampling)
Elapsed Time: 12.313 seconds (Warm-up)
              3.601 seconds (Sampling)
              15.914 seconds (Total)

SAMPLING FOR MODEL 'hbayes' NOW (CHAIN 2).
Iteration: 1000 / 1000 [100%]  (Sampling)
Elapsed Time: 11.74 seconds (Warm-up)
              3.534 seconds (Sampling)
              15.274 seconds (Total)

SAMPLING FOR MODEL 'hbayes' NOW (CHAIN 3).
Iteration: 1000 / 1000 [100%]  (Sampling)
Elapsed Time: 12.26 seconds (Warm-up)
              6.143 seconds (Sampling)
              18.403 seconds (Total)

SAMPLING FOR MODEL 'hbayes' NOW (CHAIN 4).
Iteration: 1000 / 1000 [100%]  (Sampling)
Elapsed Time: 16.227 seconds (Warm-up)
              3.512 seconds (Sampling)
              19.739 seconds (Total)

> plot(d.fit)
> d.fit
# 略

f:id:TJO:20140501155537p:plain

結構個体差がばらついているのが、プロットから見て取れますね。緑本p.231の図10.4みたいなのを描いてみたかったんですが、StanのMCサンプルから描くコードがすぐには思い付かなかったので一旦置いときます。WinBUGSでやる場合は久保先生のコード例があるのでご参考までに。

ちなみにd.fitとただ打てば個々のパラメータから推定された平均やらパーセンタイル点やらが出ます。また、extract()してdensity()して最大値を計算すれば、想定される個々のパラメータの解も得られます。緑本p.231の図10.3に出てくる $\beta, s, r_1, r_2$ の事後分布を{coda}パッケージを用いてプロットしたものを下に示しておきます。

f:id:TJO:20140501180857p:plain

うーん、sだけ何かおかしい気がするんですが（汗）他は大体緑本の通りになっているかと思います。

後日判明したポイント

BUGSから移植するに当たり、一部事前分布の指定が間違ってました。。。Stanコードを以下のように修正すれば、sも正確に推定されます。

data {
	int<lower=0> N; // サンプルサイズ
	int<lower=0> y[N]; // 種子8個当たりの生存数（目的変数）
}
parameters {
	real beta; // 全個体共通のロジスティック偏回帰係数
	real r[N]; // 個体差
	real<lower=0> s; // 個体差のばらつき
}
transformed parameters {
	real q[N];
	for (i in 1:N)
		q[i] <- inv_logit(beta+r[i]); // 生存確率を個体差でロジット変換
}
model {
	for (i in 1:N)
		y[i] ~ binomial(8,q[i]); // 二項分布で生存確率をモデリング
	beta~normal(0,1.0e+2); // ロジスティック偏回帰係数の無情報事前分布
	for (i in 1:N)
		r[i]~normal(0,s); // 個体差の階層事前分布
	s~uniform(0,1.0e+4); // r[i]を表現するための無情報事前分布
}

これでsの事後分布の最頻値は3当たりに落ち付くと思います。

個人的に階層ベイズをやる時にイメージしていること

表現がおかしいかもしれませんが、実は僕が階層ベイズをやる時は実データに対して階層分布たちという「クッション」を置く、というイメージを持ってやってます。

というのは、一般に通常のGLMであまり正確に固定効果モデルの未知母数が推定できないケースというのは、個体差も含めた何かしらの「固定効果モデルの範囲を逸脱した大きなばらつきがある」ということであり、このばらつきを何とかして緩和してやることが必要だからです。

もちろん、その「ばらつき」を同定することも大事なんですが、もしその個体差なりを反映した「ばらつき」ではなく、あくまでもシステム全体の共通的・普遍的なパラメータにしか興味がないのであれば、あえて階層モデルにそういう「ばらつき」を吸収させ、その上で真に求める共通・普遍パラメータを求めるという取り組みが必要になると思うのです。

例えば、上の緑本第10章のケースだと「生存確率を二項分布でモデリング」がまさにその「クッション」に当たるわけです。生存確率の実測値は様々な「ばらつき」によって上下動するので、これを大きな二項分布という確率パラメータなるクッションで包み込んでしまえば、サンプリング時にあまり振り回されずに済む。。。というイメージです。

そういう意味では、不要なばらつきを吸収させる「クッション」としての階層モデルというイメージを持って、階層モデルを使うことが僕は多かったりします。もしかしたらおかしな理解かもしれませんが、そうすることで求めたいパラメータのMCサンプルの収束が安定したりすることもあるので、個人的には結構有用なイメージの持ち方かもと思ってます。。。