##load the data set
data(Zambia)
##load map
g = system.file("demodata/zambia.graph", package="INLA")
# add one column for the unstructured spatial effect
Zambia$distr.unstruct = Zambia$district
##define formulas for the three models
##MOD1
formula.mod1 = hazstd ~ agc +
f(district, model="besag", graph.file = g, param = c(1,0.01)) +
f(distr.unstruct, model="iid", param = c(1,0.01)) +
edu1 + edu2 + tpr + sex + bmi
##MOD2
formula.mod2 = hazstd ~ f(agc, model = "rw2") +
f(district, model="besag", graph.file = g, param = c(1,0.01)) +
f(distr.unstruct, model="iid", param = c(1,0.01)) +
edu1 + edu2 + tpr + sex + bmi
##MOD3
formula.mod3 = hazstd ~f(agc, model = "rw2") +
f(district, bmi, model = "besag", graph.file = g, param = c(1,0.01), constr = FALSE) +
edu1 + edu2 + tpr + sex
##run the three models
mod1 = inla(formula.mod1, data = Zambia, control.family = list(initial = 1),
control.inla = list(h=1e-4),
control.compute = list(dic = TRUE, cpo=TRUE),
verbose = TRUE)
mod2 = inla(formula.mod2, data = Zambia, control.family = list(initial = 1),
control.inla = list(h = 1e-4),
control.compute = list(dic = TRUE,cpo=TRUE),
verbose = TRUE)
mod3 = inla(formula.mod3,data = Zambia, control.family = list(initial = 1),
control.inla = list(h = 1e-4),
control.compute = list(dic = TRUE,cpo=TRUE),
verbose = TRUE)
##compute the log-score
log.score1 = -mean(log(mod1$cpo))
log.score2 = -mean(log(mod2$cpo))
log.score3 = -mean(log(mod3$cpo))
##pit histogram
par(mfrow=c(3,1))
hist(mod1$pit,br=30)
hist(mod3$pit,br=30)
hist(mod2$pit,br=30)