Extracting and visualizing tidy draws from rstanarm models

イントロダクション

この vignette は、tidybayes と ggdist パッケージを使用して、 tidy データフレームを抽出し、可視化する方法について説明する。モデル変数の事後分布からのドロー、平均、および rstanarm からの予測値である。より一般的な tidybayes の紹介と、汎用のベイズモデリング言語 (Stan や JAGS など) での使用については、vignette("tidybayes") を参照。

セットアップ

このvignetteを実行するには、以下のライブラリが必要である。

library(magrittr)
library(dplyr)
library(purrr)
library(forcats)
library(tidyr)
library(modelr)
library(ggdist)
library(tidybayes)
library(ggplot2)
library(cowplot)
library(rstan)
library(rstanarm)
library(RColorBrewer)

theme_set(theme_tidybayes() + panel_border())

これらのオプションは、Stanの動作を高速化するためのものである。

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

サンプルデータセット

tidybayes を実証するために、5つの条件からそれぞれ10個のオブザベーションを持つ単純なデータセットを使用することにする。

set.seed(5)
n = 10
n_condition = 5
ABC =
  tibble(
    condition = rep(c("A","B","C","D","E"), n),
    response = rnorm(n * 5, c(0,1,2,1,-1), 0.5)
  )

データのスナップショットはこのような感じである。

head(ABC, 10)

これは典型的な整頓されたフォーマットのデータフレーム: 1行に1つの観測値である。グラフィカルに

ABC %>%
  ggplot(aes(y = condition, x = response)) +
  geom_point()

#モデル {#model}

大域平均に向かって収縮する階層的モデルを当てはめよう。

m = stan_lmer(response ~ (1|condition), data = ABC, 
  prior = normal(0, 1, autoscale = FALSE),
  prior_aux = student_t(3, 0, 1, autoscale = FALSE),
  adapt_delta = .99)

結果はこのようになる。

m

## stan_lmer
##  family:       gaussian [identity]
##  formula:      response ~ (1 | condition)
##  observations: 50
## ------
##             Median MAD_SD
## (Intercept) 0.6    0.5   
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 0.6    0.1   
## 
## Error terms:
##  Groups    Name        Std.Dev.
##  condition (Intercept) 1.13    
##  Residual              0.56    
## Num. levels: condition 5 
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

spread_draws を使ってフィットからドローをtidyフォーマットで抽出する

さて、結果が出たので、楽しいことが始まります：整頓されたフォーマットで描画を取り出すことである！まず、 関数を使って、生のモデル変数名のリストを取得する。まず、get_variables() 関数を使って、生のモデル変数名のリストを取得し、モデルからどの変数を抽出できるかを知ることにする。

get_variables(m)

##  [1] "(Intercept)"                              "b[(Intercept) condition:A]"              
##  [3] "b[(Intercept) condition:B]"               "b[(Intercept) condition:C]"              
##  [5] "b[(Intercept) condition:D]"               "b[(Intercept) condition:E]"              
##  [7] "sigma"                                    "Sigma[condition:(Intercept),(Intercept)]"
##  [9] "accept_stat__"                            "stepsize__"                              
## [11] "treedepth__"                              "n_leapfrog__"                            
## [13] "divergent__"                              "energy__"

ここで、(Intercept) は大域的な平均値、b のパラメータは各条件におけるその平均値からのオフセットである。これらのパラメータが与えられると

• b [(Intercept) condition:A]
• b [(Intercept) condition:B]
• b [(Intercept) condition:C]
• b [(Intercept) condition:D]
• b [(Intercept) condition:E]

各行が b [(Intercept) condition:A] , b [(Intercept) condition:B] , ...:C] , ...:D] , or ...:E] の どちらか一方からの draw であるデータフレームが必要だろう。 , そして、行がどの chain/iteration/draw から来たか、そしてどの条件 ( A to E ) のものであるかのインデックスを持つ列があるデータフレームが欲しいと思うだろう。これにより、条件ごとにグループ化された量を簡単に計算したり、ggplotを使用して条件ごとにプロットを生成したり、あるいはドローとオリジナルデータをマージしてデータと後置をプロットしたりすることができるようになる。

