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 `Chiro_forestier' from data source
## `C:\Rprojects\SDMs_PNA_Corse\data\donnees_brutes\taxa\Chiro_forestier.shp'
## using driver `ESRI Shapefile'
## replacing null geometries with empty geometries
## Simple feature collection with 15400 features and 69 fields (with 27 geometries empty)
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 8.544579 ymin: 41.37692 xmax: 9.546355 ymax: 43.02367
## Geodetic CRS: WGS 84# Simplification du nom d'espèce en binomial
chiroforet$species <- simplify_species_name(chiroforet$nom_valide)
# Retrait d'une araignée qui s'est glissée par erreur dans l'extraction des
# données du groupe des chiroptères forestiers
chiroforet <- chiroforet[-which(chiroforet$species %in%
c("Enoplognatha thoracica")), ]
# Dates d'échantillonnage
chiroforet$year <- as.numeric(strtrim(chiroforet$date_fin, 4))
chiroforet$month <- as.numeric(substr(chiroforet$date_fin, 6, 7))
# Visualisation de la temporalité des occurrences
ggplot(chiroforet) +
geom_boxplot(aes(x = species,
y = year))+
coord_flip() +
scale_y_continuous(breaks = breaks_pretty()) +
theme_minimal()
Filtre temporel
Les espèces présentent une saisonnalité forte ; nous ne modéliserons ici que les habitats de la période estivale pour laquelle il existe suffisamment de données. Le Groupe Chiroptères Corses suggère une temporalité allant du 15 mai au 15 septembre ; nous ne retiendrons donc que les données des mois de mai à septembre.
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.
Le nombre de données supprimées en fixant un seuil à l’année 2000 est très faible :
# Les données avant 2000 représentent un % modéré du jeu de données :
100 * length(which(chiroforet$year < 2000)) / nrow(chiroforet)## [1] 3.967351La couverture temporelle sur l’année est plus complète après 2000 :
ggplot(chiroforet) +
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 ne change pas de manière majeure avec ou sans les données pré-2000 :
p_chiroforet_all <- ggplot() +
geom_sf(data = corse) +
geom_sf(data = chiroforet, aes(col = year)) +
scale_color_continuous(type = "viridis") +
theme_minimal(base_size = 15) +
ggtitle("Toutes données\nchiroptères forestiers")
p_chiroforet_post2000 <- ggplot() +
geom_sf(data = corse) +
geom_sf(data = chiroforet[chiroforet$year >= 2000, ], aes(col = year)) +
scale_color_continuous(type = "viridis") +
theme_minimal(base_size = 15) +
ggtitle("Données post-2000\nchiroptères forestiers")
ggarrange(p_chiroforet_all,
p_chiroforet_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 des espèces.
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
chiroforet_r <- rasterize(chiroforet,
env_corse)
names(chiroforet_r) <- "chiroforet"
plot(chiroforet_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 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_chiroforet <- c(env_corse,
chiroforet_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_chiroforet_df <- values(env_chiroforet)
# On regarde le nombre d'occurrences pour lesquelles il y a des données
# manquantes :
length(which(is.na(env_chiroforet_df[, "bio1"]) &
!is.na(env_chiroforet_df[, "chiroforet"])))## [1] 9On 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_chiroforet_df[, 1])), ]
# Et ensuite sur le tableau avec variables env et présences d'espèces
env_chiroforet_df <- env_chiroforet_df[which(!is.na(env_chiroforet_df[, 1])), ]
# Comparaison du nombre d'occurrences :
# Avant rasterisation
nrow(chiroforet)## [1] 14354# Après rasterisation et élimination des données env manquantes
length(which(env_chiroforet_df[, "chiroforet"] == 1))## [1] 623Il s’agit donc du nombre d’occurrences que l’on va pouvoir utiliser
pour calibrer nos 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_chiroforet_df[, "chiroforet"])), ],
# Ensuite, on récupère la colonne qui indique présence pour chaque cellule
occurrence = env_chiroforet_df[which(!is.na(env_chiroforet_df[, "chiroforet"])),
"chiroforet"])
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_amphib_df
# Par conséquent, son nombre de lignes est égal au nombre total de pixels
# disponibles sur la Corse
nrow(env_chiroforet_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.
Il existe un biais d’échantillonnage sur les occurrences de chiroptères, qui peut influencer les résultats des modèles. Pour corriger ce biais, deux options sont possibles : prendre en compte l’accessibilité (distance aux routes) ou prendre en compte la pression d’échantillonnage sur tout le groupe. C’est cette méthode que nous allons utiliser ici, en intégrant le biais d’échantillonnage dans la génération des points de background. Cela va ainsi éviter aux modèles d’expliquer le biais d’échantillonnage avec les variables environnementales.
Pour cela, il faut générer une couche qui reflète la probabilité d’échantillonnage. Nous avons généré cette couche en utilisant toutes les occurrences de chiroptères de Corse. Il s’agit de la méthode appelée “biais du groupe cible” (target-group sampling bias). Nous utilisons ensuite cette couche pour échantillonner les points de background en Corse.
Etant donné que la Corse ne possède qu’un nombre limité de cellules, nous réduisons le nombre de backgrounds tirés à 6000 afin que le biais produise un effet (sinon quasiment toutes les cellules de la Corse sont échantillonnées, ce qui supprime l’effet du biais).
# On réduit également le nombre de background pour avoir un effet du biais
background <- spatSample(env_corse[["occ_density_chiro"]],
method = "weights",
size = 6000,
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
Les chauves-souris sont dépendantes des conditions climatologiques pour la sélection de leurs abris et sites de nichage et pour leur activité. Les études macroécologiques suggèrent que les chiroptères préfèrent des températures douces à chaudes et un degré d’humidité élevé (McCain 2007), et les variables de température et de précipitation sont systématiquement utilisées comme des prédicteurs pertinents dans les modèles de distribution de chiroptères. Nous utiliserons donc des variables reflétant les limites potentielles qui empêchent l’occurrence des chiroptères, avec l’hypothèse que les espèces préfèrent des conditions intermédiaires, douces et humides : effet limitant des températures trop froides ou trop chaudes (températures les plus chaudes et les plus froides de l’année), effet limitant de l’humidité relative (humidité relative minimale et maximale).
Noms des variables retenues :
températures les plus chaudes de l’année (bio5)
températures les plus froides sur la saison estivale (tasmin_chiro)
humidité relative minimale (cmi_min) et maximale (cmi_max).
Occupation du sol
Les chiroptères forestiers vivent et nichent dans les forêts en utilisant notamment les anfractuosités du bois mort. Ainsi, s’il est attendu que ces espèces dépendent des forêts, il est possible que toutes les forêts ne permettent pas de leur offrir un habitat réellement favorable (Miller et al., 2003). Par exemple, les forêts gérées dans lesquelles les arbres morts ne sont pas laissés sur pied ou les forêts de plantation monospécifiques avec une rotation de coupe relativement courte ne permettent probablement pas d’avoir les conditions de microhabitat requises pour la plupart des chiroptères forestiers. Ainsi, il apparait nécessaire d’effectuer une distinction dans le type de couvert forestier pour identifier plus clairement les habitats favorables aux chiroptères forestiers. Cependant, il est à attendre un fort degré d’incertitude qui peut rendre la relation moins nette ou moins évidente, à la fois à cause de l’incertitude sur la nature de la gestion des forêts et sur le décalage temporel entre moment d’observation de l’espèce et date des couches géographiques d’occupation du sol. Nous utiliserons l’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).
Ensuite, la proximité aux zones humides a également été démontrée comme ayant un effet positif sur de nombreuses espèces de chiroptères, car les zones humides fournissent des ressources conséquentes pour l’alimentation (Sirami et al., 2013). Nous testerons donc un effet positif de la distance aux zones humides sur la probabilité d’observer les chiroptères. Enfin, il a été démontré un effet négatif de l’homogénéisation du territoire sur la probabilité d’observation des chiroptères (Put et al., 2019). Nous testerons donc si l’augmentation de la diversité d’occupation augmente la probabilité d’observer les chiroptères, en utilisant par exemple des indices de Simpson ou Shannon sur le recouvrement des différentes classes d’occupation du sol. Nous étudierons également l’effet de la distance aux cavités, car celles-ci peuvent également abriter les différentes espèces de ce groupe.
Enfin, la connectivité peut beaucoup jouer sur la présence des chauves-souris en leur permettant de se déplacer ; nous utiliserons donc la variable de continuité spatiale développée par l’IUCN (Guetté et al. 2021) en supposant un effet positif de la connectivité sur la probabilité d’observer les espèces.
Noms des variables retenues :
couvert forestier (forets)
distance aux forêts (dist_forets)
intégrité biophysique des sols (integrite)
distance aux zones humides (dist_moy_pzh)
homogénéisation du paysage (simpson_landscapediv)
connectivité (connectivite)
Biais d’échantillonnage
Pour prendre en compte le biais potentiel d’échantillonnage, nous utiliserons la probabilité d’échantillonnage basée sur la méthode dite du “biais d’échantillonnage du groupe cible”. Cette méthode consiste à utiliser toutes les occurrences du groupe cible dans la zone d’étude pour spatialiser le biais d’échantillonnage, en utilisant une fonction kernel bi-dimensionnelle. Nous avons construit cette cette variable de biais d’échantillonnage en utilisant les occurrences de tous les chiroptères sur toute la Corse. Elle sera utilisée pour simuler les points de background avec le même biais d’échantillonnage.
Variables anthropogéniques
La pollution lumineuse affecte de manière différente les chauves-souris : certaines espèces sont très négativement impactées, tandis que d’autres y semblent insensibles voire en bénéficient (Azam et al., 2016). Nous testerons l’effet de la pollution avec un a priori négatif pour ce groupe qui ne semble pas inféodé aux habitats artificiels. Néanmoins, certaines des espèces de ce groupe sont connues pour occuper des habitats artificiels (e.g., Barbastella barbastellus), ce qui peut créer de l’incertitude sur l’effet de la pollution lumineuse et générer un effet favorable de l’emprise urbaine, que nous testerons également.
Noms des variables retenues :
pollution lumineuse (pollum)
zones artificielles (zones_artif)
Autres variables et commentaires
Les chauves-souris régissent leur activité en fonction du vent, préférant les périodes où le vent est faible (Barré et al. 2023). Ainsi, il est attendu qu’elles aient moins de chance d’occuper les zones fortement ventées ou sujettes aux rafales excessives. Nous testerons donc s’il existe un effet négatif des variables correspondant aux vitesses de vent.
Noms des variables retenues :
- Vitesse du vent mensualisée maximale sur l’année (sfcWind_max)
Constitution du jeu de variables finales pour les chiroptères forestiers
Préparation des rasters
Etude de la colinéarité et réduction du nombre de variables
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.
# Modifier pour les corrélations 2 à 2 ? e.g. corrplot?
var_groups <- removeCollinearity(env_chiroforet,
plot = TRUE,
multicollinearity.cutoff = 0.7,
method = "spearman")
Deux jeux de variables sont corrélées fortement : intégrité et connectivité ; couvert forestier et distance aux forêts. Nous conserverons ici l’intégrité plutôt que la connectivité, et le couvert forestier plutôt que la distance aux forêts.
env_chiroforet <- env_chiroforet[[-which(names(env_chiroforet) ==
"connectivite")]]
env_chiroforet <- env_chiroforet[[-which(names(env_chiroforet) ==
"forets")]]Au total cela fait 11 variables environnementales, ce qui peut être excessif pour les modèles et créer un risque d’overfitting. Nous allons donc réduire le nombre de variables à 10, en supprimant les zones artificielles car l’hypothèse liée à cette variable est la plus faible, et les effets des zones artificielles sont probablement inclus dans la pollution lumineuse et l’intégrité biophysique des sols.
Filtrage des variables non informatives
Les analyses préliminaires ont révélé que certaines de ces variables n’ont des effets que très limités : integrite, cmi_min, simpson_landscapediv. Ces variables seront donc éliminées.
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_chiroforet_l93 <- project(env_chiroforet,
"EPSG:2154") # Projection en Lambert 93 ici
# Ensuite on étudie le range d'autocorrélation spatiale
AC_range <- cv_spatial_autocor(env_chiroforet_l93,
num_sample = 10000)##
|
| | 0%
|
|=========== | 14%
|
|======================= | 29%
|
|================================== | 43%
|
|============================================= | 57%
|
|======================================================== | 71%
|
|==================================================================== | 86%
|
|===============================================================================| 100%
On obtient initialement un range médian qui est de 18 km. Cependant, les données de chiroptères forestiers présentent une répartition aggrégée qui cause des déséquilibres dans les plis de validation croisée. Pour corriger ce problème nous réduirons la taille des blocs à 10km².
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 = 10000, # Taille des blocs en metres
selection = "random", # Attribution des blocs aléatoire dans
# les plis
iteration = 500, # Nombre d'essais pour trouver des plis
# équilibrés
biomod2 = TRUE, # Formater les données pour biomod2
r = env_chiroforet, # Pour le fond de carte
progress = FALSE,
plot = FALSE) ##
## train_0 train_1 test_0 test_1
## 1 4984 495 1016 128
## 2 4998 525 1002 98
## 3 4992 521 1008 102
## 4 4998 529 1002 94
## 5 5006 525 994 98
## 6 5022 520 978 103On voit que nos plis sont plutôt équilibrés :
de 495 à 529 présences pour la calibration
de 94 à 128 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/chiroforet", recursive = T, showWarnings = FALSE)
run_data <- BIOMOD_FormatingData(
resp.name = "chiroforet", # Nom de l'espèce
resp.var = occurrences, # Présences + background
expl.var = env_chiroforet, # 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##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= chiroforet 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 623 présences et 6000 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,
method = "REML")
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, MARS et MAXNET donnaient des réponses biologiquement irréalistes en U
le modèle GAM donnait des réponses linéaires monotones pour toutes les variables, alors que la plupart des autres modèles donnaient des réponses non linéaires
Au final, trois modèles ont donc été retenus pour ce groupe : RF, GBM 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/chiroforet/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)
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") +
ylab("Réponse") +
ggtitle(model)
print(p)
}## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Le prédicteur le plus important de la répartition des chiroptères forestiers est la température estivale minimale, ce groupe préférant des zones plus fraiches avec un effet de seuil négatif dans les zones chaudes durant toute la saison estivale. On note néanmoins un pic de favorabilité aux températures les plus chaudes également, donnant une réponse bimodale qui peut paraître incohérente d’un point de vue biologique. Après échange avec les experts suite à la réunion du 13 février 2024, nous avons pu poser une hypothèse pour expliquer cette réponse bimodale. Les données d’observations de ce groupe contiennent deux types de données : en majorité, des données correspondant aux habitats de gîte (captures) et en minorité, des données correspondant aux habitats de chasse (radiotracking, ultrasons). Ainsi, le premier pic (températures fraiches) correspond probablement aux habitats de gîte en montagne, tandis que le second pic (températures plus chaudes) correspond probablement aux habitats de chasse en plaine. Les autres facteurs climatiques ont un effet moins important, avec une préférence pour les zones qui atteindront des températures élevées au plus fort de l’été, et les zones où l’humidité n’est pas excessive en période de pluie excédentaire.
Le second facteur le plus important est la distance aux forêts, ce qui est cohérent avec nos attentes initiales pour ce groupe. On observe une favorabilité maximale à distance très faible, qui décroit rapidement jusqu’à une distance de 2km où la favorabilité est plus faible. On note de très légers rebonds à distance plus grande (4 km) qui sont probablement dus aux aires de nourrissage et/ou à des patchs forestiers non détectés dans notre variable sur les forêts.
Le troisième facteur le plus important est la distance aux zones humides, avec un effet conforme à nos attentes de préférence pour la proximité aux zones humides, utilisées comme zones de nourrissage.
La pollution lumineuse joue également un rôle important avec un effet négatif dès la présence de pollution lumineuse, même en faible quantité. On observe de légers rebonds de la courbe à pollution lumineuse intermédiaire, qui pourraient potentiellement s’expliquer par les habitats de chasse de ce groupe. Cet effet de la pollution lumineuse traduit, comme attendu, une sensibilité particulière de ce groupe d’espèces à ce type de perturbation anthropique. Ces résultats suggèrent qu’un type d’action PNA à envisager est de limiter l’exposition à la pollution lumineuse, surtout dans les habitats favorables à ce groupe. Enfin, on note une légère préférence des chiroptères de forêts pour les zones qui ne sont pas sujettes à des vents excessifs.
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_chiroforet, # 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/chiroforet/proj_corse/proj_corse_chiroforet.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))
}
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()
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 de l’indice de favorabilité final représente la moyenne entre les modèles qui ont satisfait les critères de qualité (i.e., indice de Boyce supérieur à 0.75 et réponses cohérentes). Elle illustre le gradient de favorabilité de l’habitat sur toute la Corse (Figure 11), soulignant essentiellement les zones proches des forêts dans des altitudes intermédiaires. On peut noter deux limites sur cette carte : tout d’abord, les valeurs de favorabilité sont généralement intermédiaires, il n’y a pas beaucoup de zones avec des valeurs très élevées (i.e., indice > 800). Ensuite, on note qu’un certain nombre d’occurrences sont localisées dans des zones à favorabilité relativement basse, notamment en basse altitude. On peut supposer que ces deux limites sont dues au mélange des deux types d’observations qui reflètent des occupations de l’habitat différentes (gîte vs site de chasse).
La carte d’incertitude correspond à l’écart-type des indices des modèles ayant passé les critères de qualité. Elle complète l’analyse de la carte de l’indice de favorabilité finale. On note un certain nombre de cellules isolées où l’incertitude est élevée (points vert clair à jaune sur la Figure 12). Ces points reflètent des phénomènes ponctuels de surajustement de certains modèles (notamment, RF et XGBOOST) au niveau de certaines occurrences, ce qui illustre les difficultés des modèles à obtenir des réponses aux variables environnementales qui englobent toutes les occurrences. Ensuite, on note des zones plus étendues où l’incertitude est intermédiaire (zones couleur cyan). Ces zones illustrent probablement deux types d’incertitudes : d’une part, l’incertitude habituelle sur les zones à favorabilité intermédiaire, qui sont toujours difficiles à bien définir par les modèles, et d’autre part, l’incertitude sur les zones correspondant aux observations de lieux de chasse, minoritaires dans les données, et donc que les modèles ont du mal à reproduire.
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)
Cette carte illustre bien la zone qui présente le potentiel d’habitat le plus élevé pour le groupe des chiroptères forestiers, essentiellement située dans le centre et l’ouest de la Corse. Les zones à potentiel d’habitat intermédiaire sont essentiellement situées en marge des zones à fort potentiel, illustrant un gradient de favorabilité progressif, ainsi que dans les zones de basses altitude plus proches des côtes, et notamment à l’est et au nord-est de la Corse qui reflètent probablement des zones de chasse, suivant l’hypothèse émise par les experts du Groupe Chiroptères Corse.