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

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

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:

bottom of page