Bayesian analysis with brms

The marginaleffects package offers convenience functions to compute and display predictions, contrasts, and marginal effects from bayesian models estimated by the brms package. To compute these quantities, marginaleffects relies on workshorse functions from the brms package to draw from the posterior distribution. The type of draws used is controlled by using the type argument of the predictions or marginaleffects functions:

type = "response": Compute posterior draws of the expected value using the brms::posterior_epred function.
type = "link": Compute posterior draws of the linear predictor using the brms::posterior_linpred function.
type = "prediction": Compute posterior draws of the posterior predictive distribution using the brms::posterior_predict function.

The predictions and marginaleffects functions can also pass additional arguments to the brms prediction functions via the ... ellipsis. For example, if mod is a mixed-effects model, then this command will compute 10 draws from the posterior predictive distribution, while ignoring all group-level effects:

predictions(mod, type = "prediction", ndraws = 10, re_formula = NA)

See the brms documentation for a list of available arguments:

?brms::posterior_epred
?brms::posterior_linpred
?brms::posterior_predict

Logistic regression with multiplicative interactions

Load libraries and download data on passengers of the Titanic from the Rdatasets archive:

library(marginaleffects)
library(brms)
library(ggplot2)
library(ggdist)
library(magrittr)

dat <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/carData/TitanicSurvival.csv")
dat$survived <- ifelse(dat$survived == "yes", 1, 0)
dat$woman <- ifelse(dat$sex == "female", 1, 0)

Fit a logit model with a multiplicative interaction:

mod <- brm(survived ~ woman * age + passengerClass,
           family = bernoulli(link = "logit"),
           data = dat)

Adjusted predictions

We can compute adjusted predicted values of the outcome variable (i.e., probability of survival aboard the Titanic) using the predictions function. By default, this function calculates predictions for each row of the dataset:

pred <- predictions(mod)
head(pred)
#>   rowid     type predicted survived woman     age passengerClass  conf.low
#> 1     1 response 0.9364708        1     1 29.0000            1st 0.9069508
#> 2     2 response 0.8478001        1     0  0.9167            1st 0.7416961
#> 3     3 response 0.9426796        0     1  2.0000            1st 0.8946092
#> 4     4 response 0.5134935        0     0 30.0000            1st 0.4294245
#> 5     5 response 0.9373910        0     1 25.0000            1st 0.9076207
#> 6     6 response 0.2737386        1     0 48.0000            1st 0.2029320
#>   conf.high
#> 1 0.9591722
#> 2 0.9175600
#> 3 0.9704055
#> 4 0.5991923
#> 5 0.9602311
#> 6 0.3540721

To visualize the relationship between the outcome and one of the regressors, we can plot conditional adjusted predictions with the plot_cap function:

plot_cap(mod, condition = "age")

Compute adjusted predictions for some user-specified values of the regressors, using the newdata argument and the datagrid function:

pred <- predictions(mod,
                    newdata = datagrid(woman = 0:1,
                                       passengerClass = c("1st", "2nd", "3rd")))
pred
#>   rowid     type  predicted      age woman passengerClass   conf.low conf.high
#> 1     1 response 0.51524219 29.88113     0            1st 0.43099807 0.6012692
#> 2     2 response 0.20199559 29.88113     0            2nd 0.15214184 0.2607661
#> 3     3 response 0.08798361 29.88113     0            3rd 0.06563742 0.1148200
#> 4     4 response 0.93626465 29.88113     1            1st 0.90643165 0.9590551
#> 5     5 response 0.77759974 29.88113     1            2nd 0.71026516 0.8334703
#> 6     6 response 0.57183081 29.88113     1            3rd 0.49836765 0.6413378

The posteriordraws function samples from the posterior distribution of the model, and produces a data frame with drawid and draw columns.

pred <- posteriordraws(pred)
head(pred)
#>   drawid       draw rowid     type  predicted      age woman passengerClass
#> 1      1 0.53605685     1 response 0.51524219 29.88113     0            1st
#> 2      1 0.21803613     2 response 0.20199559 29.88113     0            2nd
#> 3      1 0.08772301     3 response 0.08798361 29.88113     0            3rd
#> 4      1 0.94042848     4 response 0.93626465 29.88113     1            1st
#> 5      1 0.79208379     5 response 0.77759974 29.88113     1            2nd
#> 6      1 0.56780970     6 response 0.57183081 29.88113     1            3rd
#>     conf.low conf.high
#> 1 0.43099807 0.6012692
#> 2 0.15214184 0.2607661
#> 3 0.06563742 0.1148200
#> 4 0.90643165 0.9590551
#> 5 0.71026516 0.8334703
#> 6 0.49836765 0.6413378

This “long” format makes it easy to plots results:

ggplot(pred, aes(x = draw, fill = factor(woman))) +
    geom_density() +
    facet_grid(~ passengerClass, labeller = label_both) +
    labs(x = "Predicted probability of survival", y = "", fill = "Woman")

Marginal effects

Use marginaleffects() to compute marginal effects (slopes of the regression equation) for each row of the dataset, and use summary() to compute “Average Marginal Effects”, that is, the average of all observation-level marginal effects:

mfx <- marginaleffects(mod)
summary(mfx)
#> Average marginal effects 
#>             Term  Contrast    Effect     2.5 %   97.5 %
#> 1          woman     dY/dX  0.365918  0.337130  0.39163
#> 2            age     dY/dX -0.005255 -0.007087 -0.00344
#> 3 passengerClass 2nd - 1st -0.237235 -0.310456 -0.16329
#> 4 passengerClass 3rd - 1st -0.388029 -0.452463 -0.31987
#> 
#> Model type:  brmsfit 
#> Prediction type:  response

Compute marginal effects with some regressors fixed at user-specified values, and other regressors held at their means:

marginaleffects(mod,
                newdata = datagrid(woman = 1,
                                   passengerClass = "1st"))
#>   rowid     type           term  contrast         dydx     conf.low
#> 1     1 response          woman     dY/dX  0.156651021  0.110625470
#> 2     1 response            age     dY/dX -0.000239209 -0.001409638
#> 3     1 response passengerClass 2nd - 1st -0.158062819 -0.222048775
#> 4     1 response passengerClass 3rd - 1st -0.364428580 -0.433930646
#>       conf.high      age woman passengerClass
#> 1  0.2083671441 29.88113     1            1st
#> 2  0.0008619362 29.88113     1            1st
#> 3 -0.1018392438 29.88113     1            1st
#> 4 -0.2969196694 29.88113     1            1st

Compute and plot conditional marginal effects:

plot_cme(mod, effect = "woman", condition = "age")

The posteriordraws produces a dataset with drawid and draw columns:

draws <- posteriordraws(mfx)

dim(draws)
#> [1] 16736000       13

head(draws)
#>   drawid       draw rowid     type  term contrast       dydx   conf.low
#> 1      1 0.14356982     1 response woman    dY/dX 0.15347608 0.10803140
#> 2      1 0.08189881     2 response woman    dY/dX 0.13727820 0.03838740
#> 3      1 0.04865343     3 response woman    dY/dX 0.05939775 0.02948603
#> 4      1 0.65251084     4 response woman    dY/dX 0.65397000 0.56801942
#> 5      1 0.13027082     5 response woman    dY/dX 0.13862164 0.09693914
#> 6      1 0.75171004     6 response woman    dY/dX 0.70703323 0.59064279
#>   conf.high survived woman     age passengerClass
#> 1 0.2046992        1     1 29.0000            1st
#> 2 0.3009955        1     0  0.9167            1st
#> 3 0.1022993        0     1  2.0000            1st
#> 4 0.7386564        0     0 30.0000            1st
#> 5 0.1885322        0     1 25.0000            1st
#> 6 0.8477174        1     0 48.0000            1st

We can use this dataset to plot our results. For example, to plot the posterior density of the marginal effect of age when the woman variable is equal to 0 or 1:

mfx <- marginaleffects(mod,
                       variables = "age",
                       newdata = datagrid(woman = 0:1)) |>
       posteriordraws()

ggplot(mfx, aes(x = draw, fill = factor(woman))) +
    stat_halfeye(slab_alpha = .5) +
    labs(x = "Marginal Effect of Age on Survival",
         y = "Posterior density",
         fill = "Woman")

Random effects model

This section replicates some of the analyses of a random effects model published in Andrew Heiss’ blog post: “A guide to correctly calculating posterior predictions and average marginal effects with multilievel Bayesian models.” The objective is mainly to illustrate the use of marginaleffects. Please refer to the original post for a detailed discussion of the quantities computed below.

Load libraries and download data:

library(brms)
library(ggdist)
library(patchwork)
library(marginaleffects)

