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 5.20

set.seed(135)
nobs <- 1:2000

x1 <- runif(nobs)
x2 <- 2*runif(nobs)

xb <- 2 + 0.75 * x1 - 1.25 * x2

exb <- pnorm(xb)                                                # probit inverse link


py <- rbinom(nobs, size=1, prob=exb)

probdata <- data.frame(py, x1, x2)

Code 5.22 Probit model in R using JAGS

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

attach(probdata)
require(R2jags)

set.seed(1944)
X <- model.matrix(~ x1 + x2,
                               data = probdata)


K <- ncol(X)


model.data <- list(Y = probdata$py,
                              N = nrow(probdata),
                              X =X,
                              K =K,
                              LogN = log(nrow(probdata)),
                              b0 = rep(0, K),
                              B0 = diag(0.00001, K)
)
sink("PROBIT.txt")


cat("
model{
    # Priors
    beta ~ dmnorm(b0[], B0[,])

    # Likelihood
    for (i in 1:N){
        Y[i] ~ dbern(p[i])
        probit(p[i]) <- max(-20, min(20, eta[i]))
        eta[i] <- inprod(beta[], X[i,])
        LLi[i] <- Y[i] * log(p[i]) +
        (1 - Y[i]) * log(1 - p[i])
    }

    LogL <- sum(LLi[1:N])

    # Information criteria
    AIC <- -2 * LogL + 2 * K
    BIC <- -2 * LogL + LogN * K
    }"
, fill = TRUE)

sink()

# Initial parameter values
inits <- function () {list(beta = rnorm(K, 0, 0.1)) }

# parameters and statistics to display
params <- c("beta", "LogL", "AIC", "BIC")

PROBT <- jags(data = model.data,
                           inits = inits,
                           parameters = params,
                           model.file = "PROBIT.txt",
                           n.thin = 1,
                           n.chains = 3,
                           n.burnin = 5000,
                           n.iter = 10000)

print(PROBT, intervals=c(0.025, 0.975), digits=3)

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

 

Output on screen:

Inference for Bugs model at "PROBIT.txt", fit using jags,

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

    n.sims = 15000 iterations saved

 

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

AIC              1635.902          2.452        1633.130       1642.264            1.009          310

BIC               1652.705         2.452        1649.932       1659.066            1.009          310

LogL              -814.951        1.226         -818.132        -813.565            1.009          310

beta[1]                 1.948        0.107              1.744             2.157             1.009         250

beta[2]                 0.710        0.128              0.464             0.971             1.003       1100

beta[3]                -1.253        0.069            -1.388            -1.118             1.016         150

deviance        1629.902        2.452        1627.130       1636.264             1.009         310

 

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 = 3.0 and DIC = 1632.9

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