Morris, M., Wheeler-Martin, K., Simpson, D., Mooney, S.J., Gelman, A., DiMaggil, C. (2019) Bayesian hierarchical spatial models: Implementing the Besag York Mollié model in stan. Spatial and Spati-temporal Epidemiology. 31.
最近の調査データ分析でよく使われるみんな大好きミスターPことMRPについて調べていると、空間構造を入れるという話があって、そりゃまああるだろうとは思うのだが、そこに出てくるBYM2モデルというのになじみがなくて面食らった。調べてみたら、Gelmanさんたちが「StanでBYM2モデルをどう書くか」というチュートリアル論文を出しておられたので、パラパラめくってみた。
本題は3節のStanコードの紹介なのだが、単にBYM2モデルってどんなんなのかを知りたいだけないので、申し訳ないけど3節以降はパス。
1. イントロダクション
[略]
2. モデル
2.1 条件つき自己回帰モデル
地域数を\(N\)とする。条件付き自己回帰モデル(CARモデル)では、正規確率変数 \(\phi = (\phi_1, \ldots, \phi_N)^\top\)があって、$$ \phi_i | \phi_j, j \neq i \sim N\left( \sum_{j=1}^N w_{ij} \phi_j, \sigma^2 \right)$$ と考える。同時分布は$$ \phi \sim N(0, Q^{-1}) $$ という形になる。いま、対角行列\(D\)の対角要素には地域\(i\)の隣接地域数をいれ、隣接行列\(A\)の非対角要素は隣接の時に1とする。すると$$ Q = D(I – \alpha A) $$ となる。\(\alpha\)がCARモデルのパラメータである。[CARモデルって、同時分布は常にMVNなんだっけ？？ なにか隣接行列に制約が必要だったような気もするんだけど… まあいいや]
\(\phi\) の対数確率密度は、近隣グラフの要素数を\(n\)として $$ \frac{n}{2} \log(det(Q)) – \frac{1}{2} \phi^\top Q \phi$$ に比例する。[??? \(N\)と\(n\)って同じ？ よくわからん] この計算は\(N\)が大きいときに大変になるので、MCMCで新しい提案のたびに\(\phi\)の確率密度を求めるのは大変すぎる。
2.2 Intrinsic CARモデル
CARモデルの\(\alpha\)を1にする。\(Q = D-A\)、つまり、対角要素は隣接地域数、非対角要素は隣接のとき-1となる。行列式は0になり、事前分布は非正則だが、事後分布は正則である。
対数確率密度の計算は楽になるけど、MCMCに数日かかっていたのが数時間になります、という感じ。
同時分布はこういう風にかける。\(i, j\)が隣接していることを\(i \sim j\)と書く。$$ p(\phi) \propto \exp \left( -\frac{1}{2} \sum_{i \sim j} (\phi_i – \phi_j)^2 \right) $$ つまり、隣接ペアの差が罰則項になっているわけだ。
2.3 Basag-York-Mollieモデル
BYMモデルは対数正規ポアソンモデルで、疾病リスク・マッピングのために開発された。空間的スムージングのためのICARコンポーネントと、ふつうのランダム効果コンポーネントを含んでいる。
未知の相対リスク\(\eta_i\)が疾患事例数\(y_i\)で表現されているとしよう。このモデルでは、まず$$ \eta_i = \mu + x \beta + \phi + \theta$$ と考える。\(\phi\)がICAR空間コンポーネントである。
このモデルはフィッティングが難しい。Besag & Mollie (1991)は\(\phi, \theta\)の精度パラメータ(\tau_\phi, \tau_\theta\)にガンマ事前分布を与えた。その後、どっちも同じくらい強調するようなフェアな事前分布をどう指定するかということがいろいろ提案されて… [略]
2.4 BYM2モデル
[本節、わかりにくいので逐語訳]
BYM2モデル(Reibler et al., 2016)は罰則つき複雑性というフレームワーク(Simpson et al., 2017)に従うモデルである。このフレームワークでは、パラメータにsensibleなハイバーパラメータを割り当てられるようにするため、パラメータが明確な解釈を持つようなモデルが支持される。BYM2モデルはBYMモデルと同様に、空間的誤差項と非空間的誤差項を持ち、Leroux et al.(2000)のモデル同様に、結合されたコンポーネントについての単一の精度(スケール)パラメータ\(\sigma\)と、空間的変動と非空間的変動の量の混合パラメータ\(\rho\)を持つ。\(\sigma\)が結合コンポーネントのSDであることを正当的にするためには、それぞれの\(i\)について\(Var(\phi_i) \approx Var(\theta_i) \approx 1\)とすることが重要である。そのためにスケーリングファクター\(s\)をモデルに追加する。\(s\)は分散の割合を\(\rho\)にスケールする。[なんのこっちゃ…]
スケーリングファクター\(s\)はデータに依存するので、それはモデルにデータとして入る。Rieblerらは、分散の幾何平均が1になるようにスケーリングすることを推奨している。このスケーリングファクターは近隣グラフから算出できる。Stanプログラムではtransformed dataブロックで求められるが、RのINLA::inla.scale.model関数で求めてStanにデータとして渡すことにする。
BYM2モデルでは、BYMモデルにおける\(\phi + \theta\)を次のように書き換える。$$ \left( (\sqrt{\rho/s}) \phi^* + (\sqrt{1-\rho}) \theta^* \right)\sigma$$ \(\phi*\)はICARモデル。\(\theta^*\)は、結合しているサブグラフの数を\(n\)として\(\theta^* \sim N(0, n)\)である。\(\rho \in [0,1]\)は分散のうちどれだけが空間的に相関する誤差項から来ていてどれだけが独立した誤差項から来ているかをモデル化している。\(s\)はスケーリングファクターで、近隣グラフに基づいて\(Var(\phi_i) \approx 1\)になるように算出される。最後に、\(\sigma \geq 0\)は結合された誤差項の全体的なSDである。
近隣グラフが完全結合でない (\(n \gt 1\))ようなBYM2モデルでは、それぞれの結合したサブグラフはそれ自身の分散を持っており、従ってスケーリングする必要がある。Stan言語は十分にパワフルなので、つながっていないサブグラフや島領域があっても大丈夫だが、それぞれのサブグラフを追跡する必要があるのでコードは複雑になる。
3. Stanプログラム
[パス]
4. 事例: NYCにおける若者の歩道でのケガ、2005-2014
[パス]
————-
ふうん。要するにBYM2モデルというのは、地域のインシデント数をアウトカムに取ったポアソンモデルで、誤差項として普通の誤差項とICAR空間誤差項の両方が入っているような奴らしい。へええ。そんなんよう識別しはりますねえ。
細かいところはよく理解できてないんだけど、ま、この話はこのくらいでいいや…