## Basic regression model
Data <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/701/data/alcohol.txt")
y <- Data$Metabol
X <- model.matrix(~Gastric+Sex+Alcohol, Data)
jData <- list(X=X, y=y, n=length(y), p=ncol(X))
model1 <- function() {
  ## Likelihood
  for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- X[i,1]*beta[1] + X[i,2]*beta[2] + X[i,3]*beta[3] + X[i,4]*beta[4]
  }  
  ## Prior
  for (j in 1:p) {
    beta[j] ~ dnorm(0, .0001)
  }
  tau <- pow(sigma, -2)
  log.sigma ~ dunif(-10, 10)
  sigma <- exp(log.sigma)
}

require(R2jags); invisible(runif(1))
fit <- jags(model=model1, data=jData, param=c("beta", "sigma"), n.iter=11000, n.chains=3, n.burn=1000, n.thin=1)
fit
dev <- dic.samples(fit$model, n.iter=5000, type="popt")
dev

require(R2OpenBUGS)
fit <- bugs(model=model1, inits=NULL, data=jData, param=c("beta", "sigma"), n.iter=11000, n.chains=1, n.burn=1000, n.thin=1)
fit

## Skeptical prior
model2 <- function() {
  ## Likelihood
  for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- X[i,1]*beta[1] + X[i,2]*beta[2] + X[i,3]*beta[3] + X[i,4]*beta[4]
  }  
  ## Prior
  beta[1] ~ dnorm(0, 0.0001)
  for (j in 2:p) {
    beta[j] ~ dnorm(0, l[j])
    l[j] <- 10 * pow(sd(X[,j]),2)
  }
  tau <- pow(sigma, -2)
  log.sigma ~ dunif(-10, 10)
  sigma <- exp(log.sigma)
}
fit <- bugs(model=model2, inits=NULL, data=jData, param=c("beta", "sigma"), n.iter=11000, n.chains=1, n.burn=1000, n.thin=1)
fit

fit <- jags(model=model2, data=jData, param=c("beta", "sigma"), n.iter=11000, n.chains=3, n.burn=1000, n.thin=1)
fit
dev <- dic.samples(fit$model, n.iter=1000, type="popt")
dev

## Bad quadratic approx.
xx <- seq(0,15,len=99)
require(numDeriv)
D <- function(x) -2*log(dt(10, df=4, ncp=x))
X <- -2*log(dt(10, df=4, ncp=xx))
plot(xx, D(xx), type="l", lwd=3, col=pal(2)[1], ylim=c(4,22), xlab=expression(theta), ylab="Deviance")
lines(xx, D(5) + grad(D, 5)*(xx-5) + hessian(D, 5)*(xx-5)^2, lwd=3, col=pal(2)[2])

## BUGS
model3 <- function() {
  y ~ dt(mu, 1, 4)
  mu ~ dt(0, 1, 4)
}
fit <- bugs(model=model3, inits=list(list(mu=5)), data=list(y=10), param=c("mu"), n.iter=10000, n.chains=1, n.thin=1)
fit
dnplot(fit$sims.list$mu)

fit <- jags(model=model3, data=list(y=10), param=c("mu"), n.iter=10000, n.chains=2, n.thin=1)
dic.samples(fit$model, 10000, type="popt")

## Interactions
X <- model.matrix(~Gastric*Sex*Alcohol, Data)
jData <- list(X=X, y=y, n=length(y), p=ncol(X))
model1 <- function() {
  ## Likelihood
  for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- X[i,1]*beta[1] + X[i,2]*beta[2] + X[i,3]*beta[3] + X[i,4]*beta[4] + X[i,5]*beta[5] + X[i,6]*beta[6] + X[i,7]*beta[7] + X[i,8]*beta[8]
  }  
  ## Prior
  for (j in 1:p) {
    beta[j] ~ dnorm(0, .0001)
  }
  tau <- pow(sigma, -2)
  log.sigma ~ dunif(-10, 10)
  sigma <- exp(log.sigma)
}
fit <- jags(model=model1, data=jData, param=c("beta", "sigma"), n.iter=11000, n.chains=3, n.burn=1000, n.thin=1)
dev <- dic.samples(fit$model, n.iter=5000, type="popt")
dev
fit <- bugs(model=model1, inits=NULL, data=jData, param=c("beta", "sigma"), n.iter=11000, n.chains=1, n.burn=1000, n.thin=1)
fit