ベイズ統計_12_ベイズ階層モデリング2

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

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

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

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

（前回のおさらい）はじめにやったこと ▌グループごとに異なる回帰係数を推定する 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回繰り返しただけ 12 ベイズ階層モデリング(2) 4

（前回のおさらい）ベイズ階層モデリング ▌本当にやりたいこと 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以外から得られる情報 を組み込めたら最高 12 ベイズ階層モデリング(2) もちろん他の地域でも同じ 5

（前回のおさらい）いざやってみると ▌事後予測チェック（データの密度）を比較してみると グループごとに独立に推定した場合 階層ベイズモデルにした場合 ppc_dens_overlay(dat$sales, SALES_pred_H1) ppc_dens_overlay(dat$sales, SALES_pred_H2) 階層事前分布を導入したところで 明確に良くなった，という感じには見えない…？ 12 ベイズ階層モデリング(2) 6

１ 階層ベイズモデルの利点 12 ベイズ階層モデリング(2) 7

なぜあまり差が出なかったのか？ ▌10店舗あれば，ある程度の推定はできてしまうため 1 地域 店舗 1 2 … 10 beta_DIST[1]の事後分布 cf. 階層ベイズモデルの場合 EAP 𝛽0 1.3117 𝛽DIST -0.1917 𝛽FLOOR 0.0554 𝛽ITEMS -0.0011 多少のズレはあれど，事後分布の幅はさほど変わらない と言えるくらいにはデータも情報を持っている，ということです 12 ベイズ階層モデリング(2) 8

利点①｜推定の安定性 ▌SALES_predの事後分布のSD 階層ベイズモデル SALES_predに限らず，階層性をもつパラメータは いずれも同じような傾向が見られます 𝑦=𝑥 階層事前分布を通して 情報が与えられる パラメータの事後分布は よりシャープな形になることが多い ただし真値に近づけるか（＝精度が向上するか）は 場合によります。 （例：あるグループが他のグループとは特異的なパラ メータを持つ場合，むしろずれてしまう可能性もある） 階層ベイズモデルではない 12 ベイズ階層モデリング(2) 9

利点①｜推定の安定性（別の例） あるいは単純に10店舗もない地域，とか 地域３の担当者がサボったため，ここだけ10店舗ではなく 3店舗しかデータが得られませんでした。 試しにsalesとdistの散布図を見てみると… 問題を簡単にするため これ以降は説明変数をdistのみとした 単回帰分析モデルを考えます ▲単回帰直線 この地域では，駅から遠いほど 売上が伸びるようです 12 ベイズ階層モデリング(2) 10

利点①｜推定の安定性（別の例） あるいは単純に10店舗もない地域，とか 地域３の担当者がサボったため，ここだけ10店舗ではなく 3店舗しかデータが得られませんでした。 実際に階層ベイズではないモデルで回帰分析をすると… 問題を簡単にするため これ以降は説明変数をdistのみとした 単回帰分析モデルを考えます ３店舗のみ 10店舗全て beta_DISTのEAP beta_DISTのEAP 0.776 -0.199 点線が真値かどうかは また別の話です データがとても少ない場合 回帰係数は大きく歪んでしまうことがある！ 12 ベイズ階層モデリング(2) 11

利点①｜推定の安定性（別の例） あるいは単純に10店舗もない地域，とか 地域３の担当者がサボったため，ここだけ10店舗ではなく 3店舗しかデータが得られませんでした。 階層ベイズモデルで回帰分析をすると… 問題を簡単にするため これ以降は説明変数をdistのみとした 単回帰分析モデルを考えます 3店舗のみ｜not階層 beta_DISTのEAP 0.776 10店舗全て beta_DISTのEAP -0.199 beta_DISTのEAP -0.093 3店舗のみ｜階層 階層ベイズモデルならば データが少ない場合でもある程度安定する 別の見方をすると「過学習の抑制」 12 ベイズ階層モデリング(2) 12

