rm(list = ls())
# enter hyperparameters of priors
# and data size
n=300
a=c(1.1,3,5.9,1.3,1.1,1.5)
b=c(8.5,11,1.7,12,4.9,5.5)
iterates=80000
burn_in=5000
library(R2OpenBUGS)
BUGS_path <-
"c:\\Program Files (x86)\\OpenBUGS\\OpenBUGS323\\OpenBUGS.exe"
setwd("c:\\E-drive\\Books\\LOGLIN3\\BAYES\\")
working_dir <- "c:\\E-drive\\Books\\LOGLIN3\\BAYES\\"

# Enter data
Trauma <- read.table(
url("http://stat.unm.edu/~fletcher/LLM/DATA/TRAUMAa.DAT"),
#"C:\\E-drive\\Books\\LOGLIN3\\DATA\\TRAUMAa.DAT",
		  sep="",col.names=c("ID","Death","ISS","TI","RTS","AGE"))
attach(Trauma)
Death=1-Death

#Enter Xtilde
Xtp=c(
 1, 25, 7.84, 60, 0,  0,
 1, 25, 3.34, 10, 0,  0,
 1, 41, 3.34, 60, 1, 60,
 1, 41, 7.84, 10, 1, 10,
 1, 33, 5.74, 35, 0,  0,
 1, 33, 5.74, 35, 1, 35)
Xtilde=t(matrix(Xtp,6,6))
Xtilde

# Find inverse of Xtilde
Xtinv=solve(Xtilde)
Xtinv

data <- list( "n", "Death", "a", "b" ,"Xtinv","ISS","RTS","AGE","TI")

inits <- function() {
  list(ptilde=c(0.5,0.5,0.5,0.5,0.5,0.5))
}

parameters <- list( "beta" )


Trauma.sim <- bugs( data, inits, parameters,
                   model.file="trauma_model.txt",
                   n.chains=1, n.iter=iterates+burn_in,
                   n.thin=1, n.burnin=burn_in,
                   OpenBUGS.pgm=BUGS_path,
                   working.directory=working_dir,
                   debug=F )
Trauma.sim$summary
# save p x t matrix of gamma posterior iterates
beta <- t(Trauma.sim$sims.list$beta)
save(beta,file="trauma_samp.Rda")
#save(beta,file="trauma-samp-thin.Rda")