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”
```R
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](https://groups.google.com/forum/#!topic/stan-users/bSb2egD2-zY) and [here](https://andrewgelman.com/2012/11/16/stantastic/) but the code he provides seems to be outdated… and anyway I need to to know what I am doing wrong.
---
### Author:  [@matsutakehoyo](https://discourse.mc-stan.org/u/matsutakehoyo)
#### Post date: [November 30, 2018, 2:38am UTC](https://discourse.mc-stan.org/t/zero-inflated-gamma-model/6788/2 "2018-11-30T02:38:15Z")
…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.
So instead of this
```stan
// 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
```stan
// 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?
---
### Author:  [@Max\_Mantei](https://discourse.mc-stan.org/u/Max_Mantei)
#### Post date: [November 30, 2018, 11:38am UTC](https://discourse.mc-stan.org/t/zero-inflated-gamma-model/6788/3 "2018-11-30T11:38:51Z")
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.
---
### Author:  [@ldeschamps](https://discourse.mc-stan.org/u/ldeschamps)
#### Post date: [November 30, 2018, 12:31pm UTC](https://discourse.mc-stan.org/t/zero-inflated-gamma-model/6788/4 "2018-11-30T12:31:09Z")
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
---
### Author:  [@Max\_Mantei](https://discourse.mc-stan.org/u/Max_Mantei)
#### Post date: [November 30, 2018, 1:06pm UTC](https://discourse.mc-stan.org/t/zero-inflated-gamma-model/6788/5 "2018-11-30T13:06:52Z")
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
---
### Author:  [@ldeschamps](https://discourse.mc-stan.org/u/ldeschamps)
#### Post date: [November 30, 2018, 1:11pm UTC](https://discourse.mc-stan.org/t/zero-inflated-gamma-model/6788/6 "2018-11-30T13:11:15Z")
Thank you Max, you are perfectly right, and it was exactly what I was looking for!
Cheers,
Lucas
---
### Author:  [@matsutakehoyo](https://discourse.mc-stan.org/u/matsutakehoyo)
#### Post date: [December 1, 2018, 12:15am UTC](https://discourse.mc-stan.org/t/zero-inflated-gamma-model/6788/7 "2018-12-01T00:15:59Z")
@Max_Mantei thank you
I will try `gamma_lpdf(y | inverse(va), inverse(va)/mu)`, any tricks to speed up sampling are appreciated!
---
### Author:  [@Xwang1020](https://discourse.mc-stan.org/u/Xwang1020)
#### Post date: [January 22, 2019, 7:25am UTC](https://discourse.mc-stan.org/t/zero-inflated-gamma-model/6788/8 "2019-01-22T07:25:34Z")
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
---
### Author:  [@Max\_Mantei](https://discourse.mc-stan.org/u/Max_Mantei)
#### Post date: [January 22, 2019, 10:59am UTC](https://discourse.mc-stan.org/t/zero-inflated-gamma-model/6788/9 "2019-01-22T10:59:00Z")
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.
---
### Author:  [@matsutakehoyo](https://discourse.mc-stan.org/u/matsutakehoyo)
#### Post date: [January 23, 2019, 2:14am UTC](https://discourse.mc-stan.org/t/zero-inflated-gamma-model/6788/10 "2019-01-23T02:14:59Z")
@Xwang1020 in my original post I was referring to [Stan User’s Guide](https://mc-stan.org/users/documentation/). 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.
---
### Author:  [@Xwang1020](https://discourse.mc-stan.org/u/Xwang1020)
#### Post date: [January 23, 2019, 3:36am UTC](https://discourse.mc-stan.org/t/zero-inflated-gamma-model/6788/11 "2019-01-23T03:36:31Z")
Thanks a lot !
---
### Author:  [@Zhenkai\_Ran](https://discourse.mc-stan.org/u/Zhenkai_Ran)
#### Post date: [June 9, 2020, 7:16am UTC](https://discourse.mc-stan.org/t/zero-inflated-gamma-model/6788/12 "2020-06-09T07:16:31Z")
Do we need to assign a prior distribution for the Variance term?
---
### Author:  [@matsutakehoyo](https://discourse.mc-stan.org/u/matsutakehoyo)
#### Post date: [June 12, 2020, 1:54am UTC](https://discourse.mc-stan.org/t/zero-inflated-gamma-model/6788/13 "2020-06-12T01:54:46Z")
I believe Stan will just use a uniform distribution for prior if you don’t specify it.