## Data; OLS
hills <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/701/data/hills.txt")
fit.ols <- lm(time~climb+dist, hills)
plot(fit.ols$fitted.values, fit.ols$residuals, pch=19, xlab="Fitted value", ylab="Residual")
plot(fit.ols$fitted.values, fit.ols$residuals, pch=19, xlab="Fitted value", ylab="Residual", cex=10*influence(fit.ols)$hat)

## Bayes setup
require(R2jags); invisible(runif(1))
pm <- pv <- matrix(NA, nrow=3, ncol=3) ## Posterior means, variances

## Normal errors
Data <- c(hills, n=nrow(hills))
model.z <- function() {
  for (i in 1:n) {
    time[i] ~ dnorm(beta[1] + beta[2]*climb[i] + beta[3]*dist[i], sigma^(-2))
  }
  for (j in 1:3) {
    beta[j] ~ dnorm(0, .0001)
  }
  log.sigma ~ dunif(-10,10)
  sigma <- exp(log.sigma)
}
fit <- jags(model=model.z, data=Data, param="beta", n.chains=1, n.iter=10000, n.burn=0, n.thin=1, DIC=FALSE)
attach.jags(fit)
pm[1,] <- apply(beta, 2, mean)
pv[1,] <- apply(beta, 2, var)
detach.jags()

## Thick-tailed errors
Data <- c(hills, n=nrow(hills))
model.t <- function() {
  for (i in 1:n) {
    time[i] ~ dt(beta[1] + beta[2]*climb[i] + beta[3]*dist[i], sigma^(-2), 5)
  }
  for (j in 1:3) {
    beta[j] ~ dnorm(0, .0001)
  }
  log.sigma ~ dunif(-10,10)
  sigma <- exp(log.sigma)
}
fit <- jags(model=model.t, data=Data, param="beta", n.chains=1, n.iter=10000, n.burn=0, n.thin=1, DIC=FALSE)
attach.jags(fit)
pm[2,] <- apply(beta, 2, mean)
pv[2,] <- apply(beta, 2, var)
detach.jags()

## Thick-tailed errors, prior on df
Data <- c(hills, n=nrow(hills))
model.t2 <- function() {
  for (i in 1:n) {
    time[i] ~ dt(beta[1] + beta[2]*climb[i] + beta[3]*dist[i], sigma^(-2), nu)
  }
  for (j in 1:3) {
    beta[j] ~ dnorm(0, .0001)
  }
  log.sigma ~ dunif(-10,10)
  sigma <- exp(log.sigma)
  nu ~ dgamma(.001,.001)
}
fit <- jags(model=model.t2, data=Data, param=c("beta", "nu"), n.chains=1, n.iter=10000, n.burn=0, n.thin=1, DIC=FALSE)
attach.jags(fit)
pm[3,] <- apply(beta, 2, mean)
pv[3,] <- apply(beta, 2, var)
dnplot(nu, xlab=expression(nu))
detach.jags()

## Summarizing the posteriors for beta
pm
pm[,2]*100
sqrt(pv)
100*sqrt(pv[,2])