Tweet
Tweet

Slide 1

Slide 1 text

不確実性を考慮した階層ベイズに よる一般化線形混合モデリング yusumi

Slide 2

Slide 2 text

目次 2 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ

Slide 3

Slide 3 text

目次 3 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ

Slide 4

Slide 4 text

背景 線形モデルの発展過程 4 最小二乗法 線形モデル 最尤推定法 一般化線形モデル 一般化線形混合モデル 階層ベイズモデル 正規分布以外の 確率分布を扱いたい 個体差・場所差といった ランダム効果を扱いたい もっと自由で 現実的な 統計モデルを！ MCMCによる事後分布の推定 GLMM GLM

Slide 5

Slide 5 text

背景 生存種子数を推定する実験 5 植物個体 𝑖 種子数 𝑛𝑖 = 8 生存数 𝑦𝑖 = 3 100個体分実験 生存数 個 体 数 個体差あり 種子 抽出 調査 この分布を推定できる 統計モデルが欲しい 生存数の観測分布 可視化

Slide 6

Slide 6 text

観測分布を推定できる統計モデルの検討 背景 6 植物個体 𝑖 100個体分実験 生存数の観測分布 生存数 個 体 数 個体差あり データ構造 8個の種子のうち𝑦𝑖 個が生存 二項分布の構造 生存種子数の ヒストグラム作成

Slide 7

Slide 7 text

背景 二項分布の確率分布 7 𝑝 𝑦𝑖 |𝑞𝑖 = 8 𝑦𝑖 𝑞𝑖 𝑦𝑖 1 − 𝑞𝑖 8−𝑦𝑖 ~ Bin 𝑦𝑖 |8, 𝑞𝑖 ‐ 𝑞𝑖 : 個体 𝑖 の生存種子確率 ‐ 𝑦𝑖 : 個体 𝑖 の生存種子数 種子数 𝑛𝑖 = 8 生存数 𝑦𝑖 = 3 調査 8個の種子のうち𝑦𝑖 個が生存する確率 ※ i. i. d : independent and identically distributed : 独立同分布 i. i. d 独立同分布に従うと仮定

Slide 8

Slide 8 text

背景 ロジスティック関数で生存確率を表現 8 𝑞𝑖 = 1 1 + exp −𝛽0 − 𝛽1 𝑥1 − ⋯ ↔ log 𝑞𝑖 1 − 𝑞𝑖 = 𝛽0 + 𝛽1 𝑥1 + ⋯ ロジット関数 logit 𝑞𝑖 線形予測子 𝑧𝑖 ‐ リンク関数に変換 𝑧 = 𝛽0 + 𝛽1 𝑥のとき

Slide 9

Slide 9 text

背景 一般化線形モデル (GLM) による推定 9 logit 𝑞𝑖 = 𝛽 植物個体 𝑖 100個体分実験 個体差あり ‐ 個体 𝑖 の生存確率 𝑞𝑖 の設計 ロジット関数 全個体共通の パラメータ (切片) 個体差の影響無視 植物個体差の影響は 考えずに生存種子数 を推定しよう GLMの立場 ※ GLM : Generalized Linear Model

Slide 10

Slide 10 text

背景 最尤推定からパラメータ 𝑞 を推定 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 = 𝑞) ‐ 𝑞𝑖 : 個体 𝑖 の生存種子確率 ‐ 𝑦𝑖 : 個体 𝑖 の生存種子数

Slide 11

Slide 11 text

背景 最尤法による生存種子数推定結果 11 最尤法により得られた二項分布 植物個体 𝑖 100個体分実験 個体差あり 個体差の影響無視 生存種子数分の ヒストグラム作成 モデルの当てはまりが悪い GLMによる単純なパラメータ推定ではダメ

Slide 12

Slide 12 text

背景 個体差無しの仮定では説明できないタスク 12 最小二乗法 線形モデル 最尤推定法 一般化線形モデル 一般化線形混合モデル 階層ベイズモデル 正規分布以外の 確率分布を扱いたい 個体差・場所差といった ランダム効果を扱いたい もっと自由で 現実的な 統計モデルを！ MCMCによる事後分布の推定 GLMM GLM GLMでは上手く 説明できないデータ

Slide 13

Slide 13 text

