In polehalieutique/demerstem: This package aims to falicitate Fish stock assessement for west african partners involved in the EU DEMERSTEM project

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

#library(devtools)
#detach("package:demerstem", unload = TRUE)
#install_github("polehalieutique/demerstem")

# Il se peut que vous ayez à installer TropFishR (estimation de paramètres biologiques) et Jags (simulation bayésienne).
# remotes::install_github("tokami/TropFishR") # Depuis le 01/01/2024 le package TropFishR n'est plus disponible sur le CRAN, il faut passer par le git
# https://sourceforge.net/projects/mcmc-jags/files/JAGS/4.x/Windows/ # Pour Jags

library(demerstem)
packageVersion("demerstem")

data(data_IA)
data(captures_aggregees)

Cet script vise à accompagner les utilisateurs du package demerstem dans la réalisation d'une évaluation de stock avec un modèle de production selon un ajustement en pseudo-équilibre.

1) Data manipulation / Manipulation des données

Nous utilisons ici les données issues de la pêcherie de Thiof dans la zone Sénégal-Gambie-Mauritanie entre 1974 et 2019. Les tableaux suivants fournissent les Captures et les Indices d'abondance sur cette période.

plot_IA <- ggplot(data = data_IA %>% mutate(title = "Abundance indices"), aes(x=year, y=AI_cor)) + 
   geom_line(size = 1.01) + 
   theme_nice() + 
   facet_grid(~title) +
  scale_color_manual(values = c("darkblue", "blue")) +  ylab("IA") +
  xlab("Year")

 Plot_Catch <- rbind(captures_aggregees %>% filter(code_pays != "MRT") %>%  group_by(year) %>%  summarise(Catch = sum(c)) %>% mutate(code_pays = "SEN-GMB"), captures_aggregees %>% filter(code_pays == "MRT")%>%  dplyr::rename(Catch = c) %>% dplyr::select(-espece))  %>% 
  #complete(code_pays, year = 1971:1989, fill = list(Catch = 0)) %>%
  mutate(title = "Catch times series") %>% 
  ggplot(aes(x = year, y = Catch, fill = code_pays), show.legend = FALSE) + 
  #geom_line(size = 1.1, col = 'black') + 
  scale_fill_manual(values = c("grey50", "grey30")) +
  geom_area(show.legend = T) + 
  geom_vline(xintercept = c(1974, 2019), linetype = "dashed") + 
  ylab("Catch (t)") +
  xlab("Year") +
  theme_nice() +
  facet_grid(~title) + #%>% filter(code_pays =="GIN")
theme(legend.title = element_text(size=11),
            legend.text = element_text(size=11)) + 
  theme(legend.position = c(0.20, 0.80))

• Captures

data_captures <- captures_aggregees %>%  group_by(year) %>%  summarise(Catch = sum(c))

data_captures

• Indices d'abondances

data_IA

Il sont représentés dans la figure suivante :

plot_grid(plotlist =list(Plot_Catch, plot_IA))

2) Calcul des Efforts de Fox

