What are the rstanarm::stan_glmer Error Term standard deviations reporting?

Modeling
1

What exactly are the posterior summaries reported for the random effect variances in the print method for stan_glmer fits? When I manually compute point estimates from the draws, I get something pretty different than the output of the print method.

For example, using data from the Pilots dataset from chapter 13 of the ARM examples in the dataframe df, I run

rstan_fit <- stan_glmer(y ~ 1 + (1 | group_id) + (1 | scenario_id), df)
print(rstan_fit)

This displays (in part):

Error terms:
 Groups      Name        Std.Dev.
 scenario_id (Intercept) 0.348   
 group_id    (Intercept) 0.067   
 Residual                0.224

In contrast, we can look directly at the draws:

rstan_draws <- as.matrix(rstan_fit)
group_id_col <- "Sigma[group_id:(Intercept),(Intercept)]"
scenario_id_col <- "Sigma[scenario_id:(Intercept),(Intercept)]"

median(sqrt(rstan_draws[, scenario_id_col]))
mean(sqrt(rstan_draws[, scenario_id_col]))

median(sqrt(rstan_draws[, group_id_col]))
mean(sqrt(rstan_draws[, group_id_col]))

The output of which is:

> median(sqrt(rstan_draws[, scenario_id_col]))
[1] 0.3168345
> mean(sqrt(rstan_draws[, scenario_id_col]))
[1] 0.3335378

> median(sqrt(rstan_draws[, group_id_col]))
[1] 0.03530919
> mean(sqrt(rstan_draws[, group_id_col]))
[1] 0.04839807

None of these manually computed posterior summaries match the printed output. The group id standard deviation is particularly different.

I wondered if the parameters are being transformed? If so, what was the transformation? Or is a different posterior summary reported? Or something else?

1 Like
2

Hi, sorry for not getting to your question earlier, it is relevant and well written.

So what gets printed is the result of VarCorr(rstan_fit) which matches with:

rstan_draws <- as.matrix(rstan_fit)
group_id_col <- "Sigma[group_id:(Intercept),(Intercept)]"
scenario_id_col <- "Sigma[scenario_id:(Intercept),(Intercept)]"

sqrt(mean(rstan_draws[, scenario_id_col]))
sqrt(mean(rstan_draws[, group_id_col]))

I am not completely sure what is the reationale for first averaging and then taking the square root, I don’t think that’s a very good quantity. I think @jonah would know better. My buest guess is that this is simply what the packages rstanarm tries to mimic report and so this is for compatibility.

Best of luck with your model!