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.13 Bayesian lognormal–logit hurdle using JAGS

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

library(R2jags)

set.seed(33559)


# Sample size
nobs <- 1000 


# Generate predictors, design matrix
x1 <- runif(nobs,0,2.5)
xc <- 0.6 + 1.25*x1
y <- rlnorm(nobs, xc, sdlog=0.4)
lndata <- data.frame(y, x1)

# Construct filter
xb <- -3 + 4.5*x1
pi <- 1/(1+exp(-(xb)))
bern <- rbinom(nobs,size=1, prob=pi)


# Add structural zeros
lndata$y <- lndata$y*bern
Xc <- model.matrix(˜ 1 + x1,data = lndata)
Xb <- model.matrix(˜ 1 + x1,data = lndata)
Kc <- ncol(Xc)
Kb <- ncol(Xb)


JAGS.data <- list(
    Y = lndata$y,                                                            # response
    Xc = Xc,                                                                   # covariates
    Xb = Xb,                                                                  # covariates
    Kc = Kc,                                                                  # number of betas
    Kb = Kb,                                                                  # number of gammas
    N = nrow(lndata),                                                    # sample size
    Zeros = rep(0, nrow(lndata)))

 

load.module(’glm’)
 

sink("ZALN.txt")
 

cat("
model{
    # Priors for both beta and gamma components
    for (i in 1:Kc) {beta[i] ˜ dnorm(0, 0.0001)}
    for (i in 1:Kb) {gamma[i] ˜ dnorm(0, 0.0001)}

 

    # Prior for sigma
    sigmaLN ˜ dgamma(1e-3, 1e-3)

 

    # Likelihood using the zero trick
    C <- 10000
   

    for (i in 1:N) {
        Zeros[i] ˜ dpois(-ll[i] + C)
       

       # LN log-likelihood
       ln1[i] <- -(log(Y[i]) + log(sigmaLN) + log(sqrt(2 * sigmaLN)))
       ln2[i] <- -0.5 * pow((log(Y[i]) - mu[i]),2)/(sigmaLN * sigmaLN)
       LN[i] <- ln1[i] + ln2[i]

       z[i] <- step(Y[i] - 1e-5)
       l1[i] <- (1 - z[i]) * log(1 - Pi[i])
       l2[i] <- z[i] * ( log(Pi[i]) + LN[i])
       ll[i] <- l1[i] + l2[i]
       mu[i] <- inprod(beta[], Xc[i,])
       logit(Pi[i]) <- inprod(gamma[], Xb[i,])
    }
}"
, fill = TRUE)

 

sink()
 

# Initial parameter values
inits <- function () {
    list(beta = rnorm(Kc, 0, 0.1),
          gamma = rnorm(Kb, 0, 0.1),
          sigmaLN = runif(1, 0, 10))}
 

# Parameter values to be displayed in output
params <- c("beta", "gamma", "sigmaLN")

 

# MCMC sampling
ZALN <- jags(data = JAGS.data,
                        inits = inits,
                        parameters = params,
                        model = "ZALN.txt",
                        n.thin = 1,
                        n.chains = 3,
                        n.burnin = 2500,
                        n.iter = 5000)

 

 

# Model results
print(ZALN, intervals = c(0.025, 0.975), digits=3)

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

Output on screen: