rm(list = ls())
# Enter data
Ache <- read.table(
url("http://stat.unm.edu/~fletcher/LLM/DATA/TAB13-5.DAT"),
#"C:\\E-drive\\Books\\LOGLIN3\\DATA\\TAB13-5.dat",
		  sep="",header=T)
attach(Ache)
Ache

library(R2jags)

setwd("c:\\E-drive\\Books\\LOGLIN3\\BAYES\\")
working_dir <- "c:\\E-drive\\Books\\LOGLIN3\\BAYES\\"

# Larger burn_in  and iterates for
# highly correlated data
iterates <- 50000
burn_in <- 10000

inits <- function() {
  list( beta = c( 0, 0, 0),tau=1 )
}

ache.dat.jags<-list("kills","age","trials","pages")
ache.params<-c("beta","tau","ExKill","r")

ache.fit <- jags(data = ache.dat.jags, inits = inits,
  parameters.to.save = ache.params, n.chains = 1, n.thin=1,
n.iter = iterates+burn_in, n.burnin = burn_in,
model.file = "ache-m.txt")

# Rename Output
beta.a <- t(ache.fit$BUGSoutput$sims.list$beta)
tau.a = ache.fit$BUGSoutput$sims.list$tau
ExKill.a = ache.fit$BUGSoutput$sims.list$ExKill
r.a = ache.fit$BUGSoutput$sims.list$r

# Save Output
save(beta.a,file="ache-beta.Rda")
save(tau.a,file="ache-tau.Rda")
save(ExKill.a,file="ache-ExKill.Rda")
save(r.a,file="ache-lambda.Rda")

#Table of Coefficients for Model Parameters
# in book Table 13.6.
rowMeans(beta.a)
apply(beta.a, 1, mean) # Redundant command
apply(beta.a, 1, sd)
t(apply(beta.a, 1, quantile, probs=c(.5,0.025,0.975)))

sigma.a=1/sqrt(tau.a)
mean(sigma.a)
apply(sigma.a, 2, mean) # Redundant command
sd(sigma.a)
apply(sigma.a, 2, sd) # Redundant command
t(apply(sigma.a, 2, quantile, probs=c(.5,0.025,0.975)))

# Standard jags Summaries
summary(ache.fit$BUGSoutput$sims.list$beta)
summary(t(beta.a))
summary(ache.fit$BUGSoutput$sims.list$tau)
summary(ache.fit$BUGSoutput$sims.ist$ExKill)
summary(ache.fit$BUGSoutput$sims.list$r)

# Code for Figure 13.15
Aplot = t(apply(ExKill.a,2 , quantile, probs=c(.5,0.025,0.975)))
plot(pages,Aplot[,3],type="l",lty=1,lwd=2,xlab="Age",
ylab="Expected Kills",ylim=c(0,.8))
lines(pages,Aplot[,1],lty=2,lwd=2)
lines(pages,Aplot[,2],lty=1,lwd=2)

# Code for Figure 13.16
AAplot = t(apply(r.a,2 , quantile, probs=c(.5,0.025,0.975)))
plot(pages,AAplot[,3],type="l",lty=1,lwd=2,xlab="Age",
ylab="Latent Variable",ylim=c(0,1.3))
lines(pages,AAplot[,1],lty=2,lwd=2)
lines(pages,AAplot[,2],lty=1,lwd=2)