Data <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/701/data/flight.txt")
source("http://web.as.uky.edu/statistics/users/pbreheny/701/S13/notes/fun.R")

## Descriptive
y <- with(Data, by(Recovered, list(Scenario, Group), mean))[,]
yy <- y[order(apply(y,1,mean)), order(apply(y,2,mean))]
col <- colorRampPalette(c("red","white"))(100)
levelplot(yy, xlab="", ylab="", col.regions=col) ## Plot this
levelplot(y, xlab="", ylab="", col.regions=col)  ## Not this

## Model
J <- nrow(y)
K <- ncol(y)
jData <- list("y", "J", "K")
model <- function() {
  ## Likelihood
  for (j in 1:J) {
    for (k in 1:K) {
      y[j,k] ~ dnorm(mu + a[j] + b[k], sigma[1]^(-2))
    }
  }
  
  ## Priors
  mu ~ dnorm(0, 0.0001)
  for (j in 1:J) {
    a[j] ~ dnorm(0, sigma[2]^(-2))
  }
  for (k in 1:K) {
    b[k] ~ dnorm(0, sigma[3]^(-2))
  }
  for (i in 1:3) {
    sigma[i] ~ dunif(0, 1)
  }
}
require(R2jags); invisible(runif(1))
fit <- jags(model=model, param=c("mu", "a", "b", "sigma"), data=jData, n.iter=20000, n.thin=1)
attach.jags(fit)

## R2jags output
plot(fit)
fit

## Convergence
max(gelman.diag(as.mcmc(fit))$psrf)
min(effectiveSize(as.mcmc(fit)))

################
## Posteriors ##
################

## Mu
psm(mu)
t.test(y) ## For comparison

## Sigma
psm(sigma, labels=c("y", "a", "b"))
prop.table(apply(sigma^2, 2, mean))

## Alpha
am <- apply(a, 2, mean)
ind <- order(am)
plot(1:J, type="n", ylim=c(-1,1), xaxt="n", xlab="Scenario", ylab=~alpha, las=1)
abline(h=0, col="gray80")
points(am[ind], pch=19)
axis(1, at=1:J, labels=abbreviate(rownames(y)))
for (j in 1:J) lines(c(j,j), quantile(a[,ind[j]], c(.025, .975)))

## Beta
bm <- apply(b, 2, mean)
ind <- order(bm)
plot(1:K, type="n", ylim=c(-1,1), xaxt="n", xlab="Group", ylab=~beta, las=1)
abline(h=0, col="gray80")
points(bm[ind], pch=19)
axis(1, at=1:K, labels=colnames(y))
for (k in 1:K) lines(c(k,k), quantile(b[,ind[k]], c(.025, .975)))

## Toledo mean
tm <- mu + a[,8]
psm(tm)