Posterior::as_draws_df() fails with stanreg object from stan_glm()

Interfaces rstanarm
1

Hello! I got an unexpected result when I attempted to use posterior::as_draws_df() on a stanreg object – an error, see reprex below, copied and pasted the offending line below.

posterior <- as_draws_df(f0)
#> Error: All variables in all chains must have the same length.

I think something is not working correctly, but hope I am wrong, and I missed something. :)

That should work, posterior::as_draws_df(f0), shouldn’t it?

# parameters to simulate linear model
N <- 30
a <- 5
b <- 5
sig <- 1

# data
x <- runif(N)
y <- a + b*x + rnorm(n = N, sd = sig)

# plot(x,y) # just checking

# model
library(rstanarm)
#> Loading required package: Rcpp
#> This is rstanarm version 2.21.4
#> - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
#> - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
#> - For execution on a local, multicore CPU with excess RAM we recommend calling
#>   options(mc.cores = parallel::detectCores())
options(mc.cores = parallel::detectCores())
f0 <- stan_glm(formula = y ~ x, data = data.frame(x, y), 
               refresh = 0)

# is it what I expected?
class(f0)
#> [1] "stanreg" "glm"     "lm"
f0
#> stan_glm
#>  family:       gaussian [identity]
#>  formula:      y ~ x
#>  observations: 30
#>  predictors:   2
#> ------
#>             Median MAD_SD
#> (Intercept) 4.5    0.3   
#> x           5.6    0.5   
#> 
#> Auxiliary parameter(s):
#>       Median MAD_SD
#> sigma 0.8    0.1   
#> 
#> ------
#> * For help interpreting the printed output see ?print.stanreg
#> * For info on the priors used see ?prior_summary.stanreg

# OK, take some draws
library(posterior)
#> This is posterior version 1.4.1
#> 
#> Attaching package: 'posterior'
#> The following objects are masked from 'package:stats':
#> 
#>     mad, sd, var
#> The following objects are masked from 'package:base':
#> 
#>     %in%, match
posterior <- as_draws_df(f0)
#> Error: All variables in all chains must have the same length.

Created on 2023-07-20 with reprex v2.0.2

Session info 
sessionInfo()
#> R version 4.2.3 (2023-03-15 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19044)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.utf8 
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] posterior_1.4.1 rstanarm_2.21.4 Rcpp_1.0.10    
#> 
#> loaded via a namespace (and not attached):
#>   [1] nlme_3.1-162         matrixStats_0.63.0   fs_1.6.2            
#>   [4] xts_0.13.1           threejs_0.3.3        rstan_2.26.22       
#>   [7] tensorA_0.36.2       R.cache_0.16.0       backports_1.4.1     
#>  [10] tools_4.2.3          utf8_1.2.3           R6_2.5.1            
#>  [13] DT_0.27              colorspace_2.1-0     withr_2.5.0         
#>  [16] tidyselect_1.2.0     gridExtra_2.3        prettyunits_1.1.1   
#>  [19] processx_3.8.1       curl_5.0.0           compiler_4.2.3      
#>  [22] cli_3.6.1            shinyjs_2.1.0        colourpicker_1.2.0  
#>  [25] checkmate_2.2.0      scales_1.2.1         dygraphs_1.1.1.6    
#>  [28] callr_3.7.3          stringr_1.5.0        digest_0.6.31       
#>  [31] StanHeaders_2.26.22  minqa_1.2.5          rmarkdown_2.21      
#>  [34] R.utils_2.12.2       base64enc_0.1-3      pkgconfig_2.0.3     
#>  [37] htmltools_0.5.5      lme4_1.1-33          styler_1.9.1        
#>  [40] fastmap_1.1.1        htmlwidgets_1.6.2    rlang_1.1.1         
#>  [43] rstudioapi_0.14      shiny_1.7.4          farver_2.1.1        
#>  [46] generics_0.1.3       zoo_1.8-12           jsonlite_1.8.4      
#>  [49] crosstalk_1.2.0      gtools_3.9.4         distributional_0.3.2
#>  [52] dplyr_1.1.2          R.oo_1.25.0          inline_0.3.19       
#>  [55] magrittr_2.0.3       loo_2.6.0            bayesplot_1.10.0    
#>  [58] Matrix_1.5-3         munsell_0.5.0        fansi_1.0.4         
#>  [61] abind_1.4-5          lifecycle_1.0.3      R.methodsS3_1.8.2   
#>  [64] stringi_1.7.12       yaml_2.3.7           MASS_7.3-58.2       
#>  [67] pkgbuild_1.4.0       plyr_1.8.8           grid_4.2.3          
#>  [70] parallel_4.2.3       promises_1.2.0.1     crayon_1.5.2        
#>  [73] miniUI_0.1.1.1       lattice_0.20-45      splines_4.2.3       
#>  [76] knitr_1.42           ps_1.7.5             pillar_1.9.0        
#>  [79] igraph_1.4.2         boot_1.3-28.1        markdown_1.6        
#>  [82] shinystan_2.6.0      reshape2_1.4.4       codetools_0.2-19    
#>  [85] stats4_4.2.3         rstantools_2.3.1     reprex_2.0.2        
#>  [88] glue_1.6.2           evaluate_0.20        V8_4.3.0            
#>  [91] RcppParallel_5.1.7   nloptr_2.0.3         vctrs_0.6.2         
#>  [94] httpuv_1.6.9         gtable_0.3.3         purrr_1.0.1         
#>  [97] ggplot2_3.4.2        xfun_0.39            mime_0.12           
#> [100] xtable_1.8-4         later_1.3.0          survival_3.5-3      
#> [103] tibble_3.2.1         shinythemes_1.2.0    ellipsis_0.3.2
1 Like
2

