1 Chargement des packages, des fonctions et des données

library(tidyverse)
library(wesanderson) # Palette de couleur
library(aspe) # Traitement de la base aspe
library(mgcv) # modele GAMMs
library(gratia)
library(patchwork)
library(lme4)
library(RVAideMemoire)
library(glmmTMB)
library(car) # fonction Anova()
library(ggh4x)
library(gstat)#semi_variogramme
library(mapview)
library(COGiter)
library(sf)
library(ggplotify)
library(ggpubr)
library(gstat)# For semivariogram

load(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/selection_pop_ope.rda")
load(file = "../processed_data/selection_pop_ope_carte.rda")

rdata_tables <- misc_nom_dernier_fichier(
  repertoire = "../../../../Liste_rouge_BFC/raw_data",
  pattern = "^tables")
load(rdata_tables) # Chargement des données

pal <- wes_palette("AsteroidCity1")

Script : 45_GAMM_question

2 Objectif

2.1 Contexte du stage

L’objectif principal de l’étude est d’évaluer les tendances démographiques des espèces de poissons d’eau douce à l’échelle de la région Bourgogne Franche-comté. Des séries chronologiques longues enregistrant plus d’une dizaine d’années ou 3 générations (critère A UICN) de données seront utilisées (sur un total disponible d’une vingtaine d’années). Ce travail s’inscrit dans la perspective de la création de la Liste rouge régionale des poissons d’eau douce de Bourgogne Franche-Comté et compte une trentaine d’espèces à évaluer.

L’OFB a mis en place depuis 1990 des réseaux de suivi temporel de l’état des milieux aquatiques à partir des peuplements de poissons. Les échantillonnages sont réalisés tous les ans ou tout les deux ans en fonction des réseaux pluriannuels des stations identiques selon deux protocoles distincts : - pêche par points (partielle) - pêche complète

L’ensemble des données sont stocké dans la base de données aspe de l’OFB.

Les tendances temporelles seront estimées par une analyse statistique de Modèles Additifs Généralisé Mixtes (GAMMs), de Modèles Linéaires Généralisé Mixtes (GLMM) ou de Zero-Inflated Poisson GLM (ZIPGLMM) en fonction de l’ajustement des modèles aux données.

2.2 Choix des variables

Pour répondre à la question initiale, les modèles vont se construire à partir des variables suivantes :

La variable réponse : valeur, effectif des poissons variable quantitative discrète bornée

Les variables explicatives :

les co-variables :

annee, variable quantitative discrète bornée
ope_surface_calculee, la surface pêchée variable quantitative continue bornée

les facteurs :

pro_libelle, les protocoles utilisés (2 modalitées) variable qualitative nominal
sta_id, les stations d’échantillonages (163 modalitées) variable qualitative nominal

2.3 Demarche des analyses statistiques

Comme les données sont appariées, un facteur aléatoire (1|sta_id) sera intégré dans tous les modèles.Les modèles seront de la forme suivante :

modele <- effectif ~ annee + protocole + surface + (1|station)

Le coefficient prédit par les modèles permettra d’estimer la tendance temporelle.

Modèles Additifs Généralisé Mixtes (GAMMs)

Les GAMMs sont utilisés en considérant que les relations entre les variables ne sont pas forcément linéaires. Par exemple, les effectifs fluctuent entre les années par un échantillonnage bi-annelles pour certaines stations et une grande variabilité annuelle des effectifs conditionnés par les conditions lors de la période de reproduction.

Modèles Linéaires Généralisé Mixtes (GLMM)

Les GLMM sont utilisés pour améliorer l’ajustement du modèle au données.

Modèle de zéro enflées

Les modèles Zero-inflated sont utilisés pour corrigé la sur-dispersion des résidus

2.4 Représentation des series temporelles par stations (163)

ope_selection %>% # Représentation graphique des stations sélectionnées et des types de pêche associées
  ggplot(
    aes(
      x = as.character(pop_libelle),
      y = annee,
      fill = pro_libelle
    ),
    legend.background = element_rect(fill = "#ffffff")
  ) +
  scale_fill_manual(values = pal) +
  geom_tile() +
  labs(
    title = "Les points de prélèvements sélectionnés et les types de pêches associées",
    subtitle = "Région Bourgogne-Franche-Comté",
    x = "",
    y = "Années",
    fill = "Type de pêche"
  ) +
  theme_light(base_size = 18) +
  theme(
    panel.grid.major = element_line(color = "#ffffff", size = 0.1),
    panel.grid.minor = element_line(color = "#ffffff"),
    panel.background = element_rect(fill = "#faf0e6"),
    legend.position = "bottom",
    legend.text = element_text(size = 18),
    # taille labels pro_libelle
    legend.title = element_text(size = 15),
    # taille titre légende
    axis.text.y = element_text(size = 16) # taille noms stations
  ) +
  coord_flip() +
  scale_y_continuous(breaks = seq(min(ope_selection$annee), max(ope_selection$annee), by = 3))

2.5 Représentation de l’effort d’echantillonnage par stations

mapview(
  list(buffer, mes_depts),
  # Visualisation de l'effort de pêche
  layer.name = c("Buffer_1_km", "Bourgogne-Franche-Comté"),
  col.regions = list("darkseagreen4", "darkseagreen1")
) +
  mapview(
    sta_id_ope_selection,
    layer.name = "Nombre operation par station",
    # size_var = "nb_ope_par_sta",
    # cex = "sta_id",
    cex = "nb_ope_par_sta",
    lwd = 0.2,
    zcol = "nb_ope_par_sta",
   # col.regions = "#D693E0",
    alpha = 0.8,
    legend = TRUE,
    labels = list(text = sta_id_ope_selection$nb_ope_par_sta, cex = 1.5)
  )

3 Exemple d’analyse avec une espèce de poisson

3.1 Resprésenation graphique - effectif total en fonction des années

# 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)

esp_ope_selection_graph <- esp_ope_selection %>%
  mef_ajouter_type_protocole() %>%
  mef_ajouter_ope_date() %>%
  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),
    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),
    indicateur = "effectif_total"
  ) %>% # creation d'une colonne pour la période (ex : 2007-2008)
  mutate(periode = {
    debut <- floor((annee - annee_debut_periode) / 2) * 2 + annee_debut_periode
    paste0(debut)
  }, periode = as.numeric(periode)) %>%
  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,
    periode,
    x,
    y,
    ope_date, 
    julian
  ) 
