wwiecek / baggr

R package for Bayesian meta-analysis models, using Stan
GNU General Public License v3.0
48 stars 12 forks source link

`example(baggr)` fails with new StanHeaders due to impossible initialization #111

Closed bgoodri closed 3 years ago

bgoodri commented 3 years ago

Currently, example(baggr) somehow works with the StanHeaders / rstan on CRAN but it really should not because baggr tries to use a uniform prior on the coefficients with both the lower bound and the upper bound equal to zero. Specifically, stan_data is

List of 16
 $ K                  : int 8
 $ theta_hat_k        : num [1:8] 1 -1 0.5 -0.5 0.7 -0.7 1.3 -1.3
 $ se_theta_k         : num [1:8] 1 1 1 1 1 1 1 1
 $ K_test             : num 0
 $ test_theta_hat_k   : num[0 (1d)] 
 $ test_se_theta_k    : num[0 (1d)] 
 $ Nc                 : num 0
 $ X                  : num[1:8, 0 ] 
 $ X_test             : num[0 , 0 ] 
 $ pooling_type       : num 1
 $ prior_hypermean_fam: num 1
 $ prior_hypermean_val: num [1:3] 0 13 0
 $ prior_hypersd_fam  : num 0
 $ prior_hypersd_val  : num [1:3] 0 9.9 0
 $ prior_beta_fam     : num 0
 $ prior_beta_val     : num [1:3] 0 0 0

the last line for prior_beta_val is a vector of three zeros and the line before that indicates that the uniform prior is being used. When all that gets passed to your

 real prior_increment_vec (int family, vector y, vector pars) {
    real inc;
    if(family == 0)
      inc = uniform_lpdf(y | pars[1], pars[2]); // pars[1] == 0 == pars[2]
    if(family == 1)
      inc = normal_lpdf(y | pars[1], pars[2]);
    if(family == 2)
      inc = cauchy_lpdf(y | pars[1], pars[2]);
    if(family == 3)
      inc = student_t_lpdf(y | pars[1], pars[2], pars[3]);
    return inc;
  }

you get an exception from the new version of StanHeaders that we are struggling to get onto CRAN:

SAMPLING FOR MODEL 'rubin' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: Exception: uniform_lpdf: Upper bound parameter is 0, but must be greater than 0  (in '/functions/prior_increment.stan' at line 16; included from 'model_rubin' at line 2)
  (in 'model_rubin' at line 80)

So, the R code should be fixed so that it specifies a non-degenerate prior. But even if you specified a lower bound that was lower than the upper bound, the Stan code is subject to exceptions and / or poor sampling because leapfrog steps can go outside those bounds. To overcome that, you need to change the declaration in the parameters block of your Stan programs to have bounds that are conditional on prior_beta_fam like:

  vector<lower = (prior_beta_fam == 0 ? prior_beta_val[1] : negative_infinity()),
         upper = (prior_beta_fam == 0 ? prior_beta_val[2] : positive_infinity())>[Nc] beta;

Even though this is not a bug with StanHeaders >= 2.26, it is a bug that prevents StanHeaders_2.26.x from passing the automatic reverse dependency checks. So please upload a new baggr to CRAN that fixes is this as soon as you possibly can.

wwiecek commented 3 years ago

Hi Ben! Thanks so much for this. The c(0,0,0) is faulty R code, I must have broken it with a recent update. Will look into conditional bounds, but will start with a quick fix to R code alone---will submit to CRAN next week with some other bug fixes. Cheers, Witold

wwiecek commented 3 years ago

@bgoodri the new version is on CRAN, I think this should work fine going forward; thanks again for letting me know! W