Understanding brms::save_pars()

matti · February 8, 2022, 3:05am

Hi!

I am trying to fit multinomial models with many levels of grouping factors, resulting in very large numbers of “random intercepts”. However I can’t wrap my head around how to use save_pars() to not save the group-level coefficients. Even when I specify save_pars = save_pars(group = FALSE), the resulting model objects contain the group level intercepts. Could someone advise me on how to not save the group-level effects (or otherwise make the resulting models smaller hehe)? Thanks!

Reproducible example:

library(brms)
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.16.7). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
#> 
#> Attaching package: 'brms'
#> The following object is masked from 'package:stats':
#> 
#>     ar
# library(tidyverse)

N <- 100
dat <- data.frame(
  id = factor(1:N),
  y1 = rbinom(N, 100, 0.1), 
  y2 = rbinom(N, 100, 0.2), 
  y3 = rbinom(N, 100, 0.3)
)
dat$size <- with(dat, y1 + y2 + y3)
fit <- brm(
  bf(cbind(y1, y2, y3) | trials(size) ~ 1 + (1 |p| id)), 
  family = multinomial(),
  data = dat, 
  save_pars = save_pars(group = FALSE),
  iter = 10, chains = 1,
  backend = "cmdstanr"
)
#> Start sampling
#> Running MCMC with 1 chain...
#> 
#> Chain 1 WARNING: No variance estimation is 
#> Chain 1          performed for num_warmup < 20 
#> Chain 1 Iteration: 1 / 10 [ 10%]  (Warmup) 
#> Chain 1 Iteration: 6 / 10 [ 60%]  (Sampling) 
#> Chain 1 Iteration: 10 / 10 [100%]  (Sampling) 
#> Chain 1 finished in 0.0 seconds.
#> 
#> Warning: 5 of 5 (100.0%) transitions ended with a divergence.
#> This may indicate insufficient exploration of the posterior distribution.
#> Possible remedies include: 
#>   * Increasing adapt_delta closer to 1 (default is 0.8) 
#>   * Reparameterizing the model (e.g. using a non-centered parameterization)
#>   * Using informative or weakly informative prior distributions

# This should be empty?
coef(fit)[[1]][1:5,,1]
#>     Estimate Est.Error       Q2.5      Q97.5
#> 1 -0.7857120         0 -0.7857120 -0.7857120
#> 2 -0.4242896         0 -0.4242896 -0.4242896
#> 3 -0.7172080         0 -0.7172080 -0.7172080
#> 4 -0.2784950         0 -0.2784950 -0.2784950
#> 5 -0.3103090         0 -0.3103090 -0.3103090

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

Session info

sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.1.2 (2021-11-01)
#>  os       macOS Monterey 12.1
#>  system   aarch64, darwin20
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       Asia/Seoul
#>  date     2022-02-08
#>  pandoc   2.17.1.1 @ /Applications/RStudio.app/Contents/MacOS/quarto/bin/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package        * version    date (UTC) lib source
#>  abind            1.4-5      2016-07-21 [1] CRAN (R 4.1.0)
#>  assertthat       0.2.1      2019-03-21 [1] CRAN (R 4.1.0)
#>  backports        1.4.1      2021-12-13 [1] CRAN (R 4.1.1)
#>  base64enc        0.1-3      2015-07-28 [1] CRAN (R 4.1.0)
#>  bayesplot        1.8.1      2021-06-14 [1] CRAN (R 4.1.0)
#>  bridgesampling   1.1-2      2021-04-16 [1] CRAN (R 4.1.0)
#>  brms           * 2.16.7     2022-02-08 [1] Github (paul-buerkner/brms@3164328)
#>  Brobdingnag      1.2-7      2022-02-03 [1] CRAN (R 4.1.1)
#>  callr            3.7.0      2021-04-20 [1] CRAN (R 4.1.0)
#>  checkmate        2.0.0      2020-02-06 [1] CRAN (R 4.1.1)
#>  cli              3.1.1      2022-01-20 [1] CRAN (R 4.1.2)
#>  cmdstanr         0.4.0.9001 2022-01-25 [1] Github (stan-dev/cmdstanr@a2a97d9)
#>  coda             0.19-4     2020-09-30 [1] CRAN (R 4.1.0)
#>  codetools        0.2-18     2020-11-04 [1] CRAN (R 4.1.2)
#>  colorspace       2.0-2      2021-06-24 [1] CRAN (R 4.1.1)
#>  colourpicker     1.1.1      2021-10-04 [1] CRAN (R 4.1.1)
#>  crayon           1.4.2      2021-10-29 [1] CRAN (R 4.1.1)
#>  crosstalk        1.2.0      2021-11-04 [1] CRAN (R 4.1.1)
#>  curl             4.3.2      2021-06-23 [1] CRAN (R 4.1.0)
#>  data.table       1.14.2     2021-09-27 [1] CRAN (R 4.1.1)
#>  DBI              1.1.2      2021-12-20 [1] CRAN (R 4.1.1)
#>  digest           0.6.29     2021-12-01 [1] CRAN (R 4.1.1)
#>  distributional   0.3.0      2022-01-05 [1] CRAN (R 4.1.1)
#>  dplyr            1.0.7      2021-06-18 [1] CRAN (R 4.1.0)
#>  DT               0.20       2021-11-15 [1] CRAN (R 4.1.1)
#>  dygraphs         1.1.1.6    2018-07-11 [1] CRAN (R 4.1.0)
#>  ellipsis         0.3.2      2021-04-29 [1] CRAN (R 4.1.0)
#>  emmeans          1.7.2      2022-01-04 [1] CRAN (R 4.1.1)
#>  estimability     1.3        2018-02-11 [1] CRAN (R 4.1.0)
#>  evaluate         0.14       2019-05-28 [1] CRAN (R 4.1.0)
#>  fansi            1.0.2      2022-01-14 [1] CRAN (R 4.1.1)
#>  farver           2.1.0      2021-02-28 [1] CRAN (R 4.1.0)
#>  fastmap          1.1.0      2021-01-25 [1] CRAN (R 4.1.0)
#>  fs               1.5.2      2021-12-08 [1] CRAN (R 4.1.1)
#>  generics         0.1.2      2022-01-31 [1] CRAN (R 4.1.1)
#>  ggplot2          3.3.5      2021-06-25 [1] CRAN (R 4.1.1)
#>  ggridges         0.5.3      2021-01-08 [1] CRAN (R 4.1.1)
#>  glue             1.6.1      2022-01-22 [1] CRAN (R 4.1.2)
#>  gridExtra        2.3        2017-09-09 [1] CRAN (R 4.1.1)
#>  gtable           0.3.0      2019-03-25 [1] CRAN (R 4.1.1)
#>  gtools           3.9.2      2021-06-06 [1] CRAN (R 4.1.0)
#>  highr            0.9        2021-04-16 [1] CRAN (R 4.1.0)
#>  htmltools        0.5.2      2021-08-25 [1] CRAN (R 4.1.1)
#>  htmlwidgets      1.5.4      2021-09-08 [1] CRAN (R 4.1.1)
#>  httpuv           1.6.5      2022-01-05 [1] CRAN (R 4.1.1)
#>  igraph           1.2.11     2022-01-04 [1] CRAN (R 4.1.1)
#>  inline           0.3.19     2021-05-31 [1] CRAN (R 4.1.0)
#>  jsonlite         1.7.3      2022-01-17 [1] CRAN (R 4.1.2)
#>  knitr            1.37       2021-12-16 [1] CRAN (R 4.1.1)
#>  later            1.3.0      2021-08-18 [1] CRAN (R 4.1.1)
#>  lattice          0.20-45    2021-09-22 [1] CRAN (R 4.1.2)
#>  lifecycle        1.0.1      2021-09-24 [1] CRAN (R 4.1.1)
#>  loo              2.4.1      2020-12-09 [1] CRAN (R 4.1.0)
#>  magrittr         2.0.2      2022-01-26 [1] CRAN (R 4.1.1)
#>  markdown         1.1        2019-08-07 [1] CRAN (R 4.1.0)
#>  Matrix           1.4-0      2021-12-08 [1] CRAN (R 4.1.1)
#>  matrixStats      0.61.0     2021-09-17 [1] CRAN (R 4.1.1)
#>  mime             0.12       2021-09-28 [1] CRAN (R 4.1.1)
#>  miniUI           0.1.1.1    2018-05-18 [1] CRAN (R 4.1.0)
#>  munsell          0.5.0      2018-06-12 [1] CRAN (R 4.1.0)
#>  mvtnorm          1.1-3      2021-10-08 [1] CRAN (R 4.1.1)
#>  nlme             3.1-155    2022-01-13 [1] CRAN (R 4.1.1)
#>  pillar           1.7.0      2022-02-01 [1] CRAN (R 4.1.1)
#>  pkgbuild         1.3.1      2021-12-20 [1] CRAN (R 4.1.1)
#>  pkgconfig        2.0.3      2019-09-22 [1] CRAN (R 4.1.0)
#>  plyr             1.8.6      2020-03-03 [1] CRAN (R 4.1.0)
#>  posterior        1.2.0      2022-01-05 [1] CRAN (R 4.1.1)
#>  prettyunits      1.1.1      2020-01-24 [1] CRAN (R 4.1.0)
#>  processx         3.5.2      2021-04-30 [1] CRAN (R 4.1.0)
#>  promises         1.2.0.1    2021-02-11 [1] CRAN (R 4.1.0)
#>  ps               1.6.0      2021-02-28 [1] CRAN (R 4.1.0)
#>  purrr            0.3.4      2020-04-17 [1] CRAN (R 4.1.0)
#>  R6               2.5.1      2021-08-19 [1] CRAN (R 4.1.1)
#>  Rcpp           * 1.0.8      2022-01-13 [1] CRAN (R 4.1.1)
#>  RcppParallel     5.1.5      2022-01-05 [1] CRAN (R 4.1.1)
#>  reprex           2.0.1      2021-08-05 [1] CRAN (R 4.1.1)
#>  reshape2         1.4.4      2020-04-09 [1] CRAN (R 4.1.0)
#>  rlang            1.0.1      2022-02-03 [1] CRAN (R 4.1.1)
#>  rmarkdown        2.11       2021-09-14 [1] CRAN (R 4.1.1)
#>  rsconnect        0.8.25     2021-11-19 [1] CRAN (R 4.1.1)
#>  rstan            2.26.6     2022-01-25 [1] local
#>  rstantools       2.1.1      2020-07-06 [1] CRAN (R 4.1.0)
#>  rstudioapi       0.13       2020-11-12 [1] CRAN (R 4.1.0)
#>  scales           1.1.1      2020-05-11 [1] CRAN (R 4.1.0)
#>  sessioninfo      1.2.2      2021-12-06 [1] CRAN (R 4.1.1)
#>  shiny            1.7.1      2021-10-02 [1] CRAN (R 4.1.1)
#>  shinyjs          2.1.0      2021-12-23 [1] CRAN (R 4.1.1)
#>  shinystan        2.5.0      2018-05-01 [1] CRAN (R 4.1.0)
#>  shinythemes      1.2.0      2021-01-25 [1] CRAN (R 4.1.0)
#>  StanHeaders      2.26.6     2022-01-25 [1] local
#>  stringi          1.7.6      2021-11-29 [1] CRAN (R 4.1.1)
#>  stringr          1.4.0      2019-02-10 [1] CRAN (R 4.1.1)
#>  tensorA          0.36.2     2020-11-19 [1] CRAN (R 4.1.0)
#>  threejs          0.3.3      2020-01-21 [1] CRAN (R 4.1.0)
#>  tibble           3.1.6      2021-11-07 [1] CRAN (R 4.1.1)
#>  tidyselect       1.1.1      2021-04-30 [1] CRAN (R 4.1.0)
#>  utf8             1.2.2      2021-07-24 [1] CRAN (R 4.1.0)
#>  V8               4.1.0      2022-02-06 [1] CRAN (R 4.1.2)
#>  vctrs            0.3.8      2021-04-29 [1] CRAN (R 4.1.0)
#>  withr            2.4.3      2021-11-30 [1] CRAN (R 4.1.1)
#>  xfun             0.29       2021-12-14 [1] CRAN (R 4.1.1)
#>  xtable           1.8-4      2019-04-21 [1] CRAN (R 4.1.0)
#>  xts              0.12.1     2020-09-09 [1] CRAN (R 4.1.0)
#>  yaml             2.2.2      2022-01-25 [1] CRAN (R 4.1.1)
#>  zoo              1.8-9      2021-03-09 [1] CRAN (R 4.1.0)
#> 
#>  [1] /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library
#> 
#> ──────────────────────────────────────────────────────────────────────────────