vdem_2015 <- read.csv("https://github.com/vincentarelbundock/marginaleffects/raw/main/data-raw/vdem_2015.csv")

head(vdem_2015)
#>   country_name country_text_id year                           region
#> 1       Mexico             MEX 2015  Latin America and the Caribbean
#> 2     Suriname             SUR 2015  Latin America and the Caribbean
#> 3       Sweden             SWE 2015 Western Europe and North America
#> 4  Switzerland             CHE 2015 Western Europe and North America
#> 5        Ghana             GHA 2015               Sub-Saharan Africa
#> 6 South Africa             ZAF 2015               Sub-Saharan Africa
#>   media_index party_autonomy_ord polyarchy civil_liberties party_autonomy
#> 1       0.837                  3     0.631           0.704           TRUE
#> 2       0.883                  4     0.777           0.887           TRUE
#> 3       0.956                  4     0.915           0.968           TRUE
#> 4       0.939                  4     0.901           0.960           TRUE
#> 5       0.858                  4     0.724           0.921           TRUE
#> 6       0.898                  4     0.752           0.869           TRUE

Fit a basic model:

mod <- brm(
  bf(media_index ~ party_autonomy + civil_liberties + (1 | region),
     phi ~ (1 | region)),
  data = vdem_2015,
  family = Beta(),
  control = list(adapt_delta = 0.9))

Posterior predictions

To compute posterior predictions for specific values of the regressors, we use the newdata argument and the datagrid function. We also use the type argument to compute two types of predictions: accounting for residual (observation-level) residual variance (prediction) or ignoring it (response).

nd = datagrid(model = mod,
              party_autonomy = c(TRUE, FALSE),
              civil_liberties = .5,
              region = "Middle East and North Africa")
p1 <- predictions(mod, type = "response", newdata = nd) %>%
    posteriordraws()
p2 <- predictions(mod, type = "prediction", newdata = nd) %>%
    posteriordraws()
pred <- rbind(p1, p2)

Extract posterior draws and plot them:

ggplot(pred, aes(x = draw, fill = party_autonomy)) +
    stat_halfeye(alpha = .5) +
    facet_wrap(~ type) +
    labs(x = "Media index (predicted)", 
         y = "Posterior density",
         fill = "Party autonomy")

Marginal effects and contrasts

As noted in the Marginal Effects vignette, there should be one distinct marginal effect for each combination of regressor values. Here, we consider only one combination of regressor values, where region is “Middle East and North Africa”, and civil_liberties is 0.5. Then, we calculate the mean of the posterior distribution of marginal effects:

mfx <- marginaleffects(mod,
                       newdata = datagrid(civil_liberties = .5,
                                          region = "Middle East and North Africa"))
mfx
#>   rowid     type            term     contrast      dydx  conf.low conf.high
#> 1     1 response  party_autonomy TRUE - FALSE 0.2522815 0.1676294 0.3332248
#> 2     1 response civil_liberties        dY/dX 0.8166910 0.6271669 1.0091115
#>   party_autonomy civil_liberties                       region
#> 1           TRUE             0.5 Middle East and North Africa
#> 2           TRUE             0.5 Middle East and North Africa

Use the posteriordraws() to extract draws from the posterio distribution of marginal effects, and plot them:

mfx <- posteriordraws(mfx)

ggplot(mfx, aes(x = draw, y = term)) +
  stat_halfeye() +
  labs(x = "Marginal effect", y = "")

Plot marginal effects, conditional on a regressor:

plot_cme(mod,
         effect = "civil_liberties",
         condition = "party_autonomy")

Continuous predictors

pred <- predictions(mod,
                    newdata = datagrid(party_autonomy = FALSE,
                                       region = "Middle East and North Africa",
                                       civil_liberties = seq(0, 1, by = 0.05))) |>
        posteriordraws()

ggplot(pred, aes(x = civil_liberties, y = draw)) +
    stat_lineribbon() +
    scale_fill_brewer(palette = "Reds") +
    labs(x = "Civil liberties",
         y = "Media index (predicted)",
         fill = "")

The slope of this line for different values of civil liberties can be obtained with:

mfx <- marginaleffects(mod,
                       newdata = datagrid(civil_liberties = c(.2, .5, .8),
                                          party_autonomy = FALSE,
                                          region = "Middle East and North Africa"),
                       variables = "civil_liberties")
