read_stan_csv parameter naming scheme differs from sampling output

[yonicd]

Summary:

The naming scheme used to reconstruct the stanfit object in read_stan_csv is different than the original rstan::sampling output

shredder::rat_data

$N
[1] 30

$T
[1] 5

$x
[1]  8 15 22 29 36

$y
# A tibble: 30 x 5
    day8 day15 day22 day29 day36
   <int> <int> <int> <int> <int>
 1   151   199   246   283   320
 2   145   199   249   293   354
 3   147   214   263   312   328
 4   155   200   237   272   297
 5   135   188   230   280   323
 6   159   210   252   298   331
 7   141   189   231   275   305
 8   159   201   248   297   338
 9   177   236   285   350   376
10   134   182   220   260   296
11   160   208   261   313   352
12   143   188   220   273   314
13   154   200   244   289   325
14   171   221   270   326   358
15   163   216   242   281   312
16   160   207   248   288   324
17   142   187   234   280   316
18   156   203   243   283   317
19   157   212   259   307   336
20   152   203   246   286   321
21   154   205   253   298   334
22   139   190   225   267   302
23   146   191   229   272   302
24   157   211   250   285   323
25   132   185   237   286   331
26   160   207   257   303   345
27   169   216   261   295   333
28   157   205   248   289   316
29   137   180   219   258   291
30   153   200   244   286   324

$xbar
[1] 22

system.file('rats.stan',package='shredder')

data {
  int<lower=0> N;
  int<lower=0> T;
  real x[T];
  real y[N,T];
  real xbar;
}
parameters {
  real alpha[N];
  real beta[N];

  real mu_alpha;
  real mu_beta;          // beta.c in original bugs model

  real<lower=0> sigmasq_y;
  real<lower=0> sigmasq_alpha;
  real<lower=0> sigmasq_beta;
}
transformed parameters {
  real<lower=0> sigma_y;       // sigma in original bugs model
  real<lower=0> sigma_alpha;
  real<lower=0> sigma_beta;

  sigma_y = sqrt(sigmasq_y);
  sigma_alpha = sqrt(sigmasq_alpha);
  sigma_beta = sqrt(sigmasq_beta);
}
model {
  mu_alpha ~ normal(0, 100);
  mu_beta ~ normal(0, 100);
  sigmasq_y ~ inv_gamma(0.001, 0.001);
  sigmasq_alpha ~ inv_gamma(0.001, 0.001);
  sigmasq_beta ~ inv_gamma(0.001, 0.001);
  alpha ~ normal(mu_alpha, sigma_alpha); // vectorized
  beta ~ normal(mu_beta, sigma_beta);  // vectorized
  for (n in 1:N)
    for (t in 1:T) 
      y[n,t] ~ normal(alpha[n] + beta[n] * (x[t] - xbar), sigma_y);

}
generated quantities {
  real alpha0;
  alpha0 = mu_alpha - xbar * mu_beta;
}

library(shredder)
library(rstan)
#> Loading required package: StanHeaders
#> Loading required package: ggplot2
#> rstan (Version 2.19.2, GitRev: 2e1f913d3ca3)
#> For execution on a local, multicore CPU with excess RAM we recommend calling
#> options(mc.cores = parallel::detectCores()).
#> To avoid recompilation of unchanged Stan programs, we recommend calling
#> rstan_options(auto_write = TRUE)

junk <- capture.output({ #verbose doesnt work in reprex
  fit <- rstan::stan(
    file  = system.file('rats.stan',package='shredder'),
    data  = shredder::rat_data,verbose = FALSE
  )  
})
#> Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
#> http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: The largest R-hat is 1.05, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#tail-ess

m <- stan_model(system.file('rats.stan',package='shredder'))

td <- file.path(tempdir(),'test')
dir.create(td)

junk <- capture.output({ #verbose doesnt work in reprex
f <- sampling(m, data = shredder::rat_data,sample_file = file.path(td,'samp'),verbose = FALSE)
})

from_csv <- rstan::read_stan_csv(list.files(td,full.names = TRUE))

class(from_csv)
#> [1] "stanfit"
#> attr(,"package")
#> [1] "rstan"

class(fit)
#> [1] "stanfit"
#> attr(,"package")
#> [1] "rstan"

