Dans deux premiers articles (Modélisation spatio-temporelle avec une équation intégro-différentielle (IDE) avec QGis et R et Modélisation spatio-temporelle avec une équation intégro-différentielle (IDE) avec QGis et R(suite)) nous avons vu comment calculer et ajuster un modèle intégro-différentiel à un ensemble de données spatio-temporelles, en utilisant R et QGis

Dans un troisième article, Comparaison et sélection de modèles intégro-différentiels (IDE): les principes, nous avons vu comment envisager la validation des modèles calculés et ajustés lors des deux premiers articles, et son corollaire, comment calculer des scores sur chaque version du modèle de manière à pouvoir choisir celui qui prédit le mieux les valeurs observées. Puis, dans l’article précédent, Comparaison et sélection de modèles intégro-différentiels (IDE): mode d’emploi des scripts QGis/R, nous avons vu le mode d’emploi des interfaces des scripts.

Dans cet article, le dernier de la série, nous allons suivre un exemple pratique avec les données des deux premiers articles.

Téléchargement des scripts et données d’exemple.

Si vous n’avez pas téléchargé les scripts, vous pouvez le faire à partir du lien ci-dessous:

(Voir les différences des versions)

Ces scripts sont à copier, avec leur fichier d’aide, dans le répertoire C:\Users\xxxxx\AppData\Roaming\QGIS\QGIS3\profiles\default\processing\rscripts

Rappel des données d’exemple

Les données sont extraites des Dénombrements du Relevé des oiseaux nicheurs (BBS). Ces données proviennent du Relevé des oiseaux nicheurs de l’Amérique du Nord.
Nous considérons en particulier les dénombrements annuels du troglodyte de Caroline (Thryothorus ludovicianus) pour la période de 1967 à 2014. L’unité d’échantillonnage du BBS est un parcours routier d’environ 39,2 km. Dans chaque unité d’échantillonnage, les observateurs bénévoles font 50 arrêts et comptent les oiseaux pendant une période de 3 minutes lorsqu’ils exécutent leurs itinéraires (habituellement en juin). Il y a plus de 4000 routes dans l’enquête nord-américaine, mais toutes ne sont pas disponibles chaque année. Aux fins des analyses présentées ici, nous considérons que les dénombrements des parcours totaux se produisent chaque année (pendant la saison de reproduction) et définissons l’emplacement spatial de chaque parcours comme étant le centroïde de la route. Ainsi, nous considérons que les données sont discrètes dans le temps, géostatistiques et irrégulières dans l’espace, et non gaussiennes dans le sens qu’elles sont des comptages.

Copiez les fichiers shapefile dans un répertoire de données et chargez la couche dans QGis.

Premier cas: modèle non variant sans validation croisée.

Nous chargeons la couche MOcarolinawren dans QGis. Si nous ouvrons la table nous voyons quatre attributs

• un identifiant
• un attribut t correspondant au rang du pas de temps, avec des valeurs comprises entre 1 et 21
• un attribut cnt qui correspond au comptage d’oiseaux
• un attribut year qui indique l’année du comptage, avec des valeurs comprises entre 1994 et 2014

Nous allons construire un modèle IDE des comptages, cnt, avec comme attribut temporel le rang t.

Nous double-cliquons sur le script R -> ide_modele_invariant

Nous renseignons les champs t avec t et val avec cnt

On garde la valeur par défaut de grid_size à 50.

Puisque notre variable temporelle est un rang on indique les unités en « days » même s’il s’agît d’années.

On demande deux périodes de prédiction supplémentaires. Elles auront comme valeur de t, 22 et 23, et elles correspondront aux années 2015 et 2016.

On laisse les valeurs par défaut pour tous les autres paramètres.

Le temps de calcul pour les 100 itérations a été de presque 4 heures et demie. Le journal du script affiche les valeurs calculés pour les différents paramètres du modèle:

Comme nous n’avons pas défini de répertoire de sortie pour les graphiques, nous trouvons le fichier Rplot.pdf dans le dossier « Documents » de l’utilisateur. Nous trouvons dans ce fichier le graphique de la corrélation des valeurs des prédictions en fonction des valeurs des observations:

Et la fenêtre cartographique montre les deux couches résultat: sous forme de grille de 50 x 50 points et sous forme d’observations.