# Creation du df par especes avec effectif total - periode
total_valeur_esp <- esp_ope_selection_graph %>%
  group_by(periode, esp_code_alternatif, indicateur) %>%
  summarise(effectif_total = sum(valeur), .groups = "drop") %>%
  select(esp_code_alternatif, periode, effectif_total, indicateur)

# Creation du df par especes avec effectif total - annee
total_valeur_esp_an <- esp_ope_selection_graph %>%
  group_by(annee, esp_code_alternatif, indicateur) %>%
  summarise(effectif_total = sum(valeur), .groups = "drop") %>%
  select(esp_code_alternatif, annee, effectif_total, indicateur)

Selection de l’espèce à traiter : ici la Truite fario “TRF”

# Selection de lespèce cible pour le graphique 
total_valeur_esp_trf <- total_valeur_esp %>%
  filter(esp_code_alternatif == "TRF") 

# Selection de lespèce cible pour le graphique 
total_valeur_esp_trf_an <- total_valeur_esp_an %>%
  filter(esp_code_alternatif == "TRF") 

# Selecetion de l'espèce pour les analyses
esp_ope_selection_trf <- esp_ope_selection_graph %>%
  filter(esp_code_alternatif == "TRF")

plot_trf_annee <- ggplot(total_valeur_esp_trf_an, aes(x = annee, y = effectif_total)) +
  facet_wrap2(
    esp_code_alternatif ~ indicateur,
    ncol = 3,
    scales = "free_y",
    axes = "all",
    trim_blank = T
  ) +
  geom_point(shape = 21, size = 2) +
  geom_line(linetype = 3, linewidth = 0.5) +
  labs(x = "Années", y = "Effectif total") +
  geom_smooth(method = "lm", color = "darkred") +
  # geom_smooth(method = "loess",
  #             se = F,
  #             color = "darkred") +  # Add loess trend lines
  theme(
    strip.background = element_rect(color = "grey85", fill = "grey85"),
    legend.position = "bottom",
    panel.background = element_rect(fill = "grey95"),
    panel.grid.major = element_line(color = "#faffff", size = 0.1),
    panel.grid.minor = element_line(color = "#faffff")
  )