tidybayes の主力は spread_draws() 関数で、この関数はこの抽出を行ってくれる。この関数には、モデル変数とそのインデックスを整頓されたデータフレームに抽出するために使用できる簡単な仕様書式が含まれている。

m %>%
  spread_draws(b[term,group]) %>%
  head(10)

## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## ℹ Please use `gather()` instead.
## ℹ The deprecated feature was likely used in the tidybayes package.
##   Please report the issue at <]8;;https://github.com/mjskay/tidybayes/issues/newhttps://github.com/mjskay/tidybayes/issues/new]8;;>.

m %>%
  spread_draws(b[t,g]) %>%
  head(10)

m %>%
  spread_draws(b[,group]) %>%
  head(10)

m %>% 
  spread_draws(b[,group]) %>%
  separate(group, c("group", "condition"), ":") %>%
  head(10)

m %>% 
  spread_draws(b[,group,condition], sep = "[, :]") %>%
  head(10)

ポイント要約と間隔

m %>%
  spread_draws(`(Intercept)`, sigma) %>%
  head(10)

m %>%
  spread_draws(`(Intercept)`, sigma) %>%
  median_qi(`(Intercept)`, sigma)

m %>%
  spread_draws(`(Intercept)`, sigma) %>%
  median_qi()

m %>%
  gather_draws(`(Intercept)`, sigma) %>%
  median_qi()

m %>%
  spread_draws(b[,group]) %>%
  median_qi()

m %>%
  spread_draws(b[,group]) %>%
  group_by(group) %>%       # this line not necessary (done by spread_draws)
  median_qi(b)                # b is not necessary (it is the only non-group column)

m %>%
  spread_draws(b[,group]) %>%
  summarise_draws()

異なるインデックスを持つ変数を1つの整頓されたフォーマットのデータフレームに結合する

spread_draws() と gather_draws() は、異なるインデックスを持つ変数を同じデータフレームに抽出することをサポートしている。同じ名前のインデックスは自動的にマッチングされ、必要に応じて値が複製され、すべてのインデックスのレベルの組み合わせごとに1つの行が作成される。例えば、各条件における平均を計算したい場合がある (これを condition_mean と呼ぶ) 。このモデルでは、その平均は切片((Intercept) )と与えられた条件での効果(b)である。

(Intercept) と b からの draw を一つのデータフレームにまとめることができる。

m %>% 
  spread_draws(`(Intercept)`, b[,group]) %>%
  head(10)

各 draw の中で、(Intercept) は b の各インデックスに対応するように必要に応じて繰り返される。 したがって、dplyr の mutate 関数を使用して、それらの合計 condition_mean (各条件の平均) を求めることができる。

m %>%
  spread_draws(`(Intercept)`, b[,group]) %>%
  mutate(condition_mean = `(Intercept)` + b) %>%
  median_qi(condition_mean)

