# 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(radonu) + ranef(M4)$county[85,c("(Intercept)")]) # \beta_{85} (second_bracket <- fixef(M4)["x"] + fixef(M4)["x:u"] * unique(radon$u) + 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
b.hat.M4


## 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

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