単回帰分析のnot 階層ベイズモデルのstanコード
▌distと切片のみの単回帰分析モデル
model_grouped_dist.stan
data {
int N;
int G;
array[N] real SALES;
vector[N] DIST;
array[N] int REGION;
}
parameters {
vector[G] beta_0;
vector[G] beta_DIST;
real <lower=0> sigma;
}

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

SALES ~ normal(yhat, sigma);

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

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

12 ベイズ階層モデリング(2)

13

単回帰分析の階層ベイズモデルのstanコード
▌distと切片のみの単回帰分析モデル
model_multilevel_dist.stan
data {
int N;
int G;
array[N] real SALES;
vector[N] DIST;
array[N] int REGION;
}
parameters {
vector[G] beta_0;
vector[G] beta_DIST;
real <lower=0> sigma;

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

real mu_beta_0;
real mu_beta_DIST;
real <lower=0> sigma_beta_0;
real <lower=0> sigma_beta_DIST;

}

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

12 ベイズ階層モデリング(2)

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

14

ベイズ階層モデルのちから ▌２つのモデルで得られるbeta_DIST[3]の事後分布 【階層モデルのないとき】 EAP 0.776 【階層モデルのあるとき】 EAP -0.093 3店舗のデータから推定するしかないので 他の地域での傾きを考慮すると 回帰係数は正のほうに寄ってしまう この地域も大体同じようなものだろう （そして情報が少ないのであまり自信はない） （データとはあまり合わないがやむなし） 12 ベイズ階層モデリング(2) 15

メリット②｜新しいデータの予測 「未開の地」だったあの地域に進出することにしたぞい 駅からの距離は近いほうが良いかね？ （実は駅から2kmのいい土地が買えそうなんだよ） ルチレベルモデルの定式 （ 考）グループで差がない 【階層モデルを知らない人】 の回帰分析 = 地域によるから，まだわからないですねェ グループごとに切片が異なる（ラン （レベル ） 切片モデル） = 資料11 p. = 8 = （レベル ） グループごとに傾きも異なる（ラン （レベル ） 切片・傾きモデル） = = （レベル ） = = ▲ ここが全体的な傾向 ベイズ階層モデリング 【階層モデルを知っている人】 この地域の情報は何もないですが， 他の地域の全体的な傾向から ◯◯だと言えそうです （駅から2kmだと売上の予測は▲▲です） 12 ベイズ階層モデリング(2) 16

予測してみよう ▌未知の地域𝑁𝐸𝑊に所属する店舗の回帰モデルの式は 𝑦𝑁𝐸𝑊 𝑖 = 𝛽0 𝑁𝐸𝑊 𝛽DIST 𝑁𝐸𝑊 DIST𝑁𝐸𝑊 𝑖 𝜀𝑁𝐸𝑊 𝑖 𝛽0 𝑁𝐸𝑊 = 𝜇𝛽0 𝑢0 𝑁𝐸𝑊 𝑢0 𝑁𝐸𝑊 ∼ 𝑁 𝜎𝛽0 𝛽DIST 𝑁𝐸𝑊 = 𝜇𝛽DIST 𝑢DIST 𝑁𝐸𝑊 𝑢DIST 𝑁𝐸𝑊 ∼ 𝑁 𝜎𝛽DIST 𝛽0 𝑁𝐸𝑊 ∼ 𝑁 𝜇𝛽0 𝜎𝛽0 𝛽DIST 𝑁𝐸𝑊 ∼ 𝑁 𝜇𝛽DIST 𝜎𝛽DIST ▌（データのない）未知の地域𝑁𝐸𝑊の回帰係数の事前分布は beta_0_new ~ normal(mu_beta_0, sigma_beta_0); beta_DIST_new ~ normal(mu_beta_DIST, sigma_beta_DIST); ベイズの定理的には，データがない場合 事後分布＝事前分布となる ▌これをもとにgenerated quantitiesブロックで回帰係数を予測する 12 ベイズ階層モデリング(2) 17