目的 個体差ありデータへのアプローチ 13 最小二乗法 線形モデル 最尤推定法 一般化線形モデル 一般化線形混合モデル 階層ベイズモデル 正規分布以外の 確率分布を扱いたい 個体差・場所差といった ランダム効果を扱いたい もっと自由で 現実的な 統計モデルを！ MCMCによる事後分布の推定 GLMM GLM 次に検討するモデル (今回のテーマ)

Slide 14

Slide 14 text

目的 実験結果からモデルを改善 14 今回検討するモデル ‐ 個体差を考慮していない ‐ データのばらつきに頑丈でない ‐ 個体差も考慮 ‐ データの不確実性を加味 生存種子数 個体 𝑖 点推定 生存種子数 個体 𝑖 分布推定 実験で用いたモデル

Slide 15

Slide 15 text

目次 15 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ

Slide 16

Slide 16 text

GLMMの導入 + ベイズモデル化 植物実験における課題 16 植物個体 𝑖 種子数 𝑛𝑖 = 8 生存数 𝑦𝑖 = 3 100個体分実験 生存数 個 体 数 個体差あり 種子 抽出 調査 生存数の観測分布 可視化 個体差を含めた生存種子数 の分布を推定したい

Slide 17

Slide 17 text

GLMMの導入 + ベイズモデル化 ベースモデルのベイズ化 17 ‐ 最尤推定 (点推定) はデータのばらつきに対して頑丈でない 不確実性を考慮 生存種子数 生存種子数 個体C 個体B 個体A 点推定 個体A 個体B 個体C 分布推定

Slide 18

Slide 18 text

GLMMの導入 + ベイズモデル化 一般化線形混合モデル (GLMM) による推定 18 logit 𝑞𝑖 = 𝛽 + 𝑟𝑖 植物個体 𝑖 100個体分実験 個体差あり ‐ 個体 𝑖 の生存確率 𝑞𝑖 の設計 ロジット関数 全個体共通の パラメータ (切片) 個体差の影響考慮 植物個体差の影響も 考えて生存種子数を 推定しよう GLMMの立場

Slide 19

Slide 19 text

GLMMの導入 + ベイズモデル化 GLMMの詳細 19 ロジット関数 logit 𝑞𝑖 = 𝛽 + 𝑟𝑖 全個体共通の パラメータ 個体 𝑖 の個体差 パラメータ 𝑟𝑖 ~𝑁 0, 𝑠2 0 𝑠 𝑠 新たに追加した パラメータ 個体 𝑖 の 生存確率 ※ GLMM : Generalized Linear Mixed Model

Slide 20

Slide 20 text

GLMMの導入 + ベイズモデル化 生存した種子が得られる確率 20 𝑝 𝒀|𝛽, 𝑟𝑖 = ෑ 𝑖=1 100 8 𝑦𝑖 𝑞𝑖 𝑦𝑖 1 − 𝑞𝑖 8−𝑦𝑖 尤度 二項分布 ‐ 個体数𝑁 = 100 ‐ 種子数 𝑛𝑖 = 8 ‐ 生存数 𝑦𝑖 ‐ 𝑟𝑖 = 𝑟1 , 𝑟2 , ⋯ , 𝑟100 ‐ 全100個体の二項分布の積

Slide 21

Slide 21 text

GLMMの導入 + ベイズモデル化 事後分布の推定 21 事後分布 ∝ 𝑝 𝒀|𝛽, 𝑟𝑖 ×事前分布 尤度 確率を更新 未定義

Slide 22

Slide 22 text

GLMMの導入 + ベイズモデル化 事前分布の設計 22 𝑝 𝛽 = 1 2𝜋 × 1002 exp −𝛽2 2 × 1002 ‐ 線形予測子の切片 𝛽 の事前分布 ‐ 個体差 𝑟𝑖 の事前分布 𝑝 𝑟𝑖 |𝑠 = 1 2𝜋𝑠2 exp −𝑟𝑖 2 2𝑠2 𝑟𝑖 ~𝑁 0, 𝑠2 0 𝑠 𝑠 𝛽~𝑁 0, 1002 平たい 𝑠の事前分布も 導入したい

Slide 23

Slide 23 text

GLMMの導入 + ベイズモデル化 超パラメータ 𝒔 の事前分布の設計 23 𝑝 𝑠 0から𝟏𝟎𝟒までの連続一様分布 0 𝟏𝟎𝟒 面積1 𝑟𝑖 ~𝑁 0, 𝑠2 𝑠の事前分布も導入したい 𝑠~𝑈 0, 104 事前分布の 事前分布を導入 𝑝 𝑟𝑖 |𝑠 : 階層事前分布 𝑝 𝑠 : 超事前分布

