Hierarchical Bayesian Poisson regression model

Modeling
1

Hello all,

I am trying to fit a hierarchical Bayesian Poisson regression model with Stan. First I tried to fit a simpler model: a Bayesian Poisson regression model given below.

\begin{array}{l} {y_{ij}} \sim {\rm{Poisson}}\left( {{\lambda _{ij}}} \right)\\ \log \left( {{\lambda _{ij}}} \right) = {\alpha _j} + {\beta _j}{t_{ij}}\\ {\alpha _j} \sim N\left( {0,100} \right)\\ {\beta _j} \sim N\left( {0,100} \right) \end{array}

However, (1) the effective sample size was too low and (2) there were many errors in evaluating the log probability at the initial value, especially in Chain 4. Do you have any idea about how to fix these issues? Thanks.

Julian

Stan code

data {
  int<lower=1> I;            // number of points in each group
  int<lower=1> J;            // number of groups
  int<lower=1> N;           // total number of points of all groups
  int<lower=1> Y[N];         // Count outcome    
  real<lower=0> time[N];     //
  int<lower=1> index_group[N]; //
}

parameters {
    real beta[J];             //
}

transformed parameters{
  real<lower=0> log_lambda[N];  
  for (i in 1:N){
    log_lambda[i]=alpha[index_group[i]]+beta[index_group[i]]*time[i];
  }
}

model {
  target +=normal_lpdf(alpha|0,100);
  target +=normal_lpdf(beta|0,100);
  for (i in 1:N){
    target += poisson_log_lpmf(Y[i]|log_lambda[i]);
  }
}

R code and data

HB_data_1d=list(Y=Y_vec,I=12,J=5,N=60,time=time_vec,index_group=index_group_vec)

n_sam=10000
n_warmup=round(n_sam/5,0)
n_chain=4

HBR_1d_fit=stan(file='1d_HBR_stan.stan',data = HB_data_1d,iter =(n_sam+n_warmup), warmup=n_warmup,
                chains = n_chain, control = list(adapt_delta = 0.95, max_treedepth = 20))

Y_vec

[1]  2444  6956 13912 18048 21996 27072 31960 41924 44932 47000 65048 91932  1170  2990  4940  7020 11960 17940 21060 23660 33020 40040 53040
[24] 73190  1016  2223  3175  5652 12002  9589 14415 17018 17653 20003 31560 34481   832  2176  2432  6976  6272 10112 20992 19136 25792 28032
[47] 32960 44032    68  1085  3189  4478  6513  8412 10108 11533 12822 16961 23813 29036

time_vec

[1]  11.25  35.25  41.25  57.25  65.25  81.25  89.25 105.25 113.25 129.25 137.25 153.25  12.00  20.00  28.00  36.00  49.00  61.50  73.00
[20]  84.50  89.50 107.00 119.50 131.50  12.00  14.25  16.50  22.00  27.50  31.75  36.00  41.00  46.00  51.50  63.50  75.50   2.75   6.25
[39]   9.75  13.75  17.75  23.00  28.00  33.75  41.75  49.75  57.75  65.75  17.00  25.00  32.00  40.00  44.75  49.50  55.00  60.00  65.00
[58]  69.00  73.00  78.00

index_group_vec

[1] 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5

Output from Stan

SAMPLING FOR MODEL '1d_HBR_stan' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -10.8695, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -1.76582, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -9.66006, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -20.8099, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -15.5109, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -21.6343, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -10.7049, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -14.3437, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -11.1223, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -11.386, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -16.8763, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -0.0814848, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -22.6518, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -15.44, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -12.1711, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -3.08141, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -20.4154, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -0.506579, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -0.0449641, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -5.52714, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -5.18335, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -13.3998, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -10.623, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -20.5497, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -10.2542, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -3.43916, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -10.6235, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -16.2179, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: log_lambda[i_0__] is -0.903154, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 1: 
Chain 1: Gradient evaluation took 0 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1:  Elapsed Time: 0.236 seconds (Warm-up)
Chain 1:                1.156 seconds (Sampling)
Chain 1:                1.392 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL '1d_HBR_stan' NOW (CHAIN 2).
Chain 2: Rejecting initial value:
Chain 2:   Error evaluating the log probability at the initial value.
Chain 2: Exception: validate transformed params: log_lambda[i_0__] is -2.29231, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 2: Rejecting initial value:
Chain 2:   Error evaluating the log probability at the initial value.
Chain 2: Exception: validate transformed params: log_lambda[i_0__] is -12.1401, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 2: Rejecting initial value:
Chain 2:   Error evaluating the log probability at the initial value.
Chain 2: Exception: validate transformed params: log_lambda[i_0__] is -4.40557, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 2: Rejecting initial value:
Chain 2:   Error evaluating the log probability at the initial value.
Chain 2: Exception: validate transformed params: log_lambda[i_0__] is -0.0288678, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 2: Rejecting initial value:
Chain 2:   Error evaluating the log probability at the initial value.
Chain 2: Exception: validate transformed params: log_lambda[i_0__] is -30.3227, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 2: Rejecting initial value:
Chain 2:   Error evaluating the log probability at the initial value.
Chain 2: Exception: validate transformed params: log_lambda[i_0__] is -9.97668, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 2: Rejecting initial value:
Chain 2:   Error evaluating the log probability at the initial value.
Chain 2: Exception: validate transformed params: log_lambda[i_0__] is -12.4269, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 2: Rejecting initial value:
Chain 2:   Error evaluating the log probability at the initial value.
Chain 2: Exception: validate transformed params: log_lambda[i_0__] is -9.82757, but must be greater than or equal to 0  (in 'model132dc39234ddd_1d_HBR_stan' at line 22)

