Zero inflated Poisson vs Zero inflated negative exponential & over-dispersion

[James Camac]

Hi all,

I've fitted two multilevel models in stan a ZIP and a ZINB. What I want to do is examine whether there is over dispersion in the ZIP model.

One approach I've taken is to plot the predictions pearsons residuals (y - pred)/ sqrt(pred) against the predictions.

I look at the standard deviation of these residuals, which if I'm understanding correctly, should be ~1 if it isn't overdispersed.

However, when I examine the ZIP the residuals are evenly spread but the standard deviation is ~ 2...

So I ran a ZINB... and found the coefficients are almost identical between the two models.

I have two questions:

1) Am I correctly calculating the residuals/predictions to make such plots?

2) I've read that WAIC could be useful in selecting between a ZIP and a ZINB. I've read the paper by Aki and Andrew and used their function which shows the ZINB has a waic score of 5041 compared to the ZIPs 8422. What I am not 100% sure of is whether I am correctly calculating the log likelihood in both models?

Any help/advice would be greatly appreciated!

ZIP model

  data{

    int<lower=0> n_obs;

    int<lower=0> n_spp;

    int<lower=1> spp[n_obs];

    real cs_TPRE[n_obs];

    real cs_LDMC[n_obs];

    int TRT[n_obs];

    real TYPE[n_obs];

    int<lower=0> y[n_obs];

  }

parameters{

  real a0_raw[n_spp];

  real a0_mu;

  real<lower=0> a0_sigma;

  real a1;

  real a2;

  real a3;

  real a4;

  real b0_raw[n_spp];

  real b0_mu;

  real<lower=0> b0_sigma;

  real b1;

  real b2;

  real b3;

  real b4;

}

transformed parameters {

  real lambda[n_obs];

  real p0[n_obs];

  real a0[n_spp];

  real b0[n_spp];

  

  for (s in 1:n_spp) {

    a0[s] <- a0_raw[s] * a0_sigma + a0_mu;

    b0[s] <- b0_raw[s] * b0_sigma + b0_mu;

  }

  

  for (i in 1:n_obs) {

    lambda[i] <- exp(b0[spp[i]] + b1 * TRT[i] + b2 * cs_LDMC[i] + b3 * TYPE[i] + b4 * cs_TPRE[i]);

    p0[i] <- inv_logit(a0[spp[i]] + a1 * TRT[i] + a2 * cs_LDMC[i] + a3 * TYPE[i] + a4 * cs_TPRE[i]);

  }

}

model {    

  for (s in 1:n_spp) {

    a0_raw[s] ~ normal(0,1);

    b0_raw[s] ~ normal(0,1);

  }

  for (i in 1:n_obs){

    if (y[i] == 0) {

      increment_log_prob(log_sum_exp(bernoulli_log(1, p0[i]),

                                     bernoulli_log(0, p0[i])

                                     + poisson_log(y[i], lambda[i])));

    } else {

      increment_log_prob(bernoulli_log(0, p0[i])

                         + poisson_log(y[i], lambda[i]));

    }

  }

}

generated quantities {

  real pred[n_obs];

  real resid[n_obs];

  real log_lik[n_obs];

  for(i in 1:n_obs) {

    if (y[i]== 0) {

      log_lik[i] <- log_sum_exp(bernoulli_log(1, p0[i]),

                                bernoulli_log(0, p0[i])

                                + poisson_log(y[i], lambda[i]));}

    else {

      log_lik[i] <- bernoulli_log(0, p0[i]) + poisson_log(y[i], lambda[i]);

    }

   pred[i] <- (1-p0[i]) * lambda[i];

   resid[i] <- (y[i] - pred[i])/sqrt(pred[i]);

  }

}

ZINB model

  data{

    int<lower=0> n_obs;

    int<lower=0> n_spp;

    int<lower=1> spp[n_obs];

    real cs_TPRE[n_obs];

    real cs_LDMC[n_obs];

    int TRT[n_obs];

    real TYPE[n_obs];

    int<lower=0> y[n_obs];

  }

