Mixture model fails on brms using cmdstanr, but not using cmdstanr directly

[fusaroli]

I am playing with the data in https://doi.org/10.1371/journal.pcbi.1005833 2 , where different kinds of primates play a “matching pennies” game (guess in which hand is the tasty treat).

I define a handful of possible strategies and want to build a mixture model trying to identify which strategy each primate is mostly likely using.

The two simplest ones (as an example here):

Strategy 1: Simple Bias: the primate just picks the left hand with a fixed rate. E.g. the primate is left handed and tends to pick left 80% of the time.
Strategy 2: Biased Nash: the primate keeps track of how many times the food was in the left hand so far, and picks accordingly (more likely to pick left, the higher the proportion of times food was in the left hand)

I build the models on a single primate and the single strategy models work nicely.
The mixture model fails when using "cmdstanr" as backend, but not when using rstan as backend, or when using cmdstanr directly. The problems only occur when I specify a formula for theta (the probability of picking strategy 1), which I'll need to do when scaling up to multilevel modeling.

SampleData1.csv

d <- readr::read_csv("SampleData1.csv")

mix <- mixture(bernoulli, bernoulli, order = "none")

mix_f0 <- bf(
  left ~ 1 ,
  mu1 ~ 1,
  mu2 ~ 1 + LeftProb
)

mix_f1 <- bf(
  left ~ 1 ,
  mu1 ~ 1,
  mu2 ~ 1 + LeftProb,
  theta1 ~ 1
)

# defining priors
mix_p1 <- c(
  prior(normal(0,0.3), class=Intercept, dpar=mu1),
  prior(normal(0,0.3), class=Intercept, dpar=mu2),
  prior(normal(0,0.3), class=b, dpar=mu2)
)

# Fitting the model without theta works
mix_m0 <- brm(
  mix_f0,
  d,
  family = mix,
  mix_p1,
  sample_prior = T,
  chains = 1,
  cores = 1,
  backend = "cmdstanr",
  threads = threading(2),
  file = paste0("mix_m0"),
  file_refit = "on_change",
  control = list(
    adapt_delta = 0.99,
    max_treedepth = 20
  )
)

# Fitting the model with theta doesn't work (using brms w cmdstanr backend)
mix_m1 <- brm(
  mix_f1,
  d,
  family = mix,
  mix_p1,
  sample_prior = T,
  chains = 1,
  cores = 1,
  backend = "cmdstanr",
  threads = threading(2),
  file = paste0("mix_m1"),
  file_refit = "on_change",
  control = list(
    adapt_delta = 0.99,
    max_treedepth = 20
  )
)

The error I get (both with and without threading) is

Compiling Stan program…
Start sampling
Running MCMC with 1 chain, with 2 thread(s) per chain…

Chain 1 num_threads = 2
Warning: Chain 1 finished unexpectedly!

Error in rstan::read_stan_csv(out$output_files()) :
csvfiles does not contain any CSV file name
In addition: Warning message:
No chains finished successfully. Unable to retrieve the fit.

What is curious is that the error doesn't happen with rstan backend:

mix_m1 <- brm(
  mix_f1,
  d,
  family = mix,
  mix_p1,
  sample_prior = T,
  chains = 1,
  cores = 1,
  #backend = "cmdstanr",
  #threads = threading(2),
  file = paste0("mix_m1"),
  file_refit = "on_change",
  control = list(
    adapt_delta = 0.99,
    max_treedepth = 20
  )
)

Nor if I fit it directly via cmdstanr

scode <- make_stancode(mix_f1,
                       d,
                       family = mix,
                       mix_p1)
sdata <- make_standata(mix_f1,
                       d,
                       family = mix,
                       mix_p1)

m <- cmdstan_model(write_stan_file(scode))

class(sdata) <- NULL # Necessary for stupid reasons
fit <- m$sample(data = as.list(sdata))
fit$output()

[fusaroli]

See also https://discourse.mc-stan.org/t/mixture-model-fails-on-cmdstanr-when-conditioning-theta-on-participant-but-not-on-rstan/22018 (where Martin Modrak confirms he can reproduce the error)

[paul-buerkner]

Interesting. Will check it out.

[fusaroli]

thanks! (I cleaned up the code above from a copy/paste error, I had forgotten commenting away the backend lines for the "the error doesn't happen with rstan backend")

[paul-buerkner]

Found and fixed it. The problem did only occur during threading because one indexing of theta1 was wrong in that case.