ベイズ統計_11_ベイズ階層モデリング1

ベイズ統計 11 ベイズ階層モデリング(1) 分寺 杏介 神戸大学大学院 経営学研究科  [email protected] ※本スライドは，クリエイティブ・コモンズ 表示-非営利 4.0 国際 ライセンス（CC BY-NC 4.0）に従って利用が可能です。

• https://creativecommons.org/licenses/by-nc/4.0/

（前回のおさらい）地域差に関する２つの仮説モデル DISTの影響が地域ごとに異なるのか同じなのか 知りたいんだよねぇ。 ま，うまいことやっといてよ！ 回帰係数には 地域差がある 𝐻1 : 𝐻0 以外 model { （事前分布は省略） for(i in 1:N){ SALES[i] ~ normal(beta_0 + beta_DIST[REGION[i]]*DIST[i], sigma); } } 地域ごとに異なる係数 v.s. 回帰係数には地域差はない 𝐻0 : 𝛽DIST,1 = 𝛽DIST,2 = ⋯ = 𝛽DIST,𝐺 model { （事前分布は省略） for(i in 1:N){ SALES[i] ~ normal(beta_0 + beta_DIST*DIST[i], sigma); } } 今回のお話はここから始まります 11 ベイズ階層モデリング(1) 2

１ 階層モデルとは またの名をマルチレベルモデル 11 ベイズ階層モデリング(1) 3

データには階層性があることが多い ▌大抵のデータに対して，階層性を考えることもできるし無視できることもある コンビニ 店舗 地域 企業 従業員 課 部 本部 学校 生徒 学級 学校 地域 企業 ▌基本的には入れ子（ネスト）の関係にあると考える 【コンビニの場合】 地域A 地域B … 店舗1 店舗2 … 店舗10 店舗11 店舗12 … 店舗20 11 ベイズ階層モデリング(1) 上位カテゴリが複数（e.g., 部署と支社） ある場合や １つの個体が複数の上位カテゴリに 所属する（e.g., SNSのグループ）場合は また別のモデルが必要になりますが この講義では最もシンプルな のような形式のみを考えていきます 4

なぜ階層性に注意して分析する必要があるのか？ ▌個体の効果と集団の効果が区別できなくなってしまうため 例 勉強時間が長い学生ほど模試の成績が高いのかを調べることにしました。 以下の図は，勉強時間と模試の得点の散布図です。 の図では明らかに 模試の得点 勉強時間が長いほど成績が低下する ことが示されている 勉強しても 意味ねーじゃねーか 勉強時間 11 ベイズ階層モデリング(1) 5

なぜ階層性に注意して分析する必要があるのか？ ▌個体の効果と集団の効果が区別できなくなってしまうため 例 勉強時間が長い学生ほど模試の成績が高いのかを調べることにしました。 以下の図は，勉強時間と模試の得点の散布図です。 ちなみにこのデータは，地元の３つの高校から得られたものです。 グループごとに見てみると 模試の得点 勉強時間が長いほど成績が向上する ことが示される 高校 学力 勉強時間 A 屈指の進学校 短め B 普通の公立 ふつう C … 長め やっぱ 勉強したほうが いいんじゃねーか 勉強時間 グループ間のベースラインの違いのせいで 全体では結果が違って見えてしまった 11 ベイズ階層モデリング(1) シンプソンのパラドックス 6

階層性にどう対処するか 高校B 模試の得点 高校A 模試の得点 模試の得点 ▌集団ごとに分析する？ 高校C 勉強時間 勉強時間 結果をまとめるのがめんどくさそう 勉強時間 グループの数が100とかなったらどうするの？ 「全体的な傾向」があるかを確認しづらい e.g., （高校ごとに違いはあれど）基本的には「勉強時間が長いほど得点が高い」か？ そして 「全体的な傾向」を確認するなら全データを使って推定したほうが精度も良い ということでマルチレベルモデルを使いましょう 11 ベイズ階層モデリング(1) データを分けるという行為は 多くの場合，単純に勿体ない 7

マルチレベルモデルの定式化 ▌（参考）グループで差がない：通常の回帰分析 𝑦𝑔𝑖 = 𝛽0 + 𝛽1 𝑥𝑔𝑖 + 𝜀𝑔𝑖 ▌グループごとに切片が異なる（ランダム切片モデル） （レベル1） 𝑦𝑔𝑖 = 𝛽0𝑔 + 𝛽1 𝑥𝑔𝑖 + 𝜀𝑔𝑖 （レベル2） 𝛽0𝑔 = 𝜇𝛽0 + 𝑢0𝑔 𝑢0𝑔 ∼ 𝑁(0, 𝜎0 ) 𝑦𝑔𝑖 = 𝜇𝛽0 + 𝛽1 𝑥𝑔𝑖 + 𝜀𝑔𝑖 + 𝑢0𝑔 ▌グループごとに傾きも異なる（ランダム切片・傾きモデル） （レベル1） （レベル2） 𝑦𝑔𝑖 = 𝛽0𝑔 + 𝛽1𝑔 𝑥𝑔𝑖 + 𝜀𝑔𝑖 𝛽0𝑔 = 𝜇𝛽0 + 𝑢0𝑔 𝑢0𝑔 ∼ 𝑁 0, 𝜎0 𝛽1𝑔 = 𝜇𝛽1 + 𝑢1𝑔 𝑢1𝑔 ∼ 𝑁(0, 𝜎1 ) 11 ベイズ階層モデリング(1) 𝑦𝑔𝑖 = 𝜇𝛽0 + 𝜇𝛽1 + 𝑢1𝑔 𝑥𝑔𝑖 + 𝜀𝑔𝑖 + 𝑢0𝑔 8

