この記事は，Stan Advent Calendar2019の6日目の記事です。
brmsパッケージ便利ですよね。みなさん使ってますか？
今回の記事ではそのbrmsパッケージで「安易に」ベイズファクターを計算すると思いも寄らない結論を下してしまう可能性を指摘します。実際、死にかけました（意味深
ただ、使い方を間違えなければbrmsでベイズファクターを計算すること自体は問題ないので、安心してくださいね。
brmsパッケージ
brmsパッケージをどのように発音するかで戦争が起こるという噂を聞いたことはありませんが、僕は「びーあーるえむえす」と読んでます。まさか「ぶらむす」なんて読んでませんよね？
さて、brmsはBayesian Regression Models using 'Stan'、の略称です。つまり、Stanを使ってベイジアンな回帰分析ができるパッケージです。 しかし、最近では回帰分析なんてちゃちなレベルじゃなくて、GLMMはもちろん、SEMやDiffusionモデルとか、心理系で使いそうなほとんどな分析が簡単なコードを書くだけでできてしまいます。すげぇ！
brmsの詳しい説明は、das_kinoさんのこちらの記事に超わかりやすく書いてるので、そちらを見てください。 ついでに、僕はStanコードを書いてると心が癒やされるタイプの人間なんで、あんまりbrmsは使ってないです（おい
ベイズファクター
ベイジアンのみなさんの多くはご存知だと思いますが、ベイズファクターとは、モデルを比較するときの指標の一種です。定義としてはそれぞれのモデルの周辺尤度の比、ということになります。
周辺尤度は、想定した確率モデル（ここでは尤度と事前分布の両方を指す）から、今回得られたデータが得られる確率です。データが真の分布から生成されたもとの考えれば、想定した確率モデルが真の分布の近似として適切であれば、周辺尤度は高くなります。
その周辺尤度の比、ということですから、どちらのモデルがデータ（あるいは真の分布）から見て適切なモデルか、ということを表しています。 ベイズファクターや周辺尤度は、汎化損失の推定値であるWAICなどに比べて、事前分布の影響を強く受けます。それはWAICなどの予測精度を表す指標は尤度を事後分布で平均した（事後）予測分布を評価の対象にしているからです。
尤度を事前分布で平均した確率モデルを評価の対象としているベイズファクターは、事前分布をどのように設定するかで、結構変わってきます。 今回の話は、このベイズファクターと事前分布の話だと思ってもらえればOKです。
bridgesamplingパッケージとベイズファクター
さて、ベイズファクターは周辺尤度を使って計算するのですが、実はこの周辺尤度の計算が難しいのです。なんせ、事前分布で積分しないといけないからです。
しかし、最近この周辺尤度をシミュレーションで計算できる方法が実用化されました。ブリッジサンプリングと呼ばれる方法です。ブリッジサンプリングについては東大の岡田さんの論文がわかりやすいです。
そして、それがRに実装されたのが、bridgesamplingパッケージです。
birdgesamplingは結構いい精度で周辺尤度を推定できます。簡単なシミュレーションをしたものに「社会科学のためのベイズ統計モデリング」の7章があります。興味ある人はぜひ。販促販促。
実はbrmsは内部にbridgesamplingパッケージを依存パッケージとして含んでいます。なので、brmsはbridgesamplingパッケージを使って、モデルの周辺尤度が計算できてしまいます。
便利ー
brmsでベイズファクターを計算してみる
使い方はやはり上で紹介した、das_kinoさんのページに書いています。 注意点として、次のような記述が。
さらに以下3点を対処する必要があります。
- 全パラメータに、無情報でない事前分布を指定
iterを多めにsave_all_pars = TRUEと引数を指定
まず全パラメータに無情報でない事前分布を指定とあります。
周辺尤度は事前分布で平均した尤度ですので、事前分布がimproperだとうまく計算ができないんですね。なので事前分布を設定しましょう、というわけです。
次に、bridgesamplingパッケージで周辺尤度を計算するためには、MCMCサンプルが大きくないと、良い精度で推定できません。それもモデルの複雑さに依存するんですが、50000くらいあると安心みたいです。おそらくESSが大きくないとダメっぽいので、単純にiterを増やせばいいというわけでもないと思いますが。
最後はbrmsの都合ですが、brmsは分析結果を軽くするために内部で作ったパラメータを全部保存するわけじゃないんですね。でもbridgesamplingのためにはそれが全部いるので、オプションを指定しようね、というわけです。
でね、何が問題かというと・・・最初のやつです。事前分布。
das_kinoさんの記事にちゃんと書いてあったんですけど、僕は当初brmsだし、回帰係数にも適切な弱情報事前分布が設定してあると思ってたんです（僕のかすか記憶では前のバージョンでは設定していたと思う）。 それがね、してないんですよ。brms。
だから、回帰係数に事前分布を設定しないとbridgesamplingがうまく行かないんです。ほえー。
実はこれ、das_kinoさんとの共同研究をしているとき、指摘してもらって気づいたんですよ。で事前分布を設定してやってみたら、あ、結果結構違うって。論文投稿する前に気づいてよかったー。感謝感謝。
で、それがどううまく行かないかを以下でお伝えします。
brmsで安易にベイズファクターを計算すると死ぬ話
ついにタイトルに来ました。
まずは、適当な回帰分析をやります。データを適当に生成します。 今回は、真のモデルが2つ説明変数で構成され、あとは正規ノイズだとします。なので、3つ目の説明変数を入れるとダメ、という設定です。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
#パッケージ library(rstan) library(brms) library(bridgesampling) library(tidyverse) library(magrittr) #パラメータの真値 alpha <- 2 beta <- c(3,4,0) #乱数でデータ生成 set.seed(123) x1=rnorm(100,0,1) x2=rnorm(100,0,1) x3 = rnorm(100,0,1) y <- alpha + beta[1]*x1 + beta[2]*x2 + beta[3]*x3 + rnorm(100,0,10) dat <- data.frame(y=y,x1=x1,x2=x2,x3=x3) |
モデル生成では3つの説明変数があるように見えますが、回帰係数の3つ目が0になってるので、結局2つで説明しているモデルとなってます。
これをbrmsで推定してみましょう。
最初は無情報事前分布、つまりbrmsのデフォルトの設定です。
|
1 2 3 |
#brm関数で推定する result0 <- brm(y ~ x1+x2+x3,data=dat,save_all_pars = TRUE,iter=10000) result1 <- brm(y ~ x1+x2,data=dat,save_all_pars = TRUE,iter=10000) |
まずはsummaryで結果を見てみましょう。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
> summary(result0) Family: gaussian Links: mu = identity; sigma = identity Formula: y ~ x1 + x2 + x3 Data: dat (Number of observations: 100) Samples: 4 chains, each with iter = 10000; warmup = 5000; thin = 1; total post-warmup samples = 20000 Population-Level Effects: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat Intercept 1.81 1.08 -0.34 3.94 23244 1.00 x1 2.43 1.18 0.14 4.77 23580 1.00 x2 4.45 1.10 2.27 6.60 24454 1.00 x3 -0.58 1.14 -2.82 1.67 23529 1.00 Family Specific Parameters: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat sigma 10.60 0.78 9.22 12.25 21491 1.00 Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat = 1). |
当たり前ではありますが、x3の回帰係数は0近くに推定され、95%CIに0を含みます。
さて、次はベイズファクター。
今回は真のモデルがresult1なので、ベイズファクタはresult１を基準に計算すれば、１より大きく（理想的には３以上）なるはずです。
では、実際に計算してみましょう。
|
1 2 |
#ベイズファクターの計算 brms::bayes_factor(result1,result0) |
計算結果は以下です。
|
1 |
Estimated Bayes factor in favor of bridge1 over bridge2: 0.31121 |
はてはて？
BF=0.31ということは、真のモデルではない、冗長なモデルが積極的に支持されるレベル（逆数を取ると3.22になる）です。
では、事前分布を適切に設定した場合のモデルです。
適切な事前分布って、なんなんでしょうか。一つの選択肢として、Jeffreysの事前分布なんかがあります。重回帰分析の場合、この論文では、Jeffreysの事前分布となるコーシー分布を使うのじゃ、と言ってます。ついでにこの論文の第一著者のRouder JNさん、ファーストネームがJeffreyなんですね。紛らわしいわ！（失礼
さて、コーシ分布は位置パラメータと尺度パラメータの二つをもちます。大抵の事前分布がそうであるように、位置パラメータは0でいいんですが、尺度パラメータはどうするか。上の論文ではモデルの残差SDを使えとあります。そんなん事前にわかるのかと思いますが、まぁ保守的にいくなら目的変数のSDを使えばいいんじゃないでしょうか。
でも今回は真のモデルがわかってるので、素直に尺度パラメータ＝10としておきます。
|
1 2 3 4 5 6 7 8 9 10 |
#事前分布を設定 result0_2 <- brm(y ~ x1+x2+x3,data=dat,save_all_pars = TRUE, prior = prior(cauchy(0,10),class=b), iter=10000) result1_2 <- brm(y ~ x1+x2,data=dat,save_all_pars = TRUE, prior = prior(cauchy(0,10),class=b), iter=10000) #ベイズファクターの計算 brms::bayes_factor(result1_2,result0_2) |
さってさて、どうなるかな？
|
1 |
Estimated Bayes factor in favor of bridge1 over bridge2: 9.90376 |
あれれーさっきとぜんぜん違うよー（コナン的な声で
ふむ。一応、真のモデルはモデル1のほうなので、こっちのほうが正しいはず。
念の為、最尤法で計算してBICからベイズファクターの近似値を出してみましょう。これはあくまで近似値なので、厳密には正しくないんですが・・・
|
1 2 3 4 5 |
bic0 <- lm(y ~ x1+x2+x3,data=dat) %>% BIC bic1 <- lm(y ~ x1+x2,data=dat) %>% BIC exp((bic0-bic1)/2) exp((bic0-bic1)/2) [1] 8.728503 |
うむー！
というわけで、無情報事前分布の結果がかなりあやしい、ということはわかるかと思います。
シミュレーションしてみるぞ
今回の結果は特定のデータだけの結果かもしれないので、一応シミュレーションをしてみます。
brmsでは毎回コンパイルしよるので、brmsで生成したStanコードと同じコードを使ってシミュレーションを回します。また、ベイズファクターもbridgesamplingから直接推定してます。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 |
#stan model model <- stan_model("brms_reg.stan") #setting N <- 100 alpha <- 2 beta <- c(3,4,0) sigma <- 10 beta_var <- 10 #non-informative prior prior_check <- 0 #start set.seed(1234) x1 <- rnorm(N,0,1) x2 <- rnorm(N,0,1) x3 <- rnorm(N,0,1) bs0_1 <- c() bs1_1 <- c() for(i in 1:100){ y <- alpha + beta[1]*x1 + beta[2]*x2 + beta[3]*x3 + rnorm(N,0,sigma) dat <- data.frame(y=y,x1=x1,x2=x2,x3=x3) X <- model.matrix(y ~ x1+x2+x3,data=dat) datastan <- list(N=N,Y=y, K=4, X=X, beta_var = beta_var, prior_check = prior_check) fit <- sampling(model,data=datastan,iter=2100, warmup=100,chains=4) bs0_1[i] <- bridge_sampler(fit) X <- model.matrix(y ~ x1+x2,data=dat) datastan <- list(N=N,Y=y, K=3, X=X, beta_var = beta_var, prior_check = prior_check) fit <- sampling(model,data=datastan,iter=2100, warmup=100,chains=4) bs1_1[i] <- bridge_sampler(fit) } #weakly informative prior prior_check <- 1 #start bs0_2 <- c() bs1_2 <- c() for(i in 1:100){ y <- alpha + beta[1]*x1 + beta[2]*x2 + beta[3]*x3 + rnorm(N,0,sigma) dat <- data.frame(y=y,x1=x1,x2=x2,x3=x3) X <- model.matrix(y ~ x1+x2+x3,data=dat) datastan <- list(N=N,Y=y, K=4, X=X, beta_var = beta_var, prior_check = prior_check) fit <- sampling(model,data=datastan,iter=2100, warmup=100,chains=4) bs0_2[i] <- bridge_sampler(fit) X <- model.matrix(y ~ x1+x2,data=dat) datastan <- list(N=N,Y=y, K=3, X=X, beta_var = beta_var, prior_check = prior_check) fit <- sampling(model,data=datastan,iter=2100, warmup=100,chains=4) bs1_2[i] <- bridge_sampler(fit) } #bayes factor (bs1_1-bs0_1) %>% exp %>% median (bs1_2-bs0_2) %>% exp %>% median |
上のシミュレーションでは、prior_checkを0にすれば無情報、 1にすればコーシー分布となるようになるようにしてます。
さて、結果は。
|
1 2 3 4 |
> (bs1_1-bs0_1) %>% exp %>% median [1] 0.2971951 > (bs1_2-bs0_2) %>% exp %>% median [1] 9.518806 |
というわけで、100回程度回しただけではありますが、無情報事前分布では真のモデルよりも冗長なモデルが明らかに選択されてしまいました。そして、弱情報事前分布を仮定したほうは、ちゃんと真のモデルが選択されていました。
今回の話の結論
というわけで、「ベイズファクターを使うときは、無情報事前分布はやめような！」という話。ベイズファクターを使うときは「そんな事前分布で大丈夫か？」と常に確認しましょう。
特に、brmsはデフォルトでは回帰係数に事前分布が設定されていないので要注意。弱情報事前分布として、とりあえずコーシー分布をおいてみよう。尺度パラメータは、目的変数のSDぐらいにしてみよう。
ただ、実際は感度分析といって、事前分布のパラメータを色々かえて、結果がどのように変わるのかを見るのが重要。今回の場合、コーシー分布の尺度パラメータを1にすると結構結果が変わります。
|
1 2 3 4 5 6 |
#コーシー分布の尺度パラメータ＝1 result0_3 <- brm(y ~ x1+x2+x3,data=dat,save_all_pars = TRUE, prior = prior(cauchy(0,1),class=b)) result1_3 <- brm(y ~ x1+x2,data=dat,save_all_pars = TRUE, prior = prior(cauchy(0,1),class=b)) brms::bayes_factor(brms::bayes_factor(result1_3,result0_3)) |
|
1 2 |
Estimated Bayes factor in favor of bridge1 over bridge2: 1.63458 |
BF=1.63というわけで、どっちも積極的に採用できない、という結果に。
今度は逆に、尺度パラメータ＝100とかにすると、
|
1 |
Estimated Bayes factor in favor of bridge1 over bridge2: 97.04573 |
BF=97となります。
あー怖いですねー
僕の印象として、BF>3だったらこっちを採用！とかそんな甘い基準だと死ぬな、という感じです。そんな小さい値でモデルの優劣を決めようと思っていると、事前分布の設定で一瞬でひっくりかえるからです。再現性とは、みたいな。
まぁBF>20兆とかなら、さすがに大丈夫だとは思いますが。
そして、そもそも、パラメータ＝0に対するモデルと自由推定したモデルのBFは結構不安定です。関数系やモデルが違ってて、パラメータの数は同じ、というときはやりやすいんですけどね。
以上です。Enjoy brms!