TP3 Analyse de réseaux biogéographiques

Avant de commencer le TP, je vais vous présenter une introduction à la biorégionalisation, ainsi que le contenu du TP. Pour aller plus loin, ce TP et le package bioregion s’appuient sur les travaux présentés dans Lenormand et al. (2019).

 

3.1 Prise en main des données

Le package R bioregion met à disposition des données d’exemple d’abondance d’espèces végétales dans le sud de la France, fournies par le conservatoire botanique méditerranéen.

• vegemat : matrice de co-occurrences sites-espèces.
  
• vegesf : données spatiales associées (nous aurons besoin de sf pour les manipuler).

Nous allons dans un premier temps calculer quelques statistiques de base pour explorer ces données.

Créez un script TP3.R et enregistrez-le dans votre répertoire de travail, que vous pouvez également nommer TP3.

Dans ce TP, je m’inspire du cours sur la bioregionalisation avec les réseaux proposé par mon collègue Boris Leroy du Muséum national d’Histoire naturelle.

3.1.1 Charger les données

# Charger les packages
library(bioregion)
library(TDLM)
library(sf)

# Définir le répertoire de travail
setwd("~/Desktop/TP3") # A adapter à votre chemin sur votre ordinateur !

# Charger les données
data(vegemat)
data(vegesf)

3.1.2 Taille et aperçu des objets

dim(vegemat)         # Dimensions et aperçu de la matrice sites-espèces
vegemat[1:10,1:10]

vegesf               # Aperçu des entités spatiales 
class(vegesf)        # Objet sf
plot(vegesf)         # Représentation spatiale de vegesf

3.1.3 Les réseaux en biogéographie

Dans le TP précédent, nous avons travaillé sur un réseau dirigé et pondéré, où chaque lien représentait un flux entre deux communes (origine → destination). En biogéographie, nous allons maintenant manipuler un réseau bipartite.

Un réseau bipartite relie deux ensembles de noeuds distincts, ici les sites (zones géographiques) et les espèces (entités biologiques). Autrement dit, un lien existe uniquement entre un site et une espèce, jamais entre deux sites ou deux espèces directement. Ce type de réseau peut être pondéré (abondance d’espèces) ou non pondéré (présence/absence).

Le réseau bipartite sites-espèces que nous considérons dans ce TP est pondéré et stocké sous forme matricielle dans l’objet vegemat. Nous pouvons ici aussi facilement convertir cette matrice en un format “trois colonnes” à l’aide de la fonction mat_to_net().

# Conversion matrice -> "réseau"
vegenet <- mat_to_net(vegemat, weight = TRUE)
vegenet[1:10,]

On peut également considérer un réseau non pondéré.

# Suppression de la colonne poids
vegenet_bin <- vegenet[,-3]
vegenet_bin[1:10,]

Et le reconvertir en matrice avec net_to_mat().

# Conversion "réseau" -> matrice  
vegemat_bin <- net_to_mat(vegenet_bin)
vegemat_bin[1:10,1:10]

3.1.4 Exercice

Tracer les histogrammes du nombre moyen d’espèces par site et du nombre moyen de sites par espèce à partir de vegemat_bin (à l’aide de la fonction apply()).

---

Solution (cliquer pour afficher)

# Histogramme du nombre d'espèces par site
species_per_site <- apply(vegemat_bin,1,sum)
hist(species_per_site, col = "steelblue3")

# Histogramme du nombre d'espèces par site
sites_per_species <- apply(vegemat_bin,2,sum)
hist(sites_per_species, col = "#CC6666")

---

3.2 Similarité entre sites

Pour comparer les sites entre eux, nous pouvons calculer un indice de similarité basé sur la composition en espèces. Par exemple, l’indice de Simpson mesure la proportion d’espèces partagées entre deux sites, en donnant un poids plus fort aux espèces présentes dans les deux zones. Il est robuste aux différences de richesse spécifique.

sim <- similarity(vegemat, metric = "Simpson")
sim[1:10,]

Pour mieux comprendre la manière dont cet indice est construit, il est possible de récupérer les trois composantes fondamentales de la comparaison entre deux sites :

• a : nombre d’espèces partagées (intersection)
  