Chain 2: 
Chain 2: Gradient evaluation took 0 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: 
Chain 2:  Elapsed Time: 0.255 seconds (Warm-up)
Chain 2:                0.724 seconds (Sampling)
Chain 2:                0.979 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL '1d_HBR_stan' NOW (CHAIN 3).

SAMPLING FOR MODEL '1d_HBR_stan' NOW (CHAIN 4).

**Warning messages:**
**1: There were 2433 divergent transitions after warmup. Increasing adapt_delta above 0.95 may help. See**
**http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup **
**2: There were 2 chains where the estimated Bayesian Fraction of Missing Information was low. See**
**http://mc-stan.org/misc/warnings.html#bfmi-low **
**3: Examine the pairs() plot to diagnose sampling problems**
** **
**4: The largest R-hat is 5.17, indicating chains have not mixed.**
**Running the chains for more iterations may help. See**
**http://mc-stan.org/misc/warnings.html#r-hat **
**5: 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 **
**6: 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**
2

I didn’t read through your stan code to see if there were errors, but N(0, 100) priors on the log scale place a lot of mass in unrealistic portions of the parameter space (I don’t know what your data is, but are responses on the scale of e^{100} realistic?), and can cause the parameters to be weakly identified (leading to computational issues in stan).

3

From a cursory look at the code, log_lambda is defined to be positive, but beta can be negative (I don’t see the definition for alpha, I presume it can be negative too). If you define them with <lower=0>, you will sample from the half-Normal, and the constraint on log_lambda should be satisfied then.

4

Thanks! I will try weakly-informative priors instead.

5

Thanks. I lower bounded alpha and beta, and the “error evaluating the log probability at the initial value” rarely occurs now. However, the effective sample size is still too low.

The warning messages are as follows:

*Warning messages:*
*1: There were 15462 divergent transitions after warmup. Increasing adapt_delta above 0.95 may help. See*
*http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup *
*2: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See*
*http://mc-stan.org/misc/warnings.html#bfmi-low *
*3: Examine the pairs() plot to diagnose sampling problems*
* *
*4: The largest R-hat is 5.19, indicating chains have not mixed.*
*Running the chains for more iterations may help. See*
*http://mc-stan.org/misc/warnings.html#r-hat *
*5: 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 *
*6: 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*

Stan code

data {
  int<lower=1> I;            // number of points in each group
  int<lower=1> J;            // number of groups
  int<lower=1> N;           // total number of points of all groups
  int<lower=1> Y[N];         // Count outcome    
  real<lower=0> time[N];     //
  int<lower=1> index_group[N]; //
}

parameters {
  real<lower=0>alpha[J];             //
  real<lower=0>beta[J];             //
}

transformed parameters{
real<lower=0> log_lambda[N];  
  for (i in 1:N){
    log_lambda[i]=alpha[index_group[i]]+beta[index_group[i]]*time[i];
  }
}

model {
  target +=cauchy_lpdf(alpha|0,10);
  target +=cauchy_lpdf(beta|0,10);
  for (i in 1:N){
    target += poisson_log_lpmf(Y[i]|log_lambda[i]);
  }
}
6

Ok, I might have given you the wrong advice before. What’s wrong is the <lower=0> on log_lambda, so you can remove the lower bounds on alpha and beta too.

3 Likes
7

Hello researchers,

Thank you for your previous advice. The Stan model works in simpler case, so I extended the previous model to the hierarchical version.

However, the effective sample size is too low, indicating unreliable posteriors. The warning messages are:

Warning messages:
1: There were 39742 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded 
2: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low 
3: Examine the pairs() plot to diagnose sampling problems
 
4: The largest R-hat is 3.78, indicating chains have not mixed.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#r-hat 
5: 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 
6: 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

The Stan code is as follows:

data {
  int<lower=1> I;            // number of points in each group
  int<lower=1> J;            // number of groups
  int<lower=1> N;
  int<lower=1> Y[N];         // Count outcome    
  int<lower=1> K;         // number of groups of coefficients
  real<lower=0> time[N];              //
  int<lower=1> intensity[N];
  int<lower=1> index_group[N]; //
}

parameters {
  real<lower=0> alpha[J];            //
    real<lower=0> beta[J];           //
    real gamma[J];            // group specific intercept
    real<lower=0> sigma2[K];
    real mu_alpha;
    real mu_beta;
    real mu_gamma;
}

transformed parameters{
  real log_lambda[N]; 
  for (i in 1:N){
    log_lambda[i]=alpha[index_group[i]]+beta[index_group[i]]*time[i]+gamma[index_group[i]]*intensity[i];
  }
}

model {
  // hyperpriors
  target +=normal_lpdf(mu_alpha|0,5);         
  target +=normal_lpdf(mu_beta|0,2);         
  target +=normal_lpdf(mu_gamma|0,1);          
  target +=inv_gamma_lpdf(sigma2|2,2);
  
  // coefficients
  target +=normal_lpdf(alpha|mu_alpha,sigma2[1]);
  target +=normal_lpdf(beta|mu_beta,sigma2[2]);
  target +=normal_lpdf(gamma|mu_gamma,sigma2[3]);
  
  // likelihood
  for (i in 1:N){
    target += poisson_log_lpmf(Y[i]|log_lambda[i]);
  }
}

The traceplots below reveal that only chains of beta have a good mix.

traceplots
traceplots.png814Ă—417 13.1 KB

Do you have any suggestions to improve the effective sample size of alpha and gamma or in general the Stan model?

Thanks a lot.

Julian

8

You still have lower bounds on alpha and beta, try without them.

1 Like