Sources

Ce cours est construit à partir du livre de Rue and Held (2005), Blangiardo et al. (2013), Krainski et al. (2018), Gómez-Rubio (2020) et de l’article Rue, Martino, and Chopin (2009)

Rappel modèle hiérarchique

Exemple illustratif - Mort subite

On s’intéresse aux données de mort subite du nourrisson en Caroline du Nord en 1974 (exemple extrait du livre de Cressie (2015)). Les données sont disponibles dans le package sf

Les données sont ici des comptages par counties, mais on peut imaginer des prix au m2 par IRIS, des taux de chomage par IRIS etc ….

Exemple illustratif - les pluies en Suisse

Rappel estimation bayésienne

Considérons un modèle de régression simple

Rappel estimation bayésienne

Considérons un modèle de régression simple, écrit dans le formalisme du modèle linéaire généralisé

Modèle

\[Y_i \overset{i.i.d} \sim \mathcal{N}(\mu_i, \sigma^2)\]

avec

\[g(\mu_i) = \eta_i = \beta_0 + \beta_1 x_i.\]

Prior

\[ \beta \sim \mathcal{N}(0, 10^6 I_2) \] \[ \log(\sigma^2) \sim logGamma(1, 10^{-5})\]

Posterior

On est surtout intéressé par

\[\pi(\beta\vert \boldsymbol{y}) = \int_{\sigma^2} \pi(\beta, \sigma^2\vert \boldsymbol{y}) d\sigma^2\]

Rappel estimation bayésienne

Les “rares” bons cas

On a une forme conjuguée et on connait \(\pi(\beta, \sigma^2\vert \boldsymbol{y})\)

Généralement

On approche \(\pi(\theta\vert \boldsymbol{y})\) par des méhodes de Monte Carlo car

\[\pi(\theta\vert \boldsymbol{y}) = \frac{\pi(\boldsymbol{y}, \theta)}{\pi(\boldsymbol{y})} =\frac{\pi(\boldsymbol{y}\vert \theta) \pi(\theta)}{\pi(\boldsymbol{y})}\propto \pi(\boldsymbol{y}\vert \theta) \pi(\theta)\]

Problème

Peut être long car

• les chaînes bougent mal
• ca converge mal

C’est souvent le cas en spatial car la forte dépendance entre les données rend la construction d’algorithme efficace difficile.

Le cadre de modélisation considéré

Un vecteur de \(n\) observations \(\mathbf{y} = (y_1, \dots, y_n)\), (dont certaines peuvent être manquantes)

Ces observations suivent une loi de probabilité de la famille exponentielle (cadre du modèle linéaire généralisé)

La moyenne \(\mu_i\) de \(y_i\) est liée au prédicteur linéaire \(\eta_i\) par une fonction adaptée.

Hypothèse d’indépendance conditionnelle

Les observations sont indépendantes conditionnellement à leur prédicteur linéaire :

\[ \eta_i = \alpha + \sum_{j=1}^{n_\beta} \beta_j z_{ji} + \sum_{k=1}^{n_f} f^{(k)}(u_{ki}) + \epsilon_i, \quad i=1,\dots,n \]

où:

• \(\alpha\) est l’intercept,
• \(\beta_j\) sont les coefficients des covariables \(z_j\),
• \(f^{(k)}\) représente les effets aléatoires,
• \(\epsilon_i\) est un terme d’erreur.

Variables latentes en INLA

Le vecteur des effets latents est :

\[ \mathbf{x} = (\eta_1, \dots, \eta_n, \alpha, \beta_1, \dots) \]

Ces effets suivent un champ aléatoire de Markov Gaussien (GMRF). Les observations sont indépendantes conditionnellement aux variables latents et aux hyperparamètres \(\theta_1\).

Implication concrète

• Les priors sur \(\beta\) doivent être gaussiens,
• Les effets aléatoires doivent être gaussiens

Remarque : un paramètre en bayésien peut être vu comme une variable latente

Matrice de Précision

\[\boldsymbol{x}\sim\mathcal{N}(0, Q(\theta_2)^{-1})\]

Le GMRF a une moyenne nulle et une matrice de précision \(Q(\theta_2)\) dépendant des hyperparamètres \(\theta_2\).

\(\theta = (\theta_1, \theta_2)\).

La matrice \(Q(\theta)\) est creuse car

• une indépendance conditionnellent entre \(x_i\) et \(x_j\) conditionnellement à \(x_{-ij}\) se traduit par un 0 dans \(Q(\theta)\).

• Cette sparsité permet d’accélérer les calculs en INLA.

Estimation

Distribution postérieure jointe

La distribution postérieure des variables latentes s’écrit :

\[ \pi(\mathbf{x}, \theta \mid \mathbf{y}) \propto \pi(\mathbf{y} \mid \mathbf{x}, \theta) \pi(\mathbf{x}, \theta) \]

Vraisemblance

Les observations étant indépendantes conditionnellement à \(\mathbf{x}\) et \(\theta\), la vraisemblance devient :

\[ \pi(\mathbf{y} \mid \mathbf{x}, \theta) = \prod_{i \in I} \pi(y_i \mid x_i, \theta) \]

avec \(I\) l’ensemble des indices des observations effectivement observées.

Distribution des variables latentes

On factorise la distribution jointe :

\[ \pi(\mathbf{x}, \theta) = \pi(\mathbf{x} \mid \theta) \pi(\theta) \]

