Test d’évaluation de la richesse floristique communale potentielle en Bretagne

1 Objectif

Le présent document vise à explorer si les données présentes dans la base du CBN Brest, couplées avec la “Carte des Grands Types de Végétation” (CGTV), permettent d’estimer la richesse potentielle en espèces de chaque commune.

La principale difficulté réside dans la grande variabilité intercommunale d’effort de prospection.

Par richesse potentielle, on entend le nombre d’espèces qui seraient répertoriées présentes sur la commune si elle était très bien prospectée.

2 Approche analytique

2.1 Le jeu de données

Le jeu de données fourni est une table avec une ligne par commune bretonne et les colonnes suivantes :

• cd_insee : code insee de la commune
• nb_obs : nombre d’observation floristique dans la commune
• nb_esp : nombre d’espèces observées actuellement dans la commune
• nb_station : non utilisé
• Shape area : superficie de la commune
• Shannon : indice de shannon calculé à partir de la CGTV : diversité relative des milieux naturels et semi-naturels
• Hémérobie : taux d’hémérobie calculé à partir de la CGTV : Impact de l’influence humaine sur les végétations
• nb_gtm : nombre de grand types de milieu dans la commune, y compris les milieux non naturels ou semi naturels. Des milieux non naturels peuvent amener de nouveaux cortèges d’espèces (murs…)

Il a été produit par le CBN en croisant ses données de relevés floristiques et les informations géographiques sur le découpage communal et la CGTV.

2.2 Analyses

L’approche comprend deux étapes :

• Modélisation le nombre d’espèces recensées dans chaque commune en fonction de l’effort de prospection (évalué à travers le nombre d’observations floristiques) et de variables connues pour influencer la richesse. Cette première phase est itérative pour obtenir un modèle qui soit à la fois satisfaisant (examen du r², des coefficients et des graphiques de diagnostic) et le plus simple possible.

• Application du modèle afin de prédire, pour chaque commune, le nombre d’espèces qui y serait observé si elle était bien prospectée. Pour le présent test, nous avons choisi un effort de prospection de n_obs = 20000 observations.

3 Chargement des données et packages

On charge les librairies R nécessaires, ainsi que les données. Quelques variables sont renommées pour respecter les règles de nommage et être plus explicites.

library(tidyverse)
library(COGiter)
library(PerformanceAnalytics)
library(mapview)
library(readxl)
library(downloadthis)

data <- readxl::read_xlsx("../raw_data/communes_cgtv_variables.xlsx") %>% 
  select(comm_id = fid,
         comm_insee = cd_insee,
         nb_esp,
         comm_surf = shape_area,
         nb_obs,
         gtm_shannon = shannon,
         gtm_heme = hemerobie,
         gtm_nat = richesse,
         gtm_tot = nb_gtm
         )

toutes_communes_var_quant <- data %>% 
  select(-comm_id) %>% 
  column_to_rownames("comm_insee")

4 Examen des variables

Avant de commencer la modélisation de nb_esp, il est essentiel d’examiner graphiquement les variables disponibles :

• Distributions (histogrammes) car on aime bien les cloches.
• Nuages de points par paires de variables afin d’identifier de potentielles redondances ou des relations non linéaires entre les variables explicatives et nb_esp.

4.1 Variables brutes

La variable que l’on cherche à prédire, nb_esp, varie de 2 à 663 espèces, pour un effort de prospection entre 2 à 17 179 relevés par commune.

chart.Correlation(toutes_communes_var_quant,
                  histogram = TRUE,
                  pch = 19)

Figure 1: Corrélation entre le nombre d’espèces recensées dans les communes et les variables explicatives. Les histogrammes sont sur la diagonale. Les nuages bivariés sont en-dessous, les coefficients de corrélation de Pearson et la significativité de la relation au-dessus.

Il apparaît que les distributions des surfaces communales comm_surf et du nombre d’observations nb_obs sont très asymétriques donc on leur applique une transformation logarithmique.

Les variables gtm_nat et gtm_tot sont très fortement liées \(\Rightarrow\) redondance. On ne conserve que gtm_tot.

En bivarié, le nombre d’espèces observées nb_esp est avant tout lié à l’effort de prospection nb_obs. Viennent ensuite les variables dérivées des GTM puis la surface.