（補足）階層性を考慮する必要について数式から考える ▌数理的には「誤差に相関が生じてしまうため」 レベル1の誤差相関 𝐶𝑜𝑟 𝜀𝑔𝑖 , 𝜀𝑔𝑗 = 0と仮定します マルチレベルモデルでは各係数を「係数の全体平均𝜇𝛽⋅ 」と「 𝛽⋅ҧ からの差𝑢⋅𝑔 」に分けて 𝑦𝑔𝑖 = 𝜇𝛽0 + 𝑢0𝑔 + 𝜇𝛽1 + 𝑢1𝑔 𝑥𝑔𝑖 + 𝜀𝑔𝑖 ここで，グループごとの違いを無視して係数を推定してしまうと 𝑦𝑔𝑖 = 𝜇𝛽0 + 𝜇𝛽1 𝑥𝑔𝑖 + 𝜀𝑔𝑖 + 𝑢0𝑔 + 𝑢1𝑔 𝑥𝑔𝑖 このとき，2つのデータの誤差共分散は 𝐶𝑜𝑣 𝜀𝑔𝑖 + 𝑢0𝑔 + 𝑢1𝑔 𝑥𝑔𝑖 , 𝜀𝑔𝑗 + 𝑢0𝑔 + 𝑢1𝑔 𝑥𝑔𝑗 ≠ 0 11 ベイズ階層モデリング(1) したがって，誤差相関も0ではなくなる 9

２ stanで階層モデル（の準備） とりあえず実装からモデル比較まで 11 ベイズ階層モデリング(1) 10

こんな問題を考えてみます 例 あるコンビニチェーンのアナリストは，各店舗の駅からの距離が利益とどう関係するかを 調べることにしました。 ただ地域によって各変数が及ぼす影響や平均的な利益は異なる気がします。 この違いがわかれば，地域ごとに異なる戦略を立てることができそうです。 【データの読み込み】 この名前はなんでもOK "data_cvs.csv"をワーキングディレクトリに配置して ▌ 中身はこんなデータ dat <- read.csv("data_cvs.csv") 今日はこれらを使います sales: その店の一日あたり平均利益（単位：千円） dist: 最寄り駅からの距離（単位：km） 被説明変数 説明変数 floor: 床面積（単位：m2） items: 取扱いアイテム数（単位：個） region: その店舗のある地域 グループ変数 neighbor: 半径1km以内にあるコンビニの数 ※データは適当に作ったので、実際とは異なります。 11 ベイズ階層モデリング(1) 11

stanコードの拡張
▌以前作成した重回帰モデル▼をもとに

資料09 p.19
model_regression.stan

data {

model {

int N;

beta_0 ~ normal(0, 100);

array[N] real SALES;

beta_DIST ~ normal(0, 100);

vector[N] DIST;

beta_FLOOR ~ normal(0, 100);

vector[N] FLOOR;

beta_ITEMS ~ normal(0, 100);

vector[N] ITEMS;

sigma ~ cauchy(0, 10);

}

SALES ~ normal(beta_0 + beta_DIST*DIST + beta_FLOOR*FLOOR + beta_ITEMS*ITEMS, sigma);
}

parameters {
real beta_0;
real beta_DIST;

array[N] real SALES_pred;

real beta_FLOOR;

SALES_pred = normal_rng(beta_0 + beta_DIST*DIST + beta_FLOOR*FLOOR + beta_ITEMS*ITEMS, sigma);

real beta_ITEMS;
real <lower=0> sigma;
}

generated quantities {

}

後ほど事後予測チェックもするため，generated quantitiesブロックもそのまま

11 ベイズ階層モデリング(1)

12

少しコードを整理してみる ▌回帰式はどうしてもコードが長くなってしまうなぁ 尤度の式に挿入されるパラメータだけ先に作っておくという手があります （data, parametersブロックはそのまま） transformed parameters { vector[N] yhat; yhat = beta_0 + beta_DIST*DIST + beta_FLOOR*FLOOR + beta_ITEMS*ITEMS; } model { （事前分布はそのまま） SALES ~ normal(yhat, sigma); } generated quantities { array[N] real SALES_pred; SALES_pred = normal_rng(yhat, sigma); } 11 ベイズ階層モデリング(1) 【メリット】 • コードが多少見やすくなる • yhatの計算回数が減る＝僅かに高速化 ▲元のコードでは２回同じ計算をしていた 1. modelブロック 2. generated quantitiesブロック • yhatの事後分布が取り出しやすい 【デメリット】 • yhatのサンプルも保存するため ファイルサイズがすこし大きくなる 13

