Suppress all output from brms in Markdown files

General
,
1

Hi all, I’ve got what might be a bug that I’m hoping someone can help with. I’m running some power simulations to determine how many people I need in my next study, and RStan/brms (I think) are causing RStudio to keep crashing, because the output gets so long.

If I only run a single model, everything is fine. If I run a smaller number of simulations (e.g., 100), still no problems. But when I run the full simulation with 1500 runs, the sampling progress for each model keeps stacking and growing longer and longer in the RMarkdown, until RStudio crashes. I can avoid this by copying and pasting the code into a script instead, but I’d rather not have to use a work-around.

I’ve set the code chunk header to have message=FALSE, warning=FALSE, results="hide", include=FALSE, refresh=0, but none of that has suppressed the output…I still get a mile-long printout of the sampling progress for each chain in each of the 1500 models.

Any ideas?

2

maybe capture all output?

output <- capture.output(fit <- brm(...))
3

So I just tried that with a single model and it does work! It produces a single model in the environment I can pull parameter estimates from. Having trouble adapting the iterative function to use it though; the command I use (parameters::model_parameters) to extract the estimates I need can’t read the captured output when its piped from update(). I get an error that says (possible class 'character' not supported'): argument is not interpreted as logical.

# works
output <- capture.output(fit <- brm(data = d1, 
                                    family = gaussian,
                                    y ~ 1 + treatment,
                                    prior = c(prior(normal(42.5, 5), class = Intercept),
                                              prior(normal(0, 10), class = b),
                                              prior(lognormal(0, 5), class = sigma)),
                                    seed = 4))

# fails
capture.output(update(fit_lin,
                        newdata = d1, 
                        seed = 1)) |>  
    parameters::model_parameters("median", ci_method = "HDI", ci=.89) |> 
    as_tibble() |> 
    select(1, 3, 5:6)

Full code for reproducability here:

pacman::p_load(brms, easystats, tidyverse, janitor, flextable)

# define the means and SD's for the idealized data
mu_struc=40
sd_struc=10

mu_UNstruc=43
sd_UNstruc=10
  
# determine the group size
n <- 50


# simulate the data
set.seed(1)

d1 <-
  tibble(group     = rep(c("structured", "unstructured"), each = n)) %>% 
  mutate(treatment = ifelse(group == "structured", 0, 1),
         y         = ifelse(group == "structured", 
                                      rnorm(n, mean = mu_struc, sd = sd_struc), 
                                      rnorm(n, mean = mu_UNstruc, sd = sd_UNstruc)))

# create practice model
fit_lin=brm(data = d1, 
      family = gaussian,
      y ~ 1 + treatment,
      prior = c(prior(normal(42.5, 5), class = Intercept),
                prior(normal(0, 10), class = b),
                prior(lognormal(0, 5), class = sigma)),
      seed = 4)


# SIMULATE
# A. Set the number of simulations
n_sim <- 1500


# B. Load function (substituting broom for easystats)
sim_d_and_fit_easystats <- function(seed, n) {
  
  set.seed(seed)
  
  d <-
    tibble(group     = rep(c("control", "treatment"), each = n)) |>  
    mutate(treatment = ifelse(group == "control", 0, 1),
           y         = ifelse(group == "control", 
                              rnorm(n, mean = mu_struc, sd = sd_struc), 
                              rnorm(n, mean = mu_UNstruc, sd = sd_UNstruc))) 
  
  update(fit_lin,
         newdata = d, 
         seed = seed) |>  
    model_parameters("median", ci_method = "HDI", ci=.89) |> as_tibble() |> select(1, 3, 5:6) |>  
    filter(Parameter == "b_treatment")
}

# C. Run simulation, with time benchmark
t1=tictoc::tic()

sim <-
  tibble(seed = 1:n_sim) |>  
  mutate(tidy = map(.x= seed, .f= sim_d_and_fit_easystats, n = 50)) |> # change sample size as necessary
  unnest(tidy) |> 
  mutate(width = CI_high - CI_low) # compute interval widths

t2=tictoc::toc()

All credit for the above goes to Solomon Kurz by the way. Absolutely amazing tutorial on his blog that I tweaked only slightly.

(In case it matters, I’m running R version 4.2.1; and brms v2.17.0)

4

didn’t try your full script but the first snippet works for me.

output <- capture.output(fit <- brm(mpg ~ cyl, data=mtcars))
output2 <- capture.output(fit2 <- update(fit, newdata=dplyr::sample_n(mtcars,10)) |> parameters::model_parameters("median", ci_method = "HDI"))

If it’s only brms output crashes the console it’s easier to debug by separating fitting & postprocessing:

output <- capture.output(fit <- brm(mpg ~ cyl, data=mtcars))
output2 <- capture.output(fit2 <- update(fit, newdata=dplyr::sample_n(mtcars,10)))
tab <- fit2 |> parameters::model_parameters("median", ci_method = "HDI")