未知の地域の係数を予測するstanコード
▌distと切片のみの単回帰分析モデル

p. 14のモデルをベースにした修正

model_multilevel_dist_pred.stan
data {
int N;
int G;
array[N] real SALES;
vector[N] DIST;
array[N] int REGION;
}

transformed parameters {
vector[N] yhat;
for(i in 1:N){
yhat[i] = beta_0[REGION[i]] + beta_DIST[REGION[i]]*DIST[i];
}
}
model {
beta_0 ~ normal(mu_beta_0, sigma_beta_0);
beta_DIST ~ normal(mu_beta_DIST, sigma_beta_DIST);
sigma ~ cauchy(0, 10);

parameters {
vector[G] beta_0;
vector[G] beta_DIST;
real <lower=0> sigma;
real mu_beta_0;
real mu_beta_DIST;
real <lower=0> sigma_beta_0;
real <lower=0> sigma_beta_DIST;

mu_beta_0 ~ normal(0, 100);
mu_beta_DIST ~ normal(0, 100);
sigma_beta_0 ~ cauchy(0, 10);
sigma_beta_DIST ~ cauchy(0, 10);
}

SALES ~ normal(yhat, sigma);

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

}

real beta_0_new;
real beta_DIST_new;
beta_0_new = normal_rng(mu_beta_0, sigma_beta_0);
beta_DIST_new = normal_rng(mu_beta_DIST, sigma_beta_DIST);

}

real SALES_pred_new;
SALES_pred_new = normal_rng(beta_0_new + beta_DIST_new*2, sigma);

「駅から2km」の店舗の予測
12 ベイズ階層モデリング(2)

18

新しい地域の回帰直線の予測 𝛽0 𝑁𝐸𝑊 の 事後分布 各iterationから回帰直線を作成して 薄ーく重ねまくったプロット 𝛽DIST 𝑁𝐸𝑊 の 事後分布 どうやら負の影響である可能性が高いです ただ，距離は僅かな影響しかないでしょう 12 ベイズ階層モデリング(2) 19

新しい地域の新しい店舗のSALESの予測 各iterationから回帰直線を作成して 薄ーく重ねまくったプロット SALESの 事後予測分布 駅から2kmに店舗を作った場合のSALESは だいたいこれくらいになります 12 ベイズ階層モデリング(2) 20

２ 階層ベイズモデルの拡張 あとは皆さんの想像力次第 12 ベイズ階層モデリング(2) 21

拡張①｜レベル2に説明変数をいれる ▌回帰係数の差異の理由を推測したい DISTがSALESに強く影響するのはどんな地域？ 【仮説のロジック 事実とは異なります】 • DISTがSALESに影響するのは電車メインの都会 • だとすると，面積が小さい・人口の多い地域ほど DISTの影響は大きくなるのでは？ もしその理由が多少なりともわかれば p. 17 未開の地域の予測 この地域は面積が◯◯であるため DISTの回帰係数は だいたい■■くらいになりそうです 予測の精度を上げることができる！ 12 ベイズ階層モデリング(2) 22

新しいデータを使って予測してみよう ▌データの読み込み "data_region.csv"をワーキングディレクトリに配置して dat_region <- read.csv("data_region.csv") この名前はなんでもOK region: 地域 population: その地域の人口 area: その地域の面積 ※データは適当に作ったので，実際とは異なります。 12 ベイズ階層モデリング(2) 23

