ベイズ統計_09_モデル評価・仮説検証

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

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

前回までのまとめ 【前回までのまとめ】 ▌ ベイズ統計学とは何かを「完全に理解した」（01-05） ベイズ統計学における「確率」の考え方とその更新 いくつかのシンプルな例について，「分布の（パラメータの）更新」の意味を確かめた ▌ ベイズ統計におけるサンプリングを「完全に理解した」(06-08) マルコフ連鎖モンテカルロ法を使えば，たいていどんな事後分布も近似できる stanでは，いろいろな設定を変えてあげることでたいてい目的が達成できる 【今回からは】 ▌ 得られた結果を使ってリアルな課題に立ち向かうには？を考えていきます 事後分布をどのように扱えば，統計モデルの評価や仮説の検証ができるのだろうか？ もう少し複雑なモデルをstanで実現するためのチュートリアル（続き） 09 モデル評価・仮説検証 2

１ モデル評価 09 モデル評価・仮説検証 3

（社会科学および）統計学がやりたいこと ▌世界の姿を簡略化して，世の中の現象を明らかにしたい 例 あるコンビニの売上は… 現実世界 人類が観察できない要因も含めた 近隣住民によって変わる 様々な条件によって変わる オーナーの方針 気候によって変わる どうしたら 売上が増えるの？ によって変わる 周辺の競合 によって変わる 景気によって変わる 09 モデル評価・仮説検証 ケースバイケース でしょうねぇ… 4

（社会科学および）統計学がやりたいこと ▌世界の姿を簡略化して，世の中の現象を明らかにしたい 例 あるコンビニの売上は… （線形）回帰分析モデル 複雑な現象を簡略化して 近隣住民 によって変わる • コンビニの売上に強く影響するもののみ考えてあとは誤差とみなす • それらが売上に及ぼす影響はすべて独立である 気候によって変わる オーナーの方針 𝑦𝑖 = 𝛽0 + 𝛽1 𝑥𝑖1 + 𝛽2 𝑥𝑖2 + 𝜀𝑖 …と仮定すると によって変わる どうしたら 売上が増えるの？ 周辺の競合 によって変わる それを使って 推論や意思決定ができるような 景気によって変わる 簡略化したモデルを考えたい 09 モデル評価・仮説検証 駅チカの物件を おさえましょう！ 5

そうなると重要なのが… ▌つくった（統計）モデルは，現実世界にどれだけ近づけているか？ 現実世界 データ 神でもないと 観察はできない 近隣住民によって変わる オーナーの方針 気候によって変わる によって変わる 周辺の競合 によって変わる 景気によって変わる 関数と事後分布 の組み合わせ 本当はこの間が どれだけ近づいているかを 評価できたら良いが… 𝑦𝑖 = 𝛽0 + 𝛽1 𝑥𝑖1 + 𝛽2 𝑥𝑖2 + 𝜀𝑖 （データで更新した）統計モデル 09 モデル評価・仮説検証 6

できる範囲で評価する ▌つくった（統計）モデルは，現実世界にどれだけ近づけているか？ 現実世界 データ 近隣住民によって変わる オーナーの方針 気候によって変わる によって変わる 周辺の競合 によって変わる この間であれば どちらも観測可能なため 直接比較することができる！ 景気によって変わる 関数と事後分布 の組み合わせ 𝑦𝑖 = 𝛽0 + 𝛽1 𝑥𝑖1 + 𝛽2 𝑥𝑖2 + 𝜀𝑖 統計モデルから作ったデータ 09 モデル評価・仮説検証 （データで更新した）統計モデル 7

ということで回帰分析を例に 例 あるコンビニチェーンのアナリストは，各店舗の利益に影響する要因を調べることにしました。 要因は色々考えられますが，とりあえずデータの中にある3つの変数が，利益にどのように 影響するかを確かめてみることにしました。 【データの読み込み】 この名前はなんでもOK “data_cvs.csv”をワーキングディレクトリに配置して ▌ 中身はこんなデータ dat <- read.csv("data_cvs.csv") 今日はここまで使います sales: その店の一日あたり平均利益（単位：千円） 被説明変数 dist: 最寄り駅からの距離（単位：km） floor: 床面積（単位：m2） 説明変数 items: 取扱いアイテム数（単位：個） region: その店舗のある地域 neighbor: 半径1km以内にあるコンビニの数 ※データは適当に作ったので、実際とは異なります。 09 モデル評価・仮説検証 8

回帰分析ってどんなんだったっけ 忘れたという人はこちらを参照 𝑦𝑖 = 𝛽0 + 𝛽1 𝑥𝑖1 + 𝛽2 𝑥𝑖2 + ⋯ + 𝑒𝑖 ▶ 「大体」の予測値 + 予測とのズレ 𝑦 𝛽0 （切片） 𝑥𝑘 の値がすべて0のとき 𝑦 の予測値（平均的な傾向）は いくつになるか 𝛽𝑘 （傾き） 𝑥𝑘 の値が1大きくなると 𝑦の予測値（平均的な傾向）は どれだけ大きくなるか 𝑥1 09 モデル評価・仮説検証 9

• https://www2.kobe-u.ac.jp/~bunji/files/lecture/stat/stat-04-regression.pdf

