Description of the simulation models in cpueTargeting

Background document: (A) the spatial multispecies simulator, simulate_bathymetric(), and (B) its extension to a longline tuna fishery with vertical targeting, simulate_tuna(). Henning Winker.

1. Origin and acknowledgement

The spatial operating model in cpueTargeting is based on the simulation model described in Appendix A of Thorson et al. (2017) and implemented as the function Generate_scenario_data in the spatial_DFA R package (github.com/James-Thorson/spatial_DFA). That simulator is itself a generalisation of the simulator developed by Ono et al. (2015) for defining spatial strata in catch and effort analyses. It tracks the biomass of several species in each cell of a spatial grid over time, with logistic population growth, habitat-driven movement and a fleet whose effort follows expected revenue, so that changes in species prices generate changes in fisher targeting.

A modified version of that code ("Run 22") was used to generate the simulated CPUE data for the fisheries-dependent CPUE standardisation guidelines prepared for the ICES WKBMSYSPICT benchmarks (Winker and Cardinale 2022). simulate_bathymetric() is a re-implementation of Run 22 as a package function. simulate_tuna() extends the same population and fleet dynamics with a vertical dimension for a longline tuna fishery. The original Run 22 scripts are kept unchanged in the package folder inst/extdata for reference.

2. Appendix A: the spatial multispecies simulator (simulate_bathymetric())

The model tracks the biomass density \(d_p(s,t)\) of species \(p = 1,\dots,P\) in grid cell \(s = 1,\dots,S\) and year \(t = 1,\dots,T\). The default (Run 22) uses \(P = 4\) species on a \(50 \times 50\) grid (\(S = 2500\)) over \(T = 30\) years.

A.1 Habitat: bathymetry

Depth \(z(s)\) is simulated as a stationary Gaussian random field with standard deviation \(\sigma_z\) (SD_O = 150 m) and Gaussian covariance \(\exp\{-(\Delta s/\phi)^2\}\) with correlation scale \(\phi\) (SpatialScale = 4 cells), as in the original code (RandomFields::RMgauss). A constant is added so that all depths are positive.

A.2 Movement and habitat preference

The probability that a fish of species \(p\) moves from cell \(s^*\) to cell \(s\) during a year is

where \(\Delta_{s^*s}\) is the distance between cell centres, \(\lambda_p\) the mean movement distance (Fish_dist_par1), and \(c_p(s^*)\) normalises each row of the movement matrix \(\mathbf{M}_p\) to sum to one. The depth preference \(h^{(z)}_p\) is a lognormal density with mean \(\mu_p(z)\) (Fish_depth_par1) and coefficient of variation \(\mathrm{CV}_p(z)\) (Fish_depth_par2); the range preference \(h^{(x)}_p\) is uniform between Fish_range_par1 and Fish_range_par2 (0 and 50 by default, which covers almost the whole grid).

A.3 Population dynamics

In each cell, biomass follows a Schaefer (logistic) model with catch removed, after which fish redistribute according to the movement matrix:

Initial biomass is the stationary distribution of the movement process scaled to \(0.5\,B_{0,p}\,\varepsilon_p\) with \(\varepsilon_p \sim \mathrm{LN}(0, 0.2)\), and the carrying capacity of each cell, \(k_p(s)\), equals this initial biomass, so stocks start at equilibrium and trends are driven by fishing.

A.4 Fishing effort

The number of sets (trips) in year \(t\) follows a logistic increase with lognormal error (CV = 0.2):

with \(U_T\) = Tot_effort (500 in Run 22) and \(\sigma_U^2 = \log(1 + 0.2^2)\). Sets are allocated to cells in proportion to expected revenue,

where \(\phi_p(t)\) is the price of species \(p\) in year \(t\) (price_fish), \(q\) the mean catchability and \(\kappa\) (Preference, 1 by default) controls how strongly effort concentrates on high-value cells. Price changes therefore shift where the fleet fishes and hence which species it targets. Each set is assigned to a vessel drawn at random from \(V\) = 30 vessels.

A.5 Catch

