In NOAA-EDAB/neusCatch: Multispecies length expansion of commercial catch

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

source("../data-raw/R/create_species_rules_table.r")

When working with mscatch, depending on what you want to do, you will need to:

Assemble data for use in the package
Aggregate the catch data based on a set of rules
Use length-weight relationships to expand the catch to weights using the length samples
Create an Age-Length Key
Calculate numbers at age

Assuming you have assembled the data (for a sample data set see, get_sample_data), the next step is to determine how the data should be aggregate and how to deal with missing length samples in YEARs/QTRs where there are catch.

There are two approaches a user can take:

Define specific rules in a species object list
Specify general rules for a more exploratory approach.

Define a species object

A species object can be defined and passed as an argument to the aggregate_landings function. This object specifies decision rules to be used in the aggregation step that answer the following questions:

Which gear types should be aggregated?
Which market codes should be aggregated?
How to deal with with missing length samples? where to borrow from if at all?
Which statistical areas define the stock?
Where is length-weight data pulled from to fit length-weight relationships? survey database, commercial database, both?
At what time scale are length weight relationships fitted? Semester, quartery, annually?
Where is the age data pulled from?
How should age length keys be calculated? Annually, by semester annually, by gear, semester, annually?
What is the maximum age of the species ?

General rules

This option is for exploring the data and assessing the impacts related to aggregation decisions to help inform building a species object. For example:

To explore the most important gear types that are associated with landing a species a user can specify a threshold value relating to the proportion of cumulative landings, over which all other gear types are considered insignificant and are grouped together.
Market codes are aggregated if the length frequency distributions are found to be significantly insignificant (using the Kolmogorov-Smirnov test). The level of the test (p-value) can be passed as an argument.

Examples

Using the exploratory approach

Assigns the name of species (speciesName = "Haddock")
Aggregates landings over QTR year to YEARly data and combines all length samples to same level of aggregation (aggregate_to = "YEAR").
Uses gear types (NEGEARs) which contribute the most to total landings. Uses a cumulative landings threshold of 90% to select gears. Other gears will be aggregated into an "other Gear" category (landingThresholdGear = 0.9)
Assigns other gears to code "998" (otherGear = "998")
Borrows length samples from nearest neighbor, whether by QTR, YEAR, or NEGEAR (borrowLengths = T).
Minimum number of length samples required when borrowing samples from other NEGEARs, QTRs, or YEARs (nLengthSamples = 1)
Assigns an output folder (outputDir=here::here("output/Haddock")) for all diagnostic plots and the logfile
Assigns a name to the log file produced (logfile = "logfile.txt"). This documents all decisions made in aggregation steps and documents the years in which length samples were borrowed and the year from which they were borrowed

 aggregatedData <- mscatch::aggregate_landings(channel,
                                               landingsData,
                                               lengthData,
                                               outputDir=here::here("output/Haddok"),
                                               speciesName = "Haddock"
                                               landingThresholdGear = 0.9,
                                               nLengthSamples = 1,
                                               aggregate_to = "YEAR",
                                               borrowLengths = T,
                                               otherGear = "998",
                                               logfile = "logfile.txt"
                                               outputPlots = T,
                                               speciesRules = NULL)

At this point, aggregateData contains length samples for every combination of YEAR, NEGEAR, MARKET_CODE. All decisions made regarding gear type selection, market code selection, and which YEAR length samples were borrowed from are recorded in the logfile. If outputPlots=T a suite of plots are also output in the outputDir.

Using a species object

A template for the species object is bundled with mscatch which can then be adapted

# Haddock species itis code = 164744
species_itis <- 164744
speciesRules <- mscatch::get_species_object(species_itis)
speciesRules$marketCodes  |>  
  flextable::flextable()

This specifies the market code aggregation; LG and XG codes are combind and relabeled as LG, SR and SK remain as they are. All other market codes are treated as unclassifieds, UN.

speciesRules$gearCodes |> 
  flextable::flextable()

All gear codes are aggregated and labelled as gear code "050"

Other related rules:

create_species_rules_table(spitis = species_itis)

This speciesRules object is then passed as an argument to aggregate_landings

Expand the catch (Length & Numbers)

To expand the catch to weight-by-length a length-weight key is required. This is achieved by obtaining parameter values that satisfy the length-weight relationship. This can be obtained in one of two ways:

Pull parameter values from Survey database (svdbs) using function get_length_weight function from the survdat package.

Pull raw length-weight data and fit own model. The bottom trawl survey data is used (for this example). A series of nested length-weight models are then fit to the data.

The data is pulled using the function get_length_weight_data from the survdat package and fit using fit_length_weight()

# Permissions and credentials are required to conect to the database
channel <- dbutils::connect_to_database("server","user")
# Haddock data are pulled (SVSPP = 74)
lengthWeightData <- survdat::get_length_weight_data(channel,year="all", species=74)
# Length-Weight models are thenn fit
fits <- mscatch::fit_length_weight(lengthWeightData$data,speciesName,speciesRules,outputDir,logfile="logfile.txt")

The fitted model is then used to expand the catch. This is achieved by using expand_landings_to_lengths()

  expandedLandings <- mscatch::expand_landings_to_lengths(aggregatedData$landings,aggregatedData$lengthData,fits$paramsH0) |> 
    dplyr::mutate(speciesName = speciesName)