# My own version of the example from http://mathstat.helsinki.fi/openbugs/ExamplesFrames.html
# Author: Mattias Villani, Statistics, Linkoping University, Sweden
#         mattias.villani@liu.se
# Date:   2012-12-05

library(BRugs)
setwd('~/Dropbox/Teaching/BayesLearn2012/Code/OpenBugs/')
nBurnin <- 1000
nIter <- 1000

# Data
seedsData <- list( r = c(10, 23, 23, 26, 17, 5, 53, 55, 32, 46, 10, 8, 10, 8, 23, 0, 3, 22, 15, 32, 3),
                   n = c(39, 62, 81, 51, 39, 6, 74, 72, 51, 79, 13, 16, 30, 28, 45, 4, 12, 41, 30, 51, 7),
                   x1 = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
                   x2 = c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1),
                   N = 21 )

# Model
LogisticRandEffectsModel <- function(){
  for (i in 1:N) {
    r[i] ~ dbin(p[i], n[i])
    b[i] ~ dnorm(0, tau)
    logit(p[i]) <- alpha0 + alpha1 * x1[i] + alpha2 * x2[i]
    + alpha12 * x1[i] * x2[i] + b[i]
  }
  alpha0 ~ dnorm(0, 1.0E-6)
  alpha1 ~ dnorm(0, 1.0E-6)
  alpha2 ~ dnorm(0, 1.0E-6)
  alpha12 ~ dnorm(0, 1.0E-6)
  tau ~ dgamma(0.001, 0.001)
  sigma <- 1 / sqrt(tau)
}

# Initial values
initVals <- list( list(alpha0 = 0, alpha1 = 0, alpha2 = 0, alpha12 = 0, tau = 10), 
                  list(alpha0 = 1, alpha1 = 2, alpha2 = 0, alpha12 = 0, tau = 10) )

# Run OpenBugs to sample from the posterior
BRugsFit(modelFile = LogisticRandEffectsModel, data = seedsData, inits = initVals, numChains = 2, 
         parametersToSave = c("alpha0", "alpha1", "alpha2", "alpha12", "tau", "sigma"),
         nBurnin, nIter, nThin = 1, coda = FALSE,
         DIC = TRUE, working.directory = NULL, digits = 5, seed=NULL,
         BRugsVerbose = getOption("BRugsVerbose"))

# Process the posterior draw. Summarize them, plot them etc
samplesHistory("*", mfrow = c(3, 2)) # Plots the MCMC trajectories (the draw plotted against the iteration number) for all parameters.
alpha0Post <- samplesSample("alpha0") # This is the vector of alpha0-draws from the posterior
hist(alpha0Post) # histogram of the draws to approximate the posterior of alpha0
samplesDensity("*") # Plots density estimates of the posterior draws for all the parameters.
plot(samplesSample("alpha1"),samplesSample("alpha2"),xlab="alpha1",ylab="alpha2", main = "Draws from joint posterior of alpha1 and alpha2")