回帰係数と地域の変数の関係を見てみる ▌とりあえずp. 14のモデルで推定した回帰係数について library(cmdstanr) model <- cmdstan_model("model_multilevel_dist.stan") stan_data <- list(N=100, G=10, SALES=dat$sales, DIST=dat$dist, REGION=dat$region) result <- model$sample(data=stan_data) beta_DIST_EAP <- result$summary("beta_DIST")$mean plot(dat_region$area, beta_DIST_EAP) の どうやら 広い地域ほど𝛽DIST の影響は弱まりそう この関係をモデルに加えて より精緻な推定を目指そう！ 地域の面積 地域4の外れ値感に引っ張られているような気もしますが そのあたりの判断はお任せします 12 ベイズ階層モデリング(2) 24

回帰式はどうなる？ ▌もとの階層線形モデルの式 レベル1 レベル ２ 𝑦𝑔𝑖 = 𝛽0𝑔 𝛽DIST 𝑔 DIST𝑔𝑖 𝛽0𝑔 = 𝜇𝛽0 𝑢0𝑔 𝛽DIST 𝑔 = 𝜇𝛽DIST 𝑢DIST 𝑔 𝜀𝑔𝑖 𝑢0𝑔 ∼ 𝑁 【plate notation】 𝜎𝛽0 𝑢DIST 𝑔 ∼ 𝑁 𝜎𝛽DIST 𝛽AREA 𝜇𝛽DIST AREA𝑔 𝛽DIST 𝑔 𝜎𝛽DIST 𝛽DIST 𝑔 ∼ 𝑁 𝜇𝛽DIST 𝜎𝛽DIST ▌レベル2に説明変数を追加すると 𝛽DIST 𝑔 = 𝜇𝛽DIST 𝛽AREA AREA𝑔 𝑢DIST𝑔 ∼ 𝑁 𝜎𝛽DIST 𝛽DIST 𝑔 ∼ 𝑁 𝜇𝛽DIST 𝑢DIST 𝑔 𝛽0𝑔 DIST𝑔𝑖 𝛽AREA AREA𝑔 𝜎𝛽DIST 12 ベイズ階層モデリング(2) 𝜇𝛽0 SALES𝑔𝑖 𝜎𝛽0 𝑖 data 𝑔 group (REGION) 𝜎𝑒 25

stanコードを拡張してみましょう ▌再びp. 14のモデルをベースにして拡張してみます ※次ページにコード例が載っているので，考えたい人は一旦ストップです 【ヒント】 • 前ページのplate notation的には，追加の登場人物は 𝛽𝐴𝑅𝐸𝐴 （パラメータ） 𝐴𝑅𝐸𝐴 （データ） の２つのみ • モデル式的には𝛽DIST 𝑔 ∼ 𝑁 𝜇𝛽DIST 𝜎𝛽DIST が 𝛽DIST 𝑔 ∼ 𝑁 𝜇𝛽DIST 𝛽AREA AREA𝑔 𝜎𝛽DIST に変わっただけ （モデルの式がきちんと理解できていれば コードの修正は割と簡単です） 12 ベイズ階層モデリング(2) 26

レベル2に説明変数を追加したstanコード
▌distと切片のみの単回帰分析モデル

p. 14のモデルをベースにした修正

model_multilevel_dist_area.stan
data {
int N;
int G;
array[N] real SALES;
vector[N] DIST;
array[N] int REGION;
vector[G] AREA;
// modelブロックでforループを使うなら
// array[G] real AREA;
// ↑でも大丈夫
}
parameters {
vector[G] beta_0;
vector[G] beta_DIST;
real <lower=0> sigma;
real beta_AREA;
real mu_beta_0;
real mu_beta_DIST;
real <lower=0> sigma_beta_0;
real <lower=0> sigma_beta_DIST;
}

transformed parameters {
vector[N] yhat;
for(i in 1:N){
yhat[i] = beta_0[REGION[i]] + beta_DIST[REGION[i]]*DIST[i];
}
}
model {
beta_0 ~ normal(mu_beta_0, sigma_beta_0);
beta_DIST ~ normal(mu_beta_DIST + beta_AREA*AREA, sigma_beta_DIST);
// forループで書くと
// for(g in 1:G){
//
beta_DIST[g] ~ normal(mu_beta_DIST + beta_AREA*AREA[g], sigma_beta_DIST);
// }
sigma ~ cauchy(0, 10);
mu_beta_0 ~ normal(0, 100);
mu_beta_DIST ~ normal(0, 100);
beta_AREA ~ normal(0, 100);
sigma_beta_0 ~ cauchy(0, 10);
sigma_beta_DIST ~ cauchy(0, 10);
}

