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

Error when generating marginal effects from brms model of the form y ~ 0 + intercept + x #113

Open
Kriaese opened this issue Jun 10, 2019 · 2 comments

Comments

Projects
None yet
2 participants
@Kriaese
Copy link

commented Jun 10, 2019

Dear Russell,

I am trying to obtain marginal effects using emmeans and tidybayes from a brms model with a prior on the real intercept and thus a model of the form y ~ 0 + intercept + x. However, this results in the following error:

Error in vcov(object)[nm, nm, drop = FALSE] : subscript out of bounds
Called from: emm_basis.brmsfit(object, trms, xlev, grid, ...)

Here’s a quick example (with mostly unchanged defaults) using the warpbreaks data set.

library(emmeans)
library(tidybayes)
library(magrittr)
library(brms)

set.seed(123)

# All seems to works here:
Model1 <- brm(breaks ~ tension,
              data = warpbreaks)

Model1 %>%
  emmeans(~ tension) %>%
  gather_emmeans_draws() %>%
  mean_qi()

# A tibble: 3 x 7
# tension .value .lower .upper .width .point .interval
# <fct>    <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
# L         36.4   30.9   42.0   0.95 mean   qi       
# M         26.2   20.6   31.7   0.95 mean   qi       
# H         21.6   16.1   27.2   0.95 mean   qi 
 
# Error occurs here:
Model2 <- brm(breaks ~ 0 + intercept + tension,
              data = warpbreaks)

Model2 %>%
  emmeans(~ tension) %>%
  gather_emmeans_draws() %>%
  mean_qi()

#Error in vcov(object)[nm, nm, drop = FALSE] : subscript out of bounds
#Called from: emm_basis.brmsfit(object, trms, xlev, grid, ...)

It seems the error has to do with how the internally generated variable intercept is handled, which might require special coding (similar to #104 ).

@rvlenth

This comment has been minimized.

Copy link
Owner

commented Jun 11, 2019

Actually, it is an oddity in the names assigned by brm vis a vis those used later:

# In the code:
    nm = eval(parse(text = "brms:::rename(colnames(X))"))
    V = vcov(object)[nm, nm, drop = FALSE]

# what we get:
Browse[2]> nm
[1] "intercept" "tensionL"  "tensionM"  "tensionH" 

Browse[2]> vcov(object)
          intercept  tensionM  tensionH
intercept  7.954110 -7.931893 -7.787033
tensionM  -7.931893 15.631325  7.807514
tensionH  -7.787033  7.807514 15.885689

So you see that I ask the package to tell me the names and to return the variance-covariance matrix, but the package is inconsistent with itself. That's kinda annoying.

I think I should probably just use the empirical cavariance of the posterior sample anyway, since we realkly don't need this result at all except in a frequestist summary.

I am traveling as I write and will not be able to do anything more with this until late June 2019.

@Kriaese

This comment has been minimized.

Copy link
Author

commented Jun 12, 2019

Aha I see - many thanks for your prompt reply! Looking forward to your fix and happy traveling!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.