Stan model runs fine in R terminal but errors when compiled in R markdown document

Interfaces RStan
,
1

Hi all,

I’m trying to fit a simple Bayesian vector autoregression model in Rstan and have no problem when I run the code chunks one by one interactively in Rstudio, but when I try to compile the file I get the error

Line 863: Error: Stan model VARModel does not contain samples. In addition: There were 22 warnings (use warnings() to see them) Execution halted.

The code leading up to this is

stan_code_var<-"
data {
  int<lower=1> P;
  int<lower=1> K;
  int<lower=0> N;
  vector[K] y[N];
}
transformed data {
  vector[K] alphamean = rep_vector(1.0, K);
  cov_matrix[K] Id = diag_matrix(rep_vector(1.0, K));
  int J = K * K;
  cov_matrix[J] BetaCov = diag_matrix(rep_vector(1.0, J));
  matrix[K,K] Betamean[P]; 
  for (p in 1:P)
    Betamean[p]=rep_matrix(0.0,K,K);
  Betamean[1,1,1] = 1.0;
  Betamean[1,2,2] = 1.0;
}
parameters {
  vector[K] alpha;
  matrix[K,K] Beta[P];
  cov_matrix[K] Sigma;
}
model {
  Sigma ~ inv_wishart(5,2.0*Id);
  alpha ~ multi_normal(alphamean, 2.0*Id);
  for (p in 1:P) {
    to_vector(Beta[p]) ~ multi_normal(to_vector(Betamean[p]),BetaCov);
  }
  for (n in (P+1):N) {
    vector[K] mu = alpha;
    for (p in 1:P){
      mu += Beta[p] * y[n-p];
    }
    y[n] ~ multi_normal(mu, Sigma);
  }
}"

stan_model_var<-stan_model(model_code = stan_code_var,model_name = "VARModel") 

options(mc.cores = parallel::detectCores()) #Use parallel computing when running MCMC

K = length(bvardata[1,])
N = length(bvardata[,1])
lags = 1 
stan_data<-list(P=lags, K = K, N = N, y=bvardata)
fit_var1<-sampling(object = stan_model_var,data = stan_data, chains = 4, iter = 2000, seed = 5678)

print(fit_var1)

plot(fit_var1)

and finally line 863 is

quietgg(stan_hist(fit_var1))

That is, it fails right at the plotting step. When I run these lines interactively, the model samples just fine, produces chains and plots histograms with no errors and the specified number of samples, but the error occurs when I recompile the full Rmd file. Removing the quietgg doesn’t change this result.

As you might have guessed from the line 863, there’s a lot of stuff going on in the code before that, mostly data cleaning stuff to go from my raw data to the bvardata object, but also a bunch of other preliminary analyses using various packages, so one or more of those may have interfered. Among these, I’m running a bunch of Julia code in the notebook through the JuliaCall library. I could also just be an idiot and running the cells in the wrong order but I’ve tried from the beginning several times.

On the model itself, I know it’s terrible (I’m trying to replicate existing analyses before moving on to more sensible modeling choices) but don’t expect that that’s the source of the problem since there are no sampling issues when run outside of the knitr compilation.

Any hints as to possible things to explore to resolve this?

For reference, I am using rstan 2.19 on R 3.5.0 in Rstudio 1.2.5033

devtools::session_info(“rstan”) produces

─ Session info ────────────────────────────────────────────────────────────────────────────────────────────────────────────
 setting  value                       
 version  R version 3.5.0 (2018-04-23)
 os       macOS  10.15.6              
 system   x86_64, darwin15.6.0        
 ui       RStudio                     
 language (EN)                        
 collate  en_US.UTF-8                 
 ctype    en_US.UTF-8                 
 tz       America/New_York            
 date     2020-09-04                  

─ Packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 package      * version   date       lib source        
 assertthat     0.2.1     2019-03-21 [1] CRAN (R 3.5.2)
 backports      1.1.6     2020-04-05 [1] CRAN (R 3.5.0)
 BH             1.72.0-3  2020-01-08 [1] CRAN (R 3.5.2)
 callr          3.4.3     2020-03-28 [1] CRAN (R 3.5.0)
 checkmate      2.0.0     2020-02-06 [1] CRAN (R 3.5.2)
 cli            2.0.2     2020-02-28 [1] CRAN (R 3.5.2)
 colorspace     1.4-1     2019-03-18 [1] CRAN (R 3.5.2)
 crayon         1.3.4     2017-09-16 [1] CRAN (R 3.5.0)
 desc           1.2.0     2018-05-01 [1] CRAN (R 3.5.0)
 digest         0.6.25    2020-02-23 [1] CRAN (R 3.5.2)
 ellipsis       0.3.1     2020-05-15 [1] CRAN (R 3.5.0)
 evaluate       0.14      2019-05-28 [1] CRAN (R 3.5.2)
 fansi          0.4.1     2020-01-08 [1] CRAN (R 3.5.2)
 farver         2.0.3     2020-01-16 [1] CRAN (R 3.5.2)
 ggplot2      * 3.3.0     2020-03-05 [1] CRAN (R 3.5.2)
 glue           1.4.1     2020-05-13 [1] CRAN (R 3.5.0)
 gridExtra    * 2.3       2017-09-09 [1] CRAN (R 3.5.0)
 gtable         0.3.0     2019-03-25 [1] CRAN (R 3.5.2)
 inline         0.3.15    2018-05-18 [1] CRAN (R 3.5.0)
 isoband        0.2.1     2020-04-12 [1] CRAN (R 3.5.0)
 labeling       0.3       2014-08-23 [1] CRAN (R 3.5.0)
 lattice        0.20-38   2018-11-04 [1] CRAN (R 3.5.0)
 lifecycle      0.2.0     2020-03-06 [1] CRAN (R 3.5.2)
 loo            2.2.0     2019-12-19 [1] CRAN (R 3.5.2)
 magrittr       1.5       2014-11-22 [1] CRAN (R 3.5.0)
 MASS           7.3-51.4  2019-03-31 [1] CRAN (R 3.5.2)
 Matrix         1.2-17    2019-03-22 [1] CRAN (R 3.5.2)
 matrixStats    0.56.0    2020-03-13 [1] CRAN (R 3.5.2)
 mgcv           1.8-28    2019-03-21 [1] CRAN (R 3.5.2)
 munsell        0.5.0     2018-06-12 [1] CRAN (R 3.5.0)
 nlme           3.1-141   2019-08-01 [1] CRAN (R 3.5.2)
 pillar         1.4.4     2020-05-05 [1] CRAN (R 3.5.0)
 pkgbuild       1.0.6     2019-10-09 [1] CRAN (R 3.5.2)
 pkgconfig      2.0.3     2019-09-22 [1] CRAN (R 3.5.2)
 pkgload        1.0.2     2018-10-29 [1] CRAN (R 3.5.0)
 praise         1.0.0     2015-08-11 [1] CRAN (R 3.5.0)
 prettyunits    1.1.1     2020-01-24 [1] CRAN (R 3.5.2)
 processx       3.4.2     2020-02-09 [1] CRAN (R 3.5.2)
 ps             1.3.2     2020-02-13 [1] CRAN (R 3.5.2)
 R6             2.4.1     2019-11-12 [1] CRAN (R 3.5.2)
 RColorBrewer   1.1-2     2014-12-07 [1] CRAN (R 3.5.0)
 Rcpp         * 1.0.5     2020-07-06 [1] CRAN (R 3.5.0)
 RcppEigen      0.3.3.7.0 2019-11-16 [1] CRAN (R 3.5.2)
 rlang          0.4.6     2020-05-02 [1] CRAN (R 3.5.0)
 rprojroot      1.3-2     2018-01-03 [1] CRAN (R 3.5.0)
 rstan        * 2.19.3    2020-02-11 [1] CRAN (R 3.5.2)
 rstudioapi     0.11      2020-02-07 [1] CRAN (R 3.5.2)
 scales         1.1.0     2019-11-18 [1] CRAN (R 3.5.2)
 StanHeaders  * 2.21.0-1  2020-01-19 [1] CRAN (R 3.5.2)
 testthat       2.3.2     2020-03-02 [1] CRAN (R 3.5.2)
 tibble       * 3.0.1     2020-04-20 [1] CRAN (R 3.5.0)
 utf8           1.1.4     2018-05-24 [1] CRAN (R 3.5.0)
 vctrs          0.3.1     2020-06-05 [1] CRAN (R 3.5.0)
 viridisLite    0.3.0     2018-02-01 [1] CRAN (R 3.5.0)
 withr          2.1.2     2018-03-15 [1] CRAN (R 3.5.0)
2

Hmmm, I have no idea why it would be fine interactively but crash when compiling the R markdown document, that’s strange!

If it’s really the plotting step that’s the issue (I doubt it but I guess it’s possible), then you could use our bayesplot package instead of the plotting functions in RStan (we’re planning to transition to that anyway in the future).

Instead of stan_hist() you’d use bayesplot’s mcmc_hist(). And instead of the default plot(fit_var1), which makes an interval plot if I recall correctly, you can use bayesplot::mcmc_intervals(). Both make ggplot objects.

In theory it shouldn’t be crashing regardless of these versions, but it’s conceivable that just updating these may help. There are newer versions of rstan, R, and RStudio.

3

I tried using bayes_plot instead of rstan plotting and while it again has no problems in R now the error message is

Error in prepare_mcmc_array(x, pars, regex_pars, transformations) : is.matrix(x) || is.array(x) is not TRUE Calls: <Anonymous> ... mcmc_hist -> .mcmc_hist -> prepare_mcmc_array -> stopifnot

which is at least telling me that it doesn’t seem to even recognize the fit object, not just that it is missing samples.

I will try the upgrades and see if any of those help. That may take me a bit longer.

4

Which version of bayesplot is that? You may need to pass it as.matrix(stanfit) instead of your stanfit object itself, although I think we allow the latter in the latest version of bayesplot.

5

Bayesplot 1.7.0 which I see is not the latest, but again, the same command (mcmc_hist(fit_var1)) runs interactively but not when compiling the R markdown file. I’ll try updating that too, though, thanks.