franzsf · February 8, 2022, 3:28am

Try using the rstan backend. Thats the only way I’ve been able to successfully use the save_pars argument.

matti · February 8, 2022, 3:32am

Hi @franzsf , whoa I did not consider this. I can confirm that it works as expected with backend = "rstan". But then I can’t use multithreading 😭.

@paul.buerkner is this (save_pars() works with rstan but not cmdstanr) the intended behaviour?

Thanks!

franzsf · February 8, 2022, 3:37am

While we’re asking for nice things, how about being able to sample_new_levels (=gaussian?) after discarding group-level effects? 😁

paul.buerkner · February 8, 2022, 9:06am

This is a known issue that has to something to do with cmdstanr not yet supporting the exclude way of not storing parameters. I think it likely has to be done from the brms side but I am still struggeling a bit how to implement this efficiently without first loading all parameters into RAM and then writing them to CSV again only for them to be loaded back in RAM by rstan for storage inside a brmsfit object.

matti · February 9, 2022, 1:31am

I see, thanks for clarifying!

Nic · May 19, 2022, 12:38pm

Thanks for clarifying! Do you already know if this will be updated in future versions of brms soon? I would also like to use threading and the chkptstanr package, both of which only work with cmdstanr, but the resulting objects are just far too large for this to work (e.g., >8GB).

Topic		Replies	Views
Brms sampling in multi-level multinomial models Modeling specification , brms , meta-analysis	8	1641	March 16, 2021
Understanding multilevel multinomial models generated by brms brms	8	1796	March 12, 2020
A group-level intercept model using a multivariate brms formula? brms	3	924	September 10, 2018
In brms: Two categorical predictors with unique intercepts Interfaces	2	431	March 3, 2019
Brms categorical model, hypothesis testing brms	4	1058	October 5, 2018

Understanding brms::save_pars()

Related Topics