4.2 Après tri et transformations

toutes_communes_var_quant <- toutes_communes_var_quant %>% 
  mutate(comm_surf = log10(comm_surf),
         nb_obs = log10(nb_obs)) %>% 
  select(-gtm_nat)

sel_communes_var_quant <- toutes_communes_var_quant %>% 
  filter(nb_obs > log10(params$nb_obs_mini))

communes_ecartees_pour_insuffisance_de_donnees <- setdiff(rownames(toutes_communes_var_quant), rownames(sel_communes_var_quant))

Nombre de communes écartées car ne disposant pas de 200 observations : 134

chart.Correlation(sel_communes_var_quant,
                  histogram = TRUE,
                  pch = 19)

Figure 2: Corrélation entre le nombre d’espèces recensées dans les communes et les variables explicatives triées et transformées, après mise à l’écart des communes insuffisamment prospectées. Les histogrammes sont sur la diagonale. Les nuages bivariés sont en-dessous, les coefficients de corrélation de Pearson et la significativité de la relation au-dessus.

5 Analyse de la richesse

5.1 Analyse préliminaire

Comme la relation bivariée entre nb_esp et nb_obs est curvi-linéaire, on introduit un terme quadratique nb_obs².

Le modèle prend donc la forme :

\(nb_esp = a.log_{10}(nb_{obs}) + b.(log_{10}(nb_{obs}))² + c.log_{10}(Surf_{comm}) + d.Shannon_{gtm} + e.Heme_{gtm}) + f.Tot_{gtm} + g\)

modele <- lm(nb_esp ~ nb_obs + I(nb_obs^2) + comm_surf + gtm_shannon + gtm_heme + gtm_tot,
             data = sel_communes_var_quant)

5.1.1 Résumé du modèle

summary(modele)

## 
## Call:
## lm(formula = nb_esp ~ nb_obs + I(nb_obs^2) + comm_surf + gtm_shannon + 
##     gtm_heme + gtm_tot, data = sel_communes_var_quant)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -174.201  -21.035   -0.277   18.786  139.502 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -423.2483    78.6126  -5.384 8.96e-08 ***
## nb_obs       117.2268    42.8655   2.735  0.00635 ** 
## I(nb_obs^2)   15.0346     7.0894   2.121  0.03418 *  
## comm_surf     -2.6452     3.8662  -0.684  0.49401    
## gtm_shannon   68.4393    26.5435   2.578  0.01006 *  
## gtm_heme      29.0728     6.0581   4.799 1.82e-06 ***
## gtm_tot        6.5004     0.5991  10.850  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35 on 1063 degrees of freedom
## Multiple R-squared:  0.8428, Adjusted R-squared:  0.8419 
## F-statistic: 949.8 on 6 and 1063 DF,  p-value: < 2.2e-16

Interprétation rapide

Le coefficient de détermination ajusté du modèle est 0.842, ce qui signifie que 84.2% de la variabilité de nb_esp est statistiquement expliqué par les variables explicatives.

Au seuil de 5%, tous les coefficients sont significativement non nuls à l’exception de la surface communale, donc le choix initial des variables semble pertinent.

5.1.2 Graphiques de diagnostic

On examine une série standard de nuages de point mettant en jeu les résidus (valeur prédite - valeur observée = erreur du modèle), les valeurs prédites et la distance de Cook (influence de l’observation sur la construction du modèle). La validation du modèle requiert les conditions ci-dessous.

• Residuals vs fitted et Scale-location (Un point par commune) \(\Rightarrow\) Absence de lien.
• Diagramme quantile-quantile des résidus (Un point par pourcentile) \(\Rightarrow\) Points sur la ligne, ce qui signifie que la distribution des résidus est gaussienne.
• Residuals vs Leverage (Un point par commune) \(\Rightarrow\) Absence de lien ; ce graphique peut aussi permettre l’identification de communes exagérément influentes.

plot(modele)

Sur les graphiques de diagnostic ci-dessus, on voit que pour certaines communes :

