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

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")

Coefficients, brms-model

Coefficients, rstanarm-model

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, brms-model

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, brms-model

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.