mfx
#>   rowid     type            term      dydx  conf.low conf.high civil_liberties
#> 1     1 response civil_liberties 0.4837312 0.3650440 0.6316934             0.2
#> 2     2 response civil_liberties 0.8034918 0.6149334 0.9927250             0.5
#> 3     3 response civil_liberties 0.8071646 0.6814830 0.9418585             0.8
#>   party_autonomy                       region
#> 1          FALSE Middle East and North Africa
#> 2          FALSE Middle East and North Africa
#> 3          FALSE Middle East and North Africa

And plotted:

mfx <- posteriordraws(mfx)

ggplot(mfx, aes(x = draw, fill = factor(civil_liberties))) +
    stat_halfeye(slab_alpha = .5) +
    labs(x = "Marginal effect of Civil Liberties on Media Index",
         y = "Posterior density",
         fill = "Civil liberties")

The marginaleffects function can use the ellipsis (...) to push any argument forward to the posterior_predict function. This can alter the types of predictions returned. For example, the re_formula=NA argument of the posterior_predict.brmsfit method will compute marginaleffects without including any group-level effects:

mfx <- marginaleffects(mod,
                       newdata = datagrid(civil_liberties = c(.2, .5, .8),
                                          party_autonomy = FALSE,
                                          region = "Middle East and North Africa"),
                       variables = "civil_liberties",
                       re_formula = NA) |>
       posteriordraws()

ggplot(mfx, aes(x = draw, fill = factor(civil_liberties))) +
    stat_halfeye(slab_alpha = .5) +
    labs(x = "Marginal effect of Civil Liberties on Media Index",
         y = "Posterior density",
         fill = "Civil liberties")

Global grand mean

pred <- predictions(mod,
                    re_formula = NA,
                    newdata = datagrid(party_autonomy = c(TRUE, FALSE))) |>
        posteriordraws()

mfx <- marginaleffects(mod,
                       re_formula = NA,
                       variables = "party_autonomy") |>
       posteriordraws()

plot1 <- ggplot(pred, aes(x = draw, fill = party_autonomy)) +
         stat_halfeye(slab_alpha = .5) +
         labs(x = "Media index (Predicted)",
              y = "Posterior density",
              fill = "Party autonomy")

plot2 <- ggplot(mfx, aes(x = draw)) +
         stat_halfeye(slab_alpha = .5)  +
         labs(x = "Contrast: Party autonomy TRUE - FALSE",
              y = "",
              fill = "Party autonomy")

# combine plots using the `patchwork` package
plot1 + plot2

Region-specific predictions and contrasts

Predicted media index by region and level of civil liberties:

pred <- predictions(mod,
                    newdata = datagrid(region = vdem_2015$region,
                                       party_autonomy = FALSE, 
                                       civil_liberties = seq(0, 1, length.out = 100))) |> 
        posteriordraws()

ggplot(pred, aes(x = civil_liberties, y = draw)) +
    stat_lineribbon() +
    scale_fill_brewer(palette = "Reds") +
    facet_wrap(~ region) +
    labs(x = "Civil liberties",
         y = "Media index (predicted)",
         fill = "")

Predicted media index by region and level of civil liberties:

pred <- predictions(mod,
                    newdata = datagrid(region = vdem_2015$region,
                                       civil_liberties = c(.2, .8),
                                      party_autonomy = FALSE)) |>
        posteriordraws()

ggplot(pred, aes(x = draw, fill = factor(civil_liberties))) +
    stat_halfeye(slab_alpha = .5) +
    facet_wrap(~ region) +
    labs(x = "Media index (predicted)",
         y = "Posterior density",
         fill = "Civil liberties")

Predicted media index by region and party autonomy:

pred <- predictions(mod,
                    newdata = datagrid(region = vdem_2015$region,
                                       party_autonomy = c(TRUE, FALSE),
                                       civil_liberties = .5)) |>
        posteriordraws()

ggplot(pred, aes(x = draw, y = region , fill = party_autonomy)) +
    stat_halfeye(slab_alpha = .5) +
    labs(x = "Media index (predicted)",
         y = "",
         fill = "Party autonomy")

TRUE/FALSE contrasts (marginal effects) of party autonomy by region:

mfx <- marginaleffects(mod,
                       variables = "party_autonomy",
                       newdata = datagrid(region = vdem_2015$region,
                                          civil_liberties = .5)) |>
        posteriordraws()

