# Zero Inflated Gamma model

I am trying to fit a Zero Inflated (or Hurdle) Gamma model, following the example in section 13.7 on the manual. I use mean and variance to parameterize the Gamma distribution because they are more intuitive for interpreting results, but I don’t think this is relevant. Almost the exact model but with the poisson_lpmg instead of gamma_lpdf works fine, what am I missing? Am I not supposed to mix lpmf and Ldp.? Any suggestions would be appreciated!

// zero inflated gamma
data {
int<lower=1> n_obs;
real y_mean;
real y_vari;
real non_zero;
real<lower=0>  y[n_obs];
}

parameters {
real t;
real m;
real v;
}

transformed parameters {
real<lower=0, upper=1> theta; // probability of zero
real<lower=0> mu; // mean
real<lower=0> va; // variance

theta = inv_logit(t);
mu = exp(m);
va = exp(v);
}

model {
t ~ normal(logit(non_zero), 1);
m ~ normal(log(y_mean),1);
v ~ normal(log(y_vari),1);

// likelihood: Zero_Inflated Gamma
for (i in 1:n_obs) {
if (y[i] == 0) {
target += log_sum_exp(bernoulli_lpmf(1 | theta),
bernoulli_lpmf(0 | theta) +
gamma_lpdf(y[i] | mu*mu/va, mu/va));
} else {
target += bernoulli_lpmf(0 | theta) +
gamma_lpdf(y[i] | mu*mu/va, mu/va);
}
}
}


I made up some data and run the model but I can not make it work. I get “… chains where the estimated Bayesian Fraction of Missing Information was low”

library(tidyverse)
library(rstan)

# total sample number
n <- 2000
# probability of zero
theta <- 0.3
# gamma pars, these are close to the values in my actual data
mu = exp(-4.5); va = exp(-11)
shape = mu^2/va; rate = mu/va

df <- tibble(y=rgamma(n*(1-theta), shape, rate)) %>%
bind_rows(tibble(y=rep(0, n*theta)))

p <- ggplot(df, aes(x=y)) +
geom_histogram()

model_data <- list(
n_obs = nrow(df),
y_mean = df %>% filter(y>0) %>% pull(y) %>% mean(),
y_vari = df %>% filter(y>0) %>% pull(y) %>% var(),
non_zero = df %>% filter(y==0) %>% nrow()/nrow(df),
y=df$y ) model_name <- "zig simulation.stan" model <- stan(file=model_name, data = model_data, cores = getOption("mc.cores", 1L)) summary(model)$summary



p.s. looking around the internet I found some discussion on this (Zero Inflated Gamma) by Richard McElreath here and here but the code he provides seems to be outdated… and anyway I need to to know what I am doing wrong.

1 Like

…ok I now realize my folly ಥ_ಥ
This isn’t really a zero inflated model, in the sense that the zeros only come from the bernoulli process, as the gamma does not produce any zeros.

// likelihood: Zero_Inflated Gamma
for (i in 1:n_obs) {
if (y[i] == 0) {
target += log_sum_exp(bernoulli_lpmf(1 | theta),
bernoulli_lpmf(0 | theta) +
gamma_lpdf(y[i] | mu*mu/va, mu/va));
} else {
target += bernoulli_lpmf(0 | theta) +
gamma_lpdf(y[i] | mu*mu/va, mu/va);
}
}


I used this

// likelihood: Gamma hurdle
for (i in 1:n_obs) {
if (y[i] == 0) {
target += bernoulli_lpmf(1 | theta);
} else {
target += bernoulli_lpmf(0 | theta) +
gamma_lpdf(y[i] | mu*mu/va, mu/va);
}
}


and that seem to work just fine!!
I think this is essentially fitting two models, bernoulli for zeros and gamma for the rest. Is there an advantage or trade-off to doing it this, or fitting two models separately?

1 Like

When you do two separate models it’s easier to vectorize… if speed is an issue. Also, you could try gamma_lpdf(y | inverse(va), inverse(va)/mu), which tends to be more stable. Even better if you have something like inverse_va in your parameter block and “sample” it directly.

3 Likes

Hello Max!

I use regularly such parametrization of the gamma distribution, and I’m very interested by your comment.

Could you please explain how you get this formulation and conclusion? Or could you suggest a reference? :)

Thank you!
Lucas

Hi Lucas,

I think I found this on the old Stan Google mailing list. I tried it and it tended to work better for my applications.

But it’s actually pretty easy to derive this formulation from the perspective of a GLM framework, i.e. that the variance function is V(\mu)=\mu^2 for the Gamma GLM. With dispersion parameter \phi the variance is given by \phi V(\mu)=\phi \mu^2. Then plug in a/b for the mean \mu and a/b^2 for the variance, i.e. a/b^2 = \phi (a/b)^2, and solve for a and b. Then you should end up with something like a=\phi^{-1} and b=\phi^{-1} / \mu, where a and b are of course the parameters of the Gamma distribution as implemented in Stan.

I’m not entirely sure if that’s what you’ve been looking for, but I hope it helps.

Cheers, Max

3 Likes

Thank you Max, you are perfectly right, and it was exactly what I was looking for!

Cheers,
Lucas

@Max_Mantei thank you
I will try gamma_lpdf(y | inverse(va), inverse(va)/mu), any tricks to speed up sampling are appreciated!

Hi there,

What manual are you talking about here, I am currently doing some research work on zero-inflated gamma model, would you mind sharing the data from the example you metioned? I really appreciate it.

Xiao

Hello Xiao,

the Gamma GLM can not really be zero-inflated - because it has zero probability density at zero. You are probably looking for a hurdle model (sometimes called two-part model, because you could estimate each part separately). In the Stan manual (on the Stan website) there is a short discussion of a Poisson hurdle model, which is even more complicated, because of the truncation of the Poisson at zero. You won’t have that in a Gamma hurdle model.

@Xwang1020 in my original post I was referring to Stan User’s Guide. There is an example with code for a Zero Inflated Poisson model. I used the previous versions of the manual so my chapter reference might be a bit off.

Thanks a lot !

Do we need to assign a prior distribution for the Variance term?

I believe Stan will just use a uniform distribution for prior if you don’t specify it.