6.3 Modeling

6.3.1 Our model:

First assumptions: Observed count is Poisson distributed (\(y_i\))

\[Y_i \sim Poisson(E_i\theta_i), i = 1, ..., n\]

With \(E_i\) is the expected count and \(\theta_i\) is the relative risk in area \(i\).

\[log(\theta_i) = \beta_0 + \beta_1 \times AFF_i + u_i + v_i\]

\[\mu | \mu_{-i} \sim N (\bar{u_{\delta_i}}, \frac{\sigma^2_u}{n_{\delta_i}})\]

\[v_i \sim N(0, \sigma²_v)\]

6.3.2 neighborhood matrix

library(spdep)
library(INLA)
nb <- poly2nb(map) # same result with map_3
nb2INLA("map.adj", nb)
g <- inla.read.graph(filename = "map.adj")

6.3.3 Using INLA

With have two random effects \(u\) and \(v\):

map$idareau <- 1:nrow(map@data) # u
map$idareav <- 1:nrow(map@data) # v
# sf map_3iareau <- 1:nrow(map_3)

Same model than chapter 5:

formula <- Y ~ AFF +
  f(idareau, model = "besag", graph = g, scale.model = TRUE) +
  f(idareav, model = "iid")

res <- inla(formula,
  family = "poisson", data = map@data,
  E = E, control.predictor = list(compute = TRUE),
  # this need to be added for latter!!!
  #control.compute=list(return.marginals.predictor=TRUE)

)
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.373, Running = 0.188, Post = 0.0238, Total = 0.585 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.306 0.120     -0.538   -0.306     -0.068 -0.308   0
## AFF          4.310 1.276      1.742    4.331      6.762  4.371   0
## 
## Random effects:
##   Name     Model
##     idareau Besags ICAR model
##    idareav IID model
## 
## Model hyperparameters:
##                           mean       sd 0.025quant 0.5quant 0.975quant    mode
## Precision for idareau     4.17     1.43       2.08     3.94       7.63    3.50
## Precision for idareav 18679.74 18489.30    1265.39 13102.10   67966.71 3476.67
## 
## Marginal log-Likelihood:  -189.63 
##  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)')

Greatly surprise by:

this was quick!
AFF increases lip cancer risk (and it is clear)

library(ggplot2)
marginal <- inla.smarginal(res$marginals.fixed$AFF) # list
marginal <- data.frame(marginal)
ggplot(marginal, aes(x = x, y = y)) + geom_line() +
  labs(x = expression(beta[1]), y = "Density") +
  geom_vline(xintercept = 0, col = "black") + theme_bw()