• b : espèces uniques au site 1
  
• c : espèces uniques au site 2

Ces quantités peuvent être obtenues via :

simabc <- similarity(vegemat, metric = "abc")
simabc[1:10,]

3.2.1 Exercices

1. Tracer l’histogramme des similarités.

---

Solution (cliquer pour afficher)

# Histogramme du nombre d'espèces par site
hist(sim$Simpson, col = "steelblue3")

---

2. A partir de la formule de Simpson que vous trouverez ici et de la fonction pmin(), reconstruisez l’indice de Simpson à partir des valeurs
a, b et c de simabc.

---

Solution (cliquer pour afficher)

# Extraire a, b et c
a <- simabc$a
b <- simabc$b
c <- simabc$c

# Reconstruction de Simpson
simpson <- 1 - pmin(b,c) / (a + pmin(b, c))

# Aperçu
sim$Simpson2 <- simpson
sim[1:10,]

---

3.3 Bioregionalisation

La biorégionalisation consiste à identifier des groupes de sites partageant des caractéristiques écologiques ou biogéographiques similaires. Nous nous baserons ici sur le réseau de similarité calculé avec l’indice de Simpson ci-dessous.

Pour cela nous pouvons appliquer un algorithme de détection de communautés, ici Louvain afin d’obtenir des régions biogéographiques cohérentes.

# Détection des communautés biogéographiques
clu <- netclu_louvain(sim, cut_weight = 0)

# Effectifs des biorégions obtenues
table(clu$clusters[,2])

# Visualisation cartographique
map_bioregions(clu, vegesf)

3.3.1 Exercice

Le paramètre cut_weight permet de retirer du réseau les liens (similarités) dont la valeur est inférieure à un seuil donné. Plus ce seuil est élevé, plus le réseau est épuré et plus les communautés détectées sont contrastées (mais potentiellement plus nombreuses). Il est important d’expérimenter plusieurs valeurs de cut_weight car cela peut fortement influencer la structure du réseau et donc la partition obtenue.

Tester différentes valeurs de cut_weight (0.5, 0.7, 0.8, 0.9) puis visualiser la biorégionalisation obtenue.

---

Solution (cliquer pour afficher)

# Détection des communautés biogéographiques avec le seuil 0.7
clu <- netclu_louvain(sim, cut_weight = 0.7, seed = 1)

# Effectifs des biorégions obtenues
table(clu$clusters[,2])

# Visualisation cartographique
map_bioregions(clu, vegesf)

---

3.4 Analyse biogéographique

3.4.1 Contributions des espèces aux biorégions

Nous allons conserver la bioregionalisation au seuil 0.7 qui contient 5 bioregions en fixant des couleurs pour chaque biorégion. L’objet sf bioreg ci-dessous contient une variable K_5 dans sa table attributaire qui permet d’identifier la bioregion assignée à chacun des sites.

# Détection des communautés biogéographiques avec le seuil 0.7
clu <- netclu_louvain(sim, cut_weight = 0.7, seed = 1)

# Ajouter des couleurs
clu <- bioregion_colors(clu, palette = "Pastel")

# Bioregions
bioreg <- map_bioregions(clu, 
                         vegesf, 
                         write_clusters = TRUE,
                         plot = TRUE)

Nous allons maintenant essayer de comprendre quelles sont les espèces qui contribuent le plus aux bioregions. Pour cela, nous allons utiliser la fonction site_species_metrics() qui va nous permettre, avec le code ci-dessous, de calculer le degré de contribution de chaque espèce à chaque bioregion.

contrib <- site_species_metrics(clu,
                                bioregion_metrics = "Rho",
                                bioregionalization_metrics = NULL,
                                comat = vegemat,
                                verbose = FALSE)

contrib <- net_to_mat(contrib$species_bioregions, weight = TRUE)

contrib[1:10,]

Les valeurs de contrib sont positives lorsqu’une espèce contribue positivement à une bioregion. Cela signifie qu’on la retrouve plus fréquemment dans cette bioregion que ce que l’on pourrait attendre dans le cas d’une distribution uniforme. A l’inverse, des valeurs négatives indiquent qu’une espèce est sous-représentée dans la bioregion considérée, c’est-à-dire qu’elle y est moins présente que ce que l’on attendrait dans le cas d’une distribution uniforme.

