9.4 Mapping exceedance probabilities

Calculate the exceedance probabilities of malaria prevalence being greater than a given threshold value.

\[P(p_1>c)=1-P(p_i\leq c)\]

Once the the marginals are found the probability that malaria prevalence exceeds 20% at this location is 0.6931.

How to calculate the marginals?

The marginal distribution of the predictions listed in the res$marginals.fitted.values[index] are indices of the stk.full:

index <- inla.stack.index(stack = stk.full, 
                          tag = "pred")$data

res2 <- 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)),
  control.compute=list(return.marginals.predictor=TRUE))

summary(res2)

## 
## 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.683, Running = 1.3, Post = 0.0703, Total = 2.05 
## 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)')

# res2$marginals.fitted.values[index]

marg <- res2$marginals.fitted.values[index][[1]]

1 - inla.pmarginal(q = 0.20, marginal = marg)

excprob <- sapply(res2$marginals.fitted.values[index],
                  FUN = function(marg){
                    1-inla.pmarginal(q = 0.20, 
                                     marginal = marg)})

head(excprob)

## fitted.APredictor.066 fitted.APredictor.067 fitted.APredictor.068 
##             0.7150007             0.7435313             0.7926329 
## fitted.APredictor.069 fitted.APredictor.070 fitted.APredictor.071 
##             0.7669259             0.6861877             0.7221892

r_excprob <- rasterize(
  x = coop, y = ra, 
  field = excprob,
  fun = mean)

pal <- colorNumeric("viridis", c(0, 1), na.color = "transparent")

leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addRasterImage(r_excprob, colors = pal, opacity = 0.5) %>%
  addLegend("bottomright",
    pal = pal,
    values = values(r_excprob), title = "P(p>0.2)"
  ) %>%
  addScaleBar(position = c("bottomleft"))

END