Schmid, V.J., Held, J. (2007) Bayesian Age-Period-Cohort Modeling and Prediction – BAMP. Journal of Statistical Software, 21(8).
ベイジアン・コウホート分析のソフトウェアBAMPの解説書。この段階ではスタンドアローンのソフトだが、その後Rパッケージに移植されている。仕事の都合で、実戦投入前の儀式として読んだ。
年齢集団のインデクスを\(i = 1, \ldots, I\), 時代のインデクスを\(j = 1, \ldots, J\), 生年コホートのインデクスを\(k = 1, \ldots, K\)とする。年齢と時代が同じ時間尺度上の値ならば \(k = k(i,j) = (I-i) + j\)である。
APCモデルの古典的文献では、APCモデルは対数線形ポワソンモデルとして捉えられるか、二項ロジットモデルとして捉えられることが多かった。後者の場合、生起数\(y_{ij}\)について$$ y_{ij} \sim B(n_{ij}, p_{ij})$$ $$ \eta_{ij} = \log \left( \frac{p_{ij}}{1 – p_{ij}} \right) = \mu + \theta_i + \phi_j + \psi_k$$ $$ \sum_i \theta_i = \sum_j \phi_j = \sum_k \psi_k = 0$$ このモデルは識別できない。任意の\(c\)について$$ \mu \rightarrow \mu – c I$$ $$ \theta_i \rightarrow \theta_i + ci$$ $$ \phi_j \rightarrow \phi_j – cj$$ $$ \psi_k \rightarrow \psi_k + ck$$ と線形変換しても\(\eta_{ij}\)は不変だからだ。従ってAPC効果の推定値は非線形トレンドに限定して解釈しなければならない。そのためBAMPは二階差分を提供する。[おっとぉ… これこれこういう制約をかけてうまいこと識別するから解釈してくれてよろしくてよ、と主張するわけではないのね]
BAMPはAPCモデルに対してベイジアン階層アプローチを用いる。パラメータ\(\theta, \phi, \psi\)についていろんな次数のランダムウォーク(RW)事前分布を用いるのである。RW1事前分布は定数トレンド、RW2事前分布は線形トレンドという仮定に基づいている。RW1事前分布は効果の一回差分を確率的に0に制約する(識別可能になる)。RW2事前分布は二階差分を確率的に0に制約する(識別は難しい)。
以下では\(\phi\)のRW事前分布について説明するが、\(\theta, \psi\)についても同様である。
精度パラメータを\(\kappa^{-1}\)とする。RW1の場合は$$ p(\phi_1) \propto const$$ $$ \phi_j | \phi_{j-1} \sim N(\phi_{j-1}, \kappa^{-1}) \ \ for \ \ j=2,\ldots, J$$ RW2の場合は$$ p(\phi_1) = p(\phi_2) \propto const$$ $$ \phi_j | \phi_{j-1}, \phi_{j-2} \sim N(2\phi_{j-1}-\phi_{j-2}, \kappa^{-1}) \ \ for \ \ j=3,\ldots, J$$ まあいずれにせよ、同時分布\(\Phi = (\phi_1, \ldots, \phi_J)’\)は$$ p(\Phi) \propto \lambda^{rg(R)/2} \exp \left( – \frac{\lambda}{2} \Phi’ R \Phi \right) $$ となる。[ここは写経したけど、学力不足の哀しさで、なにがなんだかわからない… \(\lambda^{rg(R)/2}\)ってなによ…]
BAMPは\(\kappa\)の事前分布をガンマ分布\(Ga(a,b)\)とする。\(a=1, b=0.001\)とか\(a=b=0.001\)とかを用いる。切片\(\mu\)の事前分布は\(p(\mu) \propto const\)とする。
未来についてprojectionする場合、RW1なら $$ \phi_{J+1} \sim N(\phi_J, \kappa^{-1}) $$ RW2なら $$ \phi_{J+1} \sim N(2\phi_J – \phi_{J-1}, \kappa^{-1}) $$ となる。モデル上は時代、年代、コホートのすべてについて補外が可能である(年代の補外はいらんだろうけど)。
\(n_{i,J+1}\)を決めてくれれば\(p_{i, J+1}\)もprojectionできる。$$ \eta_{i,J+1} = \mu + \theta_i + \phi_{J+1} + \psi_{k(i,J+1)}$$ $$ p_{i,J+1} = \frac{ \exp(\eta_{i,J+1})}{1+\exp(\eta_{i,J+1})} $$ とするわけね。
BAMPは事後分布をMHアルゴリズムでサンプリングする。[…中略…]
さて、APC効果で説明できないような異質性を説明するために $$ \eta_{ij} + \mu + \theta_i + \phi_j + \psi_k + z_{ij} $$ $$ z_{ij} \sim N(0, \delta^{-1}) $$ と、\(z_{ij}\)という項をいれましょうという提案もある。BAMPはこれにも対応してます。\(\delta\)の事前分布は\(Ga(a_\delta, b_\delta)\)とする。
このモデルの良いところは、APC効果の条件付き確率分布が正規分布になるという点である。
外れ値や効果の異質性のせいで効果の推定値の分散が大きくなってしまうことがある。こういうときは時間尺度に異質性を入れると良い。たとえば、すべての\(j\)について$$ \phi_j \sim (\phi_j^*, \lambda_1)$$ とし、\(\Phi^* = (\phi_1^*, \ldots, \phi_J^*)’\)の事前分布を精度\(\lambda_2\)のRWとする。\(\Phi\)と\(\Phi^*\)を縦に積んだベクトル\(\tilde{\Phi}\)を考える。その長さは\(2J\)で、\(2j-1\)番目の要素が\(\phi_j\)、\(2j\)番目の要素が\(\phi^*\)である。すると\(\tilde{\Phi}\)の事前分布は…[式は省略するけど、意外にシンプルな式になる]。
こういう修正はコホートや年齢についてもできる(年齢についてやろうと思うことはないだろうけど)。
[ふうん。要するに、状態空間モデルを考えて、\(\eta_{ij}\)のモデルは観測方程式のほうで書くってことですかね]
[スタンドアローン版BAMPの使い方。スキップ]
というわけで、ベイジアンAPCモデルは年齢・時代で層別された死力や生起率の分析に際して便利なツールである。ただし線形トレンドの解釈の際にはそのモデルがかけている確率的制約についてよく考えること。APCモデルは原因について知るための最初のステップに過ぎないのだ。
云々。
——–
年代・時代・コホートパラメータの解釈可能性について、かなり弱気というか、謙虚なのが印象に残った。いわゆるコウホート分析は、そもそもパラメータを識別できっこないという弱みを抱えているわけだが、この解説書と異なり、もっと実質的な知識を制約として導入して識別まで持っていくアプローチだと、解釈可能性についてももうちょっと強気な態度になるのかもしれませんね。