[Bug]: Very slow execution of report for BRMS models

[M-Colley]

Describe the bug

When trying to report larger BRMS models (e.g. with more data and 60,000 iterations), the report function is very slow.

For example: the duration to fit the model on a powerful PC is 140s. The generated report with report takes approx. 3 hours.

Reproducible example

I have the newest report version of 0.6.1.

The example requires custom data, which I am happy to share over email.

[robitalec]

Here's a reprex using the model example from {brms}.
report_model seems to refit the brms model, which may be related to above (but not sure why the original model fit is fast and report is slow..)

library(brms)
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.23.0). 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(report)

fit1 <- brm(
  count ~ zBase * Trt + (1|patient),
  data = epilepsy, family = poisson(),
  prior = prior(normal(0, 10), class = b) +
    prior(cauchy(0, 2), class = sd)
)
#> Compiling Stan program...
#> Start sampling
#> 
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 2.8e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 1.035 seconds (Warm-up)
#> Chain 1:                0.712 seconds (Sampling)
#> Chain 1:                1.747 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 1.4e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 0.942 seconds (Warm-up)
#> Chain 2:                0.697 seconds (Sampling)
#> Chain 2:                1.639 seconds (Total)
#> Chain 2: 
#> 
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
#> Chain 3: 
#> Chain 3: Gradient evaluation took 1.2e-05 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3: 
#> Chain 3: 
#> Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 3: 
#> Chain 3:  Elapsed Time: 1.019 seconds (Warm-up)
#> Chain 3:                0.71 seconds (Sampling)
#> Chain 3:                1.729 seconds (Total)
#> Chain 3: 
#> 
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
#> Chain 4: 
#> Chain 4: Gradient evaluation took 1.2e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4: 
#> Chain 4: 
#> Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 4: 
#> Chain 4:  Elapsed Time: 0.958 seconds (Warm-up)
#> Chain 4:                0.709 seconds (Sampling)
#> Chain 4:                1.667 seconds (Total)
#> Chain 4:


report_model(fit1)
#> Start sampling
#> 
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 1.8e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.983 seconds (Warm-up)
#> Chain 1:                0.725 seconds (Sampling)
#> Chain 1:                1.708 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 1.4e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 0.886 seconds (Warm-up)
#> Chain 2:                0.717 seconds (Sampling)
#> Chain 2:                1.603 seconds (Total)
#> Chain 2: 
#> 
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
#> Chain 3: 
#> Chain 3: Gradient evaluation took 1.2e-05 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3: 
#> Chain 3: 
#> Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 3: 
#> Chain 3:  Elapsed Time: 0.923 seconds (Warm-up)
#> Chain 3:                0.704 seconds (Sampling)
#> Chain 3:                1.627 seconds (Total)
#> Chain 3: 
#> 
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
#> Chain 4: 
#> Chain 4: Gradient evaluation took 1.2e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4: 
#> Chain 4: 
#> Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 4: 
#> Chain 4:  Elapsed Time: 1.008 seconds (Warm-up)
#> Chain 4:                0.817 seconds (Sampling)
#> Chain 4:                1.825 seconds (Total)
#> Chain 4:
#> Response residuals not available to calculate mean square error. (R)MSE
#>   is probably not reliable.
#> Bayesian poisson mixed model (estimated using MCMC sampling with 4 chains of 2000 iterations and a warmup of 1000) to predict count with zBase and Trt (formula: count ~ zBase * Trt). The model included patient as random effect (formula: ~1 | patient). Priors over parameters were set as student_t (location = 1.40, scale = 2.50) distributions
#> Bayesian poisson mixed model (estimated using MCMC sampling with 4 chains of 2000 iterations and a warmup of 1000) to predict count with zBase and Trt (formula: count ~ zBase * Trt). The model included patient as random effect (formula: ~1 | patient). Priors over parameters were set as normal (mean = 0.00, SD = 10.00) distributions
#> Bayesian poisson mixed model (estimated using MCMC sampling with 4 chains of 2000 iterations and a warmup of 1000) to predict count with zBase and Trt (formula: count ~ zBase * Trt). The model included patient as random effect (formula: ~1 | patient). Priors over parameters were set as normal (mean = 0.00, SD = 10.00) distributions
#> Bayesian poisson mixed model (estimated using MCMC sampling with 4 chains of 2000 iterations and a warmup of 1000) to predict count with zBase and Trt (formula: count ~ zBase * Trt). The model included patient as random effect (formula: ~1 | patient). Priors over parameters were set as normal (mean = 0.00, SD = 10.00) distributions
#> Bayesian poisson mixed model (estimated using MCMC sampling with 4 chains of 2000 iterations and a warmup of 1000) to predict count with zBase and Trt (formula: count ~ zBase * Trt). The model included patient as random effect (formula: ~1 | patient). Priors over parameters were set as cauchy (location = , scale = ) distributions