SALES ~ normal(yhat, sigma);

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

12 ベイズ階層モデリング(2)

27

推定＆結果を見てみる model_area <- cmdstan_model("model_multilevel_dist_area.stan") stan_data <- list(N=100, G=10, SALES=dat$sales, DIST=dat$dist, REGION=dat$region, AREA=dat$area) result_area <- model_area$sample(data=stan_data) result_area$summary("beta_AREA") mcmc_dens(result_area$draws("beta_AREA")) 結構複雑なモデルなので実は推定はシビアで あまりきれいにサンプリングできているとは言えない （とりあえず今は無視して進めます。） （余裕がある人は，どうやったらうまくサンプリング できるか考えてみてください。） 広い地域ほど 距離の影響は弱まりそうです （車社会だから？） 12 ベイズ階層モデリング(2) 28

AREAを追加した影響はほかにも 【AREAなしモデルでの結果】 ▌beta_DISTの事後分布を比べてみる 次ページと行き来して比較してみてください 12 ベイズ階層モデリング(2) 29

AREAを追加した影響はほかにも 【AREAありモデルでの結果】 ▌beta_DISTの事後分布を比べてみる 前ページと行き来して比較してみてください AREAなしモデルのほうが mu_beta_DISTのほう（EAP）へ寄っている AREAありモデルのほうは𝛽AREA AREA𝑔 の項が 追加されたことで𝜇𝛽DIST へのshrinkが抑制された？ ただそれが真値に近づいているかは別の話 もちろん「適切な説明変数を入れたら」の話です 12 ベイズ階層モデリング(2) 30

AREAを追加した影響はほかにも ▌SALES_predの事後予測分布のSDを比べてみる 𝑦=𝑥 あり AREAありモデルのほうが beta_DISTの推定に関する情報が増えた 推定がより精緻になったと言える その結果，事後予測分布についても AREAありモデルのほうが 事後SDは小さい傾向が見られる もちろん「適切な説明変数を入れたら」の話です なし 12 ベイズ階層モデリング(2) 31

拡張②｜回帰係数間の相関関係を考える ▌説明変数の影響や切片との間には関係があることも 【例】コンビニの売上について FLOORの影響が大きい地域のほうが DISTの影響も大きい，とかないかね？ Cor 𝛽FLOOR 𝛽DIST 【例】勉強時間と模試の成績について 基礎学力が高い学校ほど 勉強時間あたりの成績の上昇が大きい，とかは？ Cor 𝛽0 𝛽1 もし回帰係数の間に相関があれば，これも情報として利用できる！ • 推定や予測の精度があがる 例えばFLOORとDISTの影響はいずれも • 背後に共 の原因があるかも，と話がふくらむ その地域の都会度によって変わっているかも 12 ベイズ階層モデリング(2) 32

