Improper order of operations for VarCorr

Summary:

VarCorr uses an improper order of operations (a la Jensen's inequality) causing biased estimates.

Description:

VarCorr applies transformations to parameter estimates (variance and covariance parameters) rather than posterior draws and then estimating the parameters of interest (standard deviation and correlation parameters). This causes there to be a biased estimate for the parameters of interest.

Reproducible Steps:

library(rstanarm)
data(Gcsemv, package = "mlmRev")

seed = 1775
set.seed(seed)
M3_stanlmer <- stan_lmer(formula = course ~ gender + (1 + gender | school), 
                         data = Gcsemv[is.element(Gcsemv[,"school"],sample(unique(Gcsemv[,"school"],10))),],
                         seed = seed)

pdraws = as.data.frame(M3_stanlmer)

# Output
VarCorr(M3_stanlmer)

# Unbiased approach
mean(sqrt(pdraws[,"Sigma[school:(Intercept),(Intercept)]"]))
mean(sqrt(pdraws[,"Sigma[school:genderM,genderM]"]))
mean(pdraws[,"Sigma[school:genderM,(Intercept)]"] / sqrt(pdraws[,"Sigma[school:(Intercept),(Intercept)]"]*pdraws[,"Sigma[school:genderM,genderM]"]))

# Biased approach
estimates = colMeans(pdraws)
sqrt(estimates["Sigma[school:(Intercept),(Intercept)]"])
sqrt(estimates["Sigma[school:genderM,genderM]"])
estimates["Sigma[school:genderM,(Intercept)]"] / sqrt(estimates["Sigma[school:(Intercept),(Intercept)]"]*estimates["Sigma[school:genderM,genderM]"])

# Verify how biased approach is performed
sqrt(mean(pdraws[,"Sigma[school:(Intercept),(Intercept)]"]))
sqrt(mean(pdraws[,"Sigma[school:genderM,genderM]"]))
mean(pdraws[,"Sigma[school:genderM,(Intercept)]"]) / sqrt(mean(pdraws[,"Sigma[school:(Intercept),(Intercept)]"])*mean(pdraws[,"Sigma[school:genderM,genderM]"]))

Current Output:

VarCorr(M3_stanlmer)
Groups Name Std.Dev. Corr
school (Intercept) 9.0527
genderM 7.1805 -0.173

Expected Output:

VarCorr(M3_stanlmer)
Groups Name Std.Dev. Corr
school (Intercept) 9.0037
genderM 7.0943 -0.167

RStan Version:

packageVersion("rstan")
[1] ‘2.21.2’
packageVersion("rstanarm")
[1] ‘2.21.1’

R Version:

R version 4.0.2 (2020-06-22)

Operating System:

Microsoft 10

Start by reproducing the issue with VarCorr(M3_stanlmer) and the supplied posterior-draw calculations, then locate the VarCorr implementation in rstanarm. Compare transformation of parameter estimates with transformation of posterior draws; done means the reported standard deviations and correlations match the expected values.