Unable to save brmsfit object with saveRDS

Modeling
1

Hi all, I have simulated from the following model.

fit1 <- brm(
  bf(y ~ 0 + Intercept + a1 + a2, nl = TRUE) +
    lf(Intercept ~ 1, center = FALSE) +
    lf(as.formula(paste("a1 ~ 0 + ", paste(names(data_train)[2:16], collapse = " + "))), center = FALSE, cmc = FALSE) + 
    lf(as.formula(paste("a2 ~ 0 + ", paste(names(data_train)[17:ncol(data_train)], collapse = " + "))), center = FALSE, cmc = FALSE),
  data_train,
  family = bernoulli(link = "logit"),
  prior = prior(cauchy(0, 10), class = "b", nlpar = "Intercept") +
    prior(cauchy(0, 2.5), class = "b", nlpar = "a1") +
    prior(horseshoe(df_slab = 1, scale_slab = 2.5, par_ratio = 0.003), class = "b", nlpar = "a2"),
  warmup = 10,
  iter = 20,
  thin = 1,
  chains = 2,
  seed = 1,
  cores = 2
)

y is a binary outcome. Following recommendations in the BDA3 book, I have specified a Cauchy(0, 10) prior for the intercept and a Cauchy(0, 2.5) prior for 15 demographic predictors. Again following BDA3, the continuous predictors have mean 0 and standard deviation 0.5, whereas binary predictors have mean 0 and the two possible observed values differ by 1.
Then there is a sparsifying prior on >7000 candidate biomarkers (I expect about 20 of those 7000 to be relevant).

I would normally have used rstanarm with which I am more familiar but believe it is not possible to specify different priors on two subsets of predictors (demographic variables and candidate biomarkers) in rstanarm.

In this example I am saving a total of 20 MCMC simulations (for debugging purposes). The simulation ends and the object fit1 occupies about 280 Mbytes. I am now trying to save that as a serialised object to perform posterior predictive checks, analyse the results, etc.

Oddly, it is impossible to do so and after running

saveRDS(fit1, <path>)

followed by

fit1 <- readRDS( <path>)

on a different R session, I get the output

Error: C stack usage 7969616 is too close to the limit

I suspect the issue is with brms because I have been able to save serialized objects of much larger sizes. I am running brms on an Ubuntu machine and I have tried different versions of R (4.4, 4.3 and 3.6).

Has anybody seen this issue before?

Thanks
Víctor

2

Did you mean fit1 <- readRDS(path) instead?

You can also use the file parameter of brm to automatically save to disk. Either way, I have saved much larger models before without issue. Can you provide a reproducible example?

3

Thank you very much for your response Ax3man. Yes that is what I meant, I have edited the original post.

Here is a reproducible example. I have run it on R with the option “–max-ppsize=500000”

library(tidyverse)
library(tidymodels)
library(brms)

set.seed(1)

D <- 7350
n <- 817

data_train <- sapply(1:D, function(i) rnorm(n, 0, 0.5))
data_train <- cbind(sample(c(0, 1), size = n, replace = TRUE, prob = c(0.36, 0.64)), data_train)
colnames(data_train) <- c("y", paste0("x", 1:D))
data_train <- as_tibble(data_train)

fit1 <- brm(
  bf(y ~ 0 + Intercept + a1 + a2, nl = TRUE) +
    lf(Intercept ~ 1, center = FALSE) +
    lf(as.formula(paste("a1 ~ 0 + ", paste(names(data_train)[2:16], collapse = " + "))), center = FALSE, cmc = FALSE) + 
    lf(as.formula(paste("a2 ~ 0 + ", paste(names(data_train)[17:ncol(data_train)], collapse = " + "))), center = FALSE, cmc = FALSE),
  data_train,
  family = bernoulli(link = "logit"),
  prior = prior(cauchy(0, 10), class = "b", nlpar = "Intercept") +
    prior(cauchy(0, 2.5), class = "b", nlpar = "a1") +
    prior(horseshoe(df_slab = 1, scale_slab = 2.5, par_ratio = 0.003), class = "b", nlpar = "a2"),
  warmup = 10,
  iter = 20,
  thin = 1,
  chains = 2,
  seed = 1,
  cores = 2
)
saveRDS(fit1, "fit1.rds")

