Script : 47_ZIP_GLMM_resultats
Ce script a pour objectif d’évaluer l’impact des différents réseaux sur les tendances des indicateurs régionaux. Il a également l’objectif de regarder l’impact des changements de protocoles sur les stations au cours des années d’études grâce aux zero inflated poisson GLMM.
library(lme4)
library(aspe)
library(tidyverse)
library(aspe)
library(sjstats)
library(gt)
library(glmmTMB)
library(RVAideMemoire)
library(lmerTest) # remplace lme4 automatiquementload(file = "../processed_data/selection_especes.rda")
load(file = "../processed_data/selection_pop_ope.rda")
load(file = "../processed_data/donnees_densite_eff_abs.rda")
load(file = "../processed_data/assemblage_tab_par_ope.rda")# source(file = "../R/lm_application_modele.R")
# source(file = "../R/lm_calcul_modele.R")
source(file = "../R/zip_glmm_calcul_modele.R")
source(file = "../R/zip_glmm_application_modele.R")
source(file = "../R/glmm_calcul_modele.R")
source(file = "../R/glmm_application_modele.R")rdata_tables <- misc_nom_dernier_fichier(
repertoire = "../../../../Liste_rouge_BFC/raw_data",
pattern = "^tables")
load(rdata_tables) # Chargement des données Les modèles linéaires généralisés à effets mixtes combinent les caractéristiques des modèles linéaires généralisés (modéliser des variables non-normalement distribuées, spécialement des données binaires et de comptage) et des modèles à effets mixtes (modéliser des données groupées).
Nous utiliserons dans cette étude le package glmmTMB pour ajuster des modèles mixtes. La fonction glmmTMB de ce package estime les paramètres d’un modèle linéaire généralisé à effets mixtes. Les formules utilisées par glmmTMB suivent la forme reponse ~ predicteurs, avec une syntaxe spécifique pour les effets aléatoires.
La distribution de Poisson peut être ici utilisée (valeurs entières supérieures ou égales à 0. glmmTMB(y ~ x1 + x2 + …, data = …, family = poisson(link=log))
liste_periodes_etude <- list(c(2007, 2025),
c(2015, 2025)) # Liste des périodes (chaque période est un vecteur avec l'année de début et l'année de fin) # preparer la table pour jointure
ope_selection <- ope_selection %>%
dplyr::select(ope_id, pop_id, sta_id)
ope_biomasse_glm <- ope_indicateur %>%
mef_ajouter_type_protocole() %>%
mef_ajouter_ope_date() %>%
mef_ajouter_surf_calc() %>%
# mef_ajouter_objectif() %>% # Ajout des objectifs de pêches
inner_join(y = ope_selection, by = c("ope_id", "annee", "pop_id")) %>%
# left_join(coord_sta_id, by = "pop_id") %>%
mutate(
valeur = as.numeric(valeur) ,
annee = as.numeric(annee),
# annee_an = annee,
ope_surface_calculee = as.numeric(ope_surface_calculee) ,
pro_libelle = as.factor(pro_libelle),
pop_id = as.factor(pop_id),
sta_id = as.factor(sta_id),
# obj_libelle = as.factor(obj_libelle),
# indicateur = "effectif_total"
) %>%
# ) %>% # creation d'une colonne pour la période (ex : 2007-2008)
# mutate(annee = {
# debut <- floor((annee - annee_debut_periode) / 2) * 2 + annee_debut_periode
# paste0(debut)
# }, annee = as.numeric(annee)) %>%
# mutate(julian = lubridate::yday(ope_date)) %>% # Ajout de la date de pêche codée de 1 à 365
select(
ope_id,
valeur,
esp_code_alternatif,
annee,
ope_surface_calculee,
pro_libelle,
pop_id,
sta_id,
indicateur,
# x,
# y,
ope_date,
# julian,
# obj_libelle
) %>%
rename(espece = esp_code_alternatif) %>%
filter(indicateur == "biomasse")Afin d’exécuter l’analyse, deux fonctions ont été créer : zip_glmm_calcul_modele (en charge de calculer les valeurs des modèles ZIPGLMM) et zip_glmm_application_modele (le but étant d’excécuter la fonction zip_glm_calcul_modele pour l’ensemble des combinaisons periodes-especes-stade-indicateur). Les analyses prennent en considération les données à l’échelle de l’opération.
glmm_resultats <- glmm_application_modele(data = ope_biomasse_glm,
liste_periodes = liste_periodes_etude)glmm_resultats_annee <- glmm_resultats %>%
filter(row_name == "annee") %>%
select(periode,
esp_code_alternatif,
p_value,
Estimate,
sig,
family) %>%
mutate(trend = case_when(
Estimate > 0 ~ "\U2197",
Estimate < 0 ~ "\U2198",
TRUE ~ "."))# Filtrer les données selon les critères spécifiés
glmm_resultats_format <- glmm_resultats_annee %>%
filter((esp_code_alternatif == "ANG" & periode == "2007-2025") |
(esp_code_alternatif %in% c(c("ABL","BAF","BLN","BRB", "BRE","BRO","CCO","CHA","CHE", "EPI", "EPT", "GAR", "GOU","HOT","LOF", "LOR", "LOT","LPM" ,"LPP", "OBR", "PER", "ROT", "SAT", "SPI", "TAN", "TOX", "TRF","VAI","VAN","VAR")) & periode == "2007-2025"))Construction d’un tableau récapitulatif des résultats du modèle GLMM :
df_glmm_resultats <- glmm_resultats_format %>%
select(esp_code_alternatif, periode, trend, sig, Estimate, family) %>%
mutate(Estimate = round(Estimate, digits = 2)) %>%
rename(
"Code Espèce" = esp_code_alternatif,
"Periode" = periode,
"Tendance" = trend,
"Coefficient" = Estimate,
"Significativité" = sig,
"Famille de distribution" = family
) %>%
arrange(`Code Espèce`) # Ordonner par espèce, indicateur et stade
# Lignes positives et significatives
i_vert <- which(
df_glmm_resultats$Tendance == "\U2197" &
df_glmm_resultats$Significativité %in% c("*", "**", "***")
)
# Lignes négatives et significatives
i_rouge <- which(
df_glmm_resultats$Tendance == "\U2198" &
df_glmm_resultats$Significativité %in% c("*", "**", "***")
)
tab_glmm_resultats <- df_glmm_resultats %>%
flextable::flextable() %>%
flextable::autofit() %>%
flextable::colformat_int(j = c(1, 3), big.mark = " ") %>%
flextable::theme_vanilla() %>%
flextable::font(fontname = "calibri") %>%
flextable::fontsize(size = 13) %>%
flextable::bg(i = i_vert, bg = "darkolivegreen3") %>% # vert clair
flextable::bg(i = i_rouge, bg = "coral3") %>% # rouge clair
flextable::align(align = "center") %>%
flextable::set_caption("Synthèse des tendances des biomasses des populations de poissons d'eau douce de BFC")
tab_glmm_resultatsCode Espèce | Periode | Tendance | Significativité | Coefficient | Famille de distribution |
|---|---|---|---|---|---|
ABL | 2007-2025 | ↗ | *** | 13.73 | Gaussian |
ANG | 2007-2025 | ↘ | *** | -116.78 | Gaussian |
BAF | 2007-2025 | ↘ | * | -40.65 | Gaussian |
BLN | 2007-2025 | ↘ | ** | -30.40 | Gaussian |
BRB | 2007-2025 | ↘ | NS | -4.10 | Gaussian |
BRE | 2007-2025 | ↗ | NS | 12.85 | Gaussian |
BRO | 2007-2025 | ↗ | NS | 12.07 | Gaussian |
CCO | 2007-2025 | ↗ | NS | 55.95 | Gaussian |
CHA | 2007-2025 | ↘ | ** | -2.93 | Gaussian |
CHE | 2007-2025 | ↘ | NS | -10.26 | Gaussian |
EPI | 2007-2025 | ↘ | NS | -0.15 | Gaussian |
EPT | 2007-2025 | ↗ | NS | 0.04 | Gaussian |
GAR | 2007-2025 | ↘ | *** | -15.96 | Gaussian |
GOU | 2007-2025 | ↗ | NS | 7.89 | Gaussian |
HOT | 2007-2025 | ↘ | *** | -97.29 | Gaussian |
LOF | 2007-2025 | ↗ | ** | 3.14 | Gaussian |
LOR | 2007-2025 | ↗ | NS | 0.59 | Gaussian |
LOT | 2007-2025 | ↘ | NS | -11.93 | Gaussian |
LPP | 2007-2025 | ↗ | NS | 0.16 | Gaussian |
OBR | 2007-2025 | ↗ | NS | 7.71 | Gaussian |
PER | 2007-2025 | ↗ | NS | 2.19 | Gaussian |
ROT | 2007-2025 | ↗ | * | 3.07 | Gaussian |
SAT | 2007-2025 | ↘ | NS | -0.83 | Gaussian |
SPI | 2007-2025 | ↗ | *** | 7.22 | Gaussian |
TAN | 2007-2025 | ↘ | * | -7.10 | Gaussian |
TOX | 2007-2025 | ↗ | NS | 3.10 | Gaussian |
TRF | 2007-2025 | ↘ | *** | -42.17 | Gaussian |
VAI | 2007-2025 | ↘ | NS | -89.41 | Gaussian |
VAN | 2007-2025 | ↘ | *** | -17.28 | Gaussian |
liste_periodes_etude <- list(c(2007, 2025),
c(2015, 2025)) # Liste des périodes (chaque période est un vecteur avec l'année de début et l'année de fin) # preparer la table pour jointure
# ope_selection <- ope_selection %>%
# dplyr::select(ope_id, pop_id, sta_id)
#
# ope_occurence_glm <- ope_effectif_absences %>%
# mef_ajouter_type_protocole() %>%
# mef_ajouter_ope_date() %>%
# mef_ajouter_surf_calc() %>%
# # mef_ajouter_objectif() %>% # Ajout des objectifs de pêches
# inner_join(y = ope_selection, by = c("ope_id", "annee", "pop_id")) %>%
#
#
#
#
#
#
#
#
# # left_join(coord_sta_id, by = "pop_id") %>%
# mutate(
# valeur = as.numeric(valeur) ,
# annee = as.numeric(annee),
# # annee_an = annee,
# ope_surface_calculee = as.numeric(ope_surface_calculee) ,
# pro_libelle = as.factor(pro_libelle),
# pop_id = as.factor(pop_id),
# sta_id = as.factor(sta_id),
# # obj_libelle = as.factor(obj_libelle),
# # indicateur = "effectif_total"
# ) %>%
# # ) %>% # creation d'une colonne pour la période (ex : 2007-2008)
# # mutate(annee = {
# # debut <- floor((annee - annee_debut_periode) / 2) * 2 + annee_debut_periode
# # paste0(debut)
# # }, annee = as.numeric(annee)) %>%
# # mutate(julian = lubridate::yday(ope_date)) %>% # Ajout de la date de pêche codée de 1 à 365
# select(
# ope_id,
# valeur,
# esp_code_alternatif,
# annee,
# ope_surface_calculee,
# pro_libelle,
# pop_id,
# sta_id,
# indicateur,
# # x,
# # y,
# ope_date,
# # julian,
# # obj_libelle
# ) %>%
# rename(espece = esp_code_alternatif) %>%
# filter(indicateur == "biomasse")Afin d’exécuter l’analyse, deux fonctions ont été créer : zip_glmm_calcul_modele (en charge de calculer les valeurs des modèles ZIPGLMM) et zip_glmm_application_modele (le but étant d’excécuter la fonction zip_glm_calcul_modele pour l’ensemble des combinaisons periodes-especes-stade-indicateur). Les analyses prennent en considération les données à l’échelle de l’opération.
# glmm_resultats <- glmm_application_modele(data = ope_biomasse_glm,
# liste_periodes = liste_periodes_etude)# glmm_resultats_annee <- glmm_resultats %>%
# filter(row_name == "annee") %>%
# select(periode,
# esp_code_alternatif,
# p_value,
# Estimate,
# sig,
# family) %>%
# mutate(trend = case_when(
# Estimate > 0 ~ "\U2197",
# Estimate < 0 ~ "\U2198",
# TRUE ~ "."))# Filtrer les données selon les critères spécifiés
# glmm_resultats_format <- glmm_resultats_annee %>%
# filter((esp_code_alternatif == "ANG" & periode == "2007-2025") |
# (esp_code_alternatif %in% c(c("ABL","BAF","BLN","BRB", "BRE","BRO","CCO","CHA","CHE", "EPI", "EPT", "GAR", "GOU","HOT","LOF", "LOR", "LOT","LPM" ,"LPP", "OBR", "PER", "ROT", "SAT", "SPI", "TAN", "TOX", "TRF","VAI","VAN","VAR")) & periode == "2007-2025"))Construction d’un tableau récapitulatif des résultats du modèle GLMM :
# df_glmm_resultats <- glmm_resultats_format %>%
# select(esp_code_alternatif, periode, trend, sig, Estimate, family) %>%
# mutate(Estimate = round(Estimate, digits = 2)) %>%
# rename(
# "Code Espèce" = esp_code_alternatif,
# "Periode" = periode,
# "Tendance" = trend,
# "Coefficient" = Estimate,
# "Significativité" = sig,
# "Famille de distribution" = family
# ) %>%
# arrange(`Code Espèce`) # Ordonner par espèce, indicateur et stade
#
# # Lignes positives et significatives
# i_vert <- which(
# df_glmm_resultats$Tendance == "\U2197" &
# df_glmm_resultats$Significativité %in% c("*", "**", "***")
# )
# # Lignes négatives et significatives
# i_rouge <- which(
# df_glmm_resultats$Tendance == "\U2198" &
# df_glmm_resultats$Significativité %in% c("*", "**", "***")
# )
#
# tab_glmm_resultats <- df_glmm_resultats %>%
# flextable::flextable() %>%
# flextable::autofit() %>%
# flextable::colformat_int(j = c(1, 3), big.mark = " ") %>%
# flextable::theme_vanilla() %>%
# flextable::font(fontname = "calibri") %>%
# flextable::fontsize(size = 13) %>%
# flextable::bg(i = i_vert, bg = "darkolivegreen3") %>% # vert clair
# flextable::bg(i = i_rouge, bg = "coral3") %>% # rouge clair
# flextable::align(align = "center") %>%
# flextable::set_caption("Synthèse des tendances des biomasses des populations de poissons d'eau douce de BFC")
#
# tab_glmm_resultatsAfin d’appliquer les modèles GLMM, il est nécessaire de construire un jeu de données contenant les valeurs des effectifs, les espèces, les stades, les pop_id, mais aussi les types de réseaux, les surfaces d’opérations et les protocoles de pêches. Nous utilisons dans un premier temps, le jeu de données relatifs aux effectifs des opérations de pêches (ope_effectifs) car nous souhaitons observer les effets des réseaux et des protocoles sur les résultats des opérations.
liste_periodes_etude <- list(c(2007, 2025),
c(2015, 2025)) # Liste des périodes (chaque période est un vecteur avec l'année de début et l'année de fin) Nous construisons un jeu de données complet :
# Pour creer la colonne periode : on renseigne le début de période souhaiter : ici 2007
annee_debut_periode <- 2007
# preparer la table pour jointure
ope_selection <- ope_selection %>%
dplyr::select(ope_id, pop_id, sta_id)
# Ajout des coordonnées pour jointure avec df pour future analyse autocorrelation spatiale
coord_sta_id <- point_prelevement %>%
filter(pop_id %in% pop_serie_tempo) %>%
rename(x = pop_coordonnees_x, y = pop_coordonnees_y) %>%
select(pop_id, x, y)
ope_effectif_glm <- esp_ope_selection %>%
mef_ajouter_type_protocole() %>%
mef_ajouter_ope_date() %>%
# mef_ajouter_objectif() %>% # Ajout des objectifs de pêches
left_join(ope_selection, by = "ope_id") %>%
rename(valeur = lop_effectif, annee = annee.x) %>%
left_join(coord_sta_id, by = "pop_id") %>%
mutate(
valeur = as.numeric(valeur) ,
annee = as.numeric(annee),
# annee_an = annee,
ope_surface_calculee = as.numeric(ope_surface_calculee) ,
pro_libelle = as.factor(pro_libelle),
pop_id = as.factor(pop_id),
sta_id = as.factor(sta_id),
# obj_libelle = as.factor(obj_libelle),
indicateur = "effectif_total") %>%
# ) %>% # creation d'une colonne pour la période (ex : 2007-2008)
# mutate(annee = {
# debut <- floor((annee - annee_debut_periode) / 2) * 2 + annee_debut_periode
# paste0(debut)
# }, annee = as.numeric(annee)) %>%
mutate(julian = lubridate::yday(ope_date)) %>% # Ajout de la date de pêche codée de 1 à 365
select(
ope_id,
valeur,
esp_code_alternatif,
annee,
ope_surface_calculee,
pro_libelle,
pop_id,
sta_id,
indicateur,
x,
y,
ope_date,
julian,
# obj_libelle
) %>%
rename(espece = esp_code_alternatif)Afin d’exécuter l’analyse, deux fonctions ont été créer : zip_glmm_calcul_modele (en charge de calculer les valeurs des modèles ZIPGLMM) et zip_glmm_application_modele (le but étant d’excécuter la fonction zip_glm_calcul_modele pour l’ensemble des combinaisons periodes-especes-stade-indicateur). Les analyses prennent en considération les données à l’échelle de l’opération.
zip_glmm_resultats <- zip_glmm_application_modele(data = ope_effectif_glm,
liste_periodes = liste_periodes_etude)zip_glmm_resultats_annee <- zip_glmm_resultats %>%
filter(row_name == "annee") %>%
select(periode,
esp_code_alternatif,
p_value,
Estimate,
sig,
family) %>%
mutate(trend = case_when(
Estimate > 0 ~ "\U2197",
Estimate < 0 ~ "\U2198",
TRUE ~ "."))# Filtrer les données selon les critères spécifiés
zip_glmm_resultats_format <- zip_glmm_resultats_annee %>%
filter((esp_code_alternatif == "ANG" & periode == "2007-2025") |
(esp_code_alternatif %in% c(c("ABL","BAF","BLN","BRB", "BRE","BRO","CCO","CHA","CHE", "EPI", "EPT", "GAR", "GOU","HOT","LOF", "LOR", "LOT","LPM" ,"LPP", "OBR", "PER", "ROT", "SAT", "SPI", "TAN", "TOX", "TRF","VAI","VAN","VAR")) & periode == "2007-2025"))Construction d’un tableau récapitulatif des résultats du modèle GLMM :
df_zip_glmm_resultats <- zip_glmm_resultats_format %>%
select(esp_code_alternatif, periode, trend, sig, Estimate, family) %>%
mutate(Estimate = round(Estimate, digits = 2)) %>%
rename(
"Code Espèce" = esp_code_alternatif,
"Periode" = periode,
"Tendance" = trend,
"Coefficient" = Estimate,
"Significativité" = sig,
"Famille de distribution" = family
) %>%
arrange(`Code Espèce`) # Ordonner par espèce, indicateur et stade
# Lignes positives et significatives
i_vert <- which(
df_zip_glmm_resultats$Tendance == "\U2197" &
df_zip_glmm_resultats$Significativité %in% c("*", "**", "***")
)
# Lignes négatives et significatives
i_rouge <- which(
df_zip_glmm_resultats$Tendance == "\U2198" &
df_zip_glmm_resultats$Significativité %in% c("*", "**", "***")
)
tab_zip_glmm_resultats <- df_zip_glmm_resultats %>%
flextable::flextable() %>%
flextable::autofit() %>%
flextable::colformat_int(j = c(1, 3), big.mark = " ") %>%
flextable::theme_vanilla() %>%
flextable::font(fontname = "calibri") %>%
flextable::fontsize(size = 13) %>%
flextable::bg(i = i_vert, bg = "darkolivegreen3") %>% # vert clair
flextable::bg(i = i_rouge, bg = "coral3") %>% # rouge clair
flextable::align(align = "center") %>%
flextable::set_caption("Synthèse des tendances des effectifs des populations de poissons d'eau douce de BFC")
tab_zip_glmm_resultatsCode Espèce | Periode | Tendance | Significativité | Coefficient | Famille de distribution |
|---|---|---|---|---|---|
ABL | 2007-2025 | ↗ | *** | 0.13 | Hurdle_Negative_Binomial |
ANG | 2007-2025 | ↘ | NS | 0.00 | Negative_Binomial |
BAF | 2007-2025 | ↗ | *** | 0.13 | Hurdle_Negative_Binomial |
BLN | 2007-2025 | ↗ | NS | 0.01 | Hurdle_Negative_Binomial |
BRB | 2007-2025 | ↗ | *** | 0.04 | Zero_Inflated_Negative_Binomial |
BRO | 2007-2025 | ↘ | *** | -0.01 | Negative_Binomial |
CHA | 2007-2025 | ↗ | *** | 0.06 | Hurdle_Negative_Binomial |
CHE | 2007-2025 | ↗ | *** | 0.08 | Hurdle_Negative_Binomial |
GAR | 2007-2025 | ↗ | *** | 0.08 | Hurdle_Negative_Binomial |
GOU | 2007-2025 | ↗ | *** | 0.11 | Hurdle_Negative_Binomial |
HOT | 2007-2025 | ↗ | *** | 0.02 | Zero_Inflated_Negative_Binomial |
LOF | 2007-2025 | ↗ | *** | 0.07 | Hurdle_Negative_Binomial |
LOR | 2007-2025 | ↗ | NS | 0.03 | Negative_Binomial |
LOT | 2007-2025 | ↗ | *** | 0.01 | Negative_Binomial |
PER | 2007-2025 | ↗ | *** | 0.08 | Hurdle_Negative_Binomial |
ROT | 2007-2025 | ↗ | *** | 0.09 | Zero_Inflated_Negative_Binomial |
SAT | 2007-2025 | ↗ | ** | 0.10 | Negative_Binomial |
SPI | 2007-2025 | ↗ | *** | 0.11 | Zero_Inflated_Negative_Binomial |
TAN | 2007-2025 | ↗ | *** | 0.03 | Zero_Inflated_Negative_Binomial |
TRF | 2007-2025 | ↗ | *** | 0.04 | Hurdle_Negative_Binomial |
VAI | 2007-2025 | ↗ | *** | 0.04 | Zero_Inflated_Negative_Binomial |
VAN | 2007-2025 | ↘ | NS | -0.01 | Zero_Inflated_Negative_Binomial |
# save(zip_glmm_resultats_format,
# reg_indicateur_lrr_tendances_lm,
# reg_indicateur_lrr_tendances_mk_st,
# tab_glmm_resultats,
# tab_reg_indicateur_lrr_tendances_lm,
# tab_reg_indicateur_lrr_tendances_mk_st,
# file = "../processed_data/mk_glmm_lm_tendances.rda")