Slide 24

Slide 24 text

GLMMの導入 + ベイズモデル化 階層ベイズモデル 24 𝑠~Uniform 𝑟𝑖 ~Normal 𝛽~Normal 𝑦𝑖 ~Binomial 𝑞𝑖 ~Determistic 100 階層事前分布 超事前分布 観測変数 潜在変数 決定論的に 決まる変数 𝑝 𝛽, 𝑠, 𝑟𝑖 |𝒀 ∝ 𝑝 𝒀|𝛽, 𝑟𝑖 𝑝 𝛽 ෑ 𝑖=1 100 𝑝 𝑟𝑖 |𝑠 𝑝 𝑠 事後分布 尤度 事前分布 ‐ 個体数𝑁 = 100 ‐ 種子数 𝑛𝑖 = 8 ‐ 生存数 𝑦𝑖 ‐ 𝑟𝑖 = 𝑟1 , 𝑟2 , ⋯ , 𝑟100

Slide 25

Slide 25 text

GLMMの導入 + ベイズモデル化 補足1/2 : 事前分布の種類と選択 25 主観的な事前分布 無情報事前分布 階層事前分布 パラメータ 𝑥の事前 分布はこうなると思う パラメータ 𝑥の事前 分布は分からないわ パラメータ 𝑥の事前分布は パラメータ𝑦の事後分布 に支配されるよ

Slide 26

Slide 26 text

GLMMの導入 + ベイズモデル化 補足2/2 : 統計モデルのパラメータと事前分布 26 パラメータの種類 説明する範囲 パラメータの個数 採用すべき事前分布 全体に共通するばらつき 大域的 少数 無情報事前分布 個別のばらつき 局所的 多数 階層事前分布 生存種子数の実験の場合 ) ‐ 線形予測子𝛽 : 全個体共通のばらつき → 無情報事前分布が有効 ‐ 個体差 𝑟𝑖 : 個別のばらつき → 階層事前分布が有効

Slide 27

Slide 27 text

目次 27 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ

Slide 28

Slide 28 text

数値実験 : 個体差 MCMCサンプリングによる実験 28 ‐ Step数 ➢ Iterations : 1600回 ‐ 反復数 (サンプル列) ➢ Chains : 3回 ‐ 最初の何Stepの結果を使わないか ➢ Burn-in : 100回 ‐ 何個飛ばしてサンプリング結果を記録するか (可視化用) ➢ Skips : 3回 ‐ 正味のStep数 ➢ 1500回 ‐ 可視化 ➢ 500回 ※複数パラメータによる 事後分布の推定には ギブスサンプリングなどが 用いられる (Appendix 2)

Slide 29

Slide 29 text

数値実験 : 個体差 事後分布推定結果 29 MCMCサンプリング 𝛽の事後分布 𝑠の事後分布 𝑟𝑖 の事後分布 𝑖 = 1~100 (500回) 曲線はサンプルされたヒストグラムに カーネル密度推定を施したもの

Slide 30

Slide 30 text

数値実験 : 個体差 事後分布から個体毎の生存確率を予測 30 𝑟𝑖 の事後分布 𝑞𝑖 の分布 ‐ 𝑟𝑖 が小さい値をとりやすいと𝑞𝑖 の生存確率は小さくなる ※図は近似曲線だが実際はヒストグラム 𝑞𝑖 = 1 1 + exp −𝛽 − 𝑟𝑖 ‐ GLMM : logit 𝑞𝑖 = 𝛽 + 𝑟𝑖

Slide 31

Slide 31 text

数値実験 : 個体差 事後分布から個体毎の生存確率を予測 31 𝑟𝑖 の事後分布 𝑞𝑖 の分布 ※図は近似曲線だが実際はヒストグラム ‐ 𝑟𝑖 が大きい値をとりやすいと𝑞𝑖 の生存確率は大きくなる 𝑞𝑖 = 1 1 + exp −𝛽 − 𝑟𝑖 ‐ GLMM : logit 𝑞𝑖 = 𝛽 + 𝑟𝑖

Slide 32

Slide 32 text

数値実験 : 個体差 事後分布から個体毎の生存確率を予測 32 𝑟𝑖 の事後分布 𝑞𝑖 の分布 ※図は近似曲線だが実際はヒストグラム ‐ 𝑟𝑖 が0付近だと 𝑞𝑖 は0.5を中心とした分布となる 𝑞𝑖 = 1 1 + exp −𝛽 − 𝑟𝑖 ‐ GLMM : logit 𝑞𝑖 = 𝛽 + 𝑟𝑖

