2
votes

I'm attempting to do this in JAGS:

z[l] ~ dbeta(0.5,0.5)
y[i,l] ~ z[l]*dnorm(0,10000) + inprod(1-z[l],dnegbin(exp(eta_star[i,l]),alpha[l]))

(dnorm(0,10000) models a Dirac delta in 0: see here if you are interested in the model).

But I get:

     RUNTIME ERROR:
     Incorrect number of arguments in function dnegbin

But if I do this:

y[i,l] ~ dnegbin(exp(eta_star[i,l]),alpha[l])

It runs just fine. I wonder that I cannot multiply a value for a distribution, so I imagine that something like this could work:

z[l] ~ dbeta(0.5,0.5)
pointmass_0[l] ~ dnorm(0,10000)
y[i,l] ~ dnegbin(exp(eta_star[i,l]),alpha[l])
y_star[i,l] = z[l]*pointmass_0[l]+inprod(1-z[l],y[i,l])

If I run that I get:

ystar[1,1] is a logical node and cannot be observed
2
if by the y* in your third line you mean to introduce a new variable (I don't think y* will work, try y2 or ystar?) then yes, you have to do it that way. You can't just mix stochastic nodes into computations in that way. - Ben Bolker
@BenBolker Thanks, I edited the y* since it's not available in JAGS. Is this the only way so? At that point my response variable y_star becomes a deterministic node - Tommaso Guerrini
As far as I know, yes. - Ben Bolker
But at that point I'd get: 'ystar[2,1] is a logical node and cannot be observed' - Tommaso Guerrini

2 Answers

5
votes

You are looking to model a zero-inflated negative binomial model. You can do this in JAGS if you use the "ones trick", an pseudo-likelihood method that can be used when the distribution of your outcome variables is not one of the standard distributions in JAGS but you can still write down an expression for the likelihood.

The "ones trick" consists of creating pseudo-observations with the value 1. These are then modeled as Bernoulli random variables probability parameter Lik/C where Lik is the likelihood of your observations and C is a large constant to ensure that Lik/C << 1.

data {
   C <- 10000
   for (i in 1:N) {
      one[i,1] <- 1
   }
}
model {

   for (i in 1:N) {
      one[i,1] ~ dbern(lik[i,1]/C)

      lik[i,1] <- (y[i,1]==0)*z[1] + (1 - z[1]) * lik.NB[i,1]   
      lik.NB[i,1] <- dnegbin(y[i,1], exp(eta_star[i,1]), alpha[1])
   }

   z[l] ~ dbeta(0.5,0.5)
}

Note that the name dnegbin is overloaded in JAGS. There is a distribution that has two parameters and a function that takes three arguments and returns the likelihood. We are using the latter.

I am thinking of adding zero-inflated versions of count distributions to JAGS, since the above construction is quote awkward for the user, whereas zero-inflated distributions are quite easy to implement internally in JAGS.

2
votes

I too would like to know a better way to handle this situation.

One cheesy solution is to add a stochastic node

ystarstar[i,j] ~ dnorm(ystar[i,j],10000000) 

(i.e. a Normal distribution with a very high precision, or a Dirac delta in your terminology) to the model.