Population estimates for simulation were generated by an operating model (all estimation was completed via the R package Template Model Builder, TMB, Kristensen et al., 2016) and developed in the LIME package (Rudd and Thorson, 2018;¹) in the statistical software R (R Development Core Team, 2019). Rudd and Thorson (2018, Table 2 and Table 3 in their manuscript) described the population dynamic equations used in the operating model and functions for generating data. For natural mortality of P. argenteus and S. taty in the Bohai and Yellow Seas, the respective values (Table 1) were assumed to be constant so that in each instance the median value of natural mortality was chosen from estimates made using multiple methods (the Barefoot Ecologist’s Toolbox hosts a tool; Cope, 2017). For other biological and fisheries parameters of P. argenteus and S. taty, values were taken from local studies (Zhao, 1987; Liu et al., 1990; Tang and Ye, 1990; Chen, 1991; Jin et al., 2005a; Chen et al., 2018; Xu et al., 2019). Life history parameters for growth and length-at maturity are also listed in Table 1.

Levels of initial biomass depletion compared to carrying capacity of each simulated population were drawn from a uniform distribution between 0.35 and 0.95, and three scenarios involving combinations of fishing mortality and recruitment variability were explored (Figure 3):

Scenario (2) and (3) included a standard deviation for fishing mortality equal to 0.2 (Rudd and Thorson, 2018), a standard deviation of recruitment residuals equal to 0.737 and a first-order autoregressive coefficient equal to 0.426, the mean of the predictive distribution from a meta-analysis of recruitment variability in global fish orders (Thorson et al., 2014). Selectivity at length was assumed to follow a two-parameter logistic model (eq. 4 in the Rudd and Thorson, 2018), with an estimated parameter length at 50% selectivity and a second parameter representing the difference between length at 95 and 50% selectivity. We repeated this process running 100 iterations or replicates for each combination for each population to generate the simulation data. More information can be found on Github (see text footnote 1) and in Rudd and Thorson (2018).

A sample size of 250 individuals was assumed, which approximates the more realistic sample sizes of P. argenteus and S. taty measured annually in the Bohai Sea. Although our ensemble models were only fitted to estimated values obtained from individual length-based models with length data as the sole input, the performance of the LIME model with rich data assumptions was nevertheless tested. This helped to gauge the level of results that could be deemed acceptable when individual length-based models were fitted solely to length data. Therefore, a total of 13 data-availability scenarios for LIME were set up: (1): “Rich20” with 20 years of catch, abundance index, and length data; (2) “Rich10” with 10 years of total catch, abundance index, and length data; (3) “Rich5” with 5 years of total catch, abundance index, and length data; (4) “I_LC20” with 20 years of abundance index, and length data; (5) “I_LC10” with 10 years of abundance index, and length data; (6) “I_LC5” with 5 years of abundance index, and length data; (7) “C_LC20” with 20 years of total catch, and length data; (8) “C_LC10” with 10 s of total catch, and length data; (9) “C_LC5” with 5 years of total catch, and length data; (10) “LC20” with 20 years of length data; (11) “LC10” with 10 years of length data; (12) “LC5” with 5 years of length data; and (13) “LC1” with 1 years of length data. Only one data availability scenario was set up for LB-SPR: “LBSPR” with 1 year of length data.

Two individual data-limited length-based models that use simulated data and basic life-history and/or fishery parameters to estimate SPR were fitted. These length-based models were chosen because they can be fitted to a majority of fisheries around the world, are well established in the literature, and have been extensively simulation-tested (Rudd and Thorson, 2018; Chong et al., 2019; Pons et al., 2019, 2020; Halim et al., 2020).

Length-based spawning potential ratio is one of the prominent length-based methods for estimating reference points in data-limited fisheries and enables a rapid assessment of stock status relative to unfished levels. This model assumes equilibrium conditions (invariant recruitment and mortality) by using a static, equilibrium-based relative age-structured model (Hordyk et al., 2015; Prince et al., 2015). Since the status determination of the LB-SPR model is based on 1 year of data at a time, we only used the length data of the final year to calculate the SPR of that year. The inputs to the LB-SPR included the ratio of the natural mortality coefficient to the von Bertalanffy growth coefficient; length data in each year; the von Bertalanffy asymptotic length parameter; the lengths where 50 and 95% of the fish are mature; and length–weight parameters. LB-SPR estimates the ratio of fishing mortality to natural mortality and the lengths at 50 and 95% selectivity to best fit the predicted and observed length composition proportions, and derives SPR, outputting estimates for these four values for each year with length data (Hordyk et al., 2015). In this study, we used the LB-SPR package version 0.1.5 in R (Hordyk, 2017).

