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

MODÉLISATION SPATIALE DE LA BIOMASSE DÉRIVÉE DES RELEVÉS {#app:spatial-modeling}

Nous avons modélisé la densité prévue de biomasse dans l’espace pour chaque espèce à l’aide de modèles géostatistiques appliqués aux données des relevés au chalut de fond et à la palangre indépendants de la pêche (p. ex. figure \@ref(fig:survey-maps)). Notre approche de modélisation est conforme aux modèles récemment utilisés pour la modélisation spatiale et la normalisation de l’indice spatiotemporel des populations de poissons de fond [e.g. @shelton2014; @ward2015; @thorson2015; @thorson2016], mais elle n’a pas, à notre connaissance, été appliquée auparavant dans les documents de recherche du MPO. Il a été démontré, entre autres, que de tels modèles améliorent les estimations de l’abondance et de la distribution du sébaste [@shelton2014] et améliorent la précision de l’estimation des indices d’abondance relative du poisson de fond [@thorson2015]. Notre modèle spécifique est ajusté avec Template Model Builder (TMB) [@kristensen2016] in R [@r2018] grâce à Integrated Nested Laplace Approximations (INLA) [@rue2009; @lindgren2015] par l’intermédiaire du progiciel R sdmTMB, que nous avons codé à cet effet (annexe \@ref(app:reproducibility)).

De façon générale, ces modèles prédisent la biomasse relative ou le taux de capture dans l’espace sous la forme d’un processus continu avec un effet quadratique pour la profondeur du fond, des effets aléatoires spatiaux qui représentent un amalgame de processus spatiaux non explicitement inclus dans le modèle, et une composante d’erreur d’observation. Après avoir ajusté le modèle aux ensembles de relevés à partir de relevés au chalut ou à la palangre, nous projetons les prédictions du modèle sur un quadrillage de 2 km $\times$ 2 km à une projection UTM 9 pour obtenir des estimations de biomasse dans l’ensemble de la zone de relevé.

Tout comme les modèles de normalisation des captures commerciales par unité d’effort (annexe \@ref(app:cpue-models)), ces modèles peuvent être représentés comme des modèles linéaires généralisés à effets mixtes (MLGM) de Tweedie avec un lien logarithmique :

\begin{align} y_s &\sim \mathrm{Tweedie}(\mu_s, p, \phi), \quad 1 < p < 2,\ \mu_s &= \exp \left( \bm{X}_s \bm{\beta} + \omega_s \right), \end{align}

où $s$ représente un emplacement spatial, $y_s$ représente la densité de poissons observée pour un ensemble de relevés, $\mu_s$ représente la densité de poissons attendue, $p$ représente le paramètre de puissance de Tweedie et $\phi$ représente le paramètre de dispersion de Tweedie. Le symbole $\bm{X_s}$ représente un vecteur de prédicteurs (une intersection, une profondeur logarithmique et une profondeur logarithmique au carré) et $\bm{\beta}$ représente un vecteur correspondant de coefficients. Selon les suppositions, les effets aléatoires spatiaux $\omega_s$ sont issus d’une distribution normale multivariée avec une matrice de covariance $\bm{\Sigma}_\omega$ centrée sur zéro :

$$\bm{\omega} \sim \mathrm{MVNormal} \left( \bm{0}, \bm{\Sigma}_\omega \right).$$

Nous avons contraint les effets aléatoires spatiaux à suivre une fonction de covariance \mbox{Mat\'ern}, qui définit le taux selon lequel la corrélation spatiale diminue avec la distance (figure \@ref(fig:sdmTMB-matern-demo)). La fonction \mbox{Mat\'ern} décrit la covariance $\Phi \left( s_j, s_k \right)$ entre les emplacements spatiaux $s_j$ et $s_k$ comme suit :

$$\Phi\left( s_j,s_k \right) = \tau^2/\Gamma(\nu)2^{\nu - 1} (\kappa d_{jk})^\nu K_\nu \left( \kappa d_{jk} \right),$$

où $\tau^2$ représente la variance spatiale, $\Gamma$ représente la fonction Gamma, $K_\nu$ représente la fonction de Bessel, $d_{jk}$ représente la distance euclidienne entre les emplacements $s_j$ et $s_k$, puis $\kappa$ représentent un paramètre d’échelle estimé [e.g., @lindgren2011]. Le paramètre $\nu$ contrôle le lissage de la fonction de covariance. Nous avons établi $\nu = 1$, ce qui nous permet de tirer parti de l’approximation de l’équation différentielle partielle stochastique (EDPS) par rapport aux champs aléatoires gaussiens de Markov (CAGM) pour augmenter considérablement l’efficacité de calcul [@lindgren2011].

Notre modèle spatial entre dans la catégorie générale des modèles de « processus prédictifs » [e.g., @latimer2009; @shelton2014; @anderson2018], dans lesquels le modèle suit un nombre limité de « nœuds » qui estiment des variations spatiales inexpliquées (figure \@ref(fig:sdmTMB-spde)). En enregistrant un nombre de nœuds inférieur à celui de l’ensemble des données spatiales, nous pouvons accroître l’efficacité de calcul. Les prédictions du modèle peuvent être projetées aux emplacements des données d’origine ou à tout nouvel ensemble d’emplacements, pourvu qu’une matrice de covariance appropriée puisse être calculée [e.g., @latimer2009]. Un nombre plus élevé de nœuds permet une meilleure approximation des effets aléatoires spatiaux à un coût de calcul plus élevé. Nous ajustons les modèles spatiaux avec un nœud de moins que le nombre d’observations s’il y avait moins de 200 observations et 200 nœuds autrement. Suivant une pratique courante [e.g., @shelton2014; @thorson2017b; @anderson2018], nous avons choisi l’emplacement des nœuds à l’aide d’un algorithme de regroupement des moyennes k avec une graine aléatoire fixe pour assurer la reproductibilité des résultats.

Au lieu de modéliser directement les effets de la profondeur et de la profondeur au carré, nous avons d’abord normalisé la covariable de profondeur transformée en logarithme en soustrayant sa moyenne et en la mettant à l’échelle selon son écart type. Nous avons ensuite calculé la profondeur au carré à partir de cette variable centrée et mise à l’échelle. Cette méthode permet de s’assurer que les valeurs des covariables ne sont pas trop grandes ou trop petites pour éviter que des problèmes de calcul surviennent et que le centrage sépare les composantes linéaires et quadratiques du prédicteur.

Nous avons ajusté les quatre ensembles de données synoptiques séparément parce que seulement deux des relevés sont effectués chaque année et que les relevés sont disjoints dans l’espace et le temps. De même, nous ajustons indépendamment les relevés à la palangre sur fond dur nord et sud. Nous avons combiné les prédictions pour générer les graphiques cartographiques, mais nous avons marqué les années au cours desquelles les divers relevés ont été effectués.

À titre d’exemple, nous illustrons les composantes du modèle pour la morue du Pacifique dans le bassin de la Reine-Charlotte (figure \@ref(fig:sdmTMB-maps-combined)). Nous commençons par une couche bathymétrique et une valeur de densité de la biomasse pour chaque ensemble de relevés dans l’espace (figure \@ref(fig:sdmTMB-maps-combined)A). Après avoir ajusté le modèle, nous pouvons inspecter l’effet des prédicteurs d’effets fixes quadratiques de la profondeur du fond (figure \@ref(fig:sdmTMB-maps-combined)B) ainsi que les effets aléatoires spatiaux (figure \@ref(fig:sdmTMB-maps-combined)C). Si nous ajoutons les prédictions d’effets fixes aux effets aléatoires spatiaux dans l’espace du lien (logarithmique) et que nous élevons le résultat à une puissance, nous obtenons des prédictions de modèle qui incluent à la fois les effets fixes et aléatoires (figure \@ref(fig:sdmTMB-maps-combined)D). Nous pouvons inspecter les résidus de quantiles aléatoires dans l’espace pour vérifier s’il reste une corrélation spatiale (figure \@ref(fig:sdmTMB-resids)). Nous pouvons également examiner la relation prévue entre la profondeur et la densité de la biomasse pour toutes les espèces (figures \@ref(fig:sdmTMB-depth-all-plots1), \@ref(fig:sdmTMB-depth-all-plots2)).

```r avec $\nu = 1$. La ligne verticale en pointillés illustre $\sqrt{8\nu} / \kappa$, soit la 'plage', un point où la corrélation tombe en dessous de 0,1."} library(dplyr) library(ggplot2) dist <- seq(0, 13, length = 200) kappa <- c(0.3, 0.5, 1.1) nu <- 1.0 plyr::ldply(kappa, function(k) { data.frame( Correlation = k * dist * base::besselK(k * dist, nu), kappa = paste("kappa =", k), Distance = dist, range = sqrt(8 * 1) / k ) }) %>% mutate(Correlation = ifelse(is.na(Correlation), 1, Correlation)) %>% ggplot(aes(Distance, Correlation)) + geom_line() + facet_wrap(~kappa) + geom_vline(aes(xintercept = range), lty = 2) + coord_cartesian(expand = FALSE, ylim = c(0, 1)) + theme(strip.text.x = element_text(colour = "grey20", size = rel(1.1))) + labs(x = en2fr("Distance", french), y = en2fr("Correlation", french))

```r
library(sdmTMB)
d <- readRDS(here::here("report/data-cache/pacific-cod.rds"))
d <- d$survey_sets
dat <- gfplot:::tidy_survey_sets(d, "SYN QCS", years = 2017)
dat <- mutate(dat, density = density*1000*1000)
dat <- filter(dat, !is.na(depth))
dat <- gfplot:::scale_survey_predictors(dat)
dat <- select(dat, -X10, -Y10)

# grid_locs <- gfplot:::make_prediction_grid(filter(dat, year == 2017),
  # survey = "SYN QCS", cell_width = 1.5)$grid
# grid_locs <- rename(grid_locs, depth = akima_depth)
# grid_locs$year <- NULL

grid_locs <- gfplot::synoptic_grid %>%
  dplyr::filter(.data$survey == "SYN QCS") %>%
  dplyr::select(.data$X, .data$Y, .data$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 = 200)
options(OutDec =  ",")
plot_spde(spde)

m <- sdmTMB(
  data = dat, formula = density ~ depth_scaled + depth_scaled2,
  time = "year", spde = spde, family = tweedie(link = "log"),
  silent = TRUE, anisotropy = FALSE
)

predictions <- predict(m, newdata = grid_locs, return_tmb_object = TRUE)
plot_map <- function(.dat, column) {
  ggplot(.dat, aes_string("X", "Y", fill = column)) +
    geom_raster() +
    coord_fixed(expand = FALSE) +
    xlab(en2fr("Easting", french)) + ylab(en2fr("Northing", french)) +
    theme(legend.position = "bottom") +
    geom_point(data = dat, aes(X, Y, size = density), pch = 21, inherit.aes = FALSE) +
    scale_size_area() +
    guides(size = "none")
}

g_depth <- plot_map(grid_locs, "depth") +
  scale_fill_distiller(trans = "fourth_root_power", palette = "Blues", direction = 1,
    breaks = c(50, 200, 600)) +
  ggtitle(paste("(A) ", en2fr("Depth", french))) +
  labs(fill = paste(en2fr('Depth', french), ' (m)'))

limits <- c(0, max(exp(predictions$data$est))*1.001)
breaks <- c(1, 20, 150)
g_sdm1 <- plot_map(predictions$data, "exp(est)") +
  scale_fill_viridis_c(trans = "fourth_root_power", option = "C",
    breaks = breaks, limits = limits) +
  ggtitle(paste("(D) ", en2fr("Prediction (depth effects + spatial random effects)"))) +
  (if(french){
    labs(fill = "Estimation de la\nbiomasse (kg/m2)", parse = TRUE)
  } 
  else{
    labs(fill = "Biomass estimate\n(kg/km^2)", parse = TRUE)
    })

g_sdm2 <- plot_map(predictions$data, "exp(est_non_rf)") +
  ggtitle(paste("(B) ", en2fr("Prediction (depth effects only)", french))) +
  scale_fill_viridis_c(trans = "fourth_root_power", option = "C",
    breaks = breaks, limits = limits) +
  (if(french){
    labs(fill = "Estimation de la\nbiomasse (kg/m2)", parse = TRUE)
  } 
  else{
    labs(fill = "Biomass estimate\n(kg/km^2)", parse = TRUE)
    })

g_sdm3 <- plot_map(predictions$data, "omega_s") +
  ggtitle(paste("(C)", en2fr("Spatial random effects", french))) +
  scale_fill_gradient2(midpoint = 0,
    high = scales::muted("red"),
    mid = "white",
    low = scales::muted("blue")) +
  (if(french){
    labs(fill = "Écart par rapport à la biomasse\nprévue dans l’espace de relevé")
  } 
  else{
    labs(fill = "Deviation from expected\nbiomass in log space")
    })

\clearpage

cowplot::plot_grid(g_depth, g_sdm2, g_sdm3, g_sdm1)

dat$resids <- residuals(m) # randomized quantile residuals
# hist(dat$resids)
# qqnorm(dat$resids);abline(a = 0, b = 1)
ggplot(dat, aes(X, Y, col = resids)) + scale_colour_gradient2() +
  geom_point() + coord_fixed() +
  xlab(en2fr("Easting", french)) + ylab(en2fr("Northing", french)) +
  scale_colour_gradient2(midpoint = 0,
    high = scales::muted("red"),
    mid = "white",
    low = scales::muted("blue")) +
  labs(colour = en2fr('Residual', french))

\clearpage

g_models_dat <- readRDS(here::here('report/geostat-cache/spt-index-out-dat.rds'))

.colnames <- c(
      "Common name", "Survey", "$\\kappa$", "$\\tau_O$", 
      "$\\tau_E$", "$\\phi$", "$p$")
if (french) {
  .colnames[1] <- en2fr(.colnames[1], TRUE)
  .colnames[2] <- en2fr(.colnames[2], TRUE)
}
g_models_dat2 <- g_models_dat %>% 
  dplyr::arrange(species_name, survey) %>%
  # these have already been filtered out elsewhere for their massive CVs:
  dplyr::filter(phi < 1e3) %>% 
  dplyr::filter(tau_O < 50) %>% 
  dplyr::filter(tau_E < 50) %>% 
  mutate(species_name = gsub("Rougheye Blackspotted Rockfish Complex",
    "Rougheye/Blackspotted", species_name)) %>% 
  mutate(species_name = rosettafish::en2fr(species_name, french)) %>%
  mutate(species_name = gsub("Rougheye/Blackspotted Rockfish Complex", "Rougheye/Blackspotted Rockfish", species_name))

g_models_dat2$species_name <- purrr::map_chr(g_models_dat2$species_name , gfsynopsis:::cap)

g_models_dat2 %>%
  csasdown::csas_table(digits = 1, 
    col.names = .colnames,
    caption = "Estimations de paramètres à partir de modèles spatiotemporels géostatistiques.")