1
$\begingroup$

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?
$\endgroup$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.