• La fonction fox_effort_calculation() va calculer les efforts de Fox (cf formule), pour ajuster les modèles sous une hypothèse de pseudo-équilibre (on pondère en considérant qu'une cohorte subit plusieurs efforts au cours de sa vie, et on donne plus de poids aux efforts récents)

Efox,i = (k x Ei + (k-1) x Ei-1 + ... + Ei-k+1)/(k+(k-1)+...1)

• Elle prend en argument :

• La série d'IA moyen calculé avec mean_ai()

• Les captures totales

• k : nombre de groupes d'âges exploités (1, 3 ou 5)

• Structure des arguments :

• Pour data_IA, Year et 1 seule série d'IA doivent être sélectionnés (dans cet ordre)

• Pour catch_table, Year et Catch, dans cet ordre.

• Les deux séries de données doivent être de même tailles

#----------------------------
#      Calcul des Efox 
#----------------------------

k <- 5
table_Efox_thiof <- fox_effort_calculation(data = data_IA, catch_table = data_captures, k)

# Effort MRT
Table_effort <- data_captures %>%  dplyr::rename(Year = 'year') %>% 
  full_join(table_Efox_thiof %>% dplyr::select(Year,E), by = 'Year') %>%  dplyr::rename(Effort = 'E') 


Table_effort_2 <- Table_effort %>% mutate(title = "Theoretical effort time series - E. aeneus (Thiof) \n - MRT & SEN & GMB -") %>% dplyr::select(Year, Effort, title) %>% pivot_longer(cols = c(2), names_to = "Study_case", values_to = "Effort") %>% drop_na() 

plot_Effort <- ggplot(data = Table_effort_2, aes(x=Year, y=Effort, color = Study_case, linetype = Study_case)) + 
   geom_line(size = 1.01) + 
   theme_nice() + 
   facet_grid(~title) +
  scale_color_manual(values = c("darkblue", "blue")) 

plot_Effort

3) Application d'un modèle de production

3.1) Fonction de production de Fox

• La fonction fox_model() permet ensuite d'ajuster un modèle global de Fox, on a :

Ue(E) = a * exp(-b * E)

On cherche à ajuster cette relation sur nos données (U = PUE = IA, par rapport à l'Effort effectif = Efox) en testant toutes les valeurs possibles des paramètres a et b pour lesquels on obtient une estimation quand le modèle s'ajuste le mieux aux données. On peut alors reconstruire la courbe des rendements : Ye(E) = E * Be(E)

• En argument :

• table_Efox : table issue de la fonction fox_effort_calculation()

• graph_param : des paramètres graphiques comprenant la taille de l'axe des abscisses et du titre

• logarithmic : par défaut = TRUE. Transformation log (par défaut = TRUE, si ne fonctionne pas essayer log=FALSE).

• a_start : l'ajustement peut être sensible aux valeurs d'init choisis

• b_start : idem

mE_range <- 1.5
title <- "Fox model - Thiof"
graph_param <- c(mE_range, title) #param de mise en forme du graph
table_Efox <- table_Efox_thiof
Fox_model_thiof <- fox_model(table_Efox, graph_param, log = T, a_start = 1, b_start = 100)

Fox_model_thiof[[1]]

plot_grid(plotlist = list(Fox_model_thiof[[2]][[1]], Fox_model_thiof[[2]][[3]], Fox_model_thiof[[2]][[2]], Fox_model_thiof[[2]][[4]]))
plot_grid(plotlist = list(Fox_model_thiof[[2]][[5]]))

3.2) Fonction de production Pella-Tomlinson (modèle généralisé)

Le modèle généralisé : IA ~ (a+b*Efox)^(1/m-1). Le fonctionnement est semblable à la foonction avec le modèle de Fox. On ajoute néanmoins un argument supplémentaire, les estimations sont donc plus sensibles aux paramètres d'initialisation.

#----------------------------
#      Modèle de pella-tomlinson 
#----------------------------

lengthEfox <- 1.5
title <- "Generalised model - Thiof"
graph_param <- c(lengthEfox, title)

a_init = 1/sqrt(table_Efox$IA[1])
b_init = 1/sqrt(tail(table_Efox$IA,1)) - 1/(sqrt(table_Efox$IA[1])) 

a <- c(0.0000005, a_init, 0.005)
b <- c(0.0000001, b_init, 0.01)
m <- c(0.5, 0.7, 0.9)

warning_control=T
table_Efox <- table_Efox_thiof
Gene_model_thiof <- generalised_model(table_Efox, graph_param, a, b, m, warning_control, IC = F)

Gene_model_thiof[[1]]

plot_grid(plotlist = list(Gene_model_thiof[[2]][[1]], Gene_model_thiof[[2]][[3]], Gene_model_thiof[[2]][[2]], Gene_model_thiof[[2]][[4]]))
plot_grid(plotlist = list(Gene_model_thiof[[2]][[5]]))