回帰係数間の相関をモデルに加える ▌distのみの ルチレベルモデル（p. 14のコード） レベル1 𝑦𝑔𝑖 = 𝛽0𝑔 𝛽DIST 𝑔 DIST𝑔𝑖 レベル 𝛽0𝑔 = 𝜇𝛽0 𝑢0𝑔 ２ 𝛽DIST 𝑔 = 𝜇𝛽DIST 𝑢DIST 𝑔 𝜀𝑔𝑖 𝑢0𝑔 ∼ 𝑁 𝛽0𝑔 ∼ 𝑁 𝜇𝛽0 𝜎𝛽0 𝜎𝛽0 𝑢DIST 𝑔 ∼ 𝑁 𝜎𝛽DIST 𝛽DIST 𝑔 ∼ 𝑁 𝜇𝛽DIST 𝜎𝛽DIST 2つの係数は独立に正規分布に従うと仮定していた ▌2つの係数に相関関係を許容すると レベル1 𝑦𝑔𝑖 = 𝛽0𝑔 𝛽DIST 𝑔 DIST𝑔𝑖 レベル 𝛽0𝑔 = 𝜇𝛽0 𝑢0𝑔 ２ 𝛽DIST 𝑔 = 𝜇𝛽DIST 𝑢DIST 𝑔 𝐮𝑔 = 𝑢 𝑢0𝑔 DIST 𝑔 ∼𝑁 𝜷𝑔 = 𝜎𝛽20 係数の共分散を Cov で表すことにします Cov Cov 𝜎𝛽2DIST 𝛽0𝑔 𝛽DIST 𝑔 ∼𝑁 𝜇𝛽0 𝜇𝛽DIST 𝜎𝛽20 Cov Cov 𝜎𝛽2DIST 2つの係数が多変量正規分布に従うと仮定したらよい！ 12 ベイズ階層モデリング(2) 33

stanで多変量正規分布を扱う場合 ▌相関のあるパラメータを1つのベクトルとして用意する // parametersブロック vector[G] beta_0; vector[G] beta_DIST; // parametersブロック array[G] vector[2] beta; beta array beta[g] vector beta[g][1]=𝛽0𝑔 beta[g][2]=𝛽DIST 𝑔 ▌平均・共分散行列もそれぞれ1つにまとめて用意すると良い // parametersブロック real mu_beta_0; real mu_beta_DIST; real <lower=0> sigma_beta_0; real <lower=0> sigma_beta_DIST; // parametersブロック vector[2] mu_beta; cov_matrix[2] Sigma_beta; 12 ベイズ階層モデリング(2) cov_matrix型で指定すると 共分散行列の性質を満たさない サンプリングが行われなくなる （ので，<lower=0>などは不要） 34

stanで多変量正規分布を扱う場合（続き） ▌modelブロックもそれに合わせて 【事前分布】 mu_beta_0 ~ normal(0, 100); mu_beta_DIST ~ normal(0, 100); sigma_beta_0 ~ cauchy(0, 10); sigma_beta_DIST ~ cauchy(0, 10); mu_beta ~ normal(0, 100); Sigma_beta ~ inv_wishart(3, diag_matrix(rep_vector(1,2))); 【尤度に関するところ】 transformed parameters { vector[N] yhat; for(i in 1:N){ yhat[i] = beta_0[REGION[i]] + beta_DIST[REGION[i]]*DIST[i]; } } // 以下はmodelブロック beta_0 ~ normal(mu_beta_0, sigma_beta_0); beta_DIST ~ normal(mu_beta_DIST, sigma_beta_DIST); transformed parameters { vector[N] yhat; for(i in 1:N){ yhat[i] = beta[REGION[i]][1] + beta[REGION[i]][2]*DIST[i]; } } // 以下はmodelブロック beta ~ multi_normal(mu_beta, Sigma_beta); 12 ベイズ階層モデリング(2) 35

（補足）inv_wishartってなんですか mu_beta ~ normal(0, 100); Sigma_beta ~ inv_wishart(3, identity_matrix(2)); ▌ 多変量正規分布の共分散行列の共役事前分布として知られている 逆ガンマ分布（カイ二乗分布）を多変量に一般化したもの ▌ きちんと共分散行列の性質を満たすような空間からサンプリングしてくれる ↑の理由と合わせて，共分散行列の事前分布としてよく用いられてきた ただ逆ガン 分布のときと同様に批判もあったりします そしてサンプリングの効率の観点からは 逆ウィシャート分布は使わない方法のほうが良かったりします ▌ 2つのパラメータの設定 1つめを「パラメータ数+1」，2つ目を単位行列にする事が多い 相関（非対角成分）の周辺分布が一様分布＝無情報になるため identity_matrix(2)で2×2の単位行列を作っている 12 ベイズ階層モデリング(2) 36