names(from_csv@sim$samples[[1]])
#>  [1] "alpha.1"       "alpha.2"       "alpha.3"       "alpha.4"      
#>  [5] "alpha.5"       "alpha.6"       "alpha.7"       "alpha.8"      
#>  [9] "alpha.9"       "alpha.10"      "alpha.11"      "alpha.12"     
#> [13] "alpha.13"      "alpha.14"      "alpha.15"      "alpha.16"     
#> [17] "alpha.17"      "alpha.18"      "alpha.19"      "alpha.20"     
#> [21] "alpha.21"      "alpha.22"      "alpha.23"      "alpha.24"     
#> [25] "alpha.25"      "alpha.26"      "alpha.27"      "alpha.28"     
#> [29] "alpha.29"      "alpha.30"      "beta.1"        "beta.2"       
#> [33] "beta.3"        "beta.4"        "beta.5"        "beta.6"       
#> [37] "beta.7"        "beta.8"        "beta.9"        "beta.10"      
#> [41] "beta.11"       "beta.12"       "beta.13"       "beta.14"      
#> [45] "beta.15"       "beta.16"       "beta.17"       "beta.18"      
#> [49] "beta.19"       "beta.20"       "beta.21"       "beta.22"      
#> [53] "beta.23"       "beta.24"       "beta.25"       "beta.26"      
#> [57] "beta.27"       "beta.28"       "beta.29"       "beta.30"      
#> [61] "mu_alpha"      "mu_beta"       "sigmasq_y"     "sigmasq_alpha"
#> [65] "sigmasq_beta"  "sigma_y"       "sigma_alpha"   "sigma_beta"   
#> [69] "alpha0"        "lp__"

names(fit@sim$samples[[1]])
#>  [1] "alpha[1]"      "alpha[2]"      "alpha[3]"      "alpha[4]"     
#>  [5] "alpha[5]"      "alpha[6]"      "alpha[7]"      "alpha[8]"     
#>  [9] "alpha[9]"      "alpha[10]"     "alpha[11]"     "alpha[12]"    
#> [13] "alpha[13]"     "alpha[14]"     "alpha[15]"     "alpha[16]"    
#> [17] "alpha[17]"     "alpha[18]"     "alpha[19]"     "alpha[20]"    
#> [21] "alpha[21]"     "alpha[22]"     "alpha[23]"     "alpha[24]"    
#> [25] "alpha[25]"     "alpha[26]"     "alpha[27]"     "alpha[28]"    
#> [29] "alpha[29]"     "alpha[30]"     "beta[1]"       "beta[2]"      
#> [33] "beta[3]"       "beta[4]"       "beta[5]"       "beta[6]"      
#> [37] "beta[7]"       "beta[8]"       "beta[9]"       "beta[10]"     
#> [41] "beta[11]"      "beta[12]"      "beta[13]"      "beta[14]"     
#> [45] "beta[15]"      "beta[16]"      "beta[17]"      "beta[18]"     
#> [49] "beta[19]"      "beta[20]"      "beta[21]"      "beta[22]"     
#> [53] "beta[23]"      "beta[24]"      "beta[25]"      "beta[26]"     
#> [57] "beta[27]"      "beta[28]"      "beta[29]"      "beta[30]"     
#> [61] "mu_alpha"      "mu_beta"       "sigmasq_y"     "sigmasq_alpha"
#> [65] "sigmasq_beta"  "sigma_y"       "sigma_alpha"   "sigma_beta"   
#> [69] "alpha0"        "lp__"

R.version.string
#> [1] "R version 3.6.1 (2019-07-05)"

packageVersion("rstan")
#> [1] '2.19.2'

details::details(sessioninfo::session_info(),summary = 'R Session Info')

R Session Info

─ Session info ──────────────────────────────────────────────────────────
 setting  value                       
 version  R version 3.6.1 (2019-07-05)
 os       macOS Mojave 10.14.5        
 system   x86_64, darwin15.6.0        
 ui       X11                         
 language (EN)                        
 collate  en_US.UTF-8                 
 ctype    en_US.UTF-8                 
 tz       America/New_York            
 date     2020-01-17                  

