#setwd("c:/Users/gustaf/Documents/Workspace")
library("arm")    ### includes the R -> BUGS stuff

### useful functions
logit <- function(p) { log(p)-log(1-p) }
expit <- function(x) {1/(1+exp(-x)) }
trpl <- function(x) {c(quantile(x, .025), mean(x), quantile(x, .975))}

### choice of prior on exposure prevlances
tmp.sd <- logit(.95)/1.96
tmp.rho <- 1 - (log(6)/1.95)^2 / (2*tmp.sd^2)

### three pre-cursors to calling BUGS
### first, data to pass from R to BUGS
dattobugs <- list(
  wtot1=101, n1=580, wtot2=122, n2=564,
  xwtot1=c(168, 12, 16, 21), m1=217,
  xwtot2=c(143, 22, 17, 29), m2=211,
  mur=0, varr=tmp.sd^2, rhor=tmp.rho,
  muSN=1.675, varSN=.648^2, rhoSN=0.95,
  muSP=1.675, varSP=.648^2, rhoSP=0.95)

### second, scheme to generate initial values for MCMC
  inittobugs <- function() {
    list(lgtSP=c(2,2), lgtSN=c(2.1,2.1), lgtr=c(-0.25,-0.5))
  }

### third, parameters to monitor
  paramstobugs <- c("SN", "SP", "r")

### now call BUGS, model in file "val.txt"
bugsfit1 <- bugs(dattobugs, inittobugs, paramstobugs,
            "val.txt",
            n.chains=1, n.thin=10,
            n.iter=100000, bugs.directory="c:/WinBUGS14/",
            debug=F)

### can add debug=T above to see more output as it goes

### make output `live'
attach.bugs(bugsfit1)

### summary of inference on OR
print(trpl(exp(logit(r[,2])-logit(r[,1]))))