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)


# 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:

​

​

bottom of page