回帰分析の統計モデル ▌推定・検定の際に必要だったいくつかの仮定 𝑦𝑖 = 𝑦𝑖 = 𝛽0 + 𝛽1 𝑥𝑖1 + 𝛽2 𝑥𝑖2 + ⋯ + 「大体」の予測値 1 予測値は確率的に変動せず 説明変数の値のみで決まる 𝑦ො𝑖 = 𝛽0 + 𝛽1 𝑥𝑖1 + 𝛽2 𝑥𝑖2 + ⋯ ▶ 𝐱 𝑖 が決まれば𝑦ො𝑖 も決まる + 𝑒𝑖 予測とのズレ 2 誤差が平均0の正規分布に従う 𝑒𝑖 ∼ 𝑁 0, 𝜎𝑒2 【比較】資料05でのベイズ推定のときは 全員が同じ正規分布に従うと仮定していた 𝑦𝑖 ∼ 𝑁 𝜇, 𝜎 2 個人(𝑥𝑖 )ごとに 異なる分布から値が発生する 𝑦𝑖 = 𝑦ො𝑖 + 𝑒𝑖 𝑦𝑖 ∼ 𝑁(𝛽0 + 𝛽1 𝑥𝑖1 + 𝛽2 𝑥𝑖2 + ⋯ , 𝜎𝑒2 ) 09 モデル評価・仮説検証 10

stanで回帰分析モデル ▌まずは登場人物を確認しておく data { 【plate notation】 int N; array[N] real SALES; variables array[N] real DIST; array[N] real FLOOR; array[N] real ITEMS; } 𝛽𝑘 𝛽0 𝜎𝑒 𝑥𝑖𝑘 SALES𝑖 parameters { real beta_0; real beta_DIST; 𝑖 data real beta_FLOOR; real beta_ITEMS; real <lower=0> sigma; } 09 モデル評価・仮説検証 11

stanで回帰分析モデル ▌尤度と事前分布を書き足す data { model { beta_0 ~ normal(0, 100); array[N] real SALES; beta_DIST ~ normal(0, 100); array[N] real DIST; beta_FLOOR ~ normal(0, 100); array[N] real FLOOR; beta_ITEMS ~ normal(0, 100); array[N] real ITEMS; sigma ~ cauchy(0, 10); } 【plate notation】 variables int N; 𝛽𝑘 𝛽0 𝜎𝑒 𝑥𝑖𝑘 SALES𝑖 尤度はどうなる？ parameters { real beta_0; } real beta_DIST; 𝑖 data real beta_FLOOR; real beta_ITEMS; real <lower=0> sigma; } 09 モデル評価・仮説検証 12

回帰分析モデルの尤度の書き方 𝑦𝑖 ∼ 𝑁(𝛽0 + 𝛽1 𝑥𝑖1 + 𝛽2 𝑥𝑖2 + ⋯ , 𝜎𝑒2 ) ▌愚直に書けば SALES[1] ~ normal(beta_0 + beta_DIST*DIST[1] + beta_FLOOR*FLOOR[1] + beta_ITEMS*ITEMS[1], sigma); SALES[2] ~ normal(beta_0 + beta_DIST*DIST[2] + beta_FLOOR*FLOOR[2] + beta_ITEMS*ITEMS[2], sigma); SALES[3] ~ normal(beta_0 + beta_DIST*DIST[3] + beta_FLOOR*FLOOR[3] + beta_ITEMS*ITEMS[3], sigma); ----------------------------（以下省略）------------------------------ ▌forループを使って書けば for (i in 1:N){ SALES[i] ~ normal(beta_0 + beta_DIST*DIST[i] + beta_FLOOR*FLOOR[i] + beta_ITEMS*ITEMS[i], sigma); } ただstanではこのあたりもベクトル化して書けるのです 09 モデル評価・仮説検証 13

尤度をベクトル化してあげる 資料05 p.44 for (i in 1:N){ SALES[i] ~ normal(beta_0 + beta_DIST*DIST[i] + beta_FLOOR*FLOOR[i] + beta_ITEMS*ITEMS[i], sigma); } 𝑦𝑖 ∼ 𝑁(𝛽0 + 𝛽1 𝑥𝑖1 + 𝛽2 𝑥𝑖2 + ⋯ , 𝜎𝑒2 ) ▌今回の場合，最終的にこれでOK SALES ~ normal(beta_0 + beta_DIST*DIST + beta_FLOOR*FLOOR + beta_ITEMS*ITEMS, sigma); SALESとyhatが同じ長さなら ▲のコードは log 𝑃(SALES1 |𝑦ො1 , 𝜎) + log 𝑃(SALES2 |𝑦ො2 , 𝜎) + ⋯ + log 𝑃(SALES𝑛 |𝑦ො𝑛 , 𝜎) を計算してくれるのです DIST1 FLOOR1 ITEMS1 𝐲ො = 𝛽0 + 𝛽DIST DIST2 + 𝛽FLOOR FLOOR 2 + 𝛽ITEMS ITEMS2 DIST3 FLOOR 3 ITEMS3 線形代数的に計算できるならば normal()関数の平均パラメータは SALESと同じ長さのベクトルで表すことができる！ ベクトルならば，定数倍したり定数を足したりできますね 09 モデル評価・仮説検証 14

ただしちょっと厄介なことが DIST1 FLOOR1 ITEMS1 𝐲ො = 𝛽0 + 𝛽DIST DIST2 + 𝛽FLOOR FLOOR 2 + 𝛽ITEMS ITEMS2 DIST3 FLOOR 3 ITEMS3 ▌stanのarray型は線形代数の演算ができない arrayは複数の要素をまとめて格納できる箱だが，ただそれだけ ▶ arrayの定数倍や，array同士の四則演算はできない ▌各説明変数をarrayではなくベクトルとして与えてあげよう data { int N; } data { 修正 int N; array[N] real SALES; array[N] real SALES; array[N] real DIST; vector[N] DIST; array[N] real FLOOR; vector[N] FLOOR; array[N] real ITEMS; vector[N] ITEMS; } 09 モデル評価・仮説検証 vector型は 線形代数の演算ができるが 実数しか取れない ＝int型の値をいれることはできない 15

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; あとは予測のパートをいれるだけ real beta_FLOOR; real beta_ITEMS; real <lower=0> sigma; } 09 モデル評価・仮説検証 16

新しいブロックと新しい関数 ▌generated quantitiesブロック サンプリングされたパラメータ等を用いて四則演算するためのブロック。 ここでの計算内容はmodelブロック（＝事後分布）には影響しない。 ▌***_rng()関数 特定の確率分布から乱数を生成する関数 generated quantities { array[N] real SALES_pred; SALES_pred = normal_rng(beta_0 + beta_DIST*DIST + beta_FLOOR*FLOOR + beta_ITEMS*ITEMS, sigma); } 【ポイント】 • rng()関数は乱数を作る関数なので，左辺とは（分布を表す） ~ではなく（代入を表す） =でつなぐ • 一部のrng()関数は▲のようにベクトル化して書ける（できない場合はforループを使う） 09 モデル評価・仮説検証 17

generated quantitiesの仕組み 【サンプリングされたパラメータ】 iteration beta_0 beta_DIST beta_FLOOR beta_ITEMS sigma 【乱数生成分布】 iteration 店舗1 … 1 0.817 -0.102 0.058 -0.0013 0.28 1 𝑁(6.60, 0.28) … 2 0.124 -0.106 0.061 -0.0012 0.31 2 𝑁(6.44, 0.31) … 3 -0.173 -0.098 0.062 -0.0011 0.25 3 𝑁(6.43, 0.25) … 4 -0.208 -0.099 0.061 -0.0010 0.27 4 𝑁(6.43, 0.27) … ︙ ︙ ︙ ︙ ︙ ︙ ︙ ︙ … 【生成された乱数】 iteration 店舗1 … 1 6.74 … 2 5.81 … 3 6.33 … 4 6.49 … SALES_pred[1]の事後予測分布 ︙ ︙ … 乱数一回一回について 計算した結果を集めている 09 モデル評価・仮説検証 18

stanで回帰分析モデル
▌ようやくやりたいことがすべて揃いました
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は必ず最後のブロック
（modelブロックで得られた値を使うから？）

09 モデル評価・仮説検証

19

ということで推定しましょう library(cmdstanr) model_reg <- cmdstan_model("model_regression.stan") Rから与えるデータについては サイズさえあっていれば arrayかvectorかの区別は不要 stan_data <- list(N=100, SALES=dat$sales, DIST=dat$dist, FLOOR=dat$floor, ITEMS=dat$items) result <- model_reg$sample(data = stan_data) # 少し時間がかかるので，parallel_chainsを変えても良いかも result$summary() generated quantitiesも 全部出てくるから 見るのが大変だわね… 09 モデル評価・仮説検証 20

とりあえず収束くらいはチェックしておこう どこまで細かくチェックするかはお任せしますが，最低限の確認として… ▌ 𝑅෠ (Gelman-Rubin統計量) ▌相対ESS library(bayesplot) mcmc_neff(neff_ratio(result)) mcmc_rhat(result$summary()$rhat) mcmc_rhat(result$summary(NULL,rhat_basic)$rhat_basic)でもOK みんなほぼ1 あまり効率は よくなさそう… 09 モデル評価・仮説検証 21

事後予測チェック ▌つくった（統計）モデルは，現実世界にどれだけ近づけているか？ 現実世界 データ 近隣住民によって変わる オーナーの方針 気候によって変わる によって変わる dat$salesとSALES_predの分布が 似ていればとりあえず モデルはデータにフィットしている，と言える 周辺の競合 によって変わる 景気によって変わる 関数と事後分布 の組み合わせ 𝑦𝑖 = 𝛽0 + 𝛽1 𝑥𝑖1 + 𝛽2 𝑥𝑖2 + 𝜀𝑖 統計モデルから作ったデータ 09 モデル評価・仮説検証 （データで更新した）統計モデル 22

いよいよ事後予測チェック ▌bayesplotパッケージに大体あります （準備） ここで紹介する以外にもいろいろな関数があるので 是非試してみてください # この後の関数が使いやすい形(matrix)でSALES_predだけ抽出しておく この後の関数で 大変なことになるので 間引いておくと良い SALES_pred <- result$draws("SALES_pred", format="matrix")[seq(1,4000,by=500),] ▌分布の比較 ppc_dens_overlay(dat$sales, SALES_pred) 第１引数が実際のデータ，第２引数が事後予測で生成されたデータ 【SALES_predの形】 各行から薄青線が１つずつ引かれる 黒い線は 実データの密度 だいたい同じだったら 良い感じ 09 モデル評価・仮説検証 iteration 店舗1 店舗2 … 1 6.74 8.78 … 501 5.81 8.59 … 1001 6.33 8.62 … 1501 6.49 8.27 … ︙ ︙ ︙ … 23

説明変数との関係を見たり （準備） # この後の関数が使いやすい形(matrix)でSALES_predだけ抽出しておく SALES_pred <- result$draws("SALES_pred", format="matrix") ▌分布の比較 ppc_intervals(dat$sales, SALES_pred, x=dat$dist) 第3引数xには特定の説明変数を入れる 【店舗1の実数値とSALES_pred[1]の分布】 説明変数の分散説明率が高ければ， 𝑥に応じて𝑦もきれいに ො 右上がり（右下がり）に並んで表示されるはず 青い線は 予測分布の確信区間 （見にくいので先頭10店舗分だけ） 青い線に黒い点が 載っていたら良い感じ 黒い点は 実データの値 09 モデル評価・仮説検証 24

