In PaulRegular/MSP: Multispecies production model incorporating covariance

Methods {-}

Model formulation {-}

Trends in fish populations have frequently been described using state-space production models of the form

\begin{align} B_y =& \left ( B_{y-1} + g(B_{y - 1}) - C_{y - 1} \right ) e^{\delta_y}, & (#eq:pro) \ \ I_{y, i} =& q_i B_y e^{\varepsilon_{y, i}}, & (#eq:obs) \end{align}

where $B_y$ is biomass at the start of year $y$, $C_y$ is the catch through year $y$, $g(B)$ is production as a function of biomass, $I_{y,s}$ is an index of relative abundance in year $y$ from survey $i$, $q_i$ is the time-invariant catchability coefficient for survey index $i$, $\delta$ is process error, and $\varepsilon$ is observation error. Statistical challenges aside, the most difficult aspect of this model to parameterize is the production function as it needs to capture changes caused by growth, recruitment, and natural mortality. Schaefer [-@schaefer1954] proposed a solution by applying the logistic equation to describe self-limiting growth,

\begin{align} g(B) = rB \left( 1 - \frac{B}{K} \right), & (#eq:logistic) \end{align}

where $r$ is the maximum per-capita rate of change and $K$ is the carrying capacity. That is, a populations' intrinsic ability to grow ($rB$) is limited by the size of the current population relative to the maximum biomass the system can support ($1 - B/K$). While this formulation offers an elegant description of single-species population dynamics, it assumes that density-dependent effects are solely caused by intraspecific competition and ignores the potential effects of other species inhabiting the same ecological area, competing for the same resources. We present an extension of equation \@ref(eq:logistic) that attempts to account for intra and interspecific competition by assuming that density-dependent effects are incurred when the total biomass of multiple species, represented by $s$, exceeds the capacity of the system, $K_{\Sigma}$,

\begin{align} g(B_s) = r_s B_s \left (1 - \frac{\sum_{s} B_s}{K_{\Sigma}} \right ). & (#eq:ms-logistic) \end{align}

While intrinsic rates of growth may vary across species, this formulation implies that the growth of all species is ultimately limited by the finite amount of energy in a region (i.e., as the total population of all species in the system increases towards $K_{\Sigma}$, year-over-year growth of all species slows). Combining equations \@ref(eq:pro), \@ref(eq:obs), and \@ref(eq:ms-logistic), our model becomes

\begin{align} B_{y, s} =& \left ( B_{y-1, s} + r_s B_{y-1, s} \left (1 - \frac{\sum_{s} B_{y-1, s}}{K_{\Sigma}} \right ) - C_{y - 1, s} \right ) e^{\delta_{y, s}}, & (#eq:ms-pro) \ \ I_{y, i, s} =& q_{i, s} B_{y, s} e^{\varepsilon_{y, i, s}}. & (#eq:ms-obs) \end{align}

The inclusion of multiple species in the model permits the estimation of covariance. While covarying changes may be apparent in the observations, we assume that most covariance stems from ecosystem processes. We therefore treat observation errors as independent and normally distributed deviations such that $\varepsilon_{y, i, s} \sim N(0, \tau_{i, s})$, where the standard deviation parameter $\tau_{i, s}$ represents species and survey specific levels of observation error. A more flexible error structure is used to describe the process errors to account for temporal dependencies driven by as ecological processes. For instance, species interactions may drive positive or negative population responses resulting from direct or indirect associations. Deviations from expected production may also display temporal dependence if the factors contributing to the process errors change gradually over time. Such inertia may cause positive or negative process errors to persist across several years. A first-order autoregressive (AR1) process was therefore applied to account for temporal dependence. Both sources of dependence are modeled using a multivariate AR1 process where

\begin{align} \delta_{1, s} \sim& N(0, \Sigma_s) & \ \delta_{2, s} =& \phi \delta_{1} + \sigma_s \varepsilon_{2}, \quad \varepsilon_{2} \sim N(0, \Sigma_s) & (#eq:ar1) \ \delta_{y, s} =& \phi \delta_{y - 1} + \sigma_s \varepsilon_{y}, \quad \varepsilon_{y} \sim N(0, \Sigma_s) & \end{align}

and

\begin{align} \Sigma_s = \begin{bmatrix} \sigma^2_{1} & \sigma_{1} \sigma_{2} \rho_{1,2} & \dots & \sigma_{1} \sigma_{s} \rho_{1,s} & \ \sigma_{2} \sigma_{1} \rho_{2,1} & \sigma^2_{2} & \dots & \sigma_{2} \sigma_{s} \rho_{2,s} & \ \vdots & \vdots & \ddots & \vdots & \ \sigma_{s} \sigma_{1} \rho_{s,1} & \sigma_{s} \sigma_{2} \rho_{s,2} & \dots & \sigma^2_{s} & \end{bmatrix}. & (#eq:sigma) \end{align}

The degree of temporal correlation is controlled by $\phi$, where values between 0 to 1 represent low to high correlation, and species-to-species correlations are described by $\rho_{s, s}$, where values between -1 to 1 represent negative to positive correlation. This is a flexible structure that allows for the testing of alternate hypotheses that process errors are independent through time or across species (i.e, $\phi = 0$ or $\rho_{s, s} = 0$). The possibility that process errors are similarly correlated across all species may also be tested by estimating only one $\rho$ parameter. Finally, the magnitude of the process error deviations is controlled by the species-specific standard deviation parameters, $\sigma_s$.

Minor extensions of the formulation also permit the fitting of covariates which may describe an underlying linear effect. Two options were implemented, one that affects the process errors by substituting $e^{\delta{y,s}}$ in equation \@ref(eq:ms-pro) with $e^{\beta_{\delta} X^{\delta}{y,s} + \delta{y, s}}$, and another that affects the carrying capacity by substituting $K_{\Sigma}$ in equation \@ref(eq:ms-pro) with $K_{\Sigma} e^{\beta_K X^K_{y}}$, where $\beta$ parameters capture linear effects of covariates included in design matrices, $X$. The idea is that some factors may affect positive or negative changes in the populations while others may affect change in the total carrying capacity of the system. A covariate option for intrinsic rates of increase was not implemented since one goal of this model is to obtain estimates of this vital rate, which is not expected to change rapidly as it is shaped by natural selection [@hutchings2021]. The formulation was also modified to fit the single-species Schaefer production function by dropping the summation of biomass in equation \@ref(eq:ms-pro) and estimating species-specific carrying capacities, $K_s$, rather than a system level carrying capacity, $K_{\Sigma}$ (i.e., apply equation \@ref(eq:logistic) indexed by species).

Statistical framework {-}

This model was implemented using template model builder [TMB\; @kristensen2015], which is a R [@R] package that enables the fitting of complex nonlinear random effects such as the latent $B$ variable in state-space production models (equation \@ref(eq:pro)). Such variables are not directly measured but are inferred indirectly via observed values. Data fitting is accomplished using a combination of Laplace approximation and automatic differentiation to evaluate the joint likelihood [@kristensen2015]. Like the production model described by Pedersen et al. [-@pedersen2017], both frequentist and Bayesian inference of model parameters are possible. In development, we found that estimation of parameters was generally more successful when vaguely informative priors were specified, and in some cases, not estimatable when unconstrained.

Case study {-}

Data {-}

The multispecies production model described above requires two basic inputs for each species: a time-series of catch ($C_{y,s}$ in equation \@ref(eq:ms-pro)), and an index of population size ($I_{y,s}$ in equation \@ref(eq:ms-obs)). The Northwest Atlantic Fisheries Organization (NAFO) and Fisheries and Oceans Canada (DFO) have been collecting and curating such information for multiple fish populations along the shelves of Newfoundland and Labrador (NL) since the 1970s. The communities inhabiting these shelves can be divided into several regions with distinct productivity [i.e. ecosystem production units; @pepin2014]. For our case study, we tallied catch data and calculated survey indices of multiple demersal fish populations from three regions (Figure \@ref(fig:epu-map)): 1) the Northeast NL Shelf (NAFO divisions 2J3K), 2) the Grand Bank (NAFO divisions 3LNO), and 3) Southern NL (NAFO sub-division 3Ps).

max_area <- aggregate(multispic::index, survey_area ~ region, FUN = "max")
tot_area <- sum(max_area$survey_area) * 3.43 # convert from sq nautical miles to sq km
tot_area_text <- format(signif(tot_area, 1), big.mark = ",", scientific = 6) |> 
    paste("km^2^")

Catch data were extracted from NAFO's STATLANT 21A database (https://www.nafo.int/Data/STATLANT-21A, accessed 2022-01-21) and aggregated by region, species, and year. Survey indices were derived from the standardized, stratified random bottom-trawl surveys conducted each spring and fall by DFO; this is perhaps the largest fisheries-independent survey conducted in the world, which aims to cover more than 500,000 km^2^ annually (roughly the size of Sweden or the Yukon, Canada) to depth up to 1500 m. Since the inception of this program in 1971, survey platforms and protocols have undergone a series of changes that affect the continuity of the data collected in each region and season. A Yankee then Engel otter trawl, with nets designed to catch large demersal fish, were used between 1971 to 1994. Starting in the fall of 1995 survey gear was changed to a Campelen shrimp trawl with a small mesh codend, which allowed a broader range of species and size groups to be captured [@chadwick2007]. Within each era of the survey (Yankee, Engel, or Campelen) and for each season and region, samples used in this study were limited to strata that were covered most years (> 80%) and to species found across more than 10% of these core strata. This often resulted in the exclusion of strata >750 m as these areas have been inconsistently covered by the survey. Stratified analyses [@smith1981] were then conducted on the remaining species to obtain indices of total biomass. To minimize bias introduced by inconsistent survey coverage, indices from years where more than 20% of the biomass was likely missed, inferred from time averaged percent occupancy within strata, were excluded from our analysis [sensu NAFO guidelines; page 10, @nafo2019]. Finally, species were ranked by cumulative commercial catch and limited to the seven most commonly caught species, or species group when catch was not consistently distinguished by species, within each region. On the Northeast NL Shelf, the included species were Redfish spp. (Sebastes fasciatus and Sebastes mentella combined), Wolffish spp. (Anarhichas lupus and Anarhichas minor combined), Witch Flounder (Glyptocephalus cynoglossus), American Plaice (Hippoglossoides platessoides), Greenland Halibut (Reinhardtius hippoglossoides), Atlantic Cod (Gadus morhua), and Skate spp. (Amblyraja radiata and Malacoraja senta combined). On the Grand Bank, the included species were Redfish spp., Yellowtail Flounder (Myzopsetta ferruginea), American Plaice, Greenland Halibut, Haddock (Melanogrammus aeglefinus), and Atlantic Cod. Finally, along Southern NL, the species included were Redfish spp., Witch Flounder, American Plaice, White hake (Urophycis tenuis), Haddock, Atlantic Cod, and Skate spp.

Priors {-}

For simplicity, all priors were normally distributed and, in most cases, upper and lower inflection points ($\mu \pm \sigma$; ~68% of the total area under the curve) of each normal prior were defined using values in log or logit space. Upper and lower values were based on previous research or knowledge to impose fairly generic and vaguely informative priors.

Intrinsic growth rate, $r_s$, and carrying capacity, $K$ {-}

To capture a broad range of intrinsic growth rates [@thorson2020], 0.01 and 1 were chosen as the lower, $r_l$, and upper, $r_u$, values; this translates to a normal prior with a mean of -1.15 and a standard deviation of 1.15 on the log-scale. The log-scale $K_{\Sigma}$ prior was informed by total levels of catch and $r_l$ and $r_u$,

\begin{align} K_{\Sigma, l} = \frac{\max(\Sigma_s C_{y,s})}{r_u}, K_{\Sigma, u} = \frac{4 \max(\Sigma_s C_{y,s})}{r_l}, & (#eq:k-prior) \end{align}

where $K_{\Sigma, l}$ and $K_{\Sigma, u}$ are the lower and upper values. On the lower end, this prior imposes the assumption that the fishery is unlikely to have caught the equivalent of $K_{\Sigma}$ and, on the upper end, it assumes that the maximum observed catch represents a portion of $K_{\Sigma}$ [sensu @froese2017]. The division by the lower and upper values for $r$ also accounts for the potential range of productivity. Note that the maximum time-series catches are species specific when estimating species-specific carrying capacities.

Starting biomass, $B_{1,s}$ {-}

Like $K_{\Sigma}$, catch was used to constrain the plausible range of the biomass at the beginning of the time series, $B_{1,s}$, which we will denote as $B0$ to simplify notation. Specifically,

\begin{align} B0_l = \frac{C_{1,s}}{r_u}, B0_u = \frac{4 C_{1,s}}{r_l}, & (#eq:b0-prior) \end{align}

where $B0_l$ and $B0_u$ are the lower and upper values that are log transformed to define the normal prior for $\log(B0)$. Again, adjusting for the potential range of productivity from $r_l$ to $r_u$, this assumes that the fishery did not catch all of the biomass in the first year and it assumes that the catch represents a portion of the biomass in the first year. A potential flaw with the upper value for this prior is that a lack of market demand may contradict the assumption that landings are coarsely proportional to stock size. However, the upper range chosen was considered reasonable for this case study as there was an active fishery for each species examined here at the start of the time series.

Process error variance, $\sigma_s^2$, temporal correlation, $\phi$, and species-to-species correlation, $\rho_{s, s}$ {-}

Considerable process variability may arise from variable recruitment, natural mortality, and/or growth. For instance, species in the Scorpaenidae family, which includes the Sebastes genus, have notoriously variable recruitment which frequently results in "spasmodic" stock dynamics [@cadigan2022; @licandeo2020]. Also, there is evidence that heightened and variable natural mortality contributed to the collapse and slow recovery of Atlantic cod along the Northeast NL Shelf [@cadigan2015; @regular2022]. To account for a wide range of possible standard deviations, $\sigma_s$, a vague prior with 0.01 and 1 as the lower and upper values were chosen. In log-space, this translates to a normal distribution with a mean of -2.3 and sd of 2.3.

There is little knowledge to determine the potential level of temporal or species-to-species dependence present in the focal systems; vague priors were therefore defined for $\phi$ and $\rho_{s,s}$. For $\phi$, 0.1 and 0.9 were logit transformed [$\log(x / (1 - x))$] to define the upper and lower inflection points for a normal prior for the logit of $\phi$. For $\rho$, the logit transformation was shifted [$\log((1 + x) / (1 - x))$] to capture negative and positive lower (-0.9) and upper (0.9) values. When these parameters are estimated, these priors give most credence to moderate temporal correlation (0.5) and no species-to-species correlation (0) but still allow the possibility of high levels of correlation (0.9).

Catchability, $q_{i, s}$, and observation error variance, $\tau_{i, s}^2$ {-}

The catchability of the spring and fall surveys conducted by DFO likely changed over time given the shift in gear from a Yankee to Engel to Campelen trawl as well as spatial shifts in survey coverage of the strata in each region. Survey catchability $q_{i, s}$ is therefore indexed by season and gear, $i$, and species, $s$. Moreover, the lower inflection points for the $q_{i, s}$ priors were informed by the average survey coverage by gear and season. Survey coverage within each region was computed by dividing the average spatial coverage of strata across years by the total area of all strata (i.e., average area covered / area of survey domain). Survey coverage was then multiplied by 0.2 for deep-water species (Greenland Halibut, Atlantic Halibut, Witch Flounder, Redfish spp., White Hake, Silver Hake, and Monkfish) and 0.5 for the remainder. The lower range was widened, especially for deep-water species, to account for potential gear selectivity issues [e.g., escapement under the footgear\; @walsh1992] and availability issues (e.g., portion of the stock in deeper water than covered by the survey). The upper inflection point was set to 1 as it is possible that the survey indices represent overestimates of the true population size in some instances.

A prior for observation error variance was informed by unbiased design-based estimates of survey variance associated with annual biomass estimates [@smith1981]. These estimates were used to calculate the coefficient of variation (CV) for each survey, year, and species; specifically,

\begin{align} CV_{i,s,y} = \frac { \sqrt{Var(I_{i,s,y})} }{ I_{i,s,y} } & (#eq:cv), \end{align}

where $I$ represents the design-based indices of biomass used in equation \@ref(eq:obs). These CVs were log transformed and survey, $i$, and species, $s$, specific means and standard deviations were used to inform the prior for $\tau_{i, s}^2$ in log-space. Rather than treat the estimates of survey CV as a perfect indicator of total observation error, the prior was widened by multiplying the standard deviation of log CV by two to account for observation variance introduced by potential distributional shifts outside of the survey domain. Another two times multiplier was applied for deep-water species as they have more scope for shifting outside the survey domain.

Model selection {-}

The degree to which density-dependence of the demersal fish community inhabiting the NL shelves is driven by intra versus interspecific competition is unknown. Nor is the degree to which species interactions affect population dynamics. It is, however, well known that the community has been fished for more than 500 years and this history was punctuated by a collapse of several stocks, most notably cod, in the early 1990s [@lear1998]. This collapse was not isolated to commercial stocks, as the biomass of several non-commercial species collapsed at the same time. This was followed by a reorganization of the community, which implies that the system experienced a regime shift [@pedersen2017]. Given this context, and the structure of the multispecies surplus production model described above, a series of hypotheses were tested:

Full: For this hypothesis, it is posited that population dynamics within regions are governed by an overarching system-level carrying capacity that affects all focal species as the aggregate biomass approaches system limits (i.e., inter-specific density-dependent effects are assumed). For this hypothesis and all others, annual reported landings of each species are accounted for in the production equation \@ref(eq:ms-pro). Residual variations not explained by intrinsic growth, density-dependent effects or landings are described by process errors that are 1) correlated across species and said correlations are assumed to be unstructured, meaning relationships can be of differing strengths — positive, neutral, or negative — for each species-to-species pair; and 2) assumed to be temporally correlated, following an AR1 structure. Finally, a shift covariate was applied to the carrying capacity parameter to enable the estimation of different system limits before and after the community-wide collapse to assess support for a regime shift in each systems' capacity for the focal demersal species.

Just shift: Model structure is the same as the full model, except the process errors are assumed to be independent and identically distributed random variables (iid) to assess the hypothesis that temporal and species correlations are explained by the shift.

Just correlation: Model structure is the same as the full model, except population dynamics are assumed to be affected by a common and time-invariant carrying capacity. The shift covariate was not applied.

Shared correlation: Model structure is the same as the just correlation model, except species by species correlations in the process errors are assumed to be the same across all pairs. This structure implies that there is a common but unknown environmental variable affecting the population dynamics of all species. The shift covariate was not applied.

Just species correlation: Model structure is the same as the shared correlation model, except the process errors at each time step are assumed to be iid. This structure implies that there is a common but unknown environmental process affecting all species, but the process is noisy with no temporal dependence. The shift covariate was not applied.

Just temporal correlation: Model structure is the same as the shared correlation model, except species by species correlations are assumed to be iid. This structure implies that environmental processes affect each species differently, however, there may be carry-over effects from one year to the next. The shift covariate was not applied.

No correlation: Model structure is the same as the shared correlation model, except both temporal and species correlations are assumed to be iid. That is, population dynamics are thought to be affected by a common time-invariant carrying capacity and residual variations not explained by intrinsic growth, inter-specific density-dependent effects or landings are noisy and independent across time and species. The shift covariate was not applied.

Single-species: Model structure is the same as the no correlation model, except population dynamics are assumed to be governed by species-specific carrying capacities (i.e., intra-specific density-dependent effects are assumed). This formulation is analogous to standard state-space Schaefer production models. The shift covariate was not applied.

The predictive ability of each of these models was tested using two cross-validation approaches: 1) leave-one-out cross-validation (LOO-CV), and 2) hindcast cross-validation (Hindcast-CV). LOO-CV is a form of exhaustive cross-validation where the model is repeatedly conditioned on a training set missing one observation (i.e., one fold) until the number of model folds equal the number of observations in the data. The missing observations are predicted at each fold, permitting assessments of the models' ability to predict the actual value that was left-out at each fold. The hindcast-CV approach is similar, however it focuses on the models' ability to predict the future. Under this approach, the model is repeatedly conditioned on a training set missing observations from the terminal year such that each fold excludes an increasing number of years worth of data from the tail of the time series [@kell2016]. We folded back 20 years and, for each fold, predicted survey indices were compared to the observed survey indices (e.g., observed indices from 2020 were compared with predicted survey indices for 2020 from the model conditioned on data from 1978-2019). For both approaches, we denote predicted survey indices as $\hat{I_j}$ and the left out observations $I_j$, where $j$ represents all unique combinations of years, species, and survey indices present in the left out data. LOO-CV and Hindcast-CV prediction error scores for each model for each region were calculated by taking the mean squared error,

\begin{align} \text{CV Score} = \frac{1}{n} \sum_{j=1}^{n} \left ( log(I_j) - log(\hat{I}_j) \right ) ^ 2 & (#eq:cv). \end{align}

These scores were also averaged across methods (LOO-CV, Hindcast-CV) and region (Northeast NL Shelf, Grand Bank, Southern NL) to obtain an overall score of predictive ability of each model. A well-fitting model will result in predicted values that are close to the excluded values and, therefore, result in lower scores. Ultimately, these scores enable a model-comparison approach to hypothesis testing. These scores are comparable even for the single-species model as it compares species-specific data points to predictions.

Code and Data Availability {-}

All code used for the analyses in this study is publicly available in a research compendium on GitHub: https://github.com/PaulRegular/multispic. The compendium is organized as an R package, enabling straightforward replication of the results presented here and allowing the same methods to be applied to datasets from other regions.

PaulRegular/MSP documentation built on March 1, 2025, 1:49 p.m.