search

paul-buerkner / brms

brms R package for Bayesian generalized multivariate non-linear multilevel models using Stan
https://paul-buerkner.github.io/brms/
GNU General Public License v2.0
1.28k stars 184 forks source link

Duplicated/redundant computations with horseshoe prior? #1167

Closed mbjoseph closed 3 years ago

mbjoseph commented 3 years ago

I was playing with using a horseshoe prior in brms and looking at the generated stan code -- it looks like there's a redundant/duplicated function call in the transformed parameters block to the horseshoe function. Is this intentional? 

library(brms)
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.15.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
#> 
#> Attaching package: 'brms'
#> The following object is masked from 'package:stats':
#> 
#>     ar

n <- 100
y <- rpois(n, lambda = 1)
x <- rnorm(n)

d <- data.frame(x = x, y = y)

fit <- brm(y ~ x, data = d,
           family = poisson,
           prior = set_prior("horseshoe(1)"),
           backend = "cmdstanr", 
           silent = 2, refresh = 0)
#> Running MCMC with 4 sequential chains...
#> 
#> Chain 1 finished in 0.1 seconds.
#> Chain 2 finished in 0.2 seconds.
#> Chain 3 finished in 0.2 seconds.
#> Chain 4 finished in 0.2 seconds.
#> 
#> All 4 chains finished successfully.
#> Mean chain execution time: 0.2 seconds.
#> Total execution time: 1.1 seconds.
#> 
#> Warning: 718 of 4000 (18.0%) transitions ended with a divergence.
#> This may indicate insufficient exploration of the posterior distribution.
#> Possible remedies include: 
#>   * Increasing adapt_delta closer to 1 (default is 0.8) 
#>   * Reparameterizing the model (e.g. using a non-centered parameterization)
#>   * Using informative or weakly informative prior distributions
brms::stancode(fit)
#> // generated with brms 2.15.0
#> functions {
#>   /* Efficient computation of the horseshoe prior
#>    * see Appendix C.1 in https://projecteuclid.org/euclid.ejs/1513306866
#>    * Args:
#>    *   z: standardized population-level coefficients
#>    *   lambda: local shrinkage parameters
#>    *   tau: global shrinkage parameter
#>    *   c2: slap regularization parameter
#>    * Returns:
#>    *   population-level coefficients following the horseshoe prior
#>    */
#>   vector horseshoe(vector z, vector lambda, real tau, real c2) {
#>     int K = rows(z);
#>     vector[K] lambda2 = square(lambda);
#>     vector[K] lambda_tilde = sqrt(c2 * lambda2 ./ (c2 + tau^2 * lambda2));
#>     return z .* lambda_tilde * tau;
#>   }
#> }
#> data {
#>   int<lower=1> N;  // total number of observations
#>   int Y[N];  // response variable
#>   int<lower=1> K;  // number of population-level effects
#>   matrix[N, K] X;  // population-level design matrix
#>   // data for the horseshoe prior
#>   real<lower=0> hs_df;  // local degrees of freedom
#>   real<lower=0> hs_df_global;  // global degrees of freedom
#>   real<lower=0> hs_df_slab;  // slab degrees of freedom
#>   real<lower=0> hs_scale_global;  // global prior scale
#>   real<lower=0> hs_scale_slab;  // slab prior scale
#>   int prior_only;  // should the likelihood be ignored?
#> }
#> transformed data {
#>   int Kc = K - 1;
#>   matrix[N, Kc] Xc;  // centered version of X without an intercept
#>   vector[Kc] means_X;  // column means of X before centering
#>   for (i in 2:K) {
#>     means_X[i - 1] = mean(X[, i]);
#>     Xc[, i - 1] = X[, i] - means_X[i - 1];
#>   }
#> }
#> parameters {
#>   // local parameters for horseshoe prior
#>   vector[Kc] zb;
#>   vector<lower=0>[Kc] hs_local;
#>   real Intercept;  // temporary intercept for centered predictors
#>   // horseshoe shrinkage parameters
#>   real<lower=0> hs_global;  // global shrinkage parameters
#>   real<lower=0> hs_slab;  // slab regularization parameter
#> }
#> transformed parameters {
#>   vector[Kc] b;  // population-level effects
#>   // compute actual regression coefficients
#>   b = horseshoe(zb, hs_local, hs_global, hs_scale_slab^2 * hs_slab);
#>   // compute actual regression coefficients
#>   b = horseshoe(zb, hs_local, hs_global, hs_scale_slab^2 * hs_slab);
#> }
#> model {
#>   // likelihood including constants
#>   if (!prior_only) {
#>     target += poisson_log_glm_lpmf(Y | Xc, Intercept, b);
#>   }
#>   // priors including constants
#>   target += std_normal_lpdf(zb);
#>   target += student_t_lpdf(hs_local | hs_df, 0, 1)
#>     - rows(hs_local) * log(0.5);
#>   target += student_t_lpdf(Intercept | 3, 0, 2.5);
#>   target += student_t_lpdf(hs_global | hs_df_global, 0, hs_scale_global)
#>     - 1 * log(0.5);
#>   target += inv_gamma_lpdf(hs_slab | 0.5 * hs_df_slab, 0.5 * hs_df_slab);
#> }
#> generated quantities {
#>   // actual population-level intercept
#>   real b_Intercept = Intercept - dot_product(means_X, b);
#> }