（補足）yhat自体は保存しなくても良いなら ▌計算上は同じことを行う別の書き方 modelブロックの中で宣言するという手もあります （data, parametersブロックはそのまま） model { vector[N] yhat; yhat = beta_0 + beta_DIST*DIST + beta_FLOOR*FLOOR + beta_ITEMS*ITEMS; （事前分布はそのまま） SALES ~ normal(yhat, sigma); } generated quantities { array[N] real SALES_pred; SALES_pred = normal_rng(beta_0 + beta_DIST*DIST + beta_FLOOR*FLOOR + beta_ITEMS*ITEMS, sigma); } 11 ベイズ階層モデリング(1) 【メリット】 • コードが多少見やすくなる • yhatを保存しないので ファイルサイズは変わらない 【デメリット】 • modelブロックで計算したyhatは generated quantitiesブロックでは使えない 14

ベースとなるコード
▌ここから階層モデルに書き換えていきます

仮説モデル 𝐻0 : 切片も傾きも地域差はない

model_multilevel_H0.stan
data {
int N;
array[N] real SALES;

transformed parameters {
vector[N] yhat;
yhat = beta_0 + beta_DIST*DIST + beta_FLOOR*FLOOR + beta_ITEMS*ITEMS;
}

vector[N] DIST;
vector[N] FLOOR;
vector[N] ITEMS;
}

parameters {

model {
beta_0 ~ normal(0, 100);
beta_DIST ~ normal(0, 100);
beta_FLOOR ~ normal(0, 100);
beta_ITEMS ~ normal(0, 100);
sigma ~ cauchy(0, 10);

real beta_0;
real beta_DIST;

SALES ~ normal(yhat, sigma);
}

real beta_FLOOR;
real beta_ITEMS;
real <lower=0> sigma;
}

generated quantities {
array[N] real SALES_pred;
SALES_pred = normal_rng(yhat, sigma);
}

11 ベイズ階層モデリング(1)

ブリッジサンプリングをしたい場合などは
target記法で書きましょう

15

グループごとに自由推定にする際のコードの書き方 ▌すでに説明していました 資料10 p.36 成したモデル 地域ごとに回帰係数を自由推定 【 】 dataブロックには グループ（とグループの数）を表す変数を追加 parametersブロックでは グループごとに変動するパラメータのサイズを変更 vector型やarray型を使う ということでp. 15のモデルを , すべての切片・傾きが グループごとに異なるように 0 書き換えてみましょう！ コード例はp. 18に modelブロックでは 添字をうまく このモデルについて 定することで 仮説 （ ）がどう 定されるか 各個体iが所属するグループごとに異なる係数を 考えてみましょう 定 モデル 11 ベイズ階層モデリング(1) 16

各モデルをplate notationで表すと 【グループに関係なく同じ係数】 【グループごとに異なる】 𝑔 group (REGION) 𝛽0𝑔 𝛽0 𝛽DIST 𝑔𝑖 F OOR𝑔𝑖 M 𝑔𝑖 𝑔𝑖 𝛽FLOOR 𝑖 data 𝛽ITEMS 𝛽DIST,𝑔 𝑔𝑖 F OOR𝑔𝑖 𝜎𝑒 M 𝑔𝑖 𝑔𝑖 𝛽FLOOR,𝑔 𝑖 data 𝛽ITEMS,𝑔 𝜎𝑒 11 ベイズ階層モデリング(1) 17

ランダム切片・傾きモデル
▌切片も含めてすべての係数がグループごとに異なる

仮説モデル 𝐻1 : 回帰係数には地域差がある

model_multilevel_H1.stan
data {
int N;

int G;
array[N] real SALES;
vector[N] DIST;

transformed parameters {
vector[N] yhat;
for(i in 1:N){
yhat[i] = beta_0[REGION[i]] + beta_DIST[REGION[i]]*DIST[i] +
beta_FLOOR[REGION[i]]*FLOOR[i] + beta_ITEMS[REGION[i]]*ITEMS[i];
}
}

vector[N] FLOOR;
vector[N] ITEMS;

array[N] int REGION;
}
parameters {

model {
beta_0 ~ normal(0, 100);
beta_DIST ~ normal(0, 100);
beta_FLOOR ~ normal(0, 100);
beta_ITEMS ~ normal(0, 100);
sigma ~ cauchy(0, 10);

vector[G] beta_0;

vector[G] beta_DIST;
vector[G] beta_FLOOR;
vector[G] beta_ITEMS;
real <lower=0> sigma;
}

SALES ~ normal(yhat, sigma);

}
generated quantities {
array[N] real SALES_pred;
SALES_pred = normal_rng(yhat, sigma);
}

