Error when converting cmdstanr output to stanfit with rstan::read_stan_csv and the latest versions of the packages

[GabrielR]

Hi,

When I try to convert my cmdstanr output to a stanfit object using the latest versions of cmdstanr and rstan, I get the following error:

mod <- cmdstan_model("stancode_longi_gp_v2.stan") 
fit <- mod$sample(
  data = standata,
  chains = 6,
  thin = 1,
  iter_warmup = 10,
  iter_sampling = 15,
  save_warmup = TRUE
)

fit_rstan <- rstan::read_stan_csv(fit$output_files())
# Error in if (max(save_warmup) == 0L) { : 
#   missing value where TRUE/FALSE needed

The .csv files generated by cmdstanr seem to be fine, though I don’t really know much about their contents.

On my old laptop with older versions of the same packages, this error doesn’t occur, and the conversion works.

Any idea what might be causing this?

Many thanks in advance,

Gabriel.

P.S.: I have extensive post-processing code written for rstan, which is why I’d like to convert the object to a stanfit

[jonah]

Welcome to the Stan forum. Unfortunately this is a known issue that hasn’t been fully resolved and may not be fixed. However, there are a few things you can try.

There’s some code that @andrjohns shared that can be used to create the same type of structure for the posterior draws that RStan would give you, which may let you use your post-processing code:

Rstan::read_stan_csv throwing error with cmdstan models (versions 2.35)

If the end goal is a named list where each element is an array of posterior samples for a given variable, you could implement the extract_samples function for cmdstanr and posterior as: extract_samples <- function(fit_obj) { vars <- fit_obj$metadata()$stan_variables draws <- posterior::as_draws_rvars(fit_obj$draws()) lapply(vars, \(var_name){ posterior::draws_of(draws[[var_name]], with_chains = FALSE) }) |> setNames(vars) }

There’s also a guide from @brock that may be helpful:

Rstan::read_stan_csv throwing error with cmdstan models (versions 2.35)

To anyone else hitting this issue with legacy code and just want to run it without rewriting, I posted a guide here: `rstan` fails to read `cmdstanr` output .csv files · Issue #1133 · stan-dev/rstan · GitHub. @andrjohns if anything there seems different from what you would recommend, let me know and I’ll update the guide.

Keeping the CmdStan CSV files fully compatible with RStan is probably not going to be a priority going forward, so I would recommend eventually writing post-processing code compatible with CmdStanR (if you want to use CmdStanR). But maybe all you need is that code from @andrjohns above? That would be the simplest thing to try. If that’s not sufficient and you end up needing to rewrite a lot your code feel free to ask for help here.

[jsocolar]

Not sure it will solve your problem, but see also brms:::read_csv_as_stanfit

[GabrielR]

Hi,

Using brms:::read_csv_as_stanfit() works and creates a stanfit object. Hopefully this function can be kept in future versions (@paul.buerkner). Thanks @jsocolar for the suggestion!

I didn’t check @brock’s guide since brms already offered a simpler solution. @andrjohns’ code is also helpful for extracting posterior samples (compatible with rstan), which can simplify converting post-processing code. Thanks @jonah!

Congrats to the whole team! 🙂

[jonah]

Great, glad that’s working for you!

[Michael_Peck]

After some mostly fruitless digging through documentation and pull requests I managed to rediscover the “draws_of() trick” to get posterior draws in the default form returned by rstan:: extract(). It seems a bit roundabout but it works without any additional packages beyond cmdstanr and posterior. Example:

class(sfit)
[1] “CmdStanMCMC” “CmdStanFit” “R6”
spost ← posterior::as_draws_rvars(sfit$draws())
spost_as_array ← lapply(spost, posterior::draws_of)
names(spost)
[1] “lp__” “kernel” “beta_st” “lprof_em” “beta_em” “beta0”
[7] “tauv” “delta” “mu_y” “spec_em” “gfit” “gflux_rep”
[13] “log_lik” “ll”
names(spost_as_array)
[1] “lp__” “kernel” “beta_st” “lprof_em” “beta_em” “beta0”
[7] “tauv” “delta” “mu_y” “spec_em” “gfit” “gflux_rep”
[13] “log_lik” “ll”
spost$kernel
rvar<750,4>[11] mean ± sd:
[1] 0.0027 ± 0.0046 0.0047 ± 0.0071 0.0069 ± 0.0096 0.0082 ± 0.0123
[5] 0.2892 ± 0.0268 0.3706 ± 0.0344 0.2687 ± 0.0338 0.0225 ± 0.0232
[9] 0.0089 ± 0.0119 0.0091 ± 0.0116 0.0085 ± 0.0105
dim(spost_as_array$kernel)
[1] 3000 11

Lost the link but somebody posted a more tidyverse like solution. A simple application of lapply() works fine for my purposes.

Of course the link is in the second post by @Jonah.