median_qi() は整頓された評価 ( vignette("tidy-evaluation ", package =" rlang") を使うので、列名だけでなく、列の式も受け取ることができる。したがって、condition_mean の計算を mutate から median_qi() に移すことで、上の例を単純化することができる。

m %>%
  spread_draws(`(Intercept)`, b[,group]) %>%
  median_qi(condition_mean = `(Intercept)` + b)

複数の確率レベルを持つ区間をプロットする

median_qi() とその姉妹関数は、 引数を設定することで、任意の数の確率区間を生成することができる。 .width =

m %>%
  spread_draws(`(Intercept)`, b[,group]) %>%
  median_qi(condition_mean = `(Intercept)` + b, .width = c(.95, .8, .5))

結果は整然とした形式で、1グループにつき1行、不確定性区間幅( .width )である。これはプロットを容易にする。たとえば、-.width を size の美学に割り当てると、すべての区間が表示され、太い線がより小さい区間に対応するようになる。ggdist::geom_pointinterval() geomは、複数の確率レベルを持つ点のプロットを作成するために、データ中の .width 列に基づいて、size 美学を自動的に適切に設定する。

m %>%
  spread_draws(`(Intercept)`, b[,group]) %>%
  median_qi(condition_mean = `(Intercept)` + b, .width = c(.95, .66)) %>%
  ggplot(aes(y = group, x = condition_mean, xmin = .lower, xmax = .upper)) +
  geom_pointinterval()

## Warning: Using the `size` aesthietic with geom_segment was deprecated in ggplot2 3.4.0.
## ℹ Please use the `linewidth` aesthetic instead.

密度を持つ区間

区間とともに密度を見るには、ggdist::stat_eye() (区間とバイオリンプロットを組み合わせた「アイプロット」) 、または ggdist::stat_halfeye() (区間＋密度プロット) を使用することができる。

m %>%
  spread_draws(`(Intercept)`, b[,group]) %>%
  mutate(condition_mean = `(Intercept)` + b) %>%
  ggplot(aes(y = group, x = condition_mean)) +
  stat_halfeye()

あるいは、密度の一部をカラーでアノテートしたいとする。fill の美学は、ggdist::stat_halfeye() を含む ggdist::geom_slabinterval() ファミリーのすべてのジオムおよび統計において、スラブ内で変化させることができる。例えば、ドメイン固有のROPE (region of practical equivalence) をアノテートしたい場合、次のようなことが可能である。

m %>%
  spread_draws(`(Intercept)`, b[,group]) %>%
  mutate(condition_mean = `(Intercept)` + b) %>%
  ggplot(aes(y = group, x = condition_mean, fill = stat(abs(x) < .8))) +
  stat_halfeye() +
  geom_vline(xintercept = c(-.8, .8), linetype = "dashed") +
  scale_fill_manual(values = c("gray80", "skyblue"))

## Warning: `stat(abs(x) < 0.8)` was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(abs(x) < 0.8)` instead.

その他、分布の可視化を行う stat_slabinterval

分布を視覚化するための様々な追加統計が、ggdist::geom_slabinterval() ファミリーの統計とジオムにはある。

を参照。 vignette("slabinterval ", package =" ggdist") を参照。

事後平均値と予測値

前の例のように条件付き平均を手動で計算するのではなく、rstanarm::posterior_epred() (事後予測値の期待値から事後ドローを与える；すなわち、条件付き平均の事後分布) に類似しているが、整然としたデータ形式を使用する add_epred_draws() を使用することができる。これを modelr::data_grid() と組み合わせて、まず欲しい予測を記述したグリッドを生成し、次にそのグリッドを条件付き平均からのドローの長大なデータフレームに変換することができる。

ABC %>%
  data_grid(condition) %>%
  add_epred_draws(m) %>%
  head(10)

この例をプロットするために、ggplot 内で描画を点要約と区間に要約する ggdist::geom_pointinterval() の代わりに ggdist::stat_pointinterval() を使用することも紹介する。

ABC %>%
  data_grid(condition) %>%
  add_epred_draws(m) %>%
  ggplot(aes(x = .epred, y = condition)) +
  stat_pointinterval(.width = c(.66, .95))

分位点プロット

アルファレベルがたまたまあなたが行おうとしている決定と一致する場合、区間は良いが、事後的な形状を得ることはより良いことである (したがって、上記の目のプロット) 。一方、密度プロットから推論するのは不正確である (ある形状の面積を別の形状の割合で推定するのは難しい知覚的作業である) 。頻度形式での確率の推論はより簡単で、 quantile dotplots ( Kay et al. 2016, Fernandes et al. 2018)の動機となった。これは任意の間隔 (プロットのドット解像度まで、下の例では100) を正確に推定することも可能である。

tidybayes のジオムの slabinterval ファミリーには dots と dotsinterval ファミリーがあり、点描画の適切なビンサイズを自動的に決定し、サンプルから分位を計算して分位点描画を作成できる。ggdist::stat_dotsinterval() はサンプルに使用するために設計された水平バリアントである。

ABC %>%
  data_grid(condition) %>%
  add_epred_draws(m) %>%
  ggplot(aes(x = .epred, y = condition)) +
  stat_dotsinterval(quantiles = 100)

つまり、事後情報を1つの標準的な点または区間として考えるのではなく、 (例えば) 100のほぼ等しい可能性のある点として表現することである。

事後予測

ここで、add_epred_draws() は rstanarm::posterior_epred() に類似しており、add_predicted_draws() は rstanarm::posterior_predict() に類似しており、事後予測分布からのドローを与えている。

tidybayes::stat_interval() を使って、データと平均の事後分布と並べて予測バンドをプロットすることができる。

grid = ABC %>%
  data_grid(condition)

means = grid %>%
  add_epred_draws(m)

preds = grid %>%
  add_predicted_draws(m)

ABC %>%
  ggplot(aes(y = condition, x = response)) +
  stat_interval(aes(x = .prediction), data = preds) +
  stat_pointinterval(aes(x = .epred), data = means, .width = c(.66, .95), position = position_nudge(y = -0.3)) +
  geom_point() +
  scale_color_brewer()

フィット/予測曲線

不確実性を伴う適合曲線の描画を実証するために、mtcars のデータセットの一部に、少し素朴なモデルを適合させてみよう。

m_mpg = stan_glm(mpg ~ hp * cyl, data = mtcars)

フィットカーブを確率バンドでプロットすることができる。

mtcars %>%
  group_by(cyl) %>%
  data_grid(hp = seq_range(hp, n = 51)) %>%
  add_epred_draws(m_mpg) %>%
  ggplot(aes(x = hp, y = mpg, color = ordered(cyl))) +
  stat_lineribbon(aes(y = .epred)) +
  geom_point(data = mtcars) +
  scale_fill_brewer(palette = "Greys") +
  scale_color_brewer(palette = "Set2")

## Warning: Using the `size` aesthietic with geom_ribbon was deprecated in ggplot2 3.4.0.
## ℹ Please use the `linewidth` aesthetic instead.

## Warning: Unknown or uninitialised column: `linewidth`.

## Warning: Using the `size` aesthietic with geom_line was deprecated in ggplot2 3.4.0.
## ℹ Please use the `linewidth` aesthetic instead.

## Warning: Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.

あるいは、適当な数のフィットラインをサンプリングして (たとえば100本) 、それらをオーバープロットすることもできる。

mtcars %>%
  group_by(cyl) %>%
  data_grid(hp = seq_range(hp, n = 101)) %>%
  # NOTE: this shows the use of ndraws to subsample within add_epred_draws()
  # ONLY do this IF you are planning to make spaghetti plots, etc.
  # NEVER subsample to a small sample to plot intervals, densities, etc.
  add_epred_draws(m_mpg, ndraws = 100) %>%
  ggplot(aes(x = hp, y = mpg, color = ordered(cyl))) +
  geom_line(aes(y = .epred, group = paste(cyl, .draw)), alpha = .1) +
  geom_point(data = mtcars) +
  scale_color_brewer(palette = "Dark2")

あるいは、 (平均値ではなく) 事後予測値をプロットすることもできる。この例では また、重なり合ったバンドを見やすくするために、alpha を使用する。

mtcars %>%
  group_by(cyl) %>%
  data_grid(hp = seq_range(hp, n = 101)) %>%
  add_predicted_draws(m_mpg) %>%
  ggplot(aes(x = hp, y = mpg, color = ordered(cyl), fill = ordered(cyl))) +
  stat_lineribbon(aes(y = .prediction), .width = c(.95, .80, .50), alpha = 1/4) +
  geom_point(data = mtcars) +
  scale_fill_brewer(palette = "Set2") +
  scale_color_brewer(palette = "Dark2")

## Warning: Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.

その他のフィットラインの例については、vignette("tidy-brms") を参照。
animated hypothetical outcome plots (HOPs)。

因子の水準を比較する

各条件の平均を比較したい場合、compare_levels() は、ある因子のレベル間のある変数の値の比較を容易にする。デフォルトでは、一対の差分をすべて計算する。

ggdist::stat_halfeye() を使って compare_levels() をデモしてみよう。また 差の平均値で再注文する。

m %>%
  spread_draws(b[,,condition], sep = "[, :]") %>%
  compare_levels(b, by = condition) %>%
  ungroup() %>%
  mutate(condition = reorder(condition, b)) %>%
  ggplot(aes(y = condition, x = b)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linetype = "dashed") 

Ordinal models

data(esoph)
m_esoph_rs = stan_polr(tobgp ~ agegp, data = esoph, prior = R2(0.25), prior_counts = rstanarm::dirichlet(1))

esoph %>%
  data_grid(agegp) %>%
  add_linpred_draws(m_esoph_rs) %>%
  ggplot(aes(x = as.numeric(agegp), y = .linpred)) +
  stat_lineribbon() +
  scale_fill_brewer(palette = "Greys")

## Warning: Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.
## Unknown or uninitialised column: `linewidth`.

get_variables(m_esoph_rs)

##  [1] "agegp.L"        "agegp.Q"        "agegp.C"        "agegp^4"        "agegp^5"        "0-9g/day|10-19"
##  [7] "10-19|20-29"    "20-29|30+"      "accept_stat__"  "stepsize__"     "treedepth__"    "n_leapfrog__"  
## [13] "divergent__"    "energy__"

thresholds = m_esoph_rs %>%
  gather_draws(`.*[|].*`, regex = TRUE) %>%
  group_by(.draw) %>%
  select(.draw, threshold = .value) %>%
  summarise_all(list) %>%
  mutate(threshold = map(threshold, ~ c(., Inf)))

head(thresholds, 10)

thresholds[1,]$threshold

## [[1]]
## [1] -0.8200002  0.6020005  1.4201284        Inf

esoph %>%
  data_grid(agegp) %>%
  add_linpred_draws(m_esoph_rs) %>%
  inner_join(thresholds, by = ".draw") %>%
  mutate(`P(Y = category)` = map2(threshold, .linpred, function(alpha, beta_x)
      # this part is logit^-1(alpha_j - beta*x) - logit^-1(alpha_j-1 - beta*x)
      plogis(alpha - beta_x) - 
      plogis(lag(alpha, default = -Inf) - beta_x)
    )) %>%
  mutate(.category = list(levels(esoph$tobgp))) %>%
  unnest(c(threshold, `P(Y = category)`, .category)) %>%
  ggplot(aes(x = agegp, y = `P(Y = category)`, color = .category)) +
  stat_pointinterval(position = position_dodge(width = .4)) +
  scale_size_continuous(guide = "none") +
  scale_color_manual(values = brewer.pal(6, "Blues")[-c(1,2)]) 

esoph_plot = esoph %>%
  data_grid(agegp) %>%
  add_linpred_draws(m_esoph_rs) %>%
  inner_join(thresholds, by = ".draw") %>%
  mutate(`P(Y = category)` = map2(threshold, .linpred, function(alpha, beta_x)
      # this part is logit^-1(alpha_j - beta*x) - logit^-1(alpha_j-1 - beta*x)
      plogis(alpha - beta_x) - 
      plogis(lag(alpha, default = -Inf) - beta_x)
    )) %>%
  mutate(.category = list(levels(esoph$tobgp))) %>%
  unnest(c(threshold, `P(Y = category)`, .category)) %>%
  ggplot(aes(x = `P(Y = category)`, y = .category)) +
  coord_cartesian(expand = FALSE) +
  facet_grid(. ~ agegp, switch = "x") +
  theme_classic() +
  theme(strip.background = element_blank(), strip.placement = "outside") +
  ggtitle("P(tobacco consumption category | age group)") +
  xlab("age group")

esoph_plot +
  stat_summary(fun = median, geom = "bar", fill = "gray65", width = 1, color = "white") +
  stat_pointinterval()

esoph_plot +
  stat_ccdfinterval() +
  expand_limits(x = 0) #ensure bars go to 0