diff --git a/BYM2_SC_scaling.txt b/BYM2_SC_scaling.txt new file mode 100644 index 0000000..508746e --- /dev/null +++ b/BYM2_SC_scaling.txt @@ -0,0 +1,162 @@ + +###SCcongen90 is dataframe including the vectors obs,expe,pov,inc###### +####################################################################### +SCcongen90<-list(obs=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), +expe=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)) +######################################################################## +SCcon<-data.frame(SCcongen90) +#graphics....using SC_poly +library(INLA) +setwd("C:/A_NEW_FOLDER/A_NEW_FOLDER/A_PROGRAMS/GITHUB repository") +#SC_poly.txt is the graph file +g=inla.read.graph("SC_poly.txt") +Q=-inla.graph2matrix(g) +diag(Q) = 0 +diag(Q) = -rowSums(Q) +n=dim(Q)[1] # n=46 +Q.scaled=inla.scale.model(Q,constr=list(A=matrix(1,1,n),e=0)) +scale=Q.scaled[1,1]/Q[1,1] +# scale of SC graph is 0.3783966 +### BYM2 model code in Nimble +library(nimble) +BYM2<-nimbleCode({ +for(i in 1: n){ +obs[i]~dpois(mu[i]) +log(mu[i])<-log(expe[i])+log(theta[i]) +theta[i]<-exp(ac0+bym2[i]) +vc[i]~dnorm(0,1) +Corr[i]<-sigma*uc[i]*sqrt(rho/scale) +UCorr[i]<-sigma*vc[i]*sqrt((1-rho)) +bym2[i]<-Corr[i]+UCorr[i] +} +for (k in 1 :nsumN){wei[k]<-1} +uc[1:n]~dcar_normal(adj[1:nsumN],wei[1:nsumN],num[1:n],1,zero_mean=1) +sigma~dunif(0,10) +tau<-pow(sigma,-2) +rho~dbeta(1,1) +ac0~dnorm(0,tau0) +tau0<-pow(sd0,-2) +sd0~dunif(0,10) +}) +nsumN=228 +LSTdata<-list(obs=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)) +LSTconsts<-list(n,scale,nsumN,expe=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)) +LSTinits<-list(vc=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, +0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), +uc=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, +0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, +0, 0, 0, 0, 0),sigma=0.1,rho=0.5,ac0=0.1,sd0=0.1) +names(LSTconsts)<-c("n","scale","nsumN","expe","pov","inc","num","adj") +LBYM2<-nimbleModel(code=BYM2,name="BYM2",constants=LSTconsts, +data=LSTdata,inits=LSTinits) +################################################################### +CLBYM2<-compileNimble(LBYM2) +##################################################################### +LBYM2Conf<-configureMCMC(LBYM2,print=TRUE) +LBYM2Conf$addMonitors(c("ac0","sigma","rho","vc","uc","Corr","UCorr","theta")) +LBYM2MCMC<-buildMCMC(LBYM2Conf) +CLBYM2MCMC<-compileNimble(LBYM2MCMC,project=LBYM2) +niter<-10000 +set.seed(1) +samples<-runMCMC(CLBYM2MCMC,niter=niter,nburnin=7000,nchains=1,summary=TRUE,samplesAsCodaMCMC=TRUE) +samples$summary +names(samples$summary) +MCorr<-samples$summary[1:46,1] +MUCorr<-samples$summary[47:92,1] +Mac0<-samples$summary[93,1] +Mrho<-samples$summary[94,1] +Msd<-samples$summary[95,1] +Msigma<-samples$summary[96,1] +Mtheta<-samples$summary[97:142,1] +Muc<-samples$summary[143:188,1] +Mvc<-samples$summary[189:234,1] +hist(Muc);x11();hist(Mvc);x11();hist(Mtheta);x11() + + +############################################### +SCcon<-data.frame(MCorr,MUCorr,Muc,Mvc,Mtheta) +#SCcon is dataframe with vectors eg SCcon<-data.frame(......Muc,Mvc,Mtheta) +library(tmap) +areaID<-as.character(seq(1:46)) +area<-paste(c("area"),areaID,sep="") +attr<-data.frame(SCcon,row.names=area) +library(maptools) +polySC<-readSplus("SC_SPLus_export_map.txt") +spg<-SpatialPolygonsDataFrame(polySC, attr, match.ID = TRUE) +#### this qtm call produces common legend ############################################# +labels<-c("Uc","Vc","Theta") +qtm(spg,fill=c("Muc","Mvc","Mtheta"),fill.palette="Blues",ncol=3,free.scales=FALSE) +### use tm_layout with legend.outside=TRUE and separate titles ########################### +qtm(spg,fill=c("Muc","Mvc","Mtheta"),fill.palette="Blues",ncol=3)+tm_layout(legend.outside=TRUE,title="",panel.labels=labels) +### separate legends +qtm(spg,fill=c("Muc","Mvc","Mtheta"),fill.palette="Blues",ncol=3)+tm_layout(legend.position=c("LEFT","BOTTOM"),main.title="BYM2 results",panel.labels=labels) +############################################################# +####### mapping the UCorr and Corr components +labels<-c("Ucorr","Corr","Theta") +qtm(spg,fill=c("MUCorr","MCorr","Mtheta"),fill.palette="Blues",ncol=3)+tm_layout(legend.position=c("LEFT","BOTTOM"),main.title="BYM2 UCorr, Corr",panel.labels=labels) +###################################################### +mcmc.out<-nimbleMCMC(code=BYM2,constants=LSTconsts, +data=LSTdata,inits=LSTinits,nchains=2,nburnin=7000,niter=10000,summary=TRUE,WAIC=TRUE) \ No newline at end of file diff --git a/BYM_ordinary_SC_NOscaling.txt b/BYM_ordinary_SC_NOscaling.txt new file mode 100644 index 0000000..1676b12 --- /dev/null +++ b/BYM_ordinary_SC_NOscaling.txt @@ -0,0 +1,149 @@ + +###SCcongen90 is dataframe including the vectors obs,expe,pov,inc###### +####################################################################### +SCcongen90<-list(obs=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), +expe=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)) +######################################################################## +SCcon<-data.frame(SCcongen90) +#graphics....using SC_poly +library(INLA) +setwd("C:/A_NEW_FOLDER/A_NEW_FOLDER/A_PROGRAMS/GITHUB repository") +#SC_poly.txt is the graph file + + +### BYM model code in Nimble +library(nimble) +BYM<-nimbleCode({ +for(i in 1: n){ +obs[i]~dpois(mu[i]) +log(mu[i])<-log(expe[i])+log(theta[i]) +theta[i]<-exp(ac0+bym[i]) +vc[i]~dnorm(0,tauV) +bym[i]<-vc[i]+uc[i] +} +for (k in 1 :nsumN){wei[k]<-1} +uc[1:n]~dcar_normal(adj[1:nsumN],wei[1:nsumN],num[1:n],tauU,zero_mean=1) +tauU<-pow(sdU,-2) +tauV<-pow(sdV,-2) +sdV~dunif(0,10) +sdU~dunif(0,10) +ac0~dnorm(0,tau0) +sd0~dunif(0,10) +}) +nsumN=228 +LSTdata<-list(obs=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)) +LSTconsts<-list(n,nsumN,expe=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)) +LSTinits<-list(vc=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, +0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), +uc=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, +0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, +0, 0, 0, 0, 0),ac0=0.1,sd0=0.1,sdV=0.1,sdU=0.1) +names(LSTconsts)<-c("n","nsumN","expe","pov","inc","num","adj") +LBYM<-nimbleModel(code=BYM,name="BYM",constants=LSTconsts, +data=LSTdata,inits=LSTinits) +################################################################### +CLBYM<-compileNimble(LBYM) +##################################################################### +LBYMConf<-configureMCMC(LBYM,print=TRUE) +LBYMConf$addMonitors(c("ac0","vc","uc","theta")) +LBYMMCMC<-buildMCMC(LBYMConf) +CLBYMMCMC<-compileNimble(LBYMMCMC,project=LBYM) +niter<-10000 +set.seed(1) +samples<-runMCMC(CLBYMMCMC,niter=niter,nburnin=7000,nchains=1,summary=TRUE,samplesAsCodaMCMC=TRUE) +samples$summary + +Mac0<-samples$summary[1,1] +Msd0<-samples$summary[2,1] +MsdU<-samples$summary[3,1] +MsdV<-samples$summary[4,1] +Mtheta<-samples$summary[5:50,1] +Muc<-samples$summary[51:96,1] +Mvc<-samples$summary[97:142,1] +hist(Muc);x11();hist(Mvc);x11();hist(Mtheta);x11() + + +############################################### +SCcon<-data.frame(Muc,Mvc,Mtheta) +#SCcon is dataframe with vectors eg SCcon<-data.frame(......Muc,Mvc,Mtheta) +library(tmap) +areaID<-as.character(seq(1:46)) +area<-paste(c("area"),areaID,sep="") +attr<-data.frame(SCcon,row.names=area) +library(maptools) +polySC<-readSplus("SC_SPLus_export_map.txt") +spg<-SpatialPolygonsDataFrame(polySC, attr, match.ID = TRUE) +#### this qtm call produces common legend ############################################# +labels<-c("Uc","Vc","Theta") +qtm(spg,fill=c("Muc","Mvc","Mtheta"),fill.palette="Blues",ncol=3,free.scales=FALSE) +### use tm_layout with legend.outside=TRUE and separate titles ########################### +qtm(spg,fill=c("Muc","Mvc","Mtheta"),fill.palette="Blues",ncol=3)+tm_layout(legend.outside=TRUE,title="",panel.labels=labels) +### separate legends ############################################## +qtm(spg,fill=c("Muc","Mvc","Mtheta"),fill.palette="Blues",ncol=3)+tm_layout(legend.position=c("LEFT","BOTTOM"),main.title="BYM results",panel.labels=labels) +############################################################# + +mcmc.out<-nimbleMCMC(code=BYM,constants=LSTconsts, +data=LSTdata,inits=LSTinits,nchains=2,nburnin=7000,niter=10000,summary=TRUE,WAIC=TRUE) \ No newline at end of file diff --git a/Georgia_MultiScale_CONV_Nimble.txt b/Georgia_MultiScale_CONV_Nimble.txt new file mode 100644 index 0000000..c204c88 --- /dev/null +++ b/Georgia_MultiScale_CONV_Nimble.txt @@ -0,0 +1,314 @@ +#### Goergia county/Ph district joint multi-scale model ######## +MultiSC<-nimbleCode({ + +for( i in 1:18){ +yph[i]~dpois(muph[i]) +log(muph[i])<-log(eph[i])+log(thph[i]) +thph[i]<-exp(aph0+uph[i]+vph[i]) +vph[i]~dnorm(0,0.0001) +} +for( j in 1:159){ +yc[j]~dpois(muc[j]) +log(muc[j])<-log(ec[j])+log(thc[j]) +thc[j]<-exp(ac0+uc[j]+vc[j]) +vc[j]~dnorm(0,0.0001) +} +for (k in 1 :nsumph){weiph[k]<-1} +uph[1:18]~dcar_normal(adj[1:nsumph],weiph[1:nsumph],num[1:18],tauph) +tauph<-1/pow(sdph,2) +aph0~dflat() +sdph~dunif(0,100) +for (o in 1 :nsumc){weic[o]<-1} +uc[1:159]~dcar_normal(adj2[1:nsumc],weic[1:nsumc],num2[1:159],tauc) +tauc<-1/pow(sdc,2) +ac0~dflat() +sdc~dunif(0,100) +}) +##################################################### +MultiSCinits<-list(uc=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),vc=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),ac0=-6.6,sdc=2.4, +uph=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),vph=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), +aph0=-1.5,sdph=0.1) + +MSCdata<-list( +yph=c(24,17,13,22,24,38,8,38,24,48,5,27,30,17,8, +29,18,24) , +yc=c(3,1,0,0,5,4,2,2,1,0,8,0,0,0,1,2,1,1,1,0,0,9,0,1,15,1,1,5,2,0,8,0,19,4,4,3,0,4,0,2,0,2,0,24,0,1,9,5,0,0,1,0,2,1,1,5,5,2,2,38,1,1,5,4,7,1,28,2,3,0,0,2,0,2,7,6,0,1,0,0,2,0,0,2,0,0,2,1,1,1,1,4,3,1,0,0,3,0,0,0,2,1,1,1,4, +7,7,1,0,4,3,2,2,3,6,1,0,0,1,0,18,3,0,0,0,6,1,1,2,1,0, +2,0,1,0,2,3,1,1,0,7,0,0,0,4,2,2,0,1,0,1,0,0,1,4,0,0,2,3)) +MSConsts<-list(nsumc = 860, +num2= c(6, 5, 5, 6, 5, 6, 6, 7, 5, 7, +6, 6, 6, 4, 5, 7, 5, 6, 6, 3, +5, 6, 2, 3, 2, 5, 3, 7, 5, 4, +5, 4, 5, 8, 6, 3, 5, 7, 5, 6, +1, 7, 5, 5, 6, 6, 6, 4, 5, 3, +4, 5, 10, 5, 5, 4, 4, 5, 4, 10, +6, 4, 4, 9, 3, 6, 7, 6, 9, 7, +3, 4, 3, 3, 6, 7, 5, 6, 6, 7, +8, 4, 6, 7, 5, 5, 7, 5, 5, 4, +4, 5, 6, 7, 4, 6, 7, 6, 7, 4, +7, 7, 5, 6, 4, 3, 6, 7, 7, 6, +5, 5, 5, 4, 4, 5, 6, 3, 2, 6, +4, 5, 4, 4, 3, 8, 3, 4, 8, 7, +6, 8, 7, 6, 6, 4, 6, 7, 4, 6, +4, 6, 6, 4, 7, 5, 6, 7, 6, 6, +7, 6, 6, 5, 4, 7, 6, 7, 7 +), +adj2 = c( +151, 138, 132, 113, 80, 3, +148, 86, 34, 32, 10, +148, 113, 80, 34, 1, +101, 100, 49, 47, 43, 19, +158, 150, 117, 84, 70, +127, 97, 78, 69, 68, 59, +147, 108, 78, 69, 67, 29, +115, 112, 110, 64, 57, 33, 28, +156, 142, 134, 77, 34, +137, 92, 86, 77, 37, 34, 2, +143, 111, 102, 84, 76, 39, +158, 143, 116, 87, 76, 45, +151, 148, 113, 63, 24, 20, +136, 92, 37, 35, +89, 54, 51, 25, 16, +124, 82, 54, 53, 51, 21, 15, +124, 121, 82, 81, 53, +126, 107, 102, 85, 79, 75, +135, 120, 49, 47, 30, 4, +63, 24, 13, +138, 132, 54, 53, 16, +110, 74, 71, 60, 48, 38, +155, 146, +148, 20, 13, +51, 15, +152, 130, 128, 106, 98, +146, 64, 57, +112, 64, 60, 58, 42, 33, 8, +109, 108, 97, 78, 7, +120, 118, 49, 19, +126, 75, 60, 56, 44, +148, 86, 50, 2, +110, 60, 48, 28, 8, +148, 134, 80, 77, 10, 9, 3, 2, +159, 137, 136, 101, 37, 14, +121, 94, 90, +137, 92, 35, 14, 10, +141, 126, 99, 74, 60, 56, 22, +145, 133, 111, 102, 11, +159, 156, 142, 129, 88, 46, +146, +112, 93, 69, 61, 58, 55, 28, +125, 101, 100, 65, 4, +122, 75, 67, 60, 31, +156, 153, 134, 116, 87, 12, +156, 129, 116, 96, 76, 40, +159, 135, 101, 88, 19, 4, +110, 60, 33, 22, +125, 100, 30, 19, 4, +92, 86, 32, +124, 25, 16, 15, +157, 109, 97, 90, 73, +140, 138, 132, 103, 83, 82, 81, 21, 17, 16, +132, 89, 21, 16, 15, +144, 105, 93, 61, 42, +126, 60, 38, 31, +115, 64, 27, 8, +69, 67, 60, 42, 28, +127, 97, 73, 6, +67, 58, 56, 48, 44, 38, 33, 31, 28, 22, +112, 105, 93, 64, 55, 42, +150, 149, 81, 70, +151, 95, 20, 13, +155, 146, 112, 105, 61, 57, 28, 27, 8, +136, 101, 43, +131, 117, 109, 108, 104, 70, +147, 122, 69, 60, 58, 44, 7, +154, 139, 127, 119, 69, 6, +154, 93, 78, 68, 67, 58, 42, 7, 6, +150, 149, 131, 117, 66, 62, 5, +115, 110, 22, +141, 130, 106, 99, +97, 59, 52, +141, 38, 22, +126, 122, 107, 44, 31, 18, +143, 116, 111, 96, 46, 12, 11, +142, 137, 34, 10, 9, +108, 97, 69, 29, 7, 6, +117, 107, 104, 102, 84, 18, +153, 138, 134, 103, 34, 3, 1, +150, 149, 121, 94, 83, 62, 53, 17, +124, 53, 17, 16, +158, 150, 140, 87, 81, 53, +158, 143, 117, 102, 79, 11, 5, +145, 126, 114, 102, 18, +92, 50, 32, 10, 2, +158, 153, 143, 140, 83, 45, 12, +159, 135, 129, 47, 40, +132, 95, 91, 54, 15, +157, 94, 52, 36, +151, 132, 95, 89, +86, 50, 37, 14, 10, +154, 144, 69, 61, 55, 42, +157, 149, 131, 121, 90, 81, 36, +151, 91, 89, 63, +133, 129, 123, 111, 76, 46, +109, 78, 73, 59, 52, 29, 6, +152, 133, 130, 129, 123, 26, +145, 141, 130, 126, 114, 72, 38, +125, 49, 43, 4, +159, 136, 65, 47, 43, 35, 4, +145, 85, 84, 79, 39, 18, 11, +153, 140, 138, 80, 53, +147, 117, 108, 107, 79, 66, +155, 64, 61, 55, +130, 72, 26, +147, 122, 104, 79, 75, 18, +147, 109, 104, 78, 66, 29, 7, +157, 131, 108, 97, 66, 52, 29, +115, 71, 48, 33, 22, 8, +133, 96, 76, 39, 11, +64, 61, 42, 28, 8, +151, 148, 13, 3, 1, +145, 126, 99, 85, +110, 71, 57, 8, +156, 76, 46, 45, 12, +104, 84, 79, 70, 66, 5, +128, 120, 30, +139, 68, +152, 135, 128, 118, 30, 19, +94, 81, 36, 17, +147, 107, 75, 67, 44, +133, 129, 98, 96, +82, 51, 17, 16, +100, 49, 43, +114, 99, 85, 75, 56, 38, 31, 18, +68, 59, 6, +152, 120, 118, 26, +152, 135, 123, 98, 96, 88, 46, 40, +145, 133, 106, 99, 98, 72, 26, +157, 149, 109, 94, 70, 66, +151, 138, 91, 89, 54, 53, 21, 1, +145, 130, 123, 111, 98, 96, 39, +156, 153, 80, 45, 34, 9, +152, 129, 120, 88, 47, 19, +101, 65, 35, 14, +159, 142, 77, 37, 35, 10, +140, 132, 103, 80, 53, 21, 1, +154, 144, 119, 68, +153, 138, 103, 87, 83, 53, +99, 74, 72, 38, +159, 156, 137, 77, 40, 9, +158, 87, 84, 76, 12, 11, +154, 139, 93, 55, +133, 130, 114, 102, 99, 85, 39, +155, 64, 41, 27, 23, +122, 108, 107, 104, 67, 7, +113, 34, 32, 24, 13, 3, 2, +157, 131, 94, 81, 70, 62, +158, 83, 81, 70, 62, 5, +132, 113, 95, 91, 63, 13, 1, +135, 129, 128, 120, 98, 26, +140, 134, 103, 87, 80, 45, +144, 139, 93, 69, 68, +146, 105, 64, 23, +142, 134, 116, 46, 45, 40, 9, +149, 131, 109, 94, 90, 52, +150, 143, 87, 84, 83, 12, 5, +142, 137, 101, 88, 47, 40, 35), +phc=c(18,18,18,16,12,4,3,1,15,15,12,11,18,15,17,18, +13,10,16,17,18,10,1,18,17,14,1,2,3,14,7,18,5,18, +16,13,15,10,12,14,1,4,16,9,11,14,16,5,16,15, +17,3,13,18,2,10,1,4,4,6,2,13,17,1,16,3,8,4,4, +12,1,14,4,10,10,12,15,3,12,18,13,13,11,12,10,15, +11,16,17,13,17,15,4,13,17,14,3,14,10,16,16,12,11, +3,2,14,8,3,3,1,12,2,18,10,1,11,12,14,4,14,13,8,14, +13,16,10,4,14,14,14,13,18,14,11,16,16,15,18,4,11, +10,15,12,4,10,1,3,18,13,12,18,14,11,4,2,11,13,12,16), +eph=c(26.7313582387354,17.8569786812963, +19.000833240556,24.8322642703346,35.6928573604747, +38.1880967220473,12.4232592469938,40.2840771546551, +31.6839976247491,33.1291602142528,6.65287755667638, +23.232040109711,20.3645029329909 +,16.4453885396069,10.949025543461,16.922717476408, +15.6506686820585,23.9598964049923), +ec=c(0.842405862334888,0.375626926592719,0.484362270840443,0.199184019993243,2.11970621276708,0.735452296043789,2.64537759886506,4.0780208537788,0.813194081624956,0.782106745171208,7.27574961919763,0.564870501143738,0.728747184259648 +,0.767430521475849,1.29108568514923,2.82946339512059,1.08730655358364,1.04852943858025,0.286116028719108,2.11506421230113,0.477938492417873,4.76283314473956,2.80606583721648,0.501617383683548,11.1838451225867,0.633281396899421 +,1.24508768053215,8.19055193324381,4.87414737813503,0.155530460055929,12.4232592469938,0.325830921594408,30.6655708558571,1.84643774089311,2.05199865041532,4.71650691786731,0.762178965393165,4.94096405150847,0.604304060657466 +,1.03286854811939,0.749847186377576,0.893889867502633,1.34172569023226,31.6839976247491,0.914380314003821,0.544098721280977,4.4863762281011,5.02728650461759,0.566933612461935,0.19229135263472,2.09410487686399,0.980352987292544 +,1.03591632620309,0.527406275161016,1.0134096572773,4.75139225470228,4.4079779980096,6.18300395395692,1.00590743430203,38.1880967220473,1.25451234814483,0.123364679049487,3.34585078028669,2.30116623098126,1.13846233649622 +,0.733904962555141,32.8594552982921,1.82763529456135,7.54559520183913,0.460026935064432,1.31612435432918,1.25605968163347,1.09574655443081,0.529375608692023,7.47905986182727,5.80264124911106,0.465888046763857,2.32287578871593 +,0.603272504998367,0.601115615892979,0.791625190571074,0.406995596408039,0.449476934005468,1.23013012347522,0.769446743900452,0.349931812902442,2.19008644205377,1.40258747411909,2.89529540172852,0.393772928414137,0.51240182921049 +,4.49134645082222,1.12181677926985,1.00890832349214,0.522248496865523,0.653396732251846,1.28062946187746,0.333473811250458,1.06672232929526,0.289070029015618,1.11773744552705,1.09851299915294,0.420593375550704,0.79767385784488 +,1.90162596865489,8.57363419391819,3.82257016147108,1.35696458065077,0.635672730472786,4.96722183192189,1.15651456053045,1.30215146403775,0.783982300915024,0.73850007412749,1.88807507840582,0.461246046297912,0.925868092934693 +,0.115674900499842,0.750738075355888,0.343742478947849,9.20264870150043,3.60205169489193,0.189478019018996,0.719088072178996,0.434566265842132,2.85487717544929,1.17165967316176,0.233553578998668,1.54137859916146,0.308857142112875 +,8.89013422568712E-02,1.07816321933254,0.421296708954635,0.605476282997351,0.513433384869588,2.06259554036788,1.88390196687583,1.25545012601673,0.47512515880215,0.330332255379566,2.86964717693184,0.440755599796724,0.489942049178295 +,0.919350536724933,1.31781235449861,2.97177118718262,3.3732338941464,1.66994794540004,0.293243140545608,0.987526988012639,1.32217302160298,0.109063566502892,0.30890403100647,1.10634344438337,4.19472730993774,0.407511374237588 +,0.496225160920078,0.477844714630683,1.03193077024748), +nsumph= 88, +num = c(3, 4, 5, 4, 6, 7, 4, 8, 6, 8, +5, 6, 5, 5, 4, 2, 2, 4), +adj = c( +10, 5, 2, +6, 5, 4, 1, +13, 12, 9, 8, 4, +8, 6, 3, 2, +10, 9, 8, 6, 2, 1, +10, 9, 8, 7, 5, 4, 2, +10, 9, 8, 6, +12, 10, 9, 7, 6, 5, 4, 3, +10, 8, 7, 6, 5, 3, +14, 12, 9, 8, 7, 6, 5, 1, +18, 15, 14, 13, 12, +14, 13, 11, 10, 8, 3, +18, 17, 12, 11, 3, +16, 15, 12, 11, 10, +18, 16, 14, 11, +15, 14, +18, 13, +17, 15, 13, 11)) + + +############################################################ +MSC<-nimbleModel(code=MultiSC,name="MSC",constants=MSConsts, +data=MSCdata,inits=MultiSCinits) +######## +CMSC<-compileNimble(MSC) +################## Running with monitors and CODA output ##### +LSTConf<-configureMCMC(MSC,print=TRUE) +LSTConf$addMonitors(c("ac0","uph","vph","uc","vc","thph","thc")) +LSTMCMC<-buildMCMC(LSTConf) +CLSTMCMC<-compileNimble(LSTMCMC,project=MSC) +niter<-20000 +set.seed(1) +samples<-runMCMC(CLSTMCMC,niter=niter,nburnin=15000,nchains=1,summary=TRUE,WAIC=TRUE,samplesAsCodaMCMC=TRUE) +#samples$summary +samples$WAIC + + +###################################################################### +sampMAT<-as.matrix(samples$samples) ### dim is 5000 by 535 +par(mfrow=c(2,2)) +plot(sampMAT[,"ac0"],type="l",xlab="iteration",ylab=expression(ac0)) +plot(sampMAT[,"aph0"],type="l",xlab="iteration",ylab=expression(aph0)) +plot(sampMAT[,"sdc"],type="l",xlab="iteration",ylab=expression(sdc)) +plot(sampMAT[,"sdph"],type="l",xlab="iteration",ylab=expression(sdph)) + + +setwd("C:/A_NEW_FOLDER/A_NEW_FOLDER/A_PROGRAMS/R_program/NIMBLE examples/CAR and other models on NIMBLE/Multi_scale/") +library(maptools) +library(sp) +source("fillmap.R") +GCpoly<-readSplus("Georgia_Counties_SPlus.txt") +plot(GCpoly) +GHDpoly<-readSplus("GeorgiaHDistrictsSPlus.txt") +plot(GHDpoly) + +fillmap(GCpoly,"",,n.col=5);x11() +fillmap(GCpoly,"",,n.col=5);x11() +fillmap(GCpoly,"",,n.col=5) \ No newline at end of file diff --git a/README.md b/README.md index 6a06eb6..239eab6 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,3 @@ # Bayesian_DM_Nimble_code This repositry is solely for new Nimble code for Bayesian disease mapping. Currently included are ICAR and convolution models, M models -and ST extensions. In addition the new MICAR and MPCAR code files for fitting MCAR models on Nimble are here +and ST extensions. In addition the new MICAR and MPCAR code files for fitting MCAR models on Nimble are here in the MCAR branch. diff --git a/SC_congen90_lognormal_Nimble_new.txt b/SC_congen90_lognormal_Nimble_new.txt new file mode 100644 index 0000000..60c9ca6 --- /dev/null +++ b/SC_congen90_lognormal_Nimble_new.txt @@ -0,0 +1,74 @@ +library(nimble) +LNCode<-nimbleCode({ +for(i in 1:N){ +log(theta[i])<-a0+v[i] +lambda[i]<-theta[i]*texp[i] +y[i]~dpois(lambda[i]) +v[i]~dnorm(0,tauv) +} +a0~dnorm(0,tau0) +tauv~dgamma(2,0.5) +tau0~dgamma(2,0.5) +}) +######################################################################## + +######################## SC data ####################################### + +LNConsts<-list(N=46,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)) +LNData<-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)) +LNInits<-list(a0=0.1,tauv=0.1,tau0=0.1,v=rep(0.0,LNConsts$N)) + +LN<-nimbleModel(code=LNCode,name="LN",constants=LNConsts, +data=LNData,inits=LNInits) +##################################################################### + +######################## Compile ###################################### +CLN<-compileNimble(LN) +CLN$v + + +mcmc.out<-nimbleMCMC(code=LNCode,constants=LNConsts, +data=LNData,inits=LNInits,nchains=2,niter=10000,summary=TRUE,WAIC=TRUE) + +mcmc.out$summary +mcmc.out$WAIC + +##################################################################### + +LNConf<-configureMCMC(LN,print=TRUE) +LNConf$addMonitors(c("a0","tauv","tau0","v","theta")) +LNMCMC<-buildMCMC(LNConf) +CLNMCMC<-compileNimble(LNMCMC,project=LN) +niter<-10000 +set.seed(1) +samples<-runMCMC(CLNMCMC,niter=niter,nburnin=7000) +plot(samples[,"a0"],type="l",xlab="iteration",ylab=expression(a0));x11() +plot(samples[,"tauv"],type="l",xlab="iteration",ylab=expression(tauv));x11() +plot(samples[,"tau0"],type="l",xlab="iteration",ylab=expression(tau0)) +boxplot(samples[,"a0"],xlab=expression(a0)) +boxplot(samples[,"tauv"],xlab=expression(tauv)) + +for(k in 1:46){x11() +plot(samples[,3+k],type="l",xlab="iteration",ylab=expression(theta)) +} +for(k in 1:46){x11() +plot(samples[,49+k],type="l",xlab="iteration",ylab=expression(v)) +} + +for (k in 1:46){ +boxplot(samples[,3+k],xlab=expression(theta))} + +samplesDF<-data.frame(samples) +plot(density(samplesDF$theta.1.)) diff --git a/SC_congen_90_CONV_Nimble.txt b/SC_congen_90_CONV_Nimble.txt new file mode 100644 index 0000000..af386a7 --- /dev/null +++ b/SC_congen_90_CONV_Nimble.txt @@ -0,0 +1,168 @@ +#### 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)) + +Lconv<-nimbleModel(code=LconvCode,name="Lconv",constants=LconvConsts, +data=LconvData,inits=LconvInits) +##################################################################### + +######################## Compile ###################################### +CLconv<-compileNimble(Lconv) +CLconv$v + + +mcmc.out<-nimbleMCMC(code=LconvCode,constants=LconvConsts, +data=LconvData,inits=LconvInits,nchains=2,nburnin=7000,niter=10000,summary=TRUE,WAIC=TRUE) + +mcmc.out$summary +mcmc.out$WAIC + +##################################################################### + +LconvConf<-configureMCMC(Lconv,print=TRUE,enableWAIC=nimbleOptions('enableWAIC'=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 ruun without CODA ######################################## + +samples2<-runMCMC(CLconvMCMC,niter=niter,nburnin=7000,summary=TRUE,WAIC=TRUE) +samples2$WAIC +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)) +x11();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)) +x11() +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) \ No newline at end of file diff --git a/STinT_models_Ohio_Nimble.txt b/STinT_models_Ohio_Nimble.txt new file mode 100644 index 0000000..e6832a6 --- /dev/null +++ b/STinT_models_Ohio_Nimble.txt @@ -0,0 +1,632 @@ +#### range of ST interaction models in Nimble ######## +library(nimble) + +###################################################### +###### Ohio Data #################################### +e=structure(.Data=c(11.08265, 12.43046, 12.58617, 13.07113, 13.64872, 13.44648, 14.20172, 14.38548, 14.95299, 15.44726, +51.54766, 57.20417, 57.04018, 58.30136, 61.00702, 60.78989, 64.1098, 64.32432, 66.04699, 67.65539, +21.18653, 23.58796, 23.70401, 24.58979, 25.74367, 25.60512, 27.00496, 26.73573, 27.48413, 28.12705, +48.02128, 53.02829, 52.80107, 54.51591, 56.9513, 56.11595, 58.9642, 58.08966, 58.91479, 59.60639, +25.54036, 28.85985, 29.14606, 30.34142, 31.9459, 31.69799, 33.39785, 33.45538, 34.54226, 35.47209, +19.65611, 21.73145, 21.82331, 22.744, 23.8938, 23.95328, 25.31747, 25.35509, 26.03825, 26.72759, +38.56623, 42.0636, 41.97363, 43.43165, 45.88792, 44.81808, 46.23249, 44.58007, 45.22339, 45.36022, +14.60163, 16.3542, 16.67166, 17.49995, 18.45894, 18.54305, 19.89741, 19.88945, 20.67724, 21.42489, +118.7359, 132.3991, 133.8755, 139.09, 145.8855, 146.4967, 155.7944, 157.4248, 163.6123, 169.2534, +11.82541, 13.03703, 13.08939, 13.6795, 14.60024, 14.62257, 15.47169, 15.4304, 15.81459, 16.15364, +15.67667, 17.15562, 17.08506, 17.82085, 18.65204, 18.47113, 19.61539, 19.92314, 20.61976, 21.28398, +69.8358, 76.60449, 76.37443, 78.53021, 82.3222, 81.24256, 85.3237, 85.31212, 87.24273, 89.11415, +57.71282, 65.82852, 67.17506, 70.67236, 74.95687, 75.22547, 80.34092, 81.57098, 85.39924, 88.77426, +16.01126, 17.66246, 17.66381, 18.29873, 19.23135, 19.14769, 20.11954, 20.21182, 20.97767, 21.57185, +52.67453, 57.85574, 57.74398, 60.01879, 62.86696, 61.75908, 64.5058, 62.88037, 64.11284, 64.84256, +16.55965, 18.39023, 18.49928, 19.18715, 20.17, 20.04829, 20.82488, 20.73979, 21.24078, 21.6668, +23.48562, 25.5304, 25.47703, 26.19371, 27.21013, 27.18946, 28.543, 28.23602, 28.82986, 29.35596, +702.4902, 762.9026, 756.0005, 779.8555, 822.0158, 809.6126, 845.5759, 844.0905, 855.8276, 869.7341, +25.70603, 28.01751, 27.7992, 28.40706, 29.93484, 29.77232, 31.3417, 31.05596, 31.78028, 32.41493, +18.48065, 20.33881, 20.22841, 20.71682, 21.59112, 21.23599, 22.69088, 22.68849, 23.27034, 23.87908, +24.805, 27.5521, 27.93598, 29.38989, 31.19252, 31.27092, 33.3397, 34.18607, 35.54015, 36.99492, +36.24261, 40.57072, 40.45171, 41.55404, 43.33169, 42.92893, 44.9724, 44.73515, 45.64945, 46.47362, +42.30873, 47.8858, 48.49051, 50.51246, 53.0288, 52.6939, 55.83695, 56.66023, 58.79983, 60.75547, +12.7047, 13.97347, 13.9019, 14.48464, 15.41407, 15.19235, 15.95782, 15.92759, 16.30998, 16.69189, +406.4846, 444.3345, 447.9958, 467.6531, 495.0504, 493.6799, 522.6015, 529.1328, 546.3256, 564.1479, +17.45652, 19.26888, 19.32098, 20.00927, 21.10864, 21.07612, 22.27337, 21.93282, 22.58947, 23.08501, +13.63535, 15.36918, 15.39472, 15.77881, 16.67008, 16.62236, 17.48306, 17.37735, 17.9514, 18.41551, +34.03191, 37.97007, 38.03563, 39.53481, 41.63126, 41.20462, 43.48262, 44.34541, 45.61567, 47.0657, +60.06467, 66.19062, 66.30437, 69.08276, 72.43422, 71.29946, 75.27091, 75.69876, 78.17459, 80.2487, +19.00359, 21.44709, 21.43185, 22.10858, 23.18576, 22.81702, 23.65325, 23.04279, 23.46055, 23.6662, +407.7868, 445.5072, 444.9867, 461.4336, 484.2658, 478.9633, 503.4001, 507.0756, 517.8945, 530.1678, +29.92659, 32.92372, 33.0029, 34.45148, 36.3201, 36.24579, 38.31725, 37.87029, 38.78331, 39.59003, +14.95705, 16.66925, 16.62521, 17.09685, 17.73578, 17.54398, 18.44368, 18.25615, 18.80826, 19.24466, +8.223136, 9.199208, 9.044219, 9.307426, 9.734615, 9.357822, 9.683545, 9.43737, 9.2791, 9.210176, +13.21699, 14.46189, 14.36684, 14.86706, 15.76502, 15.66644, 16.51372, 16.39807, 16.78344, 17.14971, +15.48276, 17.13619, 17.35708, 17.92799, 18.85522, 18.89875, 19.95381, 20.02827, 20.82716, 21.48658, +11.08913, 12.39415, 12.39784, 12.78047, 13.67559, 13.66001, 14.29243, 14.34889, 14.86173, 15.29123, +13.66543, 15.03009, 15.261, 15.9618, 16.94154, 16.78611, 17.82091, 17.68577, 18.01362, 18.35322, +25.26315, 27.883, 27.92015, 28.84994, 30.42291, 30.3825, 31.86446, 32.1276, 33.03179, 33.8531, +14.05833, 15.58909, 15.45903, 15.98089, 16.97009, 16.77393, 17.66914, 17.47493, 17.97451, 18.35927, +42.13704, 46.5668, 46.22499, 47.71302, 49.26917, 48.36682, 50.33897, 48.49604, 48.90682, 49.00401, +21.60812, 23.64268, 23.77036, 24.95947, 26.3241, 26.09414, 27.57889, 27.29623, 28.12766, 28.81045, +98.29676, 108.4951, 108.9542, 113.7492, 119.6884, 117.9353, 123.3566, 123.7011, 126.4734, 129.3168, +29.23195, 32.54986, 32.50121, 33.32862, 35.20066, 34.97511, 36.30995, 36.13476, 36.90544, 37.5362, +55.50258, 61.91143, 62.72108, 65.61343, 69.11733, 68.54123, 72.42336, 72.3067, 74.21561, 76.06302, +17.98687, 19.98131, 20.0411, 20.71045, 21.93815, 21.89872, 23.17992, 23.14502, 23.98498, 24.67195, +126.3616, 139.9924, 139.737, 144.8432, 151.641, 149.4939, 156.9894, 155.4053, 158.835, 161.6906, +218.5626, 240.4953, 239.8109, 248.0714, 260.3825, 256.0689, 269.0686, 269.1015, 274.4122, 280.0392, +15.02646, 16.88968, 17.11976, 18.04467, 19.3265, 18.94633, 20.13931, 20.34076, 20.98775, 21.63898, +133.9614, 147.3142, 146.1641, 150.3313, 158.2776, 155.0762, 161.9477, 159.0297, 161.3286, 163.2364, +31.10297, 34.60226, 34.48909, 35.66503, 37.35222, 36.75197, 38.38354, 37.69953, 38.50183, 39.11044, +51.48056, 57.77391, 58.37173, 60.84568, 64.23379, 63.78598, 67.20449, 66.97814, 69.27415, 71.1111, +10.62496, 12.03359, 12.02476, 12.55771, 13.17856, 13.06422, 13.79991, 13.74192, 14.07124, 14.38829, +17.7365, 19.55017, 19.57361, 20.45003, 21.51948, 21.33668, 22.57982, 22.33302, 22.99301, 23.53497, +41.81355, 46.1152, 46.15864, 47.72734, 49.8552, 49.26023, 51.80665, 52.06061, 53.43704, 54.81957, +7.922329, 8.841202, 8.741061, 9.016768, 9.394306, 9.198502, 9.534684, 9.128947, 9.167696, 9.14244, +264.8869, 291.6328, 290.8818, 300.4148, 314.7764, 312.2435, 328.7067, 329.3978, 337.411, 345.2845, +6.43588, 7.257283, 7.257423, 7.573024, 7.930644, 7.847051, 8.22924, 8.173478, 8.375426, 8.532823, +12.47516, 13.51164, 13.67121, 14.16481, 14.84484, 14.80125, 15.68975, 15.67493, 16.16302, 16.61811, +38.55466, 42.55919, 42.73356, 44.77462, 46.89149, 46.67626, 48.96491, 48.30727, 49.45495, 50.26619, +5.165549, 5.765423, 5.743163, 5.979708, 6.258207, 6.210151, 6.660381, 6.547146, 6.696667, 6.819485, +18.76665, 20.43444, 20.47185, 21.09605, 22.19562, 21.92915, 23.11014, 23.05906, 23.47892, 23.98129, +9.838238, 10.83785, 10.78763, 11.07418, 11.64549, 11.46604, 12.08742, 12.00175, 12.39663, 12.69068, +14.30453, 15.85453, 16.0189, 16.63488, 17.36357, 17.39849, 18.35704, 18.34967, 18.87226, 19.35897, +20.39795, 22.21271, 22.12289, 22.75408, 23.67607, 23.56273, 26.124, 26.47552, 27.2708, 28.31332, +10.36488, 11.71599, 12.09417, 12.92527, 13.78977, 13.75626, 14.58608, 14.28558, 14.77699, 15.10677, +62.13098, 69.43773, 70.19184, 72.84647, 76.47707, 75.42849, 79.18259, 79.55898, 81.5173, 83.49575, +17.51529, 19.46834, 19.43428, 20.05223, 21.26592, 21.41523, 22.6589, 22.46312, 23.19864, 23.76841, +15.89973, 16.85797, 17.00544, 17.59914, 18.50092, 18.49658, 19.38861, 19.35567, 19.87964, 20.35625, +60.26829, 66.856, 66.74379, 68.72845, 71.79167, 71.5495, 75.14473, 74.25946, 75.812, 77.07059, +30.01128, 33.21524, 33.46478, 34.83443, 36.98001, 37.29299, 39.37498, 38.87397, 40.03778, 40.8903, +29.4962, 32.23379, 32.03269, 33.0141, 34.59225, 34.50268, 36.09945, 35.85596, 36.62397, 37.28642, +39.07806, 43.08545, 42.90556, 44.69506, 47.42099, 46.62039, 48.37818, 47.44009, 48.39128, 48.90422, +29.00427, 31.49937, 31.23754, 32.38823, 34.20381, 33.9683, 35.91977, 35.55103, 36.51612, 37.27735, +19.92545, 22.00916, 22.19128, 23.07338, 24.07907, 24.00417, 25.44191, 25.45151, 26.11824, 26.75481, +174.7026, 193.1358, 193.0245, 200.3059, 210.1053, 207.67, 217.5818, 215.3065, 219.1714, 222.7557, +243.7286, 267.0859, 265.44, 274.3023, 287.5316, 283.302, 295.7671, 295.0601, 301.478, 307.4574, +111.7262, 123.315, 123.1225, 127.2642, 132.7925, 130.4624, 136.5146, 134.1485, 135.9463, 137.436, +38.89249, 43.20411, 43.4246, 45.23235, 47.70868, 47.30137, 49.71968, 48.81783, 49.89583, 50.72703, +13.60943, 15.1017, 15.23089, 16.07265, 16.88501, 16.89564, 17.95465, 18.30668, 18.94278, 19.6166, +14.08841, 15.50726, 15.41871, 15.88489, 16.62306, 16.54049, 17.44178, 17.14676, 17.56031, 17.88089, +5.246998, 5.882542, 5.80798, 6.009941, 6.37127, 6.390492, 6.649333, 6.685965, 6.812811, 6.950722, +45.4547, 50.73245, 51.18831, 53.51505, 56.3983, 56.37872, 59.76026, 60.88114, 63.63759, 66.06119, +29.50638, 32.76467, 32.82683, 34.19052, 36.21935, 35.71362, 37.41943, 36.98278, 37.74334, 38.37079, +45.04144, 49.70088, 50.01446, 52.3609, 55.26991, 54.94208, 58.23152, 58.24706, 60.07682, 61.70376, +16.69849, 18.5631, 18.56716, 19.11343, 20.07093, 20.06876, 21.27903, 20.93902, 21.5981, 22.07382, +49.07087, 54.90015, 55.5739, 57.83726, 60.39749, 60.01984, 63.66206, 63.51696, 65.07695, 66.56497, +10.55323, 11.56767, 11.5246, 11.92865, 12.53544, 12.49332, 13.16783, 12.85789, 13.16342, 13.39404),.Dim=c(10,88)) +e<-t(e) +LSTConsts<-list(L=463,m=88,e, +m=88, T=10,t=c(1,2,3,4,5,6,7,8,9,10), +sumNumNeigh=463, +adj=c( +8, 36, 66, 73, +6, 81, 69, 32, 33, +42, 70, 39, 47, 52, 85, 38, +43, 28, 78, +53, 82, 37, 64, 58, 84, +75, 19, 54, 81, 2, 33, 46, +56, 61, 30, 34, 41, +13, 14, 36, 1, +31, 68, 57, 83, +34, 79, 76, 15, 41, +12, 55, 75, 46, 80, 49, +29, 57, 55, 11, 49, +31, 83, 14, 8, +8, 13, 83, 29, 24, 36, +41, 10, 76, 50, +60, 45, 42, 38, 79, 30, +59, 51, 88, 74, 39, 70, +47, 52, 77, 67, 28, 43, +68, 57, 55, 75, 6, 54, +86, 35, 69, 63, +25, 80, 51, 59, 42, 45, +62, 72, 39, 47, +65, 25, 45, 64, 37, +36, 14, 29, 49, 65, 71, +65, 49, 80, 21, 45, 23, +86, 35, 48, +44, 40, 82, 53, +67, 77, 18, 43, 4, 78, +14, 83, 57, 12, 49, 24, +61, 60, 16, 79, 34, 7, +9, 83, 13, +33, 2, 69, 35, 87, 74, 88, +80, 46, 6, 2, 32, 88, 51, +7, 30, 79, 10, 41, +32, 69, 20, 86, 26, 48, 87, +1, 8, 14, 24, 71, 66, +82, 71, 65, 23, 64, 5, +16, 42, 3, 85, 85, 76, 79, +70, 17, 74, 72, 22, 47, 3, +44, 73, 66, 71, 82, 27, +7, 34, 10, 15, +45, 21, 59, 70, 3, 38, 16, +18, 28, 4, +73, 40, 27, +64, 23, 25, 21, 42, 16, 60, +11, 75, 6, 33, 80, +22, 39, 3, 52, 18, +26, 35, 87, 62, +24, 29, 12, 11, 80, 25, 65, +15, 76, 67, 78, +21, 80, 33, 88, 17, 59, +85, 3, 47, 18, 77, +27, 82, 5, +19, 6, 81, +57, 19, 75, 11, 12, +84, 61, 7, +83, 9, 68, 19, 55, 12, 29, +5, 64, 60, 61, 84, +21, 51, 17, 70, 42, +58, 64, 45, 16, 30, 61, +58, 60, 30, 7, 56, 84, +48, 87, 72, 22, +81, 69, 20, +37, 23, 45, 60, 58, 5, +71, 24, 49, 25, 23, 37, +73, 1, 36, 71, 40, +76, 77, 18, 28, 78, 50, +9, 57, 19, +2, 81, 63, 20, 35, 87, 32, +42, 59, 17, 39, 3, +66, 36, 24, 65, 37, 82, 40, +74, 87, 62, 22, 39, +1, 66, 40, 44, +17, 88, 32, 87, 72, 39, +55, 19, 6, 46, 11, +10, 79, 38, 85, 77, 67, 50, 15, +76, 85, 52, 18, 28, 67, +50, 67, 28, 4, +30, 16, 38, 76, 10, 34, +49, 11, 46, 33, 51, 21, 25, +54, 6, 2, 69, 63, +40, 71, 37, 5, 53, 27, +13, 31, 9, 57, 29, 14, +5, 58, 61, 56, +38, 3, 52, 77, 76, +20, 35, 26, +32, 69, 35, 48, 62, 72, 74, +51, 33, 32, 74, 17), +num=c(4,5,7,3,6,7,5,4,4,5,6,5,4,6,4,6,6,6,6, +4,6,4,5,6,6,3,4,6,6,6,3,7,7,5,7,6,6,7,7,6, +4,7,3,3,7,5,5,4,7,4,6, 5,3,3,5,3,7,5,5,6, +6,4,3,6,6,5,6,3,7,5,7,5,4,6,5,8,6,4,6, +7,5,6,6,4,5,3,7,5)) +names(LSTConsts)[3]<-'e' +y=structure(.Data=c(5, 14, 12, 15, 9, 12, 20, 12, 16, 15, +76, 46, 53, 47, 62, 69, 53, 65, 69, 58, +11, 22, 18, 20, 27, 25, 34, 19, 28, 30, +38, 58, 56, 53, 59, 56, 56, 50, 53, 65, +30, 25, 33, 24, 21, 38, 21, 15, 27, 24, +18, 17, 20, 23, 16, 35, 22, 23, 19, 22, +40, 54, 64, 47, 45, 55, 44, 50, 62, 70, +7, 21, 12, 22, 27, 24, 19, 23, 30, 23, +91, 135, 129, 138, 150, 141, 139, 164, 159, 158, +12, 7, 12, 14, 16, 17, 11, 8, 19, 10, +12, 17, 14, 20, 19, 18, 19, 19, 11, 21, +60, 99, 80, 88, 73, 87, 92, 95, 114, 105, +51, 62, 71, 59, 51, 75, 94, 80, 89, 75, +11, 11, 16, 15, 19, 17, 26, 19, 22, 24, +55, 51, 63, 64, 62, 74, 65, 68, 67, 60, +13, 22, 13, 13, 18, 15, 11, 28, 13, 19, +23, 28, 21, 37, 21, 26, 35, 30, 24, 31, +825, 917, 913, 909, 940, 923, 937, 953, 952, 993, +22, 18, 18, 23, 23, 24, 37, 28, 28, 35, +15, 12, 13, 12, 10, 16, 12, 25, 16, 22, +25, 17, 10, 24, 27, 14, 24, 23, 28, 22, +33, 41, 32, 36, 40, 42, 40, 43, 32, 54, +41, 33, 38, 48, 53, 33, 54, 51, 54, 45, +15, 14, 22, 15, 15, 14, 16, 23, 17, 18, +402, 402, 431, 451, 478, 444, 524, 499, 488, 543, +12, 9, 15, 8, 16, 18, 9, 7, 16, 11, +20, 15, 17, 12, 15, 10, 6, 16, 19, 29, +20, 22, 29, 36, 25, 20, 20, 38, 28, 30, +49, 44, 44, 36, 57, 50, 71, 50, 60, 52, +17, 22, 24, 31, 32, 19, 22, 18, 33, 27, +528, 583, 542, 554, 568, 642, 639, 572, 593, 576, +24, 29, 22, 24, 41, 27, 31, 33, 36, 29, +11, 18, 15, 18, 25, 15, 26, 19, 20, 21, +7, 8, 12, 20, 22, 7, 7, 16, 12, 14, +6, 7, 13, 12, 2, 12, 23, 16, 9, 9, +24, 25, 13, 17, 24, 24, 25, 19, 21, 24, +9, 14, 14, 15, 15, 18, 18, 10, 17, 12, +4, 3, 8, 9, 5, 6, 8, 9, 6, 8, +26, 25, 28, 23, 18, 21, 21, 33, 23, 16, +21, 22, 12, 19, 14, 14, 15, 15, 26, 18, +57, 55, 68, 65, 54, 68, 56, 71, 68, 77, +25, 15, 22, 27, 28, 26, 26, 23, 27, 37, +87, 101, 83, 94, 114, 103, 127, 111, 108, 147, +35, 39, 26, 45, 36, 41, 37, 33, 46, 49, +57, 69, 72, 57, 68, 84, 61, 76, 61, 81, +11, 19, 31, 24, 24, 20, 25, 24, 30, 25, +121, 115, 110, 136, 125, 124, 124, 137, 136, 149, +280, 243, 282, 277, 292, 243, 308, 305, 290, 305, +10, 10, 13, 14, 24, 20, 21, 12, 15, 16, +137, 142, 159, 163, 148, 147, 183, 200, 163, 187, +25, 27, 36, 41, 37, 26, 40, 50, 40, 55, +27, 38, 40, 36, 34, 41, 45, 31, 40, 42, +19, 12, 14, 12, 9, 14, 17, 15, 19, 16, +16, 12, 15, 14, 14, 9, 15, 15, 12, 12, +44, 48, 37, 49, 63, 39, 48, 47, 44, 42, +7, 4, 7, 6, 7, 7, 6, 6, 8, 11, +259, 314, 306, 325, 345, 326, 341, 374, 361, 340, +7, 8, 8, 4, 7, 6, 11, 11, 13, 11, +9, 16, 16, 6, 13, 9, 12, 10, 14, 10, +49, 54, 59, 44, 59, 62, 60, 45, 75, 66, +6, 3, 2, 5, 5, 10, 9, 7, 2, 5, +10, 15, 21, 18, 26, 25, 28, 14, 26, 26, +13, 11, 15, 8, 9, 11, 6, 11, 14, 9, +10, 21, 15, 25, 20, 23, 23, 25, 14, 14, +27, 23, 15, 23, 21, 24, 22, 23, 26, 27, +8, 14, 10, 10, 20, 17, 8, 11, 13, 23, +38, 37, 49, 66, 55, 65, 58, 44, 62, 55, +13, 24, 13, 21, 27, 21, 22, 26, 18, 18, +11, 13, 14, 9, 13, 10, 12, 16, 14, 11, +45, 62, 66, 63, 59, 66, 50, 79, 79, 67, +27, 23, 33, 28, 34, 36, 29, 25, 36, 42, +23, 27, 23, 26, 24, 25, 34, 33, 22, 32, +41, 54, 38, 56, 44, 70, 56, 71, 58, 67, +19, 30, 17, 28, 36, 39, 25, 31, 41, 40, +9, 19, 21, 19, 30, 23, 22, 21, 16, 22, +182, 197, 183, 188, 206, 198, 205, 217, 244, 224, +247, 277, 253, 282, 310, 302, 292, 293, 323, 337, +91, 114, 112, 111, 133, 114, 139, 143, 146, 146, +28, 36, 40, 54, 57, 37, 58, 49, 38, 51, +4, 17, 9, 9, 14, 11, 10, 14, 12, 6, +9, 13, 12, 17, 13, 17, 13, 15, 17, 10, +6, 5, 7, 7, 13, 3, 6, 7, 9, 4, +43, 33, 37, 40, 42, 42, 55, 40, 56, 59, +20, 26, 34, 30, 32, 39, 27, 40, 40, 40, +36, 27, 28, 30, 30, 31, 46, 38, 45, 32, +10, 17, 23, 10, 13, 17, 17, 13, 22, 17, +28, 32, 30, 30, 44, 37, 55, 40, 45, 30, +7, 7, 11, 5, 10, 16, 10, 11, 21, 16),.Dim=c(10,88)) +y<-t(y) +LSTdata<-list(y) +names(LSTdata)<-'y' +########################################################## +###### inits ############################################## +LSTinits<-list(a0=0.1,sdv=0.1,sdpsi=0.1,sdg=0.1,sdu=0.1,g=c(0,0,0,0,0,0,0,0,0,0), +v=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0), +u=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0), +psi=structure(.Data=rep(0.1,880),.Dim=c(88,10))) + +########################################################### +########## model running ################################ + +LST<-nimbleModel(code=STCode,name="LST",constants=LSTConsts, +data=LSTdata,inits=LSTinits) +##################################################################### +######################## Compile ###################################### +CLST<-compileNimble(LST) +CLST$v +################## Running with monitors and CODA output ##### +LSTConf<-configureMCMC(LST,monitors=CLST$getNodeNames(stochOnly=TRUE,includeData=FALSE),print=TRUE) +##LSTConf$addMonitors(c("a0","g","u","v")) +LSTMCMC<-buildMCMC(LSTConf) +CLSTMCMC<-compileNimble(LSTMCMC,project=LST) +niter<-20000 +set.seed(1) +samples<-runMCMC(CLSTMCMC,niter=niter,nburnin=15000,nchains=1,summary=TRUE,WAIC=TRUE,samplesAsCodaMCMC=TRUE) +#samples$summary +samples$WAIC + +############################################################## +##### models for Ohio lung cancer ########################### +### KH model ############################################## +STCode<-nimbleCode( +{ +for (i in 1:m) +{ + for (k in 1:T) + {y[i,k]~dpois(mu[i,k]) + log(mu[i,k])<-log(e[i,k])+log(theta[i,k]) +log(theta[i,k])<-a0+g[k]+u[i]+v[i]+psi[i,k] +psi[i,k]~dnorm(0,taupsi) + } +v[i]~dnorm(0,tauv) +} +g[1]~dnorm(0,taug) +for (k in 2:T){ + g[k]~dnorm(g[k-1],taug) + } + u[1:m]~dcar_normal(adj[1:L],weights[1:L],num[1:m],tauu) + for(k in 1:sumNumNeigh) { + weights[k] <- 1} +a0~dflat() +tauu<-1/pow(sdu,2) +sdu~dunif(0,10) +taug<-1/pow(sdg,2) +sdg~dunif(0,10) +taupsi<-1/pow(sdpsi,2) +sdpsi~dunif(0,2) +tauv<-1/pow(sdv,2) +sdv~dunif(0,10) +for(k in 1:T){t1[k]<-t[k]} +}) +LSTinits<-list(a0=0.1,sdv=0.1,sdpsi=0.1,sdg=0.1,sdu=0.1,g=c(0,0,0,0,0,0,0,0,0,0), +v=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0), +u=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0), +psi=structure(.Data=rep(0.1,880),.Dim=c(88,10))) +######################################################### +############## model with no interaction ############### +STCode<-nimbleCode( +{ +for (i in 1:m) +{ + for (k in 1:T) + {y[i,k]~dpois(mu[i,k]) + log(mu[i,k])<-log(e[i,k])+log(theta[i,k]) +log(theta[i,k])<-a0+g[k]+u[i]+v[i] + } +v[i]~dnorm(0,tauv) +} +g[1]~dnorm(0,taug) +for (k in 2:T){ + g[k]~dnorm(g[k-1],taug) + } + u[1:m]~dcar_normal(adj[1:L],weights[1:L],num[1:m],tauu) + for(k in 1:sumNumNeigh) { + weights[k] <- 1} +a0~dflat() +tauu<-1/pow(sdu,2) +sdu~dunif(0,10) +taug<-1/pow(sdg,2) +sdg~dunif(0,10) +tauv<-1/pow(sdv,2) +sdv~dunif(0,10) +for(k in 1:T){t1[k]<-t[k]} +}) +LSTinits<-list(a0=0.1,sdv=0.1,sdg=0.1,sdu=0.1,g=c(0,0,0,0,0,0,0,0,0,0), +v=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0), +u=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0)) +LST<-nimbleModel(code=STCode,name="LST",constants=LSTConsts, +data=LSTdata,inits=LSTinits) +##################################################################### +######################## Compile ###################################### +CLST<-compileNimble(LST) +#CLST$v +################## Running with monitors and CODA output ##### +LSTConf<-configureMCMC(LST,monitors=CLST$getNodeNames(stochOnly=TRUE,includeData=FALSE),print=TRUE) +##LSTConf$addMonitors(c("a0","u","v","theta")) +LSTMCMC<-buildMCMC(LSTConf) +CLSTMCMC<-compileNimble(LSTMCMC,project=LST) +niter<-20000 +set.seed(1) +samples<-runMCMC(CLSTMCMC,niter=niter,nburnin=15000,nchains=1,summary=TRUE,WAIC=TRUE,samplesAsCodaMCMC=TRUE) +samples$summary +samples$WAIC +############################################################ +################# model with psi and no time ############## +STCode<-nimbleCode( +{ +for (i in 1:m) +{ + for (k in 1:T) + {y[i,k]~dpois(mu[i,k]) + log(mu[i,k])<-log(e[i,k])+log(theta[i,k]) +log(theta[i,k])<-a0+u[i]+v[i]+psi[i,k] +psi[i,k]~dnorm(0,taupsi) + } +v[i]~dnorm(0,tauv) +} + u[1:m]~dcar_normal(adj[1:L],weights[1:L],num[1:m],tauu) + for(k in 1:sumNumNeigh) { + weights[k] <- 1} +a0~dflat() +tauu<-1/pow(sdu,2) +sdu~dunif(0,10) +taupsi<-1/pow(sdpsi,2) +sdpsi~dunif(0,2) +tauv<-1/pow(sdv,2) +sdv~dunif(0,10) +for(k in 1:T){t1[k]<-t[k]} +}) +LSTinits<-list(a0=0.1,sdv=0.1,sdpsi=0.1,sdu=0.1, +v=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0), +u=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0), +psi=structure(.Data=rep(0.1,880),.Dim=c(88,10))) +########################################################### +################ Space-time trend ######################### +STCode<-nimbleCode( +{ +for (i in 1:m) +{ + for (k in 1:T) + {y[i,k]~dpois(mu[i,k]) + log(mu[i,k])<-log(e[i,k])+log(theta[i,k]) +log(theta[i,k])<-a0+g[k]+u[i]+v[i]+psi[i,k] +psi[i,k]~dnorm(0,taupsi) + } +v[i]~dnorm(0,tauv) +} +for (k in 1:T){ + g[k]~dnorm(beta*t[k],taug) +} +u[1:m]~dcar_normal(adj[1:L],weights[1:L],num[1:m],tauu) +for(k in 1:sumNumNeigh) {weights[k] <- 1} +a0~dnorm(0,tau0) +tauu<-1/pow(sdu,2) +sdu~dunif(0,10) +taug<-1/pow(sdg,2) +sdg~dunif(0,10) +taupsi<-1/pow(sdpsi,2) +sdpsi~dunif(0,2) +tauv<-1/pow(sdv,2) +sdv~dunif(0,10) +beta~dnorm(0,tauB) +tauB<-1/pow(sdB,2) +sdB~dunif(0,10) +}) +LSTinits<-list(a0=0.1,sdv=0.1,sdpsi=0.1,sdg=0.1,sdu=0.1,g=c(0,0,0,0,0,0,0,0,0,0), +v=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0), +u=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0), +psi=structure(.Data=rep(0.1,880),.Dim=c(88,10))) + +######################################################## +############# g UH and ST int ######################### + +STCode<-nimbleCode( +{ +for (i in 1:m) +{ + for (k in 1:T) + {y[i,k]~dpois(mu[i,k]) + log(mu[i,k])<-log(e[i,k])+log(theta[i,k]) +log(theta[i,k])<-a0+g[k]+u[i]+v[i]+psi[i,k] +psi[i,k]~dnorm(0,taupsi) + } +v[i]~dnorm(0,tauv) +} +for (k in 1:T){ + g[k]~dnorm(0,taug) +} +u[1:m]~dcar_normal(adj[1:L],weights[1:L],num[1:m],tauu) +for(k in 1:sumNumNeigh) {weights[k] <- 1} +a0~dnorm(0,tau0) +tauu<-1/pow(sdu,2) +sdu~dunif(0,10) +taug<-1/pow(sdg,2) +sdg~dunif(0,10) +taupsi<-1/pow(sdpsi,2) +sdpsi~dunif(0,2) +tauv<-1/pow(sdv,2) +sdv~dunif(0,10) + +}) +LSTinits<-list(a0=0.1,sdv=0.1,sdpsi=0.1,sdg=0.1,sdu=0.1,g=c(0,0,0,0,0,0,0,0,0,0), +v=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0), +u=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0), +psi=structure(.Data=rep(0.1,880),.Dim=c(88,10))) + +########################################################## +########### CONV + g trend ############################ +STCode<-nimbleCode( +{ +for (i in 1:m) +{ + for (k in 1:T) + {y[i,k]~dpois(mu[i,k]) + log(mu[i,k])<-log(e[i,k])+log(theta[i,k]) +log(theta[i,k])<-a0+g[k]+u[i]+v[i] +# + } +v[i]~dnorm(0,tauv) +} +g[1]~dnorm(0,taug) +for (k in 2:T){ + g[k]~dnorm(g[k-1],taug) +} +u[1:m]~dcar_normal(adj[1:L],weights[1:L],num[1:m],tauu) +for(k in 1:sumNumNeigh) {weights[k] <- 1} +a0~dnorm(0,tau0) +tauu<-1/pow(sdu,2) +sdu~dunif(0,10) +taug<-1/pow(sdg,2) +sdg~dunif(0,10) +tauv<-1/pow(sdv,2) +sdv~dunif(0,10) +tau0~dgamma(2,0.5) +}) +LSTinits<-list(a0=0.1,sdv=0.1,tau0=0.1,sdg=0.1,sdu=0.1,g=c(0,0,0,0,0,0,0,0,0,0), +v=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0), +u=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0)) + +########################################################### +############## g RW1 +conv ############################### +########################################################### +STCode<-nimbleCode( +{ +for (i in 1:m) +{ + for (k in 1:T) + {y[i,k]~dpois(mu[i,k]) + log(mu[i,k])<-log(e[i,k])+log(theta[i,k]) +log(theta[i,k])<-a0+g[k]+u[i]+v[i] +# + } +v[i]~dnorm(0,tauv) +} +g[1]~dnorm(0,taug) +for (k in 2:T){ + g[k]~dnorm(g[k-1],taug) +} +u[1:m]~dcar_normal(adj[1:L],weights[1:L],num[1:m],tauu) +for(k in 1:sumNumNeigh) {weights[k] <- 1} +a0~dnorm(0,tau0) +tauu<-1/pow(sdu,2) +sdu~dunif(0,10) +taug<-1/pow(sdg,2) +sdg~dunif(0,10) +tauv<-1/pow(sdv,2) +sdv~dunif(0,10) +tau0~dgamma(2,0.5) +}) +LSTinits<-list(a0=0.1,sdv=0.1,tau0=0.1,sdg=0.1,sdu=0.1,g=c(0,0,0,0,0,0,0,0,0,0), +v=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0), +u=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, +0,0,0,0,0,0,0,0,0)) + + + + +LST<-nimbleModel(code=STCode,name="LST",constants=LSTConsts, +data=LSTdata,inits=LSTinits) +##################################################################### +######################## Compile ###################################### +CLST<-compileNimble(LST) +#CLST$v +################## Running with monitors and CODA output ##### +LSTConf<-configureMCMC(LST,monitors=CLST$getNodeNames(stochOnly=TRUE,includeData=FALSE),print=TRUE) +##LSTConf$addMonitors(c("a0","u","v","theta")) +LSTMCMC<-buildMCMC(LSTConf) +CLSTMCMC<-compileNimble(LSTMCMC,project=LST) +niter<-20000 +set.seed(1) +samples<-runMCMC(CLSTMCMC,niter=niter,nburnin=15000,nchains=1,summary=TRUE,WAIC=TRUE,samplesAsCodaMCMC=TRUE) +#samples$summary +samples$WAIC \ No newline at end of file