Final Results
PabloVigneron
2025-06-06
# Importation des packages
library(tidyverse)
library(gridExtra)
#library(lubridate)
library(sf)
library(mapview)
library(maptiles)
library(tidyterra)
library(kableExtra)
#library(knitr)
library(janitor)
library(RVAideMemoire)
library(emmeans)
library(car)
library(MASS)
library(lme4)
library(MuMIn)
#library(stringr)
library(stringi)
library(terra)
library(ggspatial)
library(slippymath)
library(rosm)
library(writexl)# Importation des données - Secteur J25 : Petit Tregor
base_tregor <-
readxl::read_xlsx(path = '../raw_data/export_suiviloutrelocal_telecharger_csv_2025_04_03_10h58m17.xlsx') %>%
dplyr::select(
id_dataset,
code_secteur,
code_site,
nom_site,
x_l93,
y_l93,
date_visite,
observateurs,
condition_prospection,
nom_taxon,
nom_complet_taxon,
techn_observation,
statut_observation,
nb_ep_tot,
nb_ep_w,
nb_ep_dnf,
nb_ep_df
) %>%
mutate(
annee = year(date_visite),
mois = month(date_visite),
jour = day(date_visite),
nb_ep_w = abs(nb_ep_w),
statut_observation = ifelse(str_detect(statut_observation, '^Pr'), yes = 'Présent', no = 'Absent')
) %>%
filter(nom_complet_taxon =='Lutra lutra',
!(id_dataset == 85 & code_secteur == "FR5300006"),
code_secteur != "J401",
!(annee == 2022)) %>%
dplyr::distinct(code_site, annee, .keep_all = T)
#=> Suppression de la valeur de 2022 qui est une erreur de saisie#tab_1 <- tribble(~"Comparaison", ~"n",
# "Identiques", "32")
#knitr::kable((tab_1), booktabs = TRUE,
##caption = 'A test table.')# Importation avec data 2025 - Lieue de Greve
base_greve_raw <-
readxl::read_xls(path = '../raw_data/Suivi loutre mars_25.xls',
skip = 1,
col_names = T) ## Correction des données - Mise en forme - Secteur : Lieu Greve avec data 2025
base_greve <- base_greve_raw %>%
janitor::clean_names() %>%
dplyr::select(
point_n,
itineraire,
x,
y,
commune,
cours_deau,
bilan_choix_du_site,
epreinte_oct20,
date_suivi_fev_mars21,
marquage_suivi_fev_mars21,
date_suivi_oct21,
marquage_suivi_oct21,
observateurs_suivi_oct21,
date_suivi_mars22,
marquage_suivi_mars22,
observateurs_suivi_mars22,
date_suivi_decembre22,
marquage_suivi_decembre22,
observateurs_suivi_decembre22,
date_suivi_mars23,
marquage_suivi_mars23,
observateurs_suivi_mars23,
date_suivi_octobre_24,
marquage_suivi_octobre_24,
observateurs_suivi_octobre_24,
date_suivi_mars_25,
marquage_suivi_mars_25,
observateurs_suivi_mars_25
) %>%
dplyr::rename(
code_site = point_n,
y = x,
x = y,
marquage_suivi_oct_20 = epreinte_oct20
) %>%
dplyr::mutate(date_suivi_octobre_24 = as.Date(as.numeric(date_suivi_octobre_24), origin = "1899-12-30")) %>%
dplyr::mutate(date_suivi_oct_20 = as.Date("2020-10-01")) %>%
dplyr::mutate(across(
c(
date_suivi_fev_mars21,
date_suivi_oct21,
date_suivi_mars22,
date_suivi_decembre22,
date_suivi_mars23,
date_suivi_oct_20,
date_suivi_mars_25
),
~ as.POSIXct(.x)
)) %>%
mutate(across(starts_with("marquage_suivi_"), as.character)) %>%
tidyr::pivot_longer(
cols = c(
date_suivi_oct_20,
marquage_suivi_oct_20,
date_suivi_fev_mars21,
marquage_suivi_fev_mars21,
date_suivi_oct21,
marquage_suivi_oct21,
observateurs_suivi_oct21,
date_suivi_mars22,
marquage_suivi_mars22,
observateurs_suivi_mars22,
date_suivi_decembre22,
marquage_suivi_decembre22,
observateurs_suivi_decembre22,
date_suivi_mars23,
marquage_suivi_mars23,
observateurs_suivi_mars23,
date_suivi_octobre_24,
marquage_suivi_octobre_24,
observateurs_suivi_octobre_24,
date_suivi_mars_25,
marquage_suivi_mars_25,
observateurs_suivi_mars_25
),
names_to = c(".value", "session"),
names_pattern = "(date_suivi|marquage_suivi|observateurs_suivi)_(.*)"
) %>%
rename(
date_visite = date_suivi,
statut_presence = marquage_suivi,
observateurs = observateurs_suivi
) %>%
dplyr::mutate(
annee = year(date_visite),
mois = month(date_visite),
jour = day(date_visite),
code_site = as.factor(code_site),
statut_presence = case_when(
tolower(statut_presence) == "oui" ~ 1,
suppressWarnings(as.numeric(statut_presence)) >= 1 ~ 1,
tolower(statut_presence) %in% c("non", "peut_etre") ~ 0,
TRUE ~ 0
),
statut_observation = ifelse(statut_presence == 1, "Présent", "Absent"),
statut_presence = as.factor (statut_presence),
jour = ifelse(jour == 1, NA, jour),
annee = ifelse(annee == 2002, 2021, annee),
y = str_replace_all(y, ",", ""),
x = str_replace_all(x, " ", ""),
passage = ifelse(mois <= 6, 1, 2),
date_visite = ifelse(
date_visite == as.Date("2002-02-24"),
"2021-02-24",
as.character(date_visite)
),
observateurs_clean = observateurs %>%
str_to_lower() %>%
str_replace_all(" et ", "_") %>%
str_replace_all("[[:space:]]+", "_") %>%
str_replace_all("[^a-z0-9_]", "") %>%
stri_trans_general("latin-ascii")
) %>%
dplyr::mutate(observateurs_clean = as.factor(
ifelse(
observateurs_clean == "maxime_chapelle_romain_salaun",
"maxime_chapelle_romain_salun",
observateurs_clean
)
)) %>%
dplyr::filter(!is.na(code_site), !is.na(itineraire))# Export jeu de donnée pour carte SIG
readr::write_csv(base_greve,
path = "../processed_data/base_Lieue_de_greve.csv"
)# Export jeu de donnée pour carte SIG
readr::write_csv(base_j25,
path = "../processed_data/base_Petit_tregor.csv"
)Sites d’études
# Représentation cartographique du secteur d'étude et des sites correspondants
sites_geo_1 <- base_greve %>%
mutate(
code_secteur = "J22",
y = str_replace_all(y, ",", "."),
x = str_replace_all(x, ",", "."),
y = as.numeric(y),
x = as.numeric(x)
) %>%
filter(!is.na(x) & !is.na(y)) %>%
rename(nom_site = commune) %>%
dplyr::select(nom_site, code_site, code_secteur, statut_observation, annee, x, y) %>%
st_as_sf(coords = c("x", "y"), crs = 4326) %>%
st_transform(2154)
sites_geo_2 <- base_tregor %>%
filter(code_secteur == "J25") %>%
dplyr::select(nom_site:y_l93, code_site, code_secteur, statut_observation, annee) %>%
st_as_sf(coords = c("x_l93", "y_l93"), crs = 2154)
common_cols <- dplyr::intersect(names(sites_geo_1), names(sites_geo_2))
sites_geo <- bind_rows(
sites_geo_1[, common_cols],
sites_geo_2[, common_cols]
)
mapview(sites_geo, zcol = "code_secteur", col.regions = c("darkseagreen3", "indianred1"))base_j25 <- base_j25 %>%
mutate(annee = year(date_visite))
all_combos <- expand_grid(code_site = unique(base_j25$code_site),
annee = unique(base_j25$annee))
obs_complete <- all_combos %>%
left_join(base_j25 %>% dplyr::select(code_site, annee, statut_observation),
by = c("code_site", "annee")) %>%
mutate(
statut_final = case_when(
statut_observation == "Présent" ~ "Présent",
statut_observation == "Absent" ~ "Absent",
TRUE ~ "Non prospecté"
)
)
ggplot(obs_complete, aes(x = annee, y = code_site, fill = statut_final)) +
geom_tile(color = "white", linewidth = 0.3) +
scale_fill_manual(
values = c(
"Présent" = "darkseagreen3",
"Absent" = "indianred1",
"Non prospecté" = "grey80"
)
) +
labs(x = "Année",
y = "Site",
title = "Suivi des sites de prospection - Secteur Petit Tregor (J25)",
fill = "Statut d'observation") +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.border = element_rect(color = "black", fill = NA),
legend.background = element_rect(fill = "grey95", color = "black")
)
Figure 1: Statut d’observation des sites prospectée en fonction des années
Figure 2: Statut d’observation des sites prospectée en fonction des années, 1er passage (printemps)
base_greve_winter <- base_greve %>%
dplyr::filter(passage == 2)
annees_completes <- 2020:2024
all_combos <- expand_grid(
code_site = unique(base_greve_winter$code_site),
annee = annees_completes
)
obs_complete <- all_combos %>%
left_join(
base_greve_winter %>%
dplyr::select(code_site, annee, statut_observation),
by = c("code_site", "annee")
) %>%
dplyr::mutate(
statut_final = case_when(
statut_observation == "Présent" ~ "Présent",
statut_observation == "Absent" ~ "Absent",
TRUE ~ "Non prospecté"
)
)
ggplot(obs_complete, aes(x = annee, y = code_site, fill = statut_final)) +
geom_tile(color = "white", linewidth = 0.3) +
scale_fill_manual(values = c(
"Présent" = "darkseagreen3",
"Absent" = "indianred1",
"Non prospecté" = "grey80"
)) +
labs(x = "Année",
y = "Site",
title = "Suivi des sites de prospection - Secteur Lieu greve, 2e passage (hiver)",
fill = "Statut d'observation") +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.border = element_rect(color = "black", fill = NA),
legend.background = element_rect(fill = "grey95", color = "black")
)
Figure 3: Statut d’observation des sites prospectée en fonction des années, 2e passage (hiver)
1 RESULTATS
1.1 Modelisation
1.1.1 Secteur Petit Trégor
annees_completes <- 2011:2023
combi <- expand_grid(annee = annees_completes)
resume <- base_j25 %>%
group_by(annee) %>%
summarise(
n_obs = n(),
n_presences = sum(statut_presence == "1"),
pc_presences = n_presences / n_obs,
.groups = "drop"
)
resume_complet <- combi %>%
left_join(resume, by = "annee") %>%
replace_na(list(
n_obs = 0,
n_presences = 0,
pc_presences = 0
)) %>%
dplyr::mutate(annee = factor(annee, levels = annees_completes))
ggplot(resume_complet,
aes(x = factor(annee), y = pc_presences)) +
geom_bar(stat = "identity", position = position_dodge()) +
geom_text(
aes(label = scales::percent(pc_presences, accuracy = 1)),
position = position_dodge(width = 0.9),
vjust = -0.3,
size = 3
)+
labs(x = "Année",
y = "Pourcentage de présence de la loutre",
fill = "Passage",
title = "Présence de la loutre par année - Secteur Petit Trégor") +
scale_y_continuous(labels = scales::percent_format(accuracy = 1),
limits = c(0, 1)) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.border = element_rect(
color = "black",
fill = NA,
linewidth = 0.5
),
axis.text.x = element_text(angle = 45, hjust = 1),
panel.background = element_rect(fill = "grey94", color = "white"),
plot.background = element_rect(fill = "white", color = NA),
legend.background = element_rect(fill = "grey94", color = "black"),
panel.grid.major = element_line(color = "white"),
panel.grid.minor = element_line(color = "white")
)
Figure 4: Pourcentages de présence de la loutre en fonction des passages et des années
| Modèle | Fonction de lien | Estimate coefficient, annee | p values |
|---|---|---|---|
| statut_presence ~ annee + (1 | code_site) | Logit | 0.02853 | 0.411 NS |
| statut_presence ~ annee + (1 | code_site) | Probit | 0.013919 | 1.39e-09 *** |
| statut_presence ~ annee + (1 | code_site) | cauchit | 8.067e-02 | <2e-16 *** |
base_j25_year_numeric <- base_j25 %>%
mutate(annee = as.numeric(annee),
statut_presence = as.factor(statut_presence))mod3.0 <- glmer(statut_presence~annee+(1|code_site), family=binomial(link = "logit"), data = base_j25_year_numeric)#plotresid(mod3.0)#summary(mod3.0)mod3.1 <- glmer(statut_presence~annee+(1|code_site), family=binomial(link = "probit"), data = base_j25_year_numeric)#plotresid(mod3.1)#summary(mod3.1)mod3.2 <- glmer(statut_presence~annee+(1|code_site), family=binomial(link = "cauchit"), data = base_j25_year_numeric)#plotresid(mod3.2)#summary(mod3.2)1.2 Modelisation
1.2.1 Secteur de la Lieue de grève
base_greve_passage <- base_greve %>%
mutate(
passage_label = ifelse(mois <= 6, "Passage printemps", "Passage hiver"),
passage_label = factor(passage_label, levels = c("Passage printemps", "Passage hiver"))
)
annees <- sort(unique(base_greve_passage$annee))
passages <- c("Passage printemps", "Passage hiver")
combi <- expand_grid(annee = annees, passage_label = passages)
resume <- base_greve_passage %>%
group_by(annee, passage_label) %>%
summarise(
n_obs = n(),
n_presences = sum(statut_presence == '1'),
pc_presences = n_presences / n_obs,
.groups = "drop"
)
resume_complet <- combi %>%
left_join(resume, by = c("annee", "passage_label")) %>%
replace_na(list(
n_obs = 0,
n_presences = 0,
pc_presences = 0
))
ggplot(resume_complet, aes(x = factor(annee), y = pc_presences, fill = passage_label)) +
geom_bar(stat = "identity", position = position_dodge()) +
geom_text(aes(label = scales::percent(pc_presences, accuracy = 1)),
position = position_dodge(width = 0.9), vjust = -0.3, size = 3) +
scale_fill_manual(
values = c("Passage printemps" = "gray50", "Passage hiver" = "grey20")
) +
labs(
x = "Année",
y = "Pourcentage de présence de la loutre",
fill = "Passage",
title = "Présence de la loutre par année et par passage - Secteur Lieue de Grève"
) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1)) +
theme_minimal()+
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.border = element_rect(
color = "black",
fill = NA,
linewidth = 0.5
),
axis.text.x = element_text(angle = 45, hjust = 1),
panel.background = element_rect(fill = "grey94", color = "white"),
plot.background = element_rect(fill = "white", color = NA),
legend.background = element_rect(fill = "grey94",color = "black"),
panel.grid.major = element_line(color = "white"),
panel.grid.minor = element_line(color = "white")
)
Figure 5: Pourcentages de présence de la loutre en fonction des passages et des années
| Modèle | Fonction de lien | Estimate coefficient, annee | p values |
|---|---|---|---|
| statut_presence ~ annee + passage + (1 | code_site) | Logit | 0.001335 | 0.3071 NS |
| statut_presence ~ annee + passage + (1 | code_site) | Probit | 0.002219 | 0.5060 NS |
| statut_presence ~ annee + passage + (1 | code_site) | cauchit | -0.02917 | 0.7001 NS |
greve_year_factor_season <- base_greve%>%
mutate(annee = as.numeric(annee),
passage = as.factor(passage))mod5.0<- glmer(statut_presence~annee + passage + (1|code_site), family=binomial(link = "logit"), data = greve_year_factor_season)#plotresid(mod5.0)#summary(mod5.0)mod5.1<- glmer(statut_presence~annee + passage + (1|code_site), family=binomial(link = "probit"), data = greve_year_factor_season)#plotresid(mod5.1)#summary(mod5.1)mod5.2<- glmer(statut_presence~annee + passage + (1|code_site), family=binomial(link = "cauchit"), data = greve_year_factor_season)#plotresid(mod5.2)#summary(mod5.2)mod6.0<- glmer(statut_presence~annee + (1|code_site), family=binomial(link = "logit"), data = greve_year_factor_season)#plotresid(mod6.0)#summary(mod6.0)#mod6.1<- glmer(statut_presence~annee + (1|code_site), family=binomial(link = "probit"), data = greve_year_factor_season)#plotresid(mod6.1)#summary(mod6.1)mod6.2<- glmer(statut_presence~annee + (1|code_site), family=binomial(link = "cauchit"), data = greve_year_factor_season)#plotresid(mod6.2)#summary(mod6.2)| Modèle | Fonction de lien | Estimate coefficient, annee | p values |
|---|---|---|---|
| statut_presence ~ annee + (1 | code_site) | Logit | 0.042909 | <2e-16 *** |
| statut_presence ~ annee + (1 | code_site) | cauchit | 0.01462 | 0.836 NS |
1.3 Simulation (Secteur Petit Trégor)
library(tidyverse)
library(lme4)
load(file = "../processed_data/base_j25_year_numeric.rda")
source(file = "../R/simulation.R")set.seed(123)df <- base_j25_year_numeric %>%
mutate(code_site = as.factor(code_site),
statut_presence = as.integer(statut_presence) - 1,
date_annee = lubridate::yday(date_visite))df_rand <- df %>%
group_by(code_site) %>%
mutate(annee_perm = gtools::permute(annee)) %>%
ungroup() %>%
mutate(annee_perm_index = annee_perm - min(annee_perm))mod <- glmer(statut_presence~annee_perm+(1|code_site), family=binomial(link = "logit"), data = df_rand)
#summary(mod)
pente_obs <- summary(mod)$coefficients[2,1]df_trend1 <- df_rand %>%
ajouter_presences_annuelles(n_pres_suppl_par_an = 1) #%>%
# cbind(statut_presence_rand = df_rand$statut_presence) %>%
# mutate(presence_ajoutee = (statut_presence_rand != statut_presence_sim))2 Test : le modèle permet-il de détecter la tendance ?
On fait tourner le modèle avec statut_presence_sim en variable dépendante et la variable annee_perm en prédicteur, avec l’effet aléatoire du code_site.
2.1 Test initial
mod_tend1 <- glmer(statut_presence_sim ~
annee_perm +
(1|code_site),
family=binomial(link = "logit"),
data = df_trend1)
summary(mod_tend1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: statut_presence_sim ~ annee_perm + (1 | code_site)
## Data: df_trend1
##
## AIC BIC logLik deviance df.resid
## 513.7 525.8 -253.9 507.7 404
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4232 -0.6662 -0.5253 1.0032 2.2171
##
## Random effects:
## Groups Name Variance Std.Dev.
## code_site (Intercept) 0.5593 0.7478
## Number of obs: 407, groups: code_site, 60
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.645e+02 1.956e+00 -84.13 <2e-16 ***
## annee_perm 8.121e-02 9.736e-04 83.41 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## annee_perm -0.997
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.121856 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?Avec p<1E-16 et une pente positive on a bien détecté la tendance à l’augmentation du signal au fil des années.
2.2 Test avec une tendance plus faible
set.seed(124)
df_trend1 <- df_rand %>%
ajouter_presences_annuelles(n_pres_suppl_par_an = 0.1)
mod_trend1 <- glmer(statut_presence_sim ~
annee_perm +
(1|code_site),
family=binomial(link = "logit"),
data = df_trend1)
summary(mod_trend1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: statut_presence_sim ~ annee_perm + (1 | code_site)
## Data: df_trend1
##
## AIC BIC logLik deviance df.resid
## 440.9 453.0 -217.5 434.9 404
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1949 -0.5090 -0.3726 -0.3383 2.3383
##
## Random effects:
## Groups Name Variance Std.Dev.
## code_site (Intercept) 0.9314 0.9651
## Number of obs: 407, groups: code_site, 60
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -19.599419 69.731050 -0.281 0.779
## annee_perm 0.009035 0.034574 0.261 0.794
##
## Correlation of Fixed Effects:
## (Intr)
## annee_perm -1.000
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 1 negative eigenvaluesTendance non détectée. Il faut dire qu’elle est très faible :
table(df_trend1$annee_perm[df_trend1$statut_presence == 1])##
## 2011 2012 2014 2016 2017 2018 2019 2023
## 11 8 10 11 16 14 14 14table(df_trend1$annee_perm[df_trend1$statut_presence_sim == 1])##
## 2011 2012 2014 2016 2017 2018 2019 2023
## 11 8 10 11 16 14 14 15my_point_plot(df_trend1)
2.3 Avec un peu plus de tendance
set.seed(125)
df_trend1 <- df_rand %>%
ajouter_presences_annuelles(n_pres_suppl_par_an = 0.2)
mod_trend1 <- glmer(statut_presence_sim ~
annee_perm +
(1|code_site),
family=binomial(link = "logit"),
data = df_trend1)
summary(mod_trend1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: statut_presence_sim ~ annee_perm + (1 | code_site)
## Data: df_trend1
##
## AIC BIC logLik deviance df.resid
## 455.3 467.3 -224.6 449.3 404
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1440 -0.5519 -0.4125 0.8539 2.3055
##
## Random effects:
## Groups Name Variance Std.Dev.
## code_site (Intercept) 0.7748 0.8802
## Number of obs: 407, groups: code_site, 60
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.913e+01 1.649e-03 -23731.7 <2e-16 ***
## annee_perm 1.877e-02 8.392e-05 223.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## annee_perm -0.010
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0438417 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?On a ici détection de la tendance, sans avoir ajouté énormément de présences :
table(df_trend1$annee_perm[df_trend1$statut_presence == 1])##
## 2011 2012 2014 2016 2017 2018 2019 2023
## 11 8 10 11 16 14 14 14table(df_trend1$annee_perm[df_trend1$statut_presence_sim == 1])##
## 2011 2012 2014 2016 2017 2018 2019 2023
## 11 8 10 12 17 15 15 16my_point_plot(df_trend1)
3 Généralisation
Pour chacun des tests ci-dessous, on enchaîne 100 fois les opérations suivantes :
- Permutations des années pour gommer le signal
- Ajout de présences
- Calage du modèle
- Récupération de la pente associée à l’année et de la p-value associée
- Histogramme des pentes
3.1 En l’absence de tendance
On s’attend pour ce modèle “neutre” à ce que les pentes soient non significatives dans 95% des cas et que la distribution soit centrée sur zéro et symétrique.
set.seed(NULL)
test2 <- tester_tendance(df = df,
n_pres_suppl_par_an = 0,
n_permutations = 100)
my_histo(test2) +
geom_vline(xintercept = pente_obs, col = "black", lwd = 1)
Figure 6: Histogramme des pentes (n=100) obtenues par le modèle après permutation des années. Les couleurs permettent de distinguer les pentes selon leur signe et leur significativité au seuil de 5%.
pc_detec_aug_0 <- nrow(test2[test2$tendance == "Augmentation",]) / nrow(test2)On voit que la pente estimée sur le jeu de données non permuté (ligne verticale noire), de valeur très faible et non significative, est bien dans la gamme des valeurs NS après permutation.
On va ensuite recommencer avec des tendances de plus en plus nettes, c’est-à-dire en rajoutant de plus en plus de présences chaque année.
3.2 Tendance faible
test2 <- tester_tendance(df = df,
n_pres_suppl_par_an = 0.1,
n_permutations = 100)
my_histo(test2)
pc_detec_aug_0.1 <- nrow(test2[test2$tendance == "Augmentation",]) / nrow(test2)Pas vraiment de différence avec l’absence de tendance. Pas très surprenant car ça revient à ne rajouter qu’une présence en dernière année.
3.3 Tendance modérée
test2 <- tester_tendance(df = df,
n_pres_suppl_par_an = 0.2,
n_permutations = 100)
my_histo(test2)
pc_detec_aug_0.2 <- nrow(test2[test2$tendance == "Augmentation",]) / nrow(test2)Là ça change. La plupart des pentes sont significatices. Bizarrement on en a plus de négatives qu’avec les jeu de données “neutres”.
3.4 Tendance forte
test2 <- tester_tendance(df = df,
n_pres_suppl_par_an = 0.5,
n_permutations = 100)
my_histo(test2)
pc_detec_aug_0.5 <- nrow(test2[test2$tendance == "Augmentation",]) / nrow(test2)3.5 Tendance très forte
test2 <- tester_tendance(df = df,
n_pres_suppl_par_an = 1,
n_permutations = 100)
my_histo(test2)
pc_detec_aug_1 <- nrow(test2[test2$tendance == "Augmentation",]) / nrow(test2)4 Synthèse
On calcule un taux de détection pour différentes valeurs de la tendance injectée. Pour limiter le temps de calcul on ne recommence pas pour les valeurs déjà calculées.
tendances <- c(0.3, 0.4, 0.6, 0.7, 0.8, 0.9, 1.5, 2)
resultat <- map(
.x = tendances,
.f = get_pc_sig,
df = df,
n_permutations = 100
)
pc_detec_aug <- c(
pc_detec_aug_0,
pc_detec_aug_0.1,
pc_detec_aug_0.2,
resultat[1:2],
pc_detec_aug_0.5,
resultat[3:6],
pc_detec_aug_1,
resultat[7:8]
) %>%
unlist()
tendance <- c(seq(0, 1, 0.1), 1.5, 2)
df_resume <- data.frame(tendance,
pc_detec_aug)
ggplot(data = df_resume,
aes(x = tendance,
y = pc_detec_aug)) +
geom_bar(stat = "identity",
fill = "darkgreen") +
geom_smooth(col = "red", se = FALSE) +
scale_y_continuous(limits = c(0, 1), labels = scales::percent) +
labs(x = "Tendance (nb de présences annuelles ajoutées)",
y = "Pourcentage de détection de la tendance à l'augmentation")