# Decov prior just for a random intercept is supposed to be gamma but is empirically different?

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.

``````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)
 "group sigma expected, mean= 1  sd= 1"
> paste("group sigma expected 10,50,90 percentiles: ")
 "group sigma expected 10,50,90 percentiles: "
> qgamma(c(0.1, 0.5, 0.9), shape=k, scale=theta)
 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