## ANALYSIS OF THE NORMAL MODEL - pag 129-141

# define the working directory assuming - the user is already in ~/WebMaterial/ch2
working.dir <- paste(getwd(),"/",sep="")

# Load data from .dta format
library(foreign)
data <- read.dta("../ch2/phbirths.dta")
attach(data)

# Arrange the data and define relevant variables
y <- grams
N <- length(y)
X <- gestate
k <- 5000

# Run JAGS
library(R2jags)

## 4. Uncentred model with blocking (for chapter 4)
## NB Continues the analysis from chapter 2 (check modelNormal.R in the directory /ch2)
X <- gestate
m0 <- c(0,0)
prec <- 1/100000*matrix(c(1,0,0,1),nrow=2,ncol=2)
dataJags <- list("m0","prec","N","y","X","k") 
filein <- "modelNormalBlocking.txt"
params <- c("alpha","beta","sigma")
inits <- function(){
	list(lsigma=rnorm(1),coef=rnorm(2,0,1))
}

n.iter <- 10000
n.burnin <- 9500
n.thin <- floor((n.iter-n.burnin)/500)
m.4 <- jags(dataJags, inits, params, model.file=filein,
	n.chains=2, n.iter, n.burnin, n.thin,
	DIC=TRUE, working.directory=working.dir, progress.bar="text")
print(m.4,digits=3,intervals=c(0.025, 0.975))



## 5. Predictive distribution (for chapter 4)
# 5a. Uncentered version with thinning
X.star <- 28
dataJags2 <- list("N","y","X","k","X.star") 
filein <- paste(working.dir,"modelNormal2.txt",sep="")
params <- c("alpha","beta","sigma","y.star")
inits <- function(){
	list(alpha=rnorm(1),beta=rnorm(1),lsigma=rnorm(1),y.star=runif(1))
}

n.iter <- 50000
n.burnin <- 9500
n.thin <- floor((n.iter-n.burnin)/500)
m.5 <- jags(dataJags2, inits, params, model.file=filein,
	n.chains=2, n.iter, n.burnin, n.thin,
	DIC=TRUE, working.directory=working.dir, progress.bar="text")
print(m.5,digits=3,intervals=c(0.025, 0.975))

# 5b. Centred model with missing data
y <- c(y,NA)
N <- length(y)
X <- c(X,28)
# still call centred variable X so can still use same model file
#Xbar <- mean(X)
#X <- X - Xbar
dataJags3 <- list("N","y","X","k") 
filein <- "../ch2/modelNormal.txt"
params3 <- c("alpha","beta","sigma","y[1116]")
inits <- function(){
list(alpha=rnorm(1),beta=rnorm(1),lsigma=rnorm(1),
  y=c(rep(NA,(N-1)),rnorm(1)))
}
n.iter <- 50000
n.burnin <- 9500
n.thin <- floor((n.iter-n.burnin)/500)

m.6 <- jags(dataJags3, inits, params3, model.file=filein,
	n.chains=2,n.iter=n.iter,n.burnin=n.burnin,n.thin=n.thin,DIC=TRUE)
print(m.6,digits=3,intervals=c(0.025, 0.975))