I’m attempting to do a multilevel mediation analysis. The dv is a likert scale response which I am modeling with a cumulative distribution (logit link). The mediator is categorical (Yes / No). When calculating the ab path, I am wondering if I can multiply the a and b coefficients directly (as below).

Counterfactual_mod <- bf(CounterfactualResponse ~ Condition + (1|Subject), family = “bernoulli”)

model1 <- brm(Knowledge_Judgment_mod +
Counterfactual_mod,
set_rescor(FALSE),
data=data1,
iter=3000,
warmup=1500,
control = list(adapt_delta = .90),
prior = c(set_prior(“normal(0,1)”, class = “b”)))

a <- posterior_samples(model1, pars = “CounterfactualResponse_ConditionSpecific”)
b <- posterior_samples(model1, pars = “KnowledgeResponse_CounterfactualResponse”)
a_b <- a*b

Maybe this is of interest for you: I have drafted a function, mediation(), that prints a compact summary especially for mediation analysis (similar to the output from the mediation package).

If the mediator is normally distributed, and the outcome is negative binomial distributed, does taking the product of the coefficients still work for calculating the indirect effects? Thanks!

Ah, gotcha, so in your original example, it should work because the logit link applied to both. Makes sense. Any chance you know of a conversion for calculating the indirect effects in my case? Thanks again.

Hi! Thanks a lot for providing this function! I have two questions related to it

In my mediation model, the mediator is normally distributed and the dependent variable can be considered cumulative.
bfy ~bf(y~x+m+(1+x+m|ID|Subject), family=“cumulative”)
bfm ~bf(m~x+(1+x|ID|Subject), family=“gaussian”)
As these families do not have the same link (cumulative having a logit-link) - can I still use the mediation()-function to calculate the effects (especially the indirect effect)? And if not, do you know a workaround for that?

Is it possible to include covariates into the model als illustrated below? (Please note that I am not including the interaction terms with the covariate as random effect as well as this does not add much to the predictive capacity of the model)
bfy ~bf(y~xcov1+mcov1+(1+x+m+cov1|ID|Subject), family=“cumulative”)
bfm ~bf(m~x+(1+x|ID|Subject), family=“gaussian”)
Thanks a lot for your help!