事後予測チェック ▌とりあえず作ったモデルが「データにはフィットしている」ことが確認できる 手元のデータにすらフィットしないモデルはダメに決まっているので… ▌ただし「手元にないデータにフィットするか（＝汎化性能）」については未知 データから予測したパラメータ（の事後分布）を使っている ▶ そのパラメータは手元のデータに過剰適合しているかもしれない ▌手元にない（未来の）データへのフィットを確認するためには パラメータ推定 80店舗分の データ データ 20店舗分の データ 予測データ作成 09 モデル評価・仮説検証 25

きちんと検証するなら ▌クロスバリデーション (cross-validation) 全てのデータ ① ② ③ パラメータ推定（学習）用データ 学習用データ テスト用データ 1. 2. 3. 4. テスト用データ テスト用データ 学習用データ 学習用データ 学習用データでパラメータを推定する 推定したパラメータによるモデルで学習用データを予測する（尤度などでも良い） 事後予測の精度が良ければGoodなモデルだということに！ このプロセスを①から③の分割パターンで繰り返す ※特に「テスト用データが１つ」をデータの数だけ繰り返す場合をLOOCVと呼ぶ（汎化誤差の近似） 09 モデル評価・仮説検証 26

（おまけ）回帰分析モデルをもっと簡潔に・効率よく書くには ▌勘の良い方はすでに気づいているかも DIST1 FLOOR1 ITEMS1 𝐲ො = 𝛽0 + 𝛽DIST DIST2 + 𝛽FLOOR FLOOR 2 + 𝛽ITEMS ITEMS2 DIST3 FLOOR 3 ITEMS3 𝛽0 1 DIST1 𝛽DIST 𝜷= , 𝐗 = 1 DIST2 𝛽FLOOR 1 DIST3 𝛽ITEMS FLOOR1 FLOOR 2 FLOOR 3 ITEMS1 ITEMS2 ITEMS3 とおくと 𝐲ො = 𝐗𝜷 betaをベクトルの形で，また説明変数（＋切片）を行列の形で与えたら良い！ 09 モデル評価・仮説検証 27

（おまけ）ベクトルと行列を駆使した書き方
【stanコード】

generated quantitiesは省略します

【Rコード】

model_regression_vec.stan
data {
int N;
int K; //説明変数の数

model <- cmdstan_model("model_regression_vec.stan")

array[N] real SALES;

# 先に切片項(1)をくっつけた説明変数の行列を用意しておく

matrix[N, K+1] X; //説明変数+切片項
}

dat_X <- cbind(1,dat[,c("dist","floor","items")])
stan_data <- list(N=100, K=3, SALES=dat$sales, X=dat_X)

parameters {
vector[K+1] beta; //回帰係数+切片

result <- model_reg_vec$sample(data = stan_data)

real <lower=0> sigma;
}
model {
beta ~ normal(0, 100); // betaの全要素に同じ事前分布
sigma ~ cauchy(0, 10);
SALES ~ normal(X*beta, sigma);

}

09 モデル評価・仮説検証

28

（さらにおまけ）stanに時々ある「効率化を極めた」関数
【stanコード】

「たまにこういう関数も用意されているよ」の例として

【Rコード】

model_regression_vec.stan
data {
int N;

model <- cmdstan_model("model_regression_vec.stan")

int K; //説明変数の数

stan_data <- list(N=100, K=3,

vector[N] SALES;
matrix[N, K] X; //説明変数

SALES=dat$sales, X=dat[,c("dist","floor","items")])

}

result <- model_reg_vec$sample(data = stan_data)
parameters {
real beta_0; // 切片
vector[K] beta; // 回帰係数
real <lower=0> sigma;

}
model {
beta_0 ~ normal(0, 100); // 切片の事前分布
beta ~ normal(0, 100); // betaの全要素に同じ事前分布

sigma ~ cauchy(0, 10);
SALES ~ normal_id_glm(X, beta_0, beta, sigma);
}

・切片を説明変数と別で与えるので「1をくっつける」作業が不要
・単純に計算が効率化されている
09 モデル評価・仮説検証

29

２ ベイズ統計的仮説検証 あえて「仮説検定」とは言わないことにします 09 モデル評価・仮説検証 30

事例 例 あなたはある講義を２つのクラスで担当しています。いま開発中のスペシャル指導法の効果を検 証するため，クラスAではスペシャル指導法を，クラスBではいままでどおりの指導法をしました。 実際に指導法の効果を確かめるため， ２つのクラスで同じテストをしました。さて，スペシャル指 導法には効果があるでしょうか？つまり，２つのクラスのテストの得点には差があるでしょうか？ 【データの読み込み】 この名前はなんでもOK "data_test.csv"をワーキングディレクトリに配置して ▌ 中身はこんなデータ dat <- read.csv("data_test.csv") ID:出席番号 出席番号で横並びになっているだけで 各行が1人ということではありません score_A: クラスAの各生徒のテスト得点 score_B: クラスBの各生徒のテスト得点 ▶ 標本における平均値差は3.175点 09 モデル評価・仮説検証 31

（標本理論）統計的仮説検定 1 仮説を設定する ▌まずは母集団分布を仮定する 母集団＝「その指導法を受けた（潜在的な）すべての生徒」 クラスA，Bは潜在的母集団の中からのランダムサンプリングに相当する，と考える テストの得点は正規分布に従う（と考えて良い）はず ▌2つの仮説を立てる 今回は両側検定にしておきます 2つの指導法によるテストの平均点の差を 𝛿𝜇 = 𝜇𝐴 − 𝜇𝐵 とおくと… 【実際に検証したいこと】 スペシャル指導法の方が 平均点が高い 𝛿𝜇 > 0 現実には どちらか 一方だけが 正しい alternative hypothesis 【◀の逆】 スペシャル指導法のほうが 平均点が高いことはない 𝛿𝜇 ≤ 0 null hypothesis 帰無仮説 対立仮説 09 モデル評価・仮説検証 32

