-
Notifications
You must be signed in to change notification settings - Fork 15
/
SCcongen_INLA_models_plusBYM.txt
77 lines (53 loc) · 3.09 KB
/
SCcongen_INLA_models_plusBYM.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
##############################################
setwd("working directory")
##############################################
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))
region<-seq(1:46);region2<-region
library(INLA)
#UH model
formulaUH = obs~ f(region, model = "iid",param=c(2,1))
resultUH = inla(formulaUH,family="poisson",
data=SCcongen90,control.compute=list(dic=TRUE,cpo=TRUE),E=expe)
summary(resultUH)
#UH + poverty covariate
formulaUHPov = obs~ 1+pov+f(region, model = "iid",param=c(2,1))
resultUHPov = inla(formulaUHPov,family="poisson",
data=SCcongen90,control.compute=list(dic=TRUE,cpo=TRUE),E=expe)
summary(resultUHPov)
#UH+CH + poverty covariate
formulaUHCHPov = obs~ 1+pov+f(region, model = "iid",param=c(2,1))+f(region2,model="besag",graph="SC.poly.txt",param=c(2,1))
resultUHCHPov = inla(formulaUHCHPov,family="poisson",
data=SCcongen90,control.compute=list(dic=TRUE,cpo=TRUE),E=expe)
summary(resultUHCHPov)
#BYM+covariate
formulaUHCHPov2 = obs~ 1+pov+f(region2,model="bym",graph="SC.poly.txt",param=c(2,1))
resultUHCHPov2 = inla(formulaUHCHPov,family="poisson",
data=SCcongen90,control.compute=list(dic=TRUE,cpo=TRUE),E=expe)
#############################################
#### graphics ###############################
#############################################
#graphics setup#
setwd("C:/A_NEW_FOLDER/A_NEW_FOLDER/A_COURSES/A_courses/Bayesian Disease Mapping BDM courses/BDM_NEW_COURSEs_2008-present/NEW COURSE 2015/MUSC/BDMI participant files/BDMI participant files.zip/BDMI participant files")
library(maptools)
polySC<-readShapePoly("SC_county_alphasort.shp")
plot(polySC)
source("fillmap.R")
#### UH only model ###
UH_effect<-resultUH$summary.random$region[,2]
fillmap(polySC,"UH effect",UH_effect*100,n.col=5)
#### UH + Poverty model #######
UH_effect<-resultUHPov$summary.random$region[,2]
fillmap(polySC,"UH effect",UH_effect*100,n.col=5)
#### UH + CH + Poverty model #####
UH_effect<-resultUHCHPov$summary.random$region[,2]
CH_effect<-resultUHCHPov$summary.random$region2[,2]
fillmap(polySC,"UH effect",UH_effect*100,n.col=5)
fillmap(polySC,"CH effect",CH_effect*100,n.col=5)