11 ベイズ階層モデリング(1)

ブリッジサンプリングをしたい場合などは
target記法で書きましょう

18

ランダム切片モデルなら
▌切片だけグループごとに異なる
data {
int N;

int G;
array[N] real SALES;
vector[N] DIST;
vector[N] FLOOR;
vector[N] ITEMS;

array[N] int REGION;
}

transformed parameters {
vector[N] yhat;
for(i in 1:N){
yhat[i] = beta_0[REGION[i]] + beta_DIST*DIST[i] + beta_FLOOR*FLOOR[i] + beta_ITEMS*ITEMS[i];
}
}
model {
beta_0 ~ normal(0, 100);
beta_DIST ~ normal(0, 100);
beta_FLOOR ~ normal(0, 100);
beta_ITEMS ~ normal(0, 100);
sigma ~ cauchy(0, 10);

parameters {
SALES ~ normal(yhat, sigma);

vector[G] beta_0;

real beta_DIST;
real beta_FLOOR;
real beta_ITEMS;
real <lower=0> sigma;
}

}
generated quantities {
array[N] real SALES_pred;
SALES_pred = normal_rng(yhat, sigma);
}

11 ベイズ階層モデリング(1)

ブリッジサンプリングをしたい場合などは
target記法で書きましょう

19

推定してみる ▌推定 library(cmdstanr) stan_data <- list(N=100, G=10, SALES = dat$sales, DIST=dat$dist, FLOOR=dat$floor, ITEMS = dat$items, REGION = dat$region) # 係数がグループ間で同じモデル model_H0 <- cmdstan_model("model_multilevel_H0.stan") result_H0 <- model_H0$sample(data = stan_data) # 係数がグループ間で異なるモデル model_H1 <- cmdstan_model("model_multilevel_H1.stan") result_H1 <- model_H1$sample(data = stan_data) ▌（事後予測チェック用に）MCMCサンプルを取り出す SALES_pred_H0 <- result_H0$draws("SALES_pred", format="matrix")[seq(1,4000,by=500),] SALES_pred_H1 <- result_H1$draws("SALES_pred", format="matrix")[seq(1,4000,by=500),] 11 ベイズ階層モデリング(1) 20

事後予測チェック もちろんやり方は色々と考えられます（ベイズファクターや情報量規準を使うのも当然OK） ▌とりあえずデータの密度を比較してみると ppc_dens_overlay(dat$sales, SALES_pred_H0) ppc_dens_overlay(dat$sales, SALES_pred_H1) こちらのほうがyの密度に近いような気はするけど… 11 ベイズ階層モデリング(1) 21

もう少し詳細に ▌グループごとに異なるかも，と考えているのだから プロットもグループごとにやってもいいのでは？ グループごとに密度を出すので 各グループに1-2個体などの場合には使えませんが ppc_dens_overlay_grouped(dat$sales, SALES_pred_H0, group = dat$region) 地域によっては あまりうまくフィットしていない ところもありそう 11 ベイズ階層モデリング(1) 22

もう少し詳細に ▌グループごとに異なるかも，と考えているのだから プロットもグループごとにやってもいいのでは？ ppc_dens_overlay_grouped(dat$sales, SALES_pred_H1, group = dat$region) 前ページと行ったり来たりしながら 見比べてみてください 1サンプルを除いては モデルH1のほうが フィットしている？ 明らかにモデルH1のほうが よくフィットしている地域が見られる 11 ベイズ階層モデリング(1) 23

別の見方 ▌点推定値（EAP）とのズレを見るのもあり ppc_scatter_avg(dat$sales, result_H0$draws("SALES_pred", format="matrix")) EAPは正確に出したいので，間引きせず すべてのMCMCサンプルを使用 EAPによる事後予測が うまく行っているほど 点線の上に近づくはず 11 ベイズ階層モデリング(1) 24

別の見方 ▌点推定値（EAP）とのズレを見るのもあり ppc_scatter_avg(dat$sales, result_H1$draws("SALES_pred", format="matrix")) EAPは正確に出したいので，間引きせず すべてのMCMCサンプルを使用 ちなみにこいつにも _grouped()関数があります 前ページと見比べると 明らかにモデルH1のほうが 点線の近くにあるように見える 11 ベイズ階層モデリング(1) 25

グループごとに異なるかも，と考えているのだから あるかもしれない疑問 プロットもグループごとにやってもいいのでは？ 前ページと行ったり来たりしながら 見比べてみてください サンプルを除いては モデル のほうが フィットしている？ 明らかにモデル のほうが よくフィットしている地域が見られる 事後予測チェック的には モデルH1のほうが フィットしそうです 回帰係数は 地域ごとに異なるはず！ モデル 実際には回帰係数に差が無いときには モデルH0とモデルH1は同じ推定値になるはずだから 事後予測チェックでモデルH1のほうが当てはまりが悪くなる，ってことはないのでは？ 11 ベイズ階層モデリング(1) 26