Nous pouvons affecter une symbologie de type remplissage graduée,puis en appliquant un filtre par pas de temps, nous voyons les prédictions du modèle sur l’ensemble de la zone.

Ci-dessous les images représentent les prédictions du modèle pour les observations des pas de temps 20 et 21.

Qu’est(ce qu’on peut dire de ces deux images? Il semble avoir une progression de la population vers le NE et une légère diminution dans le coin SE. Une des caractéristiques de la modélisation IDE est de pouvoir séparer la partie croissance de la population de la partie dispersion des nouveaux individus.

C’est une bonne raison de tenter de modéliser ces observations avec un noyau variant spatialement. Mais avant, il est important de valider le modèle, c’est à dire accumuler des critères objectifs , des scores.

Deuxième cas: modèle non variant avec validation croisée.

Nous gardons les mêmes paramètres que pour le premier cas, mais nous allons cocher la case « effectuer la validation du modèle » et remplir les paramètres de validation.

La validation croisée implique de mettre de côté une partie des données des observations, calculer le modèle sur le reste, puis voir les différences des prédictions entre le modèle calculé et les données que nous avions laissé de côte.

Il faut donc disposer d’assez d’observations pour pouvoir s’en passer d’une partie sans pénaliser l’ajustement du modèle.

D’autre part, la modélisation peut servir pour deux objectifs différents:

• on peut chercher à modéliser le futur, c’est à dire, disposer d’un modèle qui nous permette d’anticiper l’état de notre processus, avant de disposer des observations correspondantes
• on peut chercher aussi à modéliser notre processus pour disposer de données aux endroits où on n’a pas d’observations ‘(« boucher des trous ») ou bien pour rendre continues des observations ponctuelles (interpolation, krigeage)

Selon que l’on ai l’un ou l’autre, ou les deux, objectifs, on utilisera une méthode de validation croisée différente.

Ici nous allons utiliser les deux, mais ce n’est pas obligatoire.

Validation croisée d’un pas de temps

Pour la première méthode de validation, modéliser le futur, la démarche plus appropriée est de prendre le dernier pas de temps (dans notre cas t=21) et de le mettre de côté. Une fois l’ajustement du modèle réalisé avec les pas de temps 1 à 20, on comparera les prédictions du modèle pour le pas 21 avec les observations laissées de côté.

Il faut faire attention à un problème éventuel: le script peut réaliser la validation croisée seulement avec des pas de temps qui font partie du domaine temporel de l’ajustement. Autrement dit, si vous voulez laisser de côte le dernier pas de temps des observation, il faut absolument que le paramètre forecast soit au moins égal à 1.

Dans notre cas le domaine temporel va de 1 à 21. Si on laisse de côté le pas 21, le domaine temporel ira de 1 à 20. Le modèle ne pourra pas utiliser le pas 21 qui est en dehors du domaine.

Par contre si on laisse de côté le pas 21 et qu’on demande forecast=1, le domaine du modèle sera de 1 à 21 (1 à 20 à cause des observations et 21 à cause de forecast).

Dernier élément de choix, rappelez-vous qu’IDE est basé sur la modélisation du passage d’un état au suivant. Même si toute l’étendue temporelle est utilisée pour ajuster le modèle, essayez de valider avec le pas temporel le plus proche de la période qui vous intéresse. Dans notre exemple, on suit au cours du temps la colonisation d’un territoire par une espèce extérieure. Si nous voulons prédire l’expansion à venir, il est plus approprié de valider l’ajustement du modèle aux dernières observations disponibles, plutôt qu’aux observation d’il y a 10 ou 15 ans.

Dans notre exemple on fixe « pas de temps de validation » = 21

Validation croisée d’un échantillon aléatoire

Outre le développement futur du processus d’expansion de notre oiseau, si on s’intéresse à la dynamique d’expansion, il est important de pouvoir interpoler et extrapoler les observations passées, et ce depuis le début des observations.

Pour effectuer une validation croisée appropriée, on prélèvera de manière aléatoire un certain nombre d’observations, réparties sur l’ensemble du territoire et sur l’ensemble du domaine temporel.

Ce « certain » nombre est directement dépendant du nombre d’observations. Sans notre cas, nous allons extraire 10% des observations, mais le choix sera adapté à chaque cas.

On laisse alors « pourcentage points aléatoires de validation » à sa valeur par défaut ’10).

Répertoire résultats

