Error running the same model twice: cmdstan with ulam function (Rethinking package)

Modeling
1

Thanks in advance for any help with this hopefully minor issue!

I am trying to fit models from McElreath’s Statistical Rethinking book, using ‘ulam’ with cmdstan. The issue seems to be that when I run the same exact model code twice in the same R session (i.e. run the same compiled RStan program) something hangs and the chains finish unexpectedly. However, if you run essentially the same model but with different parameter names (i.e. change the program but not the model structure), it works fine again.

The problem is that this happens even if you run the same model structure with different numbers of chains, etc., so this would be very annoying if I can’t figure it out!

I expect this is some minor settings issue. Is there an easy fix to this?

Running:
R v 4.1.1 on Mac ARM/M1
Rstan v 2.21.2
Cmdstanr v 0.4.0
Rethinking v 2.13

To illustrate with code (data attached):

library(“rstan”)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
library(cmdstanr)
library(posterior)
library(bayesplot)
color_scheme_set(“brightblue”)
library(rethinking)

Model formula/code, which works fine:

m9.1 ← ulam(
alist(
log_gdp_std ~ dnorm(mu, sigma),
mu ← a[cid] + b[cid]*( rugged_std - 0.215 ),
a[cid] ~ dnorm(1, 0.1),
b[cid] ~ dnorm(0, 0.3),
sigma ~ dexp( 1 )
), data=dd_slim, chains=1, cmdstan=TRUE )
precis(m9.1, depth=2)

It runs fine, and this is the output that I get:

m9.1 ← ulam(

  • alist(
  • log_gdp_std ~ dnorm(mu, sigma),
  • mu ← a[cid] + b[cid]*( rugged_std - 0.215 ),
  • a[cid] ~ dnorm(1, 0.1),
  • b[cid] ~ dnorm(0, 0.3),
  • sigma ~ dexp( 1 )
  • ), data=dd_slim, chains=1, cmdstan=TRUE )
    Compiling Stan program…
    Running MCMC with 1 chain, with 1 thread(s) per chain…

Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 1 Iteration: 1000 / 1000 [100%] (Sampling)
Chain 1 finished in 0.1 seconds.

precis(m9.1, depth=2)
mean sd 5.5% 94.5% n_eff Rhat4
a[1] 0.89 0.01 0.86 0.91 610 1
a[2] 1.05 0.01 1.04 1.07 710 1
b[1] 0.13 0.08 0.01 0.25 777 1
b[2] -0.14 0.06 -0.23 -0.05 712 1
sigma 0.11 0.01 0.10 0.12 734 1

However, if I run a different model with the same formula, it doesn’t work:

m9.2 ← ulam(
alist(
log_gdp_std ~ dnorm(mu, sigma),
mu ← a[cid] + b[cid]*( rugged_std - 0.215 ),
a[cid] ~ dnorm(1, 0.1),
b[cid] ~ dnorm(0, 0.3),
sigma ~ dexp( 1 )
), data=dd_slim, chains=1, cmdstan=TRUE )
precis(m9.2, depth=2)

Output:

m9.2 ← ulam(

  • alist(
  • log_gdp_std ~ dnorm(mu, sigma),
  • mu ← a[cid] + b[cid]*( rugged_std - 0.215 ),
  • a[cid] ~ dnorm(1, 0.1),
  • b[cid] ~ dnorm(0, 0.3),
  • sigma ~ dexp( 1 )
  • ), data=dd_slim, chains=1, cmdstan=TRUE )
    Compiling Stan program…
    Running MCMC with 1 chain, with 1 thread(s) per chain…
    Warning: Chain 1 finished unexpectedly!
    Error in rstan::read_stan_csv(cmdstanfit$output_files()) :
    csvfiles does not contain any CSV file name
    In addition: Warning message:
    No chains finished successfully. Unable to retrieve the fit.

But, now, if I make a third model that just changes the parameter names, it works fine again:

m9.3 ← ulam(
alist(
log_gdp_std ~ dnorm(mu, sigma),
mu ← x[cid] + y[cid]*( rugged_std - 0.215 ),
x[cid] ~ dnorm(1, 0.1),
y[cid] ~ dnorm(0, 0.3),
sigma ~ dexp( 1 )
), data=dd_slim, chains=1, cmdstan=TRUE )
precis(m9.3, depth=2)

Output:

m9.3 ← ulam(

  • alist(
  • log_gdp_std ~ dnorm(mu, sigma),
  • mu ← x[cid] + y[cid]*( rugged_std - 0.215 ),
  • x[cid] ~ dnorm(1, 0.1),
  • y[cid] ~ dnorm(0, 0.3),
  • sigma ~ dexp( 1 )
  • ), data=dd_slim, chains=1, cmdstan=TRUE )
    Compiling Stan program…
    Running MCMC with 1 chain, with 1 thread(s) per chain…

Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 1 Iteration: 1000 / 1000 [100%] (Sampling)
Chain 1 finished in 0.1 seconds.

precis(m9.3, depth=2)
mean sd 5.5% 94.5% n_eff Rhat4
x[1] 0.89 0.02 0.86 0.91 962 1
x[2] 1.05 0.01 1.03 1.07 722 1
y[1] 0.14 0.08 0.01 0.25 615 1
y[2] -0.14 0.05 -0.23 -0.06 899 1
sigma 0.11 0.01 0.10 0.12 1208 1

This must be some setting issue somewhere

2

Have you checked Richard’s GitHub issues?

If you can’t find anything there then perhaps add an issue?

3

I did check there, thanks! I thought it would be more of a Stan issue, but it just occurred to me that the problem could be at the interface of rethinking/Stan (which seems to be at the root of others’ issues).

And, lo and behold the Stan code on it’s own works fine, even if the ulam code doesn’t. Weird, but I guess that ‘solves’ it.

4

Hmm, if that’s the case then it sounds as if rethinking package might be the problem here. Could also be some interaction with RStudio perhaps? Anyways, I’d file an issue so Richard is notified about it.

5

Definitely not Rstudio because I use terminal. Only other possibility is that it has to do with it being a Mac with Apple silicon/M1. I’ll let Richard know, though!

1 Like
6

Hi @baileyhouse, did you manage to resolve this issue? I couldn’t track down an issue on the rethinking github, did you post there in the end? I’m experiencing the same issue on a new M1 mac.

7

Seems to run fine a second time with force_recompile = TRUE. So with rethinking:

model1 <- ulam(
   flist = flist,
   data = data,
   etc...
   control = list(adapt_delta = .95, force_recompile = TRUE)
)

Note that you have to specify adapt_delta again, as a new control list clears the ulam default.