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