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 5.6   Log-gamma synthetic data generated in R

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

set.seed(33559)


nobs <- 1000
r <- 20                                                                            # shape
beta1 <- 1
beta2 <- 0.66
beta3 <- -1.25


x1 <- runif(nobs)
x2 <- runif(nobs)


xb <- beta1 + beta2 * x1 + beta3 * x2


exb <- exp(xb)
py <- rgamma(nobs,shape = r, rate= r/exb)
LG <- data.frame(py, x1, x2)

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

Code 5.7  Log-gamma model in R using JAGS

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

library(R2jags)

X <- model.matrix(~ x1 + x2, data =LG)
K <- ncol(X)                                                                           # number of columns
model.data <- list(Y = LG$py,                                               # response
                             X = X,                                                       # covariates
                             N = nrow(LG),                                          # sample size
                             b0 = rep(0,K),
                            B0 = diag(0.0001, K))

sink("LGAMMA.txt")

cat("
model{
    # Diffuse priors for model betas
    beta ~ dmnorm(b0[], B0[,])

    # Diffuse prior for shape parameter
    r ~ dgamma(0.01, 0.01)
    
    # Likelihood
    C <- 10000
    for (i in 1:N){
        Y[i] ~ dgamma(r, lambda[i])
        lambda[i] <- r / mu[i]
        log(mu[i]) <- eta[i]
        eta[i] <- inprod(beta[], X[i,])
    }
}
"
,fill = TRUE)

sink()

inits <- function () {
  list(
    beta = rnorm(K,0,0.01),
    r = 1 )
}

params <- c("beta", "r")

# JAGS MCMC
LGAM <- jags(data = model.data,
                          inits = inits,
                          parameters = params,
                          model.file = "LGAMMA.txt",
                          n.thin = 1,
                          n.chains = 3,
                          n.burnin = 3000,
                          n.iter = 5000)

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

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

 
GET SOURCE

Output on screen:

Inference for Bugs model at "LGAMMA.txt", fit using jags,

    3 chains, each with 5000 iterations (first 3000 discarded)

    n.sims = 6000 iterations saved

 

                      mu.vect        sd.vect            2.5%          97.5%       Rhat       n.eff

beta[1]               1.010          0.020           0.972           1.050      1.051          47

beta[2]               0.659          0.025           0.609           0.707      1.015        140

beta[3]              -1.256          0.027         -1.312          -1.199      1.045          53

r                       19.669          0.883         17.984         21.457      1.001      3100

deviance      1216.067          2.917     1212.393     1223.227      1.009        590

 

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 = 4.2 and DIC = 1220.3

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