I have run a heirarchical model using stan_lmer but despite increasing adapt_delta up to 0.999 I continue to get the same one divergent draw… As suggested, I ran pairs, etc. and indeed I found a funnel. So I read about how to reparameterize stan code (including the famous 8 schools problem). But I would rather not have to figure out how to submit my problem in native STAN code. (After all, that is why I have been using rstanarm.) I also discovered how to recover the stan code from my model, and as you folks have been saying, that model specification is not intended for human reading. One suggestion was to use brms, but that code was only marginally more comprehensible to me.
So I have two questions:
a) Is it possible to obtain the fully compiled stan code (with expansion of all the headers, etc) to look at what is actually sent by rstanarm to STAN? I might be able to work with that.
b) Given that the funnel is a common enough issue with hierarchical models, has anyone considered giving the user the option of submitting a non-centered version of the model? I recognize that for a big enough sample, using the centered code is more efficient. But where one hits a funnel and reparameterization is the only available fix, it would be nice to be able to run a model without divergences.
I would submit the model and the data, but that is confidential.
rstanarm Version: 2.21.1 on R 4.1.2, RStudio 2021.09.0 Build 351 Windows 10 up-t-date.
Thanks for any suggestions.
Larry Hunsicker
a) Yes, you can extract the raw Stan code with rstan::get_stancode(mod$stanfit) where mod is your model fit with rstanarm.
b) I don’t know about non-centering in rstanarm, but you should be able to achieve this with little additional effort in brms. The formula syntax is very similar, so stan_glm(y ~ x + z, family = poisson(), data = mydata)
becomes brm(y ~ x + z, family = poisson(), data = mydata).
To fit the model with a non-centered parameterization, just wrap the formula in the helper function bf and specify the center = FALSE: brm(bf(y ~ x + z, center = FALSE), family = poisson(), data = mydata).
In short, you don’t have to dig into the raw Stan code made by brms to adjust the parameterization strategy. You’ll pay a small cost to compile the model each time you fit, but in most cases I’ve encountered this is on the order of seconds to a minute. The payout is much greater flexibility to fit different models that rstanarm can’t (because it relies on pre-compiled models).
Thanks all! Interesting that rstanarm routinely uses a non-centered model. I have not used brms because rstanarm is so easy to use, and generally works extremely well. But I’ll try brms with and without the code to use a non-centered model. In fact, I think that the assertion that the non-centered model is equivalent to the centered one is not quite right. The non-centered models actually permit the variances, as well as the intercepts, to vary among the schools (or whatever the lower level is). And I think that it is more realistic to assume that the “real” variances may vary meaningfully. If the estimated sd of the variance multiplier is within a credible interval, one could then use a centered approach and save a df. If the estimated variances of the lower units actually are all pretty close to uniform, they shouldn’t create negative correlations between the lower units that seem to cause the funnel. In my case, I suspect that my funnel results from some real correlations among the top level covariates. That might require a different sort of reparameterization to fix.
I wish a Happy New Year to all of you folks! Certainly one better than the last two years!