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

# Data from code 10.14

require(R2jags)
require(jagstools)


# Data
path_to_data = "https://raw.githubusercontent.com/astrobayes/BMAD/master/data/Section_10p7/GCs.csv"

# Read data
GC_dat = read.csv(file=path_to_data,header = T,dec=".")

# Prepare data to JAGS
N <- nrow(GC_dat)
x <- GC_dat$MV_T
y <- GC_dat$N_GC
X <- model.matrix(~ x, data=GC_dat)
K = ncol(X)

JAGS_data <- list(
  Y = y,
  X = X,
  N = N,
  K = K)

Code 10.16 NB-P model in R using JAGS, for modeling the relationship between globular cluster population and host galaxy visual magnitude

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

# Fit
model.NBP <- "model{
    # Priors for regression coefficients
    # Diffuse normal priors betas
    for (i in 1:K) { beta[i] ~ dnorm(0, 1e-5)}

    # Prior for size
    theta ~ dgamma(1e-3,1e-3)
    
    # Uniform prior for Q
    Q ~ dunif(0,3)

    # Likelihood
    for (i in 1:N){
        eta[i]<-inprod(beta[], X[i,])
        mu[i] <- exp(eta[i])
        theta_eff[i]<- theta*(mu[i]^Q)
        p[i]<-theta_eff[i]/(theta_eff[i]+mu[i])
        Y[i]~dnegbin(p[i],theta_eff[i])

        # Discrepancy
        expY[i] <- mu[i] # mean
        varY[i] <- mu[i] + pow(mu[i],2-Q)/theta #variance
        PRes[i] <- ((Y[i] - expY[i])/sqrt(varY[i]))^2
    }
    Dispersion <- sum(PRes)/(N-4)#
}"

# Set initial values
inits <- function () {
    list(
        beta = rnorm(K, 0, 0.1),
        theta = runif(1,0.1,5),
        Q = runif(1,0,1)
    )
}

# Identify parameters
params <- c("Q","beta","theta","Dispersion")

# Start JAGS
NBP_fit <- jags(data = JAGS_data,
                inits = inits,
                parameters = params,
                model = textConnection(model.NBP),
                n.thin = 1,
                n.chains = 3,
                n.burnin = 5000,
                n.iter = 20000)

# Output
# Plot posteriors

MyBUGSHist(out,c("Dispersion",uNames("beta",K),"theta"))

# Screen output
print(NBP_fit, intervals=c(0.025, 0.975), digits=3)

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

 

Output on screen:

Inference for Bugs model at "3", fit using jags,

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

  n.sims = 45000 iterations saved

 

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

Dispersion       1.925       0.208          1.547           2.355     1.002       1700

Q                     0.018        0.015         0.001           0.056     1.001        5500

beta[1]         -11.797        0.328      -12.454        -11.172     1.012          180

beta[2]           -0.883        0.016        -0.916          -0.852     1.012          180

theta                0.996        0.104         0.775           1.184      1.001      32000

deviance    5193.038       2.920    5189.225     5200.284      1.001      45000

 

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.3 and DIC = 5197.3

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

© 2017 by Emille E. O. Ishida