Code for the Zambia Example

##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)