Rstanarm: ranef vs b parameters

Interfaces rstanarm
,
1

Question

What’s the difference between ranef and the b parameters?

According to the https://cran.r-project.org/web/packages/rstanarm/vignettes/glmer.html#glms-with-group-specific-terms “b as the random effects” and

Consequently, b is not a “parameter” to be estimated because parameters are unknown constants that > are fixed in repeated sampling.

I want to make sure I’m interpreting the summary output correctly.

Reproducible Example

I was trying to recreate a radon model found in Gelman and Hill’s text (pg. 282) but had a question about the posterior summary results.

The model in question has a varying intercept and slope with a single group-level predictor:

\begin{aligned} y_i &\sim N(\alpha_{j[i]} + \beta_{j[i]}x_i, \sigma_y^2), \text{ for } i=1,\ldots,n \\ \begin{pmatrix} \alpha_j \\ \beta_j \end{pmatrix} &\sim N \left( \begin{pmatrix} \gamma_0^\alpha + \gamma_1^\alpha u_j \\ \gamma_0^\beta + \gamma_1^\beta u_j \end{pmatrix}, \begin{pmatrix} \sigma_\alpha^2 & \rho \sigma_\alpha \sigma_\beta \\ \rho \sigma_\alpha\sigma_\beta & \sigma_\beta^2 \end{pmatrix}, \text{ for } j=1,\ldots,J \right) \end{aligned}

also written

y_i = \left[ \gamma_0^\alpha +\gamma_1^\alpha u_{j[i]} +\eta_{j[i]}^\alpha \right] + \left[ \gamma_0^\beta + \gamma_1^\beta u_{j[i]} +\eta_{j[i]}^\beta \right]x_i +\epsilon_i

If we look at county 85, the b values are not the same as the values produced by ranef. What are the b parameter values?

library(rstanarm)
demo(package="rstanarm",topic="ARM_Ch12_13")

tail(ranef(M4)$county,8)
tail(posterior_summary(M4)[,1:3],8)

which yields

:    (Intercept) x           
: 78  0.02177865 -0.009344675
: 79 -0.11318990 -0.087323466
: 80 -0.01616251 -0.089506652
: 81  0.10016236  0.156978232
: 82  0.01636097  0.004331113
: 83 -0.05364101 -0.188520895
: 84  0.06537049  0.019689862
: 85 -0.03264432 -0.004255915
:                                       Estimate     Est.Error  Q2.5        
: b[(Intercept) county:84]               0.071846727 0.11608273 -0.140401534
: b[x county:84]                         0.027510480 0.23361550 -0.446632534
: b[(Intercept) county:85]              -0.038354974 0.14350999 -0.347426522
: b[x county:85]                        -0.004709328 0.25448657 -0.523348587
: sigma                                  0.752818158 0.01826484  0.717829578
: Sigma[county:(Intercept),(Intercept)]  0.021750010 0.01406291  0.002438374
: Sigma[county:x,(Intercept)]            0.007244463 0.01625235 -0.028200743

I was able to reproduce the calculations for the varying intercept and slope using the ranef function as well which would seem to confirm that the b parameters are not the same as ranef and that ranef behavior aligns with that in the text:

# \alpha_{85}
(first_bracket <- fixef(M4)["(Intercept)"] +  fixef(M4)["u"]*unique(radon$u)[85] + ranef(M4)$county[85,c("(Intercept)")])

# \beta_{85}
(second_bracket <- fixef(M4)["x"] + fixef(M4)["x:u"] * unique(radon$u)[85] + ranef(M4)$county[85, c("x")])

or

a.hat.M4 <- coef(M4)$county[,1] + coef(M4)$county[,3]*unique(radon$u)
b.hat.M4 <- coef(M4)$county[,2] + coef(M4)$county[,4]*unique(radon$u)

# Just show for county 85
a.hat.M4[85]
b.hat.M4[85]

Aside

Aside: Is there a nicer way of retrieving the alpha and beta vectors in the presence of a group-level predictor so they can be plotted? I’d really like to do something like the output of mcmc_intervals.

I’d be grateful for any help you could provide. Thank you and thanks to the people working on this project.

1 Like
2

Hi Alex, welcome to the Stan community!

I’m sorry, I don’t really have much time to look into this. I just wanted to say that one reason why I switched to Stan (the language) is, that I really disliked the way group level predictors work in the lme4 syntax.

I don’t actually se what the difference is (I did not run the code), but could it be a difference in posterior point estimates used? Like, could it be that ranef returns the posterior mean, and coef returns the posterior median or something? I’m tagging @jonah , although he might be too busy to have a look.

For a quick and (kind of) dirty solution I think you could go with broom::tidy(). The better way would probably include the posterior and/or the tidybayes packages.

Cheers,
Max

1 Like