Tommaso Guerrini - 5 months ago 51

R Question

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)`

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`

Answer

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.