4) Conclusion/Résumé

L'ajustement est réalisé sur la relation entre les $IA$ et $mE_{fox}$ qui permet d'analyser la réponse du stock aux variations de la pression de pêche depuis les années 75 (début de la disponibilité des indices d'abondance). Il ressort assez clairement qu'une augmentation de l'effort de pêche entraîne une forte diminution de l'abondance du stock selon une relation non linéaire. Le modèle de Fox ne semble pas adapté pour expliquer les fortes abondances du début de la série temporelle. Ainsi, sur la base des AIC, le modèle généralisé (ou de Pella-Tomlinson) sera ainsi privilégié pour réaliser le diagnostic. Néanmoins, ce modèle a également du mal à expliquer correctement les fortes abondances observées à la fin des années 70. De même, les années les plus récentes sont très mal expliquées. Cet écart entre les prédictions et les observations est lié aux très fortes captures réalisées en Mauritanie entre 2016 et 2018. Nous pouvons identifier une histoire de l'exploitation du stock qui aurait débutée dans les années 1970 jusqu'à atteindre une première phase de surexploitation au début des années 2000. L'effort de pêche a ensuite diminué durant une dizaine d'année avant de repartir à la hausse et d'exploser dans la fin des années 2010, menant à un état de surexploitation très avancé.

# Tableau 2 : points de référence
set_flextable_defaults(font.family = "Palatino Linotype")
Tableau_indicateurs <- data.frame("Scénario" = c(NA, "", ""), 
                                  "Modèle" = c(NA,"Fox", "Pella-Tomlinson"),
                                     "Indicateurs de l'état du stock" = c("MSY", 
                                                                        round(as.numeric(Fox_model_thiof[[1]][,3][2])),
                                                                        round(as.numeric(Gene_model_thiof[[1]][,4][2]))), 
                                     "Indicateurs de l'état du stock" = c("E_MSY",
                                                                          round(1/as.numeric(Fox_model_thiof[[1]][,4][2]), 2),
                                                                          round(1/as.numeric(Gene_model_thiof[[1]][,5][2]), 2)),
                                     "Indicateurs de l'état du stock" = c("B/B0",
                                                                          as.numeric(Fox_model_thiof[[1]][,6][2]),
                                                                          round(as.numeric(Gene_model_thiof[[1]][,7][2]),2)),
                                     "Indicateurs de l'état du stock" = c("E/E_MSY", 
                                                                          as.numeric(Fox_model_thiof[[1]][,4][2]),
                                                                          round(as.numeric(Gene_model_thiof[[1]][,5][2]),2)))


names(Tableau_indicateurs) <- c("Scénario", "Modèle", "Indicateurs de l'état du stock", "A_________", "B_________", "D___________")
ft <- flextable(Tableau_indicateurs) #, Fox_model_combinedimist[[1]], Fox_model_commercialsimist_log[[1]], Fox_model_commercialsimist[[1]]))
ft <- ft %>% set_caption("Tableau 2 - Indicateurs de l’état du stock selon les scénarios et différents modèles.", style = "Table Caption") %>% align(align = "center", part = "all")%>% padding(padding = 0, part = "body")%>% autofit() %>% fit_to_width(7.5)

ft <- bg(ft, bg = "#C5E0B3", part = "header")
ft <- bold(ft, part = "header", bold = T)
ft <- ft %>% merge_at(i=1, j=c(3:6), part = "head") %>% hline(i = 1, part = "header", border = fp_border(color="black", width = 1)) %>%  hline_top(part= "header", border = fp_border(color="black", width = 1)) %>% align(j = 1, align = "left", part ="all") %>% hline(i = 3, part = "body", border = fp_border(color="black", width = 1)) %>% bold(i = 1, j = c(3:6), part = "body", bold = T)
ft

polehalieutique/demerstem documentation built on Aug. 4, 2024, 5:12 a.m.