# More support for Bayesian analysis in the sj!-packages #rstats #rstan #brms

**R – Strenge Jacke!**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Another quick preview of my R-packages, especially **sjPlot**, which now also support `brmsfit`

-objects from the great brms-package. To demonstrate the new features, I load all my „core“-packages at once, using the *strengejacke*-package, which is only available from GitHub. This package simply loads four packages (*sjlabelled, sjmisc, sjstats and sjPlot*).

First, I fit two sample models, one with **brms** and one with **rstanarm**:

<span style="color:#999999;">#> install pkg "strengejacke" via GitHub: #> devtools::install_github("strengejacke/strengejacke")</span> library(strengejacke) library(ggplot2) library(lme4) library(glmmTMB) library(brms) library(rstanarm) data(Owls) m1 <- brm( SiblingNegotiation ~ FoodTreatment + ArrivalTime + SexParent + (1 | Nest), data = Owls, family = zero_inflated_poisson(link = "log", link_zi = "logit") ) m2 <- stan_glmer.nb( SiblingNegotiation ~ FoodTreatment + ArrivalTime + SexParent + (1 | Nest), data = Owls ) m3 <- stan_lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)

Next, some summary statistics, at first simply the *High Density Interval*, obtained from the `hdi()`

-function. You can define the probabilites of the interval, however, here I simply use the defaults (which is the 90% HDI).

hdi(m1) <span style="color:#999999;">#> # A tibble: 6 x 3 #> term hdi.low hdi.high #> #> 1 b_Intercept 3.798693e+00 4.584998e+00 #> 2 b_FoodTreatmentSatiated -2.795021e-01 -1.586102e-01 #> 3 b_ArrivalTime -9.768530e-02 -6.616026e-02 #> 4 b_SexParentMale -7.766154e-02 3.794376e-02 #> 5 zi 2.304234e-01 2.900676e-01 #> 6 lp__ -2.020651e+03 -2.002387e+03</span> hdi(m2) <span style="color:#999999;">#> # A tibble: 5 x 3 #> term hdi.low hdi.high #> #> 1 (Intercept) 3.82688888 5.93003148 #> 2 FoodTreatmentSatiated -0.86327579 -0.51376279 #> 3 ArrivalTime -0.15414266 -0.07134095 #> 4 SexParentMale -0.09371495 0.24440447 #> 5 reciprocal_dispersion 0.77345433 1.00889671</span>

Next, I implemented an own tidyr-function, similar to **broom’s** `tidy()`

, called `tidy_stan()`

. Unlike broom’s tidy-function, `tidy_stan()`

computes the HDI instead of credibility intervals, and you may return fixed or random effects only, or both. Furthermore, `tidy_stan()`

also returns the ratio of effective numbers of samples, *n_eff*, and Rhat statistics.

tidy_stan(m1) <span style="color:#999999;">#> # A tibble: 6 x 7 #> term estimate std.error hdi.low hdi.high n_eff Rhat #> #> 1 b_Intercept 4.206 0.241 3.819 4.581 1.000 1.000 #> 2 b_FoodTreatmentSatiated -0.220 0.037 -0.279 -0.162 1.000 1.000 #> 3 b_ArrivalTime -0.082 0.009 -0.096 -0.065 1.000 1.000 #> 4 b_SexParentMale -0.021 0.037 -0.077 0.036 1.000 1.000 #> 5 zi 0.259 0.019 0.230 0.288 1.000 0.999 #> 6 lp__ -2011.050 5.537 -2020.651 -2002.925 0.186 1.001</span> tidy_stan(m2) <span style="color:#999999;">#> # A tibble: 5 x 7 #> term estimate std.error hdi.low hdi.high n_eff Rhat #> #> 1 (Intercept) 4.873 0.652 3.847 5.880 1 1.000 #> 2 FoodTreatmentSatiated -0.693 0.112 -0.860 -0.522 1 1.000 #> 3 ArrivalTime -0.115 0.025 -0.154 -0.074 1 1.000 #> 4 SexParentMale 0.071 0.102 -0.088 0.241 1 0.999 #> 5 reciprocal_dispersion 0.883 0.073 0.776 1.005 1 1.000</span>

Finally, the new and generic `plot_model()`

function in **sjPlot** also supports brmsfit or stanreg objects, and plots estimates, random effects and marginal effects.

The coefficients-plot was already shown in a previous post. It shows the 50% and 89% HDI and the posterior median. The style can be changed using the `bpe`

and `bpe.style`

arguments.

theme_set(theme_sjplot()) plot_model(m1, axis.lim = c(.6, 1.4)) plot_model(m2, bpe = "mean", bpe.style = "dot", colors = "grey30")

Random effects are displayed in a similar way. For models with random slopes and intercepts, random effects are displayed in a grid layout (use `grid = FALSE`

to create a separate plot for each random effect).

plot_model(m1, type = "re") plot_model(m3, type = "re")

**Random effects, rstanarm-model**

Marginal effects in **sjPlot** are based on the ggeffects-package. By default, predictions are computed with `rstantools::posterior_linpred()`

.

plot_model(m1, type = "pred", terms = c("ArrivalTime", "FoodTreatment")) plot_model(m2, type = "pred", terms = c("ArrivalTime", "FoodTreatment"))

**Marginal effects, rstanarm-model**

The work on sjPlot for the next CRAN release is almost done… I hope to submit the sjPlot-update in the course of the next week! In the meantime, you may try out the new features using the GitHub-versions of my R-packages – at least *sjstats* and *ggeffects* are required.

Tagged: brms, data visualization, ggplot, R, rstan, rstanarm, rstats, sjPlot, Stan

**leave a comment**for the author, please follow the link and comment on their blog:

**R – Strenge Jacke!**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.