Slide 33

Slide 33 text

数値実験 : 個体差 生存種子数の予測分布 33 𝑝 𝑦|𝛽, 𝑠 = න −∞ ∞ 𝑝 𝑦|𝛽, 𝑟 𝑝 𝑟|𝑠 𝑑𝑟 生存種子数 𝑦 の確率分布 二項分布 個体差の事後分布 ‐ 二項分布を事後分布で平均化

Slide 34

Slide 34 text

数値実験 : 個体差 実装上の予測分布 34 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑠𝑖𝑗 = ෍ 𝑠=−10𝑠𝑖𝑗 𝑠=10𝑠𝑖𝑗 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑟𝑖𝑗 𝑝 𝑟𝑖𝑗 |𝑠 生存種子数が𝑦𝑖𝑗𝑘 となる確率 二項分布 超パラメータ 𝑠 の 事後分布を標準偏 差に持つ正規分布 ‐ MCMCの事後分布サンプルから計算 ‐ 𝑖 ∈ 1,2,3 : MCMCの反復 index ‐ 𝑗 ∈ 1,2, ⋯ , 1500 : MCMCの sampling index ‐ 𝑘 ∈ 0,1, ⋯ , 8 : 生存種子数のindex ‐ 𝑦𝑖𝑗𝑘 ∈ 0,1, ⋯ , 8 : 生存種子数が 𝑘 個 𝑝 𝑟𝑖𝑗 |𝑠 ~ 𝑁 0, 𝑠2

Slide 35

Slide 35 text

数値実験 : 個体差 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𝑠𝑖𝑗 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑟𝑖𝑗 𝑝 𝑟𝑖𝑗 |𝑠

Slide 36

Slide 36 text

数値実験 : 個体差 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 : 生存種子数が 𝑘 個

Slide 37

Slide 37 text

数値実験 : 個体差 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𝑠𝑖𝑗 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑟𝑖𝑗 𝑝 𝑟𝑖𝑗 |𝑠

Slide 38

Slide 38 text

数値実験 : 個体差 Step4. 予測分布から生存種子数を予測 38 ‐ 全ての 𝑖, 𝑗, 𝑘 について数値積分から予測分布 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑠 を計算 Step1 Step2 Step3 Step4 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑟𝑖𝑗 = 8 𝑦𝑖𝑗𝑘 𝑞𝑖𝑗 𝑦𝑖𝑗𝑘 1 − 𝑞𝑖𝑗 8−𝑦𝑖𝑗𝑘 𝑝 𝑟𝑖𝑗 |𝑠𝑖𝑗 = 1 2𝜋𝑠2 exp −𝑟𝑖𝑗 2 2𝑠2 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑠𝑖𝑗 = ෍ 𝑠=−10𝑠𝑖𝑗 𝑠=10𝑠𝑖𝑗 𝑝 𝑦𝑖𝑗𝑘 |𝛽𝑖𝑗 , 𝑟𝑖𝑗 𝑝 𝑟𝑖𝑗 |𝑠

Slide 39

Slide 39 text

数値実験 : 個体差 予測分布の別表現 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 : 生存種子数が 𝑘 個

Slide 40

Slide 40 text

数値実験 : 個体差 得られた予測分布から生存種子をサンプリング 40 𝑦11𝑘 (𝑖 = 𝑗 = 1) の予測分布可視化 生存種子数のヒストグラム e.g. ) 100個体 サンプリング

Slide 41

Slide 41 text

数値実験 : 個体差 全ての 𝑖, 𝑗 について生存種子数の分布推定 41 最尤推定で得られた生存種子数の分布 階層ベイズモデルから得られた生存種子数の分布

Slide 42

Slide 42 text

目次 42 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ

Slide 43

Slide 43 text

個体差 + 場所差の階層ベイズモデル 実験設定 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個) 各個体の生存種子数 𝑦𝑖 を観測

Slide 44

Slide 44 text

個体差+場所差の階層ベイズモデル 観測された種子数の Swarmplot と Boxplot 44

Slide 45

Slide 45 text

個体差+場所差の階層ベイズモデル 観測された種子数の Swarmplot と Boxplot 45 肥料なし 植木鉢内の個体差に散らばりがある 肥料あり

Slide 46

Slide 46 text