Dans le premier cas nous n’avions pas défini de répertoire car il n’y avait qu’un graphique en sortie. Avec la validation croisée on a un fichier pdf avec les graphiques et tableaux de validation, ainsi que deux couches shp. Si vous ne définissez pas de répertoire de sortie, tout ceci se trouvera dans votre répertoire « Documents ». Et les fichiers comportant le même nom sont écrasés à chaque exécution.

Il est donc conseillé de définir un répertoire pour les sorties du script.

Nom du fichier

Cette chaîne de caractères sera ajoutée aux noms des couches et au fichier pdf. Vous pourrez alors identifier les différentes exécutions du script dans le répertoire de sortie et, les noms n’étant pas les mêmes, il n’y aura pas écrasement des sorties précédentes.

Résultats de l’exécution du script

Le document pdf contient les corrélations entre les prédictions pour chaque bloc de données (observations, bloc pas de temps, échantillon aléatoire) et les valeurs observées.

Cette information donne juste un aperçu de l’ajustement du modèle. Les véritables informations sont les différents scores calculé qu’on trouve sous forme de tableau dans le même document pdf.

Pour l’instant nous avons ds scores pour deux types de validation du même et unique modèle. Nous ne sommes pas encore à l’étape de sélection du meilleur modèle. Celle-ci pourra avoir lieu seulement quand nous aurons ajusté plusieurs modèles .

Ce que nous pouvons déduire ici sont les performance de notre modèle par rapport aux deux types d’objectifs dont nous avons parlé plus haut: la prédiction d’un état à venir ou l’augmentation de l’information passée.

Étant donné que les scores sont meilleurs à tour de rôle pour chacun des deux lots de données de validation, on peut difficilement trancher. Rappelons que pour le biais, la valeur plus proche de 0 est la meilleure; pour le SCV c’est la valeur plus proche de 1, et pour les autres scores c’est la plus petite des deux qui est la meilleure.

Un dernier tableau donne la valeurs des critères d’information AIC et BIC, mais ces indicateurs concernent le modèle et pas la validation croisée. Ils n’ont d’utilité que pour la comparaison avec d’autres modèle.

Quel autre modèle? Comme nous l’avons vu plus haut, si on observe les résultats année par année, ont voit une dérive progressive vers le sud-ouest. On peut donc tester l’ajustement d’un modèle IDE avec un noyau variant spatialement et le comparer au modèle que nous venons de valider, qui a un noyau invariant.

Troisième cas: modèle variant avec validation croisée.

La « seule » différence entre ce type de modèle et le précédent réside dans le fait d’ajuster une dérive des paramètres thêta. On doit donc ajuster chacun des 4 paramètres (amplitude, aperture, xshift et yshift) à partir d’une série de fonctions.

Tout d’abord, on a le choix entre 4 familles de fonctions: bicarrées, gaussiennes, exponentielles et Matern32. Si vous n’avez pas une raison particulière de faire un autre choix, restez-y aux fonctions bicarrées.

Deuxième point, il faut définir une échelle pour la base de la fonction, la scale_aperture.

Par défaut, ce paramètre a une valeur de 1. Mais vous devez vérifier si cette valeur correspond à vos données.

Pour cela, vous avez un troisième script, ide_auto_basis, qui permet de voir les emprises des fonctions choisies en fonction de la scale_aperture.

Si vous exécutez ce script sur les données de test, avec basis=bisquare et scale_aperture=1 vous obtenez la sortie suivante:

Les ellipses correspondent aux emprises des 12 fonctions bicarrées proposées par le modèle.

Les points correspondent aux observations. Vous remarquerez que pour toute observation il y a au moins deux fonctions qui serviront à l’estimation des paramètres à ces endroits. Sauf dans le coin supérieur gauche, ou nous avons trois points d’observation couverts par une seule fonction. L’ajustement de ces trois points est très compliqué. Il est préférable d’avoir au moins deux fonctions par point et même trois. Pour cela, nous exécutons à nouveau le script, mais en agrandissant la scale_aperture à 1.5. Le résultat est alors:

Maintenant toutes les observations peuvent être ajustées avec au moins trois fonctions bicarrées. C’est donc cette échelle qu’on entrera dans le script d’ajustement du modèle. Juste à titre d’exemple, voici le résultat du script ide_auto_basis pour la famille de fonctions exp et notre échelle de 1.5.

