bayestestR
bayestestR copied to clipboard
easystats
Multivariate response models
Maybe we manage it to make functions like equivalence_test() also work for multivariate response models (like stanmvreg), however, this seems more like a long-term goal to me.
I think in this case we could add a column "Response" and basically have the current dataframes (with indices for all parameters) rbinded with Response?
Response | Parameter | CI | ROPE
y1 | x1 |0.9 | 0.25
y1 | x1 |0.95| 0.13
y1 | x2 |0.9 | 0.122
y1 | x2 |0.95| 0.8
y2 | x1 |0.9 | 0.54
y2 | x1 |0.95| 0.8
y2 | x2 |0.9 | 0.456
y2 | x2 |0.95| 0.4
It's a bit more difficult. We first check the rope-range ( rope_range()), and then perform the equivalence-test. The good thing: insight::get_parameters() is enough for multivariate-response models, however, rope_range() returns an array of ranges. And this must be considered in the code...
It's a bit more difficult. We first check the rope-range ( rope_range()), and then perform the equivalence-test. The good thing: insight::get_parameters() is enough for multivariate-response models, however, rope_range() returns an array of ranges. And this must be considered in the code...
rope_range() now returns a named list, so I think this issue might be accommodated at this point?
I think we should (at some point?) add/allow multivariate statistics for multivariate models.
For example, instead of testing how much of the slope of X is in the ROPE of Y1, and also testing how much of the slope of X is in the ROPE of Y2, we can test how much of the slope(s) of X are in both the ROPE of Y1 and the ROPE of Y2 = what is the % of the posterior that indicates no multivariate effect at all:
library(brms)
library(bayestestR)
library(dplyr)
model <- brm(mvbind(mpg, wt) ~ cyl + disp, data = mtcars,
chains = 1, iter = 200)
MVROPE <- rope_range(model)
samps <- insight::get_parameters(model)
# ROPE
samps %>%
mutate(
across(starts_with("b_mpg"), ~between(.x, MVROPE$mpg[1], MVROPE$mpg[2])),
across(starts_with("b_wt"), ~between(.x, MVROPE$wt[1], MVROPE$wt[2])),
) %>%
summarise(
cyl = across(ends_with("cyl")) %>% apply(1, all) %>% mean(),
disp = across(ends_with("disp")) %>% apply(1, all) %>% mean()
) %>% t() %>% `colnames<-`("% in ROPE")
#> % in ROPE
#> cyl 0.06
#> disp 1.00
Similarly, we can also do this with pd: we can ask how much of the joint posterior of the effect(s) of X has a direction = what is the probability that the multivariate effect has a direction:
samps %>%
mutate(
across(starts_with("b_mpg"), ~sign(.x) == sign(median(.x)) & sign(median(.x)) != 0),
across(starts_with("b_wt"), ~sign(.x) == sign(median(.x)) & sign(median(.x)) != 0)
) %>%
summarise(
cyl = across(ends_with("cyl")) %>% apply(1, all) %>% mean(),
disp = across(ends_with("disp")) %>% apply(1, all) %>% mean()
) %>% t() %>% `colnames<-`("p(d)")
#> p(d)
#> cyl 0.52
#> disp 0.96
