10 階層ベイズモデル -GLMMのベイズモデル化-

10.1 例題：個体差と生存種子数（個体差あり）

架空植物100個から種子8個ずつを採取．
個体$i$で生存している種子数$y_i$．

図10.1(B)：生存種子数$y_i$のヒストグラム．
二項分布と全然合わない．個体に由来するランダム効果により過分散が生じている．

二項分布で説明しようとすると，生存確率の最尤推定値は0.504．
$N=8$なので二項分布の分散は2.00くらいになると期待されるが，実際の観測データの分散は9.93．

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

生存確率が全個体で共通という仮定では説明できない．
過分散が生じる要因はこのデータだけでは特定できない．

個体ごとに由来する原因不明な差異などを組み込んだGLMがGLMM．
しかし，個体差だけでなく場所差なども考慮するGLMMではパラメータの最尤推定が困難．
$\to$ GLMMのベイズモデル化．

リンク関数と線形予測子：${\rm logit}(q_i) = \beta + r_i$．
$\beta$：全個体に共通するパラメータ
$r_i$：個体差．平均ゼロで標準偏差$s$の正規分布にしたがうと仮定．

データが得られる確率は，全個体100個の二項分布の積

p(Y | \beta, \{r_i\}) 
= \prod_i = \begin{pmatrix}
8 \\ y_i
\end{pmatrix}
q^{y_i} (1-q)^{8-y_i}

推定したい事後分布は$\propto p(Y | \beta, {r_i}) \times$事前分布なので，あとは事前分布を指定すれば良い．

$\beta$に事前の制約はないので無情報事前分布を指定．平均ゼロで標準偏差100のひらべったい正規分布．

$r_i$の事前分布は上で記したように正規分布

p(r_i|s)
= \frac{1}{\sqrt{2\pi s^2}} \exp \left(- \frac{r_i^2}{2 s^2}\right)

とする．
ベイズ統計モデルでは，$s$の事後分布を推定する．
$s$は正の値であれば何でもよいので，事前分布$p(s)$を無情報事前分布としてよい．ここでは$0$から$10^4$までの連続一様分布とする．

このように，個体差$r_i$の事前分布$p(r_i|s)$を決める$s$についても事前分布$p(s)$が設定されているとき（事前分布のパラメータにさらに事前分布が設定されているとき），$p(r_i|s)$を階層事前分布と呼ぶ．
$s$：hyper parameter
$p(s)$：hyper prior
階層事前分布を使っているベイズ統計モデルが階層ベイズモデル．

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

この例題の階層ベイズモデルの事後分布は

p(\beta,s,\{r_i\}|Y)
\propto p(Y|\beta,\{r_i\}) p(\beta) p(s) \prod_i p(r_i|s)

ここで

\begin{align}
p(Y|\beta,\{r_i\}) p(\beta) p(\{r_i\})
&= p(Y|\beta,r_1,r_2,...,r_n) p(\beta) p(r_1,r_2,...,r_n|s) p(s) \\
&= p(Y|\beta,r_1,r_2,...,r_n) p(\beta) p(r_1|s) p(r_2|s) \cdots p(r_n|s) p(s) \\
&= p(Y|\beta,\{r_i\}) p(\beta) \prod_i p(r_i|s) p(s)
\end{align}

10.3.1 階層ベイズモデルのMCMCサンプリング

例題のモデルをWinBUGSで実装．

10.3.2 階層ベイズモデルの事後分布推定と予測

MCMCサンプリングの結果．図10.3：パラメータごとの周辺事後分布．

図10.4：推定された事後分布を組み合わせて，生存種子数のヒストグラムを予測．
生存種子数$y$の確率分布は，二項分布$p(y|\beta,r)$と正規分布$p(r|s)$の無限混合分布．

p(y|\beta,s)
= \int^\infty_{-\infty} p(y|\beta,r) p(r|s) dr

$r$についての積分：事後分布$p(r|s)$にしたがう植物を無限個集めて平均を評価している．図7.8が参考になる．

${\beta,s}$については，すべてのMCMCサンプルごとに$p(y|\beta,s)$を評価して，$y$ごとにパーセンタイル点を表示．

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

図10.5：よく使われる3種類の事前分布

(A) 主観的な事前分布
多くのベイズ統計モデルは(B)と(C)で記述できるので，(A)が必要とされることは多くない．

(B) 無情報事前分布

(C) 階層事前分布¹

あるパラメータの事前分布の選択は，そのパラメータがデータ全体のどの範囲を説明しているかに依存する（表10.1）．

この章の例題では...
$\beta$：データ全体を説明する大域的なパラメータ $\to$ 無情報事前分布．
${r_i}$：似たような，ある範囲に分布するパラメータの集まり．個々の$r_i$はデータ全体のごく一部を説明． $\to$ 階層事前分布．

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

図10.6：一部に施肥処理をした場合．
植木鉢が10個あって，各植木鉢に10個体．植木鉢A-Eは無処理，植木鉢F-Jで施肥処理．

図10.7(A)：個体ごとの種子数$y_i$．処理ごとの平均も図示．
無処理では標本平均8くらい．$y_i$がポアソン分布にしたがうなら，標本標準偏差は$\sqrt{8}=2.83$くらいであるが，実際は過分散が起きてそう．

図10.7(B)：植木鉢ごとの種子数を箱ひげ図で示したもの．
過分散の原因は個体差だけでなく植木鉢の差にもありそう．
このデータだけを見ていても理由は不明．植木鉢を置いた環境や生育期間中の出来事に影響？

個体も植木鉢も反復ではなく疑似反復．
個体差と植木鉢差を同時に扱う統計モデルが必要．

個体$i$の種子数$y_i$のばらつき

p(y_i | \lambda_i)
= \frac{\lambda_i^{y_i} \exp (-\lambda_i)}{y_i !}

線形予測子と対数リンク関数を使って，平均種子数は

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

$\beta_1$：切片
$f_i$：施肥処理の有無
$\beta_2$：$f_i$の係数
$r_i$：個体$i$の効果（個体差を表す．100個）．
$r_{j(i)}$：植木鉢の効果（10個）．

事前分布は以下のようにする．
$\beta_1$と$\beta_2$は大域的な平均パラメータなので無情報事前分布．平均ゼロのひらべったい正規分布．
$s$と$s_p$は大域的なばらつきパラメータ．$0$から$10^4$の範囲をとる一様分布．
$r_i$と$r_{j(i)}$は局所的パラメータ．階層事前分布．平均ゼロで標準偏差それぞれ$s$と$s_p$．

統計モデルをWinBUGSで実装．

これにより事後分布推定から，施肥処理の効果なさそうとわかる．しかし，統計モデリングで手抜きをすると，あるはずのない「処理の効果」が「推定」されてしまうことに注意．

10.6 この章のまとめと参考文献

(A)だとひとつの確率分布で$r_i$を表そうとしている．(C)だと$r_i$は正規分布で表せるとして，その分散も確率分布を持っていると考えることで汎用性を高めている，みたいなイメージ？ ↩