Add Discard CPUE plotting code so plots can be recreated in French

NAMESPACE

@@ -6,6 +6,7 @@ export(calc_naa)
 export(calc_paa)
 export(calc_prop_catch_by_fleet)
 export(calc_vb)
+export(cpue_extraction)
 export(expand_df_by_col)
 export(export_age_comps)
 export(export_catch)
@@ -19,10 +20,14 @@ export(format_mohn)
 export(get_bash_path)
 export(get_cvs)
 export(get_os)
+export(make_re_dat)
 export(mod_paths)
 export(modify_prop_female_iscam)
 export(plot_age_mat_ogives)
+export(plot_fixed_effect_coefs)
+export(plot_random_intercepts)
 export(plot_trawl_footprint)
+export(plot_year_locality_interactions)
 export(props_all)
 export(props_comm)
 export(props_comm_data_summary)
@@ -31,6 +36,7 @@ export(props_surv_data_summary)
 export(read_psv_file)
 export(rotate_coords)
 export(rotate_df)
+export(run_discard_index_model)
 export(run_models)
 export(system_)
 export(table_prop_female)

R/cpue-extraction.R

@@ -0,0 +1,100 @@
+#' Extract CPUE for Arrowtooth Flounder
+#'
+#' @details
+#' Used for the plots in the CPUE appendix in the 2022 assessment
+#'
+#' @param fn_cpue The file name of the cpue index data used
+#' @param ret_params If `TRUE`, return `params` list instead
+#'
+#' @return A data frame, or the `params` list, if `ret_params` is `TRUE`
+#' @export
+cpue_extraction <- \(fn_cpue = file.path(drs$nongit_dir,
+                                         "cpue-figs",
+                                         "arrowtooth-cpue-to-2024-10-21.rds"),
+                     ret_params = FALSE){
+    # fn_cpue_predict = file.path(
+    #   drs$nongit_dir,
+    #   "data",
+    #   "cpue-predictions-arrowtooth-flounder-modern-3CD5ABCDE-2024.csv"
+
+  params = list(
+    species_proper = "Arrowtooth Flounder",
+    area = c("^5A|^5B|^5C|^5D|^5E|^3C|^3D"),
+    area_name = c("3CD5ABCDE"),
+    skip_single_variable_models = FALSE,
+    use_alt_year = FALSE,
+    alt_year_start_date = "02-21",
+    final_year = 2021,
+    final_date = "2021-12-31",
+    discard_only = TRUE,
+    lat_range = c(48, Inf),
+    min_positive_tows = 140,
+    min_positive_trips = 10,
+    min_yrs_with_trips = 5,
+    depth_range = c(-Inf, Inf),
+    era = "modern",
+    parallel = TRUE
+  )
+
+  if(ret_params){
+    return(params)
+  }
+
+  # The file `fn_cpue` was generated using this command:
+  # gfdata::get_cpue_index(gear = "bottom trawl", min_cpue_year = 1996) |>
+  # filter(species_common_name = "ARROWTOOTH FLOUNDER)
+  comm_cpue <- readRDS(fn_cpue)
+
+  comm_cpue$fishing_event_id_unique <-
+    paste0(comm_cpue$vessel_registration_number,
+           "-",
+           comm_cpue$trip_id,
+           "-",
+           comm_cpue$fishing_event_id)
+
+  define_fleet <- \(area, area_name){
+    out <- gfplot::tidy_cpue_index(comm_cpue,
+                                   species_common = tolower(params$species_proper),
+                                   gear = "bottom trawl",
+                                   alt_year_start_date = params$alt_year_start_date,
+                                   use_alt_year = params$use_alt_year,
+                                   year_range = c(1996, params$final_year),
+                                   lat_range = params$lat_range,
+                                   min_positive_tows = params$min_positive_tows,
+                                   min_positive_trips = params$min_positive_trips,
+                                   min_yrs_with_trips = params$min_yrs_with_trips,
+                                   depth_band_width = 25,
+                                   area_grep_pattern = area,
+                                   depth_bin_quantiles = params$depth_bin_quantiles,
+                                   min_bin_prop = 0.001,
+                                   lat_band_width = 0.1)
+    out$area <- area_name
+    out
+  }
+
+  if(params$discard_only) {
+    comm_cpue <- comm_cpue |>
+      dplyr::filter(discarded_kg > 0,
+                    landed_kg == 0)
+  }
+  if(!is.null(params$final_date)) {
+    comm_cpue <- comm_cpue |>
+      dplyr::filter(best_date <= lubridate::ymd(params$final_date))
+  }
+
+  comm_cpue <- comm_cpue |>
+    dplyr::filter(best_depth >= params$depth_range[[1]],
+                  best_depth <= params$depth_range[[2]])
+
+  dfleet <- define_fleet(params$area, params$area_name)
+
+  # dfleet_sum <- dfleet |>
+  #   bind_rows() |>
+  #   group_by(area) |>
+  #   summarise(locality = length(unique(locality)),
+  #             depth = length(unique(depth)),
+  #             latitude = length(unique(latitude)),
+  #             vessel = length(unique(vessel)),
+  #             month = length(unique(month)))
+  dfleet
+}

