**R – Strenge Jacke!**, and kindly contributed to R-bloggers)

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**:

```
#> install pkg "strengejacke" via GitHub:
#> devtools::install_github("strengejacke/strengejacke")
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)
#> # 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
hdi(m2)
#> # 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
```

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)
#> # 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
tidy_stan(m2)
#> # 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
```

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 on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...