#################
## GREAT study ##
#################

## Priors
logit <- binomial()$linkfun
inv.logit <- binomial()$linkinv
mu.delta <- c(0, log(2), log(1/.80))
sd.delta <- c(log(2)/1.96, log(2)/1.645, log(1/.80)/1.645)
mu.alpha <- c(0, 0, logit(.14)-0.5*log(1/.80))
sd.alpha <- c(1.6, 1.6, sqrt(((logit(.14)-logit(.11))/1.645)^2-sd.delta[3]^2/4))
main <- c("Skeptical", "Optimistic", "Clinical")
par(mfrow=c(3,4), mar=c(4,4,3,2))
for (i in 1:3) {
  delta <- rnorm(10000, mu.delta[i], sd.delta[i])
  dnplot(delta, xlab=expression(delta))
  mtext(main[i], line=0.5, cex=0.8)
  alpha <- rnorm(10000, mu.alpha[i], sd.alpha[i])
  dnplot(alpha, xlab=expression(alpha))
  mtext(main[i], line=0.5, cex=0.8)
  omega <- exp(delta)
  dnplot(omega, xlab=expression(omega))
  mtext(main[i], line=0.5, cex=0.8)
  L <- cbind(inv.logit(alpha+0.5*delta),
             inv.logit(alpha-0.5*delta))
  xlim <- if (i==3) c(0,0.3) else c(0,1)
  dnplot(L, xlab=expression(theta), xlim=xlim)
  mtext(main[i], line=0.5, cex=0.8)
}

## Data
y <- c(13, 23)
n <- c(163, 148)

## Reference
N <- 100000
OR <- RR <- matrix(NA, nrow=N, ncol=4)
colnames(OR) <- colnames(RR) <- c("Reference", main)
theta <- cbind(rbeta(N, y[1]+1, n[1]-y[1]+1),
               rbeta(N, y[2]+1, n[2]-y[2]+1))
RR[,1] <- theta[,2]/theta[,1]
OR[,1] <- theta[,2]/theta[,1]*(1-theta[1])/(1-theta[2])

## Model for informative priors
model <- function() {
  alpha ~ dnorm(mu.alpha, pow(sd.alpha, -2))
  delta ~ dnorm(mu.delta, pow(sd.delta, -2))
  theta[1] <- exp(alpha-0.5*delta)/(1+exp(alpha-0.5*delta))
  theta[2] <- exp(alpha+0.5*delta)/(1+exp(alpha+0.5*delta))
   for (i in 1:2) {
     y[i] ~ dbin(theta[i], n[i])
   }
}
require(R2jags)

## Sample from posteriors
rm(theta, delta, alpha)
for (j in 1:3) {
  Data <- list(y=y, n=n, mu.delta=mu.delta[j], sd.delta=sd.delta[j], mu.alpha=mu.alpha[j], sd.alpha=sd.alpha[j])
  fit <- jags(Data, model=model, param=c("theta", "delta", "alpha"), n.iter=N, n.chain=1, n.thin=1, n.burn=0)
  attach.jags(fit)
  RR[,j+1] <- theta[,2]/theta[,1]
  OR[,j+1] <- theta[,2]/theta[,1]*(1-theta[1])/(1-theta[2])
  detach.jags()
}

## Posterior
lOR <- log(OR)
par(mfrow=c(1,2), mar=c(5,4,1,2), oma=c(0,0,3,0))
dnplot(lOR, xlab=expression(delta))
dnplot(OR, xlab="Odds ratio")
toplegend(colnames(lOR), lwd=3, col=pal(4), ncol=4)

## Summary
med <- apply(lOR, 2, quantile, probs=.5)
CI <- apply(lOR, 2, quantile, probs=c(.025, .975))
B <- t(rbind(med, CI))
par(mfrow=c(1,2))
CIplot(B, xlab="Log odds")
abline(v=0, col="gray")
CIplot(B, xlab="Odds ratio", trans=exp)
abline(v=1, col="gray")

## Relative risk
par(mfrow=c(1,2), mar=c(5,4,1,2), oma=c(0,0,3,0))
dnplot(RR, xlab=expression(theta[2]/theta[1]))
dnplot(1/RR, xlab=expression(theta[1]/theta[2]))
toplegend(colnames(RR), lwd=3, col=pal(4), ncol=4)