Index selectivity and catchability

This tutorial starts off assuming that users are familiar with the setup with fishery selectivity.

Unlike fleet selectivity, index selectivity is unique to each series and no dummy fleets are used.

Is index selectivity already defined elsewhere?

First, we consider the argument vector s_selectivity, which defines where the index selectivity is defined (the terms index and survey are used interchangeably). Index selectivity may be identical to fleet selectivity, i.e., vulnerable biomass, or could be related to spawning or total biomass. If we have 4 indices with:

s_selectivity <- c("SSB", "B", 1, 2)

The first index is an index of spawning biomass as denoted by “SSB” (maturity is configured in the OM), the second index tracks total biomass as denoted by “B” (selectivity = 1 for all ages), and the third and fourth indices have the selectivity of the first and second fleets, respectively. Integers for fleets refer to the true fleet and not to selectivity blocks/dummy fleets.

In these cases, index selectivity is already defined somewhere else in the model, and no further consideration of index selectivity is needed. The RCM function call can look like this:

output <- RCM(OM, RCMdata, selectivity = selectivity, 
              s_selectivity = s_selectivity, 
              ...)

Index selectivity is independent of anything else in the model

On the other hand, if index selectivity needs to be explicitly defined, then the s_selectivity vector can indicate the functional form, using one of logistic, dome, or free. Let’s look at another situation with 5 indices with this selectivity definition:

s_selectivity <- c(2, "SSB", "B", "dome", "free")

For the fourth and fifth index, the selectivity functions are dome-shaped and free parameters, respectively. We will need to explicitly consider what the parameters defining this functions are, either as starting values to be estimated or fixed values in the model.

Selectivity parameters

Just as with the fleet selectivity parameters, the index selectivity parameters by default use OM@LFS, OM@L5, and OM@Vmaxlen for start values when s_selectivity = "logistic" or "dome". Custom start values are needed when index selectivity uses free parameters.

Custom start values are passed to the RCM in the start$ivul_par matrix with the same dimension and layout as for start$vul_par:

OM@maxage <- 5
s_selectivity <- c(2, "SSB", "B", "dome", "free")
print(ivul_par)

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0    0 55.0    1
## [2,]    0    0    0 40.0    0
## [3,]    0    0    0  0.5    0
## [4,]    0    0    0  0.0    0
## [5,]    0    0    0  0.0    0
## [6,]    0    0    0  0.0    0

The parameter slots in the first three columns (indices 1-3) are ignored. Recall that selectivity is defined elsewhere (through fishery selectivity, maturity, and total biomass, respectively) and the model will ignore these columns in ivul_par. Any value can be placed here, but users must ensure that the columns in the corresponding map matrix are set to NA so that TMB will turn off estimation of these parameters.

In column four, the first three rows are the start values for the three parameters of the dome function. The last three rows are not used and can be any value. In column five, we set up the selectivity only for age-0 animals, i.e., an index of recruits.

Recall that the map argument tells TMB what to do to the parameter of the corresponding row and column in ivul_par. Shared parameters are assigned a unique integer amongst themselves while fixed or unused parameters are assigned NA.

For the map matrix, we do the following things:

Columns 1-3 are NA because index selectivity is defined outside the ivul_par matrix. TMB turns off estimation of these parameters.
Column 4 sets up three independent parameters so that all three parameters for the dome are estimated. The last three rows in ivul_par are not used and the map should be set to NA.
Column 5 is all NA because we fix the selectivity at age schedule.

print(map_ivul_par)

##      [,1] [,2] [,3] [,4] [,5]
## [1,]   NA   NA   NA    1   NA
## [2,]   NA   NA   NA    2   NA
## [3,]   NA   NA   NA    3   NA
## [4,]   NA   NA   NA   NA   NA
## [5,]   NA   NA   NA   NA   NA
## [6,]   NA   NA   NA   NA   NA

A function call utilizing this custom set-up for index selectivity would be:

output <- RCM(OM, data, selectivity = selectivity, 
              s_selectivity = s_selectivity, 
              start = list(ivul_par = ivul_par), 
              map = list(ivul_par = map_ivul_par), 
              ...)

Index catchability

By default, each index series is time-invariant with respect to the selectivity and catchability ($q$). By default, catchability for index $s$ is solved analytically in the model,

\[ q_s = \exp\left(\dfrac{\sum_y \log(I^{\textrm{obs}}_{y,s}) - \sum_y \log(x_{y,s})}{n_s}\right) \]

where $I$ is the observed index, $x$ is the biomass or abundance vulnerable to the index, and $n_s$ is the number of data points.

This section goes through the three use cases if one or both of selectivity and catchability changes.

Case 1 - Both catchability and selectivity change

Simply break the index into two series and model them separately.

Case 2 - Catchability changes but selectivity does not

Break the index into two series and map the selectivity parameters so that the two series share the same parameters (see selectivity section above).

Case 3 - Catchability is constant but selectivity changes

Break the index into two series. Then use the map argument to identify the index series which share a common $q$. As an example with 5 index series,

RCM(..., map = list(q = c(1, 1, 2, 2, NA)))

This map argument sets up the model so that index series 1 & 2 share one $q$ parameter and series 3 and 4 share another $q$ parameter. When map is an integer, $q$ is an estimated parameter while NA for the fifth index indicates that $q$ will be calculated analytically.

Fixing catchability

RCM provides an option for fixing catchability to one:

RCMdata@abs_I <- c(0, 0, 0, 0, 1)

where $q = 1$ for the fifth index.

There are two ways to fix $q$ to an alternative value. First, rescale the index such that the new index will have $q=1$. Second, set the q prior to with a tight prior:

prior <- list()
prior$q <- matrix(NA, 5, 2)
prior$q[5, ] <- c(0.5, 0.001)
  
RCM(..., prior = prior)

where the $q$ for the fifth index has a prior mean of 0.5 and lognormal standard deviation of 0.001.

Offset catchability

If two index series have offset catchability, e.g., the $q$ of one index is twice of another index, simply re-scale one of them and share the $q$ parameter between the two.