parameters{

  real a0_raw[n_spp];

  real a0_mu;

  real<lower=0> a0_sigma;

  real a1;

  real a2;

  real a3;

  real a4;

  real b0_raw[n_spp];

  real b0_mu;

  real<lower=0> b0_sigma;

  real b1;

  real b2;

  real b3;

  real b4;

  real<lower=0> phi;

}

transformed parameters {

  real lambda[n_obs];

  real p0[n_obs];

  real a0[n_spp];

  real b0[n_spp];

  

  for (s in 1:n_spp) {

    a0[s] <- a0_raw[s] * a0_sigma + a0_mu;

    b0[s] <- b0_raw[s] * b0_sigma + b0_mu;

  }

  

  for (i in 1:n_obs) {

    lambda[i] <- exp(b0[spp[i]] + b1 * TRT[i] + b2 * cs_LDMC[i] + b3 * TYPE[i] + b4 * cs_TPRE[i]);

    p0[i] <- inv_logit(a0[spp[i]] + a1 * TRT[i] + a2 * cs_LDMC[i] + a3 * TYPE[i] + a4 * cs_TPRE[i]);

  }

}

model {    

  for (s in 1:n_spp) {

    a0_raw[s] ~ normal(0,1);

    b0_raw[s] ~ normal(0,1);

  }

  for (i in 1:n_obs){

    if (y[i] == 0) {

      increment_log_prob(log_sum_exp(bernoulli_log(1, p0[i]),

                                     bernoulli_log(0, p0[i])

                                     + neg_binomial_2_log(y[i], lambda[i], phi)));

    } else {

      increment_log_prob(bernoulli_log(0, p0[i])

                         + neg_binomial_2_log(y[i], lambda[i], phi));

    }

  }

}

generated quantities {

  real pred[n_obs];

  real resid[n_obs];

  real log_lik[n_obs];

  for(i in 1:n_obs) {

    if (y[i]== 0) {

      log_lik[i] <- log_sum_exp(bernoulli_log(1, p0[i]),

                                bernoulli_log(0, p0[i])

                                + neg_binomial_2_log(y[i], lambda[i], phi));}

    else {

      log_lik[i] <- bernoulli_log(0, p0[i]) + neg_binomial_2_log(y[i], lambda[i], phi);

    }

   pred[i] <- (1-p0[i]) * lambda[i];

   resid[i] <- (y[i] - pred[i])/sqrt(pred[i] + phi * pred[i]^2);

  }

}

[James Camac]

Sorry, I've just realised I posted up the non hierarchical ZINB model. I've edited the post above to reflect the model I'm actually talking about.

[Aki Vehtari]

I couldn't spot a error. Since you didn't mention whether you used log likelihood scale or deviance scale, it's not clear which model is now estimated to be better. The difference in WAIC seems large, unless you have big data and phi is far away from values corresponding to Poisson model. How much data you have? How does the posterior of phi look like? Do you get same WAIC values, if you fix phi so that NegBin is close to Poisson?

Aki

[James Camac]

Thanks for your reply Aki,

I was referring to the deviance scale... and yes I thought it was largish as well.

I've had a look at the observed vs predicted and the model does a reasonable job at fitting the data. So I thought I must be calculating the likelihood incorrectly.

The exact estimates I get using the function in your paper are:

ZIP

      waic        lpd     p_waic  elpd_waic      p_loo   elpd_loo 

 8422.1959 -3700.7198   510.3781 -4211.0980   451.4138 -4152.1336 

ZINB

       waic         lpd      p_waic   elpd_waic       p_loo    elpd_loo 

 4596.69227 -2234.31099    64.03514 -2298.34614    64.81982 -2299.13082 

The data I'm fitting consists of 600 observations and n_spp = 40 (which is used for the varying intercepts).

phi looks reasonable to me. Here is phi's estimated posterior

   mean se_mean   sd 2.5% 25%  50%  75% 97.5% n_eff Rhat

phi 4.23    0.01 0.32 3.65   4 4.22 4.43  4.91   622 1.01

