# Predicting from a Gamma Hurdle Model in JAGS

I fit a gamma hurdle model to invertebrate biomass data and am having trouble predicting from the model. I've been using these posts extensively to try to set it up:

Gamma hurdle model for continuous response?

Specify a Zero-inflated (Hurdle) Gamma Model in JAGS/BUGS

https://seananderson.ca/2014/05/18/gamma-hurdle/

but I can't quite figure out which coefficients and on what scale need to be multiplied together to get predictions. To complicate matters, I would very much like to predict the median biomass rather than the mean, which would presumably require some moment matching, and I'm not quite sure where that needs to happen.

Here's my model:

# JAGS model
model{
####### Priors #########
#intercept
a0~dnorm(0, 0.0001)

#habitat
ahab[1]<-0
thab[1]<-0
for(h in 2:nhab){
ahab[h]~dnorm(0, 0.0001)
thab[h] ~ dnorm(0, 0.5)
}

#week
aweek[1]<-0
tweek[1]<-0
for(w in 2:nweek){
aweek[w]~dnorm(0, 0.0001)
tweek[w] ~ dnorm(0, 0.5)
}

#prior for hab x week interaction
for(h in 1:nhab){
aint[h, 1]<-0
tint[h, 1]<-0
}
for(w in 2:nweek){
aint[1,w]<-0
tint[1,w]<-0
}
for(h in 2:nhab){
for(w in 2:nweek){
aint[h,w]~dnorm(0, 0.0001)
tint[h,w]~dnorm(0, 0.5)
}}

#year
ayear[1]<-0
ayear[2]~dnorm(0, 0.0001)

#temp, do, cond
a3~dnorm(0, 0.0001)
a4~dnorm(0, 0.0001)
a5~dnorm(0, 0.0001)

#wetland size
asize~dnorm(0, 0.0001)

sig_mu~dgamma(1,1)
tau_mu<-1/sig_site^2

#site random effect
for(s in 1:nsites){
asite[s]~dnorm(0, tau_site)
}

sig_site~dgamma(0.001,0.001)
tau_site<-1/sig_site^2

t0 ~ dnorm(0, 0.5)

tdry ~ dnorm(0, 0.5)

####### impute missing covariate data ########

for(i in 1:n){
temp[i] ~ dnorm(0, 1)
do[i] ~ dnorm(0, 1)
cond[i] ~ dnorm(0, 1)
}

# Likelihood using the zero trick
C <- 10000

for(i in 1:n){
zeroes[i] ~ dpois(-ll[i] + C)

# gamma log-likelihood
logGamma[i] <- log(dgamma(y[i], shape[i], rate[i]))

shape[i] <- pow(mu[i], 2) / pow(sig_mu, 2)
rate[i] <- mu[i] / pow(sig_mu, 2)

z[i] <- step(y[i] - 0.0001)

ll[i] <- (1 - z[i]) * log(1 - psi[i]) + z[i] * ( log(psi[i]) + logGamma[i])

log(mu[i]) <- a0 + asite[site[i]] + ahab[hab[i]] + aweek[week[i]] +
asize*size[i] +
aint[hab[i], week[i]] + a3*temp[i] + a4*do[i] + a5*cond[i] +
ayear[year[i]]

logit(psi[i]) <- t0 + thab[hab[i]] + tdry*dry[i] + tweek[week[i]] +
tint[hab[i], week[i]]
}

####### Derived Quantities ########

## predict invert energy for each habitat and week
for(i in 1:n.pred.i){
log(gamma_coef[i]) <- a0 + ahab[pred.ihab[i]] +
aweek[pred.iweek[i]] + aint[pred.ihab[i], pred.iweek[i]] +
ayear[pred.iyear[i]]# + asize*pred.isize[i]
logit(psi.pred[i]) <- t0 + thab[pred.ihab[i]] + tdry*pred.dry[i] +
tweek[pred.iweek[i]] + tint[pred.ihab[i], pred.iweek[i]]
#z.pred[i] ~ dbern(psi.pred[i])
pred.i[i]<-(psi.pred[i])*(gamma_coef[i])
}

}


What I'm wondering is:

1. Do I need to predict a 0 or 1 (z.pred that's commented out) and multiply that by gamma_coef[i]
2. Do I simply multiple psi.pred[i]*gamma_coef[i] like I have
3. Do I need to change gamma_coef[i] in some way to make sure it's only the positive mass of the gamma distribution that's getting multiplied by psi.pred[i]
4. And finally, where would I do the moment matching if I want pred.i[i] to be the median of each habitat x week rather than the mean?