Obtaining Pd values from Mediation Analysis using Bayesian Regression Models

[JulianGaviriaL]

Hello easystats team.

I am trying to implement a Mediation Analysis using Bayesian Regression Models (brms package) proposed by Daniel Lüdecke:

https://strengejacke.github.io/sjstats/articles/bayesian-statistics.html

# Fit Bayesian mediation model

f1 <- bf(job_seek ~ treat + econ_hard + sex + age)
f2 <- bf(depress2 ~ treat + job_seek + econ_hard + sex + age)
m2 <- brm(f1 + f2 + set_rescor(FALSE), data = jobs, cores = 4)
mm <-mediation(m2)

I would like to use the "p_direction" function from "bayestestR" package in order to obtain the pd values for the summary table:

p_direction(mm)

Given the strange output:

parameter |      pd
-------------------
value     |  80.00%
hdi.low   | 100.00%
hdi.high  |  80.00%

I wonder whether it would be possible to obtain such pd values from this kind of models.

Thanks in advance,

Julian,

[strengejacke]

library(bayestestR)
m2 <- insight::download_model("brms_mv_6")

mm <- mediation(m2)
mm
#> # Causal Mediation Analysis for Stan Model
#> 
#>   Treatment: treat
#>   Mediator : job_seek
#>   Response : depress2
#> 
#> Effect                 | Estimate |          89% ETI
#> ----------------------------------------------------
#> Direct Effect (ADE)    |   -0.040 | [-0.110,  0.031]
#> Indirect Effect (ACME) |   -0.015 | [-0.036,  0.004]
#> Mediator Effect        |   -0.240 | [-0.285, -0.195]
#> Total Effect           |   -0.055 | [-0.129,  0.018]
#> 
#> Proportion mediated: 28.14% [-71.11%, 127.40%]

pd(as.data.frame(mm))
#> # Probability of Direction (pd)
#> 
#> Parameter           |      pd
#> -----------------------------
#> effect_direct       |  81.65%
#> effect_indirect     |  90.35%
#> effect_mediator     | 100.00%
#> effect_total        |  87.85%
#> proportion_mediated |  83.85%

^{Created on 2020-07-08 by the reprex package (v0.3.0)}

[strengejacke]

I added an as.data.frame() method that returns the posterior samples of the effects, which you can plug in in all baytestestR functions.

[JulianGaviriaL]

Hi Daniel,

Thank you very much for your quick feedback. I apologize in advance, in case I am doing something silly: 1) I closed R and updated all the packages 2) I ran your script:

 library(bayestestR)
 library(sjstats)
 m2 <- insight::download_model("brms_mv_6")

> mm <- mediation(m2)
Warning messages:
1: Argument 'exact_match' is deprecated. Please use 'fixed' instead. 
2: Argument 'exact_match' is deprecated. Please use 'fixed' instead. 
3: Argument 'exact_match' is deprecated. Please use 'fixed' instead. 
> mm
# Causal Mediation Analysis for Stan Model

  Treatment: treat
   Mediator: job_seek
   Response: depress2

                 Estimate    HDI (90%)
  Direct effect:    -0.04 [-0.11 0.03]
Indirect effect:    -0.02 [-0.04 0.00]
   Total effect:    -0.05 [-0.13 0.02]

Proportion mediated: 28.14% [-79.57% 135.86%]

But the output is still the same in my case:

> pd(as.data.frame(mm))
# Probability of Direction (pd)

Parameter |      pd
-------------------
value     |  80.00%
hdi.low   | 100.00%
hdi.high  |  80.00%

[strengejacke]

Don't load sjstats, mediation() is re-implemented in bayestestR.

[DominiqueMakowski]

Also, do not forget to update bayestestR by running remotes::install_github("easystats/bayestestR")

[JulianGaviriaL]

1) It worked when: -I did not load sjstats -I updated remotes::install_github("easystats/bayestestR") Thanks Daniel & Dominique

2) Now I have difficulties for having the C++ compiled brm model (m2 variable in the environment), when I run my own data. Therefore I can not apply the mediation function:

> library(bayestestR)
> library(brms)

> # Fit Bayesian mediation model
> f2 <- bf(NA_AFF ~ C3a)
> f3 <- bf(C4pa ~C3a+NA_AFF)
> m2 <- brm(f2 + f3 + set_rescor(FALSE), data = dat, cores = 4)
Compiling the C++ model
> 
> mm <-mediation(m2)
Error in mediation(m2) : object 'm2' not found

[strengejacke]

This rather looks like a brms/Stan/C++/R issue, not related to bayestestR... Do you have further information from the traceback?

[JulianGaviriaL]

No, unfortunately I do not have further information from the traceback. 1) But as I guess you realized, Is interesting that one can have the brms variable in the environment, when using your insight function. See your example:

library(bayestestR)
m2 <- insight::download_model("brms_mv_6")
mm <- mediation(m2)

2) I will ask in the Rstan forum then. And I'll let you know whether I get some feedback