############## setwd or change dir  ###############################
setwd("..................................................")
library(BRugs)
modelCheck("gamma_poisson_BUGS_model.txt")
modelData("gamma_poisson_BUGS_data.txt")
modelCompile(numChains=2)
modelInits(rep("gamma_poisson_BUGS_inits.txt",2))
modelUpdate(10000)
samplesSet(c("theta","mu","a","b"))
dicSet()
modelUpdate(2000)
dicStats()
samplesStats("*")

theta<-samplesStats("theta")
###########################################################
############ log normal model ##################################
modelCheck("log_normal_BUGS_model.txt")
modelData("gamma_poisson_BUGS_data.txt")
modelCompile(numChains=2)
modelInits(rep("log_normal_BUGS_inits.txt",2))
modelUpdate(10000)
samplesSet(c("theta","a0","z1","deviance"))
dicSet()
modelUpdate(2000)
dicStats()
asd<-samplesStats("*")

library(coda)
gelman.plot(buildMCMC("deviance"))


######################################################################
############# Convolution model ######################################
modelCheck("CONV_BUGS_model.txt")
modelData("CONV_BUGS_data.txt")
modelCompile(numChains=2)
modelInits(rep("CONV_BUGS_inits.txt",2))
modelUpdate(30000)
samplesSet(c("theta","a0","v","u","deviance"))
dicSet()
modelUpdate(5000,thin=10)
dicStats()
asd<-samplesStats("*")
theta<-asd$mean[3:48]
u<-asd$mean[49:94]
v<-asd$mean[95:140]
library(coda)
gelman.plot(buildMCMC("deviance"))
######################  graphics  ###############################
source("fillmap.R")
library(maptools);SCpoly<-readSplus("SC_geobugsSPlus.txt")
fillmap(SCpoly,"theta",theta,n.col=5);x11()
fillmap(SCpoly,"correlated H",u,n.col=5);x11()
fillmap(SCpoly,"uncorrelated H",v,n.col=5)
###################################################################

##### spat MIx (SMix)   ##########################################
modelCheck("SMix_BUGS_model.txt")
modelData("SMix_BUGS_data.txt")
modelCompile(numChains=2)
modelInits(rep("SMix_BUGS_inits.txt",2))
modelUpdate(30000)
samplesSet(c("theta","a0","p","v","u","deviance"))
dicSet()
modelUpdate(5000,thin=10)
dicStats()
asd<-samplesStats("*")
library(coda)
gelman.plot(buildMCMC("deviance"))

########## Leroux model  ##########################################
modelCheck("Leroux_BUGS_model_Q1RD.txt")
modelData("Leroux_BUGS_data_Q1RDmodel.txt")
modelCompile(numChains=2)
modelInits(rep("Leroux_BUGS_inits_Q1RDmodel.txt",2))
modelUpdate(50000)
samplesSet(c("theta","a0","s","rho","deviance"))
dicSet()
modelUpdate(5000,thin=10)
dicStats()
asd<-samplesStats("*")
sLER<-samplesStats("s")
thetaLER<-samplesStats("theta")
library(coda)
gelman.plot(buildMCMC("deviance"))
LERM<-cbind(sLER,thetaLER)
#####################################################################
############ BYM2 model  ############################################
modelCheck("BYM2_BUGS_model.txt")
modelData("BYM2_BUGS_data.txt")
modelCompile(numChains=2)
modelInits(rep("BYM2_BUGS_inits.txt",2))
modelUpdate(50000)
samplesSet(c("theta","UCorr","Corr","rho","sigma","deviance"))
dicSet()
modelUpdate(5000,thin=10)
dicStats()
asd<-samplesStats("*")
Corr<-asd$mean[1:46]
UCorr<-asd$mean[47:92]
rho<-asd$mean[93]
sigma<-asd$mean[94]
theta<-asd$mean[95:140]

gelman.plot(buildMCMC("deviance"))

BYM2M<-cbind(Corr,UCorr,theta)
library(sf)
library(sp)
SCpoly<-st_read("SC_county_alphasort.shp")
SCpoly<-st_geometry(SCpoly)
library(fillmap)
fillmaps(SCpoly,"",BYM2M,n.col=4,leg.cex=0.9,
lay.m=matrix(c(1,2,3),ncol=3,nrow=1))