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 10.3 Gaussian linear mixed model, in R using JAGS, for modeling the relationship between type Ia supernovae host galaxy mass and Hubble residuals

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

library(R2jags)

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

dat <- read.csv(path_to_data, header = T)

# Prepare data to JAGS
nobs = nrow(dat)
obsx1 <- dat$LogMass
errx1 <- dat$e_LogMass
obsy <- dat$HR
erry <- dat$e_HR


type <- as.numeric(dat$Type)                                     # convert class to numeric flag 1 or 2

jags_data <- list(
         obsx1 = obsx1,
         obsy = obsy,
         errx1 = errx1,
         erry = erry,
         K = 2,
         N = nobs,
         type = type)

# Fit
NORM_errors <-" model{
    tau0 ~ dunif(1e-1,5)
    mu0 ~ dnorm(0,1)

    # Diffuse normal priors for predictors
    for (j in 1:2){
        for (i in 1:K) {
            beta[i,j] ~  dnorm(mu0, tau0)
        }
    }

    # Gamma prior for standard deviation
    tau ~ dgamma(1e-3, 1e-3)                                         # precision
    sigma <- 1 / sqrt(tau)                                                 # standard deviation

    # Diffuse normal priors for true x
    for (i in 1:N){
        x1[i] ~ dnorm(0,1e-3)
    }

    # Likelihood function
    for (i in 1:N){
        obsy[i] ~ dnorm(y[i],pow(erry[i],-2))
        y[i] ~ dnorm(mu[i],tau)
        obsx1[i] ~ dnorm(x1[i],pow(errx1[i],-2))

        mu[i] <- beta[1,type[i]] + beta[2,type[i]] * x1[i]
    }
}"

inits <- function () {
    list(beta = matrix(rnorm(4, 0, 0.01),ncol = 2))
}

params0 <- c("beta", "sigma")

# Run MCMC
NORM <- jags(
             data = jags_data,
             inits = inits,
             parameters = params0,
             model = textConnection(NORM_errors),
             n.chains = 3,
             n.iter = 40000,
             n.thin = 1,
             n.burnin = 15000)

# Output
print(NORM,justify = "left", digits=3)

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

Output on screen:

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

3 chains, each with 40000 iterations (first 15000 discarded)

n.sims = 75000 iterations saved

 

                      mu.vect     sd.vect                2.5%            25%               50%             75%              97.5%   Rhat      n.eff

beta[1,1]            0.880       0.227              0.418            0.723             0.876           1.048               1.302    1.029     160

beta[2,1]           -0.085       0.022             -0.125          -0.101            -0.085          -0.070             -0.041    1.013     230

beta[1,2]            0.221        0.181            -0.159           0.110              0.229           0.339              0.574     1.032       68

beta[2,2]           -0.021       0.017            -0.055           -0.033            -0.022          -0.011              0.015    1.032        68

sigma                 0.120       0.009             0.103             0.114             0.120            0.126             0.138     1.001 18000

deviance     -1103.264     36.488      -1173.317      -1128.154      -1103.712     -1078.989     -1030.542     1.001 23000

 

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 = 665.6 and DIC = -437.6

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

© 2017 by Emille E. O. Ishida