I’m trying to fit a simple linear mixed model with random intercepts and no slopes, and understand all priors that I assume, so that I will be able to provide informative priors if I want to. The larger context is to build a workflow for an R&D department; understanding the priors is important here, as is the ability to use informative priors.
I’m stuck with the prior for the standard deviation of the random intercept,
prior_covariance = decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
I read http://mc-stan.org/rstanarm/articles/glmer.html and many related posts on this forum (@wds15 and @bgoodri gave several useful comments). From there I understand that the standard deviation is the tau in the link, and that the distribution of tau follows a Gamma distribution with the shape and scale prior as in the arguments of decov.
However, when I checked the MCMC results of the prior using the prior_PD=TRUE argument (so that posterior=prior by ignoring the data), I found that the mean and SD of the random intercepts are much larger: mean=16.66 and sd=69.6 instead of the mean=1 and sd=1 from a gamma(shape=1, scale=1) distribution.
I can link all the other priors understand to MCMC results for the other priors (intercept, and prior_aux) as well as for regression coefficients (these last are not illustrated in the example below).
Below is the code to simulate a small dataset with a reproduction of the phenomenon. I left out any predicors in this example although there are a couple in my main problem. In this example I have 10 groups of 100 observations, so the model has a global intercept, random group terms adding to the intercept, and a residual noise term.
I hope someone can enlighten me why the MCMC results for the “Sigma[group:(Intercept),(Intercept)]” are so different from what I expect from the documentation.
Many thanks in advance.
library(rstanarm)
set.seed(37)
df <- data.frame(y=rnorm(1000), group=(0:99) %/% 10)
k <- 1 # shape
theta <- 1 # scale
mod <- stan_lmer(y ~ 1 + (1 | group), df, seed=37,
prior_intercept = normal(0, 3),
prior_aux=normal(0, 3),
prior_covariance = decov(regularization = 1, concentration = 1,
shape=k, scale=theta),
prior_PD=TRUE)
prior_summary(mod)
summary(mod, pars=c("(Intercept)","sigma", "Sigma[group:(Intercept),(Intercept)]"), digits=5)
# https://en.wikipedia.org/wiki/Gamma_distribution
paste("group sigma expected, mean=", k*theta, " sd=", sqrt(k)*theta)
paste("group sigma expected 10,50,90 percentiles: ")
qgamma(c(0.1, 0.5, 0.9), shape=k, scale=theta)
(copy of the last part of the output)
> summary(mod, pars=c("(Intercept)","sigma", "Sigma[group:(Intercept),(Intercept)]"), digits=5)
Model Info:
function: stan_lmer
family: gaussian [identity]
formula: y ~ 1 + (1 | group)
algorithm: sampling
sample: 4000 (posterior sample size)
priors: see help('prior_summary')
observations: 1000
groups: group (10)
Estimates:
mean sd 10% 50% 90%
(Intercept) -0.01731 2.93088 -3.68178 -0.02156 3.76057
sigma 2.33549 1.82016 0.33168 1.97219 4.79182
Sigma[group:(Intercept),(Intercept)] 16.66615 69.61850 0.00710 1.21325 37.12862
MCMC diagnostics
mcse Rhat n_eff
(Intercept) 0.03338 0.99965 7710
sigma 0.02612 0.99992 4857
Sigma[group:(Intercept),(Intercept)] 1.33179 1.00013 2733
For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
> # https://en.wikipedia.org/wiki/Gamma_distribution
> paste("group sigma expected, mean=", k*theta, " sd=", sqrt(k)*theta)
[1] "group sigma expected, mean= 1 sd= 1"
> paste("group sigma expected 10,50,90 percentiles: ")
[1] "group sigma expected 10,50,90 percentiles: "
> qgamma(c(0.1, 0.5, 0.9), shape=k, scale=theta)
[1] 0.1053605 0.6931472 2.3025851
R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22621)
rstanarm_2.21.3