Stanfit object fit inside function explodes in size when saved to .rds

General
1

I’m using rstanarm to fit stan_glm models inside a function. Running into a problem where the size of the saved stanfit object explodes in size when saved to .rds, but only when the model is fit inside a function. The issue seems to be that the stanfit object is storing a copy of the local environment, which is then getting saved to disk with write_rds? Manually removing large objects inside the function more or less solves the problem, but that’s a pretty clunky solution, so wondering if anyone has advice for a more elegant way to resolve this issue? Toy reprex below (warning this will write some .rds files to disk, I remove them at the end of the example but heads up anyway)

library(tidyverse)
library(rstanarm)
#> Loading required package: Rcpp
#> rstanarm (Version 2.19.3, packaged: 2020-02-11 05:16:41 UTC)
#> - Do not expect the default priors to remain the same in future rstanarm versions.
#> Thus, R scripts should specify priors explicitly, even if they are just the defaults.
#> - For execution on a local, multicore CPU with excess RAM we recommend calling
#> options(mc.cores = parallel::detectCores())
#> - bayesplot theme set to bayesplot::theme_default()
#>    * Does _not_ affect other ggplot2 plots
#>    * See ?bayesplot_theme_set for details on theme setting
library(gapminder)

# create a largeish object
test <- matrix(data = rnorm(10000), nrow = 10000/2, ncol = 10000/2)

# fit model in the global environment

a = stan_glm(lifeExp ~ gdpPercap, data = gapminder, refresh =0)

print(object.size(a), unit = "Mb")
#> 1.4 Mb


# fit model inside function , passing but not using largeish object
memfoo <- function(gap, testy, clean = FALSE){
  
  d <- testy
  
  if (clean){
    rm(d,testy)
  }
  
  a <- stan_glm(lifeExp ~ gdpPercap, data = gap, refresh = 0)
  
}

b <- memfoo(gapminder, test)

# fit model again, but removing large obects from the environment before running
d <- memfoo(gapminder, test, clean = TRUE)

print(object.size(a), unit = "Mb")
#> 1.4 Mb

print(object.size(b), unit = "Mb")
#> 1.4 Mb

print(object.size(d), unit = "Mb")
#> 1.4 Mb

# all same size in memory

# write to .rds

write_rds(a,"a.rds")

write_rds(b,"b.rds")

write_rds(d,"d.rds")

# rstan object run in function with largeish object in environment is 45 times bigger than same regression 
# fit outside function!
file.size("a.rds")
#> [1] 9026457

file.size("b.rds")
#> [1] 410011307

file.size("d.rds")
#> [1] 10011197

file.size("b.rds") / file.size("a.rds")
#> [1] 45.42328

file.remove(c("a.rds", "b.rds", "d.rds"))
#> [1] TRUE TRUE TRUE

Created on 2020-03-11 by the reprex package (v0.3.0)

1 Like
2

I noticed this post a little while ago for a tool which trims unnecessary bits from a stanfit object. Does it help with your problem?

3

Interesting, thanks! Doesn’t quite solve it though: looks like shredder is designed to reduce the size of the stanfit object. But, the problem is in the saved environment, which is in the terms / formula part of the rstanreg object. But poking around where exactly the problem was coming from using butcher led me to try a new route: specifying the environment to use for the model inside the function (global environment + only the data needed to actually fit the model), and that more or less solves it.

env <- new.env(parent = .GlobalEnv)

  env$gap <- gap
  
  a <- with(env,{stan_glm(lifeExp ~ gdpPercap, data = gap, refresh = 0)})

Easier to just add the data for the model to the environment than remove all the other objects from the environment inside the function

library(tidyverse)
library(rstanarm)
#> Loading required package: Rcpp
#> rstanarm (Version 2.19.3, packaged: 2020-02-11 05:16:41 UTC)
#> - Do not expect the default priors to remain the same in future rstanarm versions.
#> Thus, R scripts should specify priors explicitly, even if they are just the defaults.
#> - For execution on a local, multicore CPU with excess RAM we recommend calling
#> options(mc.cores = parallel::detectCores())
#> - bayesplot theme set to bayesplot::theme_default()
#>    * Does _not_ affect other ggplot2 plots
#>    * See ?bayesplot_theme_set for details on theme setting
library(gapminder)
library(butcher)

# create a largeish object
test <- matrix(data = rnorm(10000), nrow = 10000/2, ncol = 10000/2)

# fit model in the global environment

a = stan_glm(lifeExp ~ gdpPercap, data = gapminder, refresh =0)

# fit model inside function , passing but not using largeish object
memfoo <- function(gap, testy, clean = FALSE){
  
  d <- testy
  
  if (clean){
  
  
  env <- new.env(parent = .GlobalEnv)

  env$gap <- gap
  
  a <- with(env,{stan_glm(lifeExp ~ gdpPercap, data = gap, refresh = 0)})
  
  } else {
    
    a <- stan_glm(lifeExp ~ gdpPercap, data = gap, refresh = 0)
  }
}

b <- memfoo(gapminder, test)

d <- memfoo(gapminder, test, clean = TRUE)

butcher::weigh(a)
#> # A tibble: 50 x 2
#>    object              size
#>    <chr>              <dbl>
#>  1 stanfit           3.87  
#>  2 residuals         0.123 
#>  3 data.country      0.017 
#>  4 fitted.values     0.0143
#>  5 linear.predictors 0.0143
#>  6 y                 0.0143
#>  7 model.lifeExp     0.0137
#>  8 model.gdpPercap   0.0137
#>  9 data.lifeExp      0.0137
#> 10 data.gdpPercap    0.0137
#> # … with 40 more rows

butcher::weigh(b)
#> # A tibble: 50 x 2
#>    object                size
#>    <chr>                <dbl>
#>  1 terms             204.    
#>  2 formula           204.    
#>  3 stanfit             3.87  
#>  4 residuals           0.123 
#>  5 data.country        0.017 
#>  6 fitted.values       0.0143
#>  7 linear.predictors   0.0143
#>  8 y                   0.0143
#>  9 model.lifeExp       0.0137
#> 10 model.gdpPercap     0.0137
#> # … with 40 more rows

butcher::weigh(d)
#> # A tibble: 50 x 2
#>    object              size
#>    <chr>              <dbl>
#>  1 stanfit           3.87  
#>  2 residuals         0.123 
#>  3 terms             0.0700
#>  4 formula           0.0677
#>  5 data.country      0.017 
#>  6 fitted.values     0.0143
#>  7 linear.predictors 0.0143
#>  8 y                 0.0143
#>  9 model.lifeExp     0.0137
#> 10 model.gdpPercap   0.0137
#> # … with 40 more rows

write_rds(a, path = "a.rds")

write_rds(b, path = "b.rds")

write_rds(d, path = "d.rds")

file.size("a.rds")
#> [1] 9110179

file.size("b.rds")
#> [1] 410171586

file.size("d.rds")
#> [1] 9167589

file.size("d.rds") / file.size("a.rds")
#> [1] 1.006302

file.remove(c("a.rds", "b.rds", "d.rds"))
#> [1] TRUE TRUE TRUE

Created on 2020-03-12 by the reprex package (v0.3.0)