I am running a multilevel model in rstanarm (stan_lmer), and the model fails at the end of fitting with this error code:
Error in object@.MISC$summary$msd[col_idx, 1, drop = FALSE] :
subscript out of bounds
I can’t submit a replicable example as the data is not sharable, but if anyone can point me where in the package code I should look, I can take a look at debugging. :)
Updated rstanarm (and dependencies, R to the latest versions as an attempt to solve. Data looks fine by visual inspection (everything is the correct type, no NAs etc.)
Sure!
stan_lmer(formula = r ~ I + m + b +(1|sid),data = dat_cl)
r is continuous, l is numeric but constrained to be 0 or 1, m and b are numeric but constrained to be -.5 or .5 and sid is a factor. It runs fine in lmer, and when I tried simulating data with the same structure, it also ran fine.
Can you point me to where in rstanarm you calculate the summary so I can debug?
It could come from a few different places. You can look at traceback() after the error to see what sequence of functions was being called. You might have levels of sid that have spaces or other weird characters or else maybe some levels that are empty or something like that.
The traceback wasn’t that useful. I don’t think it’s just the levels of sid - if I run it with just the fixed effects or just the varying effects, it’s fine, the combination of the two runs the error.
It runs fine with brms and lmer, so it seems like this is particularly an rstanarm issue. I have a solution (I’ll just run brms). Is it helpful to keep looking into this?
It would be good if we could figure it out. This is saying that col_idx contains values that extend beyond the range of msd. That is why I thought it could be an off-by-one thing when counting up the levels, but perhaps it is due to something else.
Reproducible example (thanks to @jonah!). It’s labelling a fixed effect as “b” that is the issue. Changing the name from “b” to something else fixes the issue. Do you want a bug report on rstanarm?
d <- data.frame(x = rnorm(100), y = rnorm(100), b = rnorm(100), grp = gl(5, 20))
stan_lmer(y ~ x + b + (1|grp), data = d)