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)

Comments