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 4.6 Multivariate normal linear model in R using JAGS

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

require(R2jags)

set.seed(1056)                                                            # set seed to replicate example
nobs = 5000                                                               # number of obs in model

x1 <- runif(nobs)                                                       # random uniform variable
x2 <- runif(nobs)                                                       # random uniform variable
beta1 = 2.0                                                                 # intercept
beta2 = 3.0                                                                 # 1st coefficient
beta3 = -2.5                                                                # 2nd coefficient

xb <- beta1 + beta2*x1 + beta3*x2                           # linear predictor
y <- rnorm(nobs, xb, sd=1)                                        # create y as adjusted random normal variate

 

# Model setup
X <- model.matrix(~ 1 + x1+x2)
K <- ncol(X)

model.data <- list(Y = y,                                             # response variable
                              X = X,                                           # predictors
                              K = K,                                           # number of predictors including the intercept
                              N = nobs                                       # sample size
                              )

NORM <- "
model{
    # Diffuse normal priors for predictors
    for (i in 1:K) { beta[i] ~ dnorm(0, 0.0001) }

    # Uniform prior for standard deviation
    tau <- pow(sigma, -2)                                               # precision
    sigma ~ dunif(0, 100)                                              #  standard deviation

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

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

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

 

normOfit <- jags(data = model.data,
                             inits = inits,
                             parameters = params,
                             model = textConnection(NORM),
                             n.chains = 3,
                             n.iter = 15000,
                             n.thin = 1,
                             n.burnin = 10000)

 

print (normOfit, intervals=c(0.025, 0.975), digits=2)

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

Output on screen:

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

    each with 15000 iterations (first 10000 discarded)

    n.sims = 15000 iterations saved

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

beta[1]           2.01           0.04         1.94           2.09                 1            15000

beta[2]           3.03           0.05         2.93           3.13                 1            15000

beta[3]          -2.55           0.05       -2.65          -2.45                 1            15000

sigma             1.01          0.01          0.99           1.03                 1            11000

deviance 14317.14         2.84  14313.60   14324.26                 1             12000

 

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.0 and DIC = 14321.2

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