HSI
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:
​
​