easystats / bayestestR

:ghost: Utilities for analyzing Bayesian models and posterior distributions
https://easystats.github.io/bayestestR/
GNU General Public License v3.0
563 stars 55 forks source link

Flexible ROPE values for describe_posterior #643

Closed rduesing closed 2 weeks ago

rduesing commented 7 months ago

Hello!

First of all, many thanks for the wonderful package!

I have a question and/or feature idea. I estimated models with brms, drew posterior samples and used "describe_posterior", which really worked fine. But is there a possibiliy to set different upper and lower ROPE bounds (or also different values for pd) for different variables of the posterior?

For example: I calculated a simple t-test like model with two groups and different standard deviations for each group. Within the posterior draws, I also calculated the effect sizes, so I have 5 variables in my posterior object. If I use

post |> 
describe_posterior(centrality = "all",
                   ci_method = "hdi",
                   test = c("p_direction", "rope"),
                   rope_range = c(80, 100))

for example, the upper and lower bounds of 100 and 80 will be used for all variables, which is not useful, for this application.

It would be very helpful, if I could enter a list instead of a vector, so that each variable could get it's own ROPE bounds. Maybe something like this:

post |> 
describe_posterior(centrality = "all",
                   ci_method = "hdi",
                   test = c("p_direction", "rope"),
                   rope_range = list(
                                      c(80, 100), # ROPE for group 1
                                      c(100, 120), # ROPE for group 2
                                      c(5, 20), # ROPE for SD 1
                                      c(10, 30), # ROPE for SD 2
                                      c(-0.1, 0.1) # ROPE for standardized effect size
                                       )

As far as I found, a list is only possible for multivariate analyses, where I can submit different ROPE bounds for each dependent variable, but these would be the same for every variable/predictor within each dependent variable. Is this correct?

Please correct me, if this option is already available and I missed it. Many thanks for help! Rainer

mattansb commented 6 months ago

Not currently supported, but you can do this manually by first extracting posteriors and then processing them...

library(rstanarm)
library(bayestestR)

model <- stan_glm(mpg ~ factor(cyl) + am, data = mtcars)

posteriors <- insight::get_parameters(model, component = "all")
names(posteriors)
#> [1] "(Intercept)"  "factor(cyl)6" "factor(cyl)8" "am"           "sigma"       

ropes <- list(
  c(-10, 10),
  c(-5, 5),
  c(-5, 5),
  c(-2, 2),
  c(0, 0.5)
)

mapply(
  function(p, r) {
    describe_posterior(p, centrality = "all",
                       ci_method = "hdi",
                       test = c("p_direction", "rope"),
                       rope_range = r)
  }, 
  p = posteriors,
  r = ropes, 
  SIMPLIFY = FALSE
) |> 
  do.call(what = rbind) |> 
  datawizard::data_modify(Parameter = names(posteriors))
#> Summary of Posterior Distribution
#> 
#> Parameter    | Median |   Mean |   MAP |          95% CI |     pd |            ROPE | % in ROPE
#> -----------------------------------------------------------------------------------------------
#> (Intercept)  |  24.76 |  24.78 | 24.78 | [ 21.89, 27.32] |   100% | [-10.00, 10.00] |        0%
#> factor(cyl)6 |  -6.10 |  -6.13 | -5.88 | [ -9.35, -3.15] | 99.98% | [ -5.00,  5.00] |    22.26%
#> factor(cyl)8 | -10.03 | -10.01 | -9.98 | [-12.96, -6.97] |   100% | [ -5.00,  5.00] |        0%
#> am           |   2.58 |   2.58 |  2.52 | [ -0.06,  5.33] | 96.88% | [ -2.00,  2.00] |    31.84%
#> sigma        |   3.14 |   3.19 |  3.01 | [  2.40,  4.07] |   100% | [  0.00,  0.50] |        0%