Can I use mediation() with ordinal mediator using brms mo()?

[JonRBrauer]

I am estimating brms models with a binary outcome (y) and an ordinal mediator (m) specified as an ordered factor m with cumulative logit (or probit) link in 'm' model and as a monotonic predictor with mo(m) in 'y' model. I want to estimate (conditional) direct and indirect effects of x on y via m using bayestestR::mediation() function.

I could not tell from the documentation whether such models are supported in this package, though the similar mediate package appears to handle such models. In this case, mediation() is generating errors (e.g., "undefined columns selected"). Specifically, the error message is:

Error in [.data.frame(insight::get_parameters(model), c(coef_indirect, :
undefined columns selected

If the mediation function can/should generalize to the posterior estimates from these models, then perhaps the errors are due to brms output attaching "mo" in front of parameter names in the output?

I reproduced this error in a modified version of the vignette example, where 'job_seekf' (m) is an ordinal ordered factor and 'depress2bin' (y) is a binary factor. As you will see below, I attempted to explicitly specify the treatment, mediator, and response variables in an attempt to overcome the error. I also tried various combinations of "mojob_seekf" as mediator and/or response to no avail. Thank you for your time and any info you can provide! : )

Jon

library(bayestestR)
library(brms)

load sample data

data(jobs)

jobs <- jobs %>%
mutate(
job_seekf = if_else(job_seek >= 1 & job_seek < 2, 1, job_seek),
job_seekf = if_else(job_seek >= 2 & job_seek < 3, 2, job_seekf),
job_seekf = if_else(job_seek >= 3 & job_seek < 4, 3, job_seekf),
job_seekf = if_else(job_seek >= 4 & job_seek < 5, 4, job_seekf),
job_seekf = if_else(job_seek == 5, 5, job_seekf),
job_seekf = factor(job_seekf, ordered=TRUE, levels = c(1, 2, 3, 4, 5)),
# educf = factor(educ, ordered=TRUE, levels = c("lt-hs", "highsc", "somcol", "bach", "gradwk")),
depress2bin = factor(if_else(depress2 <2, "0", "1"))
)

Fit Bayesian mediation model in brms

f1 <- bf(job_seekf ~ 1 + treat,
family = cumulative("logit"))
f2 <- bf(depress2bin ~ 1 + treat + mo(job_seekf),
family = "bernoulli")
m2 <- brm(f1 + f2 + set_rescor(FALSE), data = jobs,
cores = 4, chains = 4, backend = "cmdstanr", seed=1138)

Model summary

summary(m2)

Mediation estimates

medm2 <- mediation(m2,
treatment="treat",
mediator="job_seekf",
response=c(m="job_seekf", y="depressbin2"),
ci=.95)
medm2