Work around, posting in case this helps someone else.

as.matrix(f0) gives the data in the format I needed (to send to bayesplot::mcmc_areas()). If this is a genuine issue, it is not preventing me from reaching my goals.

# parameters to simulate linear model
N <- 30
a <- 5
b <- 5
sig <- 1

# data
x <- runif(N)
y <- a + b*x + rnorm(n = N, sd = sig)

# plot(x,y) # just checking

# model
library(rstanarm)
#> Loading required package: Rcpp
#> This is rstanarm version 2.21.4
#> - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
#> - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
#> - For execution on a local, multicore CPU with excess RAM we recommend calling
#>   options(mc.cores = parallel::detectCores())
options(mc.cores = parallel::detectCores())
f0 <- stan_glm(formula = y ~ x, data = data.frame(x, y), 
               refresh = 0)

# is it what I expected?
class(f0)
#> [1] "stanreg" "glm"     "lm"
f0
#> stan_glm
#>  family:       gaussian [identity]
#>  formula:      y ~ x
#>  observations: 30
#>  predictors:   2
#> ------
#>             Median MAD_SD
#> (Intercept) 5.4    0.4   
#> x           3.9    0.6   
#> 
#> Auxiliary parameter(s):
#>       Median MAD_SD
#> sigma 0.9    0.1   
#> 
#> ------
#> * For help interpreting the printed output see ?print.stanreg
#> * For info on the priors used see ?prior_summary.stanreg

# OK, take some draws
library(posterior)
#> This is posterior version 1.4.1
#> 
#> Attaching package: 'posterior'
#> The following objects are masked from 'package:stats':
#> 
#>     mad, sd, var
#> The following objects are masked from 'package:base':
#> 
#>     %in%, match
posterior <- as_draws_df(f0)
#> Error: All variables in all chains must have the same length.

# Wait, what?

# this works though
my_matrix <- as.matrix(f0)

# my goal was to produce this plot, and I can do that with as.matrix()
bayesplot::mcmc_areas(x =  my_matrix)

Created on 2023-07-21 with reprex v2.0.2