R/make-re-dat.R

@@ -0,0 +1,22 @@
+#' Make a data frame from given GLMM model output showing the estimated
+#' parameter values for each parameter
+#'
+#' @param object A model as returned by [run_discard_index_model()]
+#'
+#' @return A data frame representing GLMM model parameter value outputs
+#' @export
+make_re_dat <- function(object) {
+
+  re <- glmmTMB::ranef(object)
+  plyr::ldply(re$cond, \(x) {
+    sud <- as.data.frame(x)
+    sud$par_value <- row.names(sud)
+    row.names(sud) <- NULL
+    sud
+  }) |>
+    rename(par_group = .id) |>
+    rename(est = `(Intercept)`) |>
+    as_tibble() |>
+    mutate(loc_group = gsub("^([0-9]+)[ -]*([0-9a-zA-Z-]+)$", "\\2", par_value)) |>
+    mutate(loc_year = gsub("^([0-9]+)[ -]*[0-9a-zA-Z-]+$", "\\1", par_value))
+}

R/plot-fixed-effect-coefs.R

@@ -0,0 +1,51 @@
+#' Plot fixed effect coefficients in the model run by [run_discard_index_model()]
+#'
+#' @param object A model as returned by [run_discard_index_model()]
+#'
+#' @return A [ggplot2::ggplot()] object
+#' @export
+plot_fixed_effect_coefs <- function(object){
+
+  su <- summary(object)$coefficients$cond
+  sud <- as.data.frame(su)
+  sud$param <- row.names(su)
+  row.names(sud) <- NULL
+  sud <- sud |>
+    rename(est = Estimate,
+           se = `Std. Error`) |>
+    mutate(par_value = gsub("^[A-Z_a-z]+", "", param)) |>
+    mutate(par_group = gsub("^([A-Z_a-z]+)[0-9.]+$", "\\1", param))
+
+  if(fr()){
+    sud <- sud |>
+      mutate(par_group = ifelse(par_group == "depth",
+                                "profondeur",
+                                ifelse(par_group == "latitude",
+                                       "latitude",
+                                       ifelse(par_group == "month",
+                                              "mois",
+                                              ifelse(par_group == "year_factor",
+                                                     "facteur_année",
+                                                     par_group)))))
+  }
+  sud |> ggplot(aes_string("est", "forcats::fct_rev(par_value)",
+                           yend = "forcats::fct_rev(par_value)")) +
+    ggplot2::geom_segment(aes_string(x = "est - 1.96 * se",
+                                     xend = "est + 1.96 * se"),
+                          lwd = 0.5) +
+    ggplot2::geom_segment(aes_string(x = "est - 0.67 * se",
+                                     xend = "est + 0.67 * se"),
+                          lwd = 1.25) +
+    geom_point() +
+    facet_wrap(~par_group,
+               scales = "free") +
+    theme_pbs() +
+    guides(shape = "none",
+                         colour = "none") +
+    labs(x = ifelse(fr(),
+                    "Valeur du coefficient (espace logarithmique)",
+                    "Coefficient value (log space)"),
+         y = ifelse(fr(),
+                    "Valeur prédictive",
+                    "Predictor value"))
+}

R/plot-random-intercepts.R

@@ -0,0 +1,49 @@
+#' Create a plot of random intercepts for the GLMM discard CPUE index model
+#'
+#' @param object A model as returned by [run_discard_index_model()]
+#' @param re_names A vector of names of GLMM model parameters to include
+#' in the plots
+#'
+#' @return A [ggplot2::ggplot()] object
+#' @export
+plot_random_intercepts <- function(object,
+                                   re_names = c("locality")){
+
+  re <- make_re_dat(object) |>
+    filter(par_group %in% re_names)
+
+  if(fr()){
+    re <- re |>
+      mutate(par_group =
+        ifelse(par_group == "depth",
+               "profondeur",
+               ifelse(par_group == "latitude",
+                      "latitude",
+                      ifelse(par_group == "month",
+                             "mois",
+                             ifelse(par_group == "year_factor",
+                                    "facteur_année",
+                                    ifelse(par_group == "locality",
+                                           "localité",
+                                           ifelse(par_group == "vessel",
+                                                  "navire",
+                                                  par_group)))))))
+  }
+
+  re |>
+    ggplot(aes_string("est",
+                      "forcats::fct_rev(par_value)")) +
+    geom_vline(xintercept = 0,
+               lty = 2,
+               alpha = 0.4) +
+    geom_point(bg = "white") +
+    facet_wrap(~par_group,
+               scales = "free") +
+    theme_pbs() +
+    guides(shape = "none",
+           colour = "none") +
+    labs(x = ifelse(fr(),
+                    "Valeur de l'ordonnée à l'origine aléatoire (espace logarithmique)",
+                    "Random intercept value (log space)"),
+         y = "")
+}