─ Packages ──────────────────────────────────────────────────────────────
 package     * version    date       lib source                         
 assertthat    0.2.1      2019-03-21 [1] CRAN (R 3.6.0)                 
 backports     1.1.5      2019-10-02 [1] CRAN (R 3.6.0)                 
 callr         3.3.2      2019-09-22 [1] CRAN (R 3.6.0)                 
 checkmate     1.9.4      2019-07-04 [1] CRAN (R 3.6.0)                 
 cli           2.0.1      2020-01-08 [1] CRAN (R 3.6.0)                 
 clipr         0.7.0      2019-07-23 [1] CRAN (R 3.6.0)                 
 codetools     0.2-16     2018-12-24 [1] CRAN (R 3.6.1)                 
 colorspace    1.4-1      2019-03-18 [1] CRAN (R 3.6.0)                 
 crayon        1.3.4      2017-09-16 [1] CRAN (R 3.6.0)                 
 desc          1.2.0      2019-12-01 [1] Github (r-lib/desc@61205f6)    
 details       0.2.1      2020-01-12 [1] local                          
 digest        0.6.23     2019-11-23 [1] CRAN (R 3.6.0)                 
 dplyr         0.8.3      2019-07-04 [1] CRAN (R 3.6.0)                 
 evaluate      0.14       2019-05-28 [1] CRAN (R 3.6.0)                 
 fansi         0.4.1      2020-01-08 [1] CRAN (R 3.6.0)                 
 ggplot2     * 3.2.1      2019-08-10 [1] CRAN (R 3.6.0)                 
 glue          1.3.1.9000 2020-01-07 [1] Github (tidyverse/glue@b9ffe6c)
 gridExtra     2.3        2017-09-09 [1] CRAN (R 3.6.0)                 
 gtable        0.3.0      2019-03-25 [1] CRAN (R 3.6.0)                 
 highr         0.8        2019-03-20 [1] CRAN (R 3.6.0)                 
 htmltools     0.4.0      2019-10-04 [1] CRAN (R 3.6.0)                 
 httr          1.4.1      2019-08-05 [1] CRAN (R 3.6.0)                 
 inline        0.3.15     2018-05-18 [1] CRAN (R 3.6.0)                 
 knitr         1.25       2019-09-18 [1] CRAN (R 3.6.0)                 
 lazyeval      0.2.2      2019-03-15 [1] CRAN (R 3.6.0)                 
 lifecycle     0.1.0      2019-08-01 [1] CRAN (R 3.6.0)                 
 loo           2.1.0      2019-03-13 [1] CRAN (R 3.6.0)                 
 magrittr      1.5        2014-11-22 [1] CRAN (R 3.6.0)                 
 matrixStats   0.54.0     2018-07-23 [1] CRAN (R 3.6.0)                 
 munsell       0.5.0      2018-06-12 [1] CRAN (R 3.6.0)                 
 pillar        1.4.3      2019-12-20 [1] CRAN (R 3.6.0)                 
 pkgbuild      1.0.6      2019-10-09 [1] CRAN (R 3.6.1)                 
 pkgconfig     2.0.3      2019-09-22 [1] CRAN (R 3.6.0)                 
 png           0.1-7      2013-12-03 [1] CRAN (R 3.6.0)                 
 prettyunits   1.0.2      2015-07-13 [1] CRAN (R 3.6.0)                 
 processx      3.4.1      2019-07-18 [1] CRAN (R 3.6.0)                 
 ps            1.3.0      2018-12-21 [1] CRAN (R 3.6.0)                 
 purrr         0.3.3      2019-10-18 [1] CRAN (R 3.6.0)                 
 R6            2.4.1      2019-11-12 [1] CRAN (R 3.6.0)                 
 Rcpp          1.0.3      2019-11-08 [1] CRAN (R 3.6.1)                 
 rematch2      2.1.0      2019-07-11 [1] CRAN (R 3.6.0)                 
 rlang         0.4.2      2019-11-23 [1] CRAN (R 3.6.0)                 
 rmarkdown     2.0        2019-12-12 [1] CRAN (R 3.6.0)                 
 rprojroot     1.3-2      2018-01-03 [1] CRAN (R 3.6.0)                 
 rstan       * 2.19.2     2019-07-09 [1] CRAN (R 3.6.0)                 
 scales        1.1.0      2019-11-18 [1] CRAN (R 3.6.0)                 
 sessioninfo   1.1.1      2018-11-05 [1] CRAN (R 3.6.0)                 
 shredder    * 0.1.1      2019-12-20 [1] local                          
 StanHeaders * 2.18.1-10  2019-06-14 [1] CRAN (R 3.6.0)                 
 stringi       1.4.3      2019-03-12 [1] CRAN (R 3.6.0)                 
 stringr       1.4.0      2019-02-10 [1] CRAN (R 3.6.0)                 
 tibble        2.1.3      2019-06-06 [1] CRAN (R 3.6.0)                 
 tidyselect    0.2.5      2018-10-11 [1] CRAN (R 3.6.0)                 
 withr         2.1.2      2018-03-15 [1] CRAN (R 3.6.0)                 
 xfun          0.10       2019-10-01 [1] CRAN (R 3.6.0)                 
 xml2          1.2.2      2019-08-09 [1] CRAN (R 3.6.0)                 
 yaml          2.2.0      2018-07-25 [1] CRAN (R 3.6.0)                 