ggplot(mfx, aes(x = draw, y = region , fill = party_autonomy)) +
    stat_halfeye(slab_alpha = .5) +
    labs(x = "Media index (predicted)",
         y = "",
         fill = "Party autonomy")

Hypothetical groups

We can also obtain predictions or marginal effects for a hypothetical group instead of one of the observed regions. To achieve this, we create a dataset with NA in the region column. Then we call the marginaleffects or predictions functions with the allow_new_levels argument. This argument is pushed through via the ellipsis (...) to the posterior_epred function of the brms package:

dat <- data.frame(civil_liberties = .5,
                  party_autonomy = FALSE,
                  region = "New Region")

mfx <- marginaleffects(
    mod,
    variables = "party_autonomy",
    allow_new_levels = TRUE,
    newdata = dat)

draws <- posteriordraws(mfx)

ggplot(draws, aes(x = draw)) +
     stat_halfeye() +
     labs(x = "Marginal effect of party autonomy in a generic world region", y = "")

Multinomial logit

Fit a model with categorical outcome (heating system choice in California houses) and logit link:

dat <- "https://vincentarelbundock.github.io/Rdatasets/csv/Ecdat/Heating.csv"
dat <- read.csv(dat)

mod <- brm(depvar ~ ic.gc + oc.gc,
           data = dat,
           family = categorical(link = "logit"))

Adjusted predictions

Compute predicted probabilities for each level of the outcome variable:

pred <- predictions(mod)

head(pred)
#>   rowid     type group  predicted depvar  ic.gc  oc.gc   conf.low  conf.high
#> 1     1 response    ec 0.06612210     gc 866.00 199.69 0.04523796 0.09285630
#> 2     2 response    ec 0.07662654     gc 727.93 168.66 0.05935980 0.09695905
#> 3     3 response    ec 0.10356685     gc 599.48 165.58 0.06230961 0.16012522
#> 4     4 response    ec 0.06324271     er 835.17 180.88 0.04727269 0.08334918
#> 5     5 response    ec 0.07432869     er 755.59 174.91 0.05827615 0.09383320
#> 6     6 response    ec 0.07098633     gc 666.11 135.67 0.04553055 0.10322703

Extract posterior draws and plot them:

draws <- posteriordraws(pred)

ggplot(draws, aes(x = draw, fill = group)) +
    geom_density(alpha = .2, color = "white") +
    labs(x = "Predicted probability",
         y = "Density",
         fill = "Heating system")

Use the plot_cap function to plot conditional adjusted predictions for each level of the outcome variable gear, conditional on the value of the mpg regressor:

plot_cap(mod, condition = "oc.gc") +
    facet_wrap(~ group) +
    labs(y = "Predicted probability")

Marginal effects

mfx <- marginaleffects(mod)
summary(mfx)
#> Average marginal effects 
#>    Group  Term     Effect      2.5 %     97.5 %
#> 1     ec ic.gc -1.815e-04 -4.015e-04  2.051e-05
#> 2     ec oc.gc  4.855e-04 -4.396e-04  1.433e-03
#> 3     er ic.gc  1.542e-05 -2.145e-04  2.433e-04
#> 4     er oc.gc -1.023e-03 -2.054e-03 -3.274e-06
#> 5     gc ic.gc  1.281e-05 -3.750e-04  4.005e-04
#> 6     gc oc.gc  1.060e-03 -7.152e-04  2.750e-03
#> 7     gr ic.gc  3.973e-05 -2.439e-04  3.273e-04
#> 8     gr oc.gc  9.248e-05 -1.178e-03  1.406e-03
#> 9     hp ic.gc  1.135e-04 -6.839e-05  2.960e-04
#> 10    hp oc.gc -6.157e-04 -1.449e-03  1.745e-04
#> 
#> Model type:  brmsfit 
#> Prediction type:  response

Hurdle models

This section replicates some analyses from yet another amazing blog post by Andrew Heiss.

To begin, we estimate a hurdle model in brms with random effects, using data from the gapminder package:

library(gapminder)
library(brms)
library(dplyr)
library(ggplot2)
library(ggdist)
library(cmdstanr)
library(patchwork)
library(marginaleffects)

set.seed(1024)

CHAINS <- 4
ITER <- 2000
WARMUP <- 1000
BAYES_SEED <- 1234

