Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Duplicated/redundant computations with horseshoe prior? #1167

Closed
mbjoseph opened this issue May 25, 2021 · 2 comments
Closed

Duplicated/redundant computations with horseshoe prior? #1167

mbjoseph opened this issue May 25, 2021 · 2 comments
Labels
bug
Milestone
brms 2.16.0

Comments

@mbjoseph
Copy link

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
Copy link
Owner

This is strange. Thank you. Will investigate.

@paul-buerkner paul-buerkner added this to the brms 2.15.0++ milestone May 26, 2021
paul-buerkner added a commit that referenced this issue May 26, 2021
@paul-buerkner
fix issue #1167
fb4b59b
@paul-buerkner
Copy link
Owner

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
bug
Projects
None yet
Milestone
brms 2.16.0
Development

No branches or pull requests
2 participants
@mbjoseph @paul-buerkner