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.
Le jeu de données fourni est une table avec une ligne par commune bretonne et les colonnes suivantes :
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.
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.
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")
Avant de commencer la modélisation de nb_esp
, il est essentiel d’examiner graphiquement les variables disponibles :
nb_esp
.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.
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.
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)
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.
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.
plot(modele)
Sur les graphiques de diagnostic ci-dessus, on voit que pour certaines communes :
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”.
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)
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.
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.
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.
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.
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.
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,
)
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.
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.
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.
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.