Length-based integrated mixed effects is an age-structured population dynamics mixed model and can account for variable fishing mortality and recruitment when there are only length data from a single year as well as assumed life-history parameters; including the length-at-age relationship, von Bertalanffy growth parameters, allometric length-weight parameters, natural mortality coefficient, and length at 50% maturity (Rudd and Thorson, 2018). It can also accommodate multiple years and types of data, including both length data and index data and/or catch data, in an integrated manner to improve the estimation of changes in fishing mortality over time (Rudd and Thorson, 2018). LIME estimates lengths at 50 and 95% selectivity to the fishing gear; the Dirichlet-multinomial parameter; the recruitment standard deviation and annual fishing mortality coefficient as fixed effects; and estimates of annual recruitment as random effects. Compared with LB-SPR, when fitting length data from more than 1 year, LIME does not assume equilibrium conditions when recruitment is estimable. In our model runs, when the final gradient for all parameters was less than 0.001, the model converged. For each combination of life-history type, data availability scenario, fishing mortality pattern, and recruitment dynamics, we obtained 100 iterations of generated data and ran the estimation model for each set. In this study, we used the LIME package 2.1.3 version (Rudd and Thorson, 2018).

To assess the performance of the LIME and LB-SPR models under different scenarios we compared the outputs to the simulated “truth” from the operating models, and used the relative error (RE, [estimated-true]/true) as an evaluation indicator.

The individual models combined as ensembles were intended to provide estimates of recent stock status (SPR). Therefore, the ensembles were applied to estimate SPR in the last year of the data series, which was of interest for both of management and conservation purposes. The SPR were used in the last year as the response variable and the predictions from the individual models (LIME and LB-SPR) as predictors in our ensemble models (Figure 2).

A model average for each population and four ensembles of varying complexity were compared: a random forest (RF), a support vector machine (SVM), a linear model (LM), and a boosted regression tree (generalized boosted regression modeling, GBM). The four ensemble models were fitted for each population separately (e.g., RF_PA for P. argenteus only, RF_ST for S. taty only), as well as fitting the models for a set comprising data from populations of both species (e.g., RF for both P. argenteus and S. taty).

These models can be described in terms ofθ^, the ensemble estimated SPR. The individual model estimates of SPR are represented as x^. The model average for each population was calculated as:

The linear model for each individual population or both populations combined was calculated as:

Support vector machines, based on mapping input space to a high-dimensional feature space where linear separation is easier than input space, provide a popular machine learning method and yet also represent a powerful technique for general (nonlinear) classification, regression and outlier detection with an intuitive model representation (Cortes and Vapnik, 1995; Bennett and Campbell, 2000; Chang and Lin, 2011). SVMs were initially used to train a data set (SPR in the last year of individual models fitting 80% of the iterations for simulated data under all fishing and recruitment scenarios) to obtain a model and secondly, the resultant ensemble model was used to predict information from a testing data set (SPR in the last year of individual models from fitting the remaining 20% iterations of simulated data under all fishing and recruitment scenarios). SVMs, based the pre-processing strategy in learning by mapping input space to a high-dimensional feature space, are used to find an optimal hyper-plane which maximizes the margin between itself and the nearest training examples in the new high-dimensional space and minimizes the expected generalization error (Seo, 2007). We fit SVMs with the e1071 package (Meyer et al., 2021) for R.

Generalized Boosted Regression Modeling (GBM), based on regression trees, is a type of regression model that is highly flexible, and has considerable success in predictive accuracy by maintaining a monotonic relationship between the response and each predictor (Friedman and Tibshirani, 2000; Ridgeway, 2007; Al-Mudhafar et al., 2016). GBM via the R gbm package (Ridgeway, 2007) was fitted using the default argument values.

Random forest (RF) modeling is also based on regression trees. In RF, each tree depends on the values of a random vector sampled independently and with the same distribution for all trees in the forest (Breiman, 2001). We fitted RF using the random Forest package (Liaw and Wiener, 2002) for R with the default argument values.

Repeated onefold cross-validation was used by randomly dividing the dataset into two sets, then building the ensemble model on eight-tenths of the data and evaluating its predictive performance on the remaining two-tenths. To assess the ability of the ensemble models to accurately and precisely estimate quantities of management interest i.e., SPR, the median absolute relative error (MARE) was used as an evaluation indicator to quantify precision between estimated and true SPR in the last year of data across the 20% iterations of simulated data.

As an example, the LIME and LB-SPR models were applied to stocks of species of interest (i.e., P. argenteus and S. taty in the Bohai Sea) and then the resultant SPR estimates of stock status were used as data in the previously built ensemble models (i.e., RF, RF_PA, RF_ST, SVM, SVM_PA, SVM_ST, GBM, GBM_PA, GBM_ST, LM, LM_PA, LM_ST, and Model average). The Bohai Sea is an important spawning ground for the Bohai-Yellow Sea stocks of P. argenteus and S. taty thereby playing an important role in future years’ stock recruitment (Liu et al., 1990; Jin et al., 2005a,2006, 2015). May to July is the spawning period for the Bohai-Yellow Sea stock of P. argenteus and July to November is the main feeding period. At the end of November, the stock gradually migrates toward the Yellow Sea for overwintering (Liu et al., 1990; Jin et al., 2005a,2006). In the Bohai Sea, S. taty has been the dominant species since the 1980s and plays an important role in the nearshore fish community (Jin et al., 2006). Setipinna taty spawn in the Bohai Sea from mid-May to June, and then return to the wintering ground in the Yellow Sea until the end of November (Liu et al., 1990; Jin et al., 2006, 2015).

