階層ベイズモデルについての読書メモ

岩波データサイエンス vol.1の読書メモ

階層ベイズ最初の一歩

階層ベイズモデル(階層事前分布を使った統計モデル)からベイズ統計モデルの良さを解説

階層ベイズモデルは似たようなパラメータに制約を与える特徴がある
　→ それによりパラメータが多い場合であっても統計モデルをデータにうまく当てはめられる

例題

12才の男子学生(のべ1082人)の身長を測定し、1年後に2回目の測定をしたデータを使う。最初に10調査県で50人ぐらいの男子学生の身長を測定、それからランダムに選ばれた5県に対して普通の給食(甲タイプ)とは異なる乙タイプの給食を食べさせる。１年後に身長を測定し乙タイプの給食が身長の伸びに与える効果を調査する。 正解として乙タイプ給食の効果はゼロで設定してテストデータを生成した。つまり身長増減と給食タイプは無関係ということが導出できれば良い。

まず"直線当てはめ"を使ってみる

まず直線当てはめができる条件として以下を満たしていなければいけない

それぞれ以下を表す Y{i,j}: 測定i回目における県jの平均値 Y{i,j} ~ N(μ{i,j}, σ²): 測定i回目における県jの平均値は平均μ{i,j}で標準偏差がσの正規分布に従う

測定i回目における県jでの身長の平均μ{i,j}は以下になっているものとする。
(測定1の甲・乙)　μ{1,j}　= β1 (β1は１回目の測定の平均)
(測定2の甲)　μ{2,j}　= β1 + β2 (β2は身長の変化量)
(測定2の乙)　μ{2,j}　= β1 + β2 + β3(β3は乙タイプとの比較量)
以下の説明変数を用意する A_i: 第i回目の測定 X_j: 県毎の異なる給食タイプ(甲タイプが0で乙タイプが1)

この各項についている(β1, β2, β3)を係数と呼ぶ

μ_{i,j}をデータに当てはめる

動作確認用のデータ(県、標本サイズ、平均身長、身長の標準偏差、測定回、給食タイプが含まれる)をダウンロードしてきてRのglm関数を使ってμ{i,j}のモデルへの当てはまりの良さ(尤度)を最大化するように(β1, β2, β3)を計算してみる。 そしたらβ2 = 5.647, β3 = -1.720(甲タイプの身長増加の平均は5.65cm、乙タイプの身長増加の平均は5.65 - 1.72 = 3.93cm)となった。この時p=0.045という値で有意水準5%に収まっているので乙タイプの給食では身長の伸びが悪くなるという結果が出てしまった。

モデルの不備を見直す

μ_{i,j}へのデータの当てはめを10⁵回繰り返してみる、p < 0.05が成り立つため5000回未満になるはずだがシュミレーションでは10000回以上 p < 0.05となった。このように結果が一致していなのはモデルに不備があるため。

(階層ではない)ベイズモデルにしてみる

以下のシンプルなモデルをベイズモデルにするためにはパラメータ(β1, β2, β_3)の事前分布の指定が必要になる。

事前分布とは
→ パラメータ毎に推定したい範囲を確率分布として指定したもの。どんな値でも良いとしてする場合はむ情報事前分布と呼ばれる。データのばらつきを表すσの事前分布を検討する。

以下の事前分布でベイズモデルに当てはめる - 係数{βk}と標準偏差σある値をとる時、身長平均がY_i,jとなる確率は正規分布N(μ{i,j}, σ²)から求めれる - 係数{βk}は無情報事前分布に従う - 標準偏差σは無情報事前分布に従う これで計算したらβ3は95%の区間(95%ベイズ信頼区間)が-2.411 ~ -0.0065となった

階層ベイズモデルにしてみる

より良い統計モデルを目指してみる。平均ち{Y_i,j}の不確かさ(標準誤差)を(標準偏差) / √{標本サイズ}　で近似してみる。測定回i,県j毎の標準誤差をS{i,j}とすることで身長の平均値Y{i,j}を以下のように表せられる。

それから測定i回目における県jでの身長の平均μ{i,j}を決める線形予測子の構成を考える
(測定1の甲・乙)　μ{1,j}　= β1 + r{1,j}
(測定2の甲)　μ{2,j}　= β1 + r{1,j} + β2 + r{2,j}
(測定2の乙)　μ{2,j}　= β1 + r{1,j} + β2 + r{2,j} + β3
r{1,j}とr{2,j}はそれぞれ"1回目の測定における県毎の身長差"、"2回目の測定における県毎の身長増加の差"を表している。これらを合わせると以下の式で表せられる。

事前分布について 係数{b_k}についてはベイズモデルと同様に無情報事前分布を指定 {r_i,j}の事前分布としてはr_i,j ~ N(0, σ²_i)を指定する(県毎の差(r_i,j)は平均0, 標準偏差σiの正規分布に従う)

これで計算したら95%ベイズ信頼区間が0に入っている →　甲タイプと乙タイプで差が特にないと言える

階層ベイズモデルの利点

県毎の差を表すパラメータを扱えるようになった。それにより県毎の影響を受けない正確な推定が行えるようになっている。