• les résidus sont très importants (= forte erreur du modèle, cas de communes où la richesse est sur-ou sous-estimée d’une centaine d’espèces)
• le “bras de levier” (distance de Cook) est très fort, ce qui indique une influence majeure (exagérée) sur le modèle

5.2 Suppression des communes qui ne rentrent pas dans le schéma général

Ici on cherche une règle générale qui s’appliquerait à l’ensemble des communes de la région. Il n’est donc pas dérangeant d’écarter des observations “atypiques”.

5.2.1 Identification de ces communes

On utilise une régle “à la hache” consistant à écater les observations dont la distance de Cook excède \(4/n\), \(n\) étant le nombre d’observations (ref).

dist_cook <- cooks.distance(modele) %>% 
  as.data.frame() %>% 
  rownames_to_column() %>% 
  set_names(c("comm_insee", "cooks_distance")) %>% 
  mutate(quatre_sur_n = 4 / n(),
         garder = cooks_distance < quatre_sur_n)

communes_a_conserver <- dist_cook %>% 
  filter(garder) %>% 
  pull(comm_insee)

communes_ecartees_pour_cook <- setdiff(dist_cook$comm_insee,
                                       communes_a_conserver)

5.2.2 Carte des communes écartées

communes_bzh <- COGiter::communes_geo %>% 
  mutate(dept = str_sub(DEPCOM, 1, 2)) %>% 
  filter(dept %in% c("22", "29", "35", "56"))

carte <- mapview(communes_bzh %>% filter(!(DEPCOM %in% communes_a_conserver)))

carte@map

Figure 3: Carte des communes écartées du jeu de données pour construire le modèle sur la base de la distance de Cook.

Il y a manifestement un effet géographique sur cette sélection : les communes écartées sont concentrées sur certains secteurs. Il manque peut-être une ou des variables explicatives.

5.2.3 Modèle sans ces communes

sel_communes_var_quant_propre <- sel_communes_var_quant %>% 
  filter(rownames(sel_communes_var_quant) %in% communes_a_conserver)

modele <- lm(nb_esp ~ nb_obs + I(nb_obs^2) + comm_surf + gtm_shannon + gtm_heme + gtm_tot,
             data = sel_communes_var_quant_propre)

summary(modele)

## 
## Call:
## lm(formula = nb_esp ~ nb_obs + I(nb_obs^2) + comm_surf + gtm_shannon + 
##     gtm_heme + gtm_tot, data = sel_communes_var_quant_propre)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -83.28 -18.31   0.74  18.11  81.48 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -353.0491    72.4572  -4.873 1.28e-06 ***
## nb_obs        78.3150    41.2296   1.899  0.05779 .  
## I(nb_obs^2)   21.9205     6.9490   3.154  0.00166 ** 
## comm_surf     -1.1874     3.2055  -0.370  0.71116    
## gtm_shannon   60.2112    22.0757   2.727  0.00650 ** 
## gtm_heme      26.7949     5.1520   5.201 2.41e-07 ***
## gtm_tot        5.5598     0.5043  11.024  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 27.14 on 986 degrees of freedom
## Multiple R-squared:  0.8794, Adjusted R-squared:  0.8787 
## F-statistic:  1198 on 6 and 986 DF,  p-value: < 2.2e-16

plot(modele)

Les graphiques sont nettement plus satisfaisants.

6 Modèle retenu

En dernière étape, on essaye de simplifier le modèle en éliminant les variables qui ne contribuent pas à l’améliorer, par une méthode “pas-à-pas” basée sur le critère d’Akaïke.

6.1 Elagage

summary(modele_final)

## 
## Call:
## lm(formula = nb_esp ~ nb_obs + I(nb_obs^2) + gtm_shannon + gtm_heme + 
##     gtm_tot, data = sel_communes_var_quant_propre)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -83.301 -18.235   0.744  18.118  81.752 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -364.378     65.657  -5.550 3.68e-08 ***
## nb_obs        78.773     41.193   1.912  0.05613 .  
## I(nb_obs^2)   21.807      6.939   3.143  0.00172 ** 
## gtm_shannon   64.013     19.537   3.277  0.00109 ** 
## gtm_heme      27.556      4.722   5.835 7.28e-09 ***
## gtm_tot        5.465      0.434  12.591  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 27.13 on 987 degrees of freedom
## Multiple R-squared:  0.8794, Adjusted R-squared:  0.8788 
## F-statistic:  1439 on 5 and 987 DF,  p-value: < 2.2e-16