plot_trf_periode <- ggplot(total_valeur_esp_trf, aes(x = periode, y = effectif_total)) +
  facet_wrap2(
    esp_code_alternatif ~ indicateur,
    ncol = 3,
    scales = "free_y",
    axes = "all",
    trim_blank = T
  ) +
  geom_point(shape = 21, size = 2) +
  geom_line(linetype = 3, linewidth = 0.5) +
  labs(x = "Période", y = "Effectif total") +
  geom_smooth(method = "lm", color = "darkred") +
  # geom_smooth(method = "loess",
  #             se = F,
  #             color = "darkred") +  # Add loess trend lines
  theme(
    strip.background = element_rect(color = "grey85", fill = "grey85"),
    legend.position = "bottom",
    panel.background = element_rect(fill = "grey95"),
    panel.grid.major = element_line(color = "#faffff", size = 0.1),
    panel.grid.minor = element_line(color = "#faffff")
  )

plot_trf <- ggarrange(
  plot_trf_annee,
  plot_trf_periode,
  ncol = 2,
  nrow = 1,
  labels = c(
    "Effectif en fonction de l'annee",
    "Effectif en fonction de la periode (n : n+1)"
  ),
  font.label = list(
    size = 14,
    color = "darkred",
    face = "bold",
    family = NULL
  ),
  hjust = -0.5,
  vjust = 1.9,
  align = c("none", "h", "v", "hv")
)
plot_trf

3.2 Construction du modele GAMM

mod_gamm_1.0 <- gamm(
  valeur ~ s(annee, k = 15) + s(ope_surface_calculee, k = 15)  + pro_libelle,
  random = list(sta_id =  ~ 1),
  data = esp_ope_selection_trf,
  family = poisson(link = "log"),
  method = "REML"
)

## 
##  Maximum number of PQL iterations:  20

mod_gamm_1.1 <- gamm(
  valeur ~ s(annee, k = 15) + s(ope_surface_calculee, k = 15)  + pro_libelle,
  random = list(sta_id =  ~ 1),
  data = esp_ope_selection_trf,
  family = poisson(link = "identity"),
  method = "REML"
)

## 
##  Maximum number of PQL iterations:  20

mod_gamm_1.2 <- gamm(
  valeur ~ s(annee, k = 15) + s(ope_surface_calculee, k = 15)  + pro_libelle,
  random = list(sta_id =  ~ 1),
  data = esp_ope_selection_trf,
  family = poisson(link = "sqrt"),
  method = "REML"
)

## 
##  Maximum number of PQL iterations:  20

mod_gamm_1.3 <- gam(
  valeur ~ s(annee) + s(ope_surface_calculee)  + pro_libelle+ s(x,y),
  data = esp_ope_selection_trf,
  family = poisson(link = "log"),
  method = "REML"
)

3.2.1 Autocorelation temporelle - plot

acf_mod_gamm_1.0 <- as.ggplot(~ acf(residuals(mod_gamm_1.0$lme, type = "normalized"), lag.max = 20))
acf_mod_gamm_1.1 <- as.ggplot(~ acf(residuals(mod_gamm_1.1$lme, type = "normalized"), lag.max = 20))
acf_mod_gamm_1.2 <- as.ggplot(~ acf(residuals(mod_gamm_1.2$lme, type = "normalized"), lag.max = 20))

ggarrange(
  acf_mod_gamm_1.0,
  acf_mod_gamm_1.1,
  acf_mod_gamm_1.2,
  ncol = 3,
  nrow = 1,
  labels = c("acf_mod_gamm_1.0", "acf_mod_gamm_1.1", "acf_mod_gamm_1.2"),
  font.label = list(
    size = 14,
    color = "darkred",
    face = "bold",
    family = NULL
  ),
  hjust = -0.1,
  vjust = 1.6,
  align = c("none", "h", "v", "hv")
)

3.2.2 Autocorelation spatial - semi-variogramme

ggplot(esp_ope_selection_trf, aes(x = x, y = y, size = valeur)) +
    geom_point() +
    coord_fixed()

ggplot(esp_ope_selection_trf, aes(x = annee, y = valeur)) +
  geom_point()

