LT: Shallow Dive into Bayes Factor

A presentation slide for a study meeting on statistics in Japan.
Updated: Jul.2022

Meeting web page: https://connpass.com/event/204931/

1. The Difference Between Traditional Frequentism and Bayesianism

2. What is the Bayes Factor?

3. A simple experiment in R using RStan

Transcript

1. Maxwell_110 Shallow Dive into Bayes Factor
2. 仕事: 外資系生命保険会社でデータ分析などに従事 （以前は損保で収益・リスク分析など．現在は大学に研究員としても出向中） 趣味: データ分析（Kaggle など） あてのない散歩，MTG 近況: 論文執筆や雑誌寄稿で職業作家状態 （データ分析がしたい！） @Maxwell_110
3. 最初に・・・ ✓ 仮説検定 ✓ ベイズの基本（ベイズの定理など） ✓ MCMC（Rstan） の基本的なところは前提知識とします 本 LT は 10分と限られた時間のため，ご了承ください m(_ _)m
4. 1. 従来の頻度主義とベイズ主義の違い 2. ベイズファクターとは？ 3. による簡単な実験
5. 伝統的な統計学 （頻度論） ベイズ統計学 仮説検定 パラメータ 𝜃 （母数） データ 𝑋 （標本） 1. 従来の頻度論とベイズ統計学の違い • 未知の定数 • 真値は１つ • データから 最尤法で推定 • 母集団から得られる 標本 • 確率変数 • 観測された定数 • 確率変数 • 事前分布 事後分布 （ベイズの定理） • 帰無仮説 H 0 と対立仮説 H 1 • 帰無仮説 H 0 の元で データが得られる確率 を検証する • 「帰無仮説が正しくない」 ことを示す手法 H 1 ? H 0 ? どーするの？
6. 頻度論における 95% 信頼区間 (95% CI) 何度も同じサンプルサイズの標本データを取ると， 真値が 95% の確率で信頼区間内に入る 母集団 （真値：𝜃0 ） 標本データは確率変数 𝜃 母集団の真値： 𝜃0 例えば，100 個の標本デ ータによる信頼区間があっ た時，5 回真値が入らない 真値が入っていない例
7. ベイズ統計学における 95% 信用区間 MCMC などでサンプリングして求めた事後分布は，母数の確率分布になっている． そのため，「95% の確率でその範囲内に真値がある」ということができる． 𝑓 𝜃|𝑋 = 𝑓 𝑋|𝜃 𝑓 𝜃 𝑓 𝑋 事後分布 事前分布 尤度 周辺尤度 （基準化定数，エビデンス） 事後分布の形状まで求めている （複雑な分布でも MCMC で求まる） 𝜃 ベイズ信用区間 (BCI) MAP 推定値
8. 2. ベイズファクターとは？ ✓ ベイズ統計学において仮説の評価を行うことができる （Jeffreys, 1935. Hoijtink, Klugkist, and Boelen, 2008. Hoijtink, 2011.） ✓ 母数に対する仮説を「直接」評価 頻度論の場合は，仮説の元で標本が得られる確率を評価 ✓ 頻度論における仮説検定とは似て非なるもの
9. 二群の平均に関して以下の仮説を例として考える 𝐻𝑖 : 𝜇1 ≥ 𝜇2 上式のように， 母数に対して不等式制約を用いて表される仮説 𝐻𝑖 を 「情報仮説」と呼ぶ ベイズファクターを簡単な例で考えてみよう Informative hypothesis
10. 情報仮説に対して，母数に制約を課さない仮説 𝐻𝑢 を 「無制約仮説」と呼び， この場合，ベイズファクターは 情報仮説と無情報仮説で 計算する 𝐻𝑢 : 𝜇1 𝜇2 𝐻𝑖 : 𝜇1 ≥ 𝜇2 一方で， 前頁の「情報仮説」 Unconstrained hypothesis
11. では，𝑯𝒊 と 𝑯𝒖 におけるベイズファクターとは？ 𝐵𝐹𝑖𝑢 = 𝑝𝑜𝑠𝑡𝑒𝑟𝑖𝑜𝑟 𝑜𝑑𝑑𝑠 𝑝𝑟𝑖𝑜𝑟 𝑜𝑑𝑑𝑠 = ൘ 𝑃 𝐻𝑖 |𝑋 𝑃 𝐻𝑢 |𝑋 ൘ 𝑃 𝐻𝑖 𝑃 𝐻𝑢 観察データ 𝑿 によって事前確率の比 が 事後確率の比へとどれだけ 変化したかを計算 事前確率の比（事前オッズ） （事前確率とは，任意の仮説において事前分布を仮説に与する母数 空間で積分したもの） （事後確率とは，任意の仮説において事後分布を仮説に与する母数空間 で積分したもの） 事後確率の比（事後オッズ）
12. 事前オッズ 𝑃 𝐻𝑢 = ׭ 𝑓 𝜇1 𝑓 𝜇2 𝑑𝜇1 𝑑𝜇2 = 1 𝑃 𝐻𝑖 = ׭ 𝜇1≥𝜇2 𝑓 𝜇1 𝑓 𝜇2 𝑑𝜇1 𝑑𝜇2 = 1 2 ൗ 𝑃 𝐻𝑖 𝑃 𝐻𝑢 = Τ 1 2 になる 互いに独立な正規分布: 𝑁 0, 𝐷−1 （等高線は確率密度の大きさ） 𝜇1 𝜇2 𝐻𝑢 𝜇1 ≥ 𝜇2 𝜇1 𝜇2 𝐻𝑖 事前オッズは，𝑓(𝜃) を使って計算する 互いに独立な正規分布 𝑓 𝜇1 および 𝑓 𝜇2 を考えると，
13. 事後オッズ 𝐻𝑢 (a) (b) 𝑃 𝐻𝑢 | 𝑋 = ׭ 𝑓 𝜇1 , 𝜇2 | 𝑋 𝑑𝜇1 𝑑𝜇2 = 1 𝑃 𝐻𝑖 | 𝑋 = ׭ 𝜇1≥𝜇2 𝑓 𝜇1 , 𝜇2 | 𝑋 𝑑𝜇1 𝑑𝜇2 データ 𝑋 によって事前分布が更新された結果， （a）正の相関 もしくは（b）負の相関 をもち， かつ，平均が（0, 0）ではないような事後分布が 得られた二つのケース（a）と（b）を考える 𝑃 𝐻𝑖 |𝑋 ≫ Τ 1 2 ൗ 𝑃 𝐻𝑖|𝑋 𝑃 𝐻𝑢|𝑋 ≫ Τ 1 2 𝐵𝐹𝑖𝑢 ≫ 1 (a) (b) 𝑃 𝐻𝑖 |𝑋 ≪ Τ 1 2 ൗ 𝑃 𝐻𝑖|𝑋 𝑃 𝐻𝑢|𝑋 ≪ Τ 1 2 𝐵𝐹𝑖𝑢 ≪ 1 事後オッズは，𝑓 𝜃|𝑋 を使って計算する 𝐻𝑖
14. Okada, 2013 ベイズファクターに対する基準 ✓ Jeffreys や Kass らの基準がよく知られている ✓ あくまでも大雑把な経験則で絶対的なものでは ない点に留意 ✓ p 値と似たような感覚で扱うと良いとのこと （Rosnow & Rosenthal, 1989）
15. 今回は，無制約仮説に対する情報仮説の元でベイズファクターを考えたが， 無制約仮説ではなく相補仮説（情報仮説の補空間）の元で計算されることも 増えてきている 詳しくは 岡田, ベイズ推定による情報仮説の評価：その理論と各種モデルへの応用について 専修人間科学論集心理学篇, 2016, 6, 9-17 https://core.ac.uk/download/pdf/71799133.pdf などを参照
16. 3. による簡単な実験 岡田, ベイズ統計による情報仮説の評価は分散分析にとって代わるのか? 基礎心理学研究, 2013, 32巻, 2号, p.223-231 https://www.jstage.jst.go.jp/article/psychono/32/2/32_KJ00009351488/_pdf/-char/ja ここでの実験は以下の内容の再現実験 但し，OpenBUGS ではなく RStan を使用 （一部，事前分布の設定も異なる）
17. Dataset 27 匹の鼠を 3 つの群 LD : 通常の明暗サイクル LL : 常に明るい光をつけた状態 DM : 日中は明るく，夜は薄暗い状態 に分けて，生活（食住つき）させた実験† Fonken et al., 2010 の結論は 「夜間の光は，体重 (Body Mass) を増加させる」 だが・・・，これをベイズファクターで検証してみる (少し可哀想な実験ですが，おつきあください m(_ _)m) † Fonken, L., et. al., "Light at night increases body mass by shifting time of food intake," Proceedings of the National Academy of Sciences, October 26, 2010; 107(43): 18664-18669. Lock5Data::LightatNight4Weeks
18. モデルと仮説 ✓ 3 つの群はそれぞれ平均 𝜇∎∎ の異なる正規分布に従うと仮定 （但し，分散は同一とする） ✓ 以下の 2 つの情報仮説 𝐻1 , 𝐻2 を無制約仮説 𝐻𝑢 に対して検証 𝐻1 : 𝜇𝐿𝐷 < 𝜇𝐷𝑀 < 𝜇𝐿𝐿 𝐻2 : 𝜇𝐿𝐷 < 𝜇𝐷𝑀, 𝜇𝐿𝐿 𝐻𝑢 : 𝜇𝐿𝐷, 𝜇𝐷𝑀, 𝜇𝐿𝐿
19. 𝐵𝐹𝑖𝑢 = 𝑝𝑜𝑠𝑡𝑒𝑟𝑖𝑜𝑟 𝑜𝑑𝑑𝑠 𝑝𝑟𝑖𝑜𝑟 𝑜𝑑𝑑𝑠 = ൗ 𝑃 𝐻𝑖|𝑋 𝑃 𝐻𝑢 | 𝑋 ൗ 𝑃 𝐻𝑖 𝑃 𝐻𝑢 = 𝑃 𝐻𝑖 | 𝑋 𝑃 𝐻𝑖 無制約仮説のもとでは だったので，𝑃 𝐻𝑖 と 𝑃 𝐻𝑖 | 𝑋 だけ計算すればよい 𝐻1 : 𝜇𝐿𝐷 < 𝜇𝐷𝑀 < 𝜇𝐿𝐿 𝐻2 : 𝜇𝐿𝐷 < 𝜇𝐷𝑀, 𝜇𝐿𝐿 𝑃 𝐻1 = Τ 1 6 𝑃 𝐻2 = Τ 1 3 全母数空間で事前分布を 積分するので，事前確率は 1 全母数空間で事後分布を 積分するので，事後確率は 1 𝑃 𝐻𝑖 は簡単に計算できて・・・
20. Stan model code data { int N; int d1[N]; int d2[N]; int d3[N]; real Y[N]; } parameters { real mu[3]; real sig2; } model { for (i in 1:3) { mu[i] ~ normal(0, 20); } sig2 ~ lognormal(0, 20); for (n in 1:N) { Y[n] ~ normal(mu[1] * d1[n] + mu[2] * d2[n] + mu[3] * d3[n], sig2); } } generated quantities { real f1; real f2; f1 = int_step(mu[2] - mu[1]) * int_step(mu[3] - mu[2]); f2 = int_step(mu[2] - mu[1]) * int_step(mu[3] - mu[1]); } 𝑃 𝐻𝑖 |𝑋 は MCMC sampling で求める サンプリングされた事後分布のうち，各情報仮説の条件に与 する割合を計算し推定する（Naive Monte Carlo） 本実験の場合，サンプリングされた 4000点のうちどれだけの割合が条件を満た すかを計算．但し，このやり方は事後分布の形状によっては粗い近似になってし まっている可能性がある．精緻な値を求めるためには，Bridge Samplingを使用 するのが良い（岡田, 2018など）． 事前分布は，平均は各群毎に弱情報事前分布 を分散は群に依らず同一の弱情報事前分布を 設定※ ※ 岡田, 2013 の設定から大きく変更している．岡田 2013 の設定だと収束しない（そもそも弱情報事前分布になっていなさそう）が，恐らくは歳月を経てデータセットの次元が変わったため？（未確認） f1 で情報仮説 H 1 に与するサンプリングかどうかを計算．H 1 の条件 にあてはまる時は 1 に，そうでない時は 0 になる． 同様に，f2 は H 2 の条件にあてはまる時は 1 に，そうでない時は 0 になる． 𝐻1 : 𝜇𝐿𝐷 < 𝜇𝐷𝑀 < 𝜇𝐿𝐿 𝐻2 : 𝜇𝐿𝐷 < 𝜇𝐷𝑀, 𝜇𝐿𝐿
21. Okada, 2013 ✓ 事前分布の設定などが異なるが岡田（2013）とほぼ 同じ結果 ✓ Kass & Raftery の基準に従うと H 1 の仮説が positive となり，H 2 の仮説はぎりぎり positive にはならない ✓ 仮説としては H 1 の方がエビデンスレベルが強い # 6. Compute Bayes Factor ---- mcmc.sample <- rstan::extract(stan.result) f1 <- mean(mcmc.sample[["f1"]]) f2 <- mean(mcmc.sample[["f2"]]) BF1u <- f1 / (1 / 6) BF2u <- f2 / (1 / 3) cat("BF1u:", BF1u, "¥n") # 5.601 cat("BF2u:", BF2u) # 2.823 bayes_factor.R (L84 ~ 93)
22. R Codes BF_model1.stan stan model の定義 bayes_factor.R MCMC の実行・BF の計算
23. References 1. 岡田, ベイズ統計による情報仮説の評価は分散分析にとって代わるのか? 基礎心理学研究, 2013, 32巻, 2号, pp.223-231 2. 岡田, ベイズ推定による情報仮説の評価：その理論と各種モデルへの応用について 専修人間科学論集心理学篇, 2016, 6, pp.9-17 3. 岡田, ベイズファクターによる心理学的仮説・モデルの評価 心理学評論, 2018, 61(1), pp.101-115. 4. 浜田，石田，清水 社会科学のためのベイズ統計モデリング 朝倉書店 2019 p.109 https://amzn.to/30Xiqar 5. Jeffreys H. Some Tests of Significance, Treated by the Theory of Probability. Math Proc Cambridge Philos Soc. 1935 Apr 24;31(2):203–22. https://www.cambridge.org/core/product/identifier/S030500410001330X/type/journal_article 6. Kass RE, Raftery AE. Bayes Factors. J Am Stat Assoc. 1995 Jun;90(430):773–95. http://www.tandfonline.com/doi/abs/10.1080/01621459.1995.10476572 7. Hoijtink H. Objective Bayes Factors for Inequality Constrained Hypotheses. Int Stat Rev. 2013 Aug;81(2):207–29. http://doi.wiley.com/10.1111/insr.12010