Include MAD_SD of random parameters in `mod$ses`?

[emstruong]

Summary:

I'm trying to easily get the MAD_SD for the random parameters of my model, but they seem to be missing from the ses list of the model for at least stan_glmer. Could they be added? It seems like the fixed effects and random values are included... This is related to bbolker/broom.mixed#156 (comment)

Reproducible Steps:

library(rstanarm)
#> Loading required package: Rcpp
#> This is rstanarm version 2.32.1
#> - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
#> - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
#> - For execution on a local, multicore CPU with excess RAM we recommend calling
#>   options(mc.cores = parallel::detectCores())
fit <- stan_glmer(mpg ~ wt + (1|cyl) + (1+wt|gear), data = mtcars,
                  iter = 500, chains = 2)
#> 
#> SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 4e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.4 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 500 [  0%]  (Warmup)
#> Chain 1: Iteration:  50 / 500 [ 10%]  (Warmup)
#> Chain 1: Iteration: 100 / 500 [ 20%]  (Warmup)
#> Chain 1: Iteration: 150 / 500 [ 30%]  (Warmup)
#> Chain 1: Iteration: 200 / 500 [ 40%]  (Warmup)
#> Chain 1: Iteration: 250 / 500 [ 50%]  (Warmup)
#> Chain 1: Iteration: 251 / 500 [ 50%]  (Sampling)
#> Chain 1: Iteration: 300 / 500 [ 60%]  (Sampling)
#> Chain 1: Iteration: 350 / 500 [ 70%]  (Sampling)
#> Chain 1: Iteration: 400 / 500 [ 80%]  (Sampling)
#> Chain 1: Iteration: 450 / 500 [ 90%]  (Sampling)
#> Chain 1: Iteration: 500 / 500 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.402 seconds (Warm-up)
#> Chain 1:                0.212 seconds (Sampling)
#> Chain 1:                0.614 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 1.7e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:   1 / 500 [  0%]  (Warmup)
#> Chain 2: Iteration:  50 / 500 [ 10%]  (Warmup)
#> Chain 2: Iteration: 100 / 500 [ 20%]  (Warmup)
#> Chain 2: Iteration: 150 / 500 [ 30%]  (Warmup)
#> Chain 2: Iteration: 200 / 500 [ 40%]  (Warmup)
#> Chain 2: Iteration: 250 / 500 [ 50%]  (Warmup)
#> Chain 2: Iteration: 251 / 500 [ 50%]  (Sampling)
#> Chain 2: Iteration: 300 / 500 [ 60%]  (Sampling)
#> Chain 2: Iteration: 350 / 500 [ 70%]  (Sampling)
#> Chain 2: Iteration: 400 / 500 [ 80%]  (Sampling)
#> Chain 2: Iteration: 450 / 500 [ 90%]  (Sampling)
#> Chain 2: Iteration: 500 / 500 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 0.308 seconds (Warm-up)
#> Chain 2:                0.135 seconds (Sampling)
#> Chain 2:                0.443 seconds (Total)
#> Chain 2:
#> Warning: There were 1 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
fit$ses
#>           (Intercept)                    wt  b[(Intercept) cyl:4] 
#>             3.1357937             0.9469535             2.0728322 
#>  b[(Intercept) cyl:6]  b[(Intercept) cyl:8] b[(Intercept) gear:3] 
#>             1.4681377             1.9467202             0.3974553 
#>          b[wt gear:3] b[(Intercept) gear:4]          b[wt gear:4] 
#>             0.2688166             0.3394349             0.2915114 
#> b[(Intercept) gear:5]          b[wt gear:5] 
#>             0.3644062             0.2846547

^{Created on 2024-11-15 with reprex v2.1.1}

RStanARM Version:

The version of the rstanarm package you are running (e.g., from packageVersion("rstanarm")): ‘2.32.1’

R Version:

The version of R you are running (e.g., from getRversion()): ‘4.4.1’

Operating System:

Your operating system (e.g., OS X 10.11.3): Ubuntu 22.04

[jgabry]

Thanks for pointing that out. I think for now the easiest way to get the MAD for all parameters is just by computing it from the posterior draws. You could do:

apply(as.matrix(fit), 2, mad)

Does that work for you?

[emstruong]

So that definitely works for my current purposes--though a fix would still be nice for broom.mixed functionality. Was there any particular technical reason why the MAD for the random parameters was omitted or was it just a lapse?

[jgabry]

Good question, I'm not sure, it's been so long! I'll have to take a look at the source code and refresh my memory when I have a chance. I don't have time at the moment, so I'll leave this open for now until I can get to it.

[emstruong]

OK. While I'm at it, I was hoping to clarify something--AFAICT, why is it that the print method and the summary method seem to flip flop between whether they provide the robust or non-robust computation for each parameter estimate? I did look at the documentation, AFAICT the reason for the differences doesn't appear to be there...?

Even within the print method, it seems like the robust estimates for the fixed effects are provided, yet the non-robust standard deviation is used for the error terms. (?) I'm not deep enough in HMC or MLMs to really know if this is the way it should be.

Also there may or may not be a dangling sigma argument that doesn't seem to be doing anything?

---

Edit:
I'm also not sure if

apply(as.matrix(fit), 2, mad)

Will work because now that I'm comparing the estimates, this doesn't seem like the way to get the standard-deviation/MAD of the random effects?

Regardless of whether you do

apply(as.matrix(fit), 2, mad)
apply(as.matrix(fit), 2, mean)
apply(as.matrix(fit), 2, median)
apply(as.matrix(fit), 2, sd)

None of these seem to reproduce the Std.Dev. column of the error term output. Which I guess makes sense given that Std.Dev. probably comes from VarCorr, which does a lot more than just the summary statistic.

---

I'm going to guess that the reason why the standard error/MAD_SD was missing from the random parameters is because the parameterization of the random parameters in rstanarm doesn't make it easy to get the sd, mad or 95% CI of the random parameters bbolker/broom.mixed#156 (comment)

[emstruong]

As mentioned in the other thread, perhaps we could also have an argument in summary or print to facilitate having entirely robust and non-robust outputs?