var <- variogram(valeur ~ annee + ope_surface_calculee, locations = ~ x+y, data = esp_ope_selection_trf)
var

##         np        dist     gamma dir.hor dir.ver   id
## 1  4129839    453.4066 147.71964       0       0 var1
## 2  2762738  12145.1806 131.79681       0       0 var1
## 3  4770940  21150.2871 112.84181       0       0 var1
## 4  3535669  29747.7735 205.35159       0       0 var1
## 5  4389037  37930.6049 153.84611       0       0 var1
## 6  5396952  46010.7995 213.62465       0       0 var1
## 7  4486533  55613.6147 169.80796       0       0 var1
## 8  4881343  63031.0201  68.27284       0       0 var1
## 9  4928812  71676.2461  99.59569       0       0 var1
## 10 3457271  79701.0882  88.95981       0       0 var1
## 11 3521340  87739.5354  61.03049       0       0 var1
## 12 6041220  95468.6622 145.60051       0       0 var1
## 13 5885089 104588.3740 102.31949       0       0 var1
## 14 5597444 113009.1044  85.78694       0       0 var1
## 15 5834463 121683.9966  47.35979       0       0 var1

plot(var, col= "black")

show.vgms(models = c("Exp",  # Exponential model
                     "Sph",  # Spherical model
                     "Gau",  # Gaussian model
                     "Nug"), # Nugget model
          nugget = 50)

vfit <- fit.variogram(var, vgm(c("Exp", "Gau", "Sph")))
vfit

##   model    psill    range
## 1   Nug 148.1882     0.00
## 2   Sph   0.0000 40561.33

# Did the model-fitting process converge successfully?
# This is not guaranteed for nonlinear models.
attributes(vfit)$singular

## [1] TRUE

# For convenience, use show.vgms() to calculate the
# predicted values for our fitted model
fitted <-show.vgms(min = min(var$dist), 
                   max = max(var$dist), 
                   models ="Sph",
                   range = vfit$range[2],
                   sill = vfit$psill[2],
                   nugget = vfit$psill[1],
                   plot = FALSE)

# Now we can easily plot this against the observed data
ggplot(var, aes(dist, gamma)) + 
  geom_point() + 
  geom_line(data = fitted, 
            aes(distance, semivariance),
            linetype = 2)

# plot(var, vfit, col = "black")

3.2.3 Indépendance du modèles aux data

appraise(mod_gamm_1.0$gam)

appraise(mod_gamm_1.1$gam)

appraise(mod_gamm_1.2$gam)

k.check(mod_gamm_1.0$gam)

##                         k'      edf   k-index p-value
## s(annee)                14 13.90401 0.9991631  0.4300
## s(ope_surface_calculee) 14 13.53667 0.9161402  0.0025

k.check(mod_gamm_1.1$gam)

##                         k'      edf  k-index p-value
## s(annee)                14 13.88207 1.006524  0.7150
## s(ope_surface_calculee) 14 12.04178 0.964481  0.0375

k.check(mod_gamm_1.2$gam)

##                         k'      edf   k-index p-value
## s(annee)                14 13.90506 0.9401616   0.000
## s(ope_surface_calculee) 14 12.88747 0.9591719   0.005

-> Verifier le k’/edf : edf doit être << k’. Si edf est proche de k’, cela signifie qu’il faut augmenter k. Cependant, comme on travaille sur 20 ans on ne peut pas dépasser k = 19.

3.2.4 Verifier la concurvity

concurvity(mod_gamm_1.0$gam)

##               para   s(annee) s(ope_surface_calculee)
## worst    0.5097488 0.12571435              0.63851587
## observed 0.5097488 0.03634524              0.05918667
## estimate 0.5097488 0.06501472              0.11567597

concurvity(mod_gamm_1.1$gam)

##               para   s(annee) s(ope_surface_calculee)
## worst    0.5097488 0.12571435              0.63851587
## observed 0.5097488 0.03486796              0.07327614
## estimate 0.5097488 0.06501472              0.11567597

concurvity(mod_gamm_1.2$gam)

##               para   s(annee) s(ope_surface_calculee)
## worst    0.5097488 0.12571435              0.63851587
## observed 0.5097488 0.03557414              0.06114753
## estimate 0.5097488 0.06501472              0.11567597

3.2.5 Représentation graphique des smooth