回帰係数間の相関を許容したstanコード
▌distと切片のみの単回帰分析モデル

p. 14のモデルをベースにした修正

model_multilevel_dist_cor.stan
data {
int N;
int G;
array[N] real SALES;
vector[N] DIST;
array[N] int REGION;
}
parameters {
array[G] vector[2] beta;
real <lower=0> sigma;
vector[2] mu_beta;
cov_matrix[2] Sigma_beta;
}

transformed parameters {
vector[N] yhat;
for(i in 1:N){
yhat[i] = beta[REGION[i]][1] + beta[REGION[i]][2]*DIST[i];
}
}
model {
beta ~ multi_normal(mu_beta, Sigma_beta);
mu_beta ~ normal(0, 100);
Sigma_beta ~ inv_wishart(3, identity_matrix(2));
sigma ~ cauchy(0, 10);
}

SALES ~ normal(yhat, sigma);

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

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

▌あとはいつもどおり推定してみてください
12 ベイズ階層モデリング(2)

37

結果を見てみる ▌mu_beta ▌Sigma_beta EAPで見るならば 𝜷𝑔 = 𝛽0𝑔 𝛽DIST 𝑔 ∼𝑁 7.53 − . 68 .3 2 − . 3 6 − . 3 6 . 37 12 ベイズ階層モデリング(2) −0.0306 𝛽0 𝛽DIST 間の相関は 0.312 0.137 ≃ − . 48 38

（おまけ）多変量正規分布のサンプリングを効率よくするには
▌共分散行列の代わりに「SD（ベクトル）」と「相関行列」をサンプリングする
𝚺=

𝜎𝛽20

Cov

Cov 𝜎𝛽2DIST

=

𝜎𝛽20

Cor × 𝜎𝛽0 𝜎𝛽DIST

Cor × 𝜎𝛽0 𝜎𝛽DIST

𝜎𝛽2DIST

=

// この修正に関係ある箇所だけ抜粋
parameters {
corr_matrix[2] R_beta; // 相関行列
vector<lower=0>[2] sigma_beta; // SDベクトル
}
transformed parameters {
matrix[2,2] Sigma_beta = quad_form_diag(sigma_beta, R_beta);
}

model {
R_beta ~ lkj_corr(2);
sigma_beta ~ cauchy(0, 10);
}

𝜎𝛽0

Cor 𝜎𝛽0

𝜎𝛽DIST Cor

𝜎𝛽DIST

= 𝐒′𝐑𝐒

相関行列𝐑
【transformed parametersについて】
• cov_matrix[2]で宣言してもよいが，必ず共分散行列になるので
無駄なチェックを避けるため普 のmatrixとして指定している
• quad_form_diag(matrix, vector)は
1. vectorを対角成分にもつ対角行列𝐒を作る
2. それをmatrix(𝐑)と合わせて𝐒′𝐑𝐒を計算する
▲を効率よくやってくれる関数です
【modelについて】
• lkj_corr()は相関行列の事前分布として推奨されているもの
• lkj_corr(2)は弱情報，lkj_corr(1)は無情報

12 ベイズ階層モデリング(2)

39

（おまけ）多変量正規分布のサンプリングを効率よくしたstanコード
▌全体

p. 38のモデルをベースにした修正

data {
int N;
int G;
array[N] real SALES;
vector[N] DIST;
array[N] int REGION;
}
parameters {
array[G] vector[2] beta;
real <lower=0> sigma;
vector[2] mu_beta;
corr_matrix[2] R_beta;
vector<lower=0>[2] sigma_beta;
}