gapminder <- gapminder::gapminder |> 
  filter(continent != "Oceania") |> 
  # Make a bunch of GDP values 0
  mutate(prob_zero = ifelse(lifeExp < 50, 0.3, 0.02),
         will_be_zero = rbinom(n(), 1, prob = prob_zero),
         gdpPercap = ifelse(will_be_zero, 0, gdpPercap)) |> 
  select(-prob_zero, -will_be_zero) |> 
  # Make a logged version of GDP per capita
  mutate(log_gdpPercap = log1p(gdpPercap)) |> 
  mutate(is_zero = gdpPercap == 0)

mod <- brm(
  bf(gdpPercap ~ lifeExp + year + (1 + lifeExp + year | continent),
     hu ~ lifeExp),
  data = gapminder,
  backend = "cmdstanr",
  family = hurdle_lognormal(),
  cores = 2,
  chains = CHAINS, iter = ITER, warmup = WARMUP, seed = BAYES_SEED,
  silent = 2)

Adjusted predictions

Adjusted predictions for every observation in the original data:

predictions(mod) %>% head()
#>   rowid     type predicted gdpPercap lifeExp year continent conf.low conf.high
#> 1     1 response  142.5456  779.4453  28.801 1952      Asia 103.1327  218.8240
#> 2     2 response  168.2657  820.8530  30.332 1957      Asia 124.8506  255.8359
#> 3     3 response  201.7414  853.1007  31.997 1962      Asia 152.9039  303.6894
#> 4     4 response  251.4996  836.1971  34.020 1967      Asia 196.5698  372.8495
#> 5     5 response  312.2482    0.0000  36.088 1972      Asia 249.5193  454.2419
#> 6     6 response  397.5390  786.1134  38.438 1977      Asia 324.5934  566.6799

Adjusted predictions for the hu parameter:

predictions(mod, dpar = "hu") %>% head()
#>   rowid     type predicted gdpPercap lifeExp year continent  conf.low conf.high
#> 1     1 response 0.5739495  779.4453  28.801 1952      Asia 0.4746831 0.6515670
#> 2     2 response 0.5365013  820.8530  30.332 1957      Asia 0.4416000 0.6112908
#> 3     3 response 0.4955596  853.1007  31.997 1962      Asia 0.4069162 0.5664250
#> 4     4 response 0.4455042  836.1971  34.020 1967      Asia 0.3656259 0.5106086
#> 5     5 response 0.3957129    0.0000  36.088 1972      Asia 0.3251678 0.4537019
#> 6     6 response 0.3413458  786.1134  38.438 1977      Asia 0.2823879 0.3907218

Predictions on a different scale:

predictions(mod, type = "link", dpar = "hu") %>% head()
#>   rowid type   predicted gdpPercap lifeExp year continent   conf.low
#> 1     1 link  0.29798370  779.4453  28.801 1952      Asia -0.1013542
#> 2     2 link  0.14626541  820.8530  30.332 1957      Asia -0.2346709
#> 3     3 link -0.01776198  853.1007  31.997 1962      Asia -0.3767284
#> 4     4 link -0.21885246  836.1971  34.020 1967      Asia -0.5510282
#> 5     5 link -0.42336036    0.0000  36.088 1972      Asia -0.7301227
#> 6     6 link -0.65730248  786.1134  38.438 1977      Asia -0.9326474
#>     conf.high
#> 1  0.62593415
#> 2  0.45274118
#> 3  0.26727993
#> 4  0.04244063
#> 5 -0.18572429
#> 6 -0.44427921

Plot adjusted predictions as a function of lifeExp:

plot_cap(
    mod,
    condition = "lifeExp") +
    labs(y = "mu") +
plot_cap(
    mod,
    dpar = "hu",
    condition = "lifeExp") +
    labs(y = "hu")

Predictions with more than one condition and the re_formula argument from brms:

plot_cap(
    mod,
    re_formula = NULL,
    condition = c("lifeExp", "continent"))

Extract draws with `posteriordraws()`

The posteriordraws() function extract raw samples from the posterior from objects produced by marginaleffects. This allows us to use richer geoms and summaries, such as those in the ggdist package:

predictions(
    mod,
    re_formula = NULL,
    newdata = datagrid(model = mod,
                       continent = gapminder$continent,
                       year = c(1952, 2007),
                       lifeExp = seq(30, 80, 1))) %>%
    posteriordraws() %>%
    ggplot(aes(lifeExp, draw, fill = continent, color = continent)) +
    stat_lineribbon(alpha = .25) +
    facet_grid(year ~ continent)