個体差+場所差の階層ベイズモデル 観測された種子数の Swarmplot と Boxplot 46 肥料なし 植木鉢内の個体差に散らばりがある 肥料あり 特定の植木鉢だけ種子数が少ない

Slide 47

Slide 47 text

個体差+場所差の階層ベイズモデル 個体差と植木鉢差 (場所差) の発生 47 植木鉢差 個体差 ‐ 個体差 + 植木鉢差を考慮した種子数の予測モデルを考える

Slide 48

Slide 48 text

個体差+場所差の階層ベイズモデル 種子数のばらつきの設定 48 𝑝 𝑦𝑖 |𝜆𝑖 = 𝜆𝑖 𝑦𝑖exp −𝜆𝑖 𝑦𝑖 ! ~ Poi 𝑦𝑖 |𝜆𝑖 ‐ 𝜆𝑖 : ポアソン分布の平均 (平均種子数) ‐ 𝑦𝑖 : 個体 𝑖 の種子数 ‐ ポアソン分布に従うと仮定 i. i. d

Slide 49

Slide 49 text

個体差+場所差の階層ベイズモデル GLMMの設計 49 対数リンク関数 log 𝜆𝑖 = 𝛽1 + 𝛽2 𝑓2 + 𝑟𝑖 + 𝑟𝑗 𝑖 全個体共通の パラメータ 個体 𝑖 の個体差 𝑓2 = 0で肥料なし 𝑓2 = 1で肥料あり 個体 𝑖 の 平均種子数 施肥処理 の効果 植木鉢 𝑗 の効果 新たに追加した パラメータ

Slide 50

Slide 50 text

個体差+場所差の階層ベイズモデル 事前分布の設計 50 𝛽1 , 𝛽2 は大域的なパラメータ → 無情報事前分布 𝑟𝑖 , 𝑟𝑗 𝑖 は局所的なパラメータ → 階層事前分布 𝑟𝑖 ~𝑁 0, 𝑠2 , 𝑟𝑗 𝑖 ~𝑁 0, 𝑠𝑝 2 𝛽1 , 𝛽2 ~𝑁 0, 1002 平たい 𝑠, 𝑠𝑝 は大域的なパラメータ → 無情報事前分布 𝑠, 𝑠𝑝 ~𝑈 0, 104 𝑟𝑖 , 𝑟𝑗 𝑖 の超パラメータ

Slide 51

Slide 51 text

個体差+場所差の階層ベイズモデル 階層ベイズモデル 51 𝑠~Uniform 𝑟𝑖 ~Normal 𝛽1 ~Normal 𝑦𝑖 ~Poisson 𝜆𝑖 ~Determistic 100 階層事前分布 超事前分布 𝛽2 ~Normal 𝑠𝑝 ~Uniform 𝑟 𝑗(𝑖) ~Normal 10 𝑓2(𝑖) 観測変数 潜在変数 決定論的に 決まる変数 𝑝 𝛽1 , 𝛽2 , 𝑠, 𝑠𝑝 , 𝑟𝑖 , 𝑟 𝑗(𝑖) |𝒀 ∝ 𝑝 𝒀|𝛽1 , 𝛽2 , 𝑠, 𝑠𝑝 , 𝑟𝑖 , 𝑟 𝑗(𝑖) ෑ 𝑖=1 100 𝑝 𝑟𝑖 |𝑠 ෑ 𝑗=1 10 𝑝 𝑟 𝑗(𝑖) |𝑠𝑝 𝑝 𝛽1 𝑝 𝛽2 𝑝 𝑠 𝑝 𝑠𝑝

Slide 52

Slide 52 text

目次 52 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ

Slide 53

Slide 53 text

数値実験 : 個体差+場所差 MCMCサンプリングによる実験 53 ‐ Step数 ➢ Iterations : 1600回 ‐ 反復数 (サンプル列) ➢ Chains : 3回 ‐ 最初の何Stepの結果を使わないか ➢ Burn-in : 100回 ‐ 何個飛ばしてサンプリング結果を記録するか (可視化用) ➢ Skips : 3回 ‐ 正味のStep数 ➢ 1500回 ‐ 可視化 ➢ 500回

Slide 54

Slide 54 text

数値実験 : 個体差+場所差 事後分布推定結果 1/2 54 𝛽1 の事後分布 𝛽2 の事後分布 𝑟𝑖 の事後分布

Slide 55

Slide 55 text