Created on 2021-05-25 by the reprex package (v2.0.0)

Notice the duplicated line:

#>   // compute actual regression coefficients
#>   b = horseshoe(zb, hs_local, hs_global, hs_scale_slab^2 * hs_slab);

Here's my session info in case its useful: 

sessioninfo::session_info()
─ Session info ─────────────────────────────────────────────────────────────────────────────────────────────
 setting  value                       
 version  R version 4.1.0 (2021-05-18)
 os       Ubuntu 20.04.2 LTS          
 system   x86_64, linux-gnu           
 ui       RStudio                     
 language (EN)                        
 collate  en_US.UTF-8                 
 ctype    en_US.UTF-8                 
 tz       America/Denver              
 date     2021-05-25                  

─ Packages ─────────────────────────────────────────────────────────────────────────────────────────────────
 package        * version    date       lib source                             
 abind            1.4-5      2016-07-21 [1] CRAN (R 4.1.0)                     
 assertthat       0.2.1      2019-03-21 [1] CRAN (R 4.1.0)                     
 backports        1.2.1      2020-12-09 [1] CRAN (R 4.1.0)                     
 base64enc        0.1-3      2015-07-28 [1] CRAN (R 4.1.0)                     
 bayesplot        1.8.0      2021-01-10 [1] CRAN (R 4.1.0)                     
 boot             1.3-28     2021-05-03 [4] CRAN (R 4.0.5)                     
 bridgesampling   1.1-2      2021-04-16 [1] CRAN (R 4.1.0)                     
 brms           * 2.15.0     2021-03-14 [1] CRAN (R 4.1.0)                     
 Brobdingnag      1.2-6      2018-08-13 [1] CRAN (R 4.1.0)                     
 callr            3.7.0      2021-04-20 [1] CRAN (R 4.1.0)                     
 checkmate        2.0.0      2020-02-06 [1] CRAN (R 4.1.0)                     
 cli              2.5.0      2021-04-26 [1] CRAN (R 4.1.0)                     
 clipr            0.7.1      2020-10-08 [1] CRAN (R 4.1.0)                     
 cmdstanr         0.4.0.9000 2021-05-20 [1] Github (stan-dev/cmdstanr@8cbfbc6) 
 coda             0.19-4     2020-09-30 [1] CRAN (R 4.1.0)                     
 codetools        0.2-18     2020-11-04 [4] CRAN (R 4.0.3)                     
 colorspace       2.0-1      2021-05-04 [1] CRAN (R 4.1.0)                     
 colourpicker     1.1.0      2020-09-14 [1] CRAN (R 4.1.0)                     
 crayon           1.4.1      2021-02-08 [1] CRAN (R 4.1.0)                     
 crosstalk        1.1.1      2021-01-12 [1] CRAN (R 4.1.0)                     
 curl             4.3.1      2021-04-30 [1] CRAN (R 4.1.0)                     
 data.table       1.14.0     2021-02-21 [1] CRAN (R 4.1.0)                     
 DBI              1.1.1      2021-01-15 [1] CRAN (R 4.1.0)                     
 digest           0.6.27     2020-10-24 [1] CRAN (R 4.1.0)                     
 distributional   0.2.2      2021-02-02 [1] CRAN (R 4.1.0)                     
 dplyr            1.0.6      2021-05-05 [1] CRAN (R 4.1.0)                     
 DT               0.18       2021-04-14 [1] CRAN (R 4.1.0)                     
 dygraphs         1.1.1.6    2018-07-11 [1] CRAN (R 4.1.0)                     
 ellipsis         0.3.2      2021-04-29 [1] CRAN (R 4.1.0)                     
 evaluate         0.14       2019-05-28 [1] CRAN (R 4.1.0)                     
 fansi            0.4.2      2021-01-15 [1] CRAN (R 4.1.0)                     
 farver           2.1.0      2021-02-28 [1] CRAN (R 4.1.0)                     
 fastmap          1.1.0      2021-01-25 [1] CRAN (R 4.1.0)                     
 fs               1.5.0      2020-07-31 [1] CRAN (R 4.1.0)                     
 gamm4            0.2-6      2020-04-03 [1] CRAN (R 4.1.0)                     
 generics         0.1.0      2020-10-31 [1] CRAN (R 4.1.0)                     
 ggplot2          3.3.3      2020-12-30 [1] CRAN (R 4.1.0)                     
 ggridges         0.5.3      2021-01-08 [1] CRAN (R 4.1.0)                     
 glue             1.4.2      2020-08-27 [1] CRAN (R 4.1.0)                     
 gridExtra        2.3        2017-09-09 [1] CRAN (R 4.1.0)                     
 gtable           0.3.0      2019-03-25 [1] CRAN (R 4.1.0)                     
 gtools           3.8.2      2020-03-31 [1] CRAN (R 4.1.0)                     
 highr            0.9        2021-04-16 [1] CRAN (R 4.1.0)                     
 htmltools        0.5.1.1    2021-01-22 [1] CRAN (R 4.1.0)                     
 htmlwidgets      1.5.3      2020-12-10 [1] CRAN (R 4.1.0)                     
 httpuv           1.6.1      2021-05-07 [1] CRAN (R 4.1.0)                     
 igraph           1.2.6      2020-10-06 [1] CRAN (R 4.1.0)                     
 inline           0.3.18     2021-05-18 [1] CRAN (R 4.1.0)                     
 jsonlite         1.7.2      2020-12-09 [1] CRAN (R 4.1.0)                     
 knitr            1.33       2021-04-24 [1] CRAN (R 4.1.0)                     
 later            1.2.0      2021-04-23 [1] CRAN (R 4.1.0)                     
 lattice          0.20-44    2021-05-02 [4] CRAN (R 4.0.5)                     
 lifecycle        1.0.0      2021-02-15 [1] CRAN (R 4.1.0)                     
 lme4             1.1-27     2021-05-15 [1] CRAN (R 4.1.0)                     
 loo              2.4.1      2020-12-09 [1] CRAN (R 4.1.0)                     
 magrittr         2.0.1      2020-11-17 [1] CRAN (R 4.1.0)                     
 markdown         1.1        2019-08-07 [1] CRAN (R 4.1.0)                     
 MASS             7.3-54     2021-05-03 [4] CRAN (R 4.0.5)                     
 Matrix           1.3-3      2021-05-04 [4] CRAN (R 4.0.5)                     
 matrixStats      0.58.0     2021-01-29 [1] CRAN (R 4.1.0)                     
 mgcv             1.8-35     2021-04-18 [4] CRAN (R 4.0.5)                     
 mime             0.10       2021-02-13 [1] CRAN (R 4.1.0)                     
 miniUI           0.1.1.1    2018-05-18 [1] CRAN (R 4.1.0)                     
 minqa            1.2.4      2014-10-09 [1] CRAN (R 4.1.0)                     
 munsell          0.5.0      2018-06-12 [1] CRAN (R 4.1.0)                     
 mvtnorm          1.1-1      2020-06-09 [1] CRAN (R 4.1.0)                     
 nlme             3.1-152    2021-02-04 [4] CRAN (R 4.0.3)                     
 nloptr           1.2.2.2    2020-07-02 [1] CRAN (R 4.1.0)                     
 pillar           1.6.1      2021-05-16 [1] CRAN (R 4.1.0)                     
 pkgbuild         1.2.0      2020-12-15 [1] CRAN (R 4.1.0)                     
 pkgconfig        2.0.3      2019-09-22 [1] CRAN (R 4.1.0)                     
 plyr             1.8.6      2020-03-03 [1] CRAN (R 4.1.0)                     
 posterior        0.1.5      2021-05-20 [1] Github (stan-dev/posterior@7d17795)
 prettyunits      1.1.1      2020-01-24 [1] CRAN (R 4.1.0)                     
 processx         3.5.2      2021-04-30 [1] CRAN (R 4.1.0)                     
 projpred         2.0.2      2020-10-28 [1] CRAN (R 4.1.0)                     
 promises         1.2.0.1    2021-02-11 [1] CRAN (R 4.1.0)                     
 ps               1.6.0      2021-02-28 [1] CRAN (R 4.1.0)                     
 purrr            0.3.4      2020-04-17 [1] CRAN (R 4.1.0)                     
 R6               2.5.0      2020-10-28 [1] CRAN (R 4.1.0)                     
 Rcpp           * 1.0.6      2021-01-15 [1] CRAN (R 4.1.0)                     
 RcppParallel     5.1.4      2021-05-04 [1] CRAN (R 4.1.0)                     
 reprex           2.0.0      2021-04-02 [1] CRAN (R 4.1.0)                     
 reshape2         1.4.4      2020-04-09 [1] CRAN (R 4.1.0)                     
 rlang            0.4.11     2021-04-30 [1] CRAN (R 4.1.0)                     
 rmarkdown        2.8        2021-05-07 [1] CRAN (R 4.1.0)                     
 rsconnect        0.8.18     2021-05-24 [1] CRAN (R 4.1.0)                     
 rstan            2.21.2     2020-07-27 [1] CRAN (R 4.1.0)                     
 rstantools       2.1.1      2020-07-06 [1] CRAN (R 4.1.0)                     
 rstudioapi       0.13       2020-11-12 [1] CRAN (R 4.1.0)                     
 scales           1.1.1      2020-05-11 [1] CRAN (R 4.1.0)                     
 sessioninfo      1.1.1      2018-11-05 [1] CRAN (R 4.1.0)                     
 shiny            1.6.0      2021-01-25 [1] CRAN (R 4.1.0)                     
 shinyjs          2.0.0      2020-09-09 [1] CRAN (R 4.1.0)                     
 shinystan        2.5.0      2018-05-01 [1] CRAN (R 4.1.0)                     
 shinythemes      1.2.0      2021-01-25 [1] CRAN (R 4.1.0)                     
 StanHeaders      2.21.0-7   2020-12-17 [1] CRAN (R 4.1.0)                     
 stringi          1.6.2      2021-05-17 [1] CRAN (R 4.1.0)                     
 stringr          1.4.0      2019-02-10 [1] CRAN (R 4.1.0)                     
 styler           1.4.1      2021-03-30 [1] CRAN (R 4.1.0)                     
 tensorA          0.36.2     2020-11-19 [1] CRAN (R 4.1.0)                     
 threejs          0.3.3      2020-01-21 [1] CRAN (R 4.1.0)                     
 tibble           3.1.2      2021-05-16 [1] CRAN (R 4.1.0)                     
 tidyselect       1.1.1      2021-04-30 [1] CRAN (R 4.1.0)                     
 utf8             1.2.1      2021-03-12 [1] CRAN (R 4.1.0)                     
 V8               3.4.2      2021-05-01 [1] CRAN (R 4.1.0)                     
 vctrs            0.3.8      2021-04-29 [1] CRAN (R 4.1.0)                     
 withr            2.4.2      2021-04-18 [1] CRAN (R 4.1.0)                     
 xfun             0.23       2021-05-15 [1] CRAN (R 4.1.0)                     
 xtable           1.8-4      2019-04-21 [1] CRAN (R 4.1.0)                     
 xts              0.12.1     2020-09-09 [1] CRAN (R 4.1.0)                     
 yaml             2.2.1      2020-02-01 [1] CRAN (R 4.1.0)                     
 zoo              1.8-9      2021-03-09 [1] CRAN (R 4.1.0)                     

[1] /home/max/R/x86_64-pc-linux-gnu-library/4.1
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library
paul-buerkner commented 3 years ago

This is strange. Thank you. Will investigate.

paul-buerkner commented 3 years ago

It was a problem with partial matching of list elements when using the $ operator... Should now be fixed.