A simulation framework for evaluating multi-stage sampling designs in populations with spatially structured traits

Virtual population

We simulated a virtual population by resampling the observed data from surveys to preserve the spatio-temporal and statistical features of the natural population. Since natural populations are dynamic and change over time and space in response to the variations in their environment, the virtual population encompasses biological data under different demographic and physical conditions. Thus, we generated six scenarios, which when compiled represent the average physical and demographic conditions of the virtual population. The six scenarios conforming the average virtual population included (1) warm and (2) cold scenarios, based on the average bottom temperature of the middle shelf, since the spatial distribution is highly influenced by water temperature (Stabeno et al., 2001; Ciannelli & Bailey, 2005; Barbeaux, 2017). To assess the studied trait structure, we generated two more scenarios corresponding to (3) high and (4) low age structure diversity. Diversity with respect to age structure was measured using the Simpsons diversity index in the vegan (Oksanen et al., 2018) library of R software (R Core Team, 2017). The population abundance (stage 1) and number of fish measured (stage 2) in the high and low age diversity scenarios were set to the same values to avoid confounding this effect with those of abundance variability. Finally, we incorporated (5) high and (6) low total abundance scenarios into the virtual population. In order to prevent similar confounding effects as aforementioned, the age structures of these two scenarios were constrained to be as similar as possible by keeping less than a 3% difference in the number of individuals at any given age-group.

To create every particular scenario (for instance, high total abundance), the resampling was limited to the observed data that meet the condition of that particular scenario (for instance, data on the 5 years with the highest total abundance). For each scenario, we calculated total abundance (stage 1) and percentage of individuals sub-sampled for trait analysis (stages 2–3) at any given sampling location. Further, we determined a value of the trait of interest for every individual in the virtual population using a function to estimate this trait. The parameters of the trait function can be tuned to the observed data to represent the realized spatial variability of individual traits. Populations spatially structured in abundance or traits display some degree of spatial autocorrelation across locations, which was also included in the virtual population.

The virtual population included, by site, as many fish as those historically measured on average over the years corresponding to a particular scenario (except for the high and low age diversity scenarios, which share the same values as aforementioned). Since the sex of some fish was undetermined in the observed data, we assigned sex to each fish based on the sex ratio calculated for each length interval (one cm) and site. When the sample size was lower than 10 fish, the sex ratio was calculated by a higher level of aggregation such as length interval and stratum, only length interval and finally, overall sex ratio. Further, we determined the age of every fish based on its length and the geographic location where it was caught (i.e., sampling site), using the following Eq. (1): (1) $${\text{Age}}_i=\text{α}+{\text{β}}_1{\text{Size}}_i^2+{\text{β}}_{2}r{\text{λ}}_i+e_i$$ where i corresponds to every sampled fish, α is the intercept, β₁ is the size coefficient, which varies between male (1.05 e⁻⁵) and females (1 e⁻⁵) and β₂ is the geographic coefficient (−0.031) that accounts for size at age variability in rotated latitude (rλ) across the shelf. Due to the NW—SE orientation of the EBS inner-outer shelf axis, geographical coordinates were rotated to accommodate and simplify this orientation to a N—S axis. We added a site-specific spatial auto-correlated error (e_i) to the age values obtained in Eq. (1). Spatial auto-correlation among sampling sites was calculated from a covariogram based on sampling site coordinates and average age residuals from observed data using gstat (Pebesma, 2004) and fields (Nychka et al., 2017) libraries in R software (R Core Team, 2017).

Sampling simulation

We simulated a multi-stage field sampling in every scenario of the virtual population and for each of the (sub-) sampling designs under evaluation. This was equivalent to conducting a field survey in the study area for as many times as the number of scenarios. Sampling of the virtual population (stage 1) occurred in all of the sampling sites for each of the six scenarios, resembling the same protocol used in the scientific surveys (Conner & Lauth, 2017; Puerta et al., 2018). At half of the sampling sites (randomly selected), two different sub-sampling strategies for age data collection (stage 3) were simulated to be compared: the length-stratified strategy (LSS), that collects three fish per length interval (one cm), sex and area (NW or SE of the EBS); and the RS, which randomly collects four fish at each site, regardless of the length, sex or area. Thus, the total sample size of the RS is larger than LSS. These sample sizes represent the minimum effort and were chosen to closely mimic the current sampling protocols in the Pacific cod case study. We also assessed results in which the stage 3 sample size was equal to all fish caught (i.e., no subsampling) and results where LSS and RS ended up in the same sample size. For the latter, and for each scenario, several records in the RS were randomly removed to match the total sample size obtained with the LSS.

Simulation of trait sub-sample processing