^{Created on 2025-09-15 with reprex v2.1.1}

[rempsyc]

Merci @robitalec !

[rempsyc]

Here's a profiling visualization (using profvis::profvis()) of which parts of report consume the most time and resources with the above example. Might help during future investigations.

[rempsyc]

Not sure if related but found this:

That is correct. If a stanfit or stanmodel object is assigned to a symbol (that can easily be found), then calling the same Stan program does not recompile the model and instead fishes the compiled model out of the existing stanfit or stanmodel object. Otherwise, it has to recompile because loading the same dynamic shared object off of the disk crashes R. (source)

I don't know if that's what's happening here, i.e., recompiling.

[mattansb]

The model shouldn't recompile, but I believe it is being completely resampled in standardize_parameters() because the default it method = "refit".

[rempsyc]

I tried to benchmark it using different solutions...

Single timed reprex:

library(brms)
library(report)
packageVersion("report")
#> [1] '0.6.1.4'

start_time_fit <- Sys.time()
fit1 <- brm(
  count ~ zBase * Trt + (1 | patient),
  data = epilepsy, family = poisson(),
  prior = prior(normal(0, 10), class = b) +
    prior(cauchy(0, 2), class = sd)
)
end_time_fit <- Sys.time()
end_time_fit - start_time_fit
#> Time difference of 2.065662 mins

start_time_report <- Sys.time()
report_model(fit1)
end_time_report <- Sys.time()
end_time_report - start_time_report
#> Time difference of 45.48624 secs

^{Created on 2025-09-17 with reprex v2.1.1}

So report() is like 45/2 = 22.5 times slower.

rbenchmark (took a long time, defaults to 100 replications...):

library(rbenchmark)
library(brms)
library(report)

benchmark(
  fit1 <- brm(
    count ~ zBase * Trt + (1 | patient),
    data = epilepsy, family = poisson(),
    prior = prior(normal(0, 10), class = b) +
      prior(cauchy(0, 2), class = sd)
  ),
  report_model(fit1)
)

test
1 fit1 <- brm(count ~ zBase * Trt + (1 | patient), data = epilepsy, family = poisson(), prior = prior(normal(0, 10), class = b) + prior(cauchy(0, 2), class = sd))
2 report_model(fit1)
  replications  elapsed relative user.self sys.self user.child sys.child
1          100 10223.64    2.463   3282.53    73.61         NA        NA
2          100  4150.99    1.000   4039.81    83.42         NA        NA

Similar results here but reverse? 1's supposed to be the fit and 2 report? Similar using microbenchmark

library(microbenchmark)
library(brms)
library(report)

microbenchmark(
  fit1 <- brm(
    count ~ zBase * Trt + (1 | patient),
    data = epilepsy, family = poisson(),
    prior = prior(normal(0, 10), class = b) +
      prior(cauchy(0, 2), class = sd)
  ),
  report_model(fit1)
)

Unit: seconds
                                                                                                                                                                           expr
 fit1 <- brm(count ~ zBase * Trt + (1 | patient), data = epilepsy,      family = poisson(), prior = prior(normal(0, 10), class = b) +          prior(cauchy(0, 2), class = sd))
                                                                                                                                                             report_model(fit1)
      min        lq     mean    median        uq       max neval cld
 95.87804 107.08333 110.6512 109.74569 114.11338 159.57903   100  a 
 39.29918  42.63507  44.3378  43.74013  45.95419  54.25919   100   b

Finally, using bench::mark()

library(bench)
library(brms)
library(report)

mark(
  fit1 <- brm(
    count ~ zBase * Trt + (1 | patient),
    data = epilepsy, family = poisson(),
    prior = prior(normal(0, 10), class = b) +
      prior(cauchy(0, 2), class = sd)
  ),
  report_model(fit1),
  check = FALSE
)

# A tibble: 2 × 13
  expression            min median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time result memory     time       gc      
  <bch:expr>         <bch:> <bch:>     <dbl> <bch:byt>    <dbl> <int> <dbl>   <bch:tm> <list> <list>     <list>     <list>  
1 fit1 <- brm(count…  1.77m  1.77m   0.00940  320.96MB  0.00940     1     1      1.77m <NULL> <Rprofmem> <bench_tm> <tibble>
2 report_model(fit1) 43.71s 43.71s   0.0229     2.95GB  0.0915      1     4     43.71s <NULL> <Rprofmem> <bench_tm> <tibble>

Same as other benchmarking approaches (report faster than fit). I don't know if this is an artifact of the iterations?