【10, 11章】データ解析のための統計モデリング入門実装

久保 拓弥先生が執筆した「データ解析のための統計モデリング入門」において、比較的難易度の高い最後の章を実装含めてより詳しくまとめたスライド

Transcript

1. 不確実性を考慮した階層ベイズに よる一般化線形混合モデリング yusumi
2. 目次 2 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ
3. 目次 3 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ
4. 背景 線形モデルの発展過程 4 最小二乗法 線形モデル 最尤推定法 一般化線形モデル 一般化線形混合モデル 階層ベイズモデル 正規分布以外の 確率分布を扱いたい 個体差・場所差といった ランダム効果を扱いたい もっと自由で 現実的な 統計モデルを！ MCMCによる事後分布の推定 GLMM GLM
5. 背景 生存種子数を推定する実験 5 植物個体 𝑖 種子数 𝑛𝑖 = 8 生存数 𝑦𝑖 = 3 100個体分実験 生存数 個 体 数 個体差あり 種子 抽出 調査 この分布を推定できる 統計モデルが欲しい 生存数の観測分布 可視化
6. 観測分布を推定できる統計モデルの検討 背景 6 植物個体 𝑖 100個体分実験 生存数の観測分布 生存数 個 体 数 個体差あり データ構造 8個の種子のうち𝑦𝑖 個が生存 二項分布の構造 生存種子数の ヒストグラム作成
7. 背景 二項分布の確率分布 7 𝑝 𝑦𝑖 |𝑞𝑖 = 8 𝑦𝑖 𝑞𝑖 𝑦𝑖 1 − 𝑞𝑖 8−𝑦𝑖 ~ Bin 𝑦𝑖 |8, 𝑞𝑖 ‐ 𝑞𝑖 : 個体 𝑖 の生存種子確率 ‐ 𝑦𝑖 : 個体 𝑖 の生存種子数 種子数 𝑛𝑖 = 8 生存数 𝑦𝑖 = 3 調査 8個の種子のうち𝑦𝑖 個が生存する確率 ※ i. i. d : independent and identically distributed : 独立同分布 i. i. d 独立同分布に従うと仮定
8. 背景 ロジスティック関数で生存確率を表現 8 𝑞𝑖 = 1 1 + exp −𝛽0 − 𝛽1 𝑥1 − ⋯ ↔ log 𝑞𝑖 1 − 𝑞𝑖 = 𝛽0 + 𝛽1 𝑥1 + ⋯ ロジット関数 logit 𝑞𝑖 線形予測子 𝑧𝑖 ‐ リンク関数に変換 𝑧 = 𝛽0 + 𝛽1 𝑥のとき
9. 背景 一般化線形モデル (GLM) による推定 9 logit 𝑞𝑖 = 𝛽 植物個体 𝑖 100個体分実験 個体差あり ‐ 個体 𝑖 の生存確率 𝑞𝑖 の設計 ロジット関数 全個体共通の パラメータ (切片) 個体差の影響無視 植物個体差の影響は 考えずに生存種子数 を推定しよう GLMの立場 ※ GLM : Generalized Linear Model
10. 背景 最尤推定からパラメータ 𝑞 を推定 10 𝐿 𝑞 = ෑ 𝑖=1 100 8 𝑦𝑖 𝑞𝑦𝑖 1 − 𝑞 8−𝑦𝑖 log 𝐿 𝑞 = ෍ 𝑖=1 100 log 8 𝑦𝑖 + 𝑦𝑖 log 𝑞 + 8 − 𝑦𝑖 log 1 − 𝑞 log 𝐿 𝑞 を最大にする (Appendix1) 𝑞 = 1 8 1 100 ෍ 𝑖=1 100 𝑦𝑖 (𝑞1 = 𝑞2 = ⋯ = 𝑞100 = 𝑞) ‐ 𝑞𝑖 : 個体 𝑖 の生存種子確率 ‐ 𝑦𝑖 : 個体 𝑖 の生存種子数
11. 背景 最尤法による生存種子数推定結果 11 最尤法により得られた二項分布 植物個体 𝑖 100個体分実験 個体差あり 個体差の影響無視 生存種子数分の ヒストグラム作成 モデルの当てはまりが悪い GLMによる単純なパラメータ推定ではダメ
12. 背景 個体差無しの仮定では説明できないタスク 12 最小二乗法 線形モデル 最尤推定法 一般化線形モデル 一般化線形混合モデル 階層ベイズモデル 正規分布以外の 確率分布を扱いたい 個体差・場所差といった ランダム効果を扱いたい もっと自由で 現実的な 統計モデルを！ MCMCによる事後分布の推定 GLMM GLM GLMでは上手く 説明できないデータ
13. 目的 個体差ありデータへのアプローチ 13 最小二乗法 線形モデル 最尤推定法 一般化線形モデル 一般化線形混合モデル 階層ベイズモデル 正規分布以外の 確率分布を扱いたい 個体差・場所差といった ランダム効果を扱いたい もっと自由で 現実的な 統計モデルを！ MCMCによる事後分布の推定 GLMM GLM 次に検討するモデル (今回のテーマ)
14. 目的 実験結果からモデルを改善 14 今回検討するモデル ‐ 個体差を考慮していない ‐ データのばらつきに頑丈でない ‐ 個体差も考慮 ‐ データの不確実性を加味 生存種子数 個体 𝑖 点推定 生存種子数 個体 𝑖 分布推定 実験で用いたモデル
15. 目次 15 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ
16. GLMMの導入 + ベイズモデル化 植物実験における課題 16 植物個体 𝑖 種子数 𝑛𝑖 = 8 生存数 𝑦𝑖 = 3 100個体分実験 生存数 個 体 数 個体差あり 種子 抽出 調査 生存数の観測分布 可視化 個体差を含めた生存種子数 の分布を推定したい
17. GLMMの導入 + ベイズモデル化 ベースモデルのベイズ化 17 ‐ 最尤推定 (点推定) はデータのばらつきに対して頑丈でない 不確実性を考慮 生存種子数 生存種子数 個体C 個体B 個体A 点推定 個体A 個体B 個体C 分布推定
18. GLMMの導入 + ベイズモデル化 一般化線形混合モデル (GLMM) による推定 18 logit 𝑞𝑖 = 𝛽 + 𝑟𝑖 植物個体 𝑖 100個体分実験 個体差あり ‐ 個体 𝑖 の生存確率 𝑞𝑖 の設計 ロジット関数 全個体共通の パラメータ (切片) 個体差の影響考慮 植物個体差の影響も 考えて生存種子数を 推定しよう GLMMの立場
19. GLMMの導入 + ベイズモデル化 GLMMの詳細 19 ロジット関数 logit 𝑞𝑖 = 𝛽 + 𝑟𝑖 全個体共通の パラメータ 個体 𝑖 の個体差 パラメータ 𝑟𝑖 ~𝑁 0, 𝑠2 0 𝑠 𝑠 新たに追加した パラメータ 個体 𝑖 の 生存確率 ※ GLMM : Generalized Linear Mixed Model
20. GLMMの導入 + ベイズモデル化 生存した種子が得られる確率 20 𝑝 𝒀|𝛽, 𝑟𝑖 = ෑ 𝑖=1 100 8 𝑦𝑖 𝑞𝑖 𝑦𝑖 1 − 𝑞𝑖 8−𝑦𝑖 尤度 二項分布 ‐ 個体数𝑁 = 100 ‐ 種子数 𝑛𝑖 = 8 ‐ 生存数 𝑦𝑖 ‐ 𝑟𝑖 = 𝑟1 , 𝑟2 , ⋯ , 𝑟100 ‐ 全100個体の二項分布の積
21. GLMMの導入 + ベイズモデル化 事後分布の推定 21 事後分布 ∝ 𝑝 𝒀|𝛽, 𝑟𝑖 ×事前分布 尤度 確率を更新 未定義
22. GLMMの導入 + ベイズモデル化 事前分布の設計 22 𝑝 𝛽 = 1 2𝜋 × 1002 exp −𝛽2 2 × 1002 ‐ 線形予測子の切片 𝛽 の事前分布 ‐ 個体差 𝑟𝑖 の事前分布 𝑝 𝑟𝑖 |𝑠 = 1 2𝜋𝑠2 exp −𝑟𝑖 2 2𝑠2 𝑟𝑖 ~𝑁 0, 𝑠2 0 𝑠 𝑠 𝛽~𝑁 0, 1002 平たい 𝑠の事前分布も 導入したい
23. GLMMの導入 + ベイズモデル化 超パラメータ 𝒔 の事前分布の設計 23 𝑝 𝑠 0から𝟏𝟎𝟒までの連続一様分布 0 𝟏𝟎𝟒 面積1 𝑟𝑖 ~𝑁 0, 𝑠2 𝑠の事前分布も導入したい 𝑠~𝑈 0, 104 事前分布の 事前分布を導入 𝑝 𝑟𝑖 |𝑠 : 階層事前分布 𝑝 𝑠 : 超事前分布
24. GLMMの導入 + ベイズモデル化 階層ベイズモデル 24 𝑠~Uniform 𝑟𝑖 ~Normal 𝛽~Normal 𝑦𝑖 ~Binomial 𝑞𝑖 ~Determistic 100 階層事前分布 超事前分布 観測変数 潜在変数 決定論的に 決まる変数 𝑝 𝛽, 𝑠, 𝑟𝑖 |𝒀 ∝ 𝑝 𝒀|𝛽, 𝑟𝑖 𝑝 𝛽 ෑ 𝑖=1 100 𝑝 𝑟𝑖 |𝑠 𝑝 𝑠 事後分布 尤度 事前分布 ‐ 個体数𝑁 = 100 ‐ 種子数 𝑛𝑖 = 8 ‐ 生存数 𝑦𝑖 ‐ 𝑟𝑖 = 𝑟1 , 𝑟2 , ⋯ , 𝑟100
25. GLMMの導入 + ベイズモデル化 補足1/2 : 事前分布の種類と選択 25 主観的な事前分布 無情報事前分布 階層事前分布 パラメータ 𝑥の事前 分布はこうなると思う パラメータ 𝑥の事前 分布は分からないわ パラメータ 𝑥の事前分布は パラメータ𝑦の事後分布 に支配されるよ
26. GLMMの導入 + ベイズモデル化 補足2/2 : 統計モデルのパラメータと事前分布 26 パラメータの種類 説明する範囲 パラメータの個数 採用すべき事前分布 全体に共通するばらつき 大域的 少数 無情報事前分布 個別のばらつき 局所的 多数 階層事前分布 生存種子数の実験の場合 ) ‐ 線形予測子𝛽 : 全個体共通のばらつき → 無情報事前分布が有効 ‐ 個体差 𝑟𝑖 : 個別のばらつき → 階層事前分布が有効
27. 目次 27 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ
28. 数値実験 : 個体差 MCMCサンプリングによる実験 28 ‐ Step数 ➢ Iterations : 1600回 ‐ 反復数 (サンプル列) ➢ Chains : 3回 ‐ 最初の何Stepの結果を使わないか ➢ Burn-in : 100回 ‐ 何個飛ばしてサンプリング結果を記録するか (可視化用) ➢ Skips : 3回 ‐ 正味のStep数 ➢ 1500回 ‐ 可視化 ➢ 500回 ※複数パラメータによる 事後分布の推定には ギブスサンプリングなどが 用いられる (Appendix 2)
29. 数値実験 : 個体差 事後分布推定結果 29 MCMCサンプリング 𝛽の事後分布 𝑠の事後分布 𝑟𝑖 の事後分布 𝑖 = 1~100 (500回) 曲線はサンプルされたヒストグラムに カーネル密度推定を施したもの
30. 数値実験 : 個体差 事後分布から個体毎の生存確率を予測 30 𝑟𝑖 の事後分布 𝑞𝑖 の分布 ‐ 𝑟𝑖 が小さい値をとりやすいと𝑞𝑖 の生存確率は小さくなる ※図は近似曲線だが実際はヒストグラム 𝑞𝑖 = 1 1 + exp −𝛽 − 𝑟𝑖 ‐ GLMM : logit 𝑞𝑖 = 𝛽 + 𝑟𝑖
31. 数値実験 : 個体差 事後分布から個体毎の生存確率を予測 31 𝑟𝑖 の事後分布 𝑞𝑖 の分布 ※図は近似曲線だが実際はヒストグラム ‐ 𝑟𝑖 が大きい値をとりやすいと𝑞𝑖 の生存確率は大きくなる 𝑞𝑖 = 1 1 + exp −𝛽 − 𝑟𝑖 ‐ GLMM : logit 𝑞𝑖 = 𝛽 + 𝑟𝑖
32. 数値実験 : 個体差 事後分布から個体毎の生存確率を予測 32 𝑟𝑖 の事後分布 𝑞𝑖 の分布 ※図は近似曲線だが実際はヒストグラム ‐ 𝑟𝑖 が0付近だと 𝑞𝑖 は0.5を中心とした分布となる 𝑞𝑖 = 1 1 + exp −𝛽 − 𝑟𝑖 ‐ GLMM : logit 𝑞𝑖 = 𝛽 + 𝑟𝑖
33. 数値実験 : 個体差 生存種子数の予測分布 33 𝑝 𝑦|𝛽, 𝑠 = න −∞ ∞ 𝑝 𝑦|𝛽, 𝑟 𝑝 𝑟|𝑠 𝑑𝑟 生存種子数 𝑦 の確率分布 二項分布 個体差の事後分布 ‐ 二項分布を事後分布で平均化
34. 数値実験 : 個体差 実装上の予測分布 34 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑠𝑖𝑗 = ෍ 𝑠=−10𝑠𝑖𝑗 𝑠=10𝑠𝑖𝑗 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑟𝑖𝑗 𝑝 𝑟𝑖𝑗 |𝑠 生存種子数が𝑦𝑖𝑗𝑘 となる確率 二項分布 超パラメータ 𝑠 の 事後分布を標準偏 差に持つ正規分布 ‐ MCMCの事後分布サンプルから計算 ‐ 𝑖 ∈ 1,2,3 : MCMCの反復 index ‐ 𝑗 ∈ 1,2, ⋯ , 1500 : MCMCの sampling index ‐ 𝑘 ∈ 0,1, ⋯ , 8 : 生存種子数のindex ‐ 𝑦𝑖𝑗𝑘 ∈ 0,1, ⋯ , 8 : 生存種子数が 𝑘 個 𝑝 𝑟𝑖𝑗 |𝑠 ~ 𝑁 0, 𝑠2
35. 数値実験 : 個体差 Step1. 𝑝 𝑟𝑖𝑗 |𝑠 の計算 35 ‐ 𝑠の事後分布からサンプル点を1つ取り出し計算 サンプル点 𝑝 𝑟𝑖𝑗 |𝑠 = 1 2𝜋𝑠2 exp −𝑟𝑖𝑗 2 2𝑠2 𝑠 Step3 Step4 Step2 Step1 ‐ 𝑖 ∈ 1,2,3 : MCMCの反復 index ‐ 𝑗 ∈ 1,2, ⋯ , 1500 : MCMCの sampling index ‐ 𝑘 ∈ 0,1, ⋯ , 8 : 生存種子数のindex ‐ 𝑦𝑖𝑗𝑘 ∈ 0,1, ⋯ , 8 : 生存種子数が 𝑘 個 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑠 = ෍ 𝑠=−10𝑠𝑖𝑗 𝑠=10𝑠𝑖𝑗 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑟𝑖𝑗 𝑝 𝑟𝑖𝑗 |𝑠
36. 数値実験 : 個体差 Step2. 生存確率の計算 36 Step1 Step4 Step3 ‐ 𝛽の事後分布からサンプル点を1つ取り出し生存確率 𝑞𝑖𝑗 を計算 サンプル点 𝛽𝑖𝑗 logit 𝑞𝑖𝑗 = 𝛽𝑖𝑗 + 𝑟𝑖𝑗 𝑞𝑖𝑗 = 1 1 + exp −𝛽𝑖𝑗 − 𝑟𝑖𝑗 GLMM 周辺化するパラメータ 𝑟𝑖𝑗 を定めた時の𝑞𝑖𝑗 の分布 Step2 ‐ 𝑖 ∈ 1,2,3 : MCMCの反復 index ‐ 𝑗 ∈ 1,2, ⋯ , 1500 : MCMCの sampling index ‐ 𝑘 ∈ 0,1, ⋯ , 8 : 生存種子数のindex ‐ 𝑦𝑖𝑗𝑘 ∈ 0,1, ⋯ , 8 : 生存種子数が 𝑘 個
37. 数値実験 : 個体差 Step3. 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑟𝑖𝑗 の計算 37 ‐ 二項分布から生存種子数が 𝑦𝑖𝑗𝑘 となる確率を計算 : 𝑘 ∈ 0,1, ⋯ , 8 Step4 Step2 Step1 Step3 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑟𝑖𝑗 = 8 𝑦𝑖𝑗𝑘 𝑞𝑖𝑗 𝑦𝑖𝑗𝑘 1 − 𝑞𝑖𝑗 8−𝑦𝑖𝑗𝑘 ‐ 𝑖 ∈ 1,2,3 : MCMCの反復 index ‐ 𝑗 ∈ 1,2, ⋯ , 1500 : MCMCの sampling index ‐ 𝑘 ∈ 0,1, ⋯ , 8 : 生存種子数のindex ‐ 𝑦𝑖𝑗𝑘 ∈ 0,1, ⋯ , 8 : 生存種子数が 𝑘 個 𝑞𝑖𝑗 = 1 1 + exp −𝛽𝑖𝑗 − 𝑟𝑖𝑗 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑠 = ෍ 𝑠=−10𝑠𝑖𝑗 𝑠=10𝑠𝑖𝑗 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑟𝑖𝑗 𝑝 𝑟𝑖𝑗 |𝑠
38. 数値実験 : 個体差 Step4. 予測分布から生存種子数を予測 38 ‐ 全ての 𝑖, 𝑗, 𝑘 について数値積分から予測分布 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑠 を計算 Step1 Step2 Step3 Step4 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑟𝑖𝑗 = 8 𝑦𝑖𝑗𝑘 𝑞𝑖𝑗 𝑦𝑖𝑗𝑘 1 − 𝑞𝑖𝑗 8−𝑦𝑖𝑗𝑘 𝑝 𝑟𝑖𝑗 |𝑠𝑖𝑗 = 1 2𝜋𝑠2 exp −𝑟𝑖𝑗 2 2𝑠2 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑠𝑖𝑗 = ෍ 𝑠=−10𝑠𝑖𝑗 𝑠=10𝑠𝑖𝑗 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑟𝑖𝑗 𝑝 𝑟𝑖𝑗 |𝑠
39. 数値実験 : 個体差 予測分布の別表現 39 ‐ 𝑘 についてまとめると 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑠𝑖𝑗 はカテゴリ分布となる 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑠𝑖𝑗 = ෑ 𝑙=0 8 𝑝 𝑦𝑖𝑗𝑙 |𝛽𝑖𝑗 , 𝑠𝑖𝑗 𝑚𝑙 ~ Cat 𝒎|𝒑 ‐ 𝑚𝑙 ∈ 0, 1 かつ σ𝑙=0 8 𝑚𝑙 = 1 : 1 of 𝐾 (= 9) 表現 ‐ 𝑚𝑙 = 1 if 𝑘 = 𝑙 ‐ 𝒎 = 𝑚0 , 𝑚1 , ⋯ , 𝑚8 ‐ 𝒑 = 𝑝 𝑦𝑖𝑗0 |𝛽𝑖𝑗 , 𝑠𝑖𝑗 , 𝑝 𝑦𝑖𝑗1 |𝛽𝑖𝑗 , 𝑠𝑖𝑗 , ⋯ , 𝑝 𝑦𝑖𝑗8 |𝛽𝑖𝑗 , 𝑠𝑖𝑗 ‐ 𝑖 ∈ 1,2,3 : MCMCの反復 index ‐ 𝑗 ∈ 1,2, ⋯ , 1500 : MCMCの sampling index ‐ 𝑘 ∈ 0,1, ⋯ , 8 : 生存種子数のindex ‐ 𝑦𝑖𝑗𝑘 ∈ 0,1, ⋯ , 8 : 生存種子数が 𝑘 個
40. 数値実験 : 個体差 得られた予測分布から生存種子をサンプリング 40 𝑦11𝑘 (𝑖 = 𝑗 = 1) の予測分布可視化 生存種子数のヒストグラム e.g. ) 100個体 サンプリング
41. 数値実験 : 個体差 全ての 𝑖, 𝑗 について生存種子数の分布推定 41 最尤推定で得られた生存種子数の分布 階層ベイズモデルから得られた生存種子数の分布
42. 目次 42 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ
43. 個体差 + 場所差の階層ベイズモデル 実験設定 43 Pot A Pot B Pot C Pot D Pot E Pot F Pot G Pot H Pot I Pot J 肥料なし : 𝑗 ∈ A, B, C, D, E 肥料あり : 𝑗 ∈ F, G, H, I, J 個体番号 𝑖 ∈ 1,2, ⋯ , 10 個体番号 𝑖 ∈ 11,12, ⋯ , 20 ⋯ ‐ 植物個体数 𝑁 = 100 ‐ 個体番号 𝑖 ∈ 1,2, ⋯ , 100 (100個) ‐ 植木鉢 j ∈ A, B, ⋯ , J (10個) 各個体の生存種子数 𝑦𝑖 を観測
44. 個体差+場所差の階層ベイズモデル 観測された種子数の Swarmplot と Boxplot 44
45. 個体差+場所差の階層ベイズモデル 観測された種子数の Swarmplot と Boxplot 45 肥料なし 植木鉢内の個体差に散らばりがある 肥料あり
46. 個体差+場所差の階層ベイズモデル 観測された種子数の Swarmplot と Boxplot 46 肥料なし 植木鉢内の個体差に散らばりがある 肥料あり 特定の植木鉢だけ種子数が少ない
47. 個体差+場所差の階層ベイズモデル 個体差と植木鉢差 (場所差) の発生 47 植木鉢差 個体差 ‐ 個体差 + 植木鉢差を考慮した種子数の予測モデルを考える
48. 個体差+場所差の階層ベイズモデル 種子数のばらつきの設定 48 𝑝 𝑦𝑖 |𝜆𝑖 = 𝜆𝑖 𝑦𝑖exp −𝜆𝑖 𝑦𝑖 ! ~ Poi 𝑦𝑖 |𝜆𝑖 ‐ 𝜆𝑖 : ポアソン分布の平均 (平均種子数) ‐ 𝑦𝑖 : 個体 𝑖 の種子数 ‐ ポアソン分布に従うと仮定 i. i. d
49. 個体差+場所差の階層ベイズモデル GLMMの設計 49 対数リンク関数 log 𝜆𝑖 = 𝛽1 + 𝛽2 𝑓2 + 𝑟𝑖 + 𝑟𝑗 𝑖 全個体共通の パラメータ 個体 𝑖 の個体差 𝑓2 = 0で肥料なし 𝑓2 = 1で肥料あり 個体 𝑖 の 平均種子数 施肥処理 の効果 植木鉢 𝑗 の効果 新たに追加した パラメータ
50. 個体差+場所差の階層ベイズモデル 事前分布の設計 50 𝛽1 , 𝛽2 は大域的なパラメータ → 無情報事前分布 𝑟𝑖 , 𝑟𝑗 𝑖 は局所的なパラメータ → 階層事前分布 𝑟𝑖 ~𝑁 0, 𝑠2 , 𝑟𝑗 𝑖 ~𝑁 0, 𝑠𝑝 2 𝛽1 , 𝛽2 ~𝑁 0, 1002 平たい 𝑠, 𝑠𝑝 は大域的なパラメータ → 無情報事前分布 𝑠, 𝑠𝑝 ~𝑈 0, 104 𝑟𝑖 , 𝑟𝑗 𝑖 の超パラメータ
51. 個体差+場所差の階層ベイズモデル 階層ベイズモデル 51 𝑠~Uniform 𝑟𝑖 ~Normal 𝛽1 ~Normal 𝑦𝑖 ~Poisson 𝜆𝑖 ~Determistic 100 階層事前分布 超事前分布 𝛽2 ~Normal 𝑠𝑝 ~Uniform 𝑟 𝑗(𝑖) ~Normal 10 𝑓2(𝑖) 観測変数 潜在変数 決定論的に 決まる変数 𝑝 𝛽1 , 𝛽2 , 𝑠, 𝑠𝑝 , 𝑟𝑖 , 𝑟 𝑗(𝑖) |𝒀 ∝ 𝑝 𝒀|𝛽1 , 𝛽2 , 𝑠, 𝑠𝑝 , 𝑟𝑖 , 𝑟 𝑗(𝑖) ෑ 𝑖=1 100 𝑝 𝑟𝑖 |𝑠 ෑ 𝑗=1 10 𝑝 𝑟 𝑗(𝑖) |𝑠𝑝 𝑝 𝛽1 𝑝 𝛽2 𝑝 𝑠 𝑝 𝑠𝑝
52. 目次 52 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ
53. 数値実験 : 個体差+場所差 MCMCサンプリングによる実験 53 ‐ Step数 ➢ Iterations : 1600回 ‐ 反復数 (サンプル列) ➢ Chains : 3回 ‐ 最初の何Stepの結果を使わないか ➢ Burn-in : 100回 ‐ 何個飛ばしてサンプリング結果を記録するか (可視化用) ➢ Skips : 3回 ‐ 正味のStep数 ➢ 1500回 ‐ 可視化 ➢ 500回
54. 数値実験 : 個体差+場所差 事後分布推定結果 1/2 54 𝛽1 の事後分布 𝛽2 の事後分布 𝑟𝑖 の事後分布
55. 数値実験 : 個体差+場所差 事後分布推定結果 2/2 + 𝜆𝑖 の分布 55 𝑟𝑗(𝑖) の 事後分布 𝑠 (青)と 𝑠𝑝 (オレンジ) の事後分布 𝜆𝑖 の分布
56. 数値実験 : 個体差+場所差 事後分布の推定から分かること 56 パラメータ 平均 標準偏差 下側2.5%点 上側97.5%点 𝛽1 1.35 0.520 0.358 2.49 𝛽2 -0.847 0.761 -2.24 0.714 𝑠 1.02 0.116 0.810 1.26 𝑠𝑝 1.04 0.374 0.438 1.78 ‐ パラメータ𝛽2 における事後分布の95%信用区間に0が含まれる 有意水準5%では施肥処理に効果があるとは言えない ※実験データは肥料の効果が0となるように作られているので検定結果は正しい
57. 目次 57 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ
58. 空間構造がある階層ベイズモデル 一次元空間上の個体数分布 58 個体数 𝑦𝑗 = 2 区間 𝑗 = 1 1 2 3 3 4 2 0 50 ‐ 50個の区間調査 : 𝑗 ∈ 1, 2, ⋯ , 50 ‐ 区間 𝑗 における個体数 : 𝑦𝑗 (架空データ) 一次元空間 ‐ 𝑦𝑗 の従う分布 : 𝑝 𝑦𝑗 |𝜆𝑗 = 𝜆𝑗 𝑦𝑗exp −𝜆𝑗 𝑦𝑗 ! ~ Poi 𝑦𝑗 |𝜆𝑗
59. 空間構造がある階層ベイズモデル 区間同士の隣接 59 個体数 𝑦𝑗 = 2 区間 𝑗 = 1 1 2 3 3 4 2 0 50 区間同士の幅を狭める
60. 空間構造がある階層ベイズモデル 区間相関の発生 60 区間相関 近い区間同士の 傾向は似ている 区間毎の独立性を仮定できない！ 𝑦𝑗 ~ Poi 𝑦𝑗 |𝜆𝑗 i. i. d
61. 空間構造がある階層ベイズモデル 観測個体数の分布 61 一次元空間上の架空データ 区間毎のポアソン分布の 真の平均 𝜆𝑗 観測数の変化は区間毎に滑らか 空間相関がある観測データ
62. 空間構造がある階層ベイズモデル GLMMの設計 62 対数リンク関数 log 𝜆𝑗 = 𝛽 + 𝑟𝑗 全個体共通の パラメータ 区間 𝑗 毎の場所差 区間 𝑗 の平均個体数 大域的 無情報事前分布 𝛽~𝑁 0, 1002 局所的 階層事前分布 𝑠~𝑈 0, 104
63. 空間構造がある階層ベイズモデル 近傍区間からの影響 63 注目区間 𝑦𝑗 = 3 𝑦𝑗−1 = 1 𝑦𝑗+1 = 2 近傍区間の設計 隣接区間のみの空間相関を考慮 区間 𝑗 𝑗 − 1 𝑗 + 1
64. 空間構造がある階層ベイズモデル 空間構造のある階層事前分布 64 𝑝 𝑟𝑗 |𝜇𝑗 , 𝑠 = 𝑛𝑗 2𝜋𝑠2 exp − 𝑟𝑗 − 𝜇𝑗 2 2𝑠2/𝑛𝑗 ~ 𝑁 𝜇𝑗 , 𝑠2 𝑛𝑗 ‐ 場所差 𝑟 𝑗 の事前分布 ‐ 𝜇𝑗 = 𝑟𝑗−1+𝑟𝑗+1 2 : 近傍区間から影響を考慮 ‐ 𝜇1 = 𝑟2 , 𝜇50 = 𝑟49 : 両端区間の処理 ‐ 𝑛𝑗 : 近傍区間の個数 (実験では2個)
65. 空間構造がある階層ベイズモデル 条件付き自己回帰モデル CAR 65 𝑝 𝑟𝑗 |𝑠 ∝ exp 1 2𝑠2 ෍ 𝑗~𝑗′ 𝑟𝑗 − 𝑟𝑗′ 2 ~ Gaussian 𝐂𝐀𝐑 ‐ 場所差 𝑟 𝑗 = 𝑟1 , 𝑟2 , ⋯ , 𝑟50 全体の (同時) 事前分布 ‐ 𝑗~𝑗′: 区間 𝑗 と近傍区間 𝑗′ の全ての組合せ 条件付き自己回帰 ※ CAR : Conditional Auto Regressive
66. 空間構造がある階層ベイズモデル 実装上のCARモデルの事前分布 66 𝑟 𝑗 ~ 𝑁 𝟎, 𝑠2 𝐷 − 𝛼𝑊 −1 ‐ Brook’s (1964) Lemma による設計方法を採用 ‐ 𝐷 = diag 𝑛𝑗 : 隣接区間の個数 𝑛𝑗 を成分に持つ50 × 50 対角行列 ‐ 𝛼 : 空間依存性を制御するパラメータ ‐ 𝑊 = 𝑤𝑗𝑗 = 0, 𝑤𝑗′𝑗 = 1 if 𝑗′ is neighbor of 𝑗, and 𝑤𝑗′𝑗 = 0 otherwise : 𝑤𝑖𝑗 を成分に持つ50 × 50 隣接行列 ‐ 𝛼 にも超事前分布を設計 : 𝛼 ~ 𝑈 0, 104
67. 空間構造がある階層ベイズモデル 階層ベイズモデル 67 𝑠~Uniform 𝛽~Normal 𝑦𝑗 ~Poisson 𝜆𝑗 ~Determistic 50 階層事前分布 超事前分布 観測変数 潜在変数 決定論的に 決まる変数 𝑟 𝑗 ~Gaussian CAR 𝑝 𝛽, 𝑠, 𝑟 𝑗 |𝒀 ∝ ෑ 𝑗=1 50 𝑝 𝑦𝑗 |𝜆𝑗 𝑝 𝑟 𝑗 |𝑠, 𝛼 𝑝 𝑠 𝑝 𝛽 𝑝 𝛼 事後分布 尤度 事前分布 𝛼~Uniform
68. 目次 68 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ
69. 数値実験 : 空間構造 MCMCサンプリングによる実験 69 ‐ Step数 ➢ Iterations : 60000回 ‐ 反復数 (サンプル列) ➢ Chains : 3回 ‐ 最初の何Stepの結果を使わないか ➢ Burn-in : 10000回 ‐ 何個飛ばしてサンプリング結果を記録するか (可視化用) ➢ Skips : 1回 ‐ 正味のStep数 ➢ 50000回 ‐ 可視化 ➢ 50000回
70. 数値実験 : 空間構造 事後分布および自己相関推定結果 70 MCMCサンプリング 事後分布 自己相関 実行時平均 ※ 𝑟 𝑗 の事後分布の可視化は割愛
71. 数値実験 : 空間構造 事後分布から個体数推定 71 𝜆𝑗 = exp 𝛽 + 𝑟𝑗 事後分布からサンプリングした個体数平均𝜆𝑗 の分布 ‐ 50000点のサンプリング Posterior sampling GLMM 空間相関を上手く表現できている
72. 数値実験 : 空間構造 データの一部を欠測させて再実験 72 一次元空間上の架空データ データの一部を欠測 欠測させたデータ
73. 数値実験 : 空間構造 欠測データに対する比較実験 空間相関無し 空間相関有り 73
74. 数値実験 : 空間構造 欠測データに対する比較実験 74 空間相関無し 空間相関有り
75. 数値実験 : 空間構造 階層ベイズモデル : 空間相関無し 75 𝑠~Uniform 𝛽~Normal 𝑦𝑗 ~Poisson 𝜆𝑗 ~Determistic 50 階層事前分布 超事前分布 観測変数 潜在変数 決定論的に 決まる変数 𝑟 𝑗 ~Normal 𝑝 𝛽, 𝑠, 𝑟 𝑗 |𝒀 ∝ ෑ 𝑗=1 50 𝑝 𝑦𝑗 |𝜆𝑗 ෑ 𝑖=1 50 𝑝 𝑟𝑖 |𝑠 𝑝 𝑠 𝑝 𝛽 事後分布 尤度 事前分布 𝑝 𝑟𝑗 |𝑠 ~ 𝑁 0, 𝑠2 i. i. d
76. 数値実験 : 空間構造 事後分布および自己相関推定結果 : 空間相関無し 76 MCMCサンプリング 事後分布 自己相関 実行時平均 ※ 𝑟 𝑗 の事後分布の可視化は割愛
77. 数値実験 : 空間構造 欠測データに対する比較実験 空間相関無し 空間相関有り 77
78. 数値実験 : 空間構造 階層ベイズモデル : 空間相関有り 78 𝑠~Uniform 𝛽~Normal 𝑦𝑗 ~Poisson 𝜆𝑗 ~Determistic 50 階層事前分布 超事前分布 観測変数 潜在変数 決定論的に 決まる変数 𝑟 𝑗 ~Gaussian CAR 𝑝 𝛽, 𝑠, 𝑟 𝑗 |𝒀 ∝ ෑ 𝑗=1 50 𝑝 𝑦𝑗 |𝜆𝑗 𝑝 𝑟 𝑗 |𝑠, 𝛼 𝑝 𝑠 𝑝 𝛽 𝑝 𝛼 事後分布 尤度 事前分布 𝛼~Uniform 𝑝 𝑟 𝑗 |𝑠 ~𝑁 𝜇𝑗 , 𝑠2 𝑛𝑗
79. 数値実験 : 空間構造 事後分布および自己相関推定結果 :空間相関有り 79 MCMCサンプリング 事後分布 自己相関 実行時平均 ※ 𝑟 𝑗 の事後分布の可視化は割愛
80. 数値実験 : 空間構造 欠測付き観測データの推定結果 80 空間相関無しの階層ベイズモデル 空間相関ありの階層ベイズモデル 欠測区間の分散が大きい 欠測区間の分散が小さい
81. 数値実験 : 空間構造 補足 : 二次元以上の空間構造 81 ‐ CARモデルの次元や隣接区間の個数は変化するが考え方は同じ
82. 目次 82 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ
83. まとめ 83 ‐ 階層ベイズモデルとは, 事前分布となる確率分布のパラメータにも事前分布を導入 したベイズ統計モデルである ‐ 無情報事前分布・階層事前分布による事前分布の設計方法を扱った ‐ 個体差+場所差などの複雑なデータ構造の統計モデリングでは, 階層ベイズモデル とMCMCサンプリングによるパラメータ推定がよい ‐ データの独立性を仮定できない時は, 空間相関を考慮したモデル設計が必要である
84. Appendix 1 二項分布に対する最尤法 84 log 𝐿 𝑞 = ෍ 𝑖=1 𝑛 log 𝑀 𝑦𝑖 + 𝑦𝑖 log 𝑞 + 𝑀 − 𝑦𝑖 log 1 − 𝑞 𝜕 𝜕𝑞 log 𝐿 𝑞 = 𝜕 𝜕𝑞 log 𝑞 ෍ 𝑖=1 𝑛 𝑦𝑖 + 𝜕 𝜕𝑞 log 1 − 𝑞 ෍ 𝑖=1 𝑛 𝑀 − 𝑦𝑖 = 1 𝑞 ෍ 𝑖=1 𝑛 𝑦𝑖 − 1 1 − 𝑞 ෍ 𝑖=1 𝑛 𝑀 − 𝑦𝑖 = 𝑛ത 𝑦 𝑞 − 𝑛𝑀 − 𝑛ത 𝑦 1 − 𝑞 = 0 ത 𝑦 = 1 𝑛 ෍ 𝑖=1 𝑛 𝑦𝑖 を解く. ∴ 𝑞 = ത 𝑦 𝑀 を得る. の最大化. 𝑞 𝐿 𝑞 最尤推定値
85. Appendix 2 実験環境 : 個体差および個体差+場所差 86 ‐ 言語・ライブラリ ➢ python ➢ pystan (確率的プログラミング言語Stanをpythonに移植したもの) ‐ Version ➢ Python : 3.6.6 ➢ pystan : 2.17.1.0 (最新の2.19ではコンパイルが通らない) ‐ OS : Windows10 64bit ‐ MCMC : HMC (Hamiltonian Monte Carlo​) ※ pystanの利用にはVC++コンパイラのインストールが必要
86. Appendix 3 実験環境 : 空間構造のMCMCサンプリング : 𝑟𝑗 87 ‐ 言語・ライブラリ ➢ Julia ➢ Mamba.jl ‐ Version ➢ Julia : 1.6.1 ➢ Mamba : 0.12.5 ‐ OS : Windows10 64bit ‐ MCMC : AMWG (Adaptive Metropolis within Gibbs​)
87. Appendix 4 実験環境 : 空間構造のMCMCサンプリング : 𝛼, 𝛽, 𝑠 88 ‐ MCMC : Slice (Shrinkage Slice​)
88. Appendix 5 実験環境 : 欠測付きMCMCサンプリング : 𝑦𝑗 89 ‐ MCMC : MISS (Missing Value Sampler​)
89. 参考文献 90 ‐ 久保拓弥 . “データ解析のための統計モデリング入門”. 岩波書店. ‐ 小森政嗣. "これからベイズ統計を使ってみたい人に." 日本音響学会誌 75.6 (2019): 351-357. ‐ 深澤圭太, and 角谷拓. "始めよう! ベイズ推定によるデータ解析." 日本生態学会誌 59.2 (2009): 167-170. ‐ Gelfand, Alan E., and Penelope Vounatsou. "Proper multivariate conditional autoregressive models for spatial data analysis." Biostatistics 4.1 (2003): 11-15.