Characterization of some important individual traits such as age or maturity relies on statistical estimations from technical proceedings and/or depends on subjective interpretations of the biological conditions since the true value of the trait is usually unknown. Thus, the process of extracting trait information from the collected biological sub-samples introduces a new component of variability that includes bias and error from the trait’s true value.

Ageing processes in fish mainly relies on the interpretation of annual marks (annuli) in the otoliths. Despite otolith readers being expertly trained, some degree of subjectivity is involved (Matta & Kimura, 2012). Kastelle et al. (2017) defined age-specific bias in EBS Pacific cod ageing by comparing ages derived from counting otolith marks and those (considered true ages) derived from stable oxygen isotopes (δ¹⁸O) analysis. Overall, the probability of assigning an age equal to the true age was approximately 61% (Kastelle et al., 2017). In contrast, the error in otolith reading refers to the degree to which an age estimate is reproducible by the same or a different reader (Matta & Kimura, 2012). Usually a second reader tests 20% of all specimens to calculate the age-specific inter-reader coefficient of variation (Kimura & Anderl, 2005; Matta & Kimura, 2012). Thus, in our case study simulation, we differentiate the “true age” of an individual (calculated from Eq. (1) in the virtual population) from the “otolith age,” which is recorded in the age sub-sample and derived from the true age plus the age-specific reading bias and error (Kastelle et al., 2017).

Evaluation of the (sub-) sampling strategy

We followed a (sub-) sampling strategy evaluation based on the comparison of: (i) prediction errors for the trait values derived from a statistical model, with the same covariate as those included in Eq. (1), parameterized on each of the trait sub-samples, respectively; (ii) the trait frequency distribution of the virtual population estimated from each sub-sample; (iii) the mean and modal estimated values of the trait and; (iv) the average spatial patterns of the trait recovered from each sub-sample. This approach ensures a quantitatively valid comparison of the estimates derived from each (sub-) sampling strategy.

1. We used generalized additive models (GAM; Wood, 2006) to fit “otolith age” (as described above) in relation to fish length, sex and geographical location (Puerta et al., 2018), using mgcv library (Wood, 2006) in R software (R Core Team, 2017). Model Eq. (2) was defined as follows, (2) $${\text{Age}}_i=\text{α}+s_1\left({\text{Size}}_i\right)\times {\text{Sex}}_i+s_2\left({\varphi }_i,{\text{λ}}_i\right)+{\text{ε}}_i$$ where i corresponds to every sampled fish. α is an intercept, and s represents the smooth functions for the Length and the geographic location (longitude φ, latitude λ). To reduce circularity, the model Eq. (2) is purposely different from Eq. (1) (used to predict the true age of individuals in the virtual population). Sex is included as a factor, and so, different smooth functions are fitted for males and females; and ε is a normally distributed error term. Since the RS and LSS age sub-samples are by definition not proportional to the catch, these data were weighted by the size-specific abundance value within each site in the model. This ensures that age data from a particular length interval with a higher representation in the catch (as number of fish) have larger weights in the model. Additionally, these weights correct variance structure in multistage data, removing any bias in the precision of model estimates (Zuur et al., 2009). The mean predicted errors of the models (as model “otolith age” prediction minus virtual population “true age”) parameterized on the LSS and RS sub-samples were compared as an indicator of goodness of fit and predictability skills. The mean predicted error was calculated on 300 iterations by removing 30% of the data at a time and predicting age of the deleted cases with a model fitted on the remaining data.

2. Then, we used “otolith age” predictions of the Eq. (2) model parameterized on each age (sub-) sampling strategy sample to estimate the age frequency distribution of the virtual population.

3. Mean and modal size at age obtained from the LSS and RS age sub-samples, respectively, were compared for age-groups 1 to 3. Finally, (iv) spatial patterns in average size and age over the EBS shelf derived from the LSS and RS age sub-samples were compared with those of the virtual population.

Sensitivity test

We performed a sensitivity test to disentangle and quantify (i) the contribution of the spatial term included in the model used in the sub-sampling strategy evaluation (Eq. (2)), (ii) source and magnitude of the errors in the population estimates derived from the inclusion of spatial autocorrelation in the virtual population and (iii) bias and error derived from trait sub-sample processing. The aforementioned evaluation was repeated for a Pacific cod virtual population without spatial autocorrelation and RS and LSS age sub-samples where the “otolith age” equals the “true age,” and every other combination of these errors to compare and quantify error contribution in the population estimates.

Virtual population

The combined Pacific cod virtual population included approximately 83,000 individuals across the six simulated scenarios, with an average of 13,833 individuals simulated per sampling event (i.e., survey year). The spatial distribution (Fig. 1) in abundance averaged across the six scenarios highlights northern and southern regions in the inner shelf as the hotspot in the study area, although the central middle shelf also presented remarkable abundance. Size and age exhibited similar average spatial patterns, increasing progressively from the inner to the outer continental shelf. Larger and older individuals were particularly concentrated in the northwest region of the outer shelf. However, the distributional patterns in abundance, size and age varied considerably across the six individual scenarios (Fig. S1). For example, in warm years cod abundance is higher in the inner and middle shelf of the Bering Sea, and average age and size is less structured over space than in cold years.

