Mengekstrak dan memvisualisasikan tidy draws dari model brms
Copy Right @ Prof Dr Suhartono M.Kom @ Magister Informatika UIN Maulana Malik Ibrahim@ (http://suhartono.lecturer.uin-malang.ac.id/riwayat-hidup/)
18 Januari 2021
Pengantar
vignette menjelaskan cara menggunakan paket tidybayes untuk mengekstraksi tidy data frame dari model Bbrms, dan penarikan dari distribusi posterior variabel model, kecocokan, dan prediksi dari brms :: brm. Untuk pengenalan yang lebih umum tentang tidybayes, lihat vignette(“tidybayes”).
Untuk menjalankan vignette membutuhkan package sebagai berikut :
library(magrittr)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(purrr)##
## Attaching package: 'purrr'## The following object is masked from 'package:magrittr':
##
## set_nameslibrary(forcats)
library(tidyr)##
## Attaching package: 'tidyr'## The following object is masked from 'package:magrittr':
##
## extractlibrary(modelr)
library(ggdist)
library(tidybayes)
library(ggplot2)
library(cowplot)
library(rstan)## Loading required package: StanHeaders## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)##
## Attaching package: 'rstan'## The following object is masked from 'package:tidyr':
##
## extract## The following object is masked from 'package:magrittr':
##
## extractlibrary(brms)## Loading required package: Rcpp## Loading 'brms' package (version 2.14.4). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attaching package: 'brms'## The following object is masked from 'package:rstan':
##
## loo## The following objects are masked from 'package:tidybayes':
##
## dstudent_t, pstudent_t, qstudent_t, rstudent_t## The following objects are masked from 'package:ggdist':
##
## dstudent_t, pstudent_t, qstudent_t, rstudent_t## The following object is masked from 'package:stats':
##
## arlibrary(ggrepel)
library(RColorBrewer)
library(gganimate)## No renderer backend detected. gganimate will default to writing frames to separate files
## Consider installing:
## - the `gifski` package for gif output
## - the `av` package for video output
## and restarting the R sessiontheme_set(theme_tidybayes() + panel_border())Opsi ini membantu Stan bekerja lebih cepat:
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())Untuk mendemonstrasikan tidybayes, kami akan menggunakan dataset sederhana dengan 10 observasi dari 5 kondisi masing-masing:
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)
)Data cens_df dapat dilihat sebanyak 10 sebagai berikut :
head(ABC,10)## # A tibble: 10 x 2
## condition response
## <chr> <dbl>
## 1 A -0.420
## 2 B 1.69
## 3 C 1.37
## 4 D 1.04
## 5 E -0.144
## 6 A -0.301
## 7 B 0.764
## 8 C 1.68
## 9 D 0.857
## 10 E -0.931Ini adalah tipikal dari data frame format tidy : satu pengamatan per baris. Yang dilihat secara grafis:
ABC %>%
ggplot(aes(y = condition, x = response)) +
geom_point()
Mari kita menyesuaikan model hierarki dengan penyusutan menuju rata-rata global:
m = brm(
response ~ (1|condition),
data = ABC,
prior = c(
prior(normal(0, 1), class = Intercept),
prior(student_t(3, 0, 1), class = sd),
prior(student_t(3, 0, 1), class = sigma)
),
control = list(adapt_delta = .99),
)## Compiling Stan program...## Trying to compile a simple C file## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## gcc -std=gnu99 -I"/usr/share/R/include" -DNDEBUG -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/Rcpp/include/" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/unsupported" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/BH/include" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/StanHeaders/include/src/" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/StanHeaders/include/" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppParallel/include/" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -fpic -g -O2 -fdebug-prefix-map=/build/r-base-V28x5H/r-base-3.6.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -g -c foo.c -o foo.o
## In file included from /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/Eigen/Core:88:0,
## from /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/Eigen/Dense:1,
## from /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13,
## from <command-line>:0:
## /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name ‘namespace’
## namespace Eigen {
## ^~~~~~~~~
## /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:17: error: expected ‘=’, ‘,’, ‘;’, ‘asm’ or ‘__attribute__’ before ‘{’ token
## namespace Eigen {
## ^
## In file included from /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/Eigen/Dense:1:0,
## from /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13,
## from <command-line>:0:
## /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/Eigen/Core:96:10: fatal error: complex: No such file or directory
## #include <complex>
## ^~~~~~~~~
## compilation terminated.
## /usr/lib/R/etc/Makeconf:168: recipe for target 'foo.o' failed
## make: *** [foo.o] Error 1## Start samplingHasilnya data dapat terlihat seperti ini:
m## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: response ~ (1 | condition)
## Data: ABC (Number of observations: 50)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~condition (Number of levels: 5)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.16 0.43 0.61 2.22 1.00 946 1426
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.52 0.46 -0.43 1.44 1.00 891 1235
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.56 0.06 0.46 0.70 1.00 1820 2127
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Mengekstrak draws dari fit dalam tidy-format menggunakan spread_draws
Sekarang setelah kita mendapatkan hasil, maka hasil yang didapat dimulai: mengeluarkan undian dalam tidy format! Pertama, kita akan menggunakan fungsi get_variables() untuk mendapatkan daftar nama variabel model mentah sehingga kita tahu variabel apa yang dapat kita ekstrak dari model:
get_variables(m)## [1] "b_Intercept" "sd_condition__Intercept"
## [3] "sigma" "r_condition[A,Intercept]"
## [5] "r_condition[B,Intercept]" "r_condition[C,Intercept]"
## [7] "r_condition[D,Intercept]" "r_condition[E,Intercept]"
## [9] "lp__" "accept_stat__"
## [11] "stepsize__" "treedepth__"
## [13] "n_leapfrog__" "divergent__"
## [15] "energy__"Jadi r_condition [D, Intercept] memiliki indeks D dan Intercept, dan spread_draws () memungkinkan kita mengekstrak indeks ini sebagai kolom dalam tidy data frame yang dihasilkan dari penarikan dari r_condition:
m %>%
spread_draws(r_condition[condition,term]) %>%
head(10)## # A tibble: 10 x 6
## # Groups: condition, term [1]
## condition term r_condition .chain .iteration .draw
## <chr> <chr> <dbl> <int> <int> <int>
## 1 A Intercept -0.544 1 1 1
## 2 A Intercept -0.240 1 2 2
## 3 A Intercept -0.392 1 3 3
## 4 A Intercept -0.494 1 4 4
## 5 A Intercept 0.243 1 5 5
## 6 A Intercept 0.173 1 6 6
## 7 A Intercept -0.356 1 7 7
## 8 A Intercept -0.197 1 8 8
## 9 A Intercept -0.0267 1 9 9
## 10 A Intercept -0.930 1 10 10Kita dapat memilih nama apapun yang kita inginkan untuk kolom indeks; misalnya.:
m %>%
spread_draws(r_condition[c,t]) %>%
head(10)## # A tibble: 10 x 6
## # Groups: c, t [1]
## c t r_condition .chain .iteration .draw
## <chr> <chr> <dbl> <int> <int> <int>
## 1 A Intercept -0.544 1 1 1
## 2 A Intercept -0.240 1 2 2
## 3 A Intercept -0.392 1 3 3
## 4 A Intercept -0.494 1 4 4
## 5 A Intercept 0.243 1 5 5
## 6 A Intercept 0.173 1 6 6
## 7 A Intercept -0.356 1 7 7
## 8 A Intercept -0.197 1 8 8
## 9 A Intercept -0.0267 1 9 9
## 10 A Intercept -0.930 1 10 10Tetapi nama yang lebih deskriptif dan tidak terlalu samar dari contoh sebelumnya mungkin lebih disukai.
Dalam model khusus ini, hanya ada satu istilah (Intercept), sehingga kita dapat menghilangkan indeks tersebut secara bersamaan untuk mendapatkan setiap kondisi dan nilai r_condition untuk kondisi tersebut:
m %>%
spread_draws(r_condition[condition,]) %>%
head(10)## # A tibble: 10 x 5
## # Groups: condition [1]
## condition r_condition .chain .iteration .draw
## <chr> <dbl> <int> <int> <int>
## 1 A -0.544 1 1 1
## 2 A -0.240 1 2 2
## 3 A -0.392 1 3 3
## 4 A -0.494 1 4 4
## 5 A 0.243 1 5 5
## 6 A 0.173 1 6 6
## 7 A -0.356 1 7 7
## 8 A -0.197 1 8 8
## 9 A -0.0267 1 9 9
## 10 A -0.930 1 10 10Ringkasan dan interval poin
Dengan variabel model sederhana
Misalnya, kita mungkin mengekstrak hasil imbang yang sesuai dengan distribusi posterior dari mean keseluruhan dan deviasi standar pengamatan:
m %>%
spread_draws(b_Intercept, sigma) %>%
head(10)## # A tibble: 10 x 5
## .chain .iteration .draw b_Intercept sigma
## <int> <int> <int> <dbl> <dbl>
## 1 1 1 1 0.615 0.517
## 2 1 2 2 0.724 0.634
## 3 1 3 3 0.753 0.614
## 4 1 4 4 0.567 0.673
## 5 1 5 5 0.0862 0.611
## 6 1 6 6 0.160 0.629
## 7 1 7 7 0.455 0.555
## 8 1 8 8 0.400 0.558
## 9 1 9 9 0.164 0.531
## 10 1 10 10 1.09 0.639Seperti dengan r_condition [kondisi, istilah], ini memberi kita bingkai data yang rapi. Jika kita menginginkan median dan interval kuantil 95% dari variabel, kita dapat menerapkan median_qi():
m %>%
spread_draws(b_Intercept, sigma) %>%
median_qi(b_Intercept, sigma)## # A tibble: 1 x 9
## b_Intercept b_Intercept.low… b_Intercept.upp… sigma sigma.lower sigma.upper
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.535 -0.432 1.44 0.558 0.459 0.699
## # … with 3 more variables: .width <dbl>, .point <chr>, .interval <chr>Jika Anda lebih suka memiliki daftar interval format panjang, gunakan gathering_draws () sebagai gantinya:
m %>%
gather_draws(b_Intercept, sigma) %>%
median_qi()## # A tibble: 2 x 7
## .variable .value .lower .upper .width .point .interval
## <chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 b_Intercept 0.535 -0.432 1.44 0.95 median qi
## 2 sigma 0.558 0.459 0.699 0.95 median qiRingkasan dan interval titik plot
Merencanakan ringkasan titik dan dengan satu interval langsung menggunakan ggdist :: geom_pointinterval ():
m %>%
spread_draws(b_Intercept, r_condition[condition,]) %>%
median_qi(condition_mean = b_Intercept + r_condition) %>%
ggplot(aes(y = condition, x = condition_mean, xmin = .lower, xmax = .upper)) +
geom_pointinterval()
Interval dengan beberapa tingkat probabilitas
median_qi () dan fungsi kembarnya juga dapat menghasilkan sejumlah interval probabilitas dengan menyetel argumen .width =argument:
m %>%
spread_draws(b_Intercept, r_condition[condition,]) %>%
median_qi(condition_mean = b_Intercept + r_condition, .width = c(.95, .8, .5))## # A tibble: 15 x 7
## condition condition_mean .lower .upper .width .point .interval
## <chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 A 0.197 -0.145 0.534 0.95 median qi
## 2 B 0.998 0.661 1.34 0.95 median qi
## 3 C 1.83 1.48 2.19 0.95 median qi
## 4 D 1.01 0.667 1.36 0.95 median qi
## 5 E -0.890 -1.23 -0.532 0.95 median qi
## 6 A 0.197 -0.0275 0.413 0.8 median qi
## 7 B 0.998 0.778 1.22 0.8 median qi
## 8 C 1.83 1.59 2.06 0.8 median qi
## 9 D 1.01 0.784 1.24 0.8 median qi
## 10 E -0.890 -1.11 -0.659 0.8 median qi
## 11 A 0.197 0.0767 0.311 0.5 median qi
## 12 B 0.998 0.882 1.12 0.5 median qi
## 13 C 1.83 1.71 1.95 0.5 median qi
## 14 D 1.01 0.895 1.13 0.5 median qi
## 15 E -0.890 -1.01 -0.774 0.5 median qiHasilnya dalam format rapi: satu baris per grup dan lebar interval ketidakpastian (.width). Ini memfasilitasi pembuatan plot. Misalnya, menetapkan -.width ke estetika ukuran akan menampilkan semua interval, membuat garis yang lebih tebal sesuai dengan interval yang lebih kecil. Geom ggdist :: geom_pointinterval () secara otomatis menyetel estetika ukuran dengan tepat berdasarkan kolom .width dalam data untuk menghasilkan plot titik dengan beberapa tingkat probabilitas:
m %>%
spread_draws(b_Intercept, r_condition[condition,]) %>%
median_qi(condition_mean = b_Intercept + r_condition, .width = c(.95, .66)) %>%
ggplot(aes(y = condition, x = condition_mean, xmin = .lower, xmax = .upper)) +
geom_pointinterval() 
Interval dengan kepadatan
Untuk melihat kepadatan bersama dengan intervalnya, kita dapat menggunakan ggdist :: stat_eye () (“plot mata”, yang menggabungkan interval dengan plot biola), atau ggdist :: stat_halfeye () (interval + plot kepadatan):
m %>%
spread_draws(b_Intercept, r_condition[condition,]) %>%
mutate(condition_mean = b_Intercept + r_condition) %>%
ggplot(aes(y = condition, x = condition_mean)) +
stat_halfeye()
Atau katakanlah Anda ingin memberi keterangan pada bagian kepadatan warna; estetika isian dapat bervariasi dalam slab di semua geom dan statistik dalam keluarga ggdist :: geom_slabinterval (), termasuk ggdist :: stat_halfeye (). Misalnya, jika Anda ingin memberi anotasi pada wilayah spesifik domain kesetaraan praktis (ROPE), Anda dapat melakukan sesuatu seperti ini:
m %>%
spread_draws(b_Intercept, r_condition[condition,]) %>%
mutate(condition_mean = b_Intercept + r_condition) %>%
ggplot(aes(y = condition, 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"))
Titik-titik kuantitatif
Intervalnya bagus jika level alpha kebetulan sejalan dengan keputusan apa pun yang Anda coba buat, tetapi mendapatkan bentuk posterior lebih baik (karena itu plot mata, di atas). Di sisi lain, membuat kesimpulan dari plot kepadatan tidak tepat (memperkirakan luas satu bentuk sebagai proporsi bentuk lainnya adalah tugas perseptual yang sulit). Penalaran tentang probabilitas dalam format frekuensi lebih mudah, memotivasi dotplots kuantil (Kay et al. 2016, Fernandes et al. 2018), yang juga memungkinkan estimasi yang tepat dari interval arbitrer (hingga resolusi titik plot, 100 dalam contoh di bawah) .
Dalam keluarga geom slabinterval di tidybayes adalah keluarga titik dan titik-titik, yang secara otomatis menentukan ukuran bin yang sesuai untuk dotplots dan dapat menghitung kuantil dari sampel untuk membangun dotplots kuantil. ggdist :: stat_dots () adalah varian yang dirancang untuk digunakan pada sampel:
ABC %>%
data_grid(condition) %>%
add_fitted_draws(m) %>%
ggplot(aes(x = .value, y = condition)) +
stat_dots(quantiles = 100)
Prediksi posterior
Dimana add_fitted_draws () dianalogikan dengan brms :: pas.brmsfit () (atau brms :: posterior_linpred ()), add_predicted_draws () analog dengan brms :: predict.brmsfit () (atau brms :: posterior_predict ()), memberi menarik dari distribusi prediktif posterior.
Berikut adalah contoh distribusi prediksi posterior yang diplot menggunakan ggdist :: stat_slab ():
ABC %>%
data_grid(condition) %>%
add_predicted_draws(m) %>%
ggplot(aes(x = .prediction, y = condition)) +
stat_slab()
Kita juga bisa menggunakan ggdist :: stat_interval () untuk memplot band prediktif di samping data:
ABC %>%
data_grid(condition) %>%
add_predicted_draws(m) %>%
ggplot(aes(y = condition, x = .prediction)) +
stat_interval(.width = c(.50, .80, .95, .99)) +
geom_point(aes(x = response), data = ABC) +
scale_color_brewer()
Secara keseluruhan, data, prediksi posterior, dan distribusi posterior sarana:
grid = ABC %>%
data_grid(condition)
fits = grid %>%
add_fitted_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 = .value), data = fits, .width = c(.66, .95), position = position_nudge(y = -0.3)) +
geom_point() +
scale_color_brewer()
Prediksi posterior, gaya Kruschke
Pendekatan prediksi posterior di atas terintegrasi di atas ketidakpastian parameter untuk memberikan distribusi prediksi posterior tunggal. Pendekatan lain, yang sering digunakan oleh John Kruschke dalam bukunya Doing Bayesian Data Analysis, adalah mencoba menunjukkan ketidakpastian prediktif dan ketidakpastian parameter secara bersamaan dengan menunjukkan beberapa kemungkinan distribusi prediksi yang disiratkan oleh posterior.
Kita dapat melakukan ini dengan mudah dengan menanyakan parameter distribusi untuk prediksi tertentu yang diimplikasikan oleh posterior. Kami akan melakukannya secara eksplisit di sini dengan menyetel dpar = c (“mu”, “sigma”) di add_fitted_draws (). Daripada menentukan parameter secara eksplisit, Anda juga dapat menyetel dpar = TRUE untuk mendapatkan penarikan dari semua parameter distribusi dalam model, dan ini akan berfungsi untuk distribusi respons apa pun yang didukung oleh brms. Kemudian, kita dapat memilih sejumlah kecil penarikan menggunakan sample_draws () dan kemudian menggunakan ggdist :: stat_dist_slab () untuk memvisualisasikan setiap distribusi prediksi yang tersirat oleh nilai mu dan sigma:
ABC %>%
data_grid(condition) %>%
add_fitted_draws(m, dpar = c("mu", "sigma")) %>%
sample_draws(30) %>%
ggplot(aes(y = condition)) +
stat_dist_slab(aes(dist = "norm", arg1 = mu, arg2 = sigma),
slab_color = "gray65", alpha = 1/10, fill = NA
) +
geom_point(aes(x = response), data = ABC, shape = 21, fill = "#9ECAE1", size = 2)
Untuk penjelasan lebih rinci tentang bagan ini (dan beberapa variasinya yang berguna), lihat entri blog yang sangat baik dari Solomon Kurz tentang topik tersebut.
Kami bahkan dapat menggabungkan plot distribusi prediktif bergaya Kruschke dengan setengah mata yang menunjukkan cara posterior:
ABC %>%
data_grid(condition) %>%
add_fitted_draws(m, dpar = c("mu", "sigma")) %>%
ggplot(aes(x = condition)) +
stat_dist_slab(aes(dist = "norm", arg1 = mu, arg2 = sigma),
slab_color = "gray65", alpha = 1/10, fill = NA, data = . %>% sample_draws(30), scale = .5
) +
stat_halfeye(aes(y = .value), side = "bottom", scale = .5) +
geom_point(aes(y = response), data = ABC, shape = 21, fill = "#9ECAE1", size = 2, position = position_nudge(x = -.2))
Daftar pustaka :