数値実験 : 個体差+場所差 事後分布推定結果 2/2 + 𝜆𝑖 の分布 55 𝑟𝑗(𝑖) の 事後分布 𝑠 (青)と 𝑠𝑝 (オレンジ) の事後分布 𝜆𝑖 の分布

Slide 56

Slide 56 text

数値実験 : 個体差+場所差 事後分布の推定から分かること 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となるように作られているので検定結果は正しい

Slide 57

Slide 57 text

目次 57 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ

Slide 58

Slide 58 text

空間構造がある階層ベイズモデル 一次元空間上の個体数分布 58 個体数 𝑦𝑗 = 2 区間 𝑗 = 1 1 2 3 3 4 2 0 50 ‐ 50個の区間調査 : 𝑗 ∈ 1, 2, ⋯ , 50 ‐ 区間 𝑗 における個体数 : 𝑦𝑗 (架空データ) 一次元空間 ‐ 𝑦𝑗 の従う分布 : 𝑝 𝑦𝑗 |𝜆𝑗 = 𝜆𝑗 𝑦𝑗exp −𝜆𝑗 𝑦𝑗 ! ~ Poi 𝑦𝑗 |𝜆𝑗

Slide 59

Slide 59 text

空間構造がある階層ベイズモデル 区間同士の隣接 59 個体数 𝑦𝑗 = 2 区間 𝑗 = 1 1 2 3 3 4 2 0 50 区間同士の幅を狭める

Slide 60

Slide 60 text

空間構造がある階層ベイズモデル 区間相関の発生 60 区間相関 近い区間同士の 傾向は似ている 区間毎の独立性を仮定できない！ 𝑦𝑗 ~ Poi 𝑦𝑗 |𝜆𝑗 i. i. d

Slide 61

Slide 61 text

空間構造がある階層ベイズモデル 観測個体数の分布 61 一次元空間上の架空データ 区間毎のポアソン分布の 真の平均 𝜆𝑗 観測数の変化は区間毎に滑らか 空間相関がある観測データ

Slide 62

Slide 62 text

空間構造がある階層ベイズモデル GLMMの設計 62 対数リンク関数 log 𝜆𝑗 = 𝛽 + 𝑟𝑗 全個体共通の パラメータ 区間 𝑗 毎の場所差 区間 𝑗 の平均個体数 大域的 無情報事前分布 𝛽~𝑁 0, 1002 局所的 階層事前分布 𝑠~𝑈 0, 104

Slide 63

Slide 63 text

空間構造がある階層ベイズモデル 近傍区間からの影響 63 注目区間 𝑦𝑗 = 3 𝑦𝑗−1 = 1 𝑦𝑗+1 = 2 近傍区間の設計 隣接区間のみの空間相関を考慮 区間 𝑗 𝑗 − 1 𝑗 + 1

Slide 64

Slide 64 text

空間構造がある階層ベイズモデル 空間構造のある階層事前分布 64 𝑝 𝑟𝑗 |𝜇𝑗 , 𝑠 = 𝑛𝑗 2𝜋𝑠2 exp − 𝑟𝑗 − 𝜇𝑗 2 2𝑠2/𝑛𝑗 ~ 𝑁 𝜇𝑗 , 𝑠2 𝑛𝑗 ‐ 場所差 𝑟 𝑗 の事前分布 ‐ 𝜇𝑗 = 𝑟𝑗−1+𝑟𝑗+1 2 : 近傍区間から影響を考慮 ‐ 𝜇1 = 𝑟2 , 𝜇50 = 𝑟49 : 両端区間の処理 ‐ 𝑛𝑗 : 近傍区間の個数 (実験では2個)

Slide 65

Slide 65 text

空間構造がある階層ベイズモデル 条件付き自己回帰モデル CAR 65 𝑝 𝑟𝑗 |𝑠 ∝ exp 1 2𝑠2 ෍ 𝑗~𝑗′ 𝑟𝑗 − 𝑟𝑗′ 2 ~ Gaussian 𝐂𝐀𝐑 ‐ 場所差 𝑟 𝑗 = 𝑟1 , 𝑟2 , ⋯ , 𝑟50 全体の (同時) 事前分布 ‐ 𝑗~𝑗′: 区間 𝑗 と近傍区間 𝑗′ の全ての組合せ 条件付き自己回帰 ※ CAR : Conditional Auto Regressive

Slide 66

Slide 66 text

