warmup draws omitted from stanfit objects

[wds15]

Right now, the conversion of stanfit objects into the draws format relies on the as.array function from Rstan. It would be great to document this somewhere and I would actually expect the function to convert all draws in the stanfit object, but the function drops the warmup samples without a way to influence this as the following example demonstrates:

library(rstan)
library(posterior)

ex_model_code <- '
  parameters {
    real alpha[2,3];
    real beta[2]; 
  } 
  model {
    for (i in 1:2) for (j in 1:3) 
      alpha[i, j] ~ normal(0, 1); 
    for (i in 1:2) 
      beta ~ normal(0, 2); 
  } 
'

## fit the model 
fit <- stan(model_code = ex_model_code, chains = 4) 

## we sampled and kept the warmup samples
fit

draws_post <- as_draws(fit)

## but converting to draws will omit the warmup draws! There is no way
## to get the warmup draws, since "as.array" is called. The as.array
## call will by default use "extract(x, permuted = FALSE, inc_warmup =
## FALSE)" such that warmup is always dropped.
draws_post

## I think, the default should be to return the warmup draws like so:
draws_post_alt <- as_draws(rstan::extract(fit, permuted = FALSE, inc_warmup = TRUE))

draws_post_alt

[avehtari]

This should be an issue in rstan repository. posterior package is agnostic to where the draws come fro. For example, cmdstanr includes definitions of as_draws for cmdstanr objects. posterior::as_draws uses as.array as the fallback for any R object, hoping that the related packages have that implemented, but it's not the role of posterior package to know more about the foreign objects. The correct action would be to add as_draws.stanfit function to rstan package.

(You are using a super old rstan version as that example has not worked in Stan for years)

[wds15]

So, yes, it is a thing that rstan does not return the warmup samples with the as.array function. Maybe @jgabry can weigh in here?

Given that posterior is within the Stan universe of things, it would be very desirable to have the matter fixed on either end - rstan or posterior. Conceptually, I agree in that that as.array from rstan should give back the warmup samples; not sure what was the consideration there at the time. As a user I find it very inconsistent to not get the warmup samples from as_draws even though I requested them to be sampled.

[avehtari]

As a user it is fine that you reported the issue here, you don't need to understand the modularity of the packages. I should have been more clear that my response was not to criticize you for your choice of creating the issue here. Thanks for the issue, and we package developers fix it where it is appropriate.

[wds15]

All good. Nothing to worry about. I did not feel offended at all; just wanted to make the Stan universe a bit better as usual. I did not have the design choices in mind when writing the issue.

I recall that I ran into this matter when implementing things for one of my packages:

https://github.com/Novartis/OncoBayes2/blob/0b8809c6aa7b363c16ff404615fe9108f5d41039/R/posterior.R#L201

brms is also affected by this problem which follows from this... so warmup draws are not part of the posterior draws object even if one stores the warmup samples... tagging @paul-buerkner just as FYI. The implementation in OncoBayes2 is what I made of what I would expect as a user.

[jgabry]

I'm not sure what the historical reason was, but as.array.stanfit doesn't expose an inc_warmup argument. It just calls extract internally and sets inc_warmup = FALSE:

> rstan:::as.array.stanfit
function (x, ...) 
{
    if (x@mode != 0) 
        return(numeric(0))
    out <- extract(x, permuted = FALSE, inc_warmup = FALSE, ...)
    return(out)
}

I think it makes sense to default to not including the warmup draws, but there should be a way to get them and I'm not sure why the inc_warmup argument wasn't exposed.

If as_draws() passes ... to as.array() (need to double check that) then we could just expose the inc_warmup argument to as.array and then it would be possible to do as_draws(stanfit, inc_warmup=TRUE). But we should probably just add an as_draws.stanfit method to rstan, like @avehtari said.

[jgabry]

But we should probably just add an as_draws.stanfit method to rstan, like @avehtari said.

I'm not sure when I'll have time to add this (I have a lot of work on my plate at the moment and am recovering from being sick), but it would be pretty easy for anyone to add.