Cet article est le mode d’emploi des scripts R pour QGis mis à votre disposition. Il a été mis à jour pour prendre en compte les modifications apportées à la version 1.3 des scripts.
Dans deux articles précédents (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 l’article précédent, 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.
Nous verrons ici comment utiliser les scripts avec l’option validation du modèle, et, dans le prochain article, un exemple pratique avec les mêmes 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
Les fichiers shapefile sont à copier dans un répertoire de données de votre choix.
Les deux scripts principaux
Vous avez deux scripts qui permettent de calculer et ajuster un modèle IDE à vos données spatio-temporelles.
Le choix de l’un ou l’autre se fait selon l’idée que vous vous faites du phénomène à modéliser. La modélisation IDE, comme nous ‘avons vu dans les deux premiers articles, donne une importance spéciale au pas de temps précédent celui qu’on doit modéliser.
D’un point de vue spatial, si le modèle s’applique de la même manière, quelle qu’il soit l’endroit qu’on considère, on dit qu’il a un noyau invariant spatialement.
Par contre si entre deux pas de temps on peut observer un glissement ou un changement dans l’application du modèle, dépendant de l’endroit qu’on observe, on dit que le modèle est à noyau variant spatialement.
Si vous observez les deux interfaces des scripts, vous constaterez que la première partie des paramètres est commune. Nous allons d’abord voir cette partie, puis les paramètres spécifiques à chacun des deux scripts.
Les entrées des scripts
Définition des données – paramètres t et Val
Puisque nous allons modéliser des données spatio-temporelles, il est clair que nous devons commencer par établir quelle donnée nous intéresse et quel est l’attribut temporel à prendre en compte.
Le paramètre Val correspond à l’attribut à modéliser. Il doit être forcément numérique (réel ou entier).
Le paramètre t est plus délicat. Il détermine la chronologie de nos observations de Val. Le principal problème vient de la multiplicité des types d’unités et de formats permettant d’établir une chronologie.
La manière la plus simple et efficace est d’avoir un attribut entier qui donne le rang de chaque pas de temps. Vous avez 20 années d’observations annuelles, vous les numérotez de 1 à 20, et vous avez votre paramètre t qui sera traité par les scripts sans aucun problème.
Une autre forme très courante est d’avoir un attribut de type string avec des dates de type Y-M-D ou Y/M/D ou D/M/Y, etc.
Les packages R utilisés dans les scripts chercheront à convertir la chaîne de caractère en date R. S’ils arrivent, pas de soucis. Si le script s’arrête avec un message d’erreur, il est conseillé de passer à la première méthode en ajoutant un nouvel attribut avec le rang.
La dernière manière d’avoir la chronologie dans nos données est d’avoir déjà un champ au format Date. Mais là aussi, les données source pouvant être dans une multiplicité de formats différents, les formats de Date sont aussi nombreux.
Les scripts chercheront à traduire le format DATE des données en format Date R. Même remarque que précédemment:Si le script s’arrête avec un message d’erreur, il est conseillé de passer à la première méthode en ajoutant un nouvel attribut avec le rang.
Covariable (A partir de la version 1.3)
La modélisation du paramètre Val défini plus haut peut prendre en compte un autre attribut comme co-variable. Si on étudie la distribution d’une espèce dans la zone littorale, un autre paramètre peut intervenir pour expliquer la densité des individues, par exemple la profondeur. L’espèce aura une fourchette de profondeur dans laquelle elle va se développer,. On peut donc utiliser les valeurs de profondeur pour améliorer la modélisation.
Le seul inconvénient est un problème la bibliothèque R IDE. La documentation étant très sommaire, je n’ai pas trouvé(au moins pour l’instant) comment produire la grille de sortie (grid_size x grid_size) lorsqu’on utilise une co-variable car, apparemment, la fonction recherche une valeur de la co-variable pour chaque point de la grille. Comme les valeurs sont disponibles seulement aux points d’observation, le calcule est bloqué. Les sorties sont donc limitées aux couches contenant les observations, le bloc de pas de temps et l’échantillonnage aléatoire.
Ajouter lon/lat comme co-variables (A partir de la version 1.3)
De même que dans le cas de la profondeur, dans l’exemple précédent, on peut considérer que la densité de notre paramètre à modéliser n’est pas distribuée au hasard sur le territoire. C’est le cas dans les données de test avec les comptages d’oiseaux, où on voit un gradient diagonal NW-SE. La prise en compte des coordonnées XY comme co-variable pourra aider à une meilleure modélisation de ce paramètre.
Contrairement aux covariables « attribut », il y a une valeur de XY pour chaque point de la grille. On n’a donc pas le même problème pour produire la grille en sortie quand cette option est cochée.
Le paramètre unités
Les outils de modélisation nécessitent des pas de temps de même longueur. Si nos données correspondent à des observations mensuelles, le problème est que les mois n’ont pas la même longueur… de 28 à 31 jours. Et si nos données sont annuelles, les années n’ont pas la même longueur (bissextiles ou pas). Ceci fait que les unités de pas de temps acceptables pour les modèles spatio-temporels vont de la seconde à la semaine, pas plus.
Utilisez les différentes options si vous avez déjà un attribut de type date. Si non l’option « days » est celle qui doit être appliquée pour les données de type rang.
Forecast
L’ajustement du modèle se fait avec les pas de temps des observations. Le résultat se présente sous la forme d’une grille avec le même nombre de pas de temps que les observations.Le paramètre forecast permet de demander au modèle de prédire les valeurs pour des pas de temps ultérieurs à ceux des observations. Le nombre introduit est le nombre de pas à prédire. Si nous avons 20 pas de temps dans les observations, une valeur de forecast de 1 va prédire le pas de temps numéro 21, une valeur de 2 prédira les pas de temps 21 et 22, etc.
Paramètre Itermax
L’ajustement du modèle est un processus itératif. Il prend beaucoup de temps de calcul. (beaucoup!). Parallèlement, nous savons que le travail de modélisation est aussi itératif en ce qui concerne la recherche l’utilisateur. Nous allons lancer l’ajustement un grand nombre de fois. Si la première fois on n’a aucune idée du temps de calcul et du nombre « nécessaire » d’itérations, après la première exécution nous disposons du journal du script qui trace les valeurs des paramètres du modèle obtenues après chaque itération.
Quand on a de la chance, le modèle arrive à une convergence idéale en moins d’itérations que le nombre demandé. La plupart du temps il n’atteint pas cette convergence idéale. La lecture des valeurs successives permet de se faire une idée du besoin, soit d’augmenter le nombre d’itérations si les dernières valeurs des estimations des paramètres continuent à montrer une variation, soit de diminuer le nombre d’itérations si les valeurs ne montrent pas de variation.
Les sorties des scripts
Résultat »grille » (Modifié à partir de la version 1.3)
Première sortie par défaut des scripts, elle consiste en une couche de points correspondants à la grille définie par le paramètre grid_size. Pour chaque point il y a une occurrence par pas de temps des observations et, si le paramètre forecast est différent de 0, une occurrence par forecast demandé.
Les attributs contiennent:
- le pas de temps
- la valeur de prédiction du processus Y pour le point (Predict)
- l’erreur standard du processus Y pour le point (ErreurY)
- l’erreur standard dans l’espace des observations Z (ErreurZ)
Étant donné que cette couche n’est pas produite systématiquement (elle n’est pas calculable dans le cas d’utilisation d’une co-variable attribut), elle n’est pas chargée automatiquement à la fin du script. Elle est produite dans le répertoire résultats.
Résultat « observations »
Seulesortie par défaut, elle consiste en une couche de points correspondants aux points en entrée, les observations. Même si le paramètre forecast est différent de 0, il n’y a pas de calcul pour des pas de temps non présents dans les observations.
Les attributs contiennent:
- le pas de temps
- la valeur de l’observation dans le fichier en entrée (Val)
- la valeur de prédiction du processus Y pour le point (Predict)
- l’erreur standard du processus Y pour le point (ErreurY)
- l’erreur standard dans l’espace des observations Z (ErreurZ)
Répertoire résultats et nom de fichier
Si l’option « effectuer la validation du modèle » est cochée trois autres résultats, en plus des couches par défaut, sont produits:
- un fichier pdf (Rplots.pdf) avec les graphiques et les tableaux de scores de validation
- une couche au format shapefile correspondant aux points en entrée qui auront été mis de côté en tant que bloc de pas temporel (_bloc.shp)
- une couche au format shapefile correspondant aux points en entrée qui auront été mis de côté en tant qu’échantillon aléatoire de points (_alea.shp)
- une couche au format shapefile correspondant à la grille de calcul, si elle est produite (_grille.shp)
Les scores de validation calculés pour chaque ensemble de données de validation (bloc et aléatoire) sont:
- corrélation prédictions/observations
- le Biais
- la PCV : la validation croisée prédictive
- la SCV : la validation croisée standardisée
- le CRPS : le score de probabilité classé en continu
- l’ES: le score d’énergie
- le VS: le score du Variogramme d’ordre p
- l’AIC: le critère d’information d’Akaike
- le BIC : le critère d’information Bayesien
Les attributs de ces deux couches sont les mêmes que ceux de la couche Résultat observations
Vous pouvez définir un répertoire pour ces fichiers avec le paramètre « Répertoire résultats« . Si le champ reste vide, les fichiers sont placés dans votre répertoire « Documents« .
Si les fichiers existent déjà, ils sont remplacés par les nouveaux.
Le champ « Nom du fichier » permet d’entrer une chaîne de caractères qui sera ajouté aux noms des couches bloc et alea.
Configuration de la validation croisée du modèle
Comme nous l’avons détaillé dans l’article précédent, la validation croisée du modèle implique la mise de côté de certaines observations pour le calcul et l’ajustement du modèle, puis de tester le modèle calculé sur ces observations laissées de côté.
Le champ « Pas de temps de validation » met de côté l’ensemble des observations correspondantes au pas de temps indiqué. Si le paramètre t est un rang, l’entrée doit être un entier correspondant au rang souhaité. Si ce n’est pas le cas, vous devez rentrer le contenu exact du pas de temps dans les données en entrée en respectant le format.
Le champ « Pourcentage points aléatoires de validation » est un entier qui indique la proportion de points, échantillonnés au hasard sur l’ensemble des pas de temps,excepté le pas de temps de validation s’il a été renseigné, et de l’espace de travail.
Remarque: échantillonnage au hasard sera identique à chaque exécution du script pour le même fichier d’observations en entrée. Le paramètre de R, set.seed(1), contrôle la génération des nombres au hasard, reproductibles à chaque exécution. Si vous voulez un échantillonnage différent lors d’une exécution, éditez le script et changez la valeur 1 par une autre valeur quelconque. Tant que cette valeur ne sera modifiée, l’échantillonnage sera le même à chaque exécution.
Paramètres spécifiques au modèle variant spatialement
L’explication de ces paramètres a été donnée dans l’article précédent. On ne reviendra donc pas sur le sens de ceux-ci.
Le paramètre Basis permet de sélectionne la famille de fonctions à utiliser pour les paramètres thêta, parmi les quatre familles: bisquare, Gaussian ,exp et Matern32.
Le paramètre nresolution permet de définir différentes tailles d’emprise des fonctions de base. Soyez prudent car passer de nresolution=1 à nresolution=2 multiplie par environ 10 fois le nombre de fonctions à ajuster, pour chaque paramètre. Les temps de calcul peuvent alors être extrêmement longs.
Le paramètre scale_aperture contrôle l’emprise à la base des fonctions de base. Si le type de famille est bisquare, une aperture de 1 ou 1.25 sera suffisante. Pour les autres familles de fonctions, il est conseillé d’utiliser le script ide_auto_basis pour trouver la valeur de scale_aperture approprié.
Pour les quatre paramètres thêta: amplitude, aperture, xshift et yshift vous pouvez choisir entre deux options:
- constant_basis(): le paramètre est considéré invariant et sera calculé comme si le noyau était invariant spatialement
- auto_basis: le paramètre est considéré variant et sera calculé à partir de l’ajustement des fonctions de base demandées
Dans le prochain article on verra l’exemple appliqué à nos données de test.