ベイズ推定の原則を思い出す データが増えるほど推定の精度が向上する データが多いほど 尤度関数は尖っていきますね ▌回帰係数が全グループで同じ場合について考えてみると グループごとに異なる（モデルH1） と 定した場合 全グループで同じ（モデルH0） と 定した場合 1つの回帰係数を推定するのに 1つの回帰係数を推定するのに 10店舗のデータを使用している 100店舗のデータを使用している グループごとに推定される 回帰係数は結局同じ ただし事後SDは大きめになる 回帰係数の事後SDは H1のときより小さくなる MCMCサンプルごとに計算される 事後予測密度のバラツキが大きくなる MCMCサンプルごとに計算される 事後予測密度のバラツキも小さくなる 逆に言えばEAPなどはたぶん ほぼ同じになります こちらのほうが実データの密度に近くなるはず 11 ベイズ階層モデリング(1) 27

もっとうまくできないか？ ▌ここまでやってきたこと 1 地域 店舗 1 2 2 … 10 11 12 EAP 10 … 20 EAP 𝛽0 1.3117 𝛽0 1.7616 𝛽DIST -0.1917 𝛽DIST 𝛽FLOOR 0.0554 𝛽ITEMS -0.0011 … 91 92 … 100 EAP 𝛽0 0.7723 -0.1646 𝛽DIST -0.2158 𝛽FLOOR 0.0550 𝛽FLOOR 0.0586 𝛽ITEMS -0.0013 𝛽ITEMS -0.0012 ただ「10店舗による回帰分析」を10回繰り返しただけ 11 ベイズ階層モデリング(1) 28

他の地域の情報も使えたらなあ ▌本当にやりたいこと 1 地域 店舗 1 2 2 … 10 11 12 EAP 10 … 20 EAP 𝛽0 1.3117 𝛽0 1.7616 𝛽DIST -0.1917 𝛽DIST 𝛽FLOOR 0.0554 𝛽ITEMS -0.0011 例えば … 91 92 … 100 EAP 𝛽0 0.7723 -0.1646 𝛽DIST -0.2158 𝛽FLOOR 0.0550 𝛽FLOOR 0.0586 𝛽ITEMS -0.0013 𝛽ITEMS -0.0012 地域1の回帰係数 の推定に 地域1以外から得られる情報 を組み込めたら最高 11 ベイズ階層モデリング(1) もちろん他の地域でも同じ 29

３ ベイズ階層モデリング パラメータを階層的に考える 11 ベイズ階層モデリング(1) 30

まだ階層モデルとは言えなかった ▌本来階層モデルでは，各係数の「平均」と「バラツキ」を推定していた （レベル1） 𝑦𝑔𝑖 = 𝛽0𝑔 + 𝛽1𝑔 𝑥𝑔𝑖 + 𝜀𝑔𝑖 （レベル2） 𝛽0𝑔 = 𝜇𝛽0 + 𝑢0𝑔 𝛽1𝑔 = 𝜇𝛽1 + 𝑢1𝑔 全体での 平均的な傾向 𝑢0𝑔 ∼ 𝑁 0, 𝜎0 𝑢1𝑔 ∼ 𝑁(0, 𝜎1 ) 𝛽0𝑔 ∼ 𝑁 𝜇𝛽0 , 𝜎0 𝛽1𝑔 ∼ 𝑁 𝜇𝛽1 , 𝜎1 グループ間でのばらつき 事前分布のパラメータが 他の地域の係数に基づいて 調整される ▌これをベイズ統計的に考えると 𝑃 𝑌 𝛽⋅𝑔 𝑃 𝛽⋅𝑔 𝑃 𝛽⋅𝑔 𝑌 = 𝑃 𝑌 p. 8 𝛽⋅𝑔 ∼ 𝑁 𝜇𝛽⋅ , 𝜎⋅ 他の地域から得られる情報は 事前分布の一部として 利用されることになる！ 11 ベイズ階層モデリング(1) 事後分布は尤度×事前分布 に比例するのです 31

もう少し感覚的な説明 10地域のEAPから ーネル密度推定したものです 「全体的に傾き・切片はどの程度の値に 分布しているか」がわかってくる EAP 地域1 地域2 地域10 𝛽0 1.3117 1.7616 0.7723 𝛽DIST -0.1917 -0.1646 -0.2158 𝛽FLOOR 0.0554 0.0550 𝛽ITEMS -0.0011 -0.0013 𝛽0 の分布 0.0586 -0.0012 ーネル密度 … ーネル密度 各地域ごとに推定される回帰係数から 𝛽FLOOR の分布 11 ベイズ階層モデリング(1) 32

もう少し感覚的な説明 「全体的に傾き・切片はどの程度の値に cf. beta_FLOOR ~ normal(0, 100); 分布しているか」がわかれば 各地域の係数の推定に 情報を与えることができる これを事前分布的に使う ーネル密度 事後分布 𝛽FLOOR の分布 11 ベイズ階層モデリング(1) 33

