Function for repeatedly updating models not working

General
,
1

Hi all,

I recently ran through Solomon Kurz’s amazing blog post on power simulations, and I’ve been trying to translate it from brms to rstanarm since I’m only familiar with the latter package. For some reason though, when I run the stan_glm model in the early blocks just by itself it works perfectly fine, but when I try his code that iteratively updates the model and tries to run many simulations, it never finishes. The brms simulation finished in 46 seconds, whereas the rstanarm loop kept going until I aborted it a full 6 minutes later. Anyone know what’s going wrong? I’ll post the code below so you don’t have to reference the OG article…

n <- 10
z <- 7

set.seed(3)

d <- tibble(y = rbinom(n = 50, size = 1, prob = .25))


fit1 <-
  brms::brm(data = d, 
      family = binomial,
      y | trials(1) ~ 1,
      prior(normal(0, 2), class = Intercept),
      seed = 3)


#### rstanarm ####

fit2=stan_glm(y~1, 
              family = binomial(link = "logit"), 
              data=d, 
              prior = student_t(5,location=NULL, scale=NULL, autoscale = TRUE),
              #prior_intercept = normal(), 
              #prior_PD = FALSE, 
              algorithm = c("sampling"), 
              mean_PPD = TRUE,
              adapt_delta = 0.95, 
              #QR = FALSE, 
              #sparse = FALSE,
              chains=4,iter=4000,cores=4)

sim_data_fit <- function(seed, n_player) {
  
  n_trials <- 1 # do not touch
  prob_hit <- .25 # probability of event from 0 to 1
  
  set.seed(seed)
  
  d <- tibble(y = rbinom(n    = n_player, 
                         size = n_trials, 
                         prob = prob_hit))
  
  update(fit1,
         newdata = d,
         seed = seed) %>% 
  posterior_samples() %>% 
  transmute(p = inv_logit_scaled(b_Intercept)) %>% 
  median_qi() %>% 
    select(.lower:.upper)
  
}


#### rstanarm version of the above function ####

sim_data_fit_rstanarm <- function(seed, n_player) {
  
  n_trials <- 1
  prob_hit <- .25
  
  set.seed(seed)
  
  d <- tibble(y = rbinom(n    = n_player, 
                         size = n_trials, 
                         prob = prob_hit))
  
  update(fit2,
         data = d,
         seed = seed) %>% 
    posterior_samples() %>% 
    rename("Intercept"="(Intercept)") %>% # extract posterior draws
    transmute(p = inv_logit_scaled(Intercept)) %>%# convert back from log-odds scale
    tidybayes::median_qi() %>% #
    select(.lower:.upper)
  
}

EDIT: I just tried re-running the models again and letting the rstanarm version do it’s thing to see if it eventually finishes. It did, though the time difference between the two is staggering: the brms simulations finished in 1.8 minutes (on my old laptop now) vs. 25 minutes for rstanarm. Does anyone know why rstanarm took about 12 times as long as brms here? And if there’s anything that can be done to reduce that difference?

1 Like
2

Hi,
not sure what is happening - I just tried running your code and got:

system.time( sim_data_fit(2020, 3))
#   user  system elapsed 
#   0.90    0.30    4.63 

system.time( sim_data_fit_rstanarm(2020, 3))
#   user  system elapsed 
#   1.18    0.33    8.05

So quite close together, although brms would usually be expected to be a little faster as it writes Stan code specifically for each model and thus the code could be a bit more efficient.

A possible additional cause for differences is if both packages end up using different Stan version under the hood or that you are using a version of either package that has some bug affecting performance. I am using rstanarm 2.21.1 (which uses Stan 2.21) and brms 2.14.11 with CmdStan 2.25 as backend.

What are your versions? And what are the exact arguments you have used for the sim_data_XXX calls?

Best of luck with your model!

1 Like
3

So I just ran the simulations again on my desktop with similar results:{brms} finished in 54 seconds while {rstanarm} took 12.5 minutes.

My versions are as follows:
rstan: 2.21.1 (not loaded currently)
rstanarm: 2.21.1
brms: 2.14.4
StanHeaders: 2.21.0-7

I don’t see anything called CmdStan when I search the “Packages” tab in RStudio though.

I’m on Windows 10, OS build 19042.804

sim fitting functions are below; note that the only difference between the two is the “rename” line so that Solomon Kurz’s extraction method can be translated over from brms to rstanarm:

sim_data_brms <- function(seed, n_player) {
  
  n_trials <- 1 
  prob_hit <- .25 
  
  set.seed(seed)
  
  d <- tibble(y = rbinom(n    = n_player, 
                         size = n_trials, 
                         prob = prob_hit))
  
  update(fit1,
         newdata = d,
         seed = seed) %>% 
    posterior_samples() %>% 
    transmute(p = inv_logit_scaled(b_Intercept)) %>% 
    median_qi() %>% 
    select(.lower:.upper)
  
}

