In CWetzel/XSSS: Extended Simple Stock Synthesis

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

Summary

xsss is a package for Extended Simple Stock Synthesis which was designed for data-limited assessments that use catch and indices of abundance data. XSSS is an extension to Simple Stock Synthesis which calculates a time-series of abundance based on priors for natural mortality, steepness, and stock status in a specific year (Simple Stock Synthesis can be found here: https://github.com/shcaba/SSS). XSSS applies priors for natural mortality, steepness, and stock status, but uses adaptive importance sampling to update the priors based upon index of abundance data to estimate posterior distributions for each parameter resulting in time-series estimates based on information in the index data.

Software requirements

xsss is designed for use with Stock Synthesis version 3.30.11 or greater. The executable and additional information regarding Stock Synthesis is available on NOAA Vlab (https://vlab.ncep.noaa.gov/) which is accessible to all NOAA employees. Access is also available to outside collaborators after creation of a user account. If you are not a NOAA employee and would like to join the Stock Synthesis group on Vlab please email [e-mail].

Setting up a model

xsss works using Stock Synthesis which allows the user to create a model with multiple fisheries and surveys if appropriate. Stock Synthesis is an age-based population dynamics model. The user needs to provide four input files for Stock Synthesis: starter, forecast, data, and control file. An advantage for using Stock Synthesis is the flexibility for the user to define simple of complex biology (e.g., sex specific growth, fecundity relationships) and fishery dynamics (e.g., selectivity of fisheries and surveys).

File structure

xsss will run inside the folder specified by the filepath input to the function (See below for specific information on function inputs). Inside this folder the user will need to place each of the four model files required by Stock Synthesis and the executable. The control and data file names are specified by the user when calling the function. The starter and forecast files need to be named exactly starter.ss and forecast.ss for the package to recognize the files. A "run" folder will be created inside the filepath folder and all files will be copied into the run folder preserving the original input files. The ss_summary.sso files will be saved from the final model run in the report folder inside the run folder. Saved summary files of model results, summary plots, and summary tables are avialable in the save folder.

Data file

The data file used for xsss mimics a data file that would be used for a typical Stock Synthesis except in a couple notable locations. The number of fleets specified in the model needs to be changed to N+1 to include an additional depletion survey. The depletion survey is defined similar to other surveys in Stock Synthesis where a list of fishery and survey fleet information is expected prior to the catch data where:

  library(xtable)
  library(kableExtra)

  tab = rbind( c(1, -1, 1, 1, 0, "Fishery"),
               c(3, 1007, 1, 1, 0, "Survey" ),
               c(3, 1001, 1, 1, 0, "Depletion_Survey"))

  colnames(tab) = c("Fleet_type", "Timing", "Area", "Units", "Catch_mult", "Fleet_name")

  table = tab 
  kable(table, "html")

The above example is for a model with one fishery and one survey with an index of abundance.

Specification for the depletion survey must also be defined at the top of index of abundance data section. The depletion survey is define using option 34 in the "Units" column for Stock Synthesis.

  library(xtable)
  library(kableExtra)

  tab = rbind( c(1, 1, 0, 1, "# Fishery"),
               c(2, 1, 0, 1, "# Survey" ),
               c(3, 34, 0, 1,"# Depletion_Survey"))

  colnames(tab) = c("Fleet_number", "Units", "Error_Type", "Extra_SD_Reported", "Fleet_name")

  table = tab 
  kable(table, "html")

The depletion survey should include two data points, one for the first year of the model where the stock is assumed to be a specific fraction of the unfished biomass and a second for the depletion value for a specific model year. The depletion survey should be entered at the bottom of the index of abundance data section as:

  library(xtable)
  library(kableExtra)

  tab = rbind( c(1918, 1, 3, 1, 0.001, "# Initial Depletion"),
               c(2000, 1, 3, 0.40, 0.001, "# Final Depletion"))

  colnames(tab) = c("Year", "Month", "Fleet", "Observation", "Standard_Error", "Fleet_name")

  table = tab 
  kable(table, "html")

where the above example has the depletion survey as fleet 3 where the stock is assumed in an unfished state in 1918 and depleted to 40% of unfished biomass in year 2000. The model is forced to fit the depletion survey by setting a very small standard error value (>0.001). xsss will plug in alternative depletion values for the final depletion year based on draws from either the prior or posterior distribution.

No additional edits for the depletion survey is required within the data file. However, all additional required inputs for the data file (e.g., discard, length and age bins) need to be specified. Please see example files within the example folder.

Control file

The control file for xsss is nearly identical to a control file for a full Stock Synthesis based model. The control file contains the biological parameters, stock recruitment relationship, and selectivity. xsss only estimates a single parameter, log(R0), which is solved for by the model based upon the fixed parameter values, the catch history, and the depletion survey value. The estimation phase for all other parameters in the control file need to be set to a negative value to avoid estimation.

xsss modifies the fixed parameters for natural mortality (both sexes) and steepness. The code currently matches a string in order to correctly identify the parameter lines that need to be modified. In order for the package to correctly identify these parameter lines they need to be labeled with the default Stock Synthesis labeling created in the ss_new files. The expected strings are:

• NatM_p_1_Fem_GP_1 for female natural mortality,

• NatM_p_1_Mal_GP_1 for male natural mortality, and

• SR_BH_steep for steepness.

xsss currently is only configured for a constant natural mortality value across ages for a single area based model. The parameter offset approach between sexes for biological parameters should not be selected when setting up a control file for xsss. If natural mortality is selected to be equal between sexes, the package still expects that each parameter be set up for each sex where the code will fix the parameter equal for each at the drawn MCMC value.

xsss assumes deterministic recruitment with recruitment being estimated directly from the selected stock recruitment curve. Hence, recruitment deviation estimation in Stock Synthesis needs to be turned off ("0 #do_recdev"). Additionally, all phases for recruitment deviation estimates need to be set to a negative phase.

The model required the depletion survey to have defined parameters for Q and selectivity similar to all other surveys in the model. An example set-up for the Q parameter section is:

  library(kableExtra)

  tab = rbind( c(2, 1, 0, 0, 0, 1, "# Survey"),
               c(3, 1, 2, 0, 0, 1, "# Depletion Survey"))

  colnames(tab) = c("Fleet", "Link", "Link_info", "Extra_SE", "Bias_adj", "Float", "Fleet_name")

  table = tab 
  kable(table, "html")

The "Link_info" for the depletion survey should be set to option 2, which requires to model to run through all phases when using this special survey option. Additionally, the "Float" option needs to be set to option 1 for all surveys which allows the model to analytically calculate the Q value for each survey, where the depletion survey Q value will be equal to 1. In the example above, two additional parameter lines are required for the Q parameter section for the estimation of ln(Q) parameter.

Selectivity must be defined for each fleet in the model, including the depletion survey. Since xsss uses Stock Synthesis the user is able to define the selectivity curve that is most appropriate for the fishery or survey (e.g., selectivity equal to maturity, dome-shaped). An example length-based selectivity pattern block set-up is:

  library(kableExtra)

  tab = rbind( c(1, 0, 0, 0, "# Fishery"),
               c(15, 0, 0, 1, "# Survey"),
               c(10, 0, 0, 0, "# Depletion Survey"))

  colnames(tab) = c("Pattern", "Discard", "Male", "Special", "Fleet_name")

  table = tab 
  kable(table, "html")

where this specification sets the fishery selectivity to a two-parameter logistic curve and the survey selectivity mirrors (Pattern = 15) the fishery selectivity (Special = 1, the fishery fleet number). The depletion survey is set to pattern 10 where all fish greater than the first length bin have are fully selected. Alternatively, age-based selectivity may be used for fisheries and surveys and would be defined in a similar fashion using the age-selectivity pattern block. All selectivity parameters should be fixed at the input parameter values, and hence should be considered carefully.

Starter file

The starter file is modified by xsss. The user is required to set-up an initial starter file. However, xsss performs checks and modifies input values as required (e.g., turn off report file writing for speed, set parameter jitter equal to 0).

Forecast file

xsss does not do any adjustments to the forecast file. The user needs to carefully consider the set-up of this file. Within the forecast file, the target harvest rate, the control rule, and adjustments to the estimated overfishing limit is defined by the user. Additionally, harvest may be fixed at a pre-determined value for select years during the forecast period to reflect a management system where harvest limits have already be set by management for specific years. An example set-up for a West Coast groundfish stock assessment is provided in the example folder.

Visualizing model set-up

Once all files have been set-up by the user the Stock Synthesis executable should be run to ensure that all files have been defined properly. The executable may be run from a command window using the command "ss -nohess". If the model fails to run, the echoinput file can be used to determine where the reading of the input files failed. After configuring a running model the user should confirm that only the log(R0) parameter is estimated by the model. The results from a model run can be visualized using the r4ss package which can be located here: https://github.com/r4ss. Once the r4ss package is installed the model may be read into R and visualized by:

model = SS_output("C:/My_Directory/.../")
SS_plots(model)

Call to the xsss function

The following inputs are required when calling xsss:

SSS.ais.fxn(filepath = "C:/My_Directory/.../",
            control.name = "control_file.ctl", 
            dat.name = "data_file.dat", 
            m.in = c(0.28, 0.27,0.40,0.40, 0), 
            h.in = c(0.87, 0.093, 0.20, 1, 1), 
            depl.in = c(0.60, 0.10, 0.01, 0.99, 1) #
            )

where:

• filepath is the parent directory which contains the Stock Synthesis executable and each of the four model files. The function will copy and move each of these items into a folder named run where the model will run, preserving the original files in the filepath location

• control.name is the name of the control file. This will be used to locate the file to copy and move to the run folder. Additionally, the code will identify this file by name to modify values of natural mortality and steepness for each model iteration.

• dat.name is the name of the data file. This will be used to locate the file to copy and move to the run folder. Additionally, the code will identify this file by to modify values of relative stock status (depl.in) in a fixed model year.

• m.in is a vector of median natural mortality, standard deviation, and whether female and male natural mortality should be assumed equal. The expected input to the function is c(female M median value, male M median value, female M sd, male M sd, M equal between sexes (0=No, 1=Yes)).

• h.in is a vector of mean steepness, standard deviation, and lower and upper bounds for the truncated beta distribution. The expected input to the function is c(steepness mean, steepness sd, lower bound, upper bound, distribution (1= truncated beta)). Currently, only the truncated beta distribution is available for the steepness parameter.

• depl.in is a vector of mean relative stock status (aka "depletion" from unfished conditions), the standard deviation, and the lower and upper bound for the truncated beta distribution. The expected input to the function is c(depletion in select year, sd, lower bound, upper bound, distribution (1 = truncated beta, 2 = lognormal)). If the lognormal distribution is selected the lower and upper bounds for parameters will not be used but still need to be provided.

The function has default values for other inputs to the function that can be overridden by the user. The full function call will all inputs is:

SSS.ais.fxn(filepath = "C:/My_Directory/.../", 
            control.name = "control_file.ctl", 
            dat.name = "data_file.dat", 
            tantalus = FALSE, 
            read.seed = FALSE, 
            entropy.level = 0.92, 
            Niter = 2000, 
            AIS.iter = 2000, 
            final.Niter = 5000, 
            m.in = c(0.28, 0.27,0.40,0.40, 0), 
            h.in = c(0.87, 0.093, 0.20, 1, 1), 
            depl.in = c(0.60, 0.10, 0.01, 0.99, 1), 
            fmsy.m.in = NULL, 
            bmsy.b0.in = NULL 
            ) 

where

• tantalus is an option for use with linux based machines. This has only been tested on the NWFSC system and may not work in general application.

• read.seed is an option allowing the user to use the same seeds from a previous run. The package will read the seed_list file to set seeds.

• entropy.level is the level of entropy that must be met prior to finishing adaptive importance sampling.

• Niter is the number of initial model runs for MCMC drawing from the prior distributions.

• AIS.inter is the number of model runs during the adaptive importance sampling where the priors are updated based upon the index data.

• final.Niter is the number of final model runs drawing from the posterior distributions for each parameter.

• fmsy.m.in and bmsy.b0.in inputs are currently not fully implemented.

Examples

Example function calls and model files are located in the examples folder (https://github.com/CWetzel/XSSS/tree/master/examples).

CWetzel/XSSS documentation built on May 22, 2019, 1:53 p.m.