3.4.2 Exercices

1. Utilisez la fonction order() pour identifier les espèces contribuant le plus à la première bioregion.

---

Solution (cliquer pour afficher)

contrib1 <- contrib[order(contrib[,1], decreasing = TRUE),]

contrib1[1:10,]

---

2. Existe-t-il des espèces qui contribuent positivement à la fois à la biorégion 1 à la biorégion 2 ?

---

Solution (cliquer pour afficher)

contrib12 <- contrib[contrib[,1] > 0 & contrib[,2] > 0,]

contrib12[1:10,]

---

3.4.3 Relations entre biorégions

Pour analyser les relations entre bioregions, nous allons filtrer dans contrib les espèces ayant une contribution positive et significative (supérieure à 1.96).

Ensuite, nous normalisons ces valeurs afin que les contributions de chaque espèce somment à 1. Cette étape transforme les contributions en proportions comparables entre espèces.

# Copier la table des contributions
contrib_196 <- contrib

# Filtrer uniquement les contributions significatives (> 1.96)
contrib_196[contrib_196 < 1.96] <- 0

# Normaliser les contributions par ligne (chaque espèce)
contrib_196 <- contrib_196 / apply(contrib_196, 1, sum)
contrib_196[is.na(contrib_196)] <- 0

contrib_196[1:10,]

À partir de cette table, nous pouvons calculer la matrice lambda, qui résume les relations entre bioregions.

K <- ncol(contrib_196)

lambda <- matrix(0, K, K)
colnames(lambda) <- colnames(contrib_196)
rownames(lambda) <- colnames(contrib_196)

for(k in 1:K){
  
  contrib_196_k <- contrib_196[contrib_196[, k] > 0,]
  lambda[k,] <- apply(contrib_196_k, 2, mean)
  
}

lambda

Nous allons maintenant visualiser le réseau biogéographique lambda avec le package visNetwork (Almende & Thieurmel, 2025).

Nous allons commencer par charger le package, mettre lambda au format réseau et conserver uniquement les liens supérieurs 0.01.

library(visNetwork)

lambda_net <- mat_to_net(lambda, weight = TRUE)
lambda_net <- lambda_net[lambda_net$Weight > 0.01,]

Nous pouvons maintenant visualiser le réseau dirigé pondéré. La taille des noeuds correspond au degré de spécialisation des biorégions (i.e. la diagonale de lambda), tandis que la direction et l’épaisseur des liens représentent respectivement le sens et l’intensité des contributions entre biorégions.

Le code est relativement long, mais reste simple à comprendre. N’hésitez pas à le modifier afin d’ajuster la taille des nœuds, l’épaisseur des liens ou le comportement dynamique du réseau (physique, espacement des nœuds, stabilisation).

# Noeuds pondérés
nodes <- lambda_net[lambda_net$Node1 == lambda_net$Node2, -2]

# Taille des noeuds
node_size <- nodes$Weight

# Identifiant des noeuds
node_id <- nodes$Node1

# Nodes dans visNetwork
nodes <- data.frame(
  id = node_id,
  label = paste("BR", node_id),
  value = node_size * 50       
)

# Couleur des noeuds
cols <- clu$colors$K_5$color[match(node_id,clu$colors$K_5$cluster)]

nodes$color.background <- cols
nodes$color.border <- cols

# Liens dirigées et pondérés
edges_data <- lambda_net[lambda_net$Node1 != lambda_net$Node2, ]

# Edges dans visNetwork
edges <- data.frame(
  from = edges_data$Node1,
  to = edges_data$Node2,
  width = edges_data$Weight * 20,
  arrows = "to",
  smooth = ifelse(edges_data$Node1 < edges_data$Node2, "curvedCW", "curvedCCW"),
  length = 100   # Distance idéale entre nœuds
)

# Couleur des liens
edges$color <- cols[match(edges$from, nodes$id)]

