Note sur ce template Rmarkdown
Ce fichier constitue un gabarit complet pour la modélisation des habitats potentiels d’une espèce ou d’un groupe d’espèce. Il est fourni sous licence libre CC-BY 4.0.
Vous êtes autorisé à :
Partager — copier, distribuer et communiquer le matériel par tous moyens et sous tous formats pour toute utilisation, y compris commerciale. Adapter — remixer, transformer et créer à partir du matériel pour toute utilisation, y compris commerciale. L’Offrant ne peut retirer les autorisations conférées par la licence tant que vous appliquez les termes de cette licence.
Selon les conditions suivantes :
Attribution — Vous devez créditer ce travail, intégrer un lien vers la licence et indiquer si des modifications ont été effectuées à ce code. Vous devez indiquer ces informations par tous les moyens raisonnables, sans toutefois suggérer que l’Offrant vous soutient ou soutient la façon dont vous avez utilisé son code.
Partage dans les Mêmes Conditions — Dans le cas où vous effectuez un remix, que vous transformez, ou créez à partir du matériel composant le code original, vous devez diffuser le code modifié dans les même conditions, c’est à dire avec la même licence avec laquelle le code original a été diffusé.
Pas de restrictions complémentaires — Vous n’êtes pas autorisé à appliquer des conditions légales ou des mesures techniques qui restreindraient légalement autrui à utiliser le code dans les conditions décrites par la licence.
Il a été testé fonctionnel sur la version de R R version 4.3.2 (2023-10-31 ucrt), avec les packages sf (1.0.15), terra (1.7.71), ggplot2 (3.4.4), scales (1.3.0), egg (0.4.5), virtualspecies (1.6), blockCV (3.1.3), biomod2 (4.2.5), dplyr (1.1.4), tidyterra (0.5.2), viridis (0.6.5).
Il est possible que des évolutions futures de packages (notamment, biomod2, qui est sujet à de nombreuses évolutions en 2023 et 2024) rendent certaines parties du fichier non fonctionnelles, ce qui nécessitera de corriger le code.
Pré-requis :
Chargement des packages et fonctions, chargement de données géographiques et des variables environnementales harmonisées
library(sf)
library(terra)
library(ggplot2)
library(scales)
library(egg)
library(virtualspecies)
library(blockCV)
library(biomod2)
library(dplyr)
library(tidyterra)
library(viridis)
source("scripts/functions.R")
# Shapefile de la Corse
corse <- st_read("data/corse.gpkg")## Reading layer `corse' from data source `C:\Rprojects\SDMs_PNA_Corse\data\corse.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 8.534717 ymin: 41.33323 xmax: 9.560364 ymax: 43.02755
## Geodetic CRS: WGS 84Chargement et préparation des données d’occurrence
## Reading layer `sittelle' from data source
## `C:\Rprojects\SDMs_PNA_Corse\data\donnees_brutes\taxa\sittelle.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3136 features and 69 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 8.62913 ymin: 41.59261 xmax: 9.335384 ymax: 42.52245
## Geodetic CRS: WGS 84# Simplification du nom d'espèce en binomial
sittelle$species <- simplify_species_name(sittelle$nom_valide)
# Dates d'échantillonnage
sittelle$year <- as.numeric(strtrim(sittelle$date_fin, 4))
sittelle$month <- as.numeric(substr(sittelle$date_fin, 6, 7))
# Visualisation de la temporalité des occurrences
ggplot(sittelle) +
geom_boxplot(aes(x = species,
y = year))+
coord_flip() +
scale_y_continuous(breaks = breaks_pretty()) +
theme_minimal()
Filtre temporel
Il faut établir un filtre temporel pour éliminer les données
imprécises, sachant que l’objectif est de modéliser à une résolution
assez fine, de l’ordre de 1km. Le champ precision est peu
renseigné et donc peu utile ici, il nous faut donc poser une hypothèse
sur les données qui sont imprécises. On peut considérer que les GPS ont
commencé à être largement disponibles à partir de 1990, mais leur
utilisation ne s’est généralisée qu’à partir des années 2000, notamment
grâce à leur miniaturisation. Ainsi, on peut spéculer qu’avant les
années 2000, les données étaient moins précisés car possiblement
géolocalisées en utilisant des référentiels comme les lieu-dits ou les
communes, tandis qu’à partir des années 2000 la précision s’est
améliorée grâce à la géolocalisation par satellite.
La quasi-totalité des données a été récoltée après 2000 :
# Les données avant 2000 représentent un % modéré du jeu de données :
100 * length(which(sittelle$year < 2000)) / nrow(sittelle)## [1] 0.1594388La couverture temporelle sur l’année est plus importante à partir de 2010 environ :
ggplot(sittelle) +
geom_boxplot(aes(x = species,
y = month)) +
facet_wrap(~year) +
coord_flip() +
scale_y_continuous(breaks = breaks_pretty()) +
theme_minimal()
L’emprise spatiale des données d’occurrence change de manière modérée avec ou sans les données pré-2000 :
p_sittelle_all <- ggplot() +
geom_sf(data = corse) +
geom_sf(data = sittelle, aes(col = year)) +
scale_color_continuous(type = "viridis") +
theme_minimal(base_size = 15) +
ggtitle("Toutes données\nsittelle")
p_sittelle_post2000 <- ggplot() +
geom_sf(data = corse) +
geom_sf(data = sittelle[sittelle$year >= 2000, ], aes(col = year)) +
scale_color_continuous(type = "viridis") +
theme_minimal(base_size = 15) +
ggtitle("Données post-2000\nsittelle")
ggarrange(p_sittelle_all,
p_sittelle_post2000,
nrow = 1)
On pose donc l’hypothèse raisonnable qu’un filtre à 2000 va assurer une bonne précision dans la localisation des occurrences sans perdre d’information critique sur la répartition de la sittelle
Rasterisation des occurrences
L’objectif ici est de ne garder qu’une occurrence par cellule à la résolution de nos variables environnementales afin d’éviter une forme extrême de pseudo-réplication. Par exemple, si dans une cellule donnée on a initialement 30 observations de la même espèce, alors, après rasterisation, ces 30 observations ne compteront que comme une seule occurrence. Cette étape est indispensable car elle évite de donner aux modèles, par exemple, 30 fois la même valeur de température provenant d’une seule cellule. C’est ce qu’on appelle de la pseudo-réplication et c’est très problématique pour les domaines. On s’attend donc à ce que cette étape réduise le nombre d’occurrences pour les modèles.
# On rasterise les occurrences à la résolution de nos variables
# environnementales
sittelle_r <- rasterize(sittelle,
env_corse)
names(sittelle_r) <- "sittelle" # Attention il ne faut pas nommer
# la couche "sittelle" car il y a des variables qui s'appellent
# sittelle
plot(sittelle_r)
On va ensuite éliminer les occurrences qui sont dans des zones sans
valeurs de variables environnementales (i.e., essentiellement en zones
côtières). Pour cela on va combiner les variables. Pour cela on va
combiner les variables environnementales avec les occurrences
rasterisées dans un data.frame, et supprimer les
occurrences d’espèces qui tombent sur des données environnementales
manquantes
# On crée un stack avec nos occurrences rasterisées et les variables env
env_sittelle <- c(env_corse,
sittelle_r)
# On récupère les coordonnées XY de toutes les cellules, pour préparer nos
# données finales
coorXY <- xyFromCell(env_corse,
1:ncell(env_corse))
# On transforme le raster en data.frame
env_sittelle_df <- values(env_sittelle)
env_sittelle_df[is.nan(env_sittelle_df)] <- NA
# On regarde le nombre d'occurrences pour lesquelles il y a des données
# manquantes :
length(which(is.na(env_sittelle_df[, "bio1"]) &
!is.na(env_sittelle_df[, "sittelle"])))## [1] 0On va maintenant supprimer les cellules pour lesquelles on n’a pas de données environnementales. Pour cela on va utiliser la première variable environnementale ici, car les données manquantes sont toutes les mêmes entre toutes les variables environnementales (cf. script harmonisation des données).
# On filtre d'abord sur l'objet qui contient les coordonnées
coorXY <- coorXY[-which(is.na(env_sittelle_df[, 1])), ]
# Et ensuite sur le tableau avec variables env et présences d'espèces
env_sittelle_df <- env_sittelle_df[which(!is.na(env_sittelle_df[, 1])), ]
# Comparaison du nombre d'occurrences :
# Avant rasterisation
nrow(sittelle)## [1] 3131# Après rasterisation et élimination des données env manquantes
length(which(env_sittelle_df[, "sittelle"] == 1))## [1] 552Il s’agit donc du nombre d’occurrences que l’on va pouvoir utiliser pour calibrer nos modèles. Il y a 552 occurrences ce qui est assez élevé pour la calibration des modèles.
On va maintenant formater ces occurrences en combinant coordonnées et
info sur l’occurrence dans un data.frame pour préparer la
calibration de nos modèles
P_points <- data.frame(
# D'abord on récupère les coordonnées XY qui correspondent à nos cellules de présences
coorXY[which(!is.na(env_sittelle_df[, "sittelle"])), ],
# Ensuite, on récupère la colonne qui indique présence pour chaque cellule
occurrence = env_sittelle_df[which(!is.na(env_sittelle_df[, "sittelle"])),
"sittelle"])
P_pointsGénération des points de background
Etant donné que nos observations sont des présences-seules, i.e. sans données d’absences, il nous faut générer des points de “background” pour pouvoir calibrer les modèles. Ces points de backgrounds sont des données tirées dans toute la zone d’étude qui renseignent les modèles sur comment les variables environnementales sont distribuées dans la géographie. Ces points seront fournis aux modèles comme des 0, ce qui permettra aux modèles d’identifier quels habitats apparaissent comme favorable parmi l’ensemble des habitats disponibles. Cependant, ces 0 ne sont pas interprétés comme des absences, et l’interprétation finale du modèle nécessitera des précautions particulières, comme par exemple ne pas considérer la valeur issue du modèle comme une “probabilité de présence” ; elle sera plutôt considérée comme un indice de favorabilité de l’habitat.
La littérature statistique récente suggère que les meilleures pratiques consistent à générer un grand nombre de points de background (e.g., 10000) indépendamment de la localisation des points de présence (i.e., un point de background peut être localisé au même endroit qu’un point de présence). Cela permet d’assurer une bonne représentation de l’ensemble des conditions environnementales disponibles dans le modèle. Dans le cas de la Corse, le nombre de points de background sera limité par le nombre de pixels disponibles :
# Nous avons éliminé les données manquantes du tableau env_sittelle_df
# Par conséquent, son nombre de lignes est égal au nombre total de pixels
# disponibles sur la Corse
nrow(env_sittelle_df)## [1] 13620Ainsi, nous partons sur un point de départ à 10000 backgrounds ce qui sera suffisant pour une bonne calibration des modèles. Il n’est pas nécessaire de faire plusieurs répétitions, car le nombre de points de background est déjà suffisamment élevé, les résultats de calibration ne varieraient pas entre différentes répétitions.
Les protocoles d’échantillonnages de la sittelle ont été établis pour limiter le biais d’échantillonnage lié à l’accessibilité aux routes. Ainsi, nous ne chercherons pas à corriger ce biais dans la génération des backgrounds ici.
# On réduit également le nombre de background pour avoir un effet du biais
background <- spatSample(env_corse,
method = "random",
size = 10000,
replace = FALSE, # Pas de remise
na.rm = TRUE, # Pas dans les données manquantes
xy = TRUE, # L'output inclut les coords XY
values = FALSE) # L'output exclut les variables
# On ajoute les points de background aux données de présence
P_points <- rbind.data.frame(P_points,
data.frame(background,
occurrence = 0))
# Affichage des occurrences
plot(P_points$y ~ P_points$x, pch = c(1, 16)[P_points$occurrence + 1],
asp = 1, cex = .5,
xlab = "Longitude", ylab = "Latitude")
Les points de background sont les cercles blancs, et les occurrences sont les cercles pleins.
Sélection des variables environnementales
Climat
La sittelle Corse dépend des conditions climatiques, notamment température et précipitations, pour sa survie et sa reproduction (Thibault et Villard 2005, Thibault et al. 2006). Cependant, la distribution actuelle de la sittelle Corse est probablement le reflet d’une contraction très importante d’aire de répartition à basse altitude liée à deux facteurs ; d’une part, la surexploitation de son habitat principal, le pin Laricio (Barbet-Massin et Jiguet 2011, Torre 2014), et d’autre part, l’augmentation des feux de forêts qui a favorisé la régénération du pin maritime (Guy, pers. comm. 2023).
Ainsi, il est fort probable que la recherche d’une relation en analysant la corrélation entre le climat et la distribution actuelle de la sittelle Corse résulte en l’identification de réponses climatiques qui ne reflètent pas réellement la niche écologique de l’espèce. Par conséquent, les choix de variables climatiques que nous effectuerons pour modéliser l’habitat de la sitelle Corse seront seront avant tout contraint par l’objectif d’éviter de trouver des corrélations irréalistes dûes aux à la réduction anthropogénique de l’aire de répartition de l’espèce. Ainsi, nous partirons de l’hypothèse que la limite altitudinale inférieure de la sittelle est due à une réduction anthropogénique de l’habitat et non pas à un climat trop chaud et sec - nous n’incluerons donc pas de variables climatiques liées aux températures maximales et à la sécheresse. Par contre, il est probable que la limite altitudinale supérieure de la sittelle Corse soit effectivement dûe à des caractéristiques climatiques empêchant la survie ou le développement de l’espèce (e.g., accès aux graines de pin en hiver liées à l’ouverture des cônes selon la température extérieure, ONF 2016) ; nous utiliserons donc des variables reflétant les minimums de températures annuels en supposant un effet négatif des températures trop froides en hiver. Un facteur confondant pourrait néanmoins également exister pour la limite supérieure : la limite altitudinale haute du pin laricio était contrainte par le pâturage estival, et il est possible que l’abandon de cette pratique permet aujourd’hui une reconquête du pin en altitude (Guy, pers. comm. 2023). Il existera donc une incertitude sur la limite altitudinale supérieure.
Noms des variables retenues :
- Températures minimales (bio6)
Occupation du sol
Les travaux sur l’écologie de la sittelle Corse sont formels sur le fait que le pin Laricio constitue l’habitat de la sittelle ; ainsi la distribution des pins Laricio sera utilisée comme variable prédictive avec l’hypothèse qu’il s’agit du principal facteur expliquant la répartition de la sittelle Corse. Cependant, les travaux précédents ont montré que l’exploitation forestière et la nature des parcelles de pin Laricio (notamment la taille des arbres et la présence d’arbres morts) influencent fortement la présence ou non des sittelles ; mais ces informations ne sont pas disponibles dans les variables spatialisées. Ainsi, il est possible que la distribution connue des pins Laricio surestime la répartition réelle de la sittelle Corse.
Des données LIDAR sont en cours d’acquisition sur la nature des peuplements forestiers (Guy, pers. comm. 2023), ce qui permettra dans le futur d’intégrer une couche spatialisée indiquant la nature des peuplements forestiers, avec notamment la présence de gros bois.
Les feux de forêts pourraient influencer négativement la distribution de la sittelle Corse (Guy, pers. comm. 2023) - cependant il n’existe pas de variable environnementale spatialisée sur les feux de forêts dont l’intervalle temporal correspond aux données d’occurrence d’espèces, ce qui rend difficile l’inclusion de cette variable pour le moment.
Noms des variables retenues :
- Répartition du pin laricio (laricio)
Biais d’échantillonnage
La probabilité d’observer les espèces est souvent directement liée à l’accessibilité du milieu, qui est connue pour être fortement corrélée à la distance aux routes. Nous utiliserons donc la distance aux routes comme proxy du biais d’échantillonnage afin d’éviter que les modèles ne cherchent à expliquer l’accessibilité par les autres variables environnementales.
Le plan d’échantillonnage récent pour la sittelle a été conçu pour limiter ce biais (PNA sittelle de Corse 2017-2026, Guy, pers. comm. 2023). Néanmoins, à l’échelle de la Corse entière et prenant toutes les données depuis 2000, on observe une relation forte entre proximité aux routes et densité des observations. Des tests préliminaires ont déterminé que la correction du biais de proximité aux routes améliorait la performance des modèles.
Les données de présence-absence produites par Endemys ont été échantillonnées en minimisant ce biais (Guy, pers. comm. 2023), avec notamment des points d’écoute ont été effectués loin des accès routiers. Cela permettra de tester la validité des modèles indépendamment du biais d’échantillonnage.
Variables anthropogéniques
La principale perturbation pour les sittelles est le degré d’exploitation forestière qui limite la taille des arbres et la présence de bois mort, nécessaire à l’espèce. Nous utiliserons la variable d’intégrité biophysique des sols comme indicateur de la probabilité d’exploitation forestière, car l’exploitation forestière fait partie intégrante du calcul de de l’indice d’intégrité biophysique des sols (Guetté et al. 2021). L’espèce ne semble pas sujette aux perturbations anthropiques hors régimes de feu et surexploitation forestière, à dire d’expert (Guy, pers. comm. 2023).
Noms des variables retenues :
- Intégrité biophysique des sols (integrite)
Autres variables et commentaires
L’espèce est extrêmement territoriale et sédentaire. Il est probable que le territoire hivernal soit plus large que celui occupé en période de reproduction mais il reste probablement centré sur le territoire de reproduction, ainsi il n’est pas nécessaire de séparer les occurrences estivales et hivernales (Guy, pers. comm. 2023).
Constitution du jeu de variables finales pour la sittelle Corse
Etude de la colinéarité
La colinéarité est la corrélation qui existe entre les variables environnementales. Des variables colinéaires posent des problèmes pour la calibration de nombreux modèles statistiques, donc on s’assure toujours d’éliminer les variables colinéaires avant de faire la calibration.
On étudie la colinéarité entre les variables avec le coefficient de corrélation de Spearman (car certaines variables ne sont pas distribuées normalement), en utilisant un seuil standard de 0.7.
var_groups <- removeCollinearity(env_sittelle,
plot = TRUE,
multicollinearity.cutoff = 0.7,
method = "spearman")## - No multicollinearity detected in your data at threshold 0.7
Il n’y a pas de colinéarité au seuil de 0.7 dans les données, et il n’y a que trois variables ce qui est bon pour les modèles.
Préparation de la stratégie de validation croisée des modèles
Nous ne disposons pas de jeu de données indépendant pour évaluer les modèles. Par conséquent, il nous faut utiliser une procédure de “validation croisée” qui consiste à séparer le jeu de données en deux, une partie sert à la calibration des modèles, et l’autre partie sert à l’évaluation. L’approche classique consiste à faire de découpage de manière aléatoire, mais il a été démontré qu’un découpage aléatoire est suroptimiste car les points de données de calibration sont très proches, spatialement, des points de données d’évaluation.
Pour éviter ce problème de proximité spatiale, nous allons utiliser une procédure dite de “validation croisée spatiale par blocs”. Cette validation croisée par blocs vise à réduire l’autocorrélation spatiale entre jeu de données de calibration et jeu de validation. L’autocorrélation spatiale est le fait que des points proches dans l’espace ont des valeurs de variables environnementales similaires. Eviter l’autocorrélation spatiale entre jeu de calibration et d’évaluation revient à éviter que les valeurs de variables environnementales soient similaires entre calibration et évaluation - cela permet de mieux tester la réelle capacité des modèles à prédire l’habitat favorable aux espèces.
La démarche de validation croisée par blocs est la suivante :
Définir une taille de blocs qui réduit l’autocorrélation spatiale entre calibration et évaluation
Répartir les blocs en plis (“folds”) de calibration et d’évaluation
Vérifier que les plis sont équilibrés, i.e. le nombre de points de calibration doit être similaire entre les plis. Si les plis sont déséquilibrés, recommencer les étapes 1-3 en réduisant la taille des blocs.
Définition de la taille des blocs
Il faut étudier le degré d’autocorrélation spatiale dans les variables environnementales pour avoir une idée de la taille des blocs. La taille des blocs est un compromis entre l’objectif de diminution de l’autocorrélation spatiale et les contraintes des données d’occurrences. En effet, si toutes les occurrences sont localisées dans une petite zone, il ne sera pas possible de viser des blocs trop grands, car on ne pourrait alors pas séparer les points en jeu de calibration et jeu d’évaluation.
# Pour étudier la taille des blocs à viser, il faut d'abord projeter le raster
# en mètres, sinon la fonction de calcul de l'autocorrélation échouera
env_sittelle_l93 <- project(env_sittelle,
"EPSG:2154") # Projection en Lambert 93 ici
# Ensuite on étudie le range d'autocorrélation spatiale
AC_range <- cv_spatial_autocor(env_sittelle_l93,
num_sample = 10000)##
|
| | 0%
|
|========================== | 33%
|
|===================================================== | 67%
|
|===============================================================================| 100%
On obtient initialement un range médian qui est de 14.4 km, ce qui est satisfaisant ici pour réaliser une validation croisée par blocs : il y a beaucoup de blocs, ce qui signifie que la répartition des blocs en plis sera probablement bien équilibrée.
P_points_sf <- st_as_sf(P_points,
coords = c("x", "y"),
crs = "EPSG:4326")
plis_cv <- cv_spatial(x = P_points_sf,
column = "occurrence", # Nom de la colonne des occurrences
k = 6, # Nombre de plis (folds) pour la k-fold CV
size = AC_range$range, # Taille des blocs en metres
selection = "random", # Attribution des blocs aléatoire dans
# les plis
iteration = 50, # Nombre d'essais pour trouver des plis
# équilibrés
biomod2 = TRUE, # Formater les données pour biomod2
r = env_sittelle, # Pour le fond de carte
progress = FALSE,
plot = FALSE) ##
## train_0 train_1 test_0 test_1
## 1 8299 489 1701 63
## 2 8257 490 1743 62
## 3 8138 443 1862 109
## 4 8441 430 1559 122
## 5 8462 458 1538 94
## 6 8403 450 1597 102On voit que nos plis sont plutôt équilibrés :
de 430 à 490 présences pour la calibration
de 62 à 122 présences pour l’évaluation
On peut visualiser la répartition des points de calibration (“Train”) et évaluation (“Test”) pour chaque pli sur la carte suivante :
Cette carte inclut à la fois les présences et les backgrounds.
Dernière étape, biomod2 exige un format particulier pour les plis de validation croisée, donc on va préparer ce format ici :
Calibration des modèles
Tout d’abord on prépare les données pour biomod2.
coorxy <- P_points[, c("x", "y")]
occurrences <- P_points[, "occurrence"]
dir.create("models/sittelle", recursive = T, showWarnings = FALSE)
run_data <- BIOMOD_FormatingData(
resp.name = "sittelle", # Nom de l'espèce
resp.var = occurrences, # Présences + background
expl.var = env_sittelle, # Variables environnementales prédictives
dir.name = "models", # Dossier dans lequel on va stocker les modèles
resp.xy = coorxy, # Coordonnées xy des présences et background
PA.strategy = NULL) # Pas de génération de points de background par biomod##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= sittelle Data Formating -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
##
## ! No data has been set aside for modeling evaluation
## ! No data has been set aside for modeling evaluation
## !!! Some data are located in the same raster cell.
## Please set `filter.raster = TRUE` if you want an automatic filtering.
## ! No data has been set aside for modeling evaluation
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=Biomod nous indique deux choses : que nous n’avons pas de données indépendantes pour l’évaluation, ce qui est effectivement le cas à ce stade de l’étude. Par ailleurs, que plusieurs données peuvent être dans la même cellule, ce qui est également attendu car nous avons tiré aléatoirement nos background dans toute la zone d’étude et donc ils ont pu tomber dans les mêmes cellules que des points de présence. Pas d’inquiétudes, c’est ce que l’on avait prévu.
On va pouvoir désormais préparer la calibration des modèles, en les paramétrant de manière correcte. Ce qui est important de savoir ici c’est que nos modèles vont avoir deux grosses difficultés statistiques :
déséquilibre des classes : il y a au total 552 présences et 10000 backgrounds (qui seront considérés comme des valeurs de 0 par les modèles), ce qui crée un gros déséquilibre entre les 1 et les 0. C’est ce que l’on appelle le déséquilibre des classes
chevauchement des classes : il est probable que les présences et les backgrounds se chevauchent sur les gradients de variables environnementales (d’autant plus que nous pouvons avoir parfois une présence et un background dans le même pixel), ce qui rend la distinction entre les 1 et les 0 difficile pour les modèles. C’est ce que l’on appelle le chevauchement des classes.
La solution pour bien paramétrer les modèles face au déséquilibre et au chevauchement varie selon les modèles, mais le principe général est de réduire l’importance des backgrounds lors de la calibration par rapport au présence, afin de viser un ratio équilibre 50/50 entre importance des présences et importance des backgrounds. Par exemple, on va attribuer des poids aux présences et aux backgrounds de sorte que la somme du poids des présences et des backgrounds soit égale. Cependant, cette méthode fonctionne mal sur certains modèles comme le random forest, et il faut alors le paramétrer de manière plus fine avec un rééchantillonnage à 50/50 en interne.
Par ailleurs, il est important de noter que l’évaluation des modèles avec la validation croisée n’est pas un élément validant la robustesse du modèle. Elle est plutôt à considérer comme un élément qui élimine les mauvais modèles, mais elle ne constitue pas une preuve de robustesse quand elle est bonne, car elle est limitée à la fois par la nature des données (présence-seule, pas d’absences), et par la possibilité qu’il y ait des biais dans l’échantillonnage. Ainsi, il est difficile d’utiliser la validation croisée pour identifier les meilleurs modèles ; il vaut mieux donc se baser sur des paramètres établis pour être robustes en situation de présence-seule (e.g., Valavi et al. 2021).
Préparons donc la calibration de nos modèles :
calib_summary <-
summary(run_data, calib.lines = table_cv) %>%
filter(dataset == "calibration")
iwp <- (10^6)^(1 - occurrences)
RF_param_list <- NULL
GLM_param_list <- NULL
GBM_param_list <- NULL
XGBOOST_param_list <- NULL
XGBOOST_param_list <- NULL
GAM_param_list <- NULL
MARS_param_list <- NULL
XGBOOST_param_list <- NULL
for (cvrun in 1:nrow(calib_summary)) {
prNum <- calib_summary$Presences[cvrun]
bgNum <- calib_summary$True_Absences[cvrun]
wt <- ifelse(occurrences == 1, 1, prNum / bgNum)
RF_param_list[[paste0("_",
calib_summary$PA[[cvrun]],
"_",
calib_summary$run[[cvrun]])]] <-
list(ntree = 1000,
sampsize = c("0" = prNum,
"1" = prNum),
replace = TRUE)
GLM_param_list[[paste0("_",
calib_summary$PA[[cvrun]],
"_",
calib_summary$run[[cvrun]])]] <-
list(weights = wt)
GBM_param_list[[paste0("_",
calib_summary$PA[[cvrun]],
"_",
calib_summary$run[[cvrun]])]] <-
list(interaction.depth = 5,
n.trees = 5000,
shrinkage = 0.001,
bag.fraction = 0.75,
cv.folds = 5,
weights = wt)
GAM_param_list[[paste0("_",
calib_summary$PA[[cvrun]],
"_",
calib_summary$run[[cvrun]])]] <-
list(weights = wt)
MARS_param_list[[paste0("_",
calib_summary$PA[[cvrun]],
"_",
calib_summary$run[[cvrun]])]] <-
list(weights = wt)
XGBOOST_param_list[[paste0("_",
calib_summary$PA[[cvrun]],
"_",
calib_summary$run[[cvrun]])]] <-
list(nrounds = 10000,
eta = 0.001,
max_depth = 5,
subsample = 0.75,
gamma = 0,
colsample_bytree = 0.8,
min_child_weight = 1,
weight = wt,
verbose = 0)
}
model_parameters <- bm_ModelingOptions(
data.type = "binary",
models = c("GLM", "GBM", "GAM.mgcv.gam", "MARS", "RF", "MAXNET", "XGBOOST"),
strategy = "user.defined",
user.base = "default",
user.val = list(
GLM.binary.stats.glm = GLM_param_list,
GBM.binary.gbm.gbm = GBM_param_list,
GAM.binary.mgcv.gam = GAM_param_list,
MARS.binary.earth.earth = MARS_param_list,
RF.binary.randomForest.randomForest = RF_param_list,
XGBOOST.binary.xgboost.xgboost = XGBOOST_param_list
),
bm.format = run_data,
calib.lines = table_cv
)##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Build Modeling Options -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
##
## > GLM options (datatype: binary , package: stats , function: glm )...
## > GBM options (datatype: binary , package: gbm , function: gbm )...
## > GAM options (datatype: binary , package: mgcv , function: gam )...
## > MARS options (datatype: binary , package: earth , function: earth )...
## > RF options (datatype: binary , package: randomForest , function: randomForest )...
## > MAXNET options (datatype: binary , package: maxnet , function: maxnet )...
## > XGBOOST options (datatype: binary , package: xgboost , function: xgboost )...
##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=L’étape précédente sert à définir des paramètres appropriés pour tous les modèles ; cependant, nous n’allons pas utiliser tous les modèles. Nous allons maintenant sélectionner les modèles qui seront effectivement lancés. Ce choix est basé sur les tests préliminaires qui ont mis en évidence les modèles donnant des résultats cohérents, par rapport aux modèles ne donnant pas des résultats cohérents :
les modèles GLM, GBM, MARS et MAXNET donnaient des réponses incohérentes et divergentes entre runs à la variable d’intégrité biophysique des sols
le modèle GAM donnait une représentation simplifiée de bio5 et de l’intégrité biophysique des sols, ce qui donnait une tendance à la surprédiction hors des zones à pin laricio
Au final, nous garderons donc les modèles random forests et XGBOOST.
model_runs <- BIOMOD_Modeling(
run_data,
modeling.id = "1", # ID de modélisation, on met 1 pour tous nos modèles ici
models = model_list, # Liste des modèles finaux à faire tourner
OPT.strategy = "user.defined",
OPT.user = model_parameters, # Paramètres des modèles
CV.strategy = "user.defined", # Méthode de validation croisée
CV.user.table = table_cv, # Plis générés précéemment
CV.do.full.models = FALSE,
var.import = 10, # Nombre de répétitions d'importance des variables
metric.eval = "BOYCE",
do.progress = FALSE,
nb.cpu = 16 # Nombre de coeurs à utiliser pour la modélisation
# A ajuster selon votre ordinateur, ne pas en mettre trop !
)
saveRDS(model_runs, file = "models/sittelle/model_runs.RDS")Evaluation des modèles
evals_boyce <- get_evaluations(model_runs)
ggplot(evals_boyce, aes(x = algo, y = validation)) +
geom_point(aes(col = run)) +
xlab("Algorithme") +
ylab("Indice de Boyce") +
labs(col = "Plis de\nvalidation\ncroisée") +
ylim(0, 1) +
theme_minimal()
L’indice de Boyce est un indice qui varie entre -1 et 1 (-1 = prédictions opposées à la réalité, 0 = prédiction nulles, 1 = prédictions parfaites). Ici, l’indice suggère des évaluations élevées pour tous les modèles, ce qui est encourageant : aucun modèle n’a échoué à prédire les occurrences qui n’ont pas servi à la calibration.
Il faut néamoins toujours être prudent sur l’interprétation des métriques d’évaluation car il s’agit de modèles corrélatifs et parce que l’évaluation est effectuée sur les données d’occurrence qui peuvent être biaisées. Ces métriques nous indiquent principalement qu’aucun modèle n’a donné de très mauvais résultats, c’est l’information à en retirer. En revanche, il faut se garder de la fausse impression de robustesse que peuvent donner de bonnes métriques, car les modèles peuvent faire de bonnes prédictions avec des variables qui n’ont pas de sens pour la biologie des espèces. La prochaine étape consiste donc à étudier les réponses des espèces aux variables environnementales.
Importance des variables et courbes de réponse
varimp <- get_variables_importance(model_runs)
varimp$expl.var <- reorder(varimp$expl.var,
varimp$var.imp,
median,
na.rm = TRUE)
library(dplyr)
varimp %>%
group_by(expl.var) %>%
summarise(median = median(var.imp))ggplot(varimp) +
geom_boxplot(aes(x = expl.var, y = var.imp)) +
geom_jitter(aes(x = expl.var, y = var.imp, col = algo),
alpha = .3) +
coord_flip() +
theme_minimal() +
xlab("Variable prédictive") +
ylab("Importance des variables") +
labs(col = "Algorithme")
# Variables utilisées pour la calibration
cur_vars <- model_runs@expl.var.names
# Calcul des courbes de réponse
resp <- bm_PlotResponseCurves(bm.out = model_runs,
fixed.var = "mean",
data_species = occurrences,
do.plot = FALSE,
do.progress = FALSE)$tab## No id variables; using all as measure variablescolnames(resp) <- c("Index", "Variable", "Var.value", "Model", "Response")
for (model in model_list) {
p <- ggplot(resp[grep(model, resp$Model), ], aes(x = Var.value,
y = Response)) +
geom_line(alpha = 0.2, aes(group = Model)) +
stat_smooth() +
facet_wrap(~ Variable, scales = "free_x") +
theme_bw() +
ylim(0, 1.1) +
xlab("Valeurs des variables") +
ggtitle(model)
print(p)
}## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Note: la courbe bleue est une aide pour visualiser la tendance, mais la vraie réponse des modèles correspond aux courbes grises
Comme attendu, la variable qui explique le mieux la favorabilité de l’habitat pour la sittelle Corse est la répartition des pins laricio, avec une augmentation très rapide de la favorabilité dès que le recouvrement en forêts de pins laricio est supérieur à 0% pour atteindre rapidement un plateau dès que le recouvrement atteint au moins 25% par cellule de 1km².
La seconde variable qui détermine la favorabilité de l’habitat est la température minimale annuelle (pour être précis : la valeur moyenne quotidienne de température la plus froide du mois le plus froid). Les modèles suggèrent que l’espèce occupe des habitats qui sont froids (températures minimales inférieures à 2.5°C) sans pour autant que ça soit les extrêmes de froid de la Corse (températures minimales supérieure à -7.5°C). Il est important de noter que si cette relation est détectée par les modèles, c’est qu’elle est décorrélée de la répartition des pins laricio, car l’effet de la répartition des pins laricio est déjà intégré dans les modèles.
Enfin, la sittelle occupe préférentiellement des habitats à forte intégrité biophysique, i.e. à faible degré d’exploitation humaine, ce qui correspond à notre hypothèse initiale selon laquelle la sittelle est sensible à l’exploitation forestière.
Cartes
# On ne va garder que les modèles qui ont un indice de Boyce suffisamment élevé
models_to_proj <- evals_boyce$full.name[which(evals_boyce$validation >= 0.75)]
projection_runs <- BIOMOD_Projection(
bm.mod = model_runs, # Modèles calibrés
proj.name = "corse", # Nom de la projection actuelle
new.env = env_sittelle, # Données environnementales sur lesquelles on projette les modèles
models.chosen = models_to_proj, # Modèles à projeter
build.clamping.mask = TRUE, # Le clamping mask illustre les zones où les prédictions sont en dehors des valeurs
# utilisées lors de la calibration
nb.cpu = 4)
cartes_individuelles <- rast("models/sittelle/proj_corse/proj_corse_sittelle.tif")
# Rescaling des projections qui dépassent l'intervalle 0 - 1000
cartes_individuelles[cartes_individuelles < 0] <- 0
cartes_individuelles[cartes_individuelles > 1000] <- 1000
for(i in 1:ceiling(nlyr(cartes_individuelles) / 2)) {
plot(cartes_individuelles[[(i * 2 - 1):
min(nlyr(cartes_individuelles),
(i * 2))]],
col = viridis::inferno(12))
}
Les cartes individuelles des différents modèles sont relativement convergentes dans leurs prédictions, avec quelques variations dans les zones à favorabilité intermédiaire.
Carte finale
carte_finale <- mean(cartes_individuelles)
ggplot() +
geom_spatraster(data = carte_finale) +
scale_fill_viridis(option = "inferno") +
geom_point(data = P_points[which(P_points$occurrence == 1), ],
aes(x = x, y = y),
shape = 21,
fill = "#21908CFF",
col = "white",
size = 1) +
ggtitle("Indice de favorabilité final") +
xlab("Longitude") +
ylab("Latitude") +
theme_minimal()
Comme attendu, cette carte reflète bien les données d’occurrence, illustrant les zones à plus forte favorabilité localisées au niveau des forêts de pin laricio en altitude. On observe des nuances de favorabilité au sein des zones de pin laricio, qui illustrent les préférences en termes de température et de degré d’exploitation forestière. On peut noter certaines occurrences dans des zones à faible favorabilité à l’ouest, au nord-est et au sud de l’île, qui peuvent possiblement refléter des observations d’individus errants, des occurrences dans des zones historiquement occupées par les forêts de pin laricio, ou de potentielles limites des modèles dus à d’éventuels déséquilibres dans les données d’occurrences. Cette carte a été jugée tout à fait cohérente par les experts consultés.
carte_incertitude <- app(cartes_individuelles, sd)
ggplot() +
geom_spatraster(data = carte_incertitude) +
scale_fill_continuous(type = "viridis") +
ggtitle("Incertitude\n(écart-type des probabilités)") +
theme_minimal()
La carte d’incertitude – qui correspond à l’écart-type des prédictions entre modèles – reflète en partie ces occurrences difficile à prédire, avec des incertitudes élevées dans les zones concernées. Cette carte d’incertitude illustre également des difficultés à prédire avec précision les contours des zones les plus favorables à la sittelle, comme illustré par les zones à incertitude intermédiaire au centre de l’île.
Carte de potentiel d’habitat
Pour créer la carte de potentiel d’habitat final, nous allons représenter trois catégories de potentiel d’habitat, en respectant les contraintes d’interprétation sur les modèles en présence seule. En effet, les modèles en présence seule ne peuvent pas fournir d’information sur la probabilité de présence. Par conséquent, ils ne peuvent informer sur les habitats défavorables - ils informent seulement sur les habitats favorables compte-tenu des connaissances actuelles.
Ainsi, nous ne produirons pas de carte binaire “présence-absence” qui n’aurait pas de sens dans le cadre des modèles en présence-seule et qui est également une sur-simplification de la réalité biologique, qui n’est jamais binaire. Nous allons plutôt représenter trois catégories :
- les zones à fort potentiel d’habitat
- les zones à potentiel d’habitat intermédiaire
- les zones à potentiel d’habitat faible ou méconnu
Pour établir une méthode permettant de définir ces trois catégories, on peut étudier comment les occurrences sont réparties sur le gradient de favorabilité des modèles. On peut alors utiliser les quantiles des occurrences pour identifier les seuils séparant les catégories.
favorabilite_presences <- extract(carte_finale,
P_points[which(P_points$occurrence == 1),
c("x", "y")],
ID = FALSE)
qt_favorabilite <- quantile(favorabilite_presences$mean, probs = c(.05, .25))
ggplot(favorabilite_presences) +
geom_boxplot(aes(x = mean),
col = "darkgrey") +
geom_vline(xintercept = qt_favorabilite,
col = c("#1b9e77", "#7570b3"),
linetype = 2,
linewidth = 2) +
theme_minimal() +
xlab("Indice de favorabilité") +
scale_y_continuous(breaks = 0,
labels = "Occurrences") +
xlim(0, 1000)
Dans le graphe ci-dessus, on voit la répartition des occurrences sur
l’indice de favorabilité produit par le modèle. On peut utiliser les
quantiles à 5% et 25% (représentés par les pointillés bleus) pour
séparer les catégories.
La zone à droite du quantile à 25% (le trait mauve) contient l’essentiel des occurrences du groupe d’espèces, ce qui signifie qu’au delà de ce seuil, le potentiel d’habitat est élevé.
La zone entre le quantile à 5% (trait vert) et à 25% (trait mauve) est une zone à favorabilité plus faible mais qui contient tout de même 20% des occurrences du groupe. On peut ainsi la caractériser comme zone à potentiel d’habitat intermédiaire.
La zone à gauche du quantile à 5% (trait vert) contient moins de 5% des occurrences du groupe. Il s’agit donc de valeurs de favorabilité plutôt faibles puisqu’elles ne semblent pas ou peu occupées d’après les connaissances actuelles. On peut donc qualifier cette catégorie de potentiel d’habitat faible ou méconnu.
Si l’on utilise ces seuils pour illustrer la répartition de ces trois catégories, on obtient la carte suivante :
carte_indice <- carte_finale
carte_indice[carte_finale < qt_favorabilite["5%"]] <- 0
carte_indice[carte_finale >= qt_favorabilite["5%"] &
carte_finale < qt_favorabilite["25%"]] <- 1
carte_indice[carte_finale >= qt_favorabilite["25%"]] <- 2
carte_indice <- as.factor(carte_indice)
ggplot() +
geom_spatraster(data = carte_indice) +
scale_fill_manual(values = viridis::plasma(3),
name = paste0("Potentiel d'habitat\n(% du total ",
"d'occurrences\n",
"observé dans cette classe\n",
"de favorabilité)"),
labels = c("Faible ou méconnu (< 5%)",
"Intermédiaire (5-25%)",
"Elevé (75%)"),
na.translate = F)
La plupart des occurrences de la sittelle Corse sont localisées dans les zones à très fort indice de favorabilité (75% des occurrences sont localisées dans les zones à favorabilité supérieure à 800), on peut donc qualifier cette classe de zone à potentiel d’habitat élevé. Ensuite, les zones à plus faible indice de favorabilité mais qui contiennent tout de même 20% des occurrences – et qui peuvent donc être qualifiées de zone à potentiel d’habitat intermédiaire – sont localisées sur tout le pourtour de la zone à très fort potentiel d’habitat, ainsi que dans quelques patchs isolés au nord-est et au sud-ouest.