9.2 Modeling

Steps to fit the model using the SPDE approach and the R-INLA package:

\[Y_i|P(x_i)\sim\text{Binomial}(N_i,P(x_i))\] \[logit(P(x_i))=\beta_0+\beta_1 \text{x altitude + }S(x_i)\]

\[Cov(S(x_i),S(x_j))\] \(K_v(.)\) is the modified Bessel function. \(v\) is the smoother parameter. \(\rho=\sqrt{8v}/k\) is the distance where the spatial correlation is close to 0.1.

9.2.1 Mesh construction

Build a triangulated mesh that covers The Gambia:

library(INLA)
coo <- cbind(d2$long, d2$lat)


mesh <- inla.mesh.2d(
  loc = coo, 
  max.edge = c(0.1, 5),
  cutoff = 0.01)

Number of vertices:

mesh$n

## [1] 674

plot(mesh)
points(coo, col = "red")

9.2.2 Building the SPDE model on the mesh

The smoother parameter for the stochastic partial differential equation SPDE is \(\alpha=v+d/2\) defined by \(v=1\) and \(d=2\):

?inla.spde2.matern

spde <- inla.spde2.matern(mesh = mesh, 
                          alpha = 2, 
                          constr = TRUE)

9.2.3 Index set

Index set for the SPDE model where \(s\) is the effect:

indexs <- inla.spde.make.index("s", spde$n.spde)
lengths(indexs)

##       s s.group  s.repl 
##     674     674     674

9.2.4 Projection matrix

A <- inla.spde.make.A(mesh = mesh, loc = coo)

9.2.5 Prediction data

Specify the locations where to predict the prevalence:

dp <- raster::rasterToPoints(r)
dim(dp)

## [1] 12964     3

ra <- raster::aggregate(r, fact = 5, fun = mean)

dp <- raster::rasterToPoints(ra)
dim(dp)

## [1] 627   3

coop <- dp[, c("x", "y")]

Ap <- inla.spde.make.A(mesh = mesh, loc = coop)

9.2.6 Stack with data for estimation and prediction

Construct a stack for estimation called stk.e and a stack for prediction stk.p, then put them together in a full stack stk.full:

# stack for estimation stk.e
stk.e <- inla.stack(
  tag = "est",
  data = list(y = d2$positive, 
              numtrials = d2$total),
  A = list(1, A),
  effects = list(data.frame(b0 = 1, 
                            altitude = d2$alt), 
                 s = indexs))

# stack for prediction stk.p
stk.p <- inla.stack(
  tag = "pred",
  data = list(y = NA, 
              numtrials = NA),
  A = list(1, Ap),
  effects = list(data.frame(b0 = 1, 
                            altitude = dp[, 3]),
                 s = indexs))

# stk.full has stk.e and stk.p
stk.full <- inla.stack(stk.e, stk.p)

9.2.7 Model formula

formula <- y ~ 0 + b0 + altitude + f(s, model = spde)

?inla()
?control.predictor

res <- inla(formula,
  family = "binomial", 
  Ntrials = numtrials,
  control.family = list(link = "logit"),
  data = inla.stack.data(stk.full),
  control.predictor = list(compute = TRUE, 
                           link = 1,
                           A = inla.stack.A(stk.full)))

summary(res)

## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, 
##    ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
##    scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
##    ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
##    verbose, ", " lincomb = lincomb, selection = selection, control.compute 
##    = control.compute, ", " control.predictor = control.predictor, 
##    control.family = control.family, ", " control.inla = control.inla, 
##    control.fixed = control.fixed, ", " control.mode = control.mode, 
##    control.expert = control.expert, ", " control.hazard = control.hazard, 
##    control.lincomb = control.lincomb, ", " control.update = 
##    control.update, control.lp.scale = control.lp.scale, ", " 
##    control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
##    ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
##    num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, 
##    working.directory = working.directory, ", " silent = silent, inla.mode 
##    = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = 
##    .parent.frame)") 
## Time used:
##     Pre = 0.674, Running = 1.29, Post = 0.035, Total = 1.99 
## Fixed effects:
##            mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## b0       -0.199 0.453     -1.001   -0.231      0.806 -0.279   0
## altitude -0.012 0.011     -0.032   -0.012      0.010 -0.012   0
## 
## Random effects:
##   Name     Model
##     s SPDE2 model
## 
## Model hyperparameters:
##               mean    sd 0.025quant 0.5quant 0.975quant  mode
## Theta1 for s -4.03 0.297      -4.61    -4.03      -3.44 -4.03
## Theta2 for s  2.67 0.377       1.93     2.67       3.41  2.68
## 
## Marginal log-Likelihood:  -216.61 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')