On voit que la superficie de la commune a été écartée des variables explicatives.

plot(modele_final)

On a toujours quelques observations à forte influence mais les trois autres graphiques sont Ok.

6.2 Relation richesse observée - richesse prédite

Le relation est étroite entre la richesse observée et la richesse prédite. La droite de régression se confond avec la droite identité d’équation \(y=x\), ce qui indique un bon ajustement sur l’ensemble de la gamme de richesses.

test <- modele_final$fitted.values %>% 
  cbind(modele_final$model) %>% 
  rename(nb_esp_pred = 1)

ggplot(data = test,
       aes(x = nb_esp_pred,
           y = nb_esp)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(x = "Richesse prédite",
       y = "Richesse observée")

Figure 4: Nuage de point de la richesse observée en fonction de la richesse prédite (n = 993). Chaque point représente une commune.

7 Richesse potentielle des communes

7.1 Calcul

On utilise le modèle final pour prédire ce que serait la richesse observée dans chaque commune si elle disposait de 20000 observations.

NB Au stade de la prédition, les communes qui avaient été écartées pour caler le modèle sont réintégrées au jeu de données.

new_data <- toutes_communes_var_quant %>% 
  mutate(nb_obs = log10(params$nb_obs_simu)) # car le nb_obs dans le modèle est en log

prediction <- predict(modele, newdata = new_data) %>% 
  as.data.frame() %>% 
  rownames_to_column() %>% 
  set_names(c("comm_insee", "richesse_pot"))

prediction %>%
  downloadthis::download_this(
    output_name = "richesse_potentielle",
    output_extension = ".xlsx",
    button_label = "Télécharger en Excel",
    button_type = "success",
    has_icon = TRUE,
    icon = "fa fa-save",
    csv2 = TRUE,
  )

Télécharger en Excel

7.2 Cartographie

Le modèle prédit des diversités maximales sur le littoral, en particulier Sud et Ouest.

donnees_carte <- communes_bzh %>% 
  left_join(y = prediction,
            by = c("DEPCOM" = "comm_insee"))

carte <- mapview(donnees_carte,
                 zcol = "richesse_pot")

carte@map

Figure 5: Cartographie de la richesse potentielle communale estimée pour n=20000 observations.

7.3 Distribution

richesse_mediane <- median(prediction$richesse_pot)

La richesse communale potentielle médiane, estimée pour un effort de prospection de 20000 observations, est de 597 espèces.

ggplot(data = prediction,
       aes(x = richesse_pot)) +
  geom_histogram(fill = "darkgreen") +
  labs(x = "Richesse communale potentielle",
       y = "Nombre de communes") +
  geom_vline(xintercept = richesse_mediane,
             linetype = "dashed")

Figure 6: Distribution des richesses potentielles prédites. La ligne verticale indique la médiane.

7.4 Effet de l’effort de prospection

On peut s’interroger sur l’effet du nombre d’observations retenu (ci-dessus fixé à 20000) sur la richesse potentielle estimée.

nb_obs_seq <- c(1e2, 3e2, 1e3, 3e3, 1e4, 18000, 3e4, 5e4, 7e4, 1e5)

Pour répondre à cette question, on fait tourner le modèle, en prédiction, sur l’ensemble des communes pour des efforts de prospection variant de 100 à 100 000 observations. On calcule ensuite, pour chaque nombre d’observations, la richesse communale potentielle médiane.

richesse_pot_mediane <- data.frame()

for(n_obs_simu in nb_obs_seq)
  
{
  
new_data <- toutes_communes_var_quant %>% 
  mutate(nb_obs = log10(n_obs_simu)) # car le nb_obs dans le modèle est en log

prediction <- predict(modele, newdata = new_data) %>% 
  as.data.frame() %>% 
  rownames_to_column() %>% 
  set_names(c("comm_insee", "richesse_pot"))

richesse_pot_med <- data.frame(richesse_pot = median(prediction$richesse_pot),
                               n_obs_simu)
                               
richesse_pot_mediane <- richesse_pot_mediane %>% 
  rbind(richesse_pot_med)

}

ggplot(data = richesse_pot_mediane,
       aes(x = n_obs_simu,
           y = richesse_pot)) +
  geom_point(col = "darkgreen",
             size = 3) +
#  geom_smooth(se = FALSE) +
  geom_line(col = "darkgreen") +
  scale_y_continuous(limits = c(0, NA)) +
  scale_x_continuous(breaks = seq(0, 100000, 10000),
                     labels = function(x) format(x, big.mark = " ", scientific = FALSE)) +
  labs(x = "Nombre d'observations",
       y = "Richesse communale potentielle médiane") +
  geom_hline(yintercept = richesse_mediane,
             linetype = "dashed",
             col = "red") +
  geom_vline(xintercept = params$nb_obs_simu,
             linetype = "dashed",
             col = "red") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

Figure 7: Effet de l’effort de prospection sur la richesse communale potentielle médiane. L’intersection des lignes pointillées indique les paramètres retenus pour l’ensemble de l’étude.

8 Profil des communes écartées

Certaines des communes ont été écartées car elles n’ont pas été suffisamment prospectées (n = 134) et d’autres car elles influençaient exagérément le modèle (n = 77).

On peut comparer les variables entre les communes conservées dans les analyses et ces deux lots de communes.

comparaison_conservees_ecartees <- data %>% 
  mutate(ecartee = case_when(
    comm_insee %in% communes_ecartees_pour_cook ~ "Trop influente",
    comm_insee %in% communes_ecartees_pour_insuffisance_de_donnees ~ "Données insuffisantes",
    TRUE ~ "Conservée")
    )

comparaison_conservees_ecartees <- comparaison_conservees_ecartees %>% 
  select(-comm_id) %>% 
  pivot_longer(nb_esp:gtm_tot,
               names_to = "variable",
               values_to = "valeur")

ggplot(data = comparaison_conservees_ecartees,
       aes(x = valeur,
           fill = ecartee)) +
  geom_histogram(alpha = 0.5) +
  facet_wrap(~variable,
             scales = "free") +
  labs(x = "Valeur",
       y = "Nombre de communes",
       fill = "Commune")

Figure 8: Comparaison de la distribution des variables (histogramme) entre les communes selon qu’elles sont écartées ou conservées.

Même type de graphique mais en densités. Comme la surface sous chaque courbe est 1, les communes écartées, moins nombreuses que celles qui sont conservées, apparaissent plus clairement.

ggplot(data = comparaison_conservees_ecartees,
       aes(x = valeur,
           fill = ecartee)) +
  geom_density(alpha = 0.5) +
  facet_wrap(~variable,
             scales = "free") +
  labs(x = "Valeur",
       y = "Nombre de communes",
       fill = "Commune")

Figure 9: Comparaison de la distribution des variables (diagramme de densité) entre les communes selon qu’elles sont écartées ou conservées.

Comme la distribution de nb_obs est très concentrée sur les petites valeurs, on peut la représenter en échalle log.

ggplot(data = comparaison_conservees_ecartees %>% filter(variable == "nb_obs"),
       aes(x = valeur,
           fill = ecartee)) +
  geom_histogram(alpha = 0.3,
                 col = "black") +
  scale_x_log10() +
  labs(x = "Nombre d'observations",
       y = "Nombre de communes",
       fill = "Commune")

Figure 10: Comparaison de la distribution du nombre d’observations (histogramme) entre les communes selon qu’elles sont écartées ou conservées.

Même type de graphique mais en densités. Comme la surface sous chaque courbe est 1, les communes écartées, moins nombreuses que celles qui sont conservées, apparaissent plus clairement.

ggplot(data = comparaison_conservees_ecartees %>% filter(variable == "nb_obs"),
       aes(x = valeur,
           fill = ecartee)) +
  geom_density(alpha = 0.3,
                 col = "black") +
  scale_x_log10() +
  labs(x = "Nombre d'observations",
       y = "Nombre de communes",
       fill = "Commune")

Figure 11: Comparaison de la distribution du nombre d’observations (diagramme de densité) entre les communes selon qu’elles sont écartées ou conservées.