Improper order of operations for VarCorr

[mikeguggis]

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

[jgabry]

I just noticed that this issue didn't get any attention for whatever reason, sorry about that! @bgoodri I think @mikeguggis is right about this. Do you remember if there is a reason the implementation is the way it is instead of the way @mikeguggis suggests?

[bgoodri]

What rstanarm does has always been slightly wrong for the reasons outlined. I don't remember the details of why it is that way off the top of my head. I think it had something to do with making VarCorr compatible with how lme4 does it, but it might have been something else.

On Fri, May 7, 2021 at 4:18 PM Jonah Gabry @.***> wrote:

I just noticed that this issue didn't get any attention for whatever reason, sorry about that! @bgoodri https://github.com/bgoodri I think @mikeguggis https://github.com/mikeguggis is right about this. Do you remember if there is a reason the implementation is the way it is instead of the way @mikeguggis https://github.com/mikeguggis suggests?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/stan-dev/rstanarm/issues/470#issuecomment-834748282, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAZ2XKTHR6S5X7372WEYNMTTMRDJ3ANCNFSM4R2TRAEA .

[jgabry]

On Fri, May 7, 2021 at 2:29 PM bgoodri @.***> wrote:

What rstanarm does has always been slightly wrong for the reasons outlined. I don't remember the details of why it is that way off the top of my head. I think it had something to do with making VarCorr compatible with how lme4 does it, but it might have been something else.

I don't remember either but looking at it now we don't call lme4's version of VarCorr internally so I'm not sure why we would need compatibility. Or do you mean that we tried to copy the internals of lme4's VarCorr? Either way it seems preferable for it to do the right thing, regardless of lme4. Until one of us has a chance to fix this should we just add a warning about this behavior when VarCorr is called? The current behavior of silently doing the wrong thing seems like the worst possible option.

[allisonmcampbell]

Has this issue been resolved? I'm using VarCorr, and I haven't seen any sort of warning. I just stumbled onto this thread because the VarCorr estimates looked curious, and I was having trouble reproducing them with the draws.