LT: Shallow Dive into Bayes Factor

A presentation file for study meeting of statistics in Japan.

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. Senior Data Scientist A. I. A. J Maxwell Shallow Dive into Bayes Factor
2. 自己紹介 Maxwell 仕事： 外資系生命保険会社でデータ分析 （以前は損保で収益・リスク分析など・・・ 現在は大学に研究員として出向中） 趣味： データ分析（Kaggle など） MTG 近況： 論文執筆で鬱気味 （論文が帰無仮説状態・求む！優しい査読！） Maxwell_110
3. 最初に・・・  仮説検定  ベイズの基本（ベイズの定理など）  MCMC（Rstan） の基本的なところは前提知識とします 本 LT は 10分と限られた時間のため，ご了承ください m(_ _)m
4. 本日の流れ 1. 従来の頻度主義とベイズ主義の違い 2. ベイズファクターとは？ 3. R による簡単な実験 復習
5. 伝統的な統計学 （頻度論） ベイズ統計学 仮説検定 パラメータ 𝜽 （母数） データ 𝑿 （標本） 1. 従来の頻度論とベイズ統計学の違い • 未知の定数 • 真値は１つ • データから 最尤法で推定 • 母集団から得られる 標本 • 確率変数 • 観測された定数 • 確率変数 • 事前分布 事後分布 （ベイズの定理） • 帰無仮説 H0 と対立仮説 H1 • 帰無仮説 H0 の元で データが得られる確率 を検証する • 「帰無仮説が正しくない」 ことを示す手法 H1 ? H0 ? どーするの？ 復習
6. 頻度論における 95% 信頼区間 何度も同じサンプルサイズの標本データを取ると， 真値が95%の確率で信頼区間内に⼊る 母集団 （真値：𝜃0 ） 標本データは確率変数 𝜃 母集団の真値： 𝜽𝟎 例えば， 100 個の標本データによる 信頼区間があった時， 5 回真値が⼊らない 真値が⼊っていない例
7. ベイズ統計学における 95% 信用区間 MCMC などでサンプリングして求めた事後分布は，母数の確率分布になっている． そのため，「95% の確率でその範囲内に真値がある」ということができる． 𝒇 𝜽|𝑿 = 𝒇 𝑿|𝜽 𝒇 𝜽 𝒇 𝑿 事後分布 事前分布 尤度 周辺尤度 （基準化定数，エビデンス） 事後分布の形状まで求めている （複雑な分布でも MCMC で求まる） 𝜃 ベイズ信用区間 MAP 推定値
8. 2. ベイズファクターとは？  ベイズ統計学において仮説の評価を行うことができる （Jeffreys, 1935. Hoijtink, Klugkist, and Boelen, 2008. Hoijtink, 2011.）  母数に対する仮説を「直接」評価 頻度論の場合は，仮説の元で標本が得られる確率を評価  頻度論における仮説検定とは似て非なるもの
9. 二群の平均に関して以下の仮説を例として考える 𝑯𝒊 : 𝝁𝟏 ≥ 𝝁𝟐 上式のように， 母数に対して不等式制約を用いて表される仮説を 「情報仮説」と呼ぶ ベイズファクターを簡単な例で考えてみる informative hypothesis
10. 情報仮説に対して，母数に制約を課さない仮説を 「無制約仮説」と呼び， この場合，ベイズファクターは 情報仮説と無情報仮説で 計算する 𝑯𝒖 : 𝝁𝟏 𝝁𝟐 𝑯𝒊 : 𝝁𝟏 ≥ 𝝁𝟐 一方で・・・ 前頁の「情報仮説」 unconstrained hypothesis
11. では，𝐻𝑖 と 𝐻𝑢 におけるベイズファクターとは？ 𝑩𝑭𝒊𝒖 = 𝐩𝐨𝐬𝐭𝐞𝐫𝐢𝐨𝐫 𝐨𝐝𝐝𝐬 𝐩𝐫𝐢𝐨𝐫 𝐨𝐝𝐝𝐬 = 𝑷 𝑯𝒊 |𝑿 𝑷 𝑯𝒖 |𝑿 𝑷 𝑯𝒊 𝑷 𝑯𝒖 観察データ 𝑿 によって 事前確率の比 が 事後確率の比 へとどれだけ変化したかを計算 事前確率の比（事前オッズ） （事前確率とは， 任意の仮説において事前分布を仮説に与する母数空間で積分したもの） （事後確率とは， 任意の仮説において事後分布を仮説に与する母数空間で積分したもの） 事後確率の比（事後オッズ）
12. 事前オッズ 𝒇 𝜽|𝑿 = 𝒇 𝑿|𝜽 𝒇 𝜽 𝒇 𝑿 （ベイズの定理） 𝑷 𝑯𝒖 = 𝒇 𝝁𝟏 𝒇 𝝁𝟐 𝒅𝝁𝟏 𝒅𝝁𝟐 = 𝟏 𝑷 𝑯𝒊 = 𝝁𝟏≥𝝁𝟐 𝒇 𝝁𝟏 𝒇 𝝁𝟐 𝒅𝝁𝟏 𝒅𝝁𝟐 = 𝟏 𝟐 𝑷 𝑯𝒊 𝑷 𝑯𝒖 = 𝟏 𝟐 簡単のため，互いに独立な事前分布を仮定 互いに独立な正規分布: 𝑁 0,𝐷−1 （等高線は確率密度の大きさ） 𝜇1 𝜇2 𝑯𝒖 𝑯𝒊
13. 事後オッズ 𝜇1 𝜇2 𝑯𝒖 𝑯𝒊 (a) (b) 𝒇 𝜽|𝑿 = 𝒇 𝑿|𝜽 𝒇 𝜽 𝒇 𝑿 （ベイズの定理） 𝑷 𝑯𝒖 |𝑿 = 𝒇 𝝁𝟏 , 𝝁𝟐 |𝑿 𝒅𝝁𝟏 𝒅𝝁𝟐 = 𝟏 𝑷 𝑯𝒊 |𝑿 = 𝝁𝟏≥𝝁𝟐 𝒇 𝝁𝟏 , 𝝁𝟐 |𝑿 𝒅𝝁𝟏 𝒅𝝁𝟐 データ 𝑋 によって事前分布が更新された結果， （a）正の相関 もしくは（b）負の相関 をもち， かつ，平均が（0, 0）ではないような事後分布が 得られた二つのケース（a）と （b）を考える 𝑷 𝑯𝒊 |𝑿 ≫ 𝟏 𝟐 𝑷 𝑯𝒊|𝑿 𝑷 𝑯𝒖|𝑿 ≫ 𝟏 𝟐 𝑩𝑭𝒊𝒖 ≫ 𝟏 (a) (b) 𝑷 𝑯𝒊 |𝑿 ≪ 𝟏 𝟐 𝑷 𝑯𝒊|𝑿 𝑷 𝑯𝒖|𝑿 ≪ 𝟏 𝟐 𝑩𝑭𝒊𝒖 ≪ 𝟏
14. 岡田 2013 より ベイズファクターに対する基準  Jeffreys や Kass らの基準が よく知られている  あくまでも大雑把な経験則で 絶対的なものではない点に留意  p 値と似たような感覚で扱うと 良いとのこと （Rosnow & Rosenthal, 1989）
15. 今回は，無制約仮説に対する情報仮説の元で ベイズファクターを考えたが， 無制約仮説ではなく相補仮説（情報仮説の補空間）の元で 計算されることも増えてきている 詳しくは 岡田, ベイズ推定による情報仮説の評価：その理論と各種モデルへの応用について 専修人間科学論集心理学篇, 2016, 6, 9-17 https://core.ac.uk/download/pdf/71799133.pdf など
16. 3. R による簡単な実験 岡田, ベイズ統計による情報仮説の評価は分散分析にとって代わるのか? 基礎心理学研究, 2013, 32巻, 2号, p.223-231 https://www.jstage.jst.go.jp/article/psychono/32/2/32_KJ00009351488/_pdf/-char/ja ここでの実験は以下の内容の再現実験 但し，OpenBUGS ではなく RStan を使用 （一部，事前分布の設定も異なる）
17. 使用するデータ 27 匹の鼠を 3 つのグループ LD: 通常の明暗サイクル LL: 常に明るい光をつけた状態 DM: 日中は明るく，夜は薄暗い状態 に分けて， 生活（食住つき）させた実験 結論は， 「夜間の光は体重を増加させる」だが これをベイズファクターで検証してみる 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 を無制約仮説 𝐻𝑢 に対して検証 𝑯𝟏 : 𝝁𝑳𝑫 < 𝝁𝑫𝑴 < 𝝁𝑳𝑳 𝑯𝟐 : 𝝁𝑳𝑫 < 𝝁𝑫𝑴, 𝝁𝑳𝑳 𝑯𝒖 : 𝝁𝑳𝑫, 𝝁𝑫𝑴, 𝝁𝑳𝑳
19. 𝑩𝑭𝒊𝒖 = 𝐩𝐨𝐬𝐭𝐞𝐫𝐢𝐨𝐫 𝐨𝐝𝐝𝐬 𝐩𝐫𝐢𝐨𝐫 𝐨𝐝𝐝𝐬 = 𝑷 𝑯𝒊|𝑿 𝑷 𝑯𝒖|𝑿 𝑷 𝑯𝒊 𝑷 𝑯𝒖 = 𝑷 𝑯𝒊|𝑿 𝑷 𝑯𝒊 無制約仮説のもとでは だったので， 𝑃 𝐻𝑖 と 𝑃 𝐻𝑖 |𝑋 を計算すればよい 𝑯𝟏 : 𝝁𝑳𝑫 < 𝝁𝑫𝑴 < 𝝁𝑳𝑳 𝑯𝟐 : 𝝁𝑳𝑫 < 𝝁𝑫𝑴, 𝝁𝑳𝑳 𝑷 𝑯𝟏 = 𝟏 𝟔 𝑷 𝑯𝟐 = 𝟏 𝟑 全母数空間で事前分布を 積分するので，事前確率は 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]); } 𝑃 𝐻𝑖 |𝑋 は RStan のサンプリングで求める サンプリングされた事後分布のうち， 各情報仮説の条件に与する割合を計算し推定する （本実験の場合，サンプリングされた 4000点のうちどれだけの割合が条件を満たすかを計算） 事前分布は， 平均は各群毎に弱情報事前分布を 分散は群に依らず同一の弱情報事前分布を設定※ ※ 岡田 2013 の設定から大きく変更している．岡田 2013 の設定だと収束しない（そもそも弱情報事前分布になっていなさそう）が，恐らくは歳月を経てデータセットの次元が変わったため？（未確認） f1 で情報仮説 H 1 に与するサンプリングかどうかを計算 H 1 の条件にあてはまる時は 1 に，そうでない時は 0 になる 同様に f2 は H 2 の条件にあてはまる時は 1 に，そうでない時は 0 になる 𝑯𝟏 : 𝝁𝑳𝑫 < 𝝁𝑫𝑴 < 𝝁𝑳𝑳 𝑯𝟐 : 𝝁𝑳𝑫 < 𝝁𝑫𝑴, 𝝁𝑳𝑳
21. bayes_factor.R ( from L84 to L93 ) # 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 Reference: Okada, 2013  事前分布の設定などが異なるが 岡田（2013）とほぼ同じ結果  Kass & Raftery の基準に従うと どちらの仮説も "Positive" となる  仮説としては H 1 の方がエビデンスレベルが強い
22. 今回使ったコード BF_model1.stan stan model の定義 bayes_factor.R MCMC の実行・BF の計算
23. Reference 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
24. Thank you! Any questions