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

October 11, 2017
By

(This article was first published on 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")

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.

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

To 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...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)