#### SC congenital data 1990  CONvolution model #######
library(nimble)
LconvCode<-nimbleCode({
for(i in 1:N){
log(theta[i])<-a0+v[i]+B[i]
lambda[i]<-theta[i]*texp[i]
y[i]~dpois(lambda[i])
v[i]~dnorm(0,tauv)
}
for(k in 1:L){wei[k]<-1}
B[1:N]~dcar_normal(adj[1:L],wei[1:L],num[1:N],tauB,zero_mean=1)
a0~dnorm(0,tau0)
tauv~dgamma(2,0.5)
tau0~dgamma(2,0.5)
tauB~dgamma(2,0.5)
})
########################################################################

######################## SC data  #######################################

LconvConsts<-list(N=46,L=228,texp=c(1.0807,6.3775,0.622,6.6854,0.9142,1.0744,5.6518,
8.1682,0.5749,18.0989,2.174,1.6619,1.9321,1.6148,1.6713,3.0819,1.7562,4.9952,
0.9362,1.2001,6.1293,2.5604,15.8589,2.9437,1.0399,7.276,0.9739,2.064,2.7206,2.8275,
0.9425,8.828,0.3644,1.775,1.5111,1.5111,2.5321,
4.5836,3.9647,15.0264,0.732,10.8292,5.9848,1.4357,1.9949,6.9807),
pov=c(13.6,13.8,32.3,11,24.2,19.9,12.3,13.3,17,15.4,14.1,15.7,18,24.3,21.8,
20.2,24.9,12,17.4,18.1,18.7,17.5,10.6,13.8,22.8,13.7,21,12.6,14,14,26.9,
9.4,17.8,23.1,23.1,14.4,10.8,22.1,10.1,13.6,15.9,11.2,18.3,14,26.4,10.6),
inc=c(36.786,38.534,20.727,37.205,24.3,27.607,45.822,40.161,32.247,38.458,
33.232,31.715,29.505,25.896,28.919,30.776,25.552,42.886,34.297,29.96,34.009,
35.008,41.658,34.109,27.65,34.654,27.117,39.04,33.698,32.32,25.144,45.14,
29.805,25.008,25.993,32.231,36.912,28.624,37.054,39.587,31.324,37.092,
31.948,30.801,23.748,44.619),
num = c(5, 5, 4, 5, 5, 4, 3, 6, 5, 4, 
3, 4, 5, 6, 7, 5, 4, 4, 4, 6, 
8, 5, 5, 6, 5, 3, 2, 7, 5, 7, 
5, 6, 3, 5, 4, 7, 2, 9, 3, 6, 
5, 4, 6, 7, 5, 4
),
adj = c(
33, 30, 24, 23, 4, 
41, 38, 32, 19, 6, 
25, 15, 6, 5, 
39, 37, 30, 23, 1, 
38, 25, 15, 6, 3, 
38, 5, 3, 2, 
27, 25, 15, 
45, 38, 22, 18, 14, 10, 
43, 40, 38, 32, 14, 
22, 18, 15, 8, 
46, 44, 42, 
46, 44, 29, 20, 
35, 31, 29, 28, 16, 
45, 43, 38, 21, 9, 8, 
38, 25, 18, 10, 7, 5, 3, 
35, 31, 28, 21, 13, 
35, 34, 26, 21, 
38, 15, 10, 8, 
41, 33, 24, 2, 
44, 40, 36, 29, 28, 12, 
45, 43, 35, 34, 31, 17, 16, 14, 
45, 34, 26, 10, 8, 
42, 39, 30, 4, 1, 
41, 36, 33, 30, 19, 1, 
27, 15, 7, 5, 3, 
34, 22, 17, 
25, 7, 
43, 40, 31, 29, 20, 16, 13, 
46, 28, 20, 13, 12, 
44, 42, 36, 24, 23, 4, 1, 
43, 28, 21, 16, 13, 
41, 40, 38, 36, 9, 2, 
24, 19, 1, 
45, 26, 22, 21, 17, 
21, 17, 16, 13, 
44, 41, 40, 32, 30, 24, 20, 
39, 4, 
32, 18, 15, 14, 9, 8, 6, 5, 2, 
37, 23, 4, 
43, 36, 32, 28, 20, 9, 
36, 32, 24, 19, 2, 
44, 30, 23, 11, 
40, 31, 28, 21, 14, 9, 
46, 42, 36, 30, 20, 12, 11, 
34, 22, 21, 14, 8, 
44, 29, 12, 11
),sumNumNeigh = 228)
LconvData<-list(y=c(0,7,1,5,1,1,5,16,0,17,4,0,0,1,1,7,1,3,0,0,8,2,13,7,0,8,0,3,2,4,1,11,0,1,2,3,3,8,6,14,3,11,6,0,1,5))
LconvInits<-list(a0=0.1,tauv=0.1,tau0=0.1,tauB=0.1,v=rep(0.0,LconvConsts$N),B=rep(0.0,LconvConsts$N))

