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(corrplot)
library(PerformanceAnalytics)
library(mapview)
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)
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 ⇒ 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)
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
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")
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
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")
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))
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")
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")
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")
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")