When I fix phi to be 10000 (i.e. minimising the extra variance in mu + (mu^2)/phi... I get a WAIC estimate similar to that found in the zip. The models are all converging nicely... so no issue there.

[Aki Vehtari]

You have 40+40=80 parameters. Usually p_waic and p_loo are smaller than the total number of parameters. With ZINB model this is the case, but with ZIP model p_waic and p_loo are so large that there is likely to be something wrong...

Aki

[James Camac]

Hmm but that would suggest the ZINB is possibly right... And if I can get the same waic values by setting phi as a large number then I can't see where the issue is coming from.
Could it simply be that the over dispersion in the zip model is driving these high values? After all that's really all that is changing. Or have a somehow missed something when calculating the log likelihood?

[Bob Carpenter]

Did you have a latent dispersion parameter and if so, did you
marginalize it out of the likelihood calculations?

- Bob

> On Jun 15, 2015, at 8:57 AM, James Camac <james...@gmail.com> wrote:
>
> Hmm but that would suggest the ZINB is possibly right... And if I can get the same waic values by setting phi as a large number then I can't see where the issue is coming from.
> Could it simply be that the over dispersion in the zip model is driving these high values? After all that's really all that is changing. Or have a somehow missed something when calculating the log likelihood?
>

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[James Camac]

Hi Bob,

The ZINB model contains phi (i.e.  neg_binomial_2_log(y[i], lambda[i], phi)) which accounts for over dispersion when observations are > 0.

I'm not marginalising this out - at least I don't think I am. (model code is in previous posts).

The ZIP model - the model Aki mentioned some concern about in terms of the WAIC estimates doesn't have a latent over dispersion parameter.

[Bob Carpenter]

> On Jun 17, 2015, at 7:26 PM, James Camac <james...@gmail.com> wrote:
>
> Hi Bob,
> The ZINB model contains phi (i.e. neg_binomial_2_log(y[i], lambda[i], phi)) which accounts for over dispersion when observations are > 0.
> I'm not marginalising this out - at least I don't think I am. (model code is in previous posts).

If you think of neg-binomial as a mixture of Poissons,
then neg_binomial_2 does the marginalization implicitly,
so you don't have to:

https://en.wikipedia.org/wiki/Negative_binomial_distribution#Gamma.E2.80.93Poisson_mixture

- Bob

[James Camac]

Ah ok. But how does that explain the issue associated with the WAIC estimate? 

I've set Phi as some very large number and managed to almost replicate what I found in the ZIP model... so I don't think it has much to do with the dispersion parameter. Or am I missing something obvious here.

[Bob Carpenter]

I don't know about these models, but if y's likelihood is
expressed as a normal mixture with latent boolean (discrete --- don't
try in Stan) indicator:

z[n] ~ bernoulli(theta);

y[n] ~ normal(mu[z[n]], sigma[z[n]]);

then the likelihood is usually just p(y | mu, z) and it needs to
integrate (sum) out the z[n], to get:

p(y[n] | mu, z) = theta * normal(y[n] | mu[1], sigma[1])
+ (1 - theta) * normal(y[n] | mu[2], sigma[2])

So if you have latent parameters per item, they should be marginalized
out for the likelihood.

- Bob

> On Jun 18, 2015, at 4:20 AM, James Camac <james...@gmail.com> wrote:
>
> Ah ok. But how does that explain the issue associated with the WAIC estimate?
> I've set Phi as some very large number and managed to almost replicate what I found in the ZIP model... so I don't think it has much to do with the dispersion parameter. Or am I missing something obvious here.
>
>
>

[James Camac]

Thanks Bob, I'll have a look.

Though the negative binomial (and hence the dispersion parameter) may not be the issue here. As Aki mentioned  the ZINB seems to produce 'reasonable' estimates of p_waic and p_loo but the zero inflate poisson - the model with no over dispersion parameter has very large values.

I'm really not sure why the WAIC calculations are off in the ZIP. but phi in the negative binomial appears to be important so the ZINB is probably the better model to use.

[Bob Carpenter]

Can you do some actual cross-validation on held-out data? That'll
make sure there's no fudging with latent parameters not being marginalized
out properly in the comparisons.

- Bob