####################################################################
##### quick run ####################################################
mcmc.out<-nimbleMCMC(code=LconvCode,constants=LconvConsts,
data=LconvData,inits=LconvInits,nchains=2,nburnin=7000,niter=10000,summary=TRUE,WAIC=TRUE,monitors=c("a0","theta","B"))

mcmc.out$summary
mcmc.out$WAIC

######################################################################
####### separate steps ###################################################

Lconv<-nimbleModel(code=LconvCode,name="Lconv",constants=LconvConsts,
data=LconvData,inits=LconvInits)

######################## Compile  ######################################
CLconv<-compileNimble(Lconv)

LconvConf<-configureMCMC(Lconv,print=TRUE)
LconvConf$addMonitors(c("a0","tauv","tau0","v","B","theta"))


LconvMCMC<-buildMCMC(LconvConf)
CLconvMCMC<-compileNimble(LconvMCMC,project=Lconv)
###########################################################################
niter<-10000
set.seed(1)
samples<-runMCMC(CLconvMCMC,niter=niter,nburnin=7000,nchains=2,summary=TRUE,samplesAsCodaMCMC=TRUE)  
# if specifying CODA then use library(coda); then gelman.diag or geweke.diag and plot
# the samples$samples is an mcmc.list   eg samples$samples$chain1 lists the first chain etc
# to access elements of the output eg samples$samples$chain1[,"B[1]"] will yield the sampler values for B[1]
as a vector. Diagnostics: geweke.diag(samples$samples$chain1[,"B[1]"]), plotting: plot(samples$samples$chain1[,"B[1]"])

samples$summary  # this a list as two chains run and its CODA?

samples$summary$all.chains  # this is a matrix of 142 by 5....first column is th mean
Bmean<-samples$summary$all.chains[1:46,1]   #gives B values
a0mean<-samples$summary$all.chains[47,1]
tau0mean<-samples$summary$all.chains[48,1]
tauBmean<-samples$summary$all.chains[49,1]
tauvmean<-samples$summary$all.chains[50,1]
THmean<-samples$summary$all.chains[51:96,1]
vmean<-samples$summary$all.chains[97:142,1]

#### now a single chain run without CODA  ########################################

samples2<-runMCMC(CLconvMCMC,niter=niter,nburnin=7000,summary=TRUE)
samples2$summary   # this is a matrix as only 1 chian run

#####    posterior means  ####################################################
Bmean<-samples2$summary[1:46,1]  # mean of B values
a0mean<-samples2$summary[47,1]
tau0mean<-samples2$summary[48,1]
tauBmean<-samples2$summary[49,1]
tauvmean<-samples2$summary[50,1]
THmean<-samples2$summary[51:96,1]
vmean<-samples2$summary[97:142,1]
WAIC<-samples2$WAIC

sampMAT<-as.matrix(samples2$samples)  ### dim is 3000 by 142
par(mfrow=c(2,2))
plot(sampMAT[,"a0"],type="l",xlab="iteration",ylab=expression(a0))
plot(sampMAT[,"tauv"],type="l",xlab="iteration",ylab=expression(tauv))
plot(sampMAT[,"tau0"],type="l",xlab="iteration",ylab=expression(tau0))
plot(sampMAT[,"tauB"],type="l",xlab="iteration",ylab=expression(tauB))
par(mfrow=c(2,2))
plot(density(sampMAT[,"a0"]),xlab=expression(a0))
plot(density(sampMAT[,"tauv"]),xlab=expression(tauv))
plot(density(sampMAT[,"tau0"]),xlab=expression(tau0))
plot(density(sampMAT[,"tauB"]),xlab=expression(tauB))

setwd("C:/A_NEW_FOLDER/A_NEW_FOLDER/A_PROGRAMS/R_program/NIMBLE examples")
library(maptools)
library(sp)
source("fillmap.R")
SCpoly<-readSplus("SC_SPLus_export_map.txt")
plot(SCpoly)

fillmap(SCpoly,"UH",vmean,n.col=5);x11()
fillmap(SCpoly,"Theta",THmean,n.col=5);x11()
fillmap(SCpoly,"CH",Bmean,n.col=5)