top of page

From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press

(c) 2017,  Joseph M. Hilbe, Rafael S. de Souza and Emille E. O. Ishida

you are kindly asked to include the complete citation if you used this material in a publication

â€‹

â€‹

Code 7.1 Bayesian zero-inflated Poisson model in R using JAGS

====================================================

require(MASS)
require(R2jags)
require(VGAM)

â€‹

set.seed(141)
nobs <- 1000

â€‹

x1 <- runif(nobs)
xb <- 1 + 2.0*x1
xc <- 2 - 5.0*x1

exb <- exp(xb)
exc <- 1/(1+exp(-xc))

zipy <- rzipois(n=nobs, lambda=exb, pstr0=exc)
zipdata <- data.frame(zipy,x1)

â€‹

Xc <- model.matrix(~ 1 + x1, data=zipdata)
Xb <- model.matrix(~ 1 + x1, data=zipdata)
Kc <- ncol(Xc)
Kb <- ncol(Xb)

â€‹

model.data <- list(Y = zipdata\$zipy,                      # response
Xc = Xc,                                   # covariates
Kc = Kc,                                   # number of betas
Xb = Xb,                                   # covariates
Kb = Kb,                                   # number of gammas
N = nobs)

â€‹

ZIPOIS<-"model{
# Priors - count and binary components
for (i in 1:Kc) { beta[i] ~ dnorm(0, 0.0001)}
for (i in 1:Kb) { gamma[i] ~ dnorm(0, 0.0001)}

â€‹

# Likelihood
for (i in 1:N) {
W[i] ~ dbern(1 - Pi[i])
Y[i] ~ dpois(W[i] * mu[i])
log(mu[i]) <- inprod(beta[], Xc[i,])
logit(Pi[i]) <- inprod(gamma[], Xb[i,])
}
}"

â€‹

W <- zipdata\$zipy
W[zipdata\$zipy > 0] <- 1

â€‹

inits <- function() {
list(beta = rnorm(Kc, 0, 0.1),
gamma = rnorm(Kb, 0, 0.1),
W = W)}

â€‹

params <- c("beta", "gamma")

â€‹

ZIP <- jags(data = model.data,
inits = inits,
parameters = params,
model = textConnection(ZIPOIS),
n.thin = 1,
n.chains = 3,
n.burnin = 4000,
n.iter = 5000)

â€‹

print(ZIP, intervals = c(0.025, 0.975), digits=3)

====================================================

Output on screen:

â€‹

Inference for Bugs model at "3", fit using jags, 3 chains,

each with 5000 iterations (first 4000 discarded)

n.sims = 3000 iterations saved

mu.vect     sd.vect              2.5%            97.5%          Rhat       n.eff

beta[1]               1.019       0.045            0.932              1.104         1.103          25

beta[2]               1.955       0.058            1.843              2.066         1.096          26

gamma[1]          1.913       0.166            1.588              2.241         1.013        160

gamma[2]         -4.828       0.319           -5.472            -4.208         1.014        150

deviance       3095.598     11.563      3079.409        3122.313         1.004        710

For each parameter, n.eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)

pD = 66.7 and DIC = 3162.3

DIC is an estimate of expected predictive error (lower deviance is better).

bottom of page