Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Multivariate response models #31

Open
strengejacke opened this issue Mar 11, 2019 · 9 comments · May be fixed by #687
Open

Multivariate response models #31

strengejacke opened this issue Mar 11, 2019 · 9 comments · May be fixed by #687
Assignees
Labels
Help us 👀 Help is needed to implement something

Comments

@strengejacke
Copy link
Member

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.

@DominiqueMakowski DominiqueMakowski added this to the 0.2.0 milestone Apr 2, 2019
@DominiqueMakowski
Copy link
Member

DominiqueMakowski commented May 12, 2019

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

@strengejacke
Copy link
Member Author

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...

@DominiqueMakowski DominiqueMakowski modified the milestones: 0.2.0, 0.3.0 May 19, 2019
@DominiqueMakowski DominiqueMakowski added the Help us 👀 Help is needed to implement something label Jun 21, 2019
@strengejacke
Copy link
Member Author

strengejacke commented Sep 13, 2019

Unbenannt

This suggests some despair of this issue... 😆

@DominiqueMakowski
Copy link
Member

haha 😂

@bwiernik
Copy link
Contributor

bwiernik commented Aug 9, 2021

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?

@mattansb
Copy link
Member

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

@jebyrnes
Copy link

Hey! Futzing with this with rope() on a bayesian SEM fit with brms. I like the function, but, filtering to a specific parameter is a pain and then adding ranges for those parameters is a pain. Behavior is inconsistent.

rope(rich_brms, 
     parameters = "firesev_age",
     range = list(firesev_age = c(-0.017, 0.017)))

Gives me

Error: With a multivariate model, `range` should be
  'default' or a list of named numeric vectors with
  length 2.

Which - that's the actual coefficient name.......... Ugh.

So, I see what you're trying to do with the parameters argument using regex and, while I applaud the goal, I'm finding it really hard to use. I would be happier with just "supply a vector of paramters". This would be nice as then that vector of parameters could match what is in the range list name. And probably easier to check the names of what goes into range against the names in parameters, I'm guessing.

It would also be nice to be able to specify, even for many parameters, just a vector of 2 numbers to set a range for all tests instead of having to go with a named list.

If your purpose had been to have all betas and intercepts, perhaps in the component argument, you could add betas and intercepts?

Also, can plotting MV models just use facet_wrap for each coef?

@bwiernik
Copy link
Contributor

Can you give a reproducible example?

@strengejacke strengejacke self-assigned this Dec 16, 2024
strengejacke added a commit that referenced this issue Dec 16, 2024
@strengejacke strengejacke linked a pull request Dec 16, 2024 that will close this issue
@strengejacke
Copy link
Member Author

I think this can be solved by improving the docs. See this example, which I added to the docs:

library(bayestestR)
library(brms)
model <- brm(
  bf(mvbind(mpg, disp) ~ wt + cyl) + set_rescor(rescor = TRUE),
  data = mtcars
)

r <- list(
  mpg = list(b_mpg_wt = c(-1, 1), b_mpg_cyl = c(-2, 2)),
  disp = list(b_disp_wt = c(-5, 5), b_disp_cyl = c(-4, 4))
)
rope(model, range = r)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Help us 👀 Help is needed to implement something
Projects
None yet
Development

Successfully merging a pull request may close this issue.

6 participants