Here is my session info

R version 4.3.3 (2024-02-29)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Ubuntu 20.04.6 LTS

Matrix products: default
BLAS/LAPACK: /anaconda/envs/PAPER/lib/libopenblasp-r0.3.28.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
 [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
 [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
[10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] brms_2.20.1        Rcpp_1.0.13-1      yardstick_1.3.1    workflowsets_1.1.0
 [5] workflows_1.1.4    tune_1.2.1         rsample_1.2.1      recipes_1.1.0     
 [9] parsnip_1.2.1      modeldata_1.4.0    infer_1.0.7        dials_1.3.0       
[13] scales_1.3.0       broom_1.0.7        tidymodels_1.2.0   lubridate_1.9.3   
[17] forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4        purrr_1.0.2       
[21] readr_2.1.5        tidyr_1.3.1        tibble_3.2.1       ggplot2_3.5.1     
[25] tidyverse_2.0.0   

loaded via a namespace (and not attached):
  [1] tensorA_0.36.2.1     rstudioapi_0.17.1    magrittr_2.0.3      
  [4] vctrs_0.6.5          base64enc_0.1-3      htmltools_0.5.8.1   
  [7] distributional_0.5.0 parallelly_1.38.0    StanHeaders_2.32.10 
 [10] htmlwidgets_1.6.4    plyr_1.8.9           zoo_1.8-12          
 [13] igraph_2.0.3         mime_0.12            lifecycle_1.0.4     
 [16] iterators_1.0.14     pkgconfig_2.0.3      colourpicker_1.3.0  
 [19] Matrix_1.6-5         R6_2.5.1             fastmap_1.2.0       
 [22] future_1.34.0        shiny_1.9.1          digest_0.6.37       
 [25] colorspace_2.1-1     furrr_0.3.1          ps_1.8.1            
 [28] crosstalk_1.2.1      fansi_1.0.6          timechange_0.3.0    
 [31] abind_1.4-5          compiler_4.3.3       withr_3.0.2         
 [34] backports_1.5.0      inline_0.3.19        shinystan_2.6.0     
 [37] QuickJSR_1.4.0       pkgbuild_1.4.5       MASS_7.3-60.0.1     
 [40] lava_1.8.0           gtools_3.9.5         loo_2.8.0           
 [43] tools_4.3.3          httpuv_1.6.15        future.apply_1.11.2 
 [46] threejs_0.3.3        nnet_7.3-19          glue_1.8.0          
 [49] callr_3.7.6          nlme_3.1-165         promises_1.3.0      
 [52] grid_4.3.3           checkmate_2.3.2      reshape2_1.4.4      
 [55] generics_0.1.3       gtable_0.3.6         tzdb_0.4.0          
 [58] class_7.3-22         data.table_1.15.4    hms_1.1.3           
 [61] utf8_1.2.4           foreach_1.5.2        pillar_1.9.0        
 [64] markdown_1.13        posterior_1.6.0      later_1.3.2         
 [67] splines_4.3.3        lhs_1.2.0            lattice_0.22-6      
 [70] survival_3.7-0       tidyselect_1.2.1     miniUI_0.1.1.1      
 [73] gridExtra_2.3        stats4_4.3.3         bridgesampling_1.1-2
 [76] hardhat_1.4.0        timeDate_4041.110    matrixStats_1.4.1   
 [79] DT_0.33              rstan_2.32.6         stringi_1.8.4       
 [82] DiceDesign_1.10      codetools_0.2-20     cli_3.6.3           
 [85] RcppParallel_5.1.9   rpart_4.1.23         shinythemes_1.2.0   
 [88] xtable_1.8-4         processx_3.8.4       munsell_0.5.1       
 [91] globals_0.16.3       coda_0.19-4.1        parallel_4.3.3      
 [94] rstantools_2.4.0     gower_1.0.1          dygraphs_1.1.1.6    
 [97] bayesplot_1.11.1     Brobdingnag_1.2-9    GPfit_1.0-8         
[100] listenv_0.9.1        mvtnorm_1.3-2        ipred_0.9-15        
[103] xts_0.14.1           prodlim_2024.06.25   rlang_1.1.4         
[106] shinyjs_2.1.0

Many thanks again for your help