空間構造がある階層ベイズモデル 実装上のCARモデルの事前分布 66 𝑟 𝑗 ~ 𝑁 𝟎, 𝑠2 𝐷 − 𝛼𝑊 −1 ‐ Brook’s (1964) Lemma による設計方法を採用 ‐ 𝐷 = diag 𝑛𝑗 : 隣接区間の個数 𝑛𝑗 を成分に持つ50 × 50 対角行列 ‐ 𝛼 : 空間依存性を制御するパラメータ ‐ 𝑊 = 𝑤𝑗𝑗 = 0, 𝑤𝑗′𝑗 = 1 if 𝑗′ is neighbor of 𝑗, and 𝑤𝑗′𝑗 = 0 otherwise : 𝑤𝑖𝑗 を成分に持つ50 × 50 隣接行列 ‐ 𝛼 にも超事前分布を設計 : 𝛼 ~ 𝑈 0, 104

Slide 67

Slide 67 text

空間構造がある階層ベイズモデル 階層ベイズモデル 67 𝑠~Uniform 𝛽~Normal 𝑦𝑗 ~Poisson 𝜆𝑗 ~Determistic 50 階層事前分布 超事前分布 観測変数 潜在変数 決定論的に 決まる変数 𝑟 𝑗 ~Gaussian CAR 𝑝 𝛽, 𝑠, 𝑟 𝑗 |𝒀 ∝ ෑ 𝑗=1 50 𝑝 𝑦𝑗 |𝜆𝑗 𝑝 𝑟 𝑗 |𝑠, 𝛼 𝑝 𝑠 𝑝 𝛽 𝑝 𝛼 事後分布 尤度 事前分布 𝛼~Uniform

Slide 68

Slide 68 text

目次 68 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ

Slide 69

Slide 69 text

数値実験 : 空間構造 MCMCサンプリングによる実験 69 ‐ Step数 ➢ Iterations : 60000回 ‐ 反復数 (サンプル列) ➢ Chains : 3回 ‐ 最初の何Stepの結果を使わないか ➢ Burn-in : 10000回 ‐ 何個飛ばしてサンプリング結果を記録するか (可視化用) ➢ Skips : 1回 ‐ 正味のStep数 ➢ 50000回 ‐ 可視化 ➢ 50000回

Slide 70

Slide 70 text

数値実験 : 空間構造 事後分布および自己相関推定結果 70 MCMCサンプリング 事後分布 自己相関 実行時平均 ※ 𝑟 𝑗 の事後分布の可視化は割愛

Slide 71

Slide 71 text

数値実験 : 空間構造 事後分布から個体数推定 71 𝜆𝑗 = exp 𝛽 + 𝑟𝑗 事後分布からサンプリングした個体数平均𝜆𝑗 の分布 ‐ 50000点のサンプリング Posterior sampling GLMM 空間相関を上手く表現できている

Slide 72

Slide 72 text

数値実験 : 空間構造 データの一部を欠測させて再実験 72 一次元空間上の架空データ データの一部を欠測 欠測させたデータ

Slide 73

Slide 73 text

数値実験 : 空間構造 欠測データに対する比較実験 空間相関無し 空間相関有り 73

Slide 74

Slide 74 text

数値実験 : 空間構造 欠測データに対する比較実験 74 空間相関無し 空間相関有り

Slide 75

Slide 75 text

数値実験 : 空間構造 階層ベイズモデル : 空間相関無し 75 𝑠~Uniform 𝛽~Normal 𝑦𝑗 ~Poisson 𝜆𝑗 ~Determistic 50 階層事前分布 超事前分布 観測変数 潜在変数 決定論的に 決まる変数 𝑟 𝑗 ~Normal 𝑝 𝛽, 𝑠, 𝑟 𝑗 |𝒀 ∝ ෑ 𝑗=1 50 𝑝 𝑦𝑗 |𝜆𝑗 ෑ 𝑖=1 50 𝑝 𝑟𝑖 |𝑠 𝑝 𝑠 𝑝 𝛽 事後分布 尤度 事前分布 𝑝 𝑟𝑗 |𝑠 ~ 𝑁 0, 𝑠2 i. i. d

Slide 76

Slide 76 text

数値実験 : 空間構造 事後分布および自己相関推定結果 : 空間相関無し 76 MCMCサンプリング 事後分布 自己相関 実行時平均 ※ 𝑟 𝑗 の事後分布の可視化は割愛

Slide 77

Slide 77 text

