Chapter 8 Model specification
This section describes the implementation of the modules in FIMS in milestone 1. For the first milestone, we must implement enough complexity to adequately test a very standard population model. For this reason, we implement the minimum structure that can run the model described in Li et al. 2020.
8.1 Model variables and bounds
8.2 Inherited functors from TMB
8.2.1 Atomic functions
Wherever possible, FIMS should not reinvent atomic functions with extant definitions in TMB. If there is a need for a new atomic function the development team can add it to TMB using the TMB_ATOMIC_VECTOR_FUNCTION() macro following the instructions here.
8.3 Beverton-Holt expected recruitment function
For parity with existing stock assessment models, the first recruitment option in FIMS is a Beverton-Holt [cite] parameterized with R0 and h.
\[R_t =\frac{0.8R_0hS_{t-1}}{0.2R_0\phi_0(1-h) + S_{t-1}(h-0.2)}\] Where \(R_t\) and \(S_t\) are mean recruitment and spawning biomass in time \(t\), \(h\) is Mace-Doonan steepness, and \(\phi_0\) are the unfished spawning biomass per recruit. The initial FIMS model will implement a static spawning biomass-per-recruit function, with the ability to overload the method in the future to allow for time-variation in mortality, maturity, and weight-at-age over time to account for changes in spawning biomass per recruit. Deviations are assumed to be lognormally distributed such that realized recruitment is the product of mean recruitment and the exponentiated recruitment deviation. \[R_t = R_t\mathrm{exp}(r_{dev,t}-R^2/2), r_{dev,t} \sim N(0,R^2)\]
However, true \(r_{dev,t}\) values are not known, so when using estimated recruitment deviations \(\hat{r_{dev,t}}\) the following equation is applied to calculate mean unbiased recruitment \(R*_t\) using a bias adjustment factor \(b_y=\frac{E[SD(ry)]^2}{\sigma_R^2}\) (Methot and Taylor, 2011). \[R^*_t=R_t\mathrm{exp}(\hat{r_{dev,t}}-b_y\frac{\sigma_R^2}{2})\] The recruitment function should take as input the \(R\) , \(S\) values, the \(h\), \(ln(R_0)\), and R parameters and \(\phi_0\) and return mean and realized recruitment.
8.4 Logistic function with extensions
\[x_i=\frac{1}{1+\mathrm{exp}(-ln(19)(i-\nu_1)/\nu_2)}\]
Where \(x_i\) is the quantity of interest (proportion mature, selected, etc.), \(i\) is the index (can be age or size or any other quantity), \(\nu_1\) is the index of 50% mature and \(\nu_2\) is the difference in the index at 50% and the index at 95% (\(\nu_2 = \mathrm{ln}(19)/s\) where \(s\) is the slope parameter from an alternative parameterization). Logistic functions for maturity and selectivity should inherit and extend upon the base logistic function implementation.
8.5 Catch and fishing mortality
The Baranov catch equation relates catch to instantaneous fishing and natural mortality.
\[ C_{f,a,t}=\frac{F_{f,a,t}}{F_{f,a,t}+M}(1-\mathrm{exp}(-(F_{f,a,t}+M)))N_{a,t}\]
Where \(C_{f,a,t}\) is the catch at age \(a\) at time \(t\) for fleet \(f\), \(F\) is instantaneous fishing mortality, \(M\) is assumed constant over ages and time in the minimum viable assessment model, \(N_a,t\) is the number of age \(a\) fish at time \(t\).
\[F_{a,t}=\sum_{a=0}^A s_{a,f,t}F\]
\(s_a,f\) is selectivity at age \(a\) for fleet \(f\). Selectivity-at-age is constant over time.
Catch is in metric tons and survey is in number, so calculating catch weight (\(CW_t\)) is done as follows: \[ CW_t=\sum_{a=0}^A C_{a,t}w_a \]
Survey numbers are calculated as follows
\[I_t=q\sum_{a=0}^AN_{a,t}\] Where \(I_t\) is the survey index and \(q_t\) is survey catchability at time \(t\).
8.6 Modeling loops
This tier associates the expected values for each population section associated with a data source to that data source using a likelihood function. These likelihood functions are then combined into an objective function that is passed to TMB.
The population loop will be initialized at a user-specified age, time increment, and seasonal structure, rather than assuming ages, years, or seasons follow any pre-defined structure. Population categories will be described flexibly, such that subpopulations such as unique sexes, stocks, species, or areas can be handled identically to reduce duplication. Each subpopulation will have a unique set of attributes assigned to it, such that each subpopulation can share or have a different functional process (e.g. recruitment function, size-at-age) than a different category.
Spawning time and recruitment time are user-specified and can occur more than once per year. For the purposes of replicating model comparison project outputs, in milestone 1, all processes including spawning and recruitment occur January 1, but these should be specified via the spawn_time and recruit_time inputs into FIMS to allow for future flexibility. Spawning and recruitment timing can be input as a scalar or vector to account for multiple options.
Within the population loop, matrices denoting population properties at different partitions (age, season, sex) are translated into a single, dimension-folded index. A lookup table is computed at model start so that the dimension-folded index can be mapped to its corresponding population partition or time partition (e.g. population(sex, area, age, species, time, …)) so the programmer can understand what is happening. The model steps through each specified timestep to match the data to expected values, and population processes occur in the closest specified timestep to the user-input process timing (e.g. recruitment) across a small timestep that is a predefined constant.
8.7 Expected numbers and quantities
The expected values are calculated as follows in the population.hpp file: \[ B_t=\sum_{a=0}^AN_{a,t}w_a\] where \(B_t\) is total biomass in time \(t\), \(N\) is total numbers, \(w_a\) is weight-at-age \(a\) in kilograms.
\[N_t=\sum_{a=0}^AN_{a,t}\] where \(N_t\) is the total number of fish in time \(t\).
8.8 Initial values
8.8.1 Initial N
8.8.2 Initial F
8.9 Likelihood calculations
Age composition likelihood links proportions at age from data to model using a multinomial likelihood function. The multinomial and lognormal distributions, including atomic functions are provided within TMB.
Survey index likelihood links estimated CPUE to input data CPUE in biomass using a lognormal distribution. (model.hpp)
Catch index likelihood links estimated catch to input data catch in biomass using a lognormal distribution. (model.hpp)
Age composition likelihoods link catch-at-age to expected catch-at-age using a multinomial distribution.