library(brms)
library(data.table)
library(marginaleffects)

dat <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/Stat2Data/Titanic.csv")
dat$PClass <- factor(dat$PClass)
dat$Sex <- factor(dat$Sex)
mod <- brm(Survived ~ Age + Sex, data = dat, family = bernoulli)
# marginal effects & contrasts
mfx <- marginaleffects(mod)

# extract posterior draws into a data.table
bb <- posteriordraws(mfx)
setDT(bb)

# draw weights
rudirichlet <- function(N) {
    out <- rexp(N)
    out / sum(out)
}
bb[, w := rudirichlet(.N)]

# average marginal effects
bb <- bb[, .(avg_bb = weighted.mean(draw, w)),
           by = c("term", "contrast", "drawid")]

# summary
bb[, as.list(quantile(avg_bb, c(.5, .025, .975))),
     by = c("term", "contrast")]
##    term      contrast          50%         2.5%         97.5%
## 1:  Age         dY/dX -0.001073474 -0.003213676  0.0009279447
## 2:  Sex male - female -0.545990444 -0.608615602 -0.4828708461
# fixed-X version for reference
summary(mfx)[, c("term", "contrast", "estimate", "conf.low", "conf.high")]
## Average marginal effects 
##   Term      Contrast    Effect    CI low    CI high
## 1  Age         dY/dX -0.001096 -0.003205  0.0009273
## 2  Sex male - female -0.545944 -0.608612 -0.4830281
## 
## Model type:  
## Prediction type:

GLM version:

mod <- glm(Survived ~ Age + Sex, data = dat, family = binomial)
mfx <- marginaleffects(mod)

# 1000 replicates of the unit-level marginal effects
bb <- lapply(1:1000, \(i) data.table(mfx, B = i))
bb <- rbindlist(bb)

# draws weights
bb[, w := rudirichlet(.N)]

# average marginal effects
bb <- bb[, .(avg_bb = weighted.mean(dydx, w)), by = c("term", "contrast", "B")]

# quantiles
bb[, as.list(quantile(avg_bb, c(.5, .025, .975))), by = c("term", "contrast")]
##    term      contrast          50%         2.5%        97.5%
## 1:  Age         dY/dX -0.001089035 -0.001095584 -0.001082956
## 2:  Sex male - female -0.546538370 -0.546708581 -0.546343821
# Delta method for comparison
summary(mfx)[, c("term", "contrast", "estimate", "conf.low", "conf.high")]
## Average marginal effects 
##   Term      Contrast    Effect    CI low    CI high
## 1  Age         dY/dX -0.001089 -0.003166  0.0009878
## 2  Sex male - female -0.546532 -0.608527 -0.4845368
## 
## Model type:  
## Prediction type: