Modeling time-varying growth in WHAM

R code

Giancarlo M. Correa and Jane Sullivan

Data

In the new WHAM extension (Correa et al. 2023), we can include the following data types:

  • Aggregated catch
  • Indices of abundance
  • Marginal age and length composition
  • Conditional age-at-length (CAAL)
  • Environmental covariates
  • Ageing error matrix
  • Observed mean weight-at-age

Model components

And the following model components:

  • Recruitment (abundance-at-age)
  • Selectivity
  • Natural mortality
  • Catchability
  • Environmental covariates
  • Somatic growth

Steps to implement your model in WHAM

Step 0: Install WHAM


If you want to use length information, you will need to install the developtment version of the growth branch:


remotes::install_github(repo = 'GiancarloMCorrea/wham', ref='growth')
# Then:
library(wham)

Step 1: Make your data object


You will make a list that includes all your input data and some relevant model information:


input_data = list()
input_data$ages
input_data$catch_paa # for fisheries 
input_data$index_paa # for indices
# And so on...

Step 2: Make your WHAM object


Here you will need to specify the model parameters:

my_input = wham::prepare_wham_input(model_name = "Example_1",
                         basic_info = input_data, 
                         NAA_re = list(...), # Recruitment parameters
                         M = list(...), # M parameter
                         selectivity = list(...), # Selectivity parameter
                         catchability = list(...), # Catchability parameter
                         ecov = list(...), # Environmental information
                         # Somatic growth:
                         growth = list(...), LAA = list(...), # Mean length-at-age
                         LW = list(...), # Length-weight relationship
                         WAA = list(...), # Mean weight-at-age
                         age_comp = "multinomial", # Age composition model
                         len_comp = "multinomial" # Length composition model
                         )

Step 3: Run model



Before running the model, you could also fix some parameters if you want (use your TMB knowledge).


my_model = wham::fit_wham(my_input, ...)

Step 4: Analyze results



There are some plots and tables that are generated by running:


wham::plot_wham_output(my_model)

Input data

Basic information


input_data = list()
input_data$ages = 1:10 # ages
input_data$lengths = seq(from = 2, to = 120, by = 2) # length bins
input_data$years = 1971:2020 # years
input_data$n_fleets = 1 # number of fisheries
input_data$n_indices = 1 # number of indices

Auxiliary information


Before continuing, let’s create some objects to make things easier later:

n_years = length(input_data$years)
n_ages = length(input_data$ages)
n_lengths = length(input_data$lengths)

Fishery information

# Aggregated catch:
input_data$agg_catch = matrix(my_data, ncol = input_data$n_fleets, nrow = n_years) # Obs
input_data$catch_cv = matrix(my_data, ncol = input_data$n_fleets, nrow = n_years) # Obs error
# Age comps:
input_data$catch_paa = array(my_data, dim = c(input_data$n_fleets, n_years, n_ages)) # Obs
input_data$catch_Neff = matrix(my_data, ncol = input_data$n_fleets, nrow = n_years) # Obs error
# Length comps:
input_data$catch_pal = array(my_data, dim = c(input_data$n_fleets, n_years, n_lengths)) # Obs
input_data$catch_NeffL = matrix(my_data, ncol = input_data$n_fleets, nrow = n_years) # Obs error
# CAAL data:
input_data$catch_caal = array(my_data, dim = c(input_data$n_fleets, n_years, n_lengths, n_ages)) # Obs
input_data$catch_caal_Neff = array(my_data, dim = c(n_years, input_data$n_fleets, n_lengths)) # Obs error
# Aging error matrix:
input_data$catch_aging_error = array(my_data, dim = c(input_data$n_fleets, n_ages, n_ages))

Index information

# Aggregated indices:
input_data$agg_indices = matrix(my_data, ncol = input_data$n_indices, nrow = n_years) # Obs
input_data$index_cv = matrix(my_data, ncol = input_data$n_indices, nrow = n_years) # Obs error
# Additional information:
input_data$units_indices = matrix(0L, nrow = n_years, ncol = input_data$n_indices) # 0 = numbers, 1 = biomass
input_data$fracyr_indices = matrix(0.5, ncol = input_data$n_indices, nrow = n_years) # fraccion del year
# Age comps:
input_data$index_paa = array(my_data, dim = c(input_data$n_indices, n_years, n_ages)) # Obs
input_data$index_Neff = matrix(my_data, ncol = input_data$n_indices, nrow = n_years) # Obs error
# Length comps:
input_data$index_pal = array(my_data, dim = c(input_data$n_indices, n_years, n_lengths)) # Obs
input_data$index_NeffL = matrix(my_data, ncol = input_data$n_indices, nrow = n_years) # Obs error
# CAAL data:
input_data$index_caal = array(my_data, dim = c(input_data$n_indices, n_years, n_lengths, n_ages)) # Obs
input_data$index_caal_Neff = array(my_data, dim = c(n_years, input_data$n_indices, n_lengths)) # Obs error
# Aging error matrix:
input_data$index_aging_error = array(my_data, dim = c(input_data$n_indices, n_ages, n_ages))

Selectivity pointers


We need to specify pointers for selectivity models (specified later):

input_data$selblock_pointer_fleets = matrix(1, ncol = input_data$n_fleets, nrow = n_years)
input_data$selblock_pointer_indices = matrix(2, ncol = input_data$n_indices, nrow = n_years)


This is simply a way to let WHAM know what selectivity model certain fishery or index will use (see parameters section later).

Weight-at-age

This could be empirical weight-at-age (without error) or observed mean weight-at-age (we would also need to include waa_cv):

input_data$waa = array(my_waa_data, dim = c(n_waa, n_years, n_ages))
input_data$waa_pointer_fleets = 1 #length = n_fleets. 
input_data$waa_pointer_indices = 2 #length = n_indices. 
input_data$waa_pointer_totcatch = 3 #length = 1 
input_data$waa_pointer_ssb = 4 #length = 1 
input_data$waa_pointer_jan1 = 5 #length = 1 
input_data$waa_cv = array(my_waaCV_data, dim = c(n_waa, n_years, n_ages)) # obs error

These pointers, like the selectivity pointers, is a way of indexation, which indicates what weight-at-age matrix a certain fishery or index should use.

Activate data


input_data$use_agg_catch = matrix(1, nrow = n_years, ncol = input_data$n_fleets) # 1 = use, 0 = no use
input_data$use_catch_paa = matrix(1, nrow = n_years, ncol = input_data$n_fleets) # 1 = use, 0 = no use
input_data$use_catch_pal = matrix(1, nrow = n_years, ncol = input_data$n_fleets) # 1 = use, 0 = no use
input_data$use_catch_caal = array(1, dim = c(n_years, input_data$n_fleets, n_lengths)) # 1 = use, 0 = no use
input_data$use_catch_waa = matrix(1, nrow = n_years, ncol = input_data$n_fleets) # 1 = use, 0 = no use
input_data$use_catch_aging_error = rep(1, times = input_data$n_fleets) # 1 = use, 0 = no use
input_data$use_indices = matrix(1, nrow = n_years, ncol = input_data$n_indices) # 1 = use, 0 = no use
input_data$use_index_paa = matrix(1, nrow = n_years, ncol = input_data$n_indices) # 1 = use, 0 = no use
input_data$use_index_pal = matrix(1, nrow = n_years, ncol = input_data$n_indices) # 1 = use, 0 = no use
input_data$use_index_caal = array(1, dim = c(n_years, input_data$n_indices, n_lengths)) # 1 = use, 0 = no use
input_data$use_index_waa = matrix(1, nrow = n_years, ncol = input_data$n_indices) # 1 = use, 0 = no use
input_data$use_index_aging_error = rep(1, times = input_data$n_indices) # 1 = use, 0 = no use

Additional information


input_data$maturity = matrix(my_data, nrow = n_years, ncol = n_ages) # maturity
input_data$fracyr_SSB = matrix(my_data, ncol = 1, nrow = n_years) # year fraction SSB
input_data$Fbar_ages = my_data # ages to calculate Fbar
input_data$percentSPR = 40 # % SPR reference point
input_data$percentFMSY = 100 # percent of FMSY to use for calculating catch in projections
input_data$simulate_period = c(1,0) # simulation for 1: model period, 2: projection
input_data$bias_correct_process = 1 # do process bias correction, 0 = no, 1 = yes
input_data$bias_correct_observation = 1 # do obs bias correction, 0 = no, 1 = yes

Somatic growth in WHAM

WHAM input


my_input = wham::prepare_wham_input(model_name = "Example_1",
                         basic_info = input_data, 
                         NAA_re = list(...), # NAA parameters
                         M = list(...), # M parameter
                         selectivity = list(...), # Selectivity parameter
                         catchability = list(...), # Catchability parameter
                         ecov = list(...), # Environmental information
                         # Somatic growth:
                         growth = list(...), LAA = list(...), # Mean length-at-age
                         LW = list(...), WAA = list(...), # Mean weight-at-age
                         age_comp = "multinomial", # Age composition model
                         len_comp = "multinomial" # Length composition model
                         )

WHAM input


Let’s focus only on the growth sections:

my_input = wham::prepare_wham_input(...,
                                    growth = list(model = ..., re = ...,
                                                  init_vals = ..., est_pars = ...,
                                                  SD_vals = ..., SD_est = ...), 
                                    LAA = list(re = ...,
                                               LAA_vals = ..., est_pars = ...,
                                               SD_vals = ..., SD_est = ...), 
                                    LW = list(re = ..., init_vals = ..., 
                                              est_pars = ...), 
                                    WAA = list(re = ..., WAA_vals = ...,
                                               est_pars = ...), 
                                    ...)

Parametric LAA


Corresponds to the growth argument.

  • model: (character) model to use.
    • vB_classic: von Bertalanffy (assumes \(\gamma = 1\)).
    • Richards: Richards (includes \(\gamma\)).
  • init_vals: (vector numeric) initial values for parameters \(k\), \(L_\infty\), \(L_{\tilde{a}}\), and \(\gamma\) (only if using Richards).

Parametric LAA


  • est_pars: (vector integer) index of parameter to be estimated (e.g., c(1,3) indicates that \(k\) and \(L_{\tilde{a}}\) will be estimated).

  • SD_vals: (vector numeric) initial values for parameters \(\sigma_{\tilde{a}}\) and \(\sigma_{A}\).

  • SD_est: (vector integer) index of SD parameter to estimate.

Parametric LAA


  • re: (vector character) random effects’ structure. Length equals to the number of parameters.
    • none: no deviations
    • iid_y: independent by year
    • iid_c: independent by cohort
    • ar1_y: correlated by year
    • ar1_c: correlated by cohort

Nonparametric LAA


Corresponds to the LAA argument.

  • re: (character) random effects’ structure.
    • none: no deviations
    • iid: independent over years and ages
    • 2dar1: correlated by years and ages.
    • 3dgmrf: correlated by cohorts, years and ages.

Nonparametric LAA


  • LAA_vals: (vector numeric) initial values of population mean size-at-age.

  • est_pars: (vector integer) indicates the index of the parameter of estimate.

  • SD_vals: (vector numeric) initial values of the \(\sigma_{\tilde{a}}\) and \(\sigma_{A}\) parameters.

  • SD_est: (vector integer) indicates the index of the SD parameter to estimate.

Semiparametric LAA


We have to specify both growth and LAA. From the LAA argument, we only need the re information.

Parametric WAA


Corresponds to the LW argument.

  • init_vals: (vector numeric) initial values for the \(\Omega_1\) and \(\Omega_2\) parameters.

  • est_pars: (vector integer) indicates the index of the parameter to estimate.

Parametric WAA


  • re: (vector character) random effects’ structure. Length is 2.
    • none: no deviations
    • iid_y: independent over years
    • iid_c: independent over cohorts
    • ar1_y: correlated over years
    • ar1_c: correlated over cohorts

Nonparametric WAA


Corresponds to the WAA argument.

  • re: (character) random effects’ structure.
    • none: no deviations.
    • iid: independent over years and ages.
    • 2dar1: correlated over years and ages.
    • 3dgmrf: correlated by cohorts, years and ages.

Nonparametric WAA


  • WAA_vals: (vector numeric) initial values of population mean weight-at-age.


  • est_pars: (vector integer) indicates the index of parameters to estimate.

Simulation experiment

Simulation experiment

Goal: provide guidelines for growth modeling in state-space assessment models under diverse scenarios.

  • Data type: age compositions vs length compositions vs CAAL (random vs length-stratified sampling).
  • Informative data source: fishery vs survey.
  • Data quality: data rich vs data poor.
  • Modeling approach: parametric vs nonparametric vs semiparametric vs Ecov.
  • Time-varying parameter: changes in \(k\), \(L_{\infty}\) or \(L_{\tilde{a}}\).

Preliminary results

Thanks


Collaborators: Cole Monnahan, James Thorson, Andre Punt, Tim Miller, Steve Barbeaux, Peter Hulson


Contact:
gmoron@azti.es

Find these slides at:
giancarlomcorrea.netlify.app

References

© AZTI 2023. Todos los derechos reservados.

Correa, Giancarlo M, Cole C Monnahan, Jane Y Sullivan, James T Thorson, and André E Punt. 2023. “Modelling Time-Varying Growth in State-Space Stock Assessments.” ICES Journal of Marine Science 80 (7): 2036–49. https://doi.org/10.1093/icesjms/fsad133.