## Alcohol metabolism example
## Load data
require(R2jags); invisible(runif(1))
Data <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/701/data/alcohol.txt")

## Uninformative prior
y <- Data$Metabol
X <- model.matrix(~Gastric+Sex+Alcohol, Data)
jData <- list(X=X, y=y, n=length(y), p=ncol(X))
model.ref <- function() {
  ## Likelihood
  for (i in 1:n) {
    y[i] ~ dnorm(X[i,]%*%beta, sigma^(-2))
  }
  
  ## Prior
  for (j in 1:p) {
    beta[j] ~ dnorm(0, .0001)
  }
  log.sigma ~ dunif(-10, 10)
  sigma <- exp(log.sigma)
}
fit <- jags(model=model.ref, data=jData, param=c("beta", "sigma"), n.iter=10000, n.chains=1, n.burn=0, n.thin=1, DIC=FALSE)
fit
xtable(fit$BUGSoutput$summary, digits=2)
CIplot(fit$BUGSoutput$summary[2:4, c(1, 3, 7)], labels=c("Gastric", "Sex", "Alcohol"), xlab=expression(beta))
abline(v=0, col="gray")
lm(y~Gastric+Sex+Alcohol, Data) ## Compare w/ OLS

## Skeptical prior
model.ridge <- function() {
  ## Likelihood
  for (i in 1:n) {
    y[i] ~ dnorm(X[i,]%*%beta, sigma^(-2))
  }
  
  ## Prior
  beta[1] ~ dnorm(0, 0.0001)
  for (j in 2:p) {
    beta[j] ~ dnorm(0, sd(X[,j])^2)
  }
  log.sigma ~ dunif(-10, 10)
  sigma <- exp(log.sigma)
}
fit <- jags(model=model.ridge, data=jData, param=c("beta", "sigma"), n.iter=10000, n.chains=1, n.burn=0, n.thin=1, DIC=FALSE)
fit
CIplot(fit$BUGSoutput$summary[2:4, c(1, 3, 7)], labels=c("Gastric", "Sex", "Alcohol"), xlab=expression(beta), xlim=c(-2, 3))
abline(v=0, col="gray")
require(MASS)
lm.ridge(y~Gastric+Sex+Alcohol, Data, lambda=2)

## An unidentifiable model
X <- cbind(model.matrix(~Gastric, Data), model.matrix(~0+Sex, Data))
jData <- list(X=X, y=y, n=length(y), p=ncol(X))
fit <- jags(model=model.ref, data=jData, param=c("beta", "sigma"), n.iter=10000, n.chains=1, n.burn=0, n.thin=1, DIC=FALSE)
fit
CIplot(fit$BUGSoutput$summary[1:4, c(1, 3, 7)], labels=c("Intercept", "Gastric", "Female", "Male"), xlab=expression(beta))
abline(v=0, col="gray")

## Unidentifiable model, skeptical prior
fit <- jags(model=model.ridge, data=jData, param=c("beta", "sigma"), n.iter=10000, n.chains=1, n.burn=0, n.thin=1, DIC=FALSE)
fit
CIplot(fit$BUGSoutput$summary[1:4, c(1, 3, 7)], labels=c("Intercept", "Gastric", "Female", "Male"), xlab=expression(beta))
abline(v=0, col="gray")

## Reparameterization
jData$X[,"Gastric"] <- jData$X[,"Gastric"]-mean(jData$X[,"Gastric"])
model.repar <- function() {
  ## Likelihood
  for (i in 1:n) {
    y[i] ~ dnorm(X[i,]%*%beta, sigma^(-2))
  }
  
  ## Prior
  for (j in 1:2) {
    beta[j] ~ dnorm(0, .0001)
  }
  delta ~ dnorm(0, .0001)
  beta[4] <- beta[1] + delta/2
  beta[3] <- beta[1] - delta/2
  log.sigma ~ dunif(-10, 10)
  sigma <- exp(log.sigma)
}
fit <- jags(model=model.repar, data=jData, param=c("beta", "sigma", "delta"), n.iter=10000, n.chains=1, n.burn=0, n.thin=1, DIC=FALSE)
fit
CIplot(fit$BUGSoutput$summary[c(1,3,4), c(1, 3, 7)], labels=c("Intercept", "Female", "Male"), xlab=expression(beta))
CIplot(fit$BUGSoutput$summary[c(2,5), c(1, 3, 7)], labels=c("Gastric", "Delta"), xlab=expression(beta), xlim=c(-0.5,3))
abline(v=0, col="gray")