# Visualisation du réseau
visNetwork(nodes, edges) %>%
  visEdges(smooth = TRUE) %>%
  visNodes(shadow = TRUE) %>%
  visPhysics(
    solver = "forceAtlas2Based",
    forceAtlas2Based = list(
      gravitationalConstant = -50,
      centralGravity = 0.01,
      springLength = 100,
      springConstant = 0.01,
      damping = 0.4,
      avoidOverlap = 1
    ),
    stabilization = list(iterations = 200)
  )

3.4.4 Exercices

1. Quelle est la biorégion la plus spécialisée (i.e. qui possède des espèces qui contribuent peu aux autres biorégions) ?

2. Quelle est la biorégion qui contribue le plus aux autres biorégions ?

3.5 Echelle communale

Dans cette section, nous allons travailler à l’échelle communale pour analyser la répartition des bioregions au sein des communes. Nous disposons de nos deux jeux de données : mtp et bioreg.

Nous allons commencer par charger mtp et modifier sa projection pour qu’elle corresponde à celle de bioreg.

data(mtp)
mtp_l <- st_transform(mtp, crs = st_crs(bioreg))

Assurons nous que tout est ok en visualisant les deux couches pour vérifier leur superposition.

plot(st_geometry(bioreg), border = "#CC6666")
plot(st_geometry(mtp_l), border = "steelblue3", add = TRUE)

3.5.1 Exercice

Calculer la fraction de surface de chaque polygone de mtp_l correspondant à chaque bioregion définie dans bioreg.

1. Réaliser l’intersection spatiale entre mtp_l et bioreg.
2. Calculer l’aire de chaque intersection.
3. Calculer la fraction de surface de chaque polygone correspondant à chaque bioregion (aire intersection / aire totale du polygone).
4. Agréger les fractions si un polygone est découpé en plusieurs parties pour la même bioregion.
5. Organiser le résultat sous forme de matrice contenant les fractions de surface avec en ligne les communes et en colonnes les bioregions.

---

Solution (cliquer pour afficher)

# -----------------------------------------------
# 1. Intersection spatiale
# -----------------------------------------------
# On calcule la partie de chaque polygone mtp_l qui tombe dans chaque bioregion
# st_intersection conserve les attributs de mtp_l et de bioreg$K_5
inter <- st_intersection(mtp_l, bioreg[, "K_5"])

# -----------------------------------------------
# 2. Calcul de l'aire des intersections
# -----------------------------------------------
# st_area retourne la surface (en unités du CRS) de chaque polygone d'intersection
inter$area_inter <- st_area(inter)

# -----------------------------------------------
# 3. Calcul des fractions
# -----------------------------------------------
# On divise l'aire de l'intersection par l'aire totale du polygone mtp_l
inter$Frac <- inter$area_inter / (inter$Area * 1e6)

# -----------------------------------------------
# 4. Agrégation
# -----------------------------------------------
# Si plusieurs intersections existent pour le même polygone et bioregion,
# on somme les fractions pour obtenir la fraction totale par ID et K_5
# On ne conserve que l'identifiant du polygone, la bioregion et la fraction
inter <- inter[, c("ID", "K_5", "Frac"), drop = TRUE]
inter <- aggregate(inter$Frac, list(inter$ID, inter$K_5), sum)
colnames(inter) <- c("ID", "K_5", "Frac")

# -----------------------------------------------
# 5. Conversion en matrice
# -----------------------------------------------
inter <- net_to_mat(inter, weight = TRUE)

---

3.6 Application R Shiny

Pour finir, je vais vous montrer un nouvel exemple d’application Shiny sur le même format que celle du TP2.

Les fichiers sources de cette application sont disponibles ici. Vous pouvez dézipper le dossier dans votre répertoire de travail.

Nous retrouvons la même structure d’application, composée de quatre fichiers principaux.

Comme dans le TP2, prenez le temps de parcourir ces fichiers et d’expérimenter en changeant le titre de l’application, en renommant l’onglet du tableau de bord, en modifiant la palette de couleurs de la cartes ou encore personnaliser les couleurs dans style.css.

3.6.1 Exercices

1. Si vous êtes en avance essayez d’ajouter la métrique de similarité Sorensen dans l’appli.

2. Si vous êtes très en avance ajouter un onglet permettant de visualiser le réseau lambda_net.