Hierarchical Beta-Binomial Model with Order Constraints

Modeling
1

I’m working on a hierarchical beta-binomial model with order constraints for prevalence in different age groups. Here is the description of the model and the Stan code I’ve written so far.

Model Description

The hierarchical model can be described as follows:

The number of cases y_i in each age group i follows a binomial distribution:
y_i \sim \text{Binomial}(n_i, \pi_i)
where n_i is the number of individuals in age group i and \pi_i is the prevalence.

The prevalence \pi_i follows a beta distribution with order constraints:
\pi_i \sim \text{Beta}(\alpha_i, \beta_i) I(\pi_{i-1}, \pi_{i+1})
where I(\pi_{i-1}, \pi_{i+1}) is an indicator function ensuring that \pi_{i-1} \leq \pi_i \leq \pi_{i+1} .

Hyperpriors for (\alpha) and (\beta) are specified as truncated normal distributions:
\alpha_i \sim \text{Normal}(0, 1000) T[0, \infty)
\beta_i \sim \text{Normal}(0, 1000) T[0, \infty)

Stan Code

Below is the Stan code I have written. I would appreciate any feedback or suggestions for improvements.

stan

data {
  int<lower=1> Nage;           // Number of age groups
  vector[Nage] age;            // Age vector
  int<lower=0> posi[Nage];     // Number of positive cases
  int<lower=0> ni[Nage];       // Number of trials
}

parameters {
  vector<lower=0>[Nage] alpha;
  vector<lower=0>[Nage] beta;
   vector<lower=0,upper = 1>[Nage] pi;
}

model {
  for( i in 1:Nage){
  alpha[i] ~ normal(0,sqrt(1/0.001));
  beta[i] ~ normal(0,sqrt(1/0.00001));
  }
  
  pi[1] ~ beta(alpha[1],beta[1])T[,pi[2]];
  
  
  for(k in 2:(Nage-1)){
  
  pi[k] ~ beta(alpha[k],beta[k])T[pi[k-1],pi[k+1]];

  }
  pi[Nage] ~ beta(alpha[Nage],beta[Nage])T[pi[Nage-1],0.96]; 
  
  for( i in 1:Nage){
  posi[i] ~ binomial(ni[i], pi[i]);
  }
  
}

When I run the model, I encounter the following error message:
Chain 1: Initialization between (-2, 2) failed after 100 attempts. Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model. [1] “Error : Initialization failed.” [1] “error occurred during calling the sampler; sampling not done”

Do you have any ideas to fix this problem? Thank you for your help!

2

In your current framework, stan initializes pi randomly, which most likely results in the ordered constraints being violated and hence returning log(0).

If you want to constrain these to be positive ordered you need to use the positive ordered datatype in stan. You however also have the constraints that 0<pi<1 which introduces additional complexity which isn’t trivial to deal with. You may want to consider making the transformation of pi to some unconstrained space

\pi_i=1/(1+exp(-x_i))

You can then define the an ordered datatype on x and use the parameter transformation to specify the density on x. That is on the log scale log(P(x)) = log(|d\pi/dx|) + log(P(\pi)).

Simplified code would look something like

data {
  int<lower=1> Nage;           // Number of age groups
  vector[Nage] age;            // Age vector
  array[Nage] int<lower=0> posi;     // Number of positive cases
  array[Nage] int<lower=0> ni;       // Number of trials
}

parameters {
  vector<lower=0>[Nage] alpha;
  vector<lower=0>[Nage] beta;
  ordered[Nage] pi_unconstrained;
}

transformed parameters {
  vector[Nage] pi = 1/(1+exp(-pi_unconstrained)); // positive ordered since pi_unconstrained is ordered, and constrained between [0,1]
}

model {
  alpha ~ normal(0,sqrt(1/0.001));
  beta ~ normal(0,sqrt(1/0.00001));
  
  pi ~ beta(alpha,beta);
  target += sum(-pi_unconstrained)-2*sum(log(1+exp(pi_unconstrained))); // Jacobian adjustment for parameter transformation
  
  posi ~ binomial(ni, pi);
}

You additionally constrain pi here which could also result in log(0) being returned. To account for this you would have to change the transformed parameter block to

...
transformed parameters {
  vector[Nage] pi = 0.96/(1+exp(-pi_unconstrained)); // positive ordered since pi_unconstrained is ordered, and constrained between [0,0.96]
}
...
1 Like