ggplot2でヒストグラムを描くときに外枠だけに線を引きたい
タイトルの通り、「ggplot2でヒストグラムを描くときに外枠だけ線を引きたい」という質問を受けました。図にすると以下のような状況です。
枠線を付けたいときにパッと思いつくのは、引数colourを指定することですが、そうすると下図・右のようになってしまいます。
library(ggplot2) bw <- diff(range(faithful$waiting))/20 # 左図 ------------------------------------------------- ggplot(data = faithful, mapping = aes(x = waiting)) + geom_histogram(binwidth = bw, fill = "skyblue") + labs(title = "通常のヒストグラム") # 右図 ------------------------------------------------- ggplot(data = faithful, mapping = aes(x = waiting)) + geom_histogram(binwidth = bw, fill = "skyblue", colour = "darkblue") + labs(title = "全ての棒に枠線を付けたヒストグラム")
で、どうやるんだろうと調べてみました。安直に「ggplot2 histogram only outline」というググり方をしましたが、幸いにもすぐにヒットしました。ggplot2の調べ物は英語でググるのが吉です。
参考：can I draw a histogram that has an outline around the entire shape, but not the individual bars?
ヒストグラムを描くときに、Y軸にマッピングする変数を省略するのではなく、mapping = aes(y = ..density..)としてkernel density estimateをマッピングすることを明示したうえで、geom_step()を用いてヒストグラムの輪郭だけを描くんですね。
あとはヒストグラムのレイヤを重ねてやれば、見た目上は外枠だけに線が描かれたように見えるはずです。
ただし注意が必要なのは、ヒストグラムはY軸が度数（count）である点です。そこで以下のページを参考に、ひと工夫します。ポイントは、geom_step(mapping = eval(bquote(aes(y=..count..)))というところです。
参考：データのhistogramと確率密度関数を同時に描きたい
勝った、第3部完！
dens <- density(faithful$waiting) bw <- diff(range(faithful$waiting))/20 ggplot(data = faithful, mapping = aes(x = waiting)) + geom_histogram(binwidth = bw, fill = "skyblue") + geom_step(mapping = eval(bquote(aes(y=..count..))), stat = "bin", binwidth = bw, colour = "darkblue") + xlim(range(dens$x))
ずれとるやんけ！
考えてみたらそりゃそうで、geom_step()が描く線の始点は、各階級幅の中央になるので、階級幅の半分だけずれることになります。仕方がないので、引数positionをいじって、X軸のずれを補正します。
dens <- density(faithful$waiting) bw <- diff(range(faithful$waiting))/20 ggplot(data = faithful, mapping = aes(x = waiting)) + geom_histogram(binwidth = bw, fill = "skyblue") + geom_step(mapping = eval(bquote(aes(y=..count..))), stat = "bin", binwidth = bw, position = position_nudge(x = -bw/2), #←階級幅の半分だけ左にずらしている colour = "darkblue") + xlim(range(dens$x))
完成。
ただWarning messages出てるんですよね...。うーん、もうちょっと改善の余地がありそうです。もしもっと適切な方法があれば教えていただけたら嬉しいです。
Warning messages: 1: Removed 2 rows containing missing values (geom_bar). 2: Removed 1 rows containing missing values (geom_path).
追記：
atusyさん（@Atushi776）から、「stat_bin()を使うとよい ＆ 引数pad = TRUEを書くとよい」とご提案いただきました。ありがとうございます！ 確かにこちらのほうがシンプルで良いですね。そしてWarning messagesも出ない！
bw <- diff(range(faithful$waiting))/20 ggplot(data = faithful, mapping = aes(x = waiting)) + geom_histogram(binwidth = bw, fill = "skyblue") + stat_bin(binwidth = bw, colour = "darkblue", geom = "step", position = position_nudge(x = -bw / 2), pad = TRUE)
Enjoy !
brmsパッケージを用いたベイズモデリング入門
はじめに
本記事は、Stan Advent Calendar 2018 16日目の記事です。brmsパッケージを用いたモデリングの方法や、便利機能、ちょっとしたtipsを紹介します。
brmsパッケージのコンセプトは、
ベイズで（Bayesian）、回帰モデルを（Regression Model）、扱おう、 Stanを使って（using Stan）
です。
読み方はそのまま「ビー・アール・エム・エス」です（僕は勝手にブラムスと呼んでいますが）。
その名の通り、Stanのラッパーパッケージですが、ユーザは自身でStanコードを書く必要はなく、モデルを指定すると自動的に内部でStanコードが生成され、サンプリングが実行されます。
最近は、brmsを用いて分析した論文も見かけるようになってきました。例えば難波修史さんのこの論文とか
- Dynamic Displays Enhance the Ability to Discriminate Genuine and Posed Facial Expressions of Emotion
僕もたまに使っているのですが、あまりにも多機能なので、しっかりと把握できていない部分も多いです。もし本記事に間違いがあればご指摘いただければ幸いです。
brmsパッケージの良いところ
その1
回帰モデルと一言でいっても、実はその守備範囲はかなり広いです。線形モデルも非線形モデルも、その気になればSEMもできるそうです（それならblavaanでいいと思いますが）。
- 一般線形モデル
- 一般化線形モデル
- 階層線形モデル
- その他の自由なモデリング
嬉しいのは、上記の1~3は、多くのRユーザが慣れ親しんでいるであろう記法と同じ記法でモデルを指定できることです。 例えば
- 回帰モデルでは
lm(y ~ x, data = dat) - ポアソン回帰モデルでは
glm(y ~ x, family = "poisson", data = dat) - 階層線形モデルではlme4パッケージやlmerTestパッケ―ジを用いて、
lme4::lmer(y ~ x + (1 | group), data = dat)lme4::glmer(y ~ x + (1 | group), family = "poisson", data = dat)
のように書きますよね。brmsパッケージでも、記法はほぼ同じになります。
その2
デフォルトで用意されている確率分布が、ものっっっっっっっすごく多いです。コンソールで
?brms::brmsfamily
と入力してみてください。用意されている確率分布の一覧が確認できます。 驚くのは、「ベルヌーイ分布とポアソン分布の混合分布である、ゼロ過剰ポアソン分布」等、異なる確率分布を組み合わせた、いわゆる混合分布も多く用意されていることです。通常、このような混合分布を用いたモデリングをStanで書くにはひと手間要るのですが、brmsパッケージでは確率分布を指定するだけで済みます。
また、様々な確率分布からの乱数を生成するための関数や...
確率（密度）を返す関数も用意されています。
その3
「内部で何をやっているのか」を知ることができます。
ベイズを用いる際に最も気になることの一つが、事前分布の設定でしょう。brmsでは、brms::get_prior()で指定したモデルの事前分布を知ることができます。brmsがデフォルトで採用している事前分布も分かるので、報告時に困ることはありません。
また、brms::make_stancode()で指定したモデルのStanコードを出力させることもできます。「この確率分布を使ったモデルをStanで書きたいとき、どうすればいいんだろう」と悩んだとき、brmsで同じ確率分布を使った適当なモデルでStanコードを出力させて、「なるほど、こう書けばいいのか」とヒントを得ることも可能です。
僕が驚いたのが、brms流のガンマ回帰の書き方です。parametersブロックと、modelブロックの一部を以下に抜粋します（事前分布のコードは省略しています）。
parameters { vector[Kc] b; // population-level effects real temp_Intercept; // temporary intercept real<lower=0> shape; // shape parameter } model { vector[N] mu = temp_Intercept + Xc * b; for (n in 1:N) { mu[n] = shape * exp(-(mu[n])); } target += gamma_lpdf(Y | shape, mu); }
modelブロックでは、まず切片と説明変数の線形結合をmuに代入しています。
ガンマ分布にはshapeとrate、2つのパラメータがあり、平均はshape/rateになります。brmsではガンマ回帰におけるデフォルトのリンク関数は対数リンク関数（つまり逆リンク関数は指数関数）なので、exp(mu) = shape/rateの関係があります。
これを変形すると、rate = shape * exp(-mu)になります。
ところがparametersブロックではrateに相当するパラメータはそもそも宣言しておらず、一度別の用途で使ったmuという変数を再利用して上書きしています。あたかも、「もともとのmuには興味ないでしょ」と言わんばかりに。
さらには、そのmuはmodelブロック内で宣言されているため、最終的にrate（に相当する）パラメータを取り出すことができません。なかなか思い切ったことをするなあと思ったのですが、こんな風にbrmsはサンプリング効率を高めるようなStanコードの書き方をしているため、推定がかなり速いです。
その4
他のパッケージとの連携が充実しています。
例えばbrmsはbayesplotパッケージを利用して、brms::pp_check()という関数により、容易に事後予測チェックを行うことができます（詳しくはこちら）。
その他にも、bridgesamplingパッケージを利用してbayes factorを容易に計算できたり、欠測値がある場合にmiceパッケージで多重代入法による補完を行ったあとで、補完したそれぞれのデータセットからサンプリングを行えたりします（詳しくはこちらから）。
brmsパッケージの使用例
それではさっそく、使用例を紹介していきましょう。
インストール & ロード
brmsはCRANに登録されているので、以下でOKです。なお本稿執筆時点における最新バージョンは2.6.0です。
install.packages("brms") library(brms)
上述のbayesplotパッケージやbridgesamplingパッケージ、それからStanをRから扱うためのrstanパッケージを含む、種々の依存パッケージが一緒にインストールされます。便利なものも多いので、初めてStanを使う人は、brmsをCRANからインストールするのが一番楽な気がします （Rtoolsは別途インストールする必要があります）。
2018/12/19追記：
記事を書いた直後にバージョン2.7.0にアップデートされました。"extending support for Gaussian processes"ですって...
一般線形モデル（回帰分析）
まずはベーシックな回帰分析をやってみましょう。サンプルデータは、プリセットのmtcarsです。32台の車の性能が収められたデータセットになります。
head(mtcars) mpg cyl disp hp drat wt qsec vs am gear carb Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4 Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4 Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1 Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1 Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2 Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1
モデルの指定
燃費（mpg）を、車体重量（wt）とマニュアル車かオートマ車か（am）で予測する重回帰分析を試みます。mtcars$amには0か1が格納されていますが、今回はこれらをfactor型に変換したうえで、交互作用項を含むモデルを指定することにします。
mtcars$am <- as.factor(mtcars$am)
stats::lm()を用いると以下のように書きますが...
fit <- lm(mpg ~ wt * am, data = mtcars) summary(fit)
brmsでは、どのようなモデルであっても、統一的にbrm()という関数のなかでモデルを指定します。もちろん、rstan::sampling()でサンプリングする際に指定するような諸々の引数も、brm()は対応しています。例えば、サンプリング回数（iter）、ウォームアップ期間（warmup）、乱数のシード（seed）、チェイン数（chain）など。
fit <- brm(mpg ~ wt * am, data = mtcars, iter = 2000, warmup = 1000, seed = 1234, chain = 4) summary(fit) Family: gaussian Links: mu = identity; sigma = identity Formula: mpg ~ wt * am Data: mtcars (Number of observations: 32) Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; total post-warmup samples = 4000 Population-Level Effects: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat Intercept 31.42 3.17 24.90 37.59 2139 1.00 wt -3.79 0.82 -5.37 -2.11 2150 1.00 am1 14.96 4.39 6.66 23.68 1682 1.00 wt:am1 -5.34 1.48 -8.37 -2.50 1828 1.00 Family Specific Parameters: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat sigma 2.72 0.41 2.09 3.61 1774 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).
また、サンプリング結果を格納したオブジェクトをsummary()に入れると要約された結果が返されますが、brms::waic()やbrms::loo()に入れると、WAICやlooを計算することもできます（内部ではlooパッケージを使用して計算しています）。
brms::waic(fit) Computed from 4000 by 32 log-likelihood matrix Estimate SE elpd_waic -78.9 3.9 p_waic 4.1 1.0 waic 157.8 7.9 Warning message: 4 (12.5%) p_waic estimates greater than 0.4. We recommend trying loo instead.
brms::loo(fit) Computed from 4000 by 32 log-likelihood matrix Estimate SE elpd_loo -79.0 4.0 p_loo 4.3 1.1 looic 158.1 7.9 ------ Monte Carlo SE of elpd_loo is 0.1. Pareto k diagnostic values: Count Pct. Min. n_eff (-Inf, 0.5] (good) 31 96.9% 1187 (0.5, 0.7] (ok) 1 3.1% 855 (0.7, 1] (bad) 0 0.0% <NA> (1, Inf) (very bad) 0 0.0% <NA> All Pareto k estimates are ok (k < 0.7). See help('pareto-k-diagnostic') for details.
事前分布の表示
brmsは、いくつかのパラメータについてデフォルトで事前分布を与えていることがあります。今回指定したモデルで、各パラメータにどのような事前分布が与えられていたのかを知るには、brms::get_prior()を使用します。この中にモデルを指定すると、以下のように一覧が表示されます。
brms::get_prior(mpg ~ wt * am, data = mtcars) prior class coef group resp dpar nlpar bound 1 b 2 b am 3 b wt 4 b wt:am 5 student_t(3, 19, 10) Intercept 6 student_t(3, 0, 10) sigma
Interceptの事前分布student_t(3, 19, 10)、19ってどういう根拠で設定されたんでしょうね...。brmsはこんな風に、いくつかのパラメータにデフォルトで弱情報の事前分布を与えていることがあるのですが、どういう根拠に基づいているのかは、まだよく分かっていません。
またbrmsでは基本的に（おそらく常に？）、回帰係数bには事前分布は明示されません。つまりデフォルトでは、無情報事前分布が適用されるということになります。これは後述する「bayes factorの計算」においては問題となるため、そのような場合にはユーザが事前分布を指定する必要がでてきます。
事前分布の指定
ユーザが任意の事前分布を指定したければ、brms::brm()のなかでprior =という引数にパラメータと事前分布をセットで指定します。指定方法には複数の記法があります。詳しくは、コンソール上で?brms::set_priorと入力して、ヘルプをご参照ください。
fit2 <- brm(mpg ~ wt * am, data = mtcars, iter = 2000, warmup = 1000, seed = 1234, chain = 4, prior = c(prior_string("normal(0, 10)", class = "b"), prior(student_t(3, 19, 10), class = Intercept), prior_(~student_t(3, 0, 10), class = ~sigma) ) )
Stanコードの表示
brms::make_stancode()のなかに、事前分布を含み任意のモデルを指定すると、そのStanコードをコンソールに表示してくれます。例えば先ほど事前分布を指定した、以下のモデルのStanコードを表示してみましょう。長いので、ここではmodelブロックだけ掲載します。
ちゃんと、ユーザが指定した事前分布も反映されていますね（target += normal_lpdf(b | 0, 10)という部分）。
brms::make_stancode(mpg ~ wt * am, data = mtcars, prior = c(prior_string("normal(0, 10)", class = "b"), prior(student_t(3, 19, 10), class = Intercept), prior_(~student_t(3, 0, 10), class = ~sigma) ) ) # 以下、出力 -------------------------------------------- model { vector[N] mu = temp_Intercept + Xc * b; // priors including all constants target += normal_lpdf(b | 0, 10); target += student_t_lpdf(temp_Intercept | 3, 19, 10); target += student_t_lpdf(sigma | 3, 0, 10) - 1 * student_t_lccdf(0 | 3, 0, 10); // likelihood including all constants if (!prior_only) { target += normal_lpdf(Y | mu, sigma); } }
そうそう、brmsのStanコードは全て、target +=記法を用いて書かれています。
結果の可視化
plot()
brmsは、サンプリング結果の視覚化にも長けています。rstan::stan_trace()やrstan::stan_dens()のように、描きたいグラフに応じて個別の関数が用意されていることが多いと思いますが、brmsはbayesplotパッケージのbayesplot::mcmc_combo()を用いて結果を可視化します。bayesplotパッケージについては、Stan Advent Calendar 2018の2日目の記事や、11日目の記事をご覧ください。
bayesplot::mcmc_combo()は、各パラメータについて、事後分布とトレースプロットを並べて表示してくれます。
plot(fit)
繰り返しますがこのグラフはbayesplotパッケージを用いて出力されているので、同パッケージの関数を用いて見栄え（Aesthetic）を調整することができます。また、bayesplot::mcmc_combo()と同様に、bayesplotパッケージのmcmc_ シリーズの関数のうち、どの種類の関数を組み合わせて表示するかも選択することができます（mcmc_シリーズについては、Stan Advent Calendar 2018 2日目の記事を参照してください）。ただしmcmc_combo()で表示できるグラフの種類は限られているようです。
例えばこんな風に、色調をピンクに変更し、トレースプロット（bayesplot::mcmc_trace()）とチェインごとの事後分布（bayesplot::mcmc_dens_overlay()）を並べることもできます。
bayesplot::color_scheme_set("pink") plot(fit, combo = c("trace", "dens_overlay"))
marginal_effects()
※注意：brms 2.11.0より、brms::marginal_effect()がbrms::conditional_effects()に名称変更されています。詳しくはこちら
さらにbrms::marginal_effects()を用いると、「主効果」や「交互作用」を可視化することもできます。今回は交互作用項を含むモデルを指定したので、いわゆる「単純傾斜」も自動的に可視化されます。下のグラフは、左から順に「wtの主効果」「amの主効果」「wtとamの交互作用」を示します。
引数effectsに何も指定しなければ、すべての項について可視化が行われますが、特定の項だけを指定することも可能です。また、交互作用項を指定する場合には、その順番を変えることで、どちらの変数についての単純効果を可視化するかも決めることができます。
marginal_effects(fit, effects = "wt:am") #左のグラフ marginal_effects(fit, effects = "am:wt") #右のグラフ
離散変数の場合は（今回はam）、各水準における他方の変数の単純効果が可視化されます。連続変数の場合は（今回はwt）、平均±1SDのポイントにおける、他方の変数の単純効果が示されます。
mean(mtcars$wt) + sd(mtcars$wt) #平均+1SD [1] 4.195707 mean(mtcars$wt) #平均 [1] 3.21725 mean(mtcars$wt) - sd(mtcars$wt) #平均-1SD [1] 2.238793
事後予測チェック
事後予測チェックについては、Stan Advent Calendar 2018 11日目の記事も参照してください。Stanでモデルを書いた場合、事後予測チェックを行うためには、generated quantitiesブロックでモデルから乱数を生成させる手続きが必要になります。その乱数を抽出して、実際に観測されたデータとの類似性を把握するわけですが、なんとbrmsでは、サンプリング結果が格納されたオブジェクトをbrms::pp_check()に入れるだけで、事後予測チェックができます！
brms::pp_check(fit)
Stan Advent Calendar 2018 11日目の記事を読んだ方はお察しの通り、これは内部でbayesplot::ppc_dens_overlay()を使用しています。
他の事後予測チェックにも対応していますよ。type =という引数に、描画形式を指定してみます。
brms::pp_check(fit, type = "error_hist")
これはbayesplot::ppc_error_hist()を使用していますね。
一般化線形モデル
次は、正規分布以外の指数分布族を扱えるように拡張した、一般化線形モデルを例に挙げます。今度は目的変数を、オートマ車かマニュアル車かamにし、説明変数に燃費mpgと重量wtを投入してみます。
stats::glm()と同様に、目的変数が従う確率分布を、引数family =に指定しましょう。今回は、オートマ車かマニュアル車かamは2値の変数なので、family = "bernoulli"とベルヌーイ分布を指定するのがよいように思われます（余談ですが、stats::glm()では2値の目的変数に対しては二項分布（family = "binomial"）を指定しますが、brmsでは2値の目的変数に対して二項分布を指定すると、Only 2 levels detected so that family 'bernoulli' might be a more efficient choice.とサジェストされます）。
ただしbrmsでは、family = "binomial"を指定した際の記法がstats::glm()と異なるので、ここでは紹介のためあえて二項分布を適用した書き方をしてみます。
stats::glm()ではこのように書きます。
#amをfactor型に変換していたのを、いったんリセット data(mtcars) fit <- glm(cbind(am, 1 - am) ~ mpg + wt, family = "binomial", data = mtcars) #目的変数が2値の場合は、 #fit <- glm(am ~ mpg + wt, family = "binomial", data = mtcars) でOK summary(fit)
cbind()のなかで何をしているのかというと、総試行数中に占める、成功試行数（am）と失敗試行数（1 - am）の列を指定しています。今回は1台の車がオートマ車かマニュアル車かを予測するので、試行数は1であるため、失敗試行の列は1 - amとしています。
さて、brmsでは、目的変数が二項分布に従うと仮定するとき、以下のように書きます。|の左側には成功試行数を、右側にはtrials()の中に試行数を入れます。
#amをfactor型に変換していたのを、いったんリセット data(mtcars) fit <- brm(am | trials(1) ~ mpg + wt, family = "binomial", data = mtcars) #目的変数が2値の場合は、 #fit <- brm(am ~ mpg + wt, family = "bernoulli", data = mtcars) でOK summary(fit)
brms::marginal_effects()で主効果を可視化してみると、ロジットリンク関数によって、直線ではなく曲線になっていることがわかりますね。
豊富なbrms family
ゼロ過剰ポアソン分布
さて、ここからがbrmsの真骨頂です。stats::glm()ではfamilyに指定できる確率分布は以下の通りです。
- gaussian
- binomial
- Gamma
- inverse.gaussian
- poisson
- quasi
が、brmsでは豊富なfamilyが用意されているということを、冒頭で述べました。試しにちょっと珍しい分布を使ってみましょう。
Mr.Unadon氏がTwitter上で「あなたは、片方しかない靴下を何足もっていますか？」とアンケートをとった結果を拝借します。ネタ元はこちらの片方しかない靴下の数をゼロ過剰ポアソン分布でモデリングしてみたです。
dat <- data.frame(Y = c(rep(0, 746), rep(1, 142), rep(2, 142), rep(3, 154), rep(4, 23), rep(5, 22), rep(6, 9), rep(7, 66) ) ) set.seed(1234) dat$X <- 50 + 2.5*dat$Y + rnorm(n = nrow(dat), mean = 0, sd = 10)
YがUnadon氏が集めた、片方しかない靴下の所持数のデータです。今回のサンプルの中では、ほとんどの人は靴下を片方だけなくすということはないようですね。僕が初めてこのアンケートを見た時、「靴下片方だけどうやって失くすねん」と思ったのですが、最近我が家で片方しかない靴下が発見されました。
さて、靴下の所持数はカウントデータなので、ポアソン分布と相性がよさそうに思えますが、上図のとおりゼロが圧倒的に多く、通常のポアソン分布の仮定を満たしそうにありません。そこで、ゼロが通常のポアソン分布よりも過剰であるという特徴を捉えた、ゼロ過剰ポアソン分布（zero-inflated poisson）という確率分布を考えます。
ゼロ過剰ポアソン分布は、ベルヌーイ分布とポアソン分布の混合分布です（以下の図はUnadon氏の記事から拝借しました）。靴下を失くすか否かがベルヌーイ分布に従うとと考え、もし失くさないなら当然、片方しかない靴下の所持数の回答は0。もし失くすなら、所持数はポアソン分布に従う、と考えます。なお、ポアソン分布に従う確率変数は0以上の整数をとりうることに注意してください。
一方Xは適当に作った変数ですが、仮に「日常の忙しさ指標（0-100点）」だとしておきましょう。仮定としては、忙しい人ほど靴下を片方なくしやすい、といった具合です。
さて、片方しかない靴下の所持数を、忙しさで予測する、シンプルなモデルを考えてみましょう。ここで、Stanでモデルを書くとしたら、Unadon氏の記事にあるように、target +=という記法を用いて、少し工夫した書き方が必要になります。
しかしbrmsでは、family = "zero_inflated_poisson"を指定するだけです。
fit <- brm(Y ~ X, data = dat, family = "zero_inflated_poisson", seed = 1234) summary(fit) Family: zero_inflated_poisson Links: mu = log; zi = identity Formula: Y ~ X Data: dat (Number of observations: 1304) Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; total post-warmup samples = 4000 Population-Level Effects: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat Intercept -1.64 0.16 -1.98 -1.33 2039 1.00 X 0.04 0.00 0.04 0.05 2283 1.00 Family Specific Parameters: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat zi 0.47 0.02 0.43 0.50 1868 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).
事後予測チェックしてみましょう。
brms::pp_check(fit)
いい感じですね。
階層線形モデル
本記事冒頭でも言及したように、lme4パッケージやlmerTestパッケ―ジを用いると、以下のように書きますよね。
lme4::lmer(y ~ x + (1 | group), data = dat)lme4::glmer(y ~ x + (1 | group), family = "poisson", data = dat)
brmsでは、関数をbrm()に変えるだけなので、本記事では説明を省略します。
モデル比較
brmsパッケージを用いてサンプリングした結果を利用して、モデル比較を行ってみます。モデル比較には様々な観点がありますが、ここではbayes factorを指標とすることにします。bayes factorについては、北條大樹さんのこちらの記事や、清水裕士さんのスライドをご参照ください。
brmsを用いてモデル比較するために、kidneyというデータセットを用いて、試してみましょう。survivalパッケージにも同じ名前のデータセットがあるので、区別するためにbrms::kidneyと指定したほうが無難です。
話を簡単にするために、ここでは目的変数に「病気の再発までの時間time」、説明変数に「患者の年齢age」を投入した、単回帰分析を考えます。また、本当は右側打ち切り（right-censored）が起こっているため、それを考慮したモデリングをするべきですが、打ち切りが起こっていないデータだけを抽出しています。
なおbrmsでは、brm(y | cens(censored) ~ x, data = dat)のように書くことで、打ち切りデータのモデリングも容易に実現できます（cens()のなかに、打ち切りの有無を判別する変数を投入します）。
#データセットの作成 --------------------------------------------- #1回目の再発までの時間 & 非打ち切りデータのみを使用 dat <- brms::kidney %>% dplyr::filter(recur == "1" & censored == 0) #データセットの行列数 ------------------------------------------- dim(dat) #38人分のデータが抽出された [1] 38 7 #冒頭6行分を表示 ------------------------------------------------ head(dat) time censored patient recur age sex disease 1 8 0 1 1 28 male other 2 23 0 2 1 48 female GN 3 22 0 3 1 32 male other 4 447 0 4 1 31 female other 5 30 0 5 1 10 male other 6 24 0 6 1 16 female other
目的変数timeの分布はこのように左に歪んでいます。
正規分布を仮定したモデル
まずはベーシックな単回帰分析を実行してみます。目的変数が正規分布に従うと想定することになるので、family = "normal"と指定します（familyを指定しなければ自動的に正規分布を仮定することになるので、省略も可能です）。
さらに以下3点を対処する必要があります。
- 全パラメータに、無情報でない事前分布を指定
iterを多めに- 正確には、有効なMCMCサンプルの数を多めに（2109/12/6追記）
save_all_pars = TRUEと引数を指定
モデルはこのようになります。
fit_n <- brm(time ~ age, family = "normal", prior = prior(normal(0, 5), class = b), seed = 1234, iter = 100000, warmup = 5000, data = dat, save_all_pars = TRUE)
事後予測チェックを行ってみましょう。
brms::pp_check(fit_n)
実際に観測されたデータとズレがありますね。上のヒストグラムを見ればわかる通り、今回はtimeが正規分布に従うと仮定するのは、適切ではないかもしれません。
最後にbrms::bridge_sampler()を用いて、bridge samplingによる対数周辺尤度を計算してみましょう。これは内部ではbridgesampling::bridge_sampler()を利用しています。
brms::bridge_sampler(fit_n) Iteration: 1 Iteration: 2 Iteration: 3 Iteration: 4 Bridge sampling estimate of the log marginal likelihood: -249.9636 Estimate obtained in 4 iteration(s) via method "normal".
ワイブル分布を仮定したモデル
brms::kidneyは生存時間解析に適したデータセットなので、比較対象として、パラメトリックな生存時間解析でよく用いられる、ワイブル分布を適用したモデルを考えてみます。
fit_w <- brm(time ~ age, family = weibull(), prior = prior(normal(0, 5), class = b), seed = 1234, iter = 100000, warmup = 5000, data = dat, save_all_pars = TRUE)
事後予測チェック。
brms::pp_check(fit_w)
先ほどよりマシになった気がしますね。
次はbridge samplingによる対数周辺尤度の計算。
brms::bridge_sampler(fit_w) Iteration: 1 Iteration: 2 Iteration: 3 Iteration: 4 Bridge sampling estimate of the log marginal likelihood: -232.0824 Estimate obtained in 4 iteration(s) via method "normal".
正規分布 vs ワイブル分布
いよいよモデル比較です。各モデルの対数周辺尤度に、いったん-1を掛けます。つまり、自由エネルギーとして計算します。
「自由エネルギーの差」の指数をとるとbayes factorになります。
fe_n = -1 * brms::bridge_sampler(fit_n)$logml #正規モデルの自由エネルギー fe_w = -1 * brms::bridge_sampler(fit_w)$logml #ワイブルモデルの自由エネルギー exp(fe_n - fe_w) #bayes factor [1] 58278089
今回のデータに対する当てはまりの観点では、正規分布を仮定したモデルに対して、ワイブル分布を仮定したモデルを支持する程度が、58278089倍...もう圧倒的に大きいと考えられます。
ただしbrmsは、bridgesampling::bayes_factor()を利用して直接bayes factorを計算する、brms::bayes_factor()という関数を用意しています。引数x1とx2に、各モデルのサンプリング結果を格納したオブジェクトを指定してください。
brms::bayes_factor(x1 = fit_w, x2 = fit_n) Iteration: 1 Iteration: 2 Iteration: 3 Iteration: 4 Iteration: 1 Iteration: 2 Iteration: 3 Iteration: 4 Estimated Bayes factor in favor of bridge1 over bridge2: 58273598.05176
やはり今回のデータに対する当てはまりの観点では、正規分布を仮定したモデルよりも、ワイブル分布を仮定したモデルのほうが良いと考えられます。
こんな風に、brmsを用いてモデル比較をすることもできます。ただしモデル比較という手続きは非常に難しいので、実際にはもっと慎重に行う必要があると思います。 以下、北條さんの記事から引用します。
ベイズファクターの注意点
- ベイズファクターは事前分布の影響をかなり大きくうけます ｡
→ 同一データ､同一モデルで事前分布のみを変化させた場合に､
どれだけベイズファクターが変化するかを検討するととてもよくわかります｡
- bridgesamplingパッケージでも計算がうまく行かないことはありますので､
全てのモデル間のベイズファクターが計算できるわけではありません｡
- モデル比較を行うのはとても難しいです｡
本来､我々が､比較すべきモデルは無限に存在します｡
- ベイズファクターは､2つのモデル間の比として考えるので､
それ以外のモデルについては考えてはいません｡
2019/12/6追記
清水裕士さんによる、Stan Advent Calendar 2019 6日目の記事：brmsパッケージで安易にベイズファクターを使うと死ぬ話で、brmsでベイズファクターを求める際の注意点が書かれています。
その他のtips
Stanコードの保存
brms::brm()でサンプリングを行う際に、save_model = "hogehoge"と引数を追加してみましょう（もちろんhogehogeのところは任意の名前で構いません）。
fit_w <- brm(time ~ age, data = dat, save_model = "hogehoge")
するとカレントディレクトリに、hogehogeという名前のファイルができていると思います。R上でこのファイルを開くと、brms::make_stancode()で生成できるようなStanコードが記録されていることが分かります。
仮にbrmsのバージョンが変わったとしても、少なくとも今回のサンプリングを実行した時点でのStanコードは判明しているので、結果の再生可能性の向上に寄与するかもしれません。
サンプリング結果の保存
brms::brm()でサンプリングを行う際に、file = "fugafuga"と引数を追加してみましょう（もちろんfugafugaのところは任意の名前で構いません）。
fit_w <- brm(time ~ age, data = dat, file = "fugafuga")
するとカレントディレクトリに、fugafuga.rdsという名前でサンプリング結果が保存されています。
以下のようにreadRDS()という関数で、あとでロードすることができるので、「しまった、うっかりRStudioを閉じてしまった！」なんて心配も不要です。また、例えばサンプリング結果をRMarkdownでまとめたい場合などにも重宝すると思います。
fit_again <- 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()です。bridgesampling::bf()はbayes factorを計算する関数で、bridgesampling::bayes_factor()の略称なのですが、brms::bf()はbayes factorとは関係がありません。brms::bf()は、brms::brmsformula()の略称で、brmsを用いて高度なモデリングを可能にするための関数です。
分位点回帰を例に説明しましょう。brmsは、分位点回帰も容易に実行することができます。その場合、asymmetric Laplace分布を適用して、以下のように書きます。brm()のなかに、さらにbf()という関数を書いて、その中でモデルを指定します。
fit <- brm(bf(y ~ x, quantile = 0.2), data = dat, family = asym_laplace() )
これ、bridgesampling::bf()のほうが読み込まれてしまうと、当然エラーになります。したがって以下のように、brms::bf()であることを明示したほうが無難です。
fit <- brm(brms::bf(y ~ x, quantile = 0.2), data = dat, family = asym_laplace() )
あるいは、パッケージの優先順位をつけてもいいかもしれないですね。
参考記事：R でパッケージの優先順位を変えたい #rstatsj
おわりに
ここで紹介した機能は、brmsのほんの一部にすぎません。モデリングの自由度もかなり高いですし、ユーザも多いのでWeb上で色んな情報を得ることができます。 例えばこちらの記事では、brmsを用いて信号検出理論のパラメータを推定しています。
きっとこれからも、どんどん新しい機能が追加されていくと思うので、要チェックですね。
そういえば、この記事を書くためにbrmsについて調べていたら、brmsとggplot2とtidyverseを用いてStatistical Rethinkingを解説しているドキュメントを見つけました。
Statistical Rethinking with brms, ggplot2, and the tidyverse
ざっと見たところ、brmsの具体的な使用例についても紹介しているようです。参考にしてみてください。
Enjoy!!
Introduction to bayesplot (ppc_ series)
はじめに
Stan Advent Calendar 2018 11日目の記事です。また、タイトルを見て察した人もいるかもしれませんが、Stan Advent Calendar 2018の2日目の記事である、北條大樹さんによるIntroduction to bayesplot (mcmc_ series) の続編でもあります。bayesplotパッケージそのものについては、本記事では説明を省略するため、まずは北條さんの記事をご一読ください。以下引用。
Stanには、既にggplot2ライクな事後処理用の関数が標準装備されています。ここら辺に関しては、以下の資料をご覧ください。
Stanの便利な事後処理関数
これとは別にStan公式が出しているbayesplotパッケージというものがあって、これも同じようなことができます。 既に日本語文献として小杉先生のサイトで紹介されていました（こちら）。上記も併せてご覧ください。この記事を書くうえでも参考にしています。
bayesplotパッケージでは、大きく二つの関数系に分かれています
1. mcmc_ シリーズ(収束や事後分布を確認するための関数群)
2. ppc_ シリーズ(事後予測チェックを行うのに便利な関数群)
この記事では、1.のmcmc_ シリーズの説明を行います。 2.については誰かもしくは未来の自分がやってくれることでしょうｗ
別に、2.を紹介してしまっても、構わんのだろう？
本記事で採用するモデル
僕はRにデフォルトで用意されているmtcarsというデータセットが好きなので、本記事でもこのデータセットを用います。mtcarsには、32種類の車のデータ（例：燃費mpg、車体重量wt、マニュアル車かオートマ車かam）が格納されています。
head(mtcars) mpg cyl disp hp drat wt qsec vs am gear carb Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4 Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4 Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1 Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1 Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2 Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1
ここでは、燃費を車体重量で予測する、単回帰分析をやってみます。Stanコードは以下の通りです（事前分布は省略しています）。あとで事後予測チェックを行うために、generated quantitiesブロックでy_predという変数を宣言し、今回のモデルに従う乱数を格納させています。
data{ int N; vector[N] Y; vector[N] X; } parameters{ real a; real b; real<lower = 0> sigma; } model{ Y ~ normal(a + b*X, sigma); } generated quantities{ vector[N] y_pred; for(n in 1:N){ y_pred[n] = normal_rng(a + b*X[n], sigma); } }
このモデルをregression.stanという名前で保存してコンパイルしたのち、R上でサンプリングを行います。あとでbayesplotパッケージを用いるので、ここでrstanパッケージと一緒にロードしておきます。
library(rstan) library(bayesplot) stanmodel = rstan::stan_model("regression.stan") standata = list(N = nrow(mtcars), Y = mtcars$mpg, X = mtcars$wt) fit = rstan::sampling(object = stanmodel, data = standata, seed = 1234, iter = 2000, warmup = 1000, chain = 4)
診断
北條さんの記事と同様に、bayesplotパッケージの関数を用いて、収束したかどうか確認してみましょう。切片a、傾きb、標準偏差sigmaの3つのパラメータについて、bayesplot::mcmc_combo()で事後分布とトレースプロットを並べて描画します。サンプリング結果を格納したオブジェクトを、as.array()で配列型にしてから与えることをお忘れなく。
bayesplot::mcmc_combo(as.array(fit), pars = c("a", "b", "sigma"))
便利ですね。
事後予測チェック
さて、ここからが、本記事の本命です。北條さんの記事では紹介されていなかった、ppc_ シリーズの関数を紹介していきます。bayesplotパッケージには、事後予測チェックを手軽に行える関数が、たくさん用意されています。
そもそも事後予測チェック（PPC: Posterior Predictive Check）とは、モデルチェックのための手続きです。
もし何かしらの仮説に基づきモデルを立て、そのモデルが現象をとらえられているとしたら、そのモデルから生成された将来のデータは、今回観測されたデータと似ていることが期待できます。
反対に、モデルから生成した将来のデータが、今回の観測データと似ていないとすると、モデルに改良の余地があるかもしれません。
ppc_ シリーズ
bayesplotパッケージでは、事後予測チェックのための関数は、いずれもppc_ から始まる名前が与えられています。その一覧は、available_ppc()という関数で知ることができます。あるいは公式サイトをご覧ください。
bayesplot::available_ppc() bayesplot PPC module: ppc_bars ppc_bars_grouped ppc_boxplot ppc_data ppc_dens ppc_dens_overlay ppc_ecdf_overlay ppc_error_binned ppc_error_hist ppc_error_hist_grouped ppc_error_scatter ppc_error_scatter_avg ppc_error_scatter_avg_vs_x ppc_freqpoly ppc_freqpoly_grouped ppc_hist ppc_intervals ppc_intervals_data ppc_intervals_grouped ppc_loo_intervals ppc_loo_pit ppc_loo_pit_overlay ppc_loo_pit_qq ppc_loo_ribbon ppc_ribbon ppc_ribbon_data ppc_ribbon_grouped ppc_rootogram ppc_scatter ppc_scatter_avg ppc_scatter_avg_grouped ppc_stat ppc_stat_2d ppc_stat_freqpoly_grouped ppc_stat_grouped ppc_violin_grouped
ppc_シリーズの関数はいずれも、yとyrepという引数を持ちます。
yには、サンプリングのためにStanに渡した目的変数を指定します。
y <- mtcars$mpg
一方yrepには、モデルから生成した乱数を与えます。行列形式でデータが並んでいればよいので、matrix型やdata.frame型でOKです。
行数は、事後予測チェックの描画のために用いる、サンプリングのステップ数です。
列数は、目的変数yと同じ長さ（すなわち、観測されたデータの個数）です。
今回はStanの中で、y_predというベクトルにモデルから生成した乱数を格納したので、以下のようにしてy_predだけを抽出します。この時点でR上ではmatrix型になっているので、このままyrepをppc_シリーズの関数に与えることができます。
yrep <- rstan::extract(fit)$y_pred
ただし、今回は、iter = 2000, warmup = 1000, chain = 4という設定でサンプリングを行ったため、モデルから生成した乱数は、y_pred[1]からy_pred[32]までそれぞれ(2000 - 1000) * 4 = 4000個あるのですが、すべてを描画に使うとごちゃごちゃして見にくくなってしまいますし、、処理もかなり重くなってしまいます。
View(yrep)
そのため一部の関数を用いて描画する際には、一度に使用する行数を絞る（つまり一部の乱数だけを用いる）ことが必要になってきます。
ppc_dens_overlay()
事後予測チェックを行うときに、個人的に便利だなと思うのが、この関数です。ggplot2パッケージのgeom_density()で描けるような分布を、yとyrepそれぞれについて描き、重ね合わせます。10ステップだけ絞り込んで描画を行っていることに注目してください。
bayesplot::ppc_dens_overlay(y = y, yrep = yrep[1:10, ])
色の濃い線で描かれたのが実際に観測されたデータ、薄い線で描かれているのがモデルから生成された乱数です。ピタッと重なっているわけではありませんが、悪くはなさそうです。
ppc_ribbon()
今回の回帰分析のように、目的変数Yに影響を与える説明変数Xが存在する（とモデルを組んでいる）とき、事後予測チェックにおいてはこれらの関係も重要になってくるでしょう。説明変数の変動に対する目的変数の変動の関係は、モデルから生成した乱数に対しても同様に認められるかどうかをチェックしてみましょう。
bayesplot::ppc_ribbon()は、X軸の値に対応するY軸の値が複数存在するとき、それらの区間を描画します。ggplot2::geom_ribbon()と似たようなものですね。今回はモデルから生成した乱数の一部ではなく、すべてを用いて描画してみます。
bayesplot::ppc_ribbon(y = mtcars$mpg, yrep = rstan::extract(fit)$y_pred, x = mtcars$wt, prob = 0.5, prob_outer = 0.95 )
引数probには内側の色の濃い帯、引数prob_outerには外側の薄い帯の描画に使う、パーセンタイルを指定します。最も色の濃い折れ線が、実際に観測された目的変数です。
便利ですね...そういや僕、去年のStan Advent Calendar 2017では、ggfanパッケージを用いて、同様の事後予測チェックを行う記事を書いたんですが、もうbayesplotでいいよね...。
ppc_intervals()
今回使用したmtcarsデータセットは、サンプルサイズが32しかないので、上記のように帯を描くだけでなく、1点1点に対して区間を示すことも可能です。引数はそのまま、関数をppc_intervals()に変えてみましょう。
bayesplot::ppc_intervals(y = mtcars$mpg, yrep = rstan::extract(fit)$y_pred, x = mtcars$wt, prob = 0.5, prob_outer = 0.95 )
その他にも、非常に多くの事後予測チェック用の関数が用意されています。詳しくは公式サイトをご覧ください。
グラフのカスタマイズ
bayesplotパッケージの関数の内部ではggplotオブジェクトを作っているので、以下のようにggplot2の関数を用いてレイヤを追加して、グラフを加工することもできます。以下の例では、任意の切片と傾きを持つ直線を追加するggplot2の幾何学的オブジェクトである、geom_abline()を使用しています。切片と傾きには、それぞれ推定したパラメータの事後平均を用いています。
bayesplot::ppc_intervals(y = mtcars$mpg, yrep = rstan::extract(fit)$y_pred, x = mtcars$wt, prob = 0.5, prob_outer = 0.95) + geom_abline(intercept = mean(rstan::extract(fit)$a), slope = mean(rstan::extract(fit)$b), linetype = "dashed")
もちろんppc_シリーズの関数だけでなく、mcmc_シリーズの関数でも同様です。以下の例では、themeをbwに変更するとともに、X軸上の任意の位置に垂直線を描くggplot2の幾何学的オブジェクトgeom_vline()を用いて、x = 0の位置に垂直線を描画しています。
bayesplot::mcmc_areas_ridges(as.array(fit), pars = c("a", "b", "sigma"), prob = 0.95, prob_outer = 1) + theme_bw(base_size = 15) + geom_vline(xintercept = 0, linetype = "dashed")
あるいはbayesplotパッケージが用意している、グラフの見栄え（Aesthetics）を調整する関数を用いることも可能です。一覧は公式サイトの「Aesthetics」のところをご覧ください。
例えば直前のグラフは、bayesplotパッケージが用意しているvline_0()という関数を用いて実現することもできます。ついでに、color_scheme_set()という関数で、色調を紫色ベースに変えてみましょう。
bayesplot::color_scheme_set("purple") bayesplot::mcmc_areas_ridges(as.array(fit), pars = c("a", "b", "sigma"), prob = 0.95, prob_outer = 1) + theme_bw(base_size = 15) + vline_0(linetype = "dashed")
現在セットされているカラースキームは、color_scheme_view()でViewerペインに表示させることができます。ちなみにデフォルトのカラースキームは、"blue"です。
bayesplot::color_scheme_view()
また、そのカラーコードはcolor_scheme_get()でコンソールに表示されます。
bayesplot::color_scheme_get() purple 1 #e5cce5 2 #bf7fbf 3 #a64ca6 4 #800080 5 #660066 6 #400040
こんな風に、カスタマイズの余地も残されています。
ただ自分自身では、カラースキームの変更はbayesplot::color_scheme_set()を使いますが、その他のカスタマイズは普通にggplot2の関数を使っています（要はラッパー関数なので）。
他のパッケージとの連携
実はbayesplotパッケージは、brmsパッケージやrstanarmパッケージと非常に相性が良いです（brmsは、プロットにbayesplotパッケージを利用していますし）。これらは、回帰分析や一般化線形モデル、階層線形モデルをRで実行する際の記法と類似した記法で、様々なモデルのパラメータをベイズ推定することを可能にするパッケージです。特にbrmsパッケージはその守備範囲がかなり広く、拡張性も高いため、非常に便利なパッケージです。
これらのパッケージを用いてベイズ推定した際の事後予測チェックについては、Stan Advent Calendar 2018 16日目に投稿予定の、brmsパッケージ入門編の記事で触れる予定です。
おわりに
事後予測チェックは非常に重要な手続きだと思うのですが、それを行うためにggplot2による可視化や、dplyr・tidyrによるデータ整形の技術が必要になってくるとしたら、億劫になって省略してしまうこともあるかもしれません。
bayesplotパッケージは、多くの関数を用意することで可視化の手間を省いてくれる（つまりやってみようと思わせてくれる）、非常に有用なパッケージだと思います。
今回は一部の関数しか紹介できませんでしたが、ぜひ他にも色々と試してみてください。
Enjoy!!
xaringanthemerパッケージのライトな解説
たまたま最近NARUTOを読み返す機会があって、その流れでxaringanパッケージを使い始めてみたら気に入ったので、これからもたまに使っていこうかと思っています。
一応説明しておくと...xaringanとは、Yihui Xie氏により作られた、RMarkdownでHTMLスライドを生成するためのパッケージです。
詳しくは、氏の著書「R Markdown: The Definitive Guide」をご参照ください。
RMarkdownでHTMLスライドを生成するには、デフォルトのioslidesやSlidyを使うという手もあるのですが、拡張性の高さやちょっとした小技が使えるのが、xaringanの魅力かと思います。
例えばこちらのスライドはxaringanで作成されていますが、このスライドを開きながら、キーボードの「m」を押すと、ちょっとした忍術が使えます。
さて、xaringanパッケージによりHTMLスライドが生成できるわけなので、見栄えを色々と変更するには、CSSをいじればいいということになります。xaringanにはデフォルトのテーマ（theme）が設定されていますが、このテーマを変更する、いくつかのテンプレートが用意されています。例えばmetropolisというテーマを選択すると、こちらのページのようにスライドの見栄えが変わります。詳しくは「R Markdown: The Definitive Guide」をご参照ください。
ところが！ xaringanthemerというパッケージを見つけてしまいました。その名の通り、xaringanで作成するHTMLスライドのテーマ（theme）を、ユーザが任意に操作しやすくするためのパッケージです。CSSよくわかんないや...という人でも、手軽にスライドの見栄えを調整できるようになります。
あまりの便利さに驚いたんですが、僕が調べた限りでは、xaringanthemerパッケージに関する日本語の解説記事は見つからなかったので、備忘録がてら記事にまとめました。
と言っても使い方は簡単で、とりあえずパッケージをインストールしたあとで...
# install.packages("devtools") devtools::install_github("gadenbuie/xaringanthemer")
新規.Rmdファイルを立ち上げる際に、Ninja Themed Presentationを選択するだけです。
詳しくはREADMEのところを読めばわかりますが、xaringanthemerも様々なテーマのテンプレートを用意しています。そしてそのテンプレート（これ自体が関数）は、様々な「見栄え」に関する引数を持っています。文字の色からフォントから見出しのサイズから...引数の多さにびっくりします。
例えばこちらの記事に埋め込んだスライドは、mono_light()というテンプレートを使用しているため、緑を基調とした色合いになっていますが、下のように任意に設定を変更していました。必要なことは、.Rmdファイル冒頭で適当なRチャンクを作り、その中に以下を記述するだけです。
library(xaringanthemer) mono_light( base_color = "#1c5253", text_font_size = "30px", code_font_size = "20px", padding = "1em 1em 1em 1em",#margin header_h1_font_size = "45px", header_h2_font_size = "40px", header_h3_font_size = "35px", header_color = "midnightblue", header_font_google = google_font("Song Myung"), text_font_google = google_font("Song Myung", "400", "400i"), code_font_google = google_font("Song Myung"), link_color = "chocolate" )
また、任意のCSSをextra_cssという引数に与え、カスタマイズすることもできます。例えば以下の例では、指定した文字の色を赤にしたり、フォントサイズをデフォルトの50%にしたりするための準備をしています。
extra_css <- list( ".red" = list(color = "red"), ".small" = list("font-size" = "50%"), ".full-width" = list( display = "flex", width = "100%", flex = "1 1 auto" ) ) mono_light( extra_css = extra_css )
あとはMarkdown記法でドキュメントを執筆する際に、
買い物しようと街まで出かけたら、.red[財布忘れちゃったよ！] あーあ、ツイてない。.small[トホホ...]
などと書いてレンダリングすれば、「財布忘れちゃったよ」の部分が赤く、「トホホ...」の部分のフォントサイズが半分になって表示されます。
僕自身、xaringanを使い始めたばかりで、まだよくわかっていない部分も多いのですが、お役に立てば幸いです。
Enjoy!!
「R Markdownによるドキュメント生成と バージョン管理入門」という発表をしてきました
2018年8月23日に関西学院大学のCAPS（応用心理科学研究センター）で行われたセミナーで、R Markdownによるドキュメント生成と バージョン管理入門という発表をさせていただきました。
※当日の配布資料から削った＆修正したスライドがあります
主に心理学を専門とする、研究者や大学院生を対象とした発表でした。まだR Markdownを使ったことがない方々を想定して、以下のような話をしました。
また、ドキュメントを生成して終わりではなく、Gitによるバージョン管理も重要という話題提供も行いました。
これは本筋ではなく（ご依頼いただいた内容ではなく）、勝手に盛り込んでしまったものでしたが、なかには非常に関心を持っていただいた方もおられたようで、よかったです。
今回は入門編ということで、GitHub Desktopを用いた管理方法について紹介しました。
個人的な話になりますが、今回このセミナーのために初めてxaringanパッケージでスライドを作ってみました。まだ慣れていない部分もありますが、拡張性の高さにとても魅力を感じました。
自分自身の勉強にもなり、とてもいい機会になりました。関係者の皆様、ありがとうございました。
RユーザのためのRStudio[実践]入門という本が出ます。
はじめに
RユーザのためのRStudio[実践]入門−tidyverseによるモダンな分析フローの世界−という書籍の執筆に関わらせていただく機会をいただきました。
「データ分析ワークフローを一通り解説していきます。データの収集(2章)、データの整形(3章)、可視化(4章)、レポーティング(5章)など、データ分析に欠かせないこれらの要素の基礎を押さえることができます。」という触れ込みですが、私は4章：ggplot2によるデータ可視化というパートを担当しています。
書籍の概要の紹介
共著者の皆さんもすでに各ブログで紹介してくださっているので、重複した情報になりますが、構成は以下の通りです。
- RStudioの基礎 (担当: @y__mattuさん)
- RStudio（version 1.1.423）のインストールや基本的操作、データ読み込み方法など
- スクレイピングによるデータ収集 (担当: @y__mattuさん)
- Web上の様々なデータの取得方法
- dplyr/tidyrによるデータ前処理 (担当: @yutannihilation さん)
- 後の分析や可視化の効率を高めるための、tidy（整然）なデータに整形する方法
- ggplot2によるデータ可視化 (担当: @kyn02666)
- 可視化の文法の解説や、目的に応じたレイアウトの調整方法
- R Markdownによるレポート生成 (担当: @kazutan さん)
- Rによる処理をシームレスにドキュメント・レポートへ出力する方法
本書のアピールポイント
tidyverseへの準拠
著者らは本書を「宇宙本」と呼んでいますが、それは本書がtidyverseパッケージを利用した分析フローを解説しているためです。本書カバーの宇宙や惑星のイラストは、tidyverseという宇宙に広がる、様々なパッケージ群を表しています。1～2章担当の@y__mattuさんの言葉を借りると、本書には以下のような特徴があるかと思います。
- tidyverseを扱った書籍が少なく、特に入門書はほぼ皆無である
- 本書はtidyverseの入門者でも読めるように作られている
- データ取得、読み込み、前処理、可視化、レポーティングといったどんな分析をするにも必ず通るであろうプロセスを解説した本
tidyverseパッケージ群は、いまだ改良・開発が進んでいるので、賞味期限の早い内容になってしまうのではないかという懸念があるかもしれませんが、5章担当の@kazutanさんが「これは本書に限らず技術書全般で言える問題です。そこで本書は陳腐化しないようコアな機能を中心に、さらに知らないと損をするような機能を重点的に紹介しています。」とおっしゃるように、根っこの部分をしっかり紹介することを、著者一同、念頭においていました。
贅沢な内容
総ページ数が240ページ、ということは各章は40～50ページ程度ということになります。したがって、tidyverseに含まれるパッケージの中にも言及していないものもありますし（例えばpurrrパッケージやbroomパッケージ）、決してtidyverseを隅々まで徹底的に解説しているというわけではありません。それは、繰り返しになりますが本書が主眼に置いているのは、tidyverseパッケージに準拠した分析フローであり、tidyverseパッケージ群の解説そのものが目的ではないからです。本書の内容については、目次をご参照ください。
それでも、各章は「おお！」と思わず息の漏れるような内容になっていると思います（というか私は、他の章を読んで、何度もそうなりました）。 1～2章担当の@y__mattuさんの言葉をお借りすると、それは「『この内容ならこの人しかいないよね』という人物が書いています。R勉強会やWeb上で多くの記事を発信し、時にはパッケージの開発動向まで追ってしまうような頼もしい方々です。」からだと思います。
あ、いえ、私は違うんですけどね。「奴は四天王の中でも最弱」っていうか。いやそもそも四天王ですらなくて、ただのggplot2が好きな人っていうか。
とはいえ、日本社会心理学会 第5回春の方法論セミナー（R/RStudio入門 データ可視化）など、ggplot2によるデータ可視化について、聴衆の反応をリアルタイムに感じながら解説を行ってきた経験はそこそこあると思います。今回の私の担当章では、それらの経験を踏まえて、どのような解説が分かりやすいのか、どのような喩えが分かりやすいのか、どの点を重点的に解説したほうがいいのか、など、色々と工夫したつもりです。特に、図4.1はかなり頑張って作りました。
では私が担当させていただいた4章では、いったいどのような内容を紹介したのか、以下に簡潔にまとめたいと思います。
4章：ggplot2によるデータ可視化
紹介したこと
意識したことは、
- ggplot2の背後にある「可視化の文法」を紹介すること
- 可視化から図の共有までの手順を時系列的に解説すること
- 用例集を兼ねること
ということでした。可視化の文法が存在することの利点は、それを理解することで、様々なグラフを似たような記法で作成できるという、効率の良さです。したがって、その文法を伝えることを重視しました。また、種々の機能や文法を散発的に解説するのではなく、料理レシピのように手順が明確になるように意識しました。
それと同時に、各コードが用例集としても利用可能なように、事例を工夫しました（もちろん、分野によってよく用いられるグラフは異なると思われるので、一般化はできないかもしれませんが）。
紹介していないこと
第一に、かなり重要なことなのですが、最新バージョンのggplot2 version 2.3.0には準拠していません。version 2.3.0における変更点は、こちらの@yutannihilationさんの解説記事をご参照ください。本書の発売が2018年6月29日、ggplot2 version 2.3.0のリリースが6月25日なので、仕方なかったんです...。というわけで、本章では、本書の執筆時点における最新バージョンである、ggplot2 version 2.2.1に基づき解説しました。もっとも、これらの変更による影響が出ないように配慮して本章を執筆しましたので、ご安心ください。
第二に、本書のメインタイトルはRユーザのためのRStudio[実践]入門という名前ですが、4章に関しては、あまりRStudioならではの機能は解説していません。具体的には、Addinです。
ggplot2では非常に手軽に・効率的にデータを可視化することが出来ますが、「見栄えを細かく調整する」ことは結構大変です。例えば以下のようなグラフを描いたとして...
g <- ggplot(data = iris, mapping = aes(x = Species, y = Sepal.Width, fill = Species)) +
geom_boxplot()
print(g)
凡例の位置を上に、各水準を水平方向に並べようとしたら、以下のようにtheme()を追記する必要があります。
g <- ggplot(data = iris, mapping = aes(x = Species, y = Sepal.Width, fill = Species)) +
geom_boxplot() +
theme(legend.position = "top", legend.direction = "horizontal")
print(g)
ただ、このような「見栄えの調整」のためのコードは、決して頻繁に書くというわけではないと思うので、思い出すのが大変です。いえ、私の話なんですけど。
このような場合に、見栄えの調整をGUIで操作できる、ggThemeAssistというAddinが重宝します。使用方法はこちらをご参照ください。他にも、カラーピッカーのAddinであるcolourpickerなど（使用方法はこちらをご参照ください）、便利なAddinがたくさんあります。
Addinを紹介するかどうか悩んだのですが、あくまで分析フローの紹介が本書の主眼であることから、割愛することにしました。
おわりに
長くなりましたが、まだまだ語り切れない制作秘話がたくさんあります。著者一同、熱い思いを込めて執筆した本なので、ぜひ本書を手に取っていただければ幸いです。そして、勉強会などでお目にかかった際には、ぜひご感想をお聞かせいただければ幸いです。
私自身、本書を執筆する中で新たに発見したことも多々あり、このような機会をいただけたことに、心から感謝しています。
そして何より、たびたびの校正にも丁寧に応じてくださり、執筆の舵を取ってくださった編集者T氏や、共同執筆者の皆様に、厚く御礼申し上げます。
Enjoy!!