sim_data_rstanarm <- function(seed, n_player) {
  
  n_trials <- 1
  prob_hit <- .25
  
  set.seed(seed)
  
  d <- tibble(y = rbinom(n    = n_player, 
                         size = n_trials, 
                         prob = prob_hit))
  
  update(fit2,
         data = d,
         seed = seed) %>% 
    posterior_samples() %>% 
    rename("Intercept"="(Intercept)") %>% # extract posterior draws
    transmute(p = inv_logit_scaled(Intercept)) %>%# convert back from log-odds scale
    tidybayes::median_qi() %>% #
    select(.lower:.upper)
  
}
4

Just to make sure we run exactly the same thing: what are the seed and n_player arguments you use?

5

I used set.seed(3) for the seed, d <- tibble(y = rbinom(n = 50, size = 1, prob = .25)) to generate the data, and the following as the simulation and updating commands:

brms:

sim1 <-
  tibble(seed = 1:100) %>% 
  mutate(ci = map(seed, sim_data_fit, n_player = 50)) %>% 
  unnest()

rstanarm:

sim2 <-
  tibble(seed = 1:100) %>% 
  mutate(ci = map(seed, sim_data_fit_rstanarm, n_player = 50)) %>% 
  unnest()

EDIT: I should also add that I realized I had specified the chains and cores in the rstanarm model but not in the brms model and thought that may have contributed, so I made sure they were both consistent. I went back to the models that these commands are updating and re-run them with the line of each model as chains=4, iter=2500, cores=4, but it didn’t change anything.

6

So it seems you are running rstanarm with twice as much iterations than brms. Setting rstanarm to default number of iterations reduces the difference a bit.

So I also tried running the latest rstan 2.26. And I tried switching to the (default) rstan backend for brms to have a setup more like yours. And I get brms actually slower than rstan by about a factor of 2 :-O. In any case, it appears the bottleneck is somewhere in communicating between processes, because the actual rstanarm/brms fits (as reported by Stan progress messages) finish very quickly (< 1 s) but the spawned processes are alive for about 7 seconds.

I am also on Windows 10, so not sure what is causing the difference… Maybe @jonah or @mitzimorris would know a bit better?

7

I tried setting the cores and iterations to be the same for each and re-ran the script, but that didn’t seem to help; ended up getting the same results as before. Also, not sure how relevant this is, but brms has been giving me the following error lately:

Warning messages:
1: In system(paste(CXX, ARGS), ignore.stdout = TRUE, ignore.stderr = TRUE) :
'-E' not found

It doesn’t stop anything from running and both models still give me the same answers, so I’m not sure what this about.

I’ll post the full script I used here for a quick copy-paste repro:

pacman::p_load(tidyverse,brms,report,parameters,bayestestR, insight, see, rstanarm, tidybayes)

set.seed(3)

d <- tibble(y = rbinom(n = 50, size = 1, prob = .25))

fit1 <-
  brm(data = d, 
      family = binomial,
      y | trials(1) ~ 1,
      prior(normal(0, 2), class = Intercept),
      seed = 3,
      chains=4, iter=2500, cores=4)

fit2=stan_glm(y~1, 
              family = binomial(link = "logit"), 
              data=d, 
              prior = student_t(5,location=NULL, scale=NULL, autoscale = TRUE),
              prior_intercept = student_t(5,location = 0, scale = 2.5, autoscale = FALSE), 
              algorithm = c("sampling"), 
              mean_PPD = TRUE,
              adapt_delta = 0.95, 
              chains=4,iter=2500,cores=4)


sim_data_brms <- function(seed, n_player) {
  
  n_trials <- 1 
  prob_hit <- .25 
  
  set.seed(seed)
  
  d <- tibble(y = rbinom(n    = n_player, 
                         size = n_trials, 
                         prob = prob_hit))
  
  update(fit1,
         newdata = d,
         seed = seed) %>% 
    posterior_samples() %>% 
    transmute(p = inv_logit_scaled(b_Intercept)) %>% 
    median_qi() %>% 
    select(.lower:.upper)
  
}

sim_data_rstanarm <- function(seed, n_player) {
  
  n_trials <- 1
  prob_hit <- .25
  
  set.seed(seed)
  
  d <- tibble(y = rbinom(n    = n_player, 
                         size = n_trials, 
                         prob = prob_hit))
  
  update(fit2,
         data = d,
         seed = seed) %>% 
    posterior_samples() %>% 
    rename("Intercept"="(Intercept)") %>% # extract posterior draws
    transmute(p = inv_logit_scaled(Intercept)) %>%# convert back from log-odds scale
    tidybayes::median_qi() %>% #
    select(.lower:.upper)
  
}

t1 <- Sys.time()

sim1 <-
  tibble(seed = 1:100) %>% 
  mutate(ci = map(seed, sim_data_brms, n_player = 50)) %>% 
  unnest()

t2 <- Sys.time()
totaltime_brms=t2-t1
totaltime_brms

t1 <- Sys.time()

sim2 <-
  tibble(seed = 1:100) %>% 
  mutate(ci = map(seed, sim_data_rstanarm, n_player = 50)) %>% 
  unnest()

t2 <- Sys.time()
totaltime_rstanarm=t2-t1
totaltime_rstanarm

On this last run of the script:

totaltime_brms
Time difference of 54.73016 secs
totaltime_rstanarm
Time difference of 12.69326 mins