数値実験 : 空間構造 欠測データに対する比較実験 空間相関無し 空間相関有り 77

Slide 78

Slide 78 text

数値実験 : 空間構造 階層ベイズモデル : 空間相関有り 78 𝑠~Uniform 𝛽~Normal 𝑦𝑗 ~Poisson 𝜆𝑗 ~Determistic 50 階層事前分布 超事前分布 観測変数 潜在変数 決定論的に 決まる変数 𝑟 𝑗 ~Gaussian CAR 𝑝 𝛽, 𝑠, 𝑟 𝑗 |𝒀 ∝ ෑ 𝑗=1 50 𝑝 𝑦𝑗 |𝜆𝑗 𝑝 𝑟 𝑗 |𝑠, 𝛼 𝑝 𝑠 𝑝 𝛽 𝑝 𝛼 事後分布 尤度 事前分布 𝛼~Uniform 𝑝 𝑟 𝑗 |𝑠 ~𝑁 𝜇𝑗 , 𝑠2 𝑛𝑗

Slide 79

Slide 79 text

数値実験 : 空間構造 事後分布および自己相関推定結果 :空間相関有り 79 MCMCサンプリング 事後分布 自己相関 実行時平均 ※ 𝑟 𝑗 の事後分布の可視化は割愛

Slide 80

Slide 80 text

数値実験 : 空間構造 欠測付き観測データの推定結果 80 空間相関無しの階層ベイズモデル 空間相関ありの階層ベイズモデル 欠測区間の分散が大きい 欠測区間の分散が小さい

Slide 81

Slide 81 text

数値実験 : 空間構造 補足 : 二次元以上の空間構造 81 ‐ CARモデルの次元や隣接区間の個数は変化するが考え方は同じ

Slide 82

Slide 82 text

目次 82 1. 背景・目的 2. GLMMの導入 + ベイズモデル化 3. 数値実験 4. 個体差+場所差の階層ベイズモデル 5. 数値実験 : 個体差+場所差 6. 空間構造がある階層ベイズモデル 7. 数値実験 : 空間構造 8. まとめ

Slide 83

Slide 83 text

まとめ 83 ‐ 階層ベイズモデルとは, 事前分布となる確率分布のパラメータにも事前分布を導入 したベイズ統計モデルである ‐ 無情報事前分布・階層事前分布による事前分布の設計方法を扱った ‐ 個体差+場所差などの複雑なデータ構造の統計モデリングでは, 階層ベイズモデル とMCMCサンプリングによるパラメータ推定がよい ‐ データの独立性を仮定できない時は, 空間相関を考慮したモデル設計が必要である

Slide 84

Slide 84 text

Appendix 1 二項分布に対する最尤法 84 log 𝐿 𝑞 = ෍ 𝑖=1 𝑛 log 𝑀 𝑦𝑖 + 𝑦𝑖 log 𝑞 + 𝑀 − 𝑦𝑖 log 1 − 𝑞 𝜕 𝜕𝑞 log 𝐿 𝑞 = 𝜕 𝜕𝑞 log 𝑞 ෍ 𝑖=1 𝑛 𝑦𝑖 + 𝜕 𝜕𝑞 log 1 − 𝑞 ෍ 𝑖=1 𝑛 𝑀 − 𝑦𝑖 = 1 𝑞 ෍ 𝑖=1 𝑛 𝑦𝑖 − 1 1 − 𝑞 ෍ 𝑖=1 𝑛 𝑀 − 𝑦𝑖 = 𝑛ത 𝑦 𝑞 − 𝑛𝑀 − 𝑛ത 𝑦 1 − 𝑞 = 0 ത 𝑦 = 1 𝑛 ෍ 𝑖=1 𝑛 𝑦𝑖 を解く. ∴ 𝑞 = ത 𝑦 𝑀 を得る. の最大化. 𝑞 𝐿 𝑞 最尤推定値

Slide 85

Slide 85 text

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++コンパイラのインストールが必要

Slide 86

Slide 86 text

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​)

Slide 87

Slide 87 text

Appendix 4 実験環境 : 空間構造のMCMCサンプリング : 𝛼, 𝛽, 𝑠 88 ‐ MCMC : Slice (Shrinkage Slice​)

Slide 88

Slide 88 text

Appendix 5 実験環境 : 欠測付きMCMCサンプリング : 𝑦𝑗 89 ‐ MCMC : MISS (Missing Value Sampler​)

Slide 89

Slide 89 text

参考文献 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.