Average Contrasts

What happens to gdpPercap when lifeExp increases by one?

comparisons(mod) %>% summary()
#> Average contrasts 
#>                 Group    Term    Contrast Effect  2.5 % 97.5 %
#> 1 main_marginaleffect lifeExp (x + 1) - x 689.91 515.58 811.96
#> 2 main_marginaleffect    year (x + 1) - x -63.72 -84.44 -41.05
#> 
#> Model type:  brmsfit 
#> Prediction type:  response

What happens to gdpPercap when lifeExp increases by one standard deviation?

comparisons(mod, variables = list(lifeExp = "sd")) %>% summary()
#> Average contrasts 
#>                 Group    Term                Contrast Effect 2.5 % 97.5 %
#> 1 main_marginaleffect lifeExp (x + sd/2) - (x - sd/2)   4114  3718   4741
#> 
#> Model type:  brmsfit 
#> Prediction type:  response

What happens to gdpPercap when lifeExp increases from 50 to 60 and year simultaneously increases its min to its max?

comparisons(
    mod,
    variables = list(lifeExp = c(50, 60), year = "minmax"),
    interaction = TRUE) %>%
    summary()
#> Average contrasts 
#>                 Group lifeExp      year Effect 2.5 % 97.5 %
#> 1 main_marginaleffect 60 - 50 Max - Min  874.1 523.1   1404
#> 
#> Model type:  brmsfit 
#> Prediction type:  response

Plot draws from the posterior distribution of average contrasts (not the same thing as draws from the posterior distribution of contrasts):

comparisons(mod) %>%
    summary() %>%
    posteriordraws() %>%
    ggplot(aes(estimate, term)) +
    stat_dotsinterval() +
    labs(x = "Posterior distribution of average contrasts", y = "")

Marginal effects (slopes)

Average Marginal Effect of lifeExp on different scales and for different parameters:

marginaleffects(mod) |> summary()
#> Average marginal effects 
#>      Term Effect  2.5 % 97.5 %
#> 1 lifeExp 689.56 515.45 811.41
#> 2    year -63.72 -84.44 -41.05
#> 
#> Model type:  brmsfit 
#> Prediction type:  response

marginaleffects(mod, type = "link") |> summary()
#> Average marginal effects 
#>      Term    Effect    2.5 %    97.5 %
#> 1 lifeExp  0.082019  0.07419  0.088564
#> 2    year -0.009322 -0.01201 -0.006316
#> 
#> Model type:  brmsfit 
#> Prediction type:  link

marginaleffects(mod, dpar = "hu") |> summary()
#> Average marginal effects 
#>      Term    Effect     2.5 %    97.5 %
#> 1 lifeExp -0.008128 -0.009369 -0.006688
#> 
#> Model type:  brmsfit 
#> Prediction type:  response

marginaleffects(mod, dpar = "hu", type = "link") |> summary()
#> Average marginal effects 
#>      Term   Effect   2.5 %   97.5 %
#> 1 lifeExp -0.09909 -0.1132 -0.08382
#> 
#> Model type:  brmsfit 
#> Prediction type:  link

Plot Conditional Marginal Effects

plot_cme(
    mod,
    effect = "lifeExp",
    condition = "lifeExp") +
    labs(y = "mu") +

plot_cme(
    mod,
    dpar = "hu",
    effect = "lifeExp",
    condition = "lifeExp") +
    labs(y = "hu")

Or we can call marginaleffects() or comparisons() with posteriordraws() function to have even more control:

comparisons(
    mod,
    type = "link",
    variables = "lifeExp",
    newdata = datagrid(lifeExp = c(40, 70), continent = gapminder$continent)) %>%
    posteriordraws() %>%
    ggplot(aes(draw, continent, fill = continent)) +
    stat_dotsinterval() +
    facet_grid(lifeExp ~ .) +
    labs(x = "Effect of a 1 unit change in Life Expectancy")

Logistic regression with multiplicative interactions

Adjusted predictions

Marginal effects

Random effects model

Posterior predictions

Marginal effects and contrasts

Continuous predictors

Global grand mean

Region-specific predictions and contrasts

Hypothetical groups

Multinomial logit

Adjusted predictions

Marginal effects

Hurdle models

Adjusted predictions

Extract draws with posteriordraws()

Average Contrasts

Marginal effects (slopes)

Extract draws with `posteriordraws()`