brmsパッケージを用いたモデリングの方法や、便利機能、ちょっとしたtipsを紹介します。読み方はそのまま「ビー・アール・エム・エス」です(僕は勝手にブラムスと呼んでいますが)。その名の通り、Stanのラッパーパッケージですが、ユーザは自身でStanコードを書く必要はなく、モデルを指定すると自動的に内部でStanコードが生成され、サンプリングが実行されます。最近は、brmsを用いて分析した論文も見かけるようになってきました。Emotion僕もたまに使っているのですが、あまりにも多機能なので、しっかりと把握できていない部分も多いです。本記事に間違いがあればご指摘いただければ幸いです。brmsパッケージの良いところその1回帰モデルと一言でいっても、実はその守備範囲はかなり広いです。線形モデルも非線形モデルも、その気になればSEMもできるそうです(それならblavaanでいいと思いますが)。一般線形モデル一般化線形モデル階層線形モデルその他の自由なモデリング嬉しいのは、上記の1~3は、多くのRユーザが慣れ親しんでいるであろう記法と同じ記法でモデルを指定できることです。dat)のように書きますよね。brmsパッケージでも、記法はほぼ同じになります。その2デフォルトで用意されている確率分布が、ものっっっっっっっすごく多いです。コンソールで?brms::brmsfamilyと入力してみてください。用意されている確率分布の一覧が確認できます。驚くのは、「ベルヌーイ分布とポアソン分布の混合分布である、ゼロ過剰ポアソン分布」等、異なる確率分布を組み合わせた、いわゆる混合分布も多く用意されていることです。通常、このような混合分布を用いたモデリングをStanで書くにはひと手間要るのですが、brmsパッケージでは確率分布を指定するだけで済みます。様々な確率分布からの乱数を生成するための関数や...確率(密度)を返す関数も用意されています。その3「内部で何をやっているのか」を知ることができます。ベイズを用いる際に最も気になることの一つが、事前分布の設定でしょう。brmsでは、brms::get_prior()で指定したモデルの事前分布を知ることができます。brmsがデフォルトで採用している事前分布も分かるので、報告時に困ることはありません。brms::make_stancode()で指定したモデルのStanコードを出力させることもできます。「この確率分布を使ったモデルをStanで書きたいとき、どうすればいいんだろう」と悩んだとき、brmsで同じ確率分布を使った適当なモデルでStanコードを出力させて、「なるほど、こう書けばいいのか」とヒントを得ることも可能です。僕が驚いたのが、brms流のガンマ回帰の書き方です。parametersブロックと、modelブロックの一部を以下に抜粋します(事前分布のコードは省略しています)。modelブロックでは、まず切片と説明変数の線形結合をmuに代入しています。ガンマ分布にはshapeとrate、2つのパラメータがあり、平均はshape/rateになります。shape/rateの関係があります。exp(-mu)になります。ところがparametersブロックではrateに相当するパラメータはそもそも宣言しておらず、一度別の用途で使ったmuという変数を再利用して上書きしています。「もともとのmuには興味ないでしょ」と言わんばかりに。そのmuはmodelブロック内で宣言されているため、最終的にrate(に相当する)パラメータを取り出すことができません。思い切ったことをするなあと思ったのですが、こんな風にbrmsはサンプリング効率を高めるようなStanコードの書き方をしているため、推定がかなり速いです。その4他のパッケージとの連携が充実しています。brmsはbayesplotパッケージを利用して、brms::pp_check()という関数により、容易に事後予測チェックを行うことができます(詳しくはこちら)。factorを容易に計算できたり、欠測値がある場合にmiceパッケージで多重代入法による補完を行ったあとで、補完したそれぞれのデータセットからサンプリングを行えたりします(詳しくはこちらから)。brmsパッケージの使用例それではさっそく、使用例を紹介していきましょう。ロードbrmsはCRANに登録されているので、以下でOKです。本稿執筆時点における最新バージョンは2.6.0です。install.packages("brms")library(brms)上述のbayesplotパッケージやbridgesamplingパッケージ、それからStanをRから扱うためのrstanパッケージを含む、種々の依存パッケージが一緒にインストールされます。(Rtoolsは別途インストールする必要があります)。2018/12/19追記:記事を書いた直後にバージョン2.7.0にアップデートされました。processes"ですって...一般線形モデル(回帰分析)まずはベーシックな回帰分析をやってみましょう。サンプルデータは、プリセットのmtcarsです。32台の車の性能が収められたデータセットになります。燃費(mpg)を、車体重量(wt)とマニュアル車かオートマ車か(am)で予測する重回帰分析を試みます。mtcars$amには0か1が格納されていますが、今回はこれらをfactor型に変換したうえで、交互作用項を含むモデルを指定することにします。mtcars)summary(fit)brmsでは、どのようなモデルであっても、統一的にbrm()という関数のなかでモデルを指定します。rstan::sampling()でサンプリングする際に指定するような諸々の引数も、brm()は対応しています。サンプリング回数(iter)、ウォームアップ期間(warmup)、乱数のシード(seed)、チェイン数(chain)など。1).また、サンプリング結果を格納したオブジェクトをsummary()に入れると要約された結果が返されますが、brms::waic()やbrms::loo()に入れると、WAICやlooを計算することもできます(内部ではlooパッケージを使用して計算しています)。details.事前分布の表示brmsは、いくつかのパラメータについてデフォルトで事前分布を与えていることがあります。今回指定したモデルで、各パラメータにどのような事前分布が与えられていたのかを知るには、brms::get_prior()を使用します。この中にモデルを指定すると、以下のように一覧が表示されます。10)、19ってどういう根拠で設定されたんでしょうね...。brmsはこんな風に、いくつかのパラメータにデフォルトで弱情報の事前分布を与えていることがあるのですが、どういう根拠に基づいているのかは、まだよく分かっていません。brmsでは基本的に(おそらく常に?)、回帰係数bには事前分布は明示されません。デフォルトでは、無情報事前分布が適用されるということになります。factorの計算」においては問題となるため、そのような場合にはユーザが事前分布を指定する必要がでてきます。=という引数にパラメータと事前分布をセットで指定します。指定方法には複数の記法があります。詳しくは、コンソール上で?brms::set_priorと入力して、ヘルプをご参照ください。)Stanコードの表示brms::make_stancode()のなかに、事前分布を含み任意のモデルを指定すると、そのStanコードをコンソールに表示してくれます。先ほど事前分布を指定した、以下のモデルのStanコードを表示してみましょう。長いので、ここではmodelブロックだけ掲載します。+=記法を用いて書かれています。結果の可視化plot()brmsは、サンプリング結果の視覚化にも長けています。rstan::stan_trace()やrstan::stan_dens()のように、描きたいグラフに応じて個別の関数が用意されていることが多いと思いますが、brmsはbayesplotパッケージのbayesplot::mcmc_combo()を用いて結果を可視化します。2018の2日目の記事や、11日目の記事をご覧ください。bayesplot::mcmc_combo()は、各パラメータについて、事後分布とトレースプロットを並べて表示してくれます。plot(fit)繰り返しますがこのグラフはbayesplotパッケージを用いて出力されているので、同パッケージの関数を用いて見栄え(Aesthetic)を調整することができます。2日目の記事を参照してください)。mcmc_combo()で表示できるグラフの種類は限られているようです。こんな風に、色調をピンクに変更し、トレースプロット(bayesplot::mcmc_trace())とチェインごとの事後分布(bayesplot::mcmc_dens_overlay())を並べることもできます。2.11.0より、brms::marginal_effect()がbrms::conditional_effects()に名称変更されています。詳しくはこちらさらにbrms::marginal_effects()を用いると、「主効果」や「交互作用」を可視化することもできます。今回は交互作用項を含むモデルを指定したので、いわゆる「単純傾斜」も自動的に可視化されます。下のグラフは、左から順に「wtの主効果」「amの主効果」「wtとamの交互作用」を示します。引数effectsに何も指定しなければ、すべての項について可視化が行われますが、特定の項だけを指定することも可能です。交互作用項を指定する場合には、その順番を変えることで、どちらの変数についての単純効果を可視化するかも決めることができます。#右のグラフ離散変数の場合は(今回はam)、各水準における他方の変数の単純効果が可視化されます。連続変数の場合は(今回はwt)、平均±1SDのポイントにおける、他方の変数の単純効果が示されます。11日目の記事も参照してください。quantitiesブロックでモデルから乱数を生成させる手続きが必要になります。11日目の記事を読んだ方はお察しの通り、これは内部でbayesplot::ppc_dens_overlay()を使用しています。他の事後予測チェックにも対応していますよ。=という引数に、描画形式を指定してみます。"error_hist")これはbayesplot::ppc_error_hist()を使用していますね。一般化線形モデル次は、正規分布以外の指数分布族を扱えるように拡張した、一般化線形モデルを例に挙げます。今度は目的変数を、オートマ車かマニュアル車かamにし、説明変数に燃費mpgと重量wtを投入してみます。choice.とサジェストされます)。"binomial"を指定した際の記法がstats::glm()と異なるので、ここでは紹介のためあえて二項分布を適用した書き方をしてみます。stats::glm()ではこのように書きます。am)の列を指定しています。brmsでは、目的変数が二項分布に従うと仮定するとき、以下のように書きます。|の左側には成功試行数を、右側にはtrials()の中に試行数を入れます。OKsummary(fit)brms::marginal_effects()で主効果を可視化してみると、ロジットリンク関数によって、直線ではなく曲線になっていることがわかりますね。familyゼロ過剰ポアソン分布さて、ここからがbrmsの真骨頂です。stats::glm()ではfamilyに指定できる確率分布は以下の通りです。gaussianbinomialGammainverse.gaussianpoissonquasiが、brmsでは豊富なfamilyが用意されているということを、冒頭で述べました。試しにちょっと珍しい分布を使ってみましょう。Mr.Unadon氏がTwitter上で「あなたは、片方しかない靴下を何足もっていますか?」とアンケートをとった結果を拝借します。ネタ元はこちらの片方しかない靴下の数をゼロ過剰ポアソン分布でモデリングしてみたです。10)YがUnadon氏が集めた、片方しかない靴下の所持数のデータです。今回のサンプルの中では、ほとんどの人は靴下を片方だけなくすということはないようですね。僕が初めてこのアンケートを見た時、「靴下片方だけどうやって失くすねん」と思ったのですが、最近我が家で片方しかない靴下が発見されました。靴下の所持数はカウントデータなので、ポアソン分布と相性がよさそうに思えますが、上図のとおりゼロが圧倒的に多く、通常のポアソン分布の仮定を満たしそうにありません。poisson)という確率分布を考えます。ゼロ過剰ポアソン分布は、ベルヌーイ分布とポアソン分布の混合分布です(以下の図はUnadon氏の記事から拝借しました)。靴下を失くすか否かがベルヌーイ分布に従うとと考え、もし失くさないなら当然、片方しかない靴下の所持数の回答は0。失くすなら、所持数はポアソン分布に従う、と考えます。ポアソン分布に従う確率変数は0以上の整数をとりうることに注意してください。出典:片方しかない靴下の数をゼロ過剰ポアソン分布でモデリングしてみた一方Xは適当に作った変数ですが、仮に「日常の忙しさ指標(0-100点)」だとしておきましょう。仮定としては、忙しい人ほど靴下を片方なくしやすい、といった具合です。片方しかない靴下の所持数を、忙しさで予測する、シンプルなモデルを考えてみましょう。+=という記法を用いて、少し工夫した書き方が必要になります。"zero_inflated_poisson"を指定するだけです。1).事後予測チェックしてみましょう。brms::pp_check(fit)いい感じですね。階層線形モデル本記事冒頭でも言及したように、lme4パッケージやlmerTestパッケ―ジを用いると、以下のように書きますよね。dat)brmsでは、関数をbrm()に変えるだけなので、本記事では説明を省略します。モデル比較brmsパッケージを用いてサンプリングした結果を利用して、モデル比較を行ってみます。factorを指標とすることにします。factorについては、北條大樹さんのこちらの記事や、清水裕士さんのスライドをご参照ください。brmsを用いてモデル比較するために、kidneyというデータセットを用いて、試してみましょう。survivalパッケージにも同じ名前のデータセットがあるので、区別するためにbrms::kidneyと指定したほうが無難です。話を簡単にするために、ここでは目的変数に「病気の再発までの時間time」、説明変数に「患者の年齢age」を投入した、単回帰分析を考えます。本当は右側打ち切り(right-censored)が起こっているため、それを考慮したモデリングをするべきですが、打ち切りが起こっていないデータだけを抽出しています。dat)のように書くことで、打ち切りデータのモデリングも容易に実現できます(cens()のなかに、打ち切りの有無を判別する変数を投入します)。other目的変数timeの分布はこのように左に歪んでいます。正規分布を仮定したモデルまずはベーシックな単回帰分析を実行してみます。"normal"と指定します(familyを指定しなければ自動的に正規分布を仮定することになるので、省略も可能です)。以下3点を対処する必要があります。TRUEと引数を指定モデルはこのようになります。TRUE)事後予測チェックを行ってみましょう。brms::pp_check(fit_n)実際に観測されたデータとズレがありますね。上のヒストグラムを見ればわかる通り、今回はtimeが正規分布に従うと仮定するのは、適切ではないかもしれません。samplingによる対数周辺尤度を計算してみましょう。これは内部ではbridgesampling::bridge_sampler()を利用しています。"normal".ワイブル分布を仮定したモデルbrms::kidneyは生存時間解析に適したデータセットなので、比較対象として、パラメトリックな生存時間解析でよく用いられる、ワイブル分布を適用したモデルを考えてみます。TRUE)事後予測チェック。brms::pp_check(fit_w)先ほどよりマシになった気がしますね。samplingによる対数周辺尤度の計算。ワイブル分布いよいよモデル比較です。各モデルの対数周辺尤度に、いったん-1を掛けます。自由エネルギーとして計算します。factorになります。58278089今回のデータに対する当てはまりの観点では、正規分布を仮定したモデルに対して、ワイブル分布を仮定したモデルを支持する程度が、58278089倍...もう圧倒的に大きいと考えられます。factorを計算する、brms::bayes_factor()という関数を用意しています。引数x1とx2に、各モデルのサンプリング結果を格納したオブジェクトを指定してください。58273598.05176やはり今回のデータに対する当てはまりの観点では、正規分布を仮定したモデルよりも、ワイブル分布を仮定したモデルのほうが良いと考えられます。こんな風に、brmsを用いてモデル比較をすることもできます。モデル比較という手続きは非常に難しいので、実際にはもっと慎重に行う必要があると思います。以下、北條さんの記事から引用します。ベイズファクターは事前分布の影響をかなり大きくうけます。同一データ、同一モデルで事前分布のみを変化させた場合に、どれだけベイズファクターが変化するかを検討するととてもよくわかります。全てのモデル間のベイズファクターが計算できるわけではありません。モデル比較を行うのはとても難しいです。本来、我々が、比較すべきモデルは無限に存在します。ベイズファクターは、2つのモデル間の比として考えるので、それ以外のモデルについては考えてはいません。6日目の記事:brmsパッケージで安易にベイズファクターを使うと死ぬ話で、brmsでベイズファクターを求める際の注意点が書かれています。"hogehoge"と引数を追加してみましょう(もちろんhogehogeのところは任意の名前で構いません)。"hogehoge")するとカレントディレクトリに、hogehogeという名前のファイルができていると思います。R上でこのファイルを開くと、brms::make_stancode()で生成できるようなStanコードが記録されていることが分かります。仮にbrmsのバージョンが変わったとしても、少なくとも今回のサンプリングを実行した時点でのStanコードは判明しているので、結果の再生可能性の向上に寄与するかもしれません。"fugafuga"と引数を追加してみましょう(もちろんfugafugaのところは任意の名前で構いません)。"fugafuga")するとカレントディレクトリに、fugafuga.rdsという名前でサンプリング結果が保存されています。以下のようにreadRDS()という関数で、あとでロードすることができるので、「しまった、うっかりRStudioを閉じてしまった!」なんて心配も不要です。サンプリング結果をRMarkdownでまとめたい場合などにも重宝すると思います。readRDS("fugafuga.rds")関数名の前に「brms::」と付けるクセをつけるあくまで個人的な意見ですが、何かしらのパッケージ内の関数を使用する場合には、パッケージ名を併記したほうがいいと思っています(例えばdplyr::select()とか)。これはいくつかの理由がありますが、一番の理由は、たまに複数のパッケージ間で同じ名前の関数が存在することがあるからです。filter()という関数は、少なくともstatsパッケージとdplyrパッケージに存在します。他者とコードを共有するときに、何のパッケージの関数か分からないと混乱を招くと思うので、そういった意味でもパッケージ名を明示しておいた方がいいかなと思っています。文字数が増えるのでコードが見にくくなるという側面もありますが。brmsを用いる際には、特に意識してbrms::と書いたほうがよいと思います。これは何故かというと、brmsは様々なパッケージに依存しているのですが、データセット名や関数名が衝突することがあるからです。以下のようなデータセットや関数たち。brms::kidneyデータセットとsurvival::kidneyデータセットbrms::bayes_factor()とbridgesampling::bayes_factor()brms::bridge_sampler()とbridgesampling::bridge_sampler()brms::bf()とbridgesampling::bf()一番タチが悪いのが最後のbf()です。factorとは関係がありません。brms::bf()は、brms::brmsformula()の略称で、brmsを用いて高度なモデリングを可能にするための関数です。分位点回帰を例に説明しましょう。brmsは、分位点回帰も容易に実行することができます。Laplace分布を適用して、以下のように書きます。brm()のなかに、さらにbf()という関数を書いて、その中でモデルを指定します。)これ、bridgesampling::bf()のほうが読み込まれてしまうと、当然エラーになります。したがって以下のように、brms::bf()であることを明示したほうが無難です。)あるいは、パッケージの優先順位をつけてもいいかもしれないですね。#rstatsjおわりにここで紹介した機能は、brmsのほんの一部にすぎません。モデリングの自由度もかなり高いですし、ユーザも多いのでWeb上で色んな情報を得ることができます。こちらの記事では、brmsを用いて信号検出理論のパラメータを推定しています。これからも、どんどん新しい機能が追加されていくと思うので、要チェックですね。Rethinkingを解説しているドキュメントを見つけました。tidyverseざっと見たところ、brmsの具体的な使用例についても紹介しているようです。参考にしてみてください。