The application of ensemble models in this study was based on P. argenteus and S. taty length data that were collected by fixed-station bottom trawl surveys conducted in the Bohai Sea during spring (May), summer (August- September), autumn (October- November) and winter (December) of 2016–2019 by the Yellow Sea Fisheries Research Institute, Chinese Academy of Fishery Sciences (Figure 4).

The modeling results can potentially inform a harvest strategy that targets a fishing mortality rate that is expected to result in 40% of unfished spawning output (termed “SPR_40%”), which is considered risk-averse for many stocks with very low resilience (Clark, 2002). SPR_30% is considered a threshold below which a stock is considered to have been overfished (Clark, 2002; Nadon et al., 2015; Rudd and Thorson, 2018). Therefore, these values were calculated as possible fishing mortality reference points to evaluate the status of the P. argenteus and S. taty in the Bohai Sea to inform future fishery management decisions.

Simulation testing demonstrated that the RE distributions of the LIME and LB-SPR methods under varying conditions were relatively dispersed compared with those under equilibrium conditions (Figure 5). In addition, the RE distribution for P. argenteus was relatively dispersed compared with that for S. taty. When the time-period of data becomes shorter, the bias (bias is measured as MRE) of the LIME estimator for both P. argenteus and S. taty will become larger.

Our results showed that for P. argenteus, both LIME and LB-SPR can estimate unbiased SPR when length data are available and biological characteristics are correctly specified across various scenarios of fishing mortality and recruitment patterns (Figure 5). Compared with P. argenteus, the SPR estimated for S. taty by both LIME and LB-SPR had a large relative error. For P. argenteus, the mean bias in SPR of LIME and LB-SPR was 0.008 (from −0.005 to 0.017) across all data scenarios under equilibrium conditions (constant scenario). Under an endogenous scenario, the mean bias in SPR of LIME and LB-SPR for P. argenteus was 0.014 (from −0.026 to 0.085) across all data scenarios, and the SPR from LIME demonstrated relatively large deviations under scenarios with only length and catch data (C_LC). Under a one-way scenario, the mean bias in the SPR of LIME and LB-SPR for P. argenteus was 0.088 (from 0.039 to 0.163) across all data scenarios, and the SPR of LIME also demonstrated relatively large deviations under scenarios with only length and catch data (C_LC). In general, LB_SPR and LIME with length only data had a similar bias for P. argenteus.

For S. taty, the mean bias in the SPR of LIME and LB-SPR was 0.331 (from −0.176 to 0.405) across all data scenarios under equilibrium conditions. Under endogenous scenarios, the mean bias in SPR of LIME and LB-SPR for S. taty was 0.243 (from −0.146 to 0.379) across all data scenarios. But those RE distributions had relatively large deviations compared with those obtained under equilibrium conditions. Under a one-way scenario, the mean bias in SPR of LIME and LB-SPR for S. taty was 0.348 (from −0.217 to 0.576) across all data scenarios. In general, for S. taty, the LB-SPR estimate was lower than its true value, while the LIME estimate was higher than its true value. It seems that a persistent bias exists for both individual length-based assessment methods compared with the true values in the OM.

Ensemble methods, and in particular the machine-learning ensemble models (support vector machines, SVMs), generally improved the estimates of stock status beyond those of any individual model (Table 2). Compared to the LIME models with LC20, LC10, LC5, LC1, and SVM ensembles (SVM20, SVM10, SVM5, and SVM1) decreased the MARE (median absolute relative error to quantify precision) by averages of 0.257, 0.252, 0.243, 0.130, respectively for S. taty and −0.014, 0.004, 0.011, 0.036, respectively for P. argenteus. Compared to the LB-SPR model, SVM ensembles (SVM20, SVM10, SVM5, and SVM1) decreased the MARE by averages of 0.161, 0.150, 0.148, 0.136, respectively for S. taty and 0.021, 0.019, 0.011, −0.006, respectively for P. argenteus. Compared to the LB-SPR model, only SPR from individual models with the simulation of individual species (S. taty or P. argenteus) can be used to build ensemble models with higher accuracy than those fitted by both species in combination.

We applied ensemble models using the SPR from LIME and LB-SPR with the length composition data from S. taty and P. argenteus from the surveys conducted in the Bohai Sea during 2016–2019 (Table 3). The optimal ensemble model (SVM) showed that the SPR estimates for P. argenteus and S. taty in 2019 were 0.388 and 0.525, respectively. These results (> SPR_30%) indicated that there was no evidence of overfishing for P. argenteus and S. taty in the Bohai Sea. Notwithstanding this result, the SPR of P. argenteus (0.388< SPR_40%) indicated that extra monitoring attention is needed, and further protection warranted to reduce the risk of this stock becoming overfished. Most ensemble models tested produced consistent results.