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