Brms with categorical model crashes R after package upgrade to brms_2.13.4

Yesterday I upgraded brms and rstan, and my ordinal regression models stopped to work.
Here a reproducible example:

data("stemcell", package = "brglm2")
fit_sc1 = brm(formula = research ~ 1 + religion, data = stemcell, family = cumulative("probit"))

Here’s the output:

Compiling Stan program...
Start sampling

SAMPLING FOR MODEL '7c133df9b554e766be6089f9fd09b46b' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 3.5e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.35 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)

and that is: “Selection” stands for when R crashes, but the classical options are not printed…
No debug, no hints… I tried on another machine without the upgrade and all runs smoothly.
Any help is appreciated

R version 4.0.2 (2020-06-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Arch Linux

Matrix products: default
BLAS:   /usr/lib/
LAPACK: /usr/lib/

 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            

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

other attached packages:
[1] brms_2.13.4 Rcpp_1.0.5 

loaded via a namespace (and not attached):
 [1] Brobdingnag_1.2-6    jsonlite_1.7.0       gtools_3.8.2        
 [4] StanHeaders_2.21.0-5 RcppParallel_5.0.2   threejs_0.3.3       
 [7] shiny_1.5.0          assertthat_0.2.1     stats4_4.0.2        
[10] backports_1.1.8      pillar_1.4.6         lattice_0.20-41     
[13] glue_1.4.1           digest_0.6.25        promises_1.1.1      
[16] colorspace_1.4-1     htmltools_0.5.0      httpuv_1.5.4        
[19] Matrix_1.2-18        plyr_1.8.6           dygraphs_1.1.1.6    
[22] pkgconfig_2.0.3      rstan_2.21.1         purrr_0.3.4         
[25] xtable_1.8-4         mvtnorm_1.1-1        scales_1.1.1        
[28] processx_3.4.3       later_1.1.0.1        tibble_3.0.3        
[31] bayesplot_1.7.2      generics_0.0.2       ggplot2_3.3.2       
[34] ellipsis_0.3.1       DT_0.14              shinyjs_1.1         
[37] cli_2.0.2            magrittr_1.5         crayon_1.3.4        
[40] mime_0.9             ps_1.3.3             fansi_0.4.1         
[43] nlme_3.1-148         xts_0.12-0           pkgbuild_1.1.0      
[46] colourpicker_1.0     rsconnect_0.8.16     tools_4.0.2         
[49] loo_2.3.1            prettyunits_1.1.1    lifecycle_0.2.0     
[52] matrixStats_0.56.0   stringr_1.4.0        V8_3.2.0            
[55] munsell_0.5.0        callr_3.4.3          compiler_4.0.2      
[58] rlang_0.4.7          grid_4.0.2           ggridges_0.5.2      
[61] htmlwidgets_1.5.1    crosstalk_1.1.0.1    igraph_1.2.5        
[64] miniUI_0.1.1.1       base64enc_0.1-3      codetools_0.2-16    
[67] gtable_0.3.0         inline_0.3.15        abind_1.4-5         
[70] curl_4.3             markdown_1.1         reshape2_1.4.4      
[73] R6_2.4.1             gridExtra_2.3        rstantools_2.1.1    
[76] zoo_1.8-8            bridgesampling_1.0-0 dplyr_1.0.0         
[79] fastmap_1.0.1        shinystan_2.5.0      shinythemes_1.1.2   
[82] stringi_1.4.6        parallel_4.0.2       vctrs_0.3.2         
[85] tidyselect_1.1.0     coda_0.19-3

Anyone has problems with the example I posted above?
It’s only me?

Thanks for reporting this. It seems to run fine on my mac with the latest CRAN versions of brms and rstan, but I see you’re on linux so maybe this is OS specific.

Do any other linux users also get a crash here?

Thanks @jonah!
Is there AFAYK a way to have some log? Or, are you so gentle, to paste here the stancode(fit_sc1) produced by brms, so I can try it with rstan or cmdstan and see if I get more hints on what is happening?

Yeah no problem. Below I’ve pasted the Stan code brms generates.

Another thing you can try is specifying brm(..., backend = "cmdstanr") and, if you have CmdStanR installed, it will try running the model using that instead of RStan. (Of course you can also just use the Stan code directly with CmdStanR yourself, but the new backend argument that @paul.buerkner added is pretty nice!)

Here’s the Stan code:

// generated with brms 2.13.3
functions {
  /* cumulative-probit log-PDF for a single response
   * Args:
   *   y: response category
   *   mu: latent mean parameter
   *   disc: discrimination parameter
   *   thres: ordinal thresholds
   * Returns:
   *   a scalar to be added to the log posterior
   real cumulative_probit_lpmf(int y, real mu, real disc, vector thres) {
     int nthres = num_elements(thres);
     real p;
     if (y == 1) {
       p = Phi(disc * (thres[1] - mu));
     } else if (y == nthres + 1) {
       p = 1 - Phi(disc * (thres[nthres] - mu));
     } else {
       p = Phi(disc * (thres[y] - mu)) -
           Phi(disc * (thres[y - 1] - mu));
     return log(p);
  /* cumulative-probit log-PDF for a single response and merged thresholds
   * Args:
   *   y: response category
   *   mu: latent mean parameter
   *   disc: discrimination parameter
   *   thres: vector of merged ordinal thresholds
   *   j: start and end index for the applid threshold within 'thres'
   * Returns:
   *   a scalar to be added to the log posterior
   real cumulative_probit_merged_lpmf(int y, real mu, real disc, vector thres, int[] j) {
     return cumulative_probit_lpmf(y | mu, disc, thres[j[1]:j[2]]);
data {
  int<lower=1> N;  // number of observations
  int Y[N];  // response variable
  int<lower=2> nthres;  // number of thresholds
  int<lower=1> K;  // number of population-level effects
  matrix[N, K] X;  // population-level design matrix
  int prior_only;  // should the likelihood be ignored?
transformed data {
  int Kc = K;
  matrix[N, Kc] Xc;  // centered version of X
  vector[Kc] means_X;  // column means of X before centering
  for (i in 1:K) {
    means_X[i] = mean(X[, i]);
    Xc[, i] = X[, i] - means_X[i];
parameters {
  vector[Kc] b;  // population-level effects
  ordered[nthres] Intercept;  // temporary thresholds for centered predictors
transformed parameters {
  real<lower=0> disc = 1;  // discrimination parameters
model {
  // initialize linear predictor term
  vector[N] mu = Xc * b;
  // priors including all constants
  target += student_t_lpdf(Intercept | 3, 0, 2.5);
  // likelihood including all constants
  if (!prior_only) {
    for (n in 1:N) {
      target += cumulative_probit_lpmf(Y[n] | mu[n], disc, Intercept);
generated quantities {
  // compute actual thresholds
  vector[nthres] b_Intercept = Intercept + dot_product(means_X, b);

@jonah THANKS!.

I place here the solution for those who will have my same problem, or for myself, given the decline of my short-term memory.
The solution is to use cmdstanr as a backend (so the problem, i think is rstan):

data("stemcell", package = "brglm2")
fit_sc1 = brm(formula = research ~ 1 + religion, data = stemcell, family = cumulative("probit"), backend = "cmdstanr")

@jonah, @paul.buerkner if you think it’s a good idea to get our boots dirty and understand the origin of the problem, I’m happy to help.

1 Like

Thanks! I think it’s worth looking into this further but maybe not right away since I think there will be another rstan release soon and maybe that will fix this. If it’s still broken after that then we should investigate this further.

@bgoodri not sure if it’s helpful but tagging you in case it’s useful for you to have an example of rstan crashing R on linux.

Can’t reproduce this. If you can compile it with -O0 -g in your CXX14FLAGS and run R under the debugger

we might be able to figure something out.

I tried again setting the flag in Makevars. No success. I switched to cmdstanr