In pbs-assess/gfsynopsis: Data Synopsis Reports for PBS Groundfish

NORMALISATION DE L’INDICE DES CAPTURES PAR UNITÉ D’EFFORT (CPUE) {#app:cpue-models}

Nous avons cherché à produire un indice d’abondance à partir des données sur les CPUE dans la pêche commerciale au chalut, qui ont été normalisées pour la profondeur, l’emplacement de la pêche (régions spatiales définies; figure \@ref(fig:cpue-locality-map)), le mois, le navire et la latitude. Avant d’ajuster un modèle de normalisation, nous avons dû filtrer et manipuler les données disponibles sur les prises et l’effort pour générer un ensemble de données approprié pour l’ajustement du modèle. Dans les sections suivantes, nous décrivons ces décisions pour les données de 1996 à 2017 et présentons ensuite notre modèle de normalisation de l’indice. Ce modèle et la description qui suit reposent largement sur une évaluation récente de la morue du Pacifique (en partie copiée mot pour mot), où ce modèle a été élaboré et appliqué pour la première fois [@dfo2019pcod3cd5abcd].

DÉFINITION DE LA FLOTTE DE 1996 À `r as.numeric(format(Sys.Date(), "%Y")) - 1`

Les données sur la pêche commerciale du poisson de fond au chalut de fond de 1996 à nos jours ont été enregistrées à l’échelle des événements de pêche en présence d’observateurs à bord ou par surveillance vidéo. Bien que les données sur les prises et l’effort de pêche soient disponibles pour les années antérieures pour la plupart des espèces, elles ne sont pas de la même qualité et, pour la plupart des années, ne contiennent pas d’information sur la latitude ou l’ID du navire. Ces données antérieures sont probablement utiles pour l’évaluation de certaines espèces, mais, dans le présent rapport, nous avons limité la présentation des données sur les CPUE dans les pêches commerciales aux données de grande qualité de 1996 et après pour éviter les nombreuses mises en garde à prendre en compte espèce par espèce et décennie par décennie. Nous pensons qu’il est préférable de laisser la présentation des CPUE historiques aux évaluations de stocks propres à des espèces qui peuvent prendre en compte ces données de manière plus rigoureuse.

Étant donné que nous disposons de données sur les navires individuels de cette flotte moderne, et conformément aux analyses antérieures des stocks de poisson de fond du Pacifique, nous avons défini une « flotte » qui ne comprend que les navires admissibles répondant à certains critères de capture régulière de morue du Pacifique. Nous suivons l’approche utilisée dans un certain nombre d’évaluations récentes des stocks de poisson de fond de la Colombie-Britannique en exigeant que les navires aient capturé l’espèce en au moins 100 traits pour toutes les années d’intérêt et qu’ils aient dépassé un seuil de cinq sorties (sorties ayant enregistré certaines espèces) pendant au moins cinq ans, de 1996 à r as.numeric(format(Sys.Date(), "%Y")) - 1.

DÉFINITION DES PRÉDICTEURS DU MODÈLE DE NORMALISATION

Pour la profondeur et la latitude, nous avons regroupé les valeurs dans une séquence de tranches pour permettre des relations non linéaires entre ces prédicteurs et les CPUE [e.g., @maunder2004a]. Pour la profondeur, nous avons divisé la profondeur du chalut en tranches de 25 m de large Pour la latitude, nous avons utilisé des tranches de 0,1 degré de large. Afin d’obtenir suffisamment de données pour estimer un coefficient pour chaque niveau de facteur, nous avons limité la fourchette des tranches de profondeur à ceux qui se situaient dans la probabilité cumulative d’observations positives de 0,1 % à 99,9 %, puis nous avons supprimé tout niveau de facteur (pour tous les prédicteurs) contenant moins de 0,1% de ces observations positives.

Les prédicteurs qui sont traités comme des facteurs dans un modèle statistique ont besoin d’un niveau de référence ou de base, un niveau à partir duquel les autres coefficients pour cette variable estiment une différence. Le niveau de base devient alors la valeur prédictive utilisée dans la prédiction de l’indice normalisé. Nous avons choisi le niveau de facteur le plus fréquent comme niveau de base; un choix commun pour ces types de modèles [@maunder2004a]. Par exemple, nous avons établi le mois de base comme étant le mois le plus commun observé dans l’ensemble de données filtré pour seulement les traits où l’espèce a été capturée. Ce choix du niveau de base n’affecte que l’interception ou la magnitude relative de notre indice en raison de la forme de notre modèle (voir ci-dessous) et n’a aucune incidence fonctionnelle sur les séries chronologiques des CPUE dans les pêches commerciales présentées dans ce rapport puisqu’elles sont toutes mises à l’échelle à ce même maximum et affichées sans unités.

UN MODÈLE DE NORMALISATION DE L’INDICE DU MLGM TWEEDIE

Les données sur les CPUE des pêches contiennent à la fois des valeurs nulles et des valeurs continues positives. Diverses approches ont été utilisées dans la documentation sur les pêches pour modéliser ces données. Une approche a consisté à ajuster un modèle linéaire généralisé (MLG) Delta, modèle qui ajuste les valeurs nulles par rapport aux valeurs autres selon une régression logistique (un MLG binomial et un lien logit) et les valeurs positives selon une régression linéaire adaptée aux données transformées en logarithme ou un MLG Gamma avec un lien logarithmique [e.g., @maunder2004a; @thorson2013]. La probabilité que les CPUE ne soient pas nulles pour la première composante peut alors être multipliée par les CPUE attendues pour la deuxième composante en vue d’obtenir une estimation inconditionnelle des CPUE. Toutefois, cette approche pose problème sur quelques points :

L’approche du MLG Delta ajoute de la complexité en nécessitant l’ajustement et la production de rapports sur deux modèles.
Dans l’approche typique du MLG Delta, les deux modèles sont ajustés selon des liens distincts et les coefficients ne peuvent donc pas être combinés.
L’approche du MLG Delta suppose l’indépendance entre les deux composantes [e.g., @thorson2017].
Ce qui est peut-être le plus important pour nous, c’est qu’un MLG Delta dans lequel les deux modèles utilisent des liens différents donne un indice final pour lequel la tendance dépend des niveaux de référence spécifiques auxquels les prédicteurs sont fixés [e.g., @maunder2004a].

La distribution de Tweedie [@jorgensen1987] résout les problèmes ci-dessus [e.g., @candy2004; @shono2008; @foster2013; @lecomte2013; @thorson2017], mais elle n’a pas été largement utilisée, probablement principalement en raison des frais liés au calcul de la fonction de densité de probabilité de Tweedie. Récemment, la fonction de densité de Tweedie a été introduite dans le logiciel TMB [@kristensen2016] et peut être adaptée relativement rapidement aux grands ensembles de données et aux modèles comportant de nombreux paramètres d’effets fixes et aléatoires, soit au moyen de modèles TMB codés sur mesure, soit au moyen du prologiciel du MLMG TMB R [@brooks2017].

En plus d’un paramètre moyen, la distribution de Tweedie a deux autres paramètres : un paramètre de puissance $p$ et un paramètre de dispersion $\phi$. Si $1 < p < 2$, alors la distribution de Tweedie représente une distribution composée entre la distribution de Poisson ($p = 1$) et la distribution Gamma ($p = 2$) (figure\ \@ref(fig:cpue-tweedie-ex)). En fait, le Tweedie est aussi appelé la distribution composée Poisson-Gamma dans ce cas délimité. Nous notons toutefois que la distribution composée Poisson-Gamma est souvent utilisée pour faire référence à un reparamétrage dans laquelle les composantes Poisson et Gamma sont ajustées de sorte qu’elles ne sont pas censées avoir les mêmes coefficients prédictifs que dans la distribution de Tweedie [e.g., @foster2013; @lecomte2013].

Nous adaptons le MLMG de Tweedie (modèle linéaire généralisé à effets mixtes) comme suit :

\begin{align} (#eq:cpue-tweedie) y_i &\sim \mathrm{Tweedie}(\mu_i, p, \phi), \quad 1 < p < 2,\ \mu_i &= \exp \left( \bm{X}i \bm{\beta} + \alpha^\mathrm{locality}{j[i]} + \alpha^\mathrm{locality-year}{k[i]} + \alpha^\mathrm{vessel}{l[i]} \right),\ \alpha^\mathrm{locality}j &\sim \mathrm{Normal}(0, \sigma^2{\alpha \; \mathrm{locality}}),\ (#eq:cpue-locality-year) \alpha^\mathrm{locality-year}k &\sim \mathrm{Normal}(0, \sigma^2{\alpha \; \mathrm{locality-year}}),\ (#eq:cpue-vessel) \alpha^\mathrm{vessel}l &\sim \mathrm{Normal}(0, \sigma^2{\alpha \; \mathrm{vessel}}), \end{align}

où $i$ représente un trait unique, $y_i$ représente les CPUE (kg par heures de chalutage), $\bm{X_i}$ représente un vecteur de prédicteurs à effet fixe (année, profondeur, mois, latitude), $\bm{\beta}$ représente un vecteur de coefficients et $\mu_i$ représente les CPUE attendues pour un trait. Les interceptions des effets aléatoires (symboles $\alpha$) peuvent différer de l’interception globale par emplacement $j$ ((\alpha^\mathrm{locality}_j), Figure \@ref(fig:cpue-locality-map)), emplacement-année $k$ ((\alpha^\mathrm{locality-year}_k)), et par navire $l$ ((\alpha^\mathrm{vessel}_l)), puis sont limitées par des distributions normales dont les écarts-types respectifs sont marqués par des paramètres $\sigma$. En incluant les interactions emplacement-année, nous permettons aux différents emplacements d’avoir des tendances d’indice des CPUE uniques (quelque peu limitées par la distribution des effets aléatoires) tout en estimant la tendance moyenne globale (figure \@ref(fig:cpue-spaghetti)).

On peut alors calculer l’estimation normalisée des CPUE pour l’année $t$, $\mu_t$, comme suit :

\begin{equation} \mu_t = \exp \left(\bm{X}_t \bm{\beta} \right) \end{equation}

où $\bm{X_t}$ représente un vecteur de prédicteurs fixé aux niveaux de référence ($r$) avec l’année fixée à l’année d’intérêt. Comme chacune des interceptions aléatoires $\alpha$ est fixée à zéro, l’indice est prédit pour un emplacement moyen, un emplacement-année et un navire. Nous avons estimé les effets fixes avec le maximum de vraisemblance marginale tout en intégrant les effets aléatoires avec le logiciel statistique TMB au moyen du prologiciel du MLMG TMB R. Nous avons utilisé les erreurs types ($\mathrm{SE}$) telles qu’elles sont calculées par TMB avec le $\log (\mu_t)$ à l’aide de la méthode Delta généralisée. Nous avons ensuite calculé les intervalles de confiance de Wald de 95\% comme suit : $\exp (\mu_t \pm 1.96 \mathrm{SE}_t)$. Aux fins de comparaison, nous avons calculé une série chronologique non normalisée en additionnant les prises de chaque année et en les divisant par la somme de l’effort annuel (les lignes pointillées sur les pages des figures).

library(dplyr)
d_cpue <- readRDS(here::here("report/data-cache/cpue-index-dat.rds"))
d_cpue$locality_code <- paste(d_cpue$major_stat_area_code,
  d_cpue$minor_stat_area_code, d_cpue$locality_code, sep = "-")
locality_count <- d_cpue %>%
  filter(latitude >= 45) %>%
  mutate(year = lubridate::year(best_date)) %>%
  filter(year >= 1996, year <= 2017) %>%
  group_by(locality_code) %>%
  summarize(n = n()) %>%
  arrange(-n)

(ref:cpue-locality-map-cap) Les 100 principaux emplacements du MPO (selon le nombre d’événements de pêche) utilisés dans les modèles commerciaux de normalisation des CPUE. Au total, il y a r nrow(locality_count) emplacements possibles enregistrés dans l’ensemble de données.

gfplot:::plot_dfo_localities(top_n(locality_count, 75)$locality_code) +
  xlab(en2fr("Easting", french)) +
  ylab(en2fr("Northing", french))

plot_tweedie_ex <- function(df, max = 15) {
  xx <- seq(0, max, length.out = 1000)
  out <- plyr::mdply(df, function(mu, power, phi) {
    data.frame(x = xx,
      y = tweedie::dtweedie(y = xx, mu = mu, power = power, phi = phi))
  }) %>%
    mutate(phi = paste("phi =", phi)) %>%
    mutate(power = paste("p =", power)) %>%
    mutate(mu = paste("μ =", mu))

  ggplot(out, aes(x, y, colour = mu)) +
    geom_line() +
    facet_grid(power~phi, scales = "free") +
    gfplot::theme_pbs() +
    labs(colour = "μ") +
    xlab(en2fr("Value", french)) + ylab(en2fr("Density", french)) +
    coord_cartesian(expand = FALSE, xlim = c(-0.2, max(out$x))) +
    scale_colour_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2")
}

df <- expand.grid(power = c(1.2, 1.6, 1.8), mu = c(1, 3, 6), phi = c(0.5, 1, 2))
plot_tweedie_ex(df)

cpue_model <- readRDS(here::here("report/cpue-cache/petrale-sole-5AB-model.rds"))

# d_cpue <- readRDS(here::here("report/data-cache/cpue-index-dat.rds"))
d_cpue_pred <- gfplot::predict_cpue_index_tweedie(cpue_model, center = FALSE)
d_cpue_pred$formula_version <- "Full standardization"
if (file.exists(here::here("report/cpue-cache/petrale-sole-5AB-fleet.rds")))
  d_fleet <- readRDS(here::here("report/cpue-cache/petrale-sole-5AB-fleet.rds"))

if (exists("d_fleet"))
  plot_cpue_spaghetti(
    model = cpue_model,
    fleet = d_fleet,
    index_data = d_cpue_pred,
    era = "modern")

(ref:caption-cpue-quantile-residuals) Histogrammes des résidus de quantiles aléatoires [@dunn1996] pour les MLMG de normalisation des CPUE. Les histogrammes illustrent la distribution de 10 000 résidus de quantiles aléatoire sélectionnés au hasard. Les lignes en pointillés indiquent la densité de probabilité pour une distribution normale avec le même écart-type.

qr <- gfplot::qres_tweedie(cpue_model, response = "cpue")
gfplot::plot_qres_histogram(qr)

\clearpage

files <- list.files(here::here('report/cpue-cache'), pattern = "-model")
library("doParallel")
registerDoParallel(cores = parallel::detectCores())
cpue_models_summary_file <- here::here("report/cpue-cache/model-summaries.rds")
if (!file.exists(cpue_models_summary_file)) {
  cpue_models_summary <- plyr::ldply(files, function(f) {
    # cat(f, "\n")
    this_model <- readRDS(here::here('report/cpue-cache', f))
    if (length(this_model) > 1L) {
      broom_model <-
        tibble(species = gsub("([a-z-]+)-[A-Z0-9]+-model.rds", "\\1", f)) %>%
        mutate(species = gsub('-', ' ', species)) %>%
        mutate(area = gsub("([a-z-]+)-([A-Z0-9]+)-model.rds", "\\2", f)) %>%
        mutate(sigma = sigma(this_model)) %>%
        mutate(thetaf = plogis(exp(this_model$fit$par[["thetaf"]])) + 1) %>%
        mutate(vessel = glmmTMB::VarCorr(this_model)$cond$vessel[1,1]) %>%
        mutate(locality = glmmTMB::VarCorr(this_model)$cond$locality[1,1]) %>%
        mutate(year_locality = glmmTMB::VarCorr(this_model)$cond$year_locality[1,1])
    }
  }, .parallel = TRUE) %>%
    arrange(species, area)
  saveRDS(cpue_models_summary, file = cpue_models_summary_file)
} else {
  cpue_models_summary <- readRDS(cpue_models_summary_file)
}

for_table <- mutate(
  cpue_models_summary, species = gfsynopsis:::first_cap(species)) %>%
  as_tibble() %>%
  mutate(species = gsub("Rougheyeblackspotted Rockfish Complex", "Rougheye/Blackspotted", species)) %>%
  mutate(species = en2fr(species, french))

for_table$species <- purrr::map_chr(for_table$species, gfsynopsis:::first_cap)

for_table %>% 
  csasdown::csas_table(digits = 1,
    col.names = c(
      "Espèces", "Zone", "$\\phi$", "$p$",
      "$\\tau_\\mathrm{vessel}$", "$\\tau_\\mathrm{locality}$", "$\\tau_\\mathrm{year-locality}$"),
    caption = "Estimations de paramètres à partir de GLMMs de normalisation de la CPUE")