group_level = TRUE, discrepancies between freq and Bayesian models

[DominiqueMakowski]

Frequentist models have a "EFfects" and "Group" columns, quite convenient to filter out desired group-level parameters, but Bayesian models don't

library(rstanarm)

data <- lme4::sleepstudy

model <- lme4::lmer(Reaction ~ Days + (1 | Subject), data = data)

parameters::model_parameters(model, effects = "all", group_level = TRUE) |> 
  as.data.frame() |> 
  head()
#>     Parameter Level Coefficient        SE   CI     CI_low   CI_high        t
#> 1 (Intercept)  <NA>  251.405105 9.7467163 0.95 232.169625 270.64058 25.79383
#> 2        Days  <NA>   10.467286 0.8042214 0.95   8.880127  12.05444 13.01543
#> 3 (Intercept)   308   40.783710 9.4756681 0.95  22.211742  59.35568       NA
#> 4 (Intercept)   309  -77.849554 9.4756681 0.95 -96.421522 -59.27759       NA
#> 5 (Intercept)   310  -63.108567 9.4756681 0.95 -81.680536 -44.53660       NA
#> 6 (Intercept)   330    4.406442 9.4756681 0.95 -14.165526  22.97841       NA
#>   df_error            p Effects   Group
#> 1      176 1.089809e-61   fixed        
#> 2      176 1.456699e-27   fixed        
#> 3       NA           NA  random Subject
#> 4       NA           NA  random Subject
#> 5       NA           NA  random Subject
#> 6       NA           NA  random Subject

model <- glmmTMB::glmmTMB(Reaction ~ Days + (1 | Subject), data = data)

parameters::model_parameters(model, effects = "all", group_level = TRUE) |> 
  as.data.frame() |> 
  head()
#>     Parameter Level Coefficient         SE   CI      CI_low   CI_high        z
#> 1 (Intercept)  <NA>   251.40511  9.5061750 0.95  232.773351 270.03687 26.44651
#> 2        Days  <NA>    10.46729  0.8017351 0.95    8.895919  12.03866 13.05580
#> 3 (Intercept)  <NA>    30.89542         NA 0.95   27.707999  34.44951       NA
#> 4 (Intercept)   308    40.63507 12.5348271 0.95   16.067256  65.20287       NA
#> 5 (Intercept)   309   -77.56589 12.6507490 0.95 -102.360905 -52.77088       NA
#> 6 (Intercept)   310   -62.87862 12.5961053 0.95  -87.566536 -38.19071       NA
#>   df_error             p   Component Effects   Group
#> 1      Inf 4.002539e-154 conditional   fixed        
#> 2      Inf  5.889024e-39 conditional   fixed        
#> 3       NA            NA  dispersion   fixed        
#> 4       NA            NA conditional  random Subject
#> 5       NA            NA conditional  random Subject
#> 6       NA            NA conditional  random Subject

model <- rstanarm::stan_lmer(Reaction ~ Days + (1 | Subject), data = data, refresh = 0, algorithm = "meanfield")

parameters::model_parameters(model, effects = "all", group_level = TRUE) |> 
  as.data.frame() |> 
  head()
#>                    Parameter Effects   Median   CI    CI_low  CI_high pd
#> 1                (Intercept)   fixed -25.8711 0.95 -31.58070 -18.1306  1
#> 2                       Days   fixed  10.6815 0.95   9.11947  12.0224  1
#> 3 b[(Intercept) Subject:308]  random 322.7530 0.95 302.47200 332.5520  1
#> 4 b[(Intercept) Subject:309]  random 193.9890 0.95 180.44300 215.8700  1
#> 5 b[(Intercept) Subject:310]  random 200.6020 0.95 185.84500 219.0960  1
#> 6 b[(Intercept) Subject:330]  random 275.8510 0.95 265.08600 292.1400  1
#>   ROPE_Percentage     Khat ESS Prior_Distribution Prior_Location Prior_Scale
#> 1               0 1.270982  34             normal       298.5079   140.82189
#> 2               0 1.369512  23             normal         0.0000    48.89151
#> 3               0 1.375473  30               <NA>             NA          NA
#> 4               0 1.311226  27               <NA>             NA          NA
#> 5               0 1.359794   7               <NA>             NA          NA
#> 6               0 1.392132  30               <NA>             NA          NA

model <- brms::brm(Reaction ~ Days + (1 | Subject), data = data, refresh = 0, algorithm = "meanfield")

parameters::model_parameters(model, effects = "all", group_level = TRUE) |> 
  as.data.frame() |> 
  head()