▌といっても難しいことはありません ーネル密度 stanモデルを書いていく この式に沿って書けばよいだけ （レベル1） 𝑦𝑔𝑖 = 𝛽0𝑔 + 𝛽1𝑔 𝑥𝑔𝑖 + 𝜀𝑔𝑖 （レベル2） 𝛽0𝑔 = 𝜇𝛽0 + 𝑢0𝑔 𝛽1𝑔 = 𝜇𝛽1 + 𝑢0𝑔 全体での 平均的な傾向 𝑢0𝑔 ∼ 𝑁 0, 𝜎0 𝑢1𝑔 ∼ 𝑁(0, 𝜎1 ) グループ間でのばらつき カーネル密 推定したこの分布は そのまま事前分布には使えませんが， これを近似する正規分布を考えて あげると，通常の階層モデルの式と おなじになる，というわけです 𝛽0𝑔 ∼ 𝑁 𝜇𝛽0 , 𝜎0 𝛽1𝑔 ∼ 𝑁 𝜇𝛽1 , 𝜎1 ▌例えばbeta_0に関しては parametersブロック modelブロック 𝛽0𝑔 ∼ 𝑁 𝜇𝛽0 , 𝜎0 の2つのパラメータmu_beta_0とsigma_beta_0を宣言 各グループのbeta_0について，事前分布を 𝛽0𝑔 ∼ 𝑁 𝜇𝛽0 , 𝜎0 に 定 11 ベイズ階層モデリング(1) 34

階層ベイズモデルのplate notation 𝜇𝛽0 𝜎𝛽0 𝛽𝑔0 𝑔 group (REGION) 𝜇𝛽DIST 𝛽DIST,𝑔 𝑔𝑖 F OOR𝑔𝑖 𝑔𝑖 M 𝑔𝑖 𝛽FLOOR,𝑔 𝜇𝛽FLOOR 𝜎𝛽FLOOR 𝑖 data 𝛽ITEMS,𝑔 𝜎𝑒 𝜎𝛽DIST 𝜇𝛽ITEMS 𝜎𝛽ITEMS 11 ベイズ階層モデリング(1) あるパラメータの 事前分布のパラメータのことを ハイパーパラメータ（超母数） と呼んだりすることもあります 次ページに stanコードの一例が 載っています 余裕があれば 少し考えてみましょう 35

階層ベイズモデル
▌各係数の事前分布を階層的に

定

仮説モデル 𝐻2 : 回帰係数には地域差があるが，その係数は
全に無秩序というわけでもない，という感じ？

model_multilevel_H2.stan
data {
int N;
int G;
array[N] real SALES;
vector[N] DIST;
vector[N] FLOOR;
vector[N] ITEMS;
array[N] int REGION;
}
parameters {
vector[G] beta_0;
vector[G] beta_DIST;
vector[G] beta_FLOOR;
vector[G] beta_ITEMS;
real <lower=0> sigma;
real mu_beta_0;
real mu_beta_DIST;
real mu_beta_FLOOR;
real mu_beta_ITEMS;
real <lower=0> sigma_beta_0;
real <lower=0> sigma_beta_DIST;
real <lower=0> sigma_beta_FLOOR;
real <lower=0> sigma_beta_ITEMS;
}

transformed parameters {
vector[N] yhat;
for(i in 1:N){
yhat[i] = beta_0[REGION[i]] + beta_DIST[REGION[i]]*DIST[i] +
beta_FLOOR[REGION[i]]*FLOOR[i] + beta_ITEMS[REGION[i]]*ITEMS[i];
}
}
model {
beta_0 ~ normal(mu_beta_0, sigma_beta_0);
beta_DIST ~ normal(mu_beta_DIST, sigma_beta_DIST);
beta_FLOOR ~ normal(mu_beta_FLOOR, sigma_beta_FLOOR);
beta_ITEMS ~ normal(mu_beta_ITEMS, sigma_beta_ITEMS);
sigma ~ cauchy(0, 10);
mu_beta_0 ~ normal(0, 100);
mu_beta_DIST ~ normal(0, 100);
mu_beta_FLOOR ~ normal(0, 100);
mu_beta_ITEMS ~ normal(0, 100);
sigma_beta_0 ~ cauchy(0, 10);
sigma_beta_DIST ~ cauchy(0, 10);
sigma_beta_FLOOR ~ cauchy(0, 10);
sigma_beta_ITEMS ~ cauchy(0, 10);
SALES ~ normal(yhat, sigma);
}
generated quantities {
array[N] real SALES_pred;
SALES_pred = normal_rng(yhat, sigma);
}

11 ベイズ階層モデリング(1)

ブリッジサンプリングをしたい場合などは
target記法で書きましょう

36

