TD3 - Ajuster un modèle hiérarchique spatial avec INLA

Author
Affiliations

Marie-Pierre Etienne (sur la base des idées expolorée dans FinistR2023 )

ENSAI - CREST

https://marieetienne.github.io/statspat

Published

March 11, 2025

# Base de R INLA
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.4     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(INLA)
Loading required package: Matrix

Attaching package: 'Matrix'

The following objects are masked from 'package:tidyr':

    expand, pack, unpack

This is INLA_24.12.11 built 2024-12-11 19:58:26 UTC.
 - See www.r-inla.org/contact-us for how to get help.
 - List available models/likelihoods/etc with inla.list.models()
 - Use inla.doc(<NAME>) to access documentation
library(sf)
Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE

Etude du nombre de cas de malaria en Gambie

Nous nous intéressons à un jeu de données concernant la prévalence de la malaria en Gambie (disponible dans le package {geoR}). Cet exemple est repris du livre “Spatial and Spatio-temporal Bayesian Models with R-INLA”.

data(gambia, package = "geoR")
# les coordonnées correspondent au village où se trouve les enfants
gambia <- gambia |> 
  mutate(village_index = factor(as.numeric(factor(paste(x,y, sep="_")))))

Le jeu de données contient pour chaque personne testée

  • les coordonées du village (x, y)
  • le résultat du test (pos),
  • l’âge en jours, *l’utilisation d’une moustiquaire netuse,
  • l’administration éventuel d’un traitement sur la moustiquaire treated,
  • une mesure par satellite de la végétation proche green,
  • la présence ou non d’un centre de santé dans le village.
  1. Calculer la prévalence (proportion de cas positif) dans chaque village
# A tibble: 65 × 4
   village_index total   pos  prev
   <fct>         <int> <dbl> <dbl>
 1 1                33    17 0.515
 2 2                63    19 0.302
 3 3                17     7 0.412
 4 4                24     8 0.333
 5 5                26    10 0.385
 6 6                18     7 0.389
 7 7                38    24 0.632
 8 8                56     7 0.125
 9 9                57     6 0.105
10 10               26     8 0.308
# ℹ 55 more rows

On transforme le jeu de données en données de type sf.

  1. Indiquer ce que signifie à votre avis crs = "+proj=utm +zone=28" dans la commande suivante
gambia_sf <- st_as_sf(as.data.frame(gambia), coords = c("x", "y"),
               crs = "+proj=utm +zone=28")
gambia_sf <- gambia_sf |> st_transform( crs = 4326)
gambia_sf |> ggplot() + geom_sf()

# 
# # Transformer les données de bord de Gambia en objet sf
# gambia.borders$group <- cumsum(is.na(gambia.borders$x))
# 
# gambia.borders <- gambia.borders |> 
#   filter(!is.na(x))
# polygons <- lapply(unique(gambia.borders$group), 
#                    function(a){
#                      df <- gambia.borders |> 
#                        filter( group == a) 
#                      df |> 
#                        bind_rows(slice(df ,1)) |> 
#                        summarise(geometry = st_sfc(st_polygon(list(cbind(x, y)))))  |> 
#                        st_sf() })
# gambia.borders_sf <- do.call('rbind', polygons)
# st_crs(gambia.borders_sf) <- 32628  #  UTM Zone 28N
# 
# gambia.borders_sf <- gambia.borders_sf    |> 
#  st_transform( crs = 4326)
# 
#  ggplot() + geom_sf(data=gambia.borders_sf) + geom_sf(data=gambia_sf)

Dans les approches INLA, le champs gaussien caché est vu comme la solution d’une Equation aux Dérivées Partielles Stochastiques (EDPS) qui va être résolue sur un maillage. On peut définir le maillage.

hull = inla.nonconvex.hull(gambia_sf,convex = -0.1)
gambia_mesh <- inla.mesh.2d(boundary = hull, max.edge= 0.2)

plot(gambia_mesh,main="",asp=1)
coo <- st_coordinates(gambia_sf)
points(coo, col = "red")

Mais pour résoudre l’EDPS, on a besoin de mettre des conditions de nullité au bord du domaine ce qui a des conséquences forte sur la solution. Pour limiter l’impact on ajoute buffer autour avec une maille plus large.

hull = inla.nonconvex.hull(gambia_sf,convex = -0.1)
gambia_mesh <- inla.mesh.2d( boundary = hull,max.edge = c(0.2, 2), offset = c(1, 0.2))
plot(gambia_mesh,main="",asp=1)
  1. Changer la valeur de max.edge. Que se passe t il ?

