## Read in data
lead <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/701/data/lead.txt")

## A frequentist analysis
x <- subset(lead, Smelter=="Near")$IQ
y <- subset(lead, Smelter=="Far")$IQ
t.test(x, y, var.equal=TRUE)

## Bayesian model
require(R2jags); invisible(runif(1))
Data <- list(x=x, y=y, nx=length(x), ny=length(y))
model <- function() {
  ## Likelihood
  for (i in 1:nx) {
    x[i] ~ dnorm(mu[1], tau)
  }
  for (i in 1:ny) {
    y[i] ~ dnorm(mu[2], tau)
  }
  
  ## Prior
  mu[1] ~ dunif(0, 200)
  mu[2] ~ dunif(0, 200)
  tau <- pow(sigma, -2)
  sigma ~ dunif(0, 100)
}
fit <- jags(Data, param=c("mu","sigma"), model=model, n.chains=1, n.iter=10000, n.burn=0, n.thin=1, DIC=FALSE)

## Examining the posterior
source("http://web.as.uky.edu/statistics/users/pbreheny/701/S13/notes/fun.R")
attach.jags(fit)
dim(mu)
head(mu)
dnplot(mu, c("Near", "Far"), xlab="IQ")
dnplot(sigma, xlab=expression(sigma))
dnplot(mu[,2]-mu[,1], xlab=expression(Delta*"IQ"))
dnplot((mu[,2]-mu[,1])/sigma, xlab=expression(Delta*"IQ (Std. Units)"))

rbind(Mean=apply(mu, 2, mean), SD=apply(mu, 2, sd), apply(mu, 2, quantile, p=c(.025,.975)))
rbind(Mean=apply(sigma, 2, mean), SD=apply(sigma, 2, sd), apply(sigma, 2, quantile, p=c(.025,.975)))
diff <- mu[,2]-mu[,1]
c(Mean=mean(diff), SD=sd(diff), quantile(diff, c(.025, .975)))
mean(diff<0)
SNR <- (mu[,2]-mu[,1])/sigma
rbind(Mean=apply(SNR, 2, mean), SD=apply(SNR, 2, sd), apply(SNR, 2, quantile, p=c(.025,.975)))