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

TENDANCES DE L’INDICE DE BIOMASSE DÉRIVÉ DES RELEVÉS {#app:survey-trend-models}

Tous les ans ou deux ans, Pêches et Océans Canada (MPO) effectue un certain nombre de relevés de recherche plurispécifiques indépendants de la pêche. Quatre relevés synoptiques au chalut de fond et deux relevés à la ligne et à l’hameçon comptent parmi ceux-ci (figure \@ref(fig:intro-maps)). Les zones de relevé correspondent sensiblement aux zones de gestion du poisson de fond (figure \@ref(fig:management-map)) pour lesquelles des quotas de pêche sont accordés dans les eaux de la Colombie-Britannique. En effet, pour les relevés synoptiques au chalut, le relevé sur la côte ouest de l’île de Vancouver correspond sensiblement aux zones de gestion 3CD, le relevé dans le bassin de la Reine-Charlotte correspond sensiblement aux zones 5AB, le relevé dans le détroit d’Hécate correspond sensiblement aux zones 5CD et les relevés sur la côte ouest de Haida Gwaii correspondent sensiblement à la zone 5E. Le relevé à la palangre sur fond dur est divisé en segments nord et sud. La zone de relevé sud correspond sensiblement aux zones de gestion 3CD et 5AB, tandis que la zone de relevé nord correspond sensiblement aux zones\ 5CDE.

RELEVÉ SYNOPTIQUE AU CHALUT DE FOND

Le MPO et la Canadian Groundfish Research and Conservation Society ont mis en œuvre un ensemble coordonné de relevés au chalut de fond qui couvrent le plateau continental et le talus supérieur de la majeure partie de la côte de la Colombie-Britannique. Les relevés suivent un plan aléatoire stratifié en fonction de la profondeur et utilisent les mêmes engins de pêche au chalut de fond et les mêmes protocoles de pêche [@sinclair2003syn]. Les relevés ont été conçus pour fournir un synopsis de toutes les espèces pouvant être capturées au chalut de fond plutôt que de se concentrer sur des espèces précises. Il y a des relevés synoptiques (SYN) à quatre endroits en tout : détroit d’Hécate (DH), côte ouest de l’île de Vancouver (COIV), bassin de la Reine-Charlotte (BRC) et côte ouest de Haïda Gwaii (COHG) (figure \@ref(fig:intro-maps)). Les relevés dans le bassin de la Reine-Charlotte et sur la côte ouest de Haida Gwaii ont été effectués à bord de bateaux de pêche commerciale affrétés, tandis que les relevés dans le détroit d’Hécate et sur la côte ouest de l’île de Vancouver ont été effectués sur le chalutier de recherche WE Ricker de la Garde côtière canadienne ou sur des bateaux de pêche commerciale affrétés lorsque le WE Ricker n’était pas disponible. Deux des relevés synoptiques sont effectués chaque année en alternance, de sorte que chaque levé est effectué tous les deux ans.

RELEVÉS À LA PALANGRE SUR FOND DUR

La Pacific Halibut Management Association (PHMA), en consultation avec le MPO, a entrepris en 2006 un relevé de recherche à la palangre à stratification aléatoire en fonction de la profondeur, effectué par des bateaux de pêche commerciale affrétés. C’est ce qu’on appelle les relevés à la palangre sur fond dur dans les eaux extérieures. Le relevé utilise des engins et des méthodes de pêche à la palangre normalisés et alterne annuellement entre le nord et le sud de la Colombie-Britannique. Il est conçu pour fournir les taux de prise de toutes les espèces et des échantillons biologiques de sébaste provenant des eaux côtières extérieures de la Colombie-Britannique pour l’évaluation des stocks.

Les relevés à la palangre sur fond dur dans les eaux intérieures sont effectués dans la zone de gestion 4B. Ces relevés sont conçus pour fournir des indices d’abondance indépendants de la pêche ainsi que des échantillons biologiques afin d’améliorer l’évaluation du sébaste aux yeux jaunes (Sebastes ruberrimus) et du sébaste à dos épineux (S. maliger) dans la zone de gestion 4B. Ils ont commencé dans le détroit de Johnstone et le passage Discovery, dans les secteurs 12 et 13 de gestion des pêches du Pacifique en 2003 et 2004, et ils alternent maintenant d’une année sur l’autre pour couvrir les parties nord (secteurs 12 et 13) et sud (secteurs 14-20, 28 et 29) des eaux intérieures. Ces relevés font également appel à des engins et des méthodes de pêche à la palangre normalisés.

RELEVÉ INDÉPENDANT DE LA PÊCHE DE LA COMMISSION INTERNATIONALE DU FLÉTAN DU PACIFIQUE

Le relevé à la ligne fixe indépendant de la pêche de la Commission internationale du flétan du Pacifique (CIFP) est la plus longue série chronologique de données de relevés à la palangre en Colombie-Britannique. Il fournit des données sur la répartition, la biomasse, l’âge, la croissance et la maturité qui sont utilisées dans l’évaluation annuelle du flétan du Pacifique (Hippoglossus stenolepis). À l’annexe \@ref(app:iphc-survey-index), nous décrivons comment nous utilisons les données du relevé pour établir un indice d’abondance cohérent pour d’autres espèces sur une série chronologique aussi longue que possible, malgré l’évolution de la conception du relevé au fil des ans.

AUTRES RENSEIGNEMENTS SUR LES RELEVÉS

Les documents suivants fournissent des précisions sur la conception des divers relevés mentionnés dans le présent rapport :

Relevé synoptique, bassin de la Reine-Charlotte (SYN BRC) : @williams2018synqcs
Relevé synoptique, côte ouest de l’île de Vancouver (SYN COIV) : @nottingham2017synwcvi
Relevé synoptique, détroit d’Hécate (SYN DH) : @wyeth2018synhs
Relevé synoptique, côte ouest de Haida Gwaii (SYN COHG) : @williams2018synwchg
Relevé à la palangre sur fond dur dans les eaux extérieures (RPFD EXT) : @doherty2019hbllout
Relevé à la palangre sur fond dur dans les eaux intérieures (RPFD INT) : @lochead2007irf
Relevé de la communauté d’espèces dans le détroit d’Hécate (RCE DH) : @choromanski2004hsmsa
Relevé indépendant de la pêche de la Commission internationale du flétan du Pacifique (RIP CIFP) : @flemming2012iphc

if (!file.exists(here::here("report", "other-surveys.csv"))) {
  other_survey_data <- gfplot::get_other_surveys()
  readr::write_csv(other_survey_data, here::here("report", "other-surveys.csv"))
} 
f <- here::here("report", "other-surveys.csv")
other_surveys <- read.csv(f, stringsAsFactors = FALSE, strip.white = TRUE)
names(other_surveys) <- gfplot:::firstup(tolower(gsub("_", " ", names(other_surveys))))
other_surveys$X <- NULL # just in case

if (french) {
  names(other_surveys)[1] <- en2fr("Survey", TRUE)
  names(other_surveys)[2] <- "Enquêtes menées depuis 2008"
}

\clearpage

csasdown::csas_table(other_surveys, caption = "Autres relevés effectués par le MPO dans la région du Pacifique qui peuvent s’appliquer à des analyses portant sur des espèces particulières. Dans le présent rapport, ces relevés ne sont présentés que dans les encadrés de disponibilité des données de relevés dans le coin inférieur droit de chaque série de pages de figures", format = "latex") %>% 
  kableExtra::column_spec(1, width = "4in") %>%
  kableExtra::column_spec(2, width = "2.5in")

Les principales tendances de l’indice de biomasse illustrées dans les pages de figures représentent les tendances de l’indice fondé sur le plan qui ont toujours été utilisées pour les évaluations des stocks de poisson de fond de la Station biologique du Pacifique. Nous avons extrait de GFBio les tendances relatives de l’indice de biomasse dérivé des relevés au chalut et à la palangre, générées selon la même approche utilisée pour tous les récents rapports d’évaluation des stocks de poisson de fond en Colombie-Britannique. Le code permettant d’effectuer les calculs, rédigé à l’origine par Norm Olsen de la Station biologique du Pacifique, est automatiquement appliqué aux données de relevés disponibles pour générer les indices dans la base de données GFBio. Pour plus de clarté, vous trouverez les équations pertinentes ci-dessous. Nous comparons également les estimations fondées sur des modèles géostatistiques des tendances de l’indice de biomasse dérivé des relevés au chalut (section \@ref(sec:geostatistical)).

ESTIMATIONS FONDÉES SUR LE PLAN DE RELEVÉ AU CHALUT {#sec:trawl-methods}

Pour tous les relevés au chalut et pour une espèce donnée, nous avons calculé la densité relative de biomasse $B$ au cours de l’année $t$ comme suit :

\begin{equation} B_t = \sum_{i = 1}^k C_{t,i} A_i \end{equation}

où $C_{t,i}$ représente la capture par unité d’effort (CPUE) moyenne en kg/km^2^ pour l’espèce au cours de l’année $t$ et à la strate $i$, $A_i$ représente la superficie de la strate $i$ en km^2^, puis $k$ représente le nombre de strates. Nous avons calculé la CPUE ($C_{t,i}$) pour une espèce donnée au cours de l’année $t$ et à la strate $i$ comme suit :

\begin{equation} C_{t,i} = \frac{\sum_{j = 1}^{n_{t,i}} \left( W_{t,j} / D_{t,j} w_{t,j}\right)}{n_{t,i}} \end{equation}

où $W_{t,j}$ représente le poids des captures (kg) de l’espèce au cours de l’année $t$, à la strate $i$, et au trait $j$; $D_{t,j}$ représente la distance parcourue en km par le trait $j$ au cours de l’année $y$; $w_{t,j}$ représente la largeur nette d’ouverture en km au cours de l’année $y$ et au trait $j$; $n_{t,i}$ représente le nombre de traits au cours de l’année $t$ et à la strate $i$.

ESTIMATIONS FONDÉES SUR LE PLAN DE RELEVÉ À LA PALANGRE {#sec:ll-methods}

Pour tous les relevés à la palangre sur fond dur et pour une espèce donnée, nous avons calculé la densité relative de biomasse $B$ au cours de l’année $t$ comme suit :

\begin{equation} B_t = \sum_{i = 1}^k C_{t,i} A_i \end{equation}

où $C_{t,i}$ représente la CPUE moyenne en pièce (poisson) par km\textsuperscript{2} pour l’espèce au cours de l’année $t$ et à la strate $i$, $A_i$ représente la superficie de la strate $i$ en km^2^, puis $k$ représente le nombre de strates. Nous avons calculé la CPUE ($C_{t,i}$) pour une espèce donnée au cours de l’année $t$ et à la strate $i$ comme suit :

\begin{equation} C_{t,i} = \frac{\sum_{j = 1}^{n_{t,i}} \left( N_{t,j} / H_{t,j} w_{t,j}\right)}{n_{t,i}} \end{equation}

où $N_{t,j}$ représente le nombre de poissons capturés pour l’espèce au cours de l’année $t$, à la strate $i$, et dans l’ensemble $j$; $H_{t,j}$ représente le nombre d’hameçons $\times$, l’espacement des hameçons en km dans l’ensemble $j$ à l’année $t$; $w_{t,j}$ représente une largeur de balayage arbitraire de 30 pieds ou 0,009144 km pour l’année $t$ et le trait $j$; $n_{t,i}$ représente le nombre d’ensembles à l’année $t$ et à la strate $i$. L’espacement des hameçons est de 8 pieds ou 0,0024384 km pour les relevés à la palangre sur fond dur dans les eaux intérieures et extérieures

Des précisions sur l’estimation de l’indice de densité de biomasse fondé sur le plan du relevé de la CIFP sont présentées à l’annexe \@ref(app:iphc-survey-index).

INTERVALLES DE CONFIANCE DE L’INDICE FONDÉ SUR LE PLAN

Nous avons calculé les intervalles de confiance d’auto-amorçage concernant $B_t$ en calculant de façon répétée $B_t$ à partir des équations ci-dessus, mais en ré-échantillonnant chaque fois, avec remplacement, à partir des traits disponibles dans chaque strate. Nous avons tiré 1 000 répliques par auto-amorçage de $B_t$, $B_t^{\mathrm{rep}}$, et calculé des intervalles de confiance corrigés et ajustés (BCa) à 95 % pour des intervalles de confiance à biais corrigé et ajusté (BCa) [@efron1987] concernant $B_t^{\mathrm{rep}}$.

TENDANCES GÉOSTATISTIQUES ET SPATIOTEMPORELLES DE L’INDICE DE BIOMASSE {#sec:geostatistical}

Les estimations fondées sur le plan décrites ci-dessus supposent que la densité moyenne de poissons est la même dans toutes les strates du relevé et que la seule source de variance entre les échantillons est la stochasticité d’échantillonnage elle-même [@petitgas2001]. Cependant, nous savons que ce n’est pas vrai. Une part importante de la variation de la densité de poissons dans une strate peut être attribuée au fait que l’habitat convient mieux ou moins bien à un poisson donné, et cette convenance peut découler de nombreux facteurs au-delà de la simple profondeur à laquelle les strates sont stratifiées [@shelton2014]. Nous savons aussi que les poissons ne perçoivent pas leur habitat en fonction de ces limites de strates exactes et que les processus écologiques ont tendance à être en corrélation spatiale avec les processus plus semblables plus rapprochés les uns des autres que ceux plus éloignés. Les estimations fondées sur le plan n’exploitent pas cette possible corrélation spatiale.

La modélisation géostatistique des données de relevés vise à résoudre ces problèmes en modélisant la densité des poissons sous la forme d’une surface spatiale lisse, possiblement le résultat de variables explicites relatives à l’habitat comme la profondeur, mais aussi le produit d’autres effets spatiaux non observés ou « latents ». Ces dernières années, il s’est opéré une transition vers une normalisation « fondée sur des modèles » des tendances de l’indice de biomasse dérivé de relevés [e.g., @shelton2014; @thorson2015; @webster2017]. Nous incluons comme point de comparaison les tendances de l’indice fondé sur des modèles dérivé des relevés synoptiques au chalut afin de permettre au lecteur de déterminer quand les deux approches peuvent varier. De grandes différences sont probablement le résultat du positionnement aléatoire des ensembles de relevés pour une année donnée. Les ensembles se retrouvent dans un habitat qui convient particulièrement bien ou non à une espèce donnée. Les auteurs des évaluations des stocks de poisson de fond de la Colombie-Britannique pourraient envisager une normalisation de l’indice fondé sur des modèles au cas par cas.

Nous utilisons le modèle géostatistique décrit à l’annexe \@ref(app:spatial-modeling), en y ajoutant les effets spatiotemporels aléatoires ($\epsilon_{s,t}$; définis pour les emplacements dans l’espace $s$ et au temps $t$) et les prédicteurs annuels de biomasse moyenne chaque année :

\begin{align} y_{s,t} &\sim \mathrm{Tweedie}(\mu_{s,t}, p, \phi), \quad 1 < p < 2,\ \mu_{s,t} &= \exp \left( \bm{X}{s,t} \bm{\beta} + \omega_s + \epsilon{s,t} \right). \end{align}

Comme pour le modèle spatial, nous supposons que les effets aléatoires spatiaux ($\omega_s$) sont tirés d’une distribution normale multidimensionnelle avec une certaine matrice de covariance $\bm{\Sigma}_\omega$ :

\begin{equation} \bm{\omega} \sim \mathrm{MVNormal} \left( \bm{0}, \bm{\Sigma}_\omega \right), \end{equation}

Nous supposons la même chose pour les effets aléatoires spatiotemporels, chaque tranche de temps comportant son propre ensemble indépendant d’effets aléatoires ($\bm{\epsilon}t$) avec matrice de covariance $\bm{\Sigma}{\epsilon,t}$ :

\begin{equation} \bm{\epsilon}t \sim \mathrm{MVNormal} \left( \bm{0}, \bm{\Sigma}{\epsilon,t} \right). \end{equation}

Les effets aléatoires spatiaux tiennent compte de facteurs spatiaux qui sont constants au fil du temps, par exemple la profondeur et le type de substrat. Les effets spatiotemporels aléatoires tiennent compte de facteurs qui varient d’une année à l’autre dans l’espace, comme la température du fond, les régimes de circulation de l’eau, les interactions entre espèces et les déplacements des espèces.

Nous utilisons ici la morue du Pacifique comme exemple pour illustrer les composantes du modèle. L’approche comprend la même génération de « nœuds » spatiaux que dans le modèle spatial décrit à l’annexe \@ref(app:spatial-modeling) (figure \@ref(fig:sdmTMB-spt-spde)). Nous utilisons 200 nœuds, un nombre qui, pour la couverture spatiale des relevés synoptiques au chalut, semble être d’une résolution suffisamment élevée pour capturer la variation spatiale et spatiotemporelle. Nous avons vérifié cette hypothèse en augmentant le nombre de nœuds pour certaines espèces en nous assurant que les tendances estimées ne variaient pas sur le plan qualitatif.

Nous pouvons projeter les prédictions du modèle à un quadrillage à petite échelle (2 km $\times$ 2 km) en utilisant la matrice de projection de covariance (figures \@ref(fig:sdmTMB-spt-plot-all-effects-sqrt), \@ref(fig:sdmTMB-spt-plot-all-effects-sqrt-depth)). Nous pouvons également examiner chaque composante du modèle. Pour les modèles sans prédicteurs d’effets fixes pour la profondeur, les prédictions des effets fixes chaque année sont constantes dans l’espace (figure \@ref(fig:sdmTMB-spt-plot-fixed-effects)). Les effets aléatoires spatiaux sont constants d’une année à l’autre \@ref(fig:sdmTMB-spt-plot-spatial-effects)) contrairement aux effets aléatoires spatiotemporels (figure \@ref(fig:sdmTMB-spt-plot-spatiotemporal-effects)). Nous pouvons examiner les résidus dans l’espace et dans le temps pour vérifier s’il semble subsister une autocorrélation spatiale ou spatiotemporelle (figure \@ref(fig:sdmTMB-spt-spatial-residuals)).

Nous pouvons alors calculer la biomasse attendue $B_t$ pour l’année $t$, comme suit :

\begin{equation} B_t = \sum_{j = 1}^{n_j} w_j \cdot \exp \left( \bm{X}{j,t} \bm{\beta} + \omega_j + \epsilon{j,t} \right), \end{equation}

où $j$ désigne une cellule de quadrillage de 2 km $\times$ 2 km dans la zone de relevé et $w_j$ représente le poids ou la superficie de cette cellule de quadrillage (4 km^2^) (figure \@ref(fig:sdmTMB-plot-index)). En d’autres termes, nous additionnons la biomasse prévue pour chaque année dans toutes les cellules de quadrillage faisant partie de la zone de relevé. Nous avons généré les écarts-types pour les estimations annuelles de biomasse du logarithme à l’aide de la méthode Delta généralisée mise en œuvre dans Template Model Builder (TMB) [@kristensen2016]. Un peu comme dans l’étude @thorson2015, nous avons constaté que les tendances de l’indice de biomasse fondé sur des modèles avaient souvent des coefficients de variation plus faibles que les tendances de l’indice fondé sur le plan et qu’ils contribuaient souvent à stabiliser les estimations de biomasse pour les années suivantes comparativement aux tendances de l’indice fondé sur le plan (p. ex. figure \@ref(fig:sdmTMB-plot-index)).

Nous avons constaté peu de variation de l’indice de biomasse prévu entre les modèles qui incluaient ou non des covariables relatives à la profondeur (p. ex. figure \@ref(fig:sdmTMB-plot-index)). Nous avons constaté que la principale différence entre les modèles avec ou sans covariables relatives à la profondeur se trouvait dans les modèles dont ces covariables présentaient des estimations de biomasse avec une résolution légèrement plus grande dans l’espace (figure \@ref(fig:sdmTMB-spt-plot-all-effects-sqrt) par rapport à (figure \@ref(fig:sdmTMB-spt-plot-all-effects-sqrt-depth)). Nous avons choisi de ne pas inclure la profondeur et la profondeur au carré comme prédicteurs dans l’indice principal de biomasse illustré dans le présent rapport. Même si presque tous les cas qui comprenaient des covariables relatives à la profondeur ont généré un indice semblable aux cas les excluant, dans quelques cas, leur inclusion a généré ce qui semblait être des écarts irréalistes dans la biomasse lorsque la forme quadratique de la relation entre la profondeur et la biomasse générait des estimations excessivement élevées ou faibles de biomasse à la limite des polygones de relevé. Ce point demeure un sujet de recherche sur lequel les auteurs enquêteront plus tard.

library(dplyr)
library(ggplot2)
library(sdmTMB)

survey <- "SYN QCS"
n_knots <- 200
bias_correct <- FALSE
# dat <- readRDS(here::here("report/data-cache/walleye-pollock.rds"))
dat <- readRDS(here::here("report/data-cache/pacific-cod.rds"))
dat <- gfplot:::tidy_survey_sets(dat$survey_sets, survey, years = seq(1, 1e6),
  density_column = "density_kgpm2")
dat <- mutate(dat, density = density * 1000 * 1000)
dat <- gfplot:::interp_survey_bathymetry(dat)$data
dat <- gfplot:::scale_survey_predictors(dat)

# grid_locs <- gfplot:::make_prediction_grid(
#   filter(dat, year == max(dat$year)), survey = survey,
#   cell_width = 2)$grid
# grid_locs <- rename(grid_locs, depth = akima_depth)
# grid_locs$year <- NULL

grid_locs <- gfplot::synoptic_grid %>%
  dplyr::filter(survey == "SYN QCS") %>%
  dplyr::select(X, Y, depth)

grid_locs$depth_scaled <-
  (log(grid_locs$depth) - dat$depth_mean[1]) / dat$depth_sd[1]
grid_locs$depth_scaled2 <- grid_locs$depth_scaled^2

# Expand the prediction grid to create a slice for each time:
original_time <- sort(unique(dat$year))
nd <- do.call("rbind",
  replicate(length(original_time), grid_locs, simplify = FALSE))
nd[["year"]] <- rep(original_time, each = nrow(grid_locs))
grid_locs <- nd

options(OutDec =  ".")
spde <- make_spde(dat$X, dat$Y, n_knots = n_knots)
options(OutDec =  ",")
plot_spde(spde)

m_st <- sdmTMB(
  formula = density ~ 0 + as.factor(year),
  data = dat, time = "year", spde = spde, family = tweedie(link = "log"),
  silent = TRUE)
predictions <- predict(m_st, newdata = grid_locs, return_tmb_object = TRUE)

m_st_depth <- sdmTMB(
  formula = density ~ 0 + as.factor(year) + depth_scaled + depth_scaled2,
  data = dat, time = "year", spde = spde, family = tweedie(link = "log"),
  silent = TRUE)
predictions_depth <- predict(m_st_depth, newdata = grid_locs, return_tmb_object = TRUE)

dat$resids <- residuals(m_st) # randomized quantile residuals
# gfplot::plot_qres_histogram(dat$resids)

dat$resids_depth <- residuals(m_st_depth) # randomized quantile residuals
# gfplot::plot_qres_histogram(dat$resids_depth)

plot_map_spt <- function(dat, column) {
  ggplot(dat, aes_string("X", "Y", fill = column)) +
    geom_raster() +
    facet_wrap(~year) +
    coord_fixed() +
    labs(x = en2fr('Easting', french), y = en2fr('Northing', french))
}

ticks <- c(100, 500, 1000, 2000)
limits <- c(0, max(exp(predictions$data$est))*1.001)
plot_map_spt(predictions$data, "exp(est)") +
  ggtitle(en2fr("Prediction (fixed effects + all random effects)", french)) +
  scale_fill_viridis_c(trans = "fourth_root_power", option = "C",
    breaks = ticks, limits = limits) +
  labs(fill = if(french){
    'Estimation de la\nbiomasse (kg/km2)'
    }
    else{
      'Biomass estimate\n(kg/km^2)'
    })

\clearpage

plot_map_spt(predictions_depth$data, "exp(est)") +
  ggtitle(en2fr("Prediction (fixed effects + all random effects)", french)) +
  scale_fill_viridis_c(trans = "fourth_root_power", option = "C",
    breaks = ticks, limits = limits) +
  labs(fill = if(french){
    'Estimation de la\nbiomasse (kg/km2)'
    }
    else{
      'Biomass estimate\n(kg/km^2)'
    })

\clearpage

plot_map_spt(predictions$data, "exp(est_non_rf)") +
  ggtitle(en2fr("Prediction (fixed effects only)", french)) +
  scale_fill_viridis_c(trans = "fourth_root_power", option = "C") +
  labs(fill = if(french){
    'Estimation moyenne\n(kg/km^2)'
    }
    else{
     'Mean estimate\nkg/km^2'
      })

\clearpage

red_blue_fill <- scale_fill_gradient2(midpoint = 0,
    high = scales::muted("red"),
    mid = "white",
    low = scales::muted("blue"))

plot_map_spt(predictions$data, "omega_s") +
  ggtitle(en2fr("Spatial random effects only", french)) +
  red_blue_fill +
  labs(fill = if(french){
    'Écart par rapport à la\nbiomasse prévue dans\nl’espace de relevé'
  }
    else{
      'Deviation from\nexpected biomass\nin log space'
    })

\clearpage

plot_map_spt(predictions$data, "epsilon_st") +
  ggtitle(en2fr("Spatiotemporal random effects only", french)) +
  red_blue_fill +
  labs(fill = if(french){
    'Écart par rapport à la\nbiomasse prévue dans\nl’espace de relevé'
  }
    else{
      'Deviation from\nexpected biomass\nin log space'
    })

\clearpage

red_blue_fill2 <- scale_colour_gradient2(midpoint = 0,
    high = scales::muted("red"),
    mid = "white",
    low = scales::muted("blue"))

ggplot(dat, aes(X, Y, col = resids)) +
  geom_point() + facet_wrap(~year) + coord_fixed() +
  labs(x = en2fr('Easting', french), y = en2fr('Northing', french), colour = en2fr("Residual", french)) +
  red_blue_fill2

\clearpage

ind <- get_index(predictions, bias_correct = FALSE)
ind_depth <- get_index(predictions_depth, bias_correct = FALSE)

\clearpage

# ind_design <- readRDS(here::here('report/data-cache/walleye-pollock.rds'))
ind_design <- readRDS(here::here('report/data-cache/pacific-cod.rds'))
ind_design <- ind_design$survey_index
ind_design <- ind_design %>%
  filter(survey_abbrev  == 'SYN QCS') %>%
  mutate(
    est = biomass / exp(mean(log(biomass), na.rm = TRUE)),
    lwr = lowerci / exp(mean(log(biomass), na.rm = TRUE)),
    upr = upperci / exp(mean(log(biomass), na.rm = TRUE)),
    type = en2fr("Design based", french)
  )

ind_geo <- ind %>%
  mutate(
    lwr = lwr / exp(mean(log(est), na.rm = TRUE)),
    upr = upr / exp(mean(log(est), na.rm = TRUE)),
    est = est / exp(mean(log(est), na.rm = TRUE)),
    type = en2fr("Geostatistical", french)
  )

ind_geo_depth <- ind_depth %>%
  mutate(
    lwr = lwr / exp(mean(log(est), na.rm = TRUE)),
    upr = upr / exp(mean(log(est), na.rm = TRUE)),
    est = est / exp(mean(log(est), na.rm = TRUE)),
    type = paste0(en2fr("Geostatistical", french), " (", tolower(en2fr("Depth", french)), ")")
  )

ind_combined <- bind_rows(ind_design, ind_geo) %>%
  bind_rows(ind_geo_depth)
ggplot(ind_combined, aes(year, est, fill = type)) +
  geom_line(aes(colour = type)) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.4) +
  xlab(en2fr('Year', french)) +
  ylab(en2fr('Biomass estimate divided by geometric mean', french)) +
  labs(fill = en2fr('Type', french), colour = en2fr('Type', french)) +
  scale_fill_brewer(palette = "Dark2") +
  scale_colour_brewer(palette = "Dark2")