En général, l’utilisation des paramètres du maillage est assez complexe.

  1. Il est naturel d’espèrer un résultat plus fin aux endroits où l’on a des données, on peut lui passer les localisations du jeu de données pour construire le maillage grâce à l’argument loc. Que constatez vous ?
gambia_mesh <- inla.mesh.2d( boundary = hull,loc = gambia_sf, max.edge = c(0.2, 2) )
plot(gambia_mesh,main="",asp=1)
coo <- st_coordinates(gambia_sf)
points(coo, col = "red")

On définit à partir du maillage le champ spatial SPDE qui correspond à un champ spatial gaussien avec une covariance Matérn. Nous considérons les lois a priori par défaut sur les paramètres de variance et de portée du champ spatial.

gambia_spde <- inla.spde2.matern(mesh = gambia_mesh, constr = TRUE)

L’argument constr = TRUE impose une moyenne nulle pour le champs spatial.

On souhaite ajuster le modèle suivant

Pour l’enfant \(j\) du village \(i\), on suppose que conditionnellement aux covariables et à un champs spatial caché, l’infection à la malaria est le résultat d’un tirage de bernoulli.

  1. Ecrire le modèle correspondant sous forme mathématique.

La mise en oeuvre dans INLA est ensuite un peu technique (le package inlabru permettrait de simplifier les étapes suivantes)

  1. Il faut créer les éléments dont INLA a besoin

6.a Un Index set

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

6.b Une matrice de projection. Cette matrice fait le lien entre les positions des observations et les noeuds du maillage. C’est une matrice très creuse.

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

6.c On doit ensuite empiler ces données dans une structure spécifique à INLA

gambia_effects <- gambia_sf |> 
  select(-pos, -village_index) |> st_drop_geometry() |> mutate(b0=1)

stk.e <- inla.stack(
  tag = "est",
  data = list( pos = gambia_sf$pos), A = list(1, A),
  effects = list(gambia_effects, s= indexs)
)
  1. On peut appeler INLA. Explorer les sorties contenues dans l’objet fit
formula = pos ~ 0 + b0 + treated +
  netuse +
  age +
  green +
  phc +
   f(s, model = gambia_spde)

fit <- inla(formula,
  family = "binomial",
  control.family = list(link = "logit"),
  data = inla.stack.data(stk.e),
  control.predictor = list(compute = TRUE, link = 1,
                           A = inla.stack.A(stk.e)),
  control.compute = list(return.marginals.predictor = TRUE)
)
summary(fit)
  1. Accéder aux distributions marginales des effets aléatoires et des hyperparamètres :
fit$summary.random 
fit$summary.hyperpar
  1. Tacer les distributions a posteriori marginales des effets, par exemple pour l’effet de l’âge ou de l’utilisation de la moustiquaire :
  1. Représenter la composante spatiale de la malaria qui n’est pas expliquée par les covariable. Pour ceci, on peut créer un jeu de données de prédiction
gambia_pred_sf <- gambia_sf |> group_by(village_index) |> 
  summarise(age_moyen = mean(age), netuse = 0,  treated = 0, phc = mean(phc), green = mean(green))
  
coo.pred <- st_coordinates(gambia_pred_sf)
gambia_pred <- gambia_pred_sf |> st_drop_geometry() |> mutate(b0= 1)
A.pred <- inla.spde.make.A(mesh = gambia_mesh, loc = coo.pred)
stk.p <- inla.stack(
  tag = "pred",
  data = list(pos = NA),
  A = list(1, A.pred),
  effects = list(gambia_pred,
    s = indexs
  )
)

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

fit <- inla(formula,
  family = "binomial",
  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(fit)

index <- inla.stack.index(stack = stk.full, tag = "pred")$data
prev_mean <- res$summary.fitted.values[index, "mean"]


gambia_pred_sf |> mutate(prev_mean_nonet= prev_mean) |>  
ggplot() + geom_sf( aes(col = prev_mean_nonet ))
  1. le package geodata contient les altitudes en Gambie. Reprenez les étapes pour ajuster un modèle dans lequel on explique la prévalence de la malaria par l’altitude. On pourra alors grouper les données par village.
library(geodata)
library(terra)

r <- elevation_30s(country = "GMB", path = tempdir()) 
gambia[, c("long", "lat")] <- st_coordinates(gambia_sf)
gambia$elev <- terra::extract(r, gambia[, c("long", "lat")])$GMB_elv_msk

References