Vous voyez que pour chaque cas il faut tester le résultat obtenu, car il n’y aura pas de message d’erreur, simplement les résultats ne seront pas bons.

Reste maintenant à décider si nous considérons que tous les paramètres thêta dérivent entre deux pas de temps,ou bien s’il y en a que certains qui le font.

L’amplitude et l’aperture correspondent à des paramètres évoluant plus dans le temps que dans l’espace. tandis que xshift et yshift, au contraire, évoluent plus dans l’espace que dans le temps.

Dans notre exemple on n’observe pas forcément un accroissement des patchs ou une augmentation sensible d’année en année des valeurs à un point donné. On observe par contre une dérive spatiale (en X et Y).

Nous allons donc laisser l’amplitude et l’aperture invariantes (option constant_basis()) spatialement et nous allons définir xshift et yshift comme variantes (option auto_basis).

Le script va ajuster chacun des deux paramètres avec 12 fonctions chacun, réparties sur l’ensemble du territoire.

Résultats de l’exécution du script

Commençons par les résultats des corrélations prédictions-observations.

Comparons avec les résultats du modèle sans noyau variant spatialement:

Globalement et pour un pas de temps donné, les résultats sont comparables. Par contre pour estimer des valeurs prises aléatoirement sur toute l’historique des observations, le modèle paraît moins performant.

On observe les différents scores:

Les différents scores confirment notre première observation. Le deuxième modèle est moins performant quant il s’agit de « boucher les trous » tout au long de la période d’observation.

En ce qui concerne la prédiction d’un pas de temps, les résultats sont moins faciles à interpréter. Si on se base sur les scores « univariés » (Biais, PCV, SCV, CRPS) le résultat de notre deuxième modèle est légèrement moins bon. Par contre si on regarde les scores « multivariés » (ES, VS) c’est le contraire.

Dernier élément, les indices d’information:

Il n’y a pas de test statistique pour comparer les indices d’information. La valeur de l’AIC indique que la log-vraisemblance des deux modèles est pratiquement la même. L’indice BIC est moins favorable pour le deuxième modèle, mais comme nous avons ajouté l’ajustement de 24 nouvelles fonctions pour modéliser le noyau variant, ceci n’a rien d’étonnant et n’indique pas forcément que le modèle global soit moins bien ajusté.

Conclusion: si le but de notre modèles est de « boucher les trous » des observations, il vaut peut-être mieux utiliser le modèle à noyau invariant. Si notre but est de prévoir l’état futur de la population, tous les deux semblent fonctionner a peu près de la même manière.

Mais ceci reste à décider d’un point de vue moins mathématique. Si on affiche les résultats du modèle pour les mêmes pas de temps (20 et 21) montrés plus haut pour le cas 2:

Finalement, c’est la connaissance du terrain et du phénomène étudié qui fera pencher la balance, entre une population plutôt hétérogène (cas 2) et une population graduée assez régulièrement.

Un dernier élément que nous n’avons pas encore vu, c’est la modélisation du noyau.

Rappelons d’abord qu’est-ce que le noyau représente. A chaque pas de temps, la population augmente. La partie noyau d’IDE représente où est-ce qu’on va retrouver les descendants nés à un endroit au moment t, au moment t+1.

L’ajustement de notre modèle sur xshift et yshift, indique que majoritairement, les nouveaux individus ont une tendance à se déplacer vers la partie SE de notre territoire.

Si on regarde la population sur les 20 ans d’observations, on remarque que celle-ci progresse vers le NW du territoire. Apparemment ceci est contradictoire avec le résultat de notre modèle. Mais il est tout à fait possible que les deux processus ne soient pas si contradictoires. Si on considère que les conditions au SE du territoire sont les meilleures, il est logique que les nouveaux venus aient un tropisme vers cet endroit. Mais comme c’est là où l’on trouve le plus d’individus, c’est là où la compétition pour les ressources et pour les partenaires de reproduction sont les plus dures. Donc, petit à petit la population grandissante est obligée de coloniser de nouveaux territoires, moins propices. C’est ce qui pourrait expliquer que la colonisation soit si lente.

Si cet article vous a intéressé et que vous pensez qu'il pourrait bénéficier à d'autres personnes, n'hésitez pas à le partager sur vos réseaux sociaux en utilisant les boutons ci-dessous. Votre partage est apprécié !

Atilio Francois