Posterior predictive checking

Modeling
13

Hey! A few thoughts…

Warm-up?

Hm, I guess these values only occur during warm-up, right? Because in the plot, that you posted, y_rep looks fine (and those come from post-warmup-up iterations). You could try running stan with the save_warmup = FALSE option and see if you still get that error.

Priors?

Also, you could tweak your priors a bit, especially those of the NegBinom. For example, cauchy(0,5) can have extremely large values and having this as a prior for overdispersion is perhaps overkill. You could try

...
parameters {
  ...
  real log_phi[N];
}
...
model{
  ...
  log_phi ~ normal(0, 1);
  ...
    ... neg_binomial_2_lpmf(y[i]|mu[1,i], exp(log_phi[i])) ...
}

this way you have the prior median for \phi at 1, which implies “no over- or under-dispersion”, but have some wiggle-room to go above or below, putting equal prior probability on over- and under-dispersion.

(Also, just as an aside… Letting \phi vary by observation strikes me as unusual. It’s not that it couldn’t be reasonable, but I think normally you’d just have one parameter \phi that measures the dispersion of the whole distribution—like the \sigma in linear regression. Maybe having just one common \phi will give you a bit more stable inference?)

In general, you’d want to do proper prior predictive checks. There’s a lot of good resources on that for you to read, if your not familiar with it. It boils down to do the same as with posterior predictive checks, but using only the prior.

RNGs?

If you still get the error, then you could try something like this (look at this thread for reference):

      real mu_hat = which == 1 ? mu[i]*exp(ci) : mu[i];
      if (log(mu_hat) > 20){ // set some reasonable threshold here.. like a really big number
        y_hat[i] = 9999999; // set something here, that you want to put into y_hat in case mu is really large. Or some value that you can "filter" out
      } else {
        y_hat[i] = 2 + neg_binomial_2_rng(mu_hat, phi[i]);
      }

… this is super hacky (and for the values I’ve put in probably even a bit ridiculous), but I think you get the idea.

Likelihood?

However, you have said, that the maximum you’d expect y to be 11. So, basically the outcome is an integer y \in [2,11]? Are values beyond 11 possible? If no, then I think you should have a look into modeling your data using a Binomial distribution. To get random integer y between (and including) 2 and 11, you could try in the model block binomial_logit(y_ast | 9, alpha + x*beta) and so on. Further down, in generated quantities it would be something like y_hat[i] = 2 + binomial_rng(9, inv_logit(alpha + x[i]*beta)). Note, that now your coefficients would be on the logit-scale, so you’d have to change your priors accordingly (like I said above, look into prior predictive checks if you haven’t already, it’s a great way of coming up with reasonable priors).

In case you want to incorporate under- or over-dispersion, you could try the Beta-Binomial distribution. But here, you’d have to re-parametrize yourself.

I hope the suggestions are helpful!

Cheers! :)

1 Like