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.5 Multivariate normal model in R using JAGS for accessing the relationship between period, luminosity, and color in early-type contact binaries

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

library(R2jags)

# Data
PLC <- read.csv("https://raw.githubusercontent.com/astrobayes/BMAD/master/data/Section_10p3/PLC.csv", header = T)

# Prepare data for JAGS
nobs = nrow(PLC)                                                                  # number of data points
x1 <- PLC$logP                                                                      # log period
x2 <- PLC$V_I                                                                       # V-I color
y <- PLC$M_V                                                                       # V magnitude
type <- as.numeric(PLC$type)                                               # type NC/GC
X <- model.matrix(~ 1 + x1+x2)                                            # covariate matrix
K <- ncol(X)                                                                           # number of covariates per type

jags_data <- list(Y = y,
                           X = X,
                           K = K,
                           type = type,
                            N = nobs)

# Fit
NORM <-"model{
    # Shared hyperprior
    tau0 ~ dgamma(0.001,0.001)
    mu0 ~ dnorm(0,1e-3)

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

    # Uniform prior for standard deviation
    for(i in 1:2) {
        tau[i] <- pow(sigma[i], -2)                         # precision
        sigma[i] ~ dgamma(1e-3, 1e-3)                 # standard deviation
    }

    # Likelihood function
    for (i in 1:N){
        Y[i]~dnorm(mu[i],tau[type[i]])
        mu[i] <- eta[i]
        eta[i] <- beta[1, type[i]] * X[i, 1] + beta[2, type[i]] * X[i, 2] +
        beta[3, type[i]] * X[i, 3]
    }
}"

# Determine initial values
inits <- function () {
  list(beta = matrix(rnorm(6,0, 0.01),ncol=2))
}

# Identify parameters
params <- c("beta", "sigma")

# Fit
jagsfit <- jags(data = jags_data,
                        inits = inits,
                        parameters = params,
                        model = textConnection(NORM),
                        n.chains = 3,
                        n.iter = 5000,
                        n.thin = 1,
                       n.burnin = 2500)

# Output
print(jagsfit,justify = "left", digits=2)

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

Output screen:

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

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

n.sims = 7500 iterations saved

 

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

beta[1,1]    -1.01         0.27    -1.55      -1.18     -1.00      -0.83      -0.47             1       7500

beta[2,1]    -3.31         0.97    -5.21      -3.94     -3.30      -2.67      -1.41             1       7500

beta[3,1]     7.24         1.29     4.68        6.39      7.26       8.10        9.73             1       7500

beta[1,2]     -0.41        0.16    -0.73      -0.52     -0.42      -0.31      -0.10             1       7500

beta[2,2]     -3.19        0.58    -4.31      -3.58     -3.19      -2.81       -2.01            1       7500

beta[3,2]      8.48        0.82     6.85       7.94       8.48       9.03       10.08            1       6900

sigma[1]      0.62        0.09     0.47      0.55        0.61       0.67         0.82            1       1500

sigma[2]      0.43        0.05     0.34      0.39        0.42       0.46         0.55            1       4300

deviance    91.96       4.34      85.5     88.78     91.27     94.39     102.15            1       7300

 

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 = 9.4 and DIC = 101.4

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