draw(mod_gamm_1.0, residuals = T)

draw(mod_gamm_1.1, residuals = T)

draw(mod_gamm_1.2, residuals = T)

draw(mod_gamm_1.0, residuals = F)

draw(mod_gamm_1.1, residuals = F)

draw(mod_gamm_1.2, residuals = F)

3.3 Construction du modele GLMM

summary(esp_ope_selection_graph)

##      ope_id          valeur         esp_code_alternatif     annee     
##  Min.   :11207   Min.   :   1.000   Length:142523       Min.   :2007  
##  1st Qu.:36802   1st Qu.:   1.000   Class :character    1st Qu.:2011  
##  Median :38261   Median :   1.000   Mode  :character    Median :2016  
##  Mean   :53164   Mean   :   7.957                       Mean   :2016  
##  3rd Qu.:84312   3rd Qu.:   1.000                       3rd Qu.:2020  
##  Max.   :96052   Max.   :3511.000                       Max.   :2025  
##                                                                       
##  ope_surface_calculee                                    pro_libelle   
##  Min.   :  52.0       Pêche complète à un ou plusieurs passages:55006  
##  1st Qu.: 937.5       Pêche partielle par points (grand milieu):87517  
##  Median : 937.5                                                        
##  Mean   : 983.8                                                        
##  3rd Qu.:1250.0                                                        
##  Max.   :2783.0                                                        
##                                                                        
##       pop_id           sta_id        indicateur           periode    
##  33562   :  2829   8884   :  2829   Length:142523      Min.   :2007  
##  35119   :  2525   9166   :  2525   Class :character   1st Qu.:2011  
##  69061   :  2250   19535  :  2250   Mode  :character   Median :2015  
##  14895981:  2058   4331694:  2058                      Mean   :2015  
##  67287   :  2035   18756  :  2035                      3rd Qu.:2019  
##  81644   :  2032   24171  :  2032                      Max.   :2025  
##  (Other) :128794   (Other):128794                                    
##        x                y              ope_date                       julian   
##  Min.   :691396   Min.   :6566885   Min.   :2007-05-31 09:30:01   Min.   :102  
##  1st Qu.:774154   1st Qu.:6638880   1st Qu.:2011-08-31 13:30:01   1st Qu.:180  
##  Median :857024   Median :6679085   Median :2016-07-07 14:00:01   Median :240  
##  Mean   :847751   Mean   :6678175   Mean   :2016-04-12 21:58:20   Mean   :225  
##  3rd Qu.:918183   3rd Qu.:6718963   3rd Qu.:2020-07-03 09:00:01   3rd Qu.:266  
##  Max.   :998675   Max.   :6798620   Max.   :2025-10-24 08:15:01   Max.   :314  
##

mod_glmm_1.0 <- glmer(
  valeur ~ annee + ope_surface_calculee + pro_libelle  + (1 |sta_id),
  data = esp_ope_selection_trf,
  family = poisson(link = "log")
)
mod_glmm_1.1 <- glmer(
  valeur ~ annee + ope_surface_calculee + pro_libelle  + (1 |sta_id),
  data = esp_ope_selection_trf,
  family = poisson(link = "identity")
)
mod_glmm_1.2 <- glmer(
  valeur ~ annee + ope_surface_calculee + pro_libelle  + (1 |sta_id),
  data = esp_ope_selection_trf,
  family = poisson(link = "sqrt")
)

mod_glmm_1.3 <- glmer(
  valeur ~ annee + ope_surface_calculee + pro_libelle + julian + (1 |sta_id),
  data = esp_ope_selection_trf,
  family = poisson(link = "log")
)

3.3.1 Indépendance du modèle aux données

res_mod_glmm_1.0 <- as.ggplot(~ plotresid(mod_glmm_1.0))
res_mod_glmm_1.1 <- as.ggplot(~ plotresid(mod_glmm_1.1))
res_mod_glmm_1.2 <- as.ggplot(~ plotresid(mod_glmm_1.2))

ggarrange(
  res_mod_glmm_1.0,
  res_mod_glmm_1.1,
  res_mod_glmm_1.2,
  ncol = 2,
  nrow = 2,
  labels = c("log", "identity", "sqrt"),
  font.label = list(
    size = 13,
    color = "darkred",
    face = "bold",
    family = NULL
  ),
  hjust = -0.1,
  vjust = 1.6,
  align = c("none", "h", "v", "hv")
)