（標本理論）平均値差の検定 2 帰無仮説が正しいときの検定統計量の分布を考える 𝑡= ▌検定統計量は推定のときと同じ考え方で 𝛿𝜇 𝑠Ƹ𝑘2 は群𝑘の不偏分散 1 1 𝑛1 𝑠Ƹ12 + 𝑛2 𝑠Ƹ22 + 𝑛1 𝑛2 𝑛1 + 𝑛2 − 2 母平均の検定の場合は，標本平均を用いて考えたら良い 𝛿𝜇 ≤ 0の中で 対立仮説に最も近いものとして 𝛿𝜇 = 0が正しい状況を 考えていきます 指導法によっては 平均点に差はない 𝛿𝜇 = 0 𝑡(78) 帰無仮説の元での分布 ▶ 帰無分布 null hypothesis 帰無仮説 こちらが正しい場合 𝛿𝜇 から算出した検定統計量 𝑡 は 自由度𝑛1 + 𝑛2 − 2の𝑡分布に従う 09 モデル評価・仮説検証 33

（標本理論）平均値差の検定 3 標本から得られた検定統計量の値の起こりやすさを考える ▌両側検定 今回はおよそ0.873 自由度78の𝑡分布において𝑡が標本から計算した値よりも極端になる確率は 極端に高い 𝑡(78) 𝑃(𝑡 ≥ 0.873) = 0.193 4 帰無仮説を棄却するかを判断する 𝑡が0.873になるような平均値差は 𝛿𝜇 ≤ 0 だとしても割と起こり得ること ▶ 帰無仮説が間違っているとは言い切れない 09 モデル評価・仮説検証 34

標本理論の仮説検定の問題点 【実際に検証したいこと】 スペシャル指導法の方が 平均点が高い >0 現実には どちらか 一方だけが 正しい 【◀の逆】 スペシャル指導法のほうが 平均点が高いことはない 0 ▌帰無仮説は噛ませ犬 極論を言えば「平均値差がちょうどゼロ」なんてありえない （特に両側検定のとき） ▌「帰無仮説が保持された」≠「帰無仮説が正しい」 一言で言えば「背理法だから」 ▌𝑝値は仮説の強さを表現したものではない 𝑝 = 0.04の時と𝑝 = 0.0001の時を比べて 「𝑝 = 0.0001 の時の方が平均値差が大きい」などとは言えない ベイズ統計的にはどうなるのでしょうか？ 09 モデル評価・仮説検証 35

直感的には ▌事後分布が「信念」を表すとしたら スペシャル指導法のほうが 平均点が高いことはない 𝑃(𝛿𝜇 ≤ 0|𝑌) スペシャル指導法の方が 平均点が高い 𝑃(𝛿𝜇 > 0|𝑌) ◀の場合 𝛿𝜇 > 0 という設定のほうが 信念が強いと考えられる 09 モデル評価・仮説検証 36

ということで，まずはstanでモデルを書こう ▌基本方針 1. まずは各群（指導法）の母平均値を推定する 2. 推定した母平均値の差を計算し，事後分布を求める ▌と言っても難しい事はありません まずは普通に「各群の母平均値を推定するstanコード」を書いてみましょう 【ヒント】 (for group A) • 尤度は正規分布に従っています 𝜇𝐴 • 標準偏差はわからないので，同じように推定しましょう • とりあえず「クラスA」についてのstanコードを書いてみましょう。 𝑥𝐴𝑖 𝜎𝐴 𝑖 data • 「クラスB」も「クラスA」と構造は全く同じです。ということは？ 次ページに回答例が載っています 09 モデル評価・仮説検証 37

クラスAのみのモデル
data {

int N_A;
array[N_A] real SCORE_A;
}

資料05 p. 105
2つの群の人数が違っても良いように，添字Aをつけておきました
※今回使用するデータに関しては，人数は同じなので無くてもOK
テスト得点は整数しか入っていないが，ここではreal型を指定する

parameters {

こちらも添字Aをつけておきました

real <lower=0, upper=100> mu_A;

real <lower=0> sigma_A;
}

mu_Aは平均点なので0-100の値しか取らない
▶ lower, upperをつけておきました（今回のデータではどちらでも良い）

model {
mu_A ~ normal(50, 100);
sigma_A ~ cauchy(0, 10);

クラスBも同じ構造＆独立なので
単純に同じように書き足していけば…

SCORE_A ~ normal(mu_A, sigma_A);

}

09 モデル評価・仮説検証

38

2クラス同時に推定するモデル data { 【plate notation】 model { int N_A; mu_A ~ normal(50, 100); array[N_A] real SCORE_A; sigma_A ~ cauchy(0, 10); int N_B; SCORE_A ~ normal(mu_A, sigma_A); array[N_B] real SCORE_B; mu_B ~ normal(50, 100); } sigma_B ~ cauchy(0, 10); SCORE_B ~ normal(mu_B, sigma_B); parameters { real <lower=0> sigma_A; real <lower=0, upper=100> mu_B; } 𝜎𝐴 𝜇𝐵 } real <lower=0, upper=100> mu_A; real <lower=0> sigma_B; 𝜇𝐴 𝑥𝐴𝑖 𝑥𝐵𝑖 𝑖 data 𝑖 data 赤：クラスAに関する部分 青：クラスBに関する部分 あとは母平均値の差𝛿𝜇 を用意するだけです 09 モデル評価・仮説検証 39 𝜎𝐵

