## Hypothetical case-control study
N <- 100000
model <- function() {
  for (i in 1:2) {
    theta[i] ~ dunif(0,1)
    y[i] ~ dbin(theta[i], n[i])
  }
}
require(R2jags); invisible(runif(1))
RR <- OR <- matrix(NA, nrow=N, ncol=2)

## Analyzing as case-control
Data <- list(y=c(35, 15), n=c(75, n=75))
fit <- jags(Data, model=model, param="theta", n.iter=N, n.chain=1, n.thin=1, n.burn=0)
attach.jags(fit)
RR[,1] <- theta[,1]/theta[,2]
OR[,1] <- theta[,1]/(1-theta[,1])*(1-theta[,2])/theta[,2]
detach.jags()

## Analyzing as cohort (even though it isn't)
Data <- list(y=c(35, 40), n=c(50, n=100))
fit <- jags(Data, model=model, param="theta", n.iter=N, n.chain=1, n.thin=1, n.burn=0)
attach.jags(fit)
RR[,2] <- theta[,1]/theta[,2]
OR[,2] <- theta[,1]/(1-theta[,1])*(1-theta[,2])/theta[,2]
detach.jags()

pdf("or-vs-rr.pdf", height=4, width=8)
par(mfrow=c(1,2), mar=c(5,4,1,2), oma=c(0,0,3,0))
dnplot(RR, xlab="Relative risk")
dnplot(OR, xlab="Odds ratio")
toplegend(c("Case-control", "Cohort"), lwd=3, col=pal(2), ncol=2)
dev.off()