Each vessel has a mean catchability \(q_v \sim \mathrm{LN}(\log q - \sigma_v^2/2,\, \sigma_v^2)\) with CV 0.2 (\(q\) = qq_original = 0.1). The catchability of set \(u\) for species \(p\) has Tweedie observation error, \(q_p(u) \sim \mathrm{Tweedie}(\text{mean } q_{v(u)},\ \text{power } 1.2,\ \text{dispersion } 0.1)\). All sets in a cell share the biomass of that cell: the total catch in cell \(s\) is

and is split among the sets in proportion to their catchability, \(c_p(u) = c_p(s,t)\, q_p(u) / \sum_{u' \in s} q_p(u')\). The catch per set is the CPUE record.

A.6 Data and truth

For each set the output (Data) contains year, vessel, true cell position, reported position, bottom depth and the catch of each species. As in the guidelines, reported positions are the centres of \(5 \times 5\) blocks of cells (obs_res), i.e. coarser than the 1 x 1 cells on which the population is simulated. The true index of abundance is the total biomass of each species at the start of each year, scaled to a mean of one.

A.7 Run 22 settings

Fleet: 30 vessels, \(U_T\) = 500 sets per year, \(q\) = 0.1, vessel CV 0.2, Tweedie power 1.2 and dispersion 0.1, \(\kappa\) = 1. Space: 50 x 50 cells, \(\sigma_z\) = 150 m, \(\phi\) = 4 cells. The species correspond loosely to petrale sole, Dover sole, darkblotched rockfish and sablefish. Thorson et al. (2017) used different parameter values (their Table B1).

A.8 Differences from the original code

simulate_bathymetric() keeps the structure and settings of Run 22. The following changes were made when turning the scripts into a package function:

• Growth rate by species. The original code computed \(r \times d\) with \(r\) a vector over species and \(d\) a cells x species matrix, so R recycled \(r\) over cells rather than species. \(r_p\) is now applied by species.
• Stable initial distribution. The original iterated the movement process 1000 times from a uniform distribution, which may not converge for short movement distances. Because the distance kernel is symmetric, the stationary distribution has the closed form \(\pi_p(s) \propto h_p(s) \sum_{s'} K_p(s,s')\,h_p(s')\) (detailed balance), which is now used. For the same reason the movement matrix is never formed explicitly, which makes the model much faster.
• Reporting blocks. Reported positions are the centres of 5 x 5 blocks aligned with the grid; the original rounding produced 11 uneven blocks per axis.
• Number of species. The initial biomass noise assumed exactly four species; it now uses \(P\).
• Area closures. In closed cells catchability is zero, and open cells keep the mean catchability (the original set it to one).
• Software. The archived RandomFields package is replaced by a circulant-embedding simulator with the same Gaussian covariance, and the Tweedie generator is implemented directly (compound Poisson-gamma). Results are therefore not identical number by number to the original runs.

3. Appendix B: extension to a longline tuna fishery (simulate_tuna())

The tuna model keeps the population, movement, effort and catch structure of Appendix A and adds what matters for longline targeting: a pelagic horizontal habitat, a vertical distribution of each species, the depth of the hooks set with a given number of hooks between floats (HBF), and a fleet that chooses both where to fish and how deep. The aim is a conceptual, easily tuned model; all values are illustrative.

B.1 Species and horizontal habitat

Four species live in a conceptual Southern Hemisphere ocean (north = towards the equator = warm): yellowfin (YFT) and bigeye (BET) prefer warm water; southern bluefin (BFT) and albacore (ALB) prefer cooler water. Bathymetry is replaced by sea surface temperature,

with mean \(\overline{\mathrm{SST}}\) = 23 °C, a south-north gradient \(g\) = 10 °C and a Gaussian random field \(\eta\) (SD 1 °C, scale 6 cells). Habitat suitability is a normal preference, \(h_p(s) = \mathrm{N}\big(\mathrm{SST}(s);\, \mathrm{SST}_p^{opt},\, \tau_p\big)\), used in the movement model of A.2 in place of the depth and range terms, with mean movement distance \(\lambda_p\) = 1.5 cells. Wide tolerances \(\tau_p\) keep the species mixed in space. Carrying capacity \(B_{0,p}\varepsilon_p\) is distributed at the stable distribution, and stocks start at a fraction \(\delta_p\) of it (B_start; 1 = unfished). By default \(\delta\) = 0.6, 0.7, 0.2 and 0.9 (YFT, BET, BFT, ALB), with \(q\) = 3. Effort follows the logistic increase of A.3 towards 800 sets per year but is cut by 25% in years 21-30 (effort can be one value per year). The true trends then differ between species: ALB is fished down to about half while targeted, BET declines steadily, BFT recovers from 20% of carrying capacity and declines again once targeted, and YFT dips and recovers.

B.2 Vertical distribution and the thermocline

The water column (0-500 m) is divided into 10-m bins \(j\) with mid-depths \(z_j\). In cell \(s\) and year \(t\), the vertical distribution of species \(p\) is a discretised normal profile,

with mean depth \(\mu_p\) and SD \(\sigma_p\). Within each pair one species is shallow (YFT, ALB) and one deep (BET, BFT). The thermocline anomaly \(a(s,t)\) shifts all profiles up (shallow thermocline) or down (deep thermocline):

where \(a_0\) and \(a_t\) are Gaussian random fields with SDs thermocline_sd and thermocline_annual_field_sd, and \(\delta(t)\) is a stationary AR(1) process with SD thermocline_annual_sd and autocorrelation \(\rho\) = 0.5. By default only the fixed pattern is used; values of about 15 m and 10 m for the annual terms give a realistic ENSO-like range of year-to-year variation.

B.3 Hook depth

Each HBF level \(h\) (5, 10, 15, 20) has a typical median hook depth \(D_h\) (110, 160, 205 and 250 m), roughly following the aggregate hook-depth distributions of the Japanese longline fleet in Fig. 2 of Bigelow et al. (2002). Rather than modelling the mechanics of the line (catenary geometry, line angle, currents), hook depth is summarised by two lognormal errors. The median hook depth of set \(u\) is

so that HBF is an imperfect proxy for fishing depth (set_cv = 0.2), and the hooks of a set are lognormally distributed around \(D(u)\) with \(\mathrm{sdlog}(\mathrm{CV}_{hook})\) (hook_cv = 0.35), where \(\mathrm{sdlog}(\mathrm{CV}) = \sqrt{\log(1 + \mathrm{CV}^2)}\). This gives the proportion of hooks in each depth bin, \(g_u(z_j)\). The expected profile for an HBF level, used by skippers and by the habitat-based standardisation, is lognormal with median \(D_h\) and the combined spread \(\sqrt{\mathrm{sdlog}(\mathrm{CV}_{set})^2 + \mathrm{sdlog}(\mathrm{CV}_{hook})^2}\).

B.4 Effective catchability

Following the habitat-based approach of Bigelow et al. (2002), the effective effort of a set depends on the overlap between where its hooks fish and where the fish are. The vertical availability of species \(p\) to set \(u\) in cell \(s\) is

scaled so that \(A_p = 1\) corresponds to fish spread evenly over a 100-m layer fished by all hooks (\(\Delta z\) = 10 m). The catchability of the set is \(q_p(u) = \tilde q_p(u)\, A_p(u)\), where \(\tilde q_p(u)\) is the vessel catchability with Tweedie error of A.5 (mean \(q\) = 3). Catch is then calculated as in A.5.

B.5 Targeting: where and how deep

Skippers evaluate the expected value of each cell \(s\) and HBF level \(h\) in year \(t\),

where \(\bar A_p(s,h,t)\) is the availability under the expected hook profile of HBF \(h\) and the current thermocline. Effort is allocated to cells in proportion to \(\big[\max_h V(s,h,t)\big]^{\kappa}\) with \(\kappa\) = 1.5 (spatial targeting). Each set then chooses its HBF with probability

where \(\tau\) (HBF_choice_temp = 0.15) controls how strictly the best gear is chosen (vertical targeting). Prices change in three 10-year blocks: albacore (shallow, cool water), then bigeye (deep, warm water), then bigeye and southern bluefin together (deep water, warm north and cool south). This loosely follows the historical shift from albacore to bigeye targeting in the Japanese longline fleet (Hoyle and Okamoto 2013) and the later interest in southern bluefin as it recovers. The first switch changes both where and how deep the fleet fishes; in the last decade the fleet splits between two deep targets in different areas. Prices of the target species are four times those of the others (tuna_settings("mild"): twice, with \(\kappa\) = 1.25 and \(\tau\) = 0.2).

B.6 Data and truth

For each set and species, Data contains year, set and vessel identifiers, true and reported positions, HBF, SST, the thermocline anomaly and the catch (CPUE). It also keeps latent truth that an analyst would not observe: the median hook depth of the set, the realised vertical availability \(A_p(u)\), and the availability a habitat model would predict from HBF, the expected hook profile and the average (fixed) thermocline (used only by the habitat-based standardisation). The true index is the total biomass at the start of each year.

B.7 Default settings

B.8 What the extension adds, in short

4. Appendix C: standardisation models

The simulated data are standardised with Tweedie generalised additive (mixed) models, following the workflow of the WKBMSYSPICT CPUE guidelines (Winker and Cardinale 2022). This appendix describes the general model, the data processing behind each way of accounting for targeting, and how each model is specified in R. All models are fitted with standardize_cpue(), which uses the mgcv package (Wood 2017).

C.1 General model

Let \(y_{i}\) be the catch of the species of interest in set \(i\), fished in year \(t(i)\) by vessel \(v(i)\) at reported position \(\big(\mathrm{lon}_i, \mathrm{lat}_i\big)\). Catches contain many zeros and are right-skewed, so they are modelled with a Tweedie distribution, \(y_i \sim \mathrm{Tweedie}(\mu_i, p, \phi)\) with \(\mathrm{Var}(y_i) = \phi\,\mu_i^{p}\) and \(1 < p < 2\) (a compound Poisson-gamma distribution with a point mass at zero; Candy 2004; Dunn and Smyth 2005). The mean is

where \(\beta_t\) are the year effects (year as a factor), \(S(\cdot)\) is a smooth spatial surface, \(T_i\) is one of the targeting terms of C.3 and \(\alpha_v\) is an optional random vessel effect (Helser et al. 2004). The power \(p\) is estimated together with the other parameters (family = tw() in mgcv), rather than in a separate step as in the guidelines.

C.2 The spatial term \(S\)

\(S(\mathrm{lon}, \mathrm{lat})\) is a two-dimensional thin-plate regression spline of the reported positions. In the simulator, positions are reported as the centres of \(5 \times 5\) blocks of cells (\(X_s, Y_s\)), so the spatial term is s(Xs, Ys, k = ...); with real data it would be, for example, s(lon, lat, k = ...) on the reported 1° or 5° squares. The basis dimension \(k\) controls the maximum flexibility; by default it grows with the number of reported positions (8 for 100 blocks, at most 60). In model names, the suffix +S adds this term (e.g. DPC+S). Settings: spatial_res (block size in cells; 1 = exact position) and spatial_k.

C.3 Targeting terms

Catch composition (used by Cluster and DPC). For each set with a positive total catch, the catches of all \(P\) species are converted to proportions by weight, which removes the effect of catch size, and fourth-root transformed to down-weight dominant species:

Sets with zero catch of all species are excluded from the analysis.

Cluster (discrete fishing tactics). The composition matrix \(\mathbf{W}\) is grouped into \(G\) clusters with k-means (20 random starts), and the cluster of each set is used as a categorical "fishing tactic" (FT) (He et al. 1997):

Following Winker et al. (2014), the number of clusters is the number of principal components used by DPC plus one, \(G = n + 1\) (by default 4; setting n_clusters). The guidelines used CLARA clustering; here k-means is used. Clustering assumes that discrete targeting strategies exist, as recommended by Hoyle et al. (2014).

DPC (continuous targeting, direct principal components). The composition matrix \(\mathbf{W}\) is centred, scaled and decomposed by principal component analysis. The scores of the first \(n\) components, \(\mathrm{PC}_{1i}, \dots, \mathrm{PC}_{ni}\), are used as continuous predictors with smooth effects (Winker et al. 2013, 2014):

where \(f_k\) are thin-plate splines with low flexibility (s(PCk, k = 3)). Because the PCs are decorrelated linear combinations of the composition, DPC treats targeting as a continuum rather than as discrete groups. The number of components \(n\) is set with n_pc (default 3; the guidelines used 2). Winker et al. (2014) select the non-trivial components with the optimal coordinates test (a non-graphical solution to Cattell's scree test, nFactors::nScree()), because model selection criteria such as AIC tend to retain too many components. With \(P\) species there are at most \(P - 1\) informative components; using more components captures more of the targeting signal but also increases the risk of absorbing real changes in relative abundance (C.6).

HBF (recorded gear covariate). The number of hooks between floats is used as a factor, one catchability multiplier per HBF level (Hoyle et al. 2014):

It needs no catch-composition information, but assumes that each HBF level has the same effect on catchability in every year and area.

Combinations. Cluster_HBF and DPC_HBF add HBF to the cluster or DPC terms, \(T_i = \gamma_{\mathrm{FT}(i)} + \gamma_{\mathrm{HBF}(i)}\) or \(T_i = \sum_k f_k(\mathrm{PC}_{ki}) + \gamma_{\mathrm{HBF}(i)}\), combining the recorded gear with the realised catch mix.

HBS (optional, habitat-based). Habitat-based standardisation (Bigelow et al. 2002) estimates no gear effect; instead the predicted vertical availability \(\hat A_{pi}\) (B.4) enters as an offset, i.e. effort is converted to effective effort: \(T_i = \log \hat A_{pi}\). It requires a habitat model and ocean data and is therefore not part of the routine comparison.

C.4 Random vessel effect

Differences in fishing power among vessels can be included as a random intercept, \(\alpha_v \sim \mathrm{N}(0, \sigma_\alpha^2)\). In mgcv this is a smooth with a random-effect basis, s(vessel, bs = "re"), estimated by REML together with the other smooths (option vessel_re = TRUE, which adds the term to every model). When predicting the index, the vessel term is excluded, i.e. the index refers to an average vessel. The guidelines fitted the same model with mgcv::gamm() and random = list(vessel = ~1), which gives equivalent estimates but requires the Tweedie power to be fixed in advance.

C.5 Model formulae in R

All models share the structure y ~ year + [S] + [targeting] + [vessel], with year and vessel as factors:

In the simulator lon, lat are the reported block centres Xs, Ys. A single model fitted directly with mgcv:

library(mgcv)
fit <- gam(y ~ year + s(lon, lat, k = 8) + s(PC1, k = 3) + s(PC2, k = 3) +
             s(PC3, k = 3) + s(vessel, bs = "re"),
           family = tw(), method = "REML", data = sets)
# faster for large data sets:
fit <- bam(y ~ year + s(lon, lat, k = 8) + HBF + s(vessel, bs = "re"),
           family = tw(), method = "fREML", discrete = TRUE, data = sets)

The same comparison with the package:

fits <- standardize_cpue(
  sim,
  models      = c("Base", "Base+S", "Cluster+S", "DPC+S", "HBF+S"),
  spatial_res = 5,      # reported blocks of 5 x 5 cells
  spatial_k   = NULL,   # basis size of S (NULL = automatic)
  n_pc        = 3,      # principal components for DPC
  n_clusters  = 4,      # fishing tactics for Cluster (default n_pc + 1)
  vessel_re   = FALSE,  # TRUE adds s(vessel, bs = "re")
  engine      = "bam")  # "gam" (REML) or "bam" (fREML, discrete)
fits$formulas           # the formulae actually fitted
# further models, e.g. an HBF x period interaction:
standardize_cpue(sim, models = "HBF+S",
  add_models = c(HBF_period = "year + s(Xs, Ys, k = 8) + HBF_period"))

C.6 From model to index, and how accuracy is measured

The annual index is the prediction for each year with all other terms held at fixed reference values: the median reported position, PC scores of zero, the most common fishing tactic and HBF level, an average vessel (random effect excluded) and no offset. Predictions are made on the log scale, back-transformed and scaled to a mean of one; approximate 95% confidence intervals are \(\exp(\hat\eta_t \pm 1.96\,\mathrm{SE}_t)\), scaled in the same way. The accuracy of an index is the root mean squared error of the log index relative to the true (simulated) biomass index \(B_t\), also scaled to a mean of one:

A caveat for composition-based methods. Cluster and DPC use the same catches on both sides of the model: the catch of the species of interest is the response, and it is also part of the composition used to define tactics or PCs. Within a year, variation in composition among sets mainly reflects targeting, which the PC or cluster effects capture; a real change in the relative abundance of one species, however, also shifts the average composition of that year and can be partly absorbed by the targeting terms, removing part of the true trend. The risk is largest with few species (with two species the composition has a single dimension, so targeting and abundance cannot be separated) and when many components are used. Model selection criteria are therefore not reliable guides to the number of components or clusters (Winker et al. 2014).

References