transformed parameters {
matrix[2,2] Sigma_beta = quad_form_diag(R_beta, sigma_beta);
vector[N] yhat;
for(i in 1:N){
yhat[i] = beta[REGION[i]][1] + beta[REGION[i]][2]*DIST[i];
}
}
model {
beta ~ multi_normal(mu_beta, Sigma_beta);
mu_beta ~ normal(0, 100);
R_beta ~ lkj_corr(2);
sigma_beta ~ cauchy(0, 10);
sigma ~ cauchy(0, 10);
}

SALES ~ normal(yhat, sigma);

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

さらに効率 するためには，
• 相関行列をコレスキー分解する
• betaをそれぞれ独立に（一変量）正規分布からサンプリングし，
といったテクニックが有効です
12 ベイズ階層モデリング(2)

を使って多変量正規分布に従うように変換する

40

扱わなかったけれども重要な話 ▌回帰係数にはグループ効果と個体の効果が混ざっている可能性 𝛽DIST 𝑔 がマイナスの場合「駅から遠い店舗ほど売上が下がる」 ただし個体の「駅からの距離」は「グループ平均」と「グループ平均からの差」の和 回帰式のうち「駅からの距離」に関する部分は グループ平均からの差 𝛽DIST 𝑔 DIST𝑔𝑖 = 𝛽DIST 𝑔 DIST𝑔𝑖 − DIST𝑔 グループ平均 𝛽DIST 𝑔 DIST𝑔 𝛽DIST 𝑔 はこの両方を考慮して推定された値になっている • 𝛽DIST 𝑔 DIST𝑔𝑖 − DIST𝑔 について𝛽DIST 𝑔 がマイナスとは グループ内で相対的に遠い店舗ほど売上が下がる • 𝛽DIST 𝑔 DIST𝑔 について𝛽DIST 𝑔 がマイナスとは 平均的に駅から遠いグループほど売上が下がる 12 ベイズ階層モデリング(2) 41

扱わなかったけれども重要な話｜中心 の必要性 ▌式に従って分けてあげたらよいのでは？ グループ平均からの差 グループ平均 𝛽DIST 𝑔 DIST𝑔𝑖 = 𝛽DIST 𝑔 DIST𝑔𝑖 − DIST𝑔 店舗レベルの効果を 表しているので レベル1の説明変数 𝛽DIST 𝑔 DIST𝑔 グループレベルで 中心 した値を グループレベルの効果を 表しているので レベル2の説明変数 ∗ DIST𝑔𝑖 とする ▌すると全体のモデル式は レベル1 𝑦𝑔𝑖 = 𝛽0 𝑁𝐸𝑊 レベル 𝛽0𝑔 = 𝜇𝛽0 ２ 𝛽DIST 𝑔 = 𝜇𝛽DIST ∗ 𝛽DIST 𝑔 DIST𝑔𝑖 𝜀𝑔𝑖 𝑢0𝑔 𝛽DIST 𝑔 DIST𝑔 𝑢0𝑔 ∼ 𝑁 𝑢DIST 𝑔 𝑢DIST 𝑔 ∼ 𝑁 グループレベルの効果を分離するためには 説明変数の中心 が必要です 12 ベイズ階層モデリング(2) 𝜎𝛽0 𝜎𝛽DIST このあたりの話は，ベイズ統計とは別に ルチレベルモデルについて勉強してください。 42

まとめと次回予告 【まとめ】 ▌ベイズ階層モデルについて深掘りしました どのような場合にメリットがあらわれるのか もっと複雑な問題に対応できるようにするためには，自由に拡張したらよい （ただしうまく推定できる範囲内で） 【次回予告】 ▌別のモデルの実装例を紹介します たぶん離散パラメータを使うモデル（混合分布・ゼロ過剰分布）を紹介します ※ で実行する場合少しテクニックがいるのですが，知っておくと表現の幅が広がると思います。 （GLMの話もするかも？） 12 ベイズ階層モデリング(2) 43