3.3.2 Homoscédasticité des résidus

overdisp.glmer(mod_glmm_1.0)

## Residual deviance: 72030.05 on 16123 degrees of freedom (ratio: 4.468)

overdisp.glmer(mod_glmm_1.1)

## Residual deviance: 72097.968 on 16123 degrees of freedom (ratio: 4.472)

overdisp.glmer(mod_glmm_1.2)

## Residual deviance: 72067.698 on 16123 degrees of freedom (ratio: 4.47)

3.3.3 Test post-hoc - fonction anova

Si les conditions d’applications du modèle sont bien respectées :

# Anova(mod_glmm_1.0)
# Anova(mod_glmm_1.1)
# Anova(mod_glmm_1.2)

# summary(mod_glmm_1.0)
# summary(mod_glmm_1.1)
# summary(mod_glmm_1.2)

Si problème de sur ou sous-dispersion : esssayer une Quasi binomiale négative.

3.4 Construction du modèle GLMM avec une Quasi binomiale négative

# mod_glmm_qbn <- glmer.nb(valeur ~ annee + ope_surface_calculee + pro_libelle+ (1|sta_id),  data = esp_ope_selection_trf,  family = poisson(link = "log")
# )

3.5 Construction des modeles GLMM Zero-inlfated

mod_zi_1.0 <- glmmTMB(
  valeur ~ annee + ope_surface_calculee + pro_libelle + (1 |
                                                           sta_id) ,
  data = esp_ope_selection_trf,
  family = poisson(link = "log"),
  ziformula = ~ 1
)
mod_zi_1.1 <- glmmTMB(
  valeur ~ annee + ope_surface_calculee + pro_libelle + (1 | sta_id) ,
  family = nbinom2(link = "log"),
  data = esp_ope_selection_trf,
  ziformula = ~ 1
)
mod_zi_1.2 <- glmmTMB(
  valeur ~ annee + ope_surface_calculee + pro_libelle + (1 | sta_id) ,
  family = truncated_poisson(link = "log"),
  data = esp_ope_selection_trf,
  ziformula = ~ 1
)
mod_zi_1.3 <- glmmTMB(
  valeur ~ annee + ope_surface_calculee + pro_libelle + (1 | sta_id) ,
  family = truncated_nbinom2(link = "log"),
  data = esp_ope_selection_trf,
  ziformula = ~ 1
)

On selectionne à l’aide du critère AIC le meilleur modèle c’est à dire le modèle avec le plus petit AIC.

AIC(mod_zi_1.0, mod_zi_1.1, mod_zi_1.2, mod_zi_1.3)

##            df       AIC
## mod_zi_1.0  6 106394.25
## mod_zi_1.1  7  57702.92
## mod_zi_1.2  6  93880.00
## mod_zi_1.3  7  25950.06

res_mod_zi_1.0 <- as.ggplot(~ plotresid(mod_zi_1.0))
res_mod_zi_1.1 <- as.ggplot(~ plotresid(mod_zi_1.1))
res_mod_zi_1.2 <- as.ggplot(~ plotresid(mod_zi_1.2))
res_mod_zi_1.3 <- as.ggplot(~ plotresid(mod_zi_1.3))

ggarrange(
  res_mod_zi_1.0,
  res_mod_zi_1.1,
  res_mod_zi_1.2,
  res_mod_zi_1.3,
  ncol = 2,
  nrow = 2,
  labels = c("poisson", "nbinom2", "truncated_poisson", "truncated_nbinom2"),
  font.label = list(
    size = 13,
    color = "darkred",
    face = "bold",
    family = NULL
  ),
  hjust = -0.1,
  vjust = 1.6,
  align = c("none", "h", "v", "hv")
)

summary(mod_zi_1.0)

