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")

## Setup
aid <- as.numeric(Data$Scenario)
uid <- as.numeric(Data$Group)
J <- length(levels(Data$Scenario))
K <- length(levels(Data$Group))
y <- Data$Recovered
n <- nrow(Data)

## Model
jData <- list("y", "J", "K", "n", "aid", "uid")
model <- function() {
  ## Likelihood
  for (i in 1:n) {
    y[i] ~ dbern(theta[i])
    logit(theta[i]) <- mu + a[aid[i]] + b[uid[i]]
  }
  
  ## Priors
  mu ~ dnorm(0, 0.0001)
  for (j in 1:J) {
    a[j] ~ dnorm(0, sigma[1]^(-2))
  }
  for (k in 1:K) {
    b[k] ~ dnorm(0, sigma[2]^(-2))
  }
  for (i in 1:2) {
    sigma[i] ~ dunif(0, 5)
  }
}
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
psm(sigma, labels=c("a", "b"))
psm(mu)
ilogit(psm(mu))

## Alpha
am <- apply(a, 2, mean)
ind <- order(am)
plot(1:J, type="n", ylim=c(-5,6), 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(levels(Data$Scenario)))
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(-5,6), xaxt="n", xlab="Group", ylab=~beta, las=1)
abline(h=0, col="gray80")
points(bm[ind], pch=19)
axis(1, at=1:K, labels=levels(Data$Group))
for (i in 1:K) lines(c(i,i), quantile(b[,ind[i]], c(.025, .975)))

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