前ページの補足 ▌一つ一つ見ていけばわかると思いますが… parametersブロック 𝛽0𝑔 ∼ 𝑁 𝜇𝛽0 , 𝜎0 の2つのパラメータmu_beta_0とsigma_beta_0を宣言 // beta_0に関するところだけ取り出すと parameters { real mu_beta_0; real <lower=0> sigma_beta_0; } modelブロック これを同様に 各係数パラメータについても宣言 各グループのbeta_0について，事前分布を 𝛽0𝑔 ∼ 𝑁 𝜇𝛽0 , 𝜎0 に // beta_0に関するところだけ取り出すと model { beta_0 ~ normal(mu_beta_0, sigma_beta_0); mu_beta_0 ~ normal(0, 100); sigma_beta_0 ~ cauchy(0, 10); } 11 ベイズ階層モデリング(1) 定 これを同様に 各係数パラメータについても宣言 37

推定しましょう ▌推定 # stan_dataはp.20と同じ stan_data <- list(N=100, G=10, SALES = dat$sales, DIST=dat$dist, FLOOR=dat$floor, ITEMS = dat$items, REGION = dat$region) # 階層ベイズモデル model_H2 <- cmdstan_model("model_multilevel_H2.stan") result_H2 <- model_H2$sample(data = stan_data) たぶん"transitions ended with a divergence"や "transitions hit the maximum treedepth limit of 10."などのwarningが出ると思います。 階層モデルはなかなか推定が難しいのです 対処法はこのあとすぐ！ ▌（事後予測チェック用に）MCMCサンプルを取り出す SALES_pred_H2 <- result_H2$draws("SALES_pred", format="matrix")[seq(1,4000,by=500),] 11 ベイズ階層モデリング(1) 38

サンプリングがうまくいかない… ▌階層モデルくらいになってくると推定が厳しくなってくる 【主な理由】パラメータのスケールがそれぞれ異なるため 例えばmu_beta_0とmu_beta_ITEMSの散布図は とはいえ頻 主義ではもはやうまく推定量を 作れなかったりしてお手上げになったりする ので，ベイズ推定ができるだけマシとも言え ます。 X座標とY座標の スケールを 合わせると mu_beta_ITEMSのスケールが 他のパラメータに比べてかなり小さい 11 ベイズ階層モデリング(1) 39

サンプリングがうまくいかない… ▌パラメータのスケールが違うと何が良くないのか？ 適切なステップサイズを決めにくくなる HMC法における適切なステップサイズ＆ステップ数は 流しそうめんの（縦に割った）竹あるいは ハーフパイプのような形をイメージしてください 効率的に尤度関数全体を往来できるように 定したい 関数が のような形だとすると… • 𝜇𝛽0 の方向で十分に移動できるようにすると 𝜇𝛽ITEMS の方向については動き過ぎになる • 𝜇𝛽ITEMS の方向で十分に移動できるようにすると 𝜇𝛽0 の方向については動かなさすぎになる 結果的に自己相関が高くなったり Divergent transitionが発生してしまう 11 ベイズ階層モデリング(1) 40

効率的なサンプリングのために ▌ではどうしたらよいのか？ stanのマニュアルも参照してください 【方法1】パラメータのスケールが近づくように，説明変数のスケールを変更しておく 例えばgをkgに変えるように，単位を変えてしまえば良い 今回のデータはこんな感じなので • floorを100で割る • itemsを1000で割る と，各説明変数のスケールが大体揃いそう！ stan_data <- list(N=100, G=10, SALES=dat$sales, DIST=dat$dist, FLOOR=dat$floor/100, ITEMS=dat$item/1000, REGION = dat$region) # 後は同じように推定するだけ もちろん，得られる回帰係数の値は変化し，解釈も変更後の単位ベースになるので気をつけてください。 11 ベイズ階層モデリング(1) 41

• https://mc-stan.org/docs/stan-users-guide/efficiency-tuning.html#standardizing-predictors

効率的なサンプリングのために ▌ではどうしたらよいのか？ stanのマニュアルも参照してください 【方法2】サンプリング自体は標準正規分布から行い，後からtransformする モデルの書き換えが必要です // beta_0に関するところだけ取り出すと parameters { vector[G] beta_0_raw; real mu_beta_0; real <lower=0> sigma_beta_0; } model { beta_0_raw ~ std_normal(); // beta_0_raw ~ normal(0,1)と同じ // そしてbeta_0 ~ normal(mu_beta_0, sigma_beta_0)と同じ transformed parameters { vector[G] beta_0; beta_0 = mu_beta_0 + beta_0_raw*sigma_beta_0; } mu_beta_0 ~ normal(0, 100); sigma_beta_0 ~ cauchy(0, 10); } この「サンプリング自体は標準正規分布から行う」は よく使われるモデル改善法だと思うので，覚えておいて損はないと思います 11 ベイズ階層モデリング(1) 42

• https://mc-stan.org/docs/stan-users-guide/efficiency-tuning.html#non-centered-parameterization

