Zero-Inflated Negative-Binomial model

Modeling
1

Had a negative-binomial model that sampled very well. But the fit showed inflated zeros. So I’ve tried to parameterize both p (binomial parameter) and lambda (nbinom) by a linear combination. However, the sampling simply freezes after repeatedly warning me that:
“Exception: neg_binomial_2_lpmf: Location parameter[1] is 0/Inf, but must be positive finite!”
(I don’t know if this warning is actually important for my freezing problem)

Code below

data{
    int rows;
    int dm[rows];
    int htn[rows];
    int hpl[rows];
    int smk[rows];
    int metSyn_alpha[6];
    int metSyn[rows];
    vector[rows] metSynV;
    int etoh[rows];
    int etohXsmk[rows];
    int cad[rows];
    vector[rows] bmi;
    int stroke_tia[rows];
    int prev[rows];
    int mRS_preop_alpha[6];
    int mRS_preop[rows];
    vector[rows] mRS_preopV;
    int gcs_alpha[12];
    int gcs[rows];
    vector[rows] gcsV;
    vector[rows] age;
    vector[rows] time;
    int sex[rows];
    int rupt[rows];
    int migraine[rows];
    int mRS_discharge_alpha[6];
    int mRS_discharge[rows];
    vector[rows] mRS_dischargeV;
    int insXrupt[rows];
    vector[rows] dist;
    int comp[rows];
    vector[rows] mean_lnInc;
    int FU[rows];
    vector[rows] sd_lnInc;
    vector[rows] sem_lnInc;
    int rowid[rows];
    int ins[rows];
    real grandSEM_lnInc;
    real grandMean_lnInc;
}
parameters{
    vector[4] a_ins;
    vector[rows] nu;
    vector[rows] stdN1;
    vector[rows] stdN2;
    vector[4] b_lnInc;
    real b_time;
    real<lower=0> phi;
    //real<lower=1, upper=1> p; // NEW
}
transformed parameters{
  vector[rows] lnInc;
  vector[rows] lnIncStd;
  
  for(i in 1:rows) lnInc[i] = ((nu[rowid[i]] + 
                                stdN2[rowid[i]] * sem_lnInc[i]) + 
                               stdN1[rowid[i]] * sd_lnInc[i]);
  for(i in 1:rows) lnIncStd[i] = (lnInc[i] - mean_lnInc[i]) / sd_lnInc[i];
}
model{
    vector[rows] p; // NEW
    vector[rows] lambda;
    phi ~ exponential( 2 );

    b_lnInc ~ normal( 0 , 0.25 );
    stdN2 ~ normal( 0 , 1 );
    stdN1 ~ normal( 0 , 1 );
    nu ~ normal( grandMean_lnInc , grandSEM_lnInc );
    a_ins ~ normal( 0 , 0.25);
    b_time ~ normal( 1 , 0.25 );
    
    for ( i in 1:rows ) { // NEW
      p[i] = a_ins[insXrupt[i]] + 
             b_lnInc[insXrupt[i]] * lnIncStd[i] + 
             b_time * time[i];
      p[i] = inv_logit(p[i]);
    }
    
    for ( i in 1:rows ) {
        lambda[i] = a_ins[insXrupt[i]] + 
                    b_lnInc[insXrupt[i]] * lnIncStd[i] + 
                    b_time * time[i]; 
        lambda[i] = exp(lambda[i]);
    }
    
    for(i in 1:rows){ // NEW
      if (FU[i] == 0)
        target += log_sum_exp(bernoulli_lpmf(1 | p),
                              bernoulli_lpmf(0 | p) 
                              + neg_binomial_2_lpmf(FU[i] | lambda, phi));
      else
        target += bernoulli_lpmf(0 | p)
                  + neg_binomial_2_lpmf(FU[i] | lambda, phi);
    
    //FU ~ neg_binomial_2( lambda , phi ); REMOVED
  }
}

Many thanks!

2

I think you need to add the i-index for lambda in the last for-loop.

2 Likes
3

Yes, that worked.