#>                  Parameter Effects   Component    Median   CI    CI_low
#> 1              b_Intercept   fixed conditional  70.40605 0.95  57.45550
#> 2                   b_Days   fixed conditional  10.71455 0.95   8.48174
#> 3 r_Subject[308,Intercept]  random conditional 247.12000 0.95 227.44300
#> 4 r_Subject[309,Intercept]  random conditional  98.98150 0.95  75.24570
#> 5 r_Subject[310,Intercept]  random conditional 124.77550 0.95 101.67000
#> 6 r_Subject[330,Intercept]  random conditional 180.86550 0.95 159.36200
#>    CI_high pd ROPE_Percentage     Khat
#> 1  81.3284  1               0 2.306278
#> 2  13.1851  1               0 2.372424
#> 3 267.6700  1               0 2.379014
#> 4 122.8430  1               0 2.341338
#> 5 148.4760  1               0 2.383970
#> 6 205.3240  1               0 2.324848

^{Created on 2022-03-26 by the reprex package (v2.0.1)}

[strengejacke]

insight::clean_parameters() comes close to what you're looking for, I guess, and was written exactly for this purpose. Mostly to support printing by components.

library(rstanarm)
#> Loading required package: Rcpp
#> This is rstanarm version 2.21.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())
data <- lme4::sleepstudy
model <- rstanarm::stan_lmer(Reaction ~ Days + (1 | Subject), data = data, refresh = 0)
insight::clean_parameters(model)
#>                                 Parameter Effects   Component
#> 1                             (Intercept)   fixed conditional
#> 2                                    Days   fixed conditional
#> 3              b[(Intercept) Subject:308]  random conditional
#> 4              b[(Intercept) Subject:309]  random conditional
#> 5              b[(Intercept) Subject:310]  random conditional
#> 6              b[(Intercept) Subject:330]  random conditional
#> 7              b[(Intercept) Subject:331]  random conditional
#> 8              b[(Intercept) Subject:332]  random conditional
#> 9              b[(Intercept) Subject:333]  random conditional
#> 10             b[(Intercept) Subject:334]  random conditional
#> 11             b[(Intercept) Subject:335]  random conditional
#> 12             b[(Intercept) Subject:337]  random conditional
#> 13             b[(Intercept) Subject:349]  random conditional
#> 14             b[(Intercept) Subject:350]  random conditional
#> 15             b[(Intercept) Subject:351]  random conditional
#> 16             b[(Intercept) Subject:352]  random conditional
#> 17             b[(Intercept) Subject:369]  random conditional
#> 18             b[(Intercept) Subject:370]  random conditional
#> 19             b[(Intercept) Subject:371]  random conditional
#> 20             b[(Intercept) Subject:372]  random conditional
#> 21 Sigma[Subject:(Intercept),(Intercept)]  random conditional
#> 22                                  sigma   fixed       sigma
#>                 Group Cleaned_Parameter
#> 1                           (Intercept)
#> 2                                  Days
#> 3  Intercept: Subject       Subject:308
#> 4  Intercept: Subject       Subject:309
#> 5  Intercept: Subject       Subject:310
#> 6  Intercept: Subject       Subject:330
#> 7  Intercept: Subject       Subject:331
#> 8  Intercept: Subject       Subject:332
#> 9  Intercept: Subject       Subject:333
#> 10 Intercept: Subject       Subject:334
#> 11 Intercept: Subject       Subject:335
#> 12 Intercept: Subject       Subject:337
#> 13 Intercept: Subject       Subject:349
#> 14 Intercept: Subject       Subject:350
#> 15 Intercept: Subject       Subject:351
#> 16 Intercept: Subject       Subject:352
#> 17 Intercept: Subject       Subject:369
#> 18 Intercept: Subject       Subject:370
#> 19 Intercept: Subject       Subject:371
#> 20 Intercept: Subject       Subject:372
#> 21    SD/Cor: Subject       (Intercept)
#> 22                                sigma

^{Created on 2022-03-26 by the reprex package (v2.0.1)}

[DominiqueMakowski]

that's awesome! should we add clean_parameters inside parameters(group_level = TRUE) to have more consistent outputs?

[strengejacke]

It's added as attributes attr(params, "parameter_info"), and used in the format.parameters_stan() method for printing. So the "printing" should be identical.

[DominiqueMakowski]

Indeed it is! Sorry I missed this :) What do you think about adding the Group / Level columns in the main output (not just as an attribute)? It'd be easier to work with

[strengejacke]

should be resolved in #686, let's wait for tests, to see whether any code breaks...

[vincentarelbundock]

Close this?