更に追加 【ベイズによる統計的推測のポイント】 知りたいパラメータを一つにまとめ上げること 一個の事後分布が得られる形に してあげましょう，ということ 今回の例では「平均点の差」なのでmu_A – mu_Bがあるとうれしい 𝛿𝜇 = 𝜇𝐴 − 𝜇𝐵 ▌generated quantitiesブロックの出番です generated quantities { } real delta_mu; 直接delta_muをサンプリングするのではなく delta_mu = mu_A - mu_B; サンプリングされたmu_Aとmu_Bの差をとっている 09 モデル評価・仮説検証 40

完成したモデル（パターン1）
model_ttest.stan
data {

【plate notation】

model {

int N_A;

mu_A ~ normal(50, 100);

array[N_A] real SCORE_A;

sigma_A ~ cauchy(0, 10);

int N_B;

SCORE_A ~ normal(mu_A, sigma_A);

array[N_B] real SCORE_B;

mu_B ~ normal(50, 100);

}

sigma_B ~ cauchy(0, 10);

SCORE_B ~ normal(mu_B, sigma_B);
parameters {

}

𝜇𝐴

𝜎𝐴

𝜇𝐵

𝜎𝐵

}

real <lower=0, upper=100> mu_A;
real <lower=0> sigma_A;

𝛿𝜇

generated quantities {

real <lower=0, upper=100> mu_B;

real delta_mu;

real <lower=0> sigma_B;

delta_mu = mu_A - mu_B;
}

𝑥𝐴𝑖

𝑥𝐵𝑖

𝑖 data

𝑖 data

他のパラメータによって決定する値は
白二重丸で表していきます

では早速推定してみましょう！
09 モデル評価・仮説検証

41

stanで推定 ▌推定 library(cmdstanr) model <- cmdstan_model("model_ttest.stan") stan_data <- list(N_A=40, N_B=40, この程度のモデルだとすぐ終わりますが • 複数chainを同時に走らせたり • サンプリング数を増やしたり してみても良いですね。 SCORE_A=dat$score_A, SCORE_B = dat$score_B) result <- model$sample(data = stan_data) ▌事後分布を出してみたりして library(bayesplot) mcmc_dens(result$draws(), pars = "delta_mu") 09 モデル評価・仮説検証 42

generated quantitiesの仕組み iteration 1 2 3 4 … mu_A 69.98 59.03 63.64 56.71 … mu_B 63.25 59.19 63.18 64.21 … iteration 1 2 3 4 … mu_A – mu_B 6.73 -0.16 0.46 -7.50 … 乱数一回一回について 計算した結果を集めている delta_muの事後分布 09 モデル評価・仮説検証 43

あとはお好きなように

事後分布があればなんとでもなる

スペシャル指導法の方が
平均点が高い
𝑃 𝛿𝜇 > 0 𝑌

▌シンプルに𝑃 𝛿𝜇 > 0 𝑌 を求めてみる
draws <- result$draws(format = "df")

> draws
# A draws_df: 1000 iterations, 4 chains, and 6 variables

mean(draws$delta_mu > 0)

lp__ mu_A sigma_A mu_B sigma_B delta_mu

[1] 0.8085

1

-255

60

18

65

15

-5.39

2

-252

62

18

61

13

1.32

3

-256

63

18

67

12

-3.99

4 -254

66

20

63

17

3.10

------------------------（以下略）----------------------------------

80.85%

09 モデル評価・仮説検証

この列のうち
0以上の値の割合を
出せば良い！

44

点帰無仮説の場合は？ ▌ある仮説に対する信念は事後分布の面積と考えられる ▶ ある（帰無）仮説が点仮説の場合どうする？ 【実際に検証したいこと】 スペシャル指導法には 何らかの効果がある 𝛿𝜇 ≠ 0 現実には どちらか 一方だけが 正しい 【◀の逆】 スペシャル指導法は 従来の指導法と同じ 𝛿𝜇 = 0 alternative hypothesis null hypothesis 帰無仮説 対立仮説 ▌連続確率変数の確率分布 ある点を取る確率（e.g., 𝑃 𝛿𝜇 = 0|𝒀 ）はゼロとして扱われる ベイズの場合，帰無仮説は自動的に否定されてしまうのか？ 09 モデル評価・仮説検証 45

（考え方1）統計的仮説検定のように考える ▌標本分布と仮説検定の関係性 詳しくはこちら 検定統計量が棄却域に含まれる 母数の95%信頼区間が ▶ 帰無仮説が棄却される 帰無仮説の値𝛿𝜇 = 0を含まない 【今回の例の場合】 𝑡(78) 𝛿𝜇 の95%信頼区間 𝑝 = 0.385 0を含んでいる 09 モデル評価・仮説検証 46

• https://www2.kobe-u.ac.jp/~bunji/files/lecture/stat/stat-11-example-hypothesis-testing.pdf#page=42.00

ということは…？ ▌両側検定で帰無仮説が棄却されるとき すなわち「有意な効果がある」と判断されるのは ギリギリ棄却されるときの 𝛿𝜇 の95%信頼区間 2.5% （両側検定の場合）標本分布を 帰無仮説の値で区切ったときに 小さい方の面積が2.5%以下なら 帰無仮説は棄却されている！ 95% 標本における平均値差が7.246点だった場合 09 モデル評価・仮説検証 47