R/plot-year-locality-interactions.R

@@ -0,0 +1,29 @@
+#' Plot year-locality interactions for the GLMM model for Discard CPUE
+#'
+#' @param object A model as returned by [run_discard_index_model()]
+#'
+#' @return A [ggplot2::ggplot()] object
+#' @export
+plot_year_locality_interactions <- function(object) {
+
+  re <- make_re_dat(object) |>
+    filter(par_group == "year_locality")
+
+  re |>
+    ggplot(aes_string("as.numeric(loc_year)",
+                      "est",
+                      group = "loc_group")) +
+    geom_hline(yintercept = 0,
+               lty = 2,
+               alpha = 0.4) +
+    geom_point(alpha = 0.7) +
+    geom_line(alpha = 0.3) +
+    facet_wrap(~loc_group) +
+    scale_x_continuous(breaks = seq(1900, 3000, 10)) +
+    theme_pbs() +
+    guides(shape = "none", colour = "none") +
+    labs(x = tr("Year"),
+         y = ifelse(fr(),
+                    "Valeur de l'intercept aléatoire\n(interaction) (espace logarithmique)",
+                    "Random intercept\n(interaction) value (log space)"))
+}

R/run-discard-index-model.R

@@ -0,0 +1,104 @@
+#' Run the discard index GLMM model
+#'
+#' @param dfleet A list, as output from [cpue_extraction()]
+#' @param params A parameter list, as returned by [cpue_extraction()] with its
+#' argument `ret_params` set to `TRUE`
+#'
+#' @return The model object, and a CSV is written with the output
+#' @export
+run_discard_index_model <- \(dfleet,
+                             params,
+                             fn_csv = "/srv/arrowtooth/arrowtooth-nongit/data"){
+
+  dfleet$year_locality <- paste(dfleet$year_factor, dfleet$locality)
+
+  formulas <- tibble::tibble(
+    formula = c(
+      "cpue ~ 0 + year_factor",
+      "cpue ~ 0 + year_factor + depth",
+      "cpue ~ 0 + year_factor + month",
+      "cpue ~ 0 + year_factor + latitude",
+      "cpue ~ 0 + year_factor + (1 | locality)",
+      "cpue ~ 0 + year_factor + (1 | vessel)",
+      "cpue ~ 0 + year_factor + depth + month + latitude + (1 | locality) + (1 | vessel)",
+      "cpue ~ 0 + year_factor + depth + month + latitude + (1 | locality) + (1 | vessel) + (1 | year_locality)"
+    ),
+    formula_version = c(
+      "Unstandardized",
+      "Depth",
+      "Month",
+      "Latitude",
+      "Locality",
+      "Vessel",
+      "Full (without\nlocality-year effects)",
+      "Full standardization"
+    )
+  )
+
+  torun <- expand.grid(formula = formulas$formula,
+                       area = params$area_name,
+                       stringsAsFactors = FALSE)
+  torun <- torun |>
+    inner_join(formulas, by = "formula")
+
+  if(params$skip_single_variable_models){
+    torun <- torun |>
+      filter(formula_version %in% c("Unstandardized",
+                                    "Full standardization"))
+  }
+
+  if(params$parallel){
+    library("doParallel")
+    registerDoParallel(cores = floor(parallel::detectCores()/2L))
+  }
+  fn_model <- file.path(fn_csv,
+                        paste0("cpue-models-arrowtooth-",
+                               params$era,
+                               "-",
+                               params$area_name,
+                               "-",
+                               lubridate::as_date(now()),
+                               ".rds"))
+  fn_predictions <- file.path(fn_csv,
+                        paste0("cpue-predictions-arrowtooth-",
+                               params$era,
+                               "-",
+                               params$area_name,
+                               "-",
+                               lubridate::as_date(now()),
+                               ".csv"))
+
+  if(!file.exists(fn_model)){
+    if(params$discard_only){
+
+      fit_func <- \(dat, formula = cpue ~ year_factor, ...){
+        glmmTMB::glmmTMB(
+          as.formula(formula),
+          data = dat,
+          family = Gamma(link = "log"),
+          control = glmmTMB::glmmTMBControl(optCtrl = list(iter.max = 2000,
+                                                           eval.max = 2000),
+                                            profile = TRUE,
+                                            collect = FALSE),
+          ...)
+      }
+    }else{
+      fit_func <- gfplot::fit_cpue_index_glmmtmb
+    }
+    system.time({
+      model <- plyr::mlply(torun, \(formula, area, formula_version){
+        message("Fitting model ", formula)
+        fit_func(dfleet, as.formula(formula))
+      },
+      .parallel = params$parallel)
+    })
+    saveRDS(model, fn_model)
+  } else {
+    model <- readRDS(fn_model)
+  }
+
+  predictions <- plyr::ldply(model, predict_cpue_index_tweedie)
+  write_csv(predictions, fn_predictions)
+
+  model
+}