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.3 Zero-inflated negative binomial synthetic data in R

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

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

set.seed(141)
nobs <- 1000

 

x1 <- runif(nobs)
x2 <- rbinom(nobs, size=1, 0.6)

xb <- 1 + 2.0*x1 + 1.5*x2
xc <- 2 - 5.0*x1 + 3*x2

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

 

zinby <- rzinegbin(n=nobs, munb = exb, size=1/alpha, pstr0=exc)
zinbdata <- data.frame(zinby,x1,x2)


Code 7.4  Bayesian zero-inflated negative binomial model  using JAGS

=========================================================
Xc <- model.matrix(~ 1 + x1+x2, data=zinbdata)
Xb <- model.matrix(~ 1 + x1+x2, data=zinbdata)

Kc <- ncol(Xc)
Kb <- ncol(Xb)


model.data <- list(Y = zinbdata$zinby,                                                   # response
                              Xc = Xc,                                                                    # covariates
                              Kc = Kc,                                                                    # number of betas
                              Xb = Xb,                                                                   # covariates
                              Kb = Kb,                                                                   # number of gammas
                              N = nrow(zinbdata))

ZINB<-"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)}

    alpha ~ dunif(0.001, 5) 

 

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


W <- zinbdata$zinby
W[zinbdata$zinby > 0] <- 1

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

 

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

 

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

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

 

# Figures for parameter trace plots and histogramss
source("https://raw.githubusercontent.com/astrobayes/BMAD/master/auxiliar_functions/CH-Figures.R")

 

out <- ZINB$BUGSoutput
MyBUGSHist(out,c(uNames("beta",Kc),"alpha",uNames("gamma",Kb)))
MyBUGSChains(out,c(uNames("beta",Kc),"alpha",uNames("gamma",Kb)))

 

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

 

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

alpha                     2.151        0.255            1.714             2.709          1.002      3000

beta[1]                  1.023        0.267            0.487             1.541          1.026           85

beta[2]                  1.967        0.356            1.309             2.700          1.036           65

beta[3]                  1.753        0.191            1.398             2.126          1.002       1400

gamma[1]             1.825        0.299            1.248             2.405          1.006       2200

gamma[2]            -4.736        0.505          -5.859            -3.782          1.034         640

gamma[3]             2.821        0.295            2.296             3.451          1.013         320

deviance         2550.517      41.743       2471.993       2635.648         1.002       3000

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 = 871.2 and DIC = 3421.7

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