（おまけ）方法1,2の組み合わせ+αで効率化したコード
data {
int N;
効率化のため，説明変数は行列形式で与えています
int G;
資料09 pp.27-28を参照
int K; //説明変数の数
array[N] real SALES;
matrix[N,K] X; // 説明変数行列
array[N] int REGION;
}
transformed data {
matrix[N, K+1] X_scaled; //標準化した説明変数行列
vector[K] sd_X; //各説明変数の標準偏差を入れておくハコ
X_scaled[,1] = rep_vector(1,N); //切片項をここでくっつける
for(k in 1:K){
sd_X[k] = sd(X[,k]); //標準偏差を求める
X_scaled[,k+1] = X[,k] / sd_X[k]; // 標準偏差で割る
}
}
stanの中で各説明変数のSDを1に変換して推定
parameters {
matrix[G, K+1] beta_raw;
real <lower=0> sigma;

vector[K+1] mu_beta;
vector<lower=0> [K+1] sigma_beta;
}
transformed parameters {
matrix[G, K+1] beta_scaled;
前ページで説明したtransform
for(k in 1:(K+1)){
beta_scaled[,k] = mu_beta[k] + beta_raw[,k] * sigma_beta[k];
}
}

model {
vector[N] yhat = rows_dot_product(X_scaled, beta_scaled[REGION,]);
to_vector(beta_raw) ~ std_normal();
sigma ~ cauchy(0, 10);
mu_beta ~ normal(0, 10);
sigma_beta ~ cauchy(0, 10);
SALES ~ normal(yhat, sigma);
}
generated quantities {
matrix[G, K+1] beta;
beta[,1] = beta_scaled[,1]; // 切片項はそのまま
for(k in 1:K){
beta[,k+1] = beta_scaled[,k+1] / sd_X[k];
}
元の説明変数のスケールでの回帰係数を復元している
}

# スケール変更や切片(1)の追加はstan内で行っている
stan_data <- list(N=100, G=10, K=3, SALES = dat$sales,
X=dat[,c("dist","floor","items")],
REGION = dat$region)

11 ベイズ階層モデリング(1)

手元の環境では，元のコードより
10倍くらい速くなりました

43

結果の解釈 result_H2$summary() ▼ p. 41に沿ってデータ変換後の結果｜ハイパーパラメータの部分だけ抜粋して表示 ▌ EAP(mean)ベースで考えてみると 𝛽DIST ∼ 𝑁 −0.15, 0.0494 𝛽𝐹𝐿𝑂𝑂𝑅 ∼ 𝑁 5.84, 0.0457 𝛽𝐼𝑇𝐸𝑀𝑆 ∼ 𝑁(−1.13, 0.0298) 地域ごとに違いはあれど，全体的な傾向としては • 駅からの距離（dist）が長いほどsalesは小さい • 床面積（floor）が大きいほどsalesは大きい • 取り扱い商品数（items）が多いほどsalesは小さい 11 ベイズ階層モデリング(1) 44

結果の比較 ▌試しにbeta_DISTについてH1とH2の結果を比べてみると mcmc_areas(result_H1$draws("beta_DIST",format="df")) 次ページのH2の結果と交互に べてください 11 ベイズ階層モデリング(1) 45

結果の比較 ▌試しにbeta_DISTについてH1とH2の結果を比べてみると mcmc_areas(result_H2$draws("beta_DIST",format="df")) 前ページのH1の結果と交互に べてください 階層モデルのほうが • 事後分布の幅が狭くなっている ▲ 階層事前分布によって他の地域の情報が追加されたため • 各事後分布がmu_beta_DISTのEAPに寄っている ▲ これはモデルの仕様上あたりまえの話 shrinkage が起こっている，とも言えます 11 ベイズ階層モデリング(1) 46

事後予測チェックもしてみよう もちろんやり方は色々と考えられます（ベイズファクターや情報量規準を使うのも当然OK） ▌とりあえずデータの密度を比較してみると ppc_dens_overlay(dat$sales, SALES_pred_H1) ppc_dens_overlay(dat$sales, SALES_pred_H2) 明確に良くなった，という感じには見えない…？ 11 ベイズ階層モデリング(1) 47

事後予測チェックもしてみよう（続き） ▌グループごとに異なるかも，と考えているのだから プロットもグループごとにやってもいいのでは？ ppc_dens_overlay_grouped(dat$sales, SALES_pred_H2, group = dat$region) p. 23と行ったり来たりしながら 見比べてみてください やはり明らかに良くなっているようには見えない…？ 11 ベイズ階層モデリング(1) 48

まとめと次回予告 【まとめ】 ▌ベイズ階層モデルの導入を行いました 通常のマルチレベルモデルと同様に「係数の分布」を考える ベイズ統計ではそれが事前分布の役割を担うことで他のグループの情報を取り込める ただし推定は結構たいへんになってくる そして今回のデータの場合特に大きな恩恵は見られなかった 【次回予告】 ▌ベイズ階層モデルについてもう少し深掘りしていきます なぜ今回はさほど改善が見られなかったのか？ そうなると，わざわざ階層モデルにするメリットはどこに現れるのか？ 11 ベイズ階層モデリング(1) 49