##  Family: poisson  ( log )
## Formula:          
## valeur ~ annee + ope_surface_calculee + pro_libelle + (1 | sta_id)
## Zero inflation:          ~1
## Data: esp_ope_selection_trf
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##  106394.2  106440.4  -53191.1  106382.2     16122 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  sta_id (Intercept) 0.1392   0.3731  
## Number of obs: 16128, groups:  sta_id, 100
## 
## Conditional model:
##                                                        Estimate Std. Error
## (Intercept)                                          -1.488e+01  2.331e+00
## annee                                                 7.726e-03  1.157e-03
## ope_surface_calculee                                 -3.717e-04  3.449e-05
## pro_libellePêche partielle par points (grand milieu) -8.347e-02  8.626e-02
##                                                      z value Pr(>|z|)    
## (Intercept)                                           -6.382 1.75e-10 ***
## annee                                                  6.679 2.40e-11 ***
## ope_surface_calculee                                 -10.777  < 2e-16 ***
## pro_libellePêche partielle par points (grand milieu)  -0.968    0.333    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -21.38     345.67  -0.062    0.951

summary(mod_zi_1.1)

##  Family: nbinom2  ( log )
## Formula:          
## valeur ~ annee + ope_surface_calculee + pro_libelle + (1 | sta_id)
## Zero inflation:          ~1
## Data: esp_ope_selection_trf
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##   57702.9   57756.7  -28844.5   57688.9     16121 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  sta_id (Intercept) 0.1296   0.36    
## Number of obs: 16128, groups:  sta_id, 100
## 
## Dispersion parameter for nbinom2 family ():  1.3 
## 
## Conditional model:
##                                                        Estimate Std. Error
## (Intercept)                                          -2.016e+01  4.677e+00
## annee                                                 1.031e-02  2.330e-03
## ope_surface_calculee                                 -2.233e-04  6.230e-05
## pro_libellePêche partielle par points (grand milieu) -1.271e-01  9.028e-02
##                                                      z value Pr(>|z|)    
## (Intercept)                                           -4.311 1.63e-05 ***
## annee                                                  4.424 9.67e-06 ***
## ope_surface_calculee                                  -3.584 0.000338 ***
## pro_libellePêche partielle par points (grand milieu)  -1.408 0.159223    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -22.49     603.53  -0.037     0.97

summary(mod_zi_1.2)

##  Family: truncated_poisson  ( log )
## Formula:          
## valeur ~ annee + ope_surface_calculee + pro_libelle + (1 | sta_id)
## Zero inflation:          ~1
## Data: esp_ope_selection_trf
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##   93880.0   93926.1  -46934.0   93868.0     16122 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  sta_id (Intercept) 1.567    1.252   
## Number of obs: 16128, groups:  sta_id, 100
## 
## Conditional model:
##                                                        Estimate Std. Error
## (Intercept)                                          -2.983e+01  4.645e+00
## annee                                                 1.462e-02  2.280e-03
## ope_surface_calculee                                 -5.330e-04  4.205e-05
## pro_libellePêche partielle par points (grand milieu) -2.358e-01  2.845e-01
##                                                      z value Pr(>|z|)    
## (Intercept)                                           -6.421 1.35e-10 ***
## annee                                                  6.413 1.42e-10 ***
## ope_surface_calculee                                 -12.675  < 2e-16 ***
## pro_libellePêche partielle par points (grand milieu)  -0.829    0.407    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -21.98     467.81  -0.047    0.963

summary(mod_zi_1.3)

##  Family: truncated_nbinom2  ( log )
## Formula:          
## valeur ~ annee + ope_surface_calculee + pro_libelle + (1 | sta_id)
## Zero inflation:          ~1
## Data: esp_ope_selection_trf
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##   25950.1   26003.9  -12968.0   25936.1     16121 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  sta_id (Intercept) 2.099    1.449   
## Number of obs: 16128, groups:  sta_id, 100
## 
## Dispersion parameter for truncated_nbinom2 family (): 1.42e-10 
## 
## Conditional model:
##                                                        Estimate Std. Error
## (Intercept)                                          -1.117e+02  1.143e+03
## annee                                                 4.413e-02  4.071e-03
## ope_surface_calculee                                 -5.982e-04  1.419e-04
## pro_libellePêche partielle par points (grand milieu) -3.033e-01  3.394e-01
##                                                      z value Pr(>|z|)    
## (Intercept)                                           -0.098    0.922    
## annee                                                 10.840  < 2e-16 ***
## ope_surface_calculee                                  -4.216 2.48e-05 ***
## pro_libellePêche partielle par points (grand milieu)  -0.894    0.371    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -24.53    1671.66  -0.015    0.988