# Base de R INLA
library(tidyverse)
library(INLA)
library(sf)
library(gstat)TD4 - Krigeage avec INLA
$$
$$
Etude des précipitations à Parana au Brésil
Comme dans l’exemple des pluies suisses nous allons chercher à prédire les précipitation à Parana au Brésil.
Présentation des données
Les données sont disponibles dan ls le pakage `geoR et contiennent
parana$coords les coordonnées en km dans un système arbitraire,
parana$borders les contours de l’état dans le même système,
parana$data les données de précipitation
Nous allons mettre ces données au standar sf et créer une grille de prédiction.
# load data
data(parana, package = "geoR")
# locations data and reprsentation
locations <-as.data.frame(parana$coords) |> st_as_sf(coords = c('east', 'north'))
locations |> ggplot() + geom_sf()
# rainfall data
pluies <- parana$data |> as.data.frame() |> rename(precip = `parana$data`)
# join rainfall and locations data
pluies_sf <- pluies |> bind_cols(locations) |> st_as_sf()
pluies_sf |> ggplot() + geom_sf(aes(geometry = geometry, col = precip)) +
scale_colour_gradient(low = "#ffd433",high = "#052265")
# parana state limit
border_sf <- st_polygon(list(parana$borders))
pluies_sf |> ggplot() + geom_sf(aes(geometry = geometry, col = precip)) +
scale_colour_gradient(low = "#ffd433",high = "#052265") + geom_sf(data = border_sf, fill = NA)
# define prediction grid
mycell <- 20 #km2
my_grid <- st_make_grid(border_sf, cellsize = mycell) |> st_as_sf() |> mutate(center = st_centroid(x))
pluies_sf |> ggplot() + geom_sf(aes(geometry = geometry, col = precip)) +
scale_colour_gradient(low = "#ffd433",high = "#052265") + geom_sf(data = border_sf, fill = NA) + geom_sf(data = my_grid, fill = NA, aes(geometry = x)) +geom_sf(data = my_grid, aes(geometry = center), col = "red") On souhaite prédire la pluviométrie en tout point avec INLA
On est dans le cadre de données géostatistiques, on va donc procéder par krigeage. Avant de kriger il faut estimer les paramètres du champs spatial. On l’avait par ajustement du variaogramme dans la version géostatistique classique.
- Ajuster un variogramme exponentiel sur les données de pluies. Discuter l’ajustement.
Il serait possible de rajouter un effet pépite pour tenir compte d’une variabilité à très faible échelle. On la le faire dans un cadre hiérarchique.
On note \(Y_i\) la précipitation au point \(s_i\)
\[Y_i \overset{ind}{\sim}\mathcal{N}(\mu_i, \sigma^2)\]
\(\$mu_i = \beta_0 + Z(s_i).\)$
On suppose que \(Z(s)\) est un champs gaussien de fonction de covariance de Matern pour pouvoir appliquer l’approche INLA-SPDE.
Quelques rappels
Matèrn model
\[Cov(Z(s_i), Z(s_j)) = \frac{\sigma^2_Z}{2^{v-1}\Gamma(v)}(k||s_i - s_j||)^vK_v(k||s_i-s_j||)\]
Ici, \(\sigma^2_z\) représente la variance marginale du champ spatial, et \(K_{\nu}(\cdot)\) est la fonction de Bessel modifiée de deuxième espèce et d’ordre $> 0 $.
La valeur de \(\nu\) détermine la différentiabilité au sens des moindres carrés du processus et est généralement fixée, car elle est mal identifiée dans les applications.
Pour \(\nu = \frac{1}{2}\), la fonction de covariance de Matérn est équivalente à la fonction de covariance exponentielle.
Le paramètre \(\kappa > 0\) est lié au portée spatiale \(\rho\), c’est-à-dire la distance à laquelle la corrélation entre deux points est approximativement nulle. Plus précisément,
\(\rho = \frac{\sqrt{8\nu}}{\kappa}\) et à cette distance, la corrélation spatiale est proche de 0.1.
SPDE et Covariance de Matern
Un champ de markov gaussien (GRF) \(z(s)\) sur un domaine \(D \subset {\mathbb{R}}^2\) avec une matrice de covariance de Matérn peut être exprimé comme la solution de l’équation différentielle stochastique aux dérivées partielles (SPDE) suivante dans un domaine continu :
\[ (\kappa^2 - \Delta)^{\alpha/2} (\tau z(s)) = W(s). \]
Ici \(W(s)\) est un bruit blanc spatial gaussien. Le paramètre $$contrôle la régularité du champ \(z(s)\) , $$contrôle la variance, et $> 0 \(est un **paramètre d'échelle**.\)$ représente le laplacien
Les paramètres de la fonction de covariance de Matérn et ceux de l’équation SPDE sont couplés de la manière suivante :
Le paramètre de régularité \(\nu\) de la fonction de covariance de Matérn est lié au paramètre \(\alpha\) de l’équation SPDE par la relation (valable sur \({\mathbb{R}}^2\)):
\[ \nu = \alpha - 1. \]
La variance marginale \(\sigma^2_Z\) est reliée aux paramètres de l’équation SPDE par :
\[ \sigma^2_Z = \frac{\Gamma(\nu)}{\Gamma(\alpha)} \frac{1}{(4\pi)^{d/2}} \frac{1}{\kappa^{2\nu} \tau^2}. \]
Pour \(\nu = \frac{1}{2}\) (ce qui correspond à la fonction de covariance exponentielle), on obtient \(\alpha = 3/2\). Dans le package R-INLA, la valeur par défaut est \(\alpha = 2\) , bien que des options dans l’intervalle \(0 \leq \alpha < 2\) soient également disponibles.
- Construire un premier maillage sans tenir compte des points de localisation. On peut choisir 50 et 150 comme distance au noeud le plus proche.
- Construire un second maillage plus raffiné aux endroits où des données sont disponibles
- Pour représenter la maillage, il est possible d’utiliser les facilités du package inlabru
Quel est le nombre de noeuds du maillage
Preparation inla
6.1 Définir l’équation SPDE à partir du maillage avec un paramètre \(alpha = 2\) et en imposant la contrainte que la moyenne du champs soit nulle
6.2 Créer les données pour un effet spatial nommé s dans inla (ìnla.spde.make.index`)
6.3 Créer la matrice de projection qui fait le lien entre les données et les points du maillage. Quelles sont les dimensions de cette matrice, à quoi correspondent elles ?
6.4 Construire la pile (stack) inla pour ajuster le modèle
6.5 Appeler inla pour ajuster un modèle spatial comme spécifier au début du TP.
L’optimisation numérique ne se passe pas très bien. pour aider l’estimation on va modifier le prior, pour contraindre la variance d’observation dans des gammes de valeurs raisonnables. La variogramme initial laisse supposer une pépite d’environ 600, donc on peut s’attendre à une estimation de la précision de l’observation d’environ \(1/600\approx0.0015\). On va spécifier une loi gamma sur la précision de moyenne 1 et variance 1 pour contraindre le problème d’optimisation. on peut aussi donner notre estimation de la prcision comme valeur initiale.
Il suffit d’ajouter control.family = list(hyper = list(prec = list(prior = "loggamma", param = c(1, 1), initial = 0.0015 ))) dans l’appel à INLA.
6.6 Predire sur la grille. Ajouter sur la pile inla la grille comme jeu de données de prédiction.