avec \(\pi(\theta)\) une distribution a priori, souvent factorisée en lois univariées.

Distribution postérieure des variables latentes

Comme \(\mathbf{x}\) est un GMRF :

\[ \pi(\mathbf{x} \mid \theta) \propto |Q(\theta)|^{1/2} \exp \left(-\frac{1}{2} \mathbf{x}^T Q(\theta) \mathbf{x} \right) \]

La distribution postérieure devient :

\[ \pi(\mathbf{x}, \theta \mid \mathbf{y}) \propto \pi(\theta) |Q(\theta)|^{1/2} \exp \left(-\frac{1}{2} \mathbf{x}^T Q(\theta) \mathbf{x} + \sum_{i \in I} \log \pi(y_i \mid x_i, \theta) \right) \]

Marginales Postérieures pour les paramètres (hyperparamètres) \(\theta\)

INLA ne cherche pas la postérieure complète mais les marginales des variables latentes et hyperparamètres :

On va vouloir calculer

\[ \pi(x_l \mid \mathbf{y}) = \int \pi(x_l \mid \theta, \mathbf{y}) \pi(\theta \mid \mathbf{y}) d\theta \] et

\[ \pi(\theta_k \mid \mathbf{y}) = \int \pi(\theta \mid \mathbf{y}) d\theta_{-k} \]

avec \(\theta_{-k}\) le vecteur \(\theta\) sans \(\theta_k\).

Approximation en INLA

L’approximation \(\tilde{\pi}(\theta \mid \mathbf{y})\) de Rue, Martino et Chopin (2009) est :

\[ \tilde{\pi}(\theta \mid \mathbf{y}) \propto \frac{\pi(\mathbf{x}, \theta, \mathbf{y})}{ \tilde{\pi}_G(\mathbf{x} \mid \theta, \mathbf{y})} \bigg|_{\mathbf{x} = \mathbf{x}^*(\theta)} \]

avec \(\tilde{\pi}_G(\mathbf{x} \mid \theta, \mathbf{y})\) une approximation gaussienne et \(\mathbf{x}^*(\theta)\) le mode de la conditionnelle.

Approximations Finales

La marginale \(\pi(\theta_k \mid \mathbf{y})\) est obtenue par intégration de \(\tilde{\pi}(\theta \mid \mathbf{y})\).

L’approximation de \(\pi(x_i \mid \mathbf{y})\) utilise l’intégration numérique, avec la méthode de Laplace donnant de meilleurs résultats.

La flexibilité de cette approche

• On rendre compte d’effets non linéaires de covariables en modélisant l’effet de la covariabble comme un processus gaussien (dont il faudra préciser la forme)

• Les temps de calcul son imbattables : quelques secondes pour un champs spatial

Et le spatial dans tout ca ?

INLA et modèles spatiaux

• Le terme d’effets aléatoires peut rendre compte des aspects spatiaux.

Les limites

Si le nombre de données est important, le vecteur \(\boldsymbol{x}\) devient trop grand et lestemps de calcul augmentent rapidement

Limiter l’impact du nombre de points : l’approche SPDE

Le nombre de points important en géostatistique.

Il se trouve qu’un champs gaussien de Matern peut s’écrire comme la solution d’une équation aux dérivées partielles stochastique

Processus Gaussiens avec Covariance de Matérn

Un champ gaussien \(U(s)\) avec covariance de Matérn est solution de l’EDPS fractionnaire suivante :

\[ (\kappa^2 - \Delta)^{\alpha/2} u(s) = W(s), \quad s \in \mathbb{R}^d, \quad \alpha = \nu + d/2, \quad \kappa > 0, \quad \nu > 0. \]

où : - \(\Delta\) est l’opérateur de Laplace, - \(W(s)\) est un bruit blanc gaussien spatial de variance unité.

Processus Gaussiens avec Covariance de Matérn

Fonction de Corrélation de Matérn

Rappel est définie par un paramètre d’échelle \(\kappa > 0\) et un paramètre de régularité \(\nu > 0\). Pour deux positions \(s_i\) et \(s_j\), la fonction de corrélation stationnaire et isotrope de Matérn est donnée par :

\[ \text{Cor}_M(U(s_i), U(s_j)) = \frac{2^{1-\nu}}{\Gamma(\nu)} (\kappa \|s_i - s_j\|)^\nu K_\nu (\kappa \|s_i - s_j\|) \]

où \(K_\nu\) est la fonction de Bessel modifiée de seconde espèce et \(\Gamma(\nu)\) est la fonction gamma.

Un champ gaussien \(U(s)\) avec covariance de Matérn est solution de l’EDPS fractionnaire suivante :

\[ (\kappa^2 - \Delta)^{\alpha/2} u(s) = W(s), \quad s \in \mathbb{R}^d, \quad \alpha = \nu + d/2, \quad \kappa > 0, \quad \nu > 0. \]

où : - \(\Delta\) est l’opérateur de Laplace, - \(W(s)\) est un bruit blanc gaussien spatial de variance unité.

Résolution de l’EDPS sur un maillage

C’est une technique classique des numériciens on se donne un maillage et on résoud l’EDPS sur les points du maillage. On interpole ensuite au sein du maillage.

La finesse du maillage contrôle la complexité du problème et sa précision

Conclusion

• On dispose d’un outil rapide pour des modèles à variable latente y compris spatiaux (et même spatio temporels)

• des contraintes dans la modélisation

• des habitudes de modélisation à acquérir.