方向性確率 (probability of direction [PD]: Makowski et al., 2019) ▌前ページと同じ考え方を事後分布に当てはめたもの あとはお好きなように 1. 事後分布を，「効果がない」を表す点 (point null) で区切る 事後分布があればなんとでもなる スペシャル指導法の方が 2. 大きい方の面積がPD 正確には「中央値がある方の面積」 【p. 44の場合】 シンプルに ( > 0| ) を求めてみる d a s a s d a s d t d a s at d 平均点が高い ( > 0| ) 【ポイント】 • 定義上，必ず0.5から1の値を取る d a s • 仮説検定の𝑝値と直接対応している d a s_d t at s c a s a d va a s ta_ 片側検定の場合は 𝑝 𝑝𝑑 = 1 − 𝑝 𝑝𝑑 = 1. − 2 . . ▲のため，非ベイジアンにも理解されやすい __ . • これがPD モデル _ s a_ _ s a_ d ta_ . （以下略） （Makowskiらの提案する基準） pd <= 95% ~ p > .1: uncertain pd > 95% ~ p < .1: possibly existing この列のうち pd > 97%: likely existing ◀これくらいは必要？ 以上の値の割合を pd > 99%: probably existing 出せば良い！ pd > 99.9%: certainly existing ・仮説検定 09 モデル評価・仮説検証 48

• https://doi.org/10.21105/joss.01541

（考え方2）実質的な効果の大きさで考える ▌帰無仮説が「点」だからややこしくなっている 厳密な 𝑃 𝛿𝜇 = 0 は定義上0になってしまうため ▌代わりに「実質的に効果のない区間」を決めてあげたら良いのでは？ 統計的仮説検定でも「実質的に効果がないのにサンプルサイズが大きいから有意」 というのが問題視されますからね 0.5点の差もないなら もう実質的に同じじゃないですかね 事後分布と「±0.5の範囲」を比較して 実質的な効果の有無を判断する 09 モデル評価・仮説検証 49

実質的等価領域 (region of practical equivalence [ROPE]) ▼これのことです。 【ROPEの幅の決め方 by Kruschke】 • 結果に大きく影響するので慎重に決める必要あり • 効果量の基準 (Cohen, 1988)に則るならば 標準化したスケールにおいて±0.1とか？ ▌事後分布がどの程度ROPEの中に含まれているか，に基づいて評価を行う ▶ 基本的には • ROPEにすっぽり入っていたら「実質的に効果なし」 • ROPEに全く入っていなかったら「実質的に効果あり」 • 中途半端だったら「確信は持てない」 09 モデル評価・仮説検証 50

• https://easystats.github.io/bayestestR/articles/region_of_practical_equivalence.html#how-to-define-the-rope-range

ROPEによる評価方法その1 （Makowskiらはこちらをオススメ） ▌事後分布全体を使った方法 スペシャル指導法には 何らかの効果がある 𝛿𝜇 ≠ 0 v.s. スペシャル指導法は 従来の指導法と同じ 𝛿𝜇 = 0 事後分布の大半が ROPEの中にある場合 事後分布がROPEの 中にも外にもある場合 事後分布の大半が ROPEの中にない場合 ▶ 実質的に効果なし ▶ 効果に確信なし ▶ 実質的に効果あり Makowskiら的には97.5%以上 Makowskiら的には2.5%以下 09 モデル評価・仮説検証 51

ROPEによる評価方法その2 この方法の場合，確信区間の幅(95%なのか90%なのか)も重要となります ▌事後分布の確信区間 (HDI) を使った方法 確信区間は「信念の大半」を表す区間と言えます スペシャル指導法には 何らかの効果がある 𝛿𝜇 ≠ 0 v.s. スペシャル指導法は 従来の指導法と同じ 𝛿𝜇 = 0 確信区間が完全に ROPEの中にある場合 確信区間がROPEの 中にも外にもある場合 確信区間が完全に ROPEの中にない場合 ▶ 実質的に効果なし ▶ 効果に確信なし ▶ 実質的に効果あり 09 モデル評価・仮説検証 52

（考え方3）point nullの信念の相対的な強さを考える ▌最も強い信念と比べてpoint nullの信念はどれくらい強いのか？ ▶ MAP-based 𝑝-value (Mills, 2017) 𝑝MAP = 尤度比検定と同じような発想 𝑃(𝛿𝜇 = 0|𝑌) MAP 𝑃 MAP 𝛿𝜇 𝑌 【右図の場合】 0.073 𝑝MAP = ≃ 0.676 0.108 ▶ 𝑝MAP > 0.05なので 帰無仮説は棄却されない ここまで来るとかなり 仮説検定感が強いですね 09 モデル評価・仮説検証 53

Rで出してみる ▌このあたりの話をまとめたbayestestRパッケージを使います いろいろ出せる関数 library(bayestestR) describe_posterior(result$draws("delta_mu"), test=c("pd","rope","p_map")) Summary of Posterior Distribution Parameter | Median | 95% CI | p (MAP) | pd | ROPE | % in ROPE ---------------------------------------------------------------------------------delta_mu | 3.10 | [-4.02, 10.11] | 0.714 | 80.85% | [-0.10, 0.10] | 1.87% ▌調整したい場合は rope_ci=1にすると describe_posterior(result$draws("delta_mu"), test=c("pd","rope","p_map"), 事後分布全体 ci=0.9, ci_method="HDI", rope_range = c(-0.3,0.3), rope_ci = 1) 確信区間の幅と算出方法 ROPEの幅 09 モデル評価・仮説検証 ROPEの判定に用いる区間の幅 54