[1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library

^{Created on 2020-01-17 by the reprex package (v0.3.0)}

[yonicd]

in addition the csv generated output returns in sim$fnames_io

> from_csv@sim$fnames_oi
 [1] "alpha[1]"      "alpha[2]"      "alpha[3]"      "alpha[4]"      "alpha[5]"      "alpha[6]"      "alpha[7]"      "alpha[8]"      "alpha[9]"      "alpha[10]"     "alpha[11]"    
[12] "alpha[12]"     "alpha[13]"     "alpha[14]"     "alpha[15]"     "alpha[16]"     "alpha[17]"     "alpha[18]"     "alpha[19]"     "alpha[20]"     "alpha[21]"     "alpha[22]"    
[23] "alpha[23]"     "alpha[24]"     "alpha[25]"     "alpha[26]"     "alpha[27]"     "alpha[28]"     "alpha[29]"     "alpha[30]"     "beta[1]"       "beta[2]"       "beta[3]"      
[34] "beta[4]"       "beta[5]"       "beta[6]"       "beta[7]"       "beta[8]"       "beta[9]"       "beta[10]"      "beta[11]"      "beta[12]"      "beta[13]"      "beta[14]"     
[45] "beta[15]"      "beta[16]"      "beta[17]"      "beta[18]"      "beta[19]"      "beta[20]"      "beta[21]"      "beta[22]"      "beta[23]"      "beta[24]"      "beta[25]"     
[56] "beta[26]"      "beta[27]"      "beta[28]"      "beta[29]"      "beta[30]"      "mu_alpha"      "mu_beta"       "sigmasq_y"     "sigmasq_alpha" "sigmasq_beta"  "sigma_y"      
[67] "sigma_alpha"   "sigma_beta"    "alpha0"        "lp__"   

which is a problem given the different names in the names(sim$samples)

> names(from_csv@sim$samples[[1]])
 [1] "alpha.1"       "alpha.2"       "alpha.3"       "alpha.4"       "alpha.5"       "alpha.6"       "alpha.7"       "alpha.8"       "alpha.9"       "alpha.10"      "alpha.11"     
[12] "alpha.12"      "alpha.13"      "alpha.14"      "alpha.15"      "alpha.16"      "alpha.17"      "alpha.18"      "alpha.19"      "alpha.20"      "alpha.21"      "alpha.22"     
[23] "alpha.23"      "alpha.24"      "alpha.25"      "alpha.26"      "alpha.27"      "alpha.28"      "alpha.29"      "alpha.30"      "beta.1"        "beta.2"        "beta.3"       
[34] "beta.4"        "beta.5"        "beta.6"        "beta.7"        "beta.8"        "beta.9"        "beta.10"       "beta.11"       "beta.12"       "beta.13"       "beta.14"      
[45] "beta.15"       "beta.16"       "beta.17"       "beta.18"       "beta.19"       "beta.20"       "beta.21"       "beta.22"       "beta.23"       "beta.24"       "beta.25"      
[56] "beta.26"       "beta.27"       "beta.28"       "beta.29"       "beta.30"       "mu_alpha"      "mu_beta"       "sigmasq_y"     "sigmasq_alpha" "sigmasq_beta"  "sigma_y"      
[67] "sigma_alpha"   "sigma_beta"    "alpha0"        "lp__" 

[yonicd]

this also does not work

> rstan::stan_rhat(object = from_csv,'alpha')
Error in if (x@stan_args[[1L]]$method == "variational") stop("Plot only available for models estimated using MCMC",  : 
  argument is of length zero

> rstan::stan_rhat(object = fit,'alpha')
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

this is because the value for from_csv@stan_args[[1L]]$method is NULL

from_csv@stan_args[[1L]]$method
NULL
from_csv@stan_args[[1L]]$method <- 'sampling'
rstan::stan_rhat(object = from_csv,'alpha')
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.