Evaluation of the (sub-) sampling strategies

Age data sub-samples obtained by the RS and LSS strategies differed in the number of records collected (stage 3) as a result of their respective intrinsic designs. While the random sub-sampling strategy accounted for ∼1,000 records per scenario including all the preselected sampling sites (stage 1), the LSS collected only ∼800, dismissing ∼5% (13–22) sampling sites per scenario as a result of the targeted number of otoliths having already been reached. As expected due to the differences in the designs, the random sub-sample showed a normal distribution for the length measurements (stage 2), where the most common lengths are more representative (Fig. 2). Naturally, the length-stratified sub-sample led to a uniform distribution where all length intervals are equally represented, with the exception of the extremes that are not commonly found in nature (Fig. 2).

1. Generalized additive models parameterized on the RS and LSS age sub-samples showed minimal differences in terms of explained variability, prediction errors (Table 1) and functional shape in size at age (Fig. 3). These negligible differences were also observed at each individual scenario, but we observed greater standard deviation in predicted errors from the LSS strategy due to the lower sample size used. Model residuals were inspected for normality, homoscedasticity and independence. None of these model residuals exhibited strong spatial autocorrelation, indicating that the model formulation was sufficient to address the spatial autocorrelation incorporated in the virtual population (Fig. S3).

2. Similar age frequency distributions were predicted from the two age sub-samples (Fig. 4). This was expected considering the negligible differences in size-at-age relationships. However, both sub-sampling strategies misrepresent the number of fish at the most representative age groups (ages 2, 3 and 4) and, particularly the age 2 group. This age accounted for ∼26% of the individuals misplaced in an age group regardless of the sub-sampling strategy used. Although shape and prediction of the age frequency distribution varied notably across the different individual scenarios (Fig. S2), the underestimation of the age 2 group was observed in all cases. The number of miss-aged individuals increased by 2% approximately in the individual scenarios comparing with the average.

3. The average point estimates of mean and modal size at age derived from the RS and LSS age sub-samples were both almost identical to the values observed in the virtual population for age-groups 1–3 (Fig. 5A). However, a mismatch between population and sub-sample point estimates was frequent and sometimes remarkable in the individual scenarios, with the RS tending to provide more accurate and precise estimates (Fig. 5). The most diverging differences occurred in cold conditions and high diversity of age structure. Spatial distribution of cod in simulated scenarios also notably diverged between warm and cold years (Fig. S1).

4. The largest and most remarkable differences between the RS and LSS sub-sampling strategies were observed in the recovery of the average spatial patterns in size and age (Fig. 6). The length-stratified sample was not able to capture the spatial patterns in the virtual population and largely overestimated these patterns in shape and particularly, in magnitude. Larger and older fish estimated by the LSS were disproportionally widely distributed in the northwest region of the outer shelf, presenting an average age of 3 years older than the age observed in the virtual population. By contrast, the RS strategy recovered very similar patterns as those in the virtual population. Despite the fact that the population distributional patterns considerably varied between the different individual scenarios, the misrepresentation described for the LSS was observed in all cases (Fig. S4).

Sensitivity test

We analyzed models parameterized in the RS and LSS, respectively, as in Eq. (2) but removing the spatial term. Results indicated that spatial location had a small contribution in the variance explained and modestly reduced the predictability skills (Table 1). However, this small reduction turned into large differences when comparing the predicted age distribution (number of fish per age group) with and without including the spatial term in the model (Eq. (2)). The effects of incorporating different errors, such as autocorrelation and reading bias and error were mainly detected in the predicted age distribution. Autocorrelation errors in the virtual population did not affect the predictive ability of the models, since as mentioned, it was properly addressed by the spatial term in the model formulation. However, the autocorrelation error contributed to ∼13% of the fish being mis-aged (Fig. 7). Removing the reading bias and error resulted in the mis-aged number of fish to approximately be 5% (with a larger reduction in RS than in LSS, Fig. 7).

Unexpectedly, this variation in the predicted number of fish at a given age did not have a notable effect on the point estimates for mean and modal size at age (Fig. 8). Nevertheless, a mismatch in the point estimates between population and age sub-sample values can be an issue associated with small sample size in the RS and LSS sub-samples. Indeed, when we increased the sub-sample size to all fish caught (stage 3, i.e., not random or length-stratified sub-sampling either), the points estimates match perfectly between the population and the age sub-sample data (results not shown). The sampling strategy evaluation provide the same results when considering LSS and RS subsamples with the same total sample size (results not shown).