Session info 
sessionInfo()
#> R version 4.2.3 (2023-03-15 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19044)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.utf8 
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] posterior_1.4.1 rstanarm_2.21.4 Rcpp_1.0.10    
#> 
#> loaded via a namespace (and not attached):
#>   [1] minqa_1.2.5          colorspace_2.1-0     ellipsis_0.3.2      
#>   [4] ggridges_0.5.4       markdown_1.6         base64enc_0.1-3     
#>   [7] fs_1.6.2             rstudioapi_0.14      farver_2.1.1        
#>  [10] rstan_2.26.22        DT_0.27              fansi_1.0.4         
#>  [13] xml2_1.3.4           codetools_0.2-19     splines_4.2.3       
#>  [16] R.methodsS3_1.8.2    knitr_1.42           shinythemes_1.2.0   
#>  [19] bayesplot_1.10.0     jsonlite_1.8.4       nloptr_2.0.3        
#>  [22] R.oo_1.25.0          shiny_1.7.4          compiler_4.2.3      
#>  [25] httr_1.4.5           backports_1.4.1      Matrix_1.5-3        
#>  [28] fastmap_1.1.1        cli_3.6.1            later_1.3.0         
#>  [31] htmltools_0.5.5      prettyunits_1.1.1    tools_4.2.3         
#>  [34] igraph_1.4.2         gtable_0.3.3         glue_1.6.2          
#>  [37] reshape2_1.4.4       dplyr_1.1.2          V8_4.3.0            
#>  [40] styler_1.9.1         vctrs_0.6.2          nlme_3.1-162        
#>  [43] crosstalk_1.2.0      tensorA_0.36.2       xfun_0.39           
#>  [46] stringr_1.5.0        ps_1.7.5             lme4_1.1-33         
#>  [49] mime_0.12            miniUI_0.1.1.1       lifecycle_1.0.3     
#>  [52] gtools_3.9.4         MASS_7.3-58.2        zoo_1.8-12          
#>  [55] scales_1.2.1         colourpicker_1.2.0   promises_1.2.0.1    
#>  [58] parallel_4.2.3       inline_0.3.19        shinystan_2.6.0     
#>  [61] yaml_2.3.7           curl_5.0.0           gridExtra_2.3       
#>  [64] ggplot2_3.4.2        loo_2.6.0            StanHeaders_2.26.22 
#>  [67] stringi_1.7.12       highr_0.10           dygraphs_1.1.1.6    
#>  [70] checkmate_2.2.0      boot_1.3-28.1        pkgbuild_1.4.0      
#>  [73] rlang_1.1.1          pkgconfig_2.0.3      matrixStats_0.63.0  
#>  [76] distributional_0.3.2 evaluate_0.20        lattice_0.20-45     
#>  [79] purrr_1.0.1          rstantools_2.3.1     htmlwidgets_1.6.2   
#>  [82] labeling_0.4.2       processx_3.8.1       tidyselect_1.2.0    
#>  [85] plyr_1.8.8           magrittr_2.0.3       R6_2.5.1            
#>  [88] generics_0.1.3       pillar_1.9.0         withr_2.5.0         
#>  [91] xts_0.13.1           survival_3.5-3       abind_1.4-5         
#>  [94] tibble_3.2.1         crayon_1.5.2         utf8_1.2.3          
#>  [97] rmarkdown_2.21       grid_4.2.3           callr_3.7.3         
#> [100] threejs_0.3.3        reprex_2.0.2         digest_0.6.31       
#> [103] xtable_1.8-4         R.cache_0.16.0       httpuv_1.6.9        
#> [106] R.utils_2.12.2       RcppParallel_5.1.7   stats4_4.2.3        
#> [109] munsell_0.5.0        shinyjs_2.1.0
2 Likes
3

A while ago I was planning on making this work with rstanarm, but apparently I never ended up doing it and neither did anyone else: Implement `posterior::as_draws()` for rstanarm objects · Issue #540 · stan-dev/rstanarm · GitHub

For now I think you can use as_draws_df(as.array(f0)) to produce what as_draws_df(f0) should produce. I’ll try to remember to get this into the next rstanarm release.

4
5

Many thanks @jonah !