モデルの書き方を変えてみる ▌p. 41の書き方ではdelta_muに直接事前分布を置くことができない つまりdelta_mu（スペシャル指導法の効果）には事前の信念を置くことができない 【例】スペシャル指導法には3点くらい平均点を上げる効果があると思う。しらんけど delta_muをパラメータとしてあげるのもあり 【p. 14の表現】 【別の表現】 𝛿𝜇 𝛿𝜇 𝜇𝐴 𝜎𝐴 𝜇𝐵 𝜎𝐵 𝜇𝐴 𝜎𝐴 𝜇𝐴𝐿𝐿 𝜇𝐵 𝑥𝐴𝑖 𝑥𝐵𝑖 𝑥𝐴𝑖 𝑥𝐵𝑖 𝑖 data 𝑖 data 𝑖 data 𝑖 data 09 モデル評価・仮説検証 全体平均 𝜇𝐴𝐿𝐿 と 差 𝛿𝜇 に分解する 𝜎𝐵 55

モデルを少し書き換えてみよう ▌p. 14での書き方 generated quantities { real delta_mu; 直接delta_muをサンプリングするのではなく delta_mu = mu_A - mu_B; サンプリングされたmu_Aとmu_Bの差をとっている } ▌今回の表現では？ 𝛿𝜇 𝜇𝐴 𝜎𝐴 𝜇𝐴𝐿𝐿 𝜇𝐵 𝜇𝐴 , 𝜇𝐵 はいずれも 𝜎𝐵 • 他のパラメータから作られるものだが • それ自体が尤度に影響している ▶ generated quantitiesブロックではダメ！ 𝑥𝐴𝑖 𝑥𝐵𝑖 𝑖 data 𝑖 data 09 モデル評価・仮説検証 56

新しいブロックを使って書き直す
▌まずはparametersブロック
parameters {

これらを使って各群の平均パラメータ

real delta_mu;
real <lower=0, upper=100> mu_all;
︙

▌新しいブロック

【transformed parametersブロック】

transformed parameters {

real <lower=0, upper=100> mu_A;
real <lower=0, upper=100> mu_B;
mu_A = mu_all + delta_mu/2;
mu_B = mu_all - delta_mu/2;
}

mu_Aとmu_Bを作ります

サンプリングされたパラメータを
別のパラメータに変換するためのブロック。
ここでの計算内容はmodelブロック（＝事後分布）に
用いることができる。

得られるサンプルの形式などは
transformed parametersもgenerated quantitiesも
(parametersブロックのときと）全く同じです
09 モデル評価・仮説検証

57

完成したモデル（パターン2）
model_ttest_v2.stan
data {

transformed parameters {

int N_A;

real <lower=0, upper=100> mu_A;

array[N_A] real SCORE_A;

real <lower=0, upper=100> mu_B;

int N_B;

mu_A = mu_all + delta_mu/2;

array[N_B] real SCORE_B;

mu_B = mu_all - delta_mu/2;

}

}

parameters {

model {

parameters と modelブロックの間
（modelブロックで使うから）

real <lower=0, upper=100> mu_all;

mu_all ~ normal(50, 100);

real delta_mu;

delta_mu ~ normal(0, 10);

real <lower=0> sigma_A;
real <lower=0> sigma_B;

sigma_A ~ cauchy(0, 10);

}

【plate notation】

𝛿𝜇

𝜇𝐴

𝜎𝐴

𝜇𝐴𝐿𝐿
𝜇𝐵

𝑥𝐴𝑖

𝑥𝐵𝑖

𝑖 data

𝑖 data

𝜎𝐵

SCORE_A ~ normal(mu_A, sigma_A);
sigma_B ~ cauchy(0, 10);
SCORE_B ~ normal(mu_B, sigma_B);
事前分布の置き方の違いなので
データが多ければ
どちらの書き方でも結局変わりません

}

09 モデル評価・仮説検証

58

ちょっとまってよ ベイズでは事前分布による恣意性が問題となる 【試しに主張の強い事前分布をおいてみると】 逆に言えば，そういう事前情報を 考慮した推論が可能，というのが ベイズ統計のメリットなわけですが あのスペシャル指導法，どうせ効果なんかないと， delta_muの 事後分布 私はそう思ってますぜ，ダンナ。 事前分布 delta_mu ~ normal(0, 0.5) パターン1であれば mu_A ~ normal(mu_all, 0.5); mu_B ~ normal(mu_all, 0.5); （どちらの平均点もほぼ同じ） 𝑝𝑑 = 0.5487 𝑅𝑂𝑃𝐸 ±0.3 = 46.17% 𝑝MAP = 0.990 09 モデル評価・仮説検証 どうやら 効果なし？ 59

ちょっとまってよ ベイズでは事前分布による恣意性が問題となる 【対処法】 1. 事前分布を無情報・弱情報分布にする 2. データをもっと増やす 3. 事後分布の代わりに「データはどちらの仮説を支持するか」を評価する モデル比較の話（ベイズファクター）へと続きます 09 モデル評価・仮説検証 60

まとめと次回予告 【まとめ】 ▌ベイズ統計におけるモデル評価の話がある程度分かりました とりあえず「モデルがデータにフィットしているか」は事後予測チェックで検討可能 ▌ベイズ統計における仮説の評価方法がある程度分かりました 今回の内容は「パラメータの事後分布を使った方法」 ただしこれでは事前分布の影響を強く受けてしまう 【次回予告】 ▌ベイズ統計的「モデル比較」を考えていきます 仮説検証や予測の目的は「より正しい設定を明らかにすること」 ▶ モデルを直接比べることで更に迫ることができるかもしれないのです 09 モデル評価・仮説検証 61