Effects of misrepresentative length samples on individual growth and stock condition estimates

INTRODUCTION

The importance of obtaining accurate estimates of individual growth parameters is reflected in a large amount of scientific literature related to fisheries, aquaculture and ecology (Brunel and Dickey-Collas 2010Brunel T., Dickey-Collas M. 2010. Effects of temperature and population density on von Bertalanffy growth parameters in Atlantic herring: a macro-ecological analysis. Mar. Ecol. Prog. Ser. 405: 15-28. https://doi.org/10.3354/meps08491 , Hutchinson and TenBrink 2011Hutchinson C.E., TenBrink T.T. 2011. Age determination of the yellow Irish lord: management implications as a result of new estimates of maximum age. N. Am. J. Fish. Manag. 31: 1116-1122. https://doi.org/10.1080/02755947.2011.646453 , Lee et al. 2020Lee L., Atkinson D., Hirst A.G., Cornell S.J. 2020. A new framework for growth curve fitting based on the von Bertalanffy growth function. Sci. Rep. 10: 7953. https://doi.org/10.1038/s41598-020-64839-y ). Often, owing to the selectivity of fishing gear, samples do not represent the complete size structure in the population (Goodyear 1995Goodyear C.P. 1995. Mean size at age: An evaluation of sampling strategies with simulated Red Grouper data. Trans. Am. Fish. Soc. 124: 746-755. https://doi.org/10.1577/1548-8659(1995)124<0746:MSAAAE>2.3.CO;2 , 2019Goodyear C.P. 2019. Modeling growth: consequences from selecting samples by size. Trans. Am. Fish. Soc. 148: 528-551. https://doi.org/10.1002/tafs.10152 , Kraak et al. 2019Kraak S.B.M., Haase S., Minto C., Santos J. 2019. The Rosa Lee phenomenon and its consequences for fisheries advice on changes in fishing mortality or gear selectivity. ICES J. Mar. Sci. 76: 2179-2192. https://doi.org/10.1093/icesjms/fsz107 ) even in data-rich scenarios (Frater and Stefansson 2020Frater P.N., Stefansson G. 2020. Comparison and evaluation of approaches aimed at correcting or reducing selectivity bias in growth parameter estimates for fishes. Fish. Res. 225: 105464. https://doi.org/10.1016/j.fishres.2019.105464 ). This can be problematic given that individual growth influences estimates of mortality, fecundity, condition factor, growth performance, structure, dynamics and variability of stocks, food webs and ecological networks (Tsoukali et al. 2016Tsoukali S., Olsson K.H., Visser A.W., MacKenzie B.R. 2016. Adult lifetime reproductive value in fish depends on size and fecundity type. Can. J. Fish. Aqua. Sci. 73: 1405-1412. https://doi.org/10.1139/cjfas-2015-0378 , Stawitz and Essington 2018Stawitz C.C., Essington T.E. 2018. Somatic growth contributes to population variation in marine fishes. J. Anim. Ecol. 88: 315-329. https://doi.org/10.1111/1365-2656.12921 , N´Dri et al. 2020N’Dri O.R., Konan Y.A., Bamba M., et al. 2020. Length-weight relationships and condition factor of twenty-four freshwater fish species from Lake Buyo, Côte D’Ivoire. J. Fish. Aqua. Sci. 15: 27-34. https://doi.org/10.3923/jfas.2020.27.34 ). More directly for management, growth parameters can also influence estimated abundance, yield and ecosystem-based management reference points (Parma and Deriso 1990Parma A.M., Deriso R.B. 1990. Dynamics of age and size composition in a population subject to size-selective mortality: Effects of phenotypic variability in growth. Can. J. Fish. Aqua. Sci. 47: 274-289. https://doi.org/10.1139/f90-030 , Jennings and Dulvy 2005Jennings S., Dulvy N.K. 2005. Reference points and reference directions for size-based indicators of community structure. ICES J. Mar. Sci. 62: 397-404. https://doi.org/10.1016/j.icesjms.2004.07.030 , Cope and Punt 2009Cope J.M., Punt A.E. 2009. Length-based reference points for data-limited situations: applications and restrictions. Marine and Coastal Fisheries: Dyna., Man., and Eco. Sci. 1: 169-186. https://doi.org/10.1577/C08-025.1 ).

The von Bertalanffy (vB) individual growth model has been used for diverse species of fishes, mammals, birds, and invertebrates (Lee et al. 2020Lee L., Atkinson D., Hirst A.G., Cornell S.J. 2020. A new framework for growth curve fitting based on the von Bertalanffy growth function. Sci. Rep. 10: 7953. https://doi.org/10.1038/s41598-020-64839-y ). This model is often expressed as L_t $={\mathrm{L}}_{\mathrm{\infty }}\left[1-\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{k}(\mathrm{t}-{\mathrm{t}}_0)\right]$ (Pauly 1979Pauly D. 1979. Gill size and temperature as governing factors in fish growth: a generalization of von Bertalanffy’s growth formula. PhD thesis, Christian-Albrechts-Universität Kiel, Germany, 156 pp.). ${\mathrm{L}}_{\mathrm{\infty }}$ is the asymptotic length, k is a constant representing catabolic stress referred to as the Brody growth coefficient (Hart and Chute 2009Hart D.R., Chute A.S. 2009. Estimating von Bertalanffy growth parameters from growth increment data using a linear mixed effect model, with an application to the sea scallop Placopecten magellanicus. ICES J. Mar. Sci. 66: 2165-2175. https://doi.org/10.1093/icesjms/fsp188 ), and t₀ is a theoretical age when length is zero. Despite the wide applicability of the vB model, it is often difficult to compare growth between different taxa (Brey 1999Brey T. 1999. Growth performance and mortality in aquatic macrobenthic invertebrates. Adv. Mar. Biol. 35: 153-223. https://doi.org/10.1016/S0065-2881(08)60005-X ), and there have been several attempts to address this problem (e.g. the index of Gallucci and Quinn II 1979Gallucci V.F., Quinn II. T.J. 1979. Reparametrizing, fitting, and testing a simple growth model. Trans. Am. Fish. Soc. 108: 14-25. https://doi.org/10.1577/1548-8659(1979)108<14:RFATAS>2.0.CO;2 ). A commonly used length-based index of growth performance is ϕ’= ${\mathrm{l}\mathrm{o}\mathrm{g}}_{10}\mathrm{k}+2\mathrm{ }{\mathrm{l}\mathrm{o}\mathrm{g}}_{10}\mathrm{ }{\mathrm{L}}_{\mathrm{\infty }}$ (Pauly and Munro 1984Pauly D., Munro J.L. 1984. Once more on the comparison of growth in fish and invertebrates. Fishbyte 2: 21.), which is a species-specific index used to compare reliability of vB parameters between and within species or stocks (Etim et al. 1999Etim L., Lebo P.E., King R.P. 1999. The dynamics of an exploited population of a siluroid catfish (Schilbe intermidius Reupell 1832) in the Cross River, Nigeria. Fisheries Research 40, 295-307. https://doi.org/10.1016/S0165-7836(98)00217-3 , Moura et al. 2017Moura P., Garaulet L.L., Vasconcelos P. et al. 2017. Age and growth of a highly successful invasive species: the Manila clam Ruditapes philippinarum (Adams & Reeve, 1850) in the Tagus Estuary (Portugal). Aquat. Invasions 12: 133-146. https://doi.org/10.3391/ai.2017.12.2.02 ). The growth performance parameter ϕ’ has found widespread use in comparing integral performance of vB growth curves (Quaas and Skonhoft 2022Quaas M., Skonhoft A. 2022. Welfare effects of changing technological efficiency in regulated openaccess fisheries. Env. Res. Econ. 82: 869-888. https://doi.org/10.1007/s10640-022-00693-y , Rodríguez-Castañeda et al. 2022Rodríguez-Castañeda J.C., Ventero A., Manuela García-Márquez G., Iglesias M. 2022. Spatial and temporal analysis (2009-2020) of the biological parameters, abundance and distribution of Trachurus mediterraneus (Steindachner, 1868) in the western Mediterranean. Fish. Res. 256, 106483. https://doi.org/10.1016/j.fishres.2022.106483 , Şimşek et al. 2022Şimşek E., Gözler Z. A., Samsun O. 2022. Age and growth parameters of African catfish (Clarias gariepinus Burchell, 1822) from Asi River, Turkey. Acta Nat. Sci. 3 32-43. https://doi.org/10.29329/actanatsci.2022.351.04 ).

In addition to growth parameters and growth performance, the condition factor (Kn) is an important index in fisheries biology and allows inferences about the fitness of an individual in a population (N’Dri et al. 2020N’Dri O.R., Konan Y.A., Bamba M., et al. 2020. Length-weight relationships and condition factor of twenty-four freshwater fish species from Lake Buyo, Côte D’Ivoire. J. Fish. Aqua. Sci. 15: 27-34. https://doi.org/10.3923/jfas.2020.27.34 ). Kn can be expressed as W₀/ $\widehat{\text{P}}$ , where $\widehat{\text{P}}$ is the weight estimated with the length-weight relationship and W₀ is observed weight; this expression is known as the relative condition factor (Le Cren 1951Le Cren E.D. 1951. The length-weight relationship and seasonal cycle in gonad weight and condition in the perch (Perca fluviatilis). J. Anim. Ecol. 20(2): 201-219. https://doi.org/10.2307/1540 ). Kn can be used to compare the status of conspecific organisms or the status between species, sexes and sizes and in different seasons of the year or between years. Individuals are considered to be in relatively good condition when Kn is greater than 1 and in poor condition if Kn is lower than 1 (Jisr et al. 2018Jisr N., Younes G., Sukhn C., El-Dakdouki M.H. 2018. Length-weight relationships and relative condition factor of fish inhabiting the marine area of the Eastern Mediterranean city, Tripoli-Lebanon. Egypt. J. Aqua. Res. 44, 299-305. https://doi.org/10.1016/j.ejar.2018.11.004 ).

In the present study we explore how simulated length structure in samples affects estimates of key stock condition parameters, with an illustrative example. We consider a sample of lengths to be biased if it consistently over- or underrepresents the entire stock size structure. Under various sampling scenarios, we compute the accuracy of vB individual growth parameters and how they subsequently affect the weight-length relationship, W_∞ (maximum theoretical weight), ϕ’ and Kn. Our aim is to show the sensitivity of these parameters to the length sampling configuration and derive conclusions to inform sampling methods.

MATERIALS AND METHODS

We simulate growth data for the brown swimming crab, Callinectes bellicosus and compare estimated parameter values with those reported in the literature. In the simulations, individual growth follows a vB model with multiplicative random impacts $\mathrm{\gamma }$ _i. We use the general model L_i $=\mathrm{L}\left({\mathrm{t}}_{\mathrm{i}}\right){\mathrm{\gamma }}_{\mathrm{i}}$ , where L_i is observed length, L(t_i) is median length-at age t_i and γ_i is the random term. Known, baseline parameter values are ${\mathrm{L}}_{\mathrm{\infty }}\text{=190.44}$ mm carapace width (CW), $\text{k}\text{=1.038}{\mathrm{y}}^{-1}$ and t₀= $\text{-0.14}$ y (Villa-Diharce et al. 2021Villa-Diharce E.R., Cisneros-Mata M.Á., Rodríguez-Félix D., Ramírez-Félix E.A., Rodríguez-Domínguez G. 2021. Molting and individual growth models of Callinectes bellicosus. Fish. Res. 239: 105897. https://doi.org/10.1016/j.fishres.2021.105897 ).

The vB model with a multiplicative error is L_i $={\mathrm{L}}_{\mathrm{\infty }}\left[1-\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{k}\left({\mathrm{t}}_{\mathrm{i}}-{\mathrm{t}}_0\right)\right)\right]{\mathrm{\gamma }}_{\mathrm{i}}$ where ${\mathrm{L}}_{\mathrm{i}}$ is the length-at-age ${\mathrm{t}}_{\mathrm{i}}$ with log-mean = 0 and arbitrary log-sd = 0.10. An additive model results from logarithmic transformation of the model ${\mathrm{y}}_{\mathrm{i}}={\mathrm{\mu }}_{\mathrm{i}}+{\mathrm{ϵ}}_{\mathrm{i}}$ , where ${\mathrm{y}}_{\mathrm{i}}=\mathrm{l}\mathrm{o}\mathrm{g}\left({\mathrm{L}}_{\mathrm{i}}\right)$ , ${\mathrm{\mu }}_{\mathrm{i}}=\mathrm{l}\mathrm{o}\mathrm{g}\left\{{\mathrm{L}}_{\mathrm{\infty }}\left[1-\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{k}\left({\mathrm{t}}_{\mathrm{i}}-{\mathrm{t}}_0\right)\right)\right]\right\}$ ${\mathrm{ϵ}}_{\mathrm{i}}~\mathrm{N}\left(0,\mathrm{\sigma }\right)$ . The loglikelihood is $\mathrm{L}\mathrm{L}\left({\mathrm{L}}_{\mathrm{\infty }},\mathrm{k},{\mathrm{t}}_0,\mathrm{\sigma }\right)=-\left(\mathrm{n}/2\right)\mathrm{l}\mathrm{o}\mathrm{g}\left(2\mathrm{\pi }\right)-\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{g}\left(\mathrm{\sigma }\right)-\left(1/2\right){\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left[\left({\mathrm{y}}_{\mathrm{i}}-{\mathrm{\mu }}_{\mathrm{i}}\right)/\mathrm{\sigma }\right]}^2$ (Burnham and Anderson 2002Burnham K.P., Anderson D.R. 2002. Model selection and multi-model inference: A practical information-theoretic approach. Springer, New York, NY, USA.). The maximum likelihood estimators are the parameter values such that they maximize the loglikelihood, that is, $\left({\widehat{\mathrm{L}}}_{\mathrm{\infty }},\widehat{\mathrm{k}},{\widehat{\mathrm{t}}}_0,\widehat{\mathrm{\sigma }}\right)=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{L}\mathrm{L}\left({\mathrm{L}}_{\mathrm{\infty }},\mathrm{k},{\mathrm{t}}_0,\mathrm{\sigma }\right)$ .

We used a stratified sampling scheme to generate random samples. Three length (mm) segments were considered: (26-80), (80-135) and (135-190); these cover the length range of the brown swimming crab. Different sampling schemes could have been used, yet our aim was to obtain simulated samples of different, non-overlapping size segments to contrast results of our parameter estimates. Thus, the most parsimonious scheme in this case is stratified sampling. With these three segments we considered five sampling configurations (see below). We obtained pairs of L_i and their corresponding $\mathrm{\gamma }$ _i to obtain L(t_i)=L_i/ $\mathrm{\gamma }$ _i. that satisfy the restriction L₀< L(t_i)< L_∞; that is, L(t_i) values smaller than 26 and greater than 190 mm were discarded. The age t_i that corresponds to length L_i was obtained by solving the vB model:

Fifty pairs of L_i and γ_i were randomly drawn for each length segment to estimate vB growth and stock condition parameter values; this was repeated 1000 times using the Monte Carlo (MC) method (Janssen 2013Janssen H. 2013. Monte-Carlo based uncertainty analysis: Sampling efficiency and sampling convergence. Reliab. Eng. Syst. Saf. 109: 123-132. https://doi.org/10.1016/j.ress.2012.08.003 ). The means of the following parameters values were estimated: 1) vB growth equation, k, ${\mathrm{L}}_{\mathrm{\infty }}$ and t₀; 2) growth performance index, ϕ’; 3) weight-length relationship, a, b; 4) maximum theoretical weight, ${\mathrm{W}}_{\mathrm{\infty }}$ ; and 5) condition factor, Kn=W₀/ $\widehat{\text{P}}$ . Reference parameters a and b (and their variability) were also estimated in a previous work (Villa-Diharce et al. 2021Villa-Diharce E.R., Cisneros-Mata M.Á., Rodríguez-Félix D., Ramírez-Félix E.A., Rodríguez-Domínguez G. 2021. Molting and individual growth models of Callinectes bellicosus. Fish. Res. 239: 105897. https://doi.org/10.1016/j.fishres.2021.105897 ).

Samples from each length segment were drawn using the following five configurations (Table 1).

Except for Kn, all mean MC-estimated parameter values were compared with the original values in terms of their relative biases, mean squared error and standard errors. To analyse the values, in equation Kn=W₀/ $\widehat{\text{P}}$ we substitute for the observed weight W₀ the expression used for its simulation. We have ${\text{W}}_{\text{0}}\text{=a}{\text{L}}^{\text{b}}\text{λ}$ , where a and b are the observed values of the parameters of the weight-length model. In this equation we assumed a multiplicative error term with lognormal distribution. Substituting terms, we have Kn $={\mathrm{W}}_0/\widehat{\mathrm{a}}{\mathrm{L}}^{\widehat{\mathrm{b}}}=\mathrm{a}{\mathrm{L}}^{\mathrm{b}}\text{λ}/\widehat{\mathrm{a}}{\mathrm{L}}^{\widehat{\mathrm{b}}}=\left(\mathrm{a}/\widehat{\mathrm{a}}\right){\mathrm{L}}^{\left(\mathrm{b}-\widehat{\mathrm{b}}\right)}\text{λ}$ . For each length configuration, random lengths of uniform distributions (within the smallest and largest length) were generated. Weight was then obtained considering multiplicative or lognormal errors with log-mean=0 and cv= 10% (or sd=0.101). This expression shows the influence of discrepancies between a and $\widehat{\mathrm{a}}$ and between b and $\widehat{\mathrm{b}}$ on the value of Kn. To better represent this influence, we tabulated magnitudes for the different sampling configurations. When a and b take their real values, then $\left(\mathrm{a}/\widehat{\mathrm{a}}\right)$ =1, (b- $\widehat{\mathrm{b}}$ )=0, and Kn= $\text{λ}$ .

For each sampling configuration, 50 Kn values were obtained as follows. With the MC-originated lengths, an observed mean weight W₀ was computed using the observed parameters of the weight-length relationship (see below). For each mean length, weight was obtained using the MC parameters a and b obtained for samples from the three length segments. Kn values were then plotted against their corresponding length to observe the behaviour of Kn depending on the sampling scheme. To better understand the behaviour of Kn, we conducted a closer analysis of the differences of the true and MC-estimated a and b values of the weight-length relationship.

As mentioned, the vB parameter values used as baselines were those obtained by Villa-Diharce et al. (2021)Villa-Diharce E.R., Cisneros-Mata M.Á., Rodríguez-Félix D., Ramírez-Félix E.A., Rodríguez-Domínguez G. 2021. Molting and individual growth models of Callinectes bellicosus. Fish. Res. 239: 105897. https://doi.org/10.1016/j.fishres.2021.105897 . Using maximum likelihood, we estimated parameters of the weight-length model $\mathrm{W}=\mathrm{a}{\mathrm{L}}^{\mathrm{b}}$ (Haddon 2011Haddon M. 2011. Modelling and Quantitative Methods in Fisheries 2nd Edition. Chapman and Hall/CRC, New York.406 pp. https://doi.org/10.1201/9781439894170 ) for combined sexes; after observing the dispersion of data we assumed a multiplicative error term (e.g. Curiel-Bernal et al. 2021Curiel-Bernal M.V., Aragón-Noriega E.A., Cisneros-Mata M.Á., et al. 2021. Using observed residual error structure yields the best estimates of individual growth parameters. Fishes 6: 35. https://doi.org/10.3390/fishes6030035 ). The statistical model is ${\mathrm{W}}_{\mathrm{i}}=\mathrm{a}{\mathrm{L}}_{\mathrm{i}}^{\mathrm{b}}{\text{λ}}_{\mathrm{i}}$ with ${\text{λ}}_{\mathrm{i}}$ lognormally distributed with log-mean zero and log-standard deviation σ. We log-transformed the model and obtained an additive model ${\mathrm{y}}_{\mathrm{i}}={\mathrm{\mu }}_{\mathrm{i}}+{\mathrm{ϵ}}_{\mathrm{i}}$ , where ${\mathrm{y}}_{\mathrm{i}}=\mathrm{l}\mathrm{o}\mathrm{g}\left({\mathrm{W}}_{\mathrm{i}}\right)$ , ${\mathrm{\mu }}_{\mathrm{i}}=\mathrm{l}\mathrm{o}\mathrm{g}\left(\mathrm{a}{\mathrm{L}}_{\mathrm{i}}^{\mathrm{b}}\right)$ and ${\mathrm{ϵ}}_{\mathrm{i}}~\mathrm{N}\left(0,\mathrm{\sigma }\right)$ . The loglikelihood function is (Burnham and Anderson 2002Burnham K.P., Anderson D.R. 2002. Model selection and multi-model inference: A practical information-theoretic approach. Springer, New York, NY, USA.) $\mathrm{L}\mathrm{L}\left(\mathrm{a},\mathrm{b},\mathrm{\sigma }\right)=-\left(1/2\right)\log \left(2\mathrm{\pi }\right)-\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{g}\left(\mathrm{\sigma }\right)-\left(1/2\right){\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left(\left({\mathrm{y}}_{\mathrm{i}}-{\mathrm{\mu }}_{\mathrm{i}}\right)/\mathrm{\sigma }\right)}^2$ .

The estimators $\widehat{\mathrm{a}}$ , $\widehat{\mathrm{b}}$ and $\widehat{\mathrm{\sigma }}$ of parameters a, b and σ, respectively, are those that maximize the loglikelihood function, that is $\left(\widehat{\text{a}}\text{,}\widehat{\text{b}}\text{,}\widehat{\text{σ}}\right)\text{=}\text{argmax}\text{LL}\left(\text{a,b,σ}\right)$ . To numerically maximize the loglikelihood we used the function nlminb( ) written in R (R Core Team 2021R Core Team. 2021. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.). The significance of b in the weight-length equation was tested using a t-test (Pauly 1984Pauly D. 1984. Fish population dynamics in tropical waters: a manual for use with programme calculators. ICLARM Studies and Reviews 8. Manila, Philippines. 325 pp.).

Using known values of the vB parameters and of the weight-length model, we estimated the maximum loglikelihood values of ${\mathrm{\varphi }}^{\mathrm{\text{'}}}$ and ${\mathrm{W}}_{\mathrm{\infty }}$ , ${\widehat{\mathrm{\varphi }}}^{\mathrm{\text{'}}}={\mathrm{l}\mathrm{o}\mathrm{g}}_{10}\widehat{\mathrm{k}}+2\mathrm{ }{\mathrm{l}\mathrm{o}\mathrm{g}}_{10}\mathrm{ }{\widehat{\mathrm{L}}}_{\mathrm{\infty }}$ (Pauly and Munro 1984Pauly D., Munro J.L. 1984. Once more on the comparison of growth in fish and invertebrates. Fishbyte 2: 21.) and ${\widehat{\mathrm{W}}}_{\mathrm{\infty }}=\widehat{\mathrm{a}}{\widehat{\mathrm{L}}}_{\mathrm{\infty }}^{\widehat{\mathrm{b}}}$ (Haddon 2011Haddon M. 2011. Modelling and Quantitative Methods in Fisheries 2nd Edition. Chapman and Hall/CRC, New York.406 pp. https://doi.org/10.1201/9781439894170 ).

We compared the quality of a parameter estimator H using the mean squared difference between the estimator H and the parameter θ, i.e. the mean square error (MSE) of parameter H. The MSE can be divided into variance and bias (Casella and Berger 1990Casella G., Berger R.L. 1990. Statistical Inference. Wadsworth & Brooks/Cole. Pacific Grove, CA, USA. 688 pp.): $\text{MSE}\left(\text{H}\right)\text{=}{\text{E}}_{\text{θ}}{\left(\text{H-θ}\right)}^{\text{2}}\text{=}{\text{V}}_{\text{θ}}\left(\text{H}\right)\text{+}{\left[{\text{Bias}}_{\text{θ}}\left(\text{H}\right)\right]}^{\text{2}}$ . For a sequence H₁, H₂, ... H_n of estimates of a parameter θ, one can obtain the terms of the previous relationship as $\left(\text{1}/\text{n}\right)\sum _{\text{i=1}}^{\text{n}}{\left({\text{H}}_{\text{i}}\text{-θ}\right)}^{\text{2}} \text{=}\left(\text{1}/\text{n}\right)\sum _{\text{i=1}}^{\text{n}}{\left({\text{H}}_{\text{i}} \text{-}\stackrel{\mathrm{-}}{\text{H}}\right)}^{\text{2}}\text{+}{\left[\stackrel{\mathrm{-}}{\text{H}} \text{- θ}\right]}^{\text{2}}$ .

We then took the square root of these quantities so that they are expressed in the same scale of the estimated parameter: square root of the mean error $\mathrm{S}\text{RME}\left(\text{H}\right)\text{=}\sqrt{\left(\text{1}/\text{n}\right)\sum _{\text{i=1}}^{\text{n}}{\left({\text{H}}_{\text{i}}\text{-θ}\right)}^{\text{2}}}$ , standard error $\text{SE}\left(\text{H}\right)\text{=}\sqrt{\left(\text{1}/\text{n}\right)\sum _{\text{i=1}}^{\text{n}}{\left({\text{H}}_{\text{i}}\text{-}\stackrel{\mathrm{-}}{\text{H}}\right)}^{\text{2}}}$ and Bias $\text{Bias}\left(\text{H}\right)\text{= (}\stackrel{-}{\mathrm{H}}\text{-θ})$ . To further appreciate the magnitudes, we took their value relative to the magnitude of the parameter to be estimated $\text{θ}$ (Lehmann and Casella 1998Lehmann E.L., Casella G. 1998. Theory of point estimation. 2nd ed. Springer. New York, NY, USA. 588 pp., Dekking et al. 2005Dekking F.M., Kraaikamp C., Lopuhaä H.P., Meester L.E. 2005. A modern introduction to probability and statistics: understanding why and how. Springer-Verlag, London, UK. 486 pp.https://doi.org/10.1007/1-84628-168-7 , Wang et al. 2021Wang K., Zhang C., Ming, S. et al. 2021. Fishing pressure and lifespan affect the estimation of growth parameters using ELEFAN. Fish. Res. 238. https://doi.org/10.1016/j.fishres.2021.105903 ). These relative values were then tabulated with columns referring to the segment sampled and rows showing the relative quantities estimated, as a percentage of the original values.

Finally, with each set of vB parameter values estimated sampling the five configurations of lengths, a plot was generated and compared with the baseline curve using the best estimates (Villa-Diharce et al. 2021Villa-Diharce E.R., Cisneros-Mata M.Á., Rodríguez-Félix D., Ramírez-Félix E.A., Rodríguez-Domínguez G. 2021. Molting and individual growth models of Callinectes bellicosus. Fish. Res. 239: 105897. https://doi.org/10.1016/j.fishres.2021.105897 ). This provided an integrated insight into biases and errors that can be committed with incomplete length sampling.

RESULTS

The weight-length relationship estimated was W=0.000017939 L^3.349, b being significantly greater than 3 (p<0.001) (CI: 3.2923, 3.4059). Using the mentioned vB parameters, we estimated the maximum weight ${\mathrm{W}}_{\mathrm{\infty }}=0.000017939\mathrm{ }{\mathrm{L}}_{\mathrm{\infty }}^{3.349}=773.94\mathrm{ }\mathrm{g}$ . The baseline value of phi prime was ϕ’= ${\mathrm{l}\mathrm{o}\mathrm{g}}_{10}1.038+2\mathrm{ }{\mathrm{l}\mathrm{o}\mathrm{g}}_{10}190.44$ = 4.58.

The mean, relative bias (RB), relative standard error (RSE) and relative root mean square error (RMSE) of the parameters estimated under the tested sampling schemes are shown, respectively, in Tables 2 to 5. The following section highlights the most relevant information contained in these five tables.

We obtained the minimum value of ${\text{L}}_{\text{∞}}$ when sampling only from the smallest length segment (26-80 mm) and the nearest to the baseline value using extreme lengths, 26-80 and 135-190 (Table 2). When sampling the whole range of lengths (Table 3), the smallest plus largest, and the largest lengths yielded relatively small bias values (-0.24 to -3.38 mm, 26-80 and 135-190 and 135-190, respectively). Sampling the smallest length range produced the highest bias, followed by sampling the central range. Small values of the RSE (Table 4), from 3.22 to 3.78%, were found for symmetric sampling (whole range, and smallest and largest). When sampling only the largest range, the variability of ${\text{L}}_{\text{∞}}$ more than doubled. The RMSE’s two components (bias and RSE) behaved similarly (Table 5).

Estimates of parameter k showed small biases when samples came from segments 26-80 and 135-190 mm (5.89 and 4.32%) (Table 3). The smallest bias occurred when samples came from the largest length segment. In general, when only one segment was sampled, the RB was always large. The size of the biases when sampling combined and separated segments was also observed in the RSE values (Table 4).

For t₀ the largest bias was negative when samples came from the largest length segment, 135-190 mm (-459.66%). Relatively small biases resulted when samples came from symmetric combinations of segments; when sampling separate segments, biases were larger (Table 3). The RSE (variability) values showed the same described pattern: lower values when sampling symmetric combinations of length segments, and notably larger values when sampling separate segments. Most of the magnitude of RMSE was due to variability of RSE.

Estimates of growth performance ϕ’ were stable and small, changing from positive to negative when samples came from single length segments. When sampling came from combined segments, the global quality index, RMSE, took very small positive values (Table 5).

For the weight-length (W-L) relationship, the bias of coefficient a was smaller when samples came from segments containing the whole range of lengths, followed by combined smallest and largest (Table 3). The RSE behaved similarly (Table 4): an increasing RMSE of coefficient a resulted when samples came from individual segments (Table 5). Estimates of the exponent b were stable, as can be seen both in the RB values and in the RSE values (Tables 3 and 4). This stability was also observed in the values of RSE and RMSE (Tables 4 and 5).

The behaviour of estimated ${\mathrm{W}}_{\mathrm{\infty }}$ was similar to that of ${\mathrm{L}}_{\mathrm{\infty }}$ : both had smaller biases and standard errors in configurations with extreme lengths (1 and 2) and in the larger lengths (last configuration). ${\mathrm{W}}_{\mathrm{\infty }}$ showed a greater variability than ${\mathrm{L}}_{\mathrm{\infty }}$ because of the variability of both ${\mathrm{L}}_{\mathrm{\infty }}$ and the scale parameter a of the W-L relationship. The value of parameter b did not influence variation of ${\mathrm{W}}_{\mathrm{\infty }}$ because of its stability (Tables 2-5).

Figure 1 shows the scatter plots of 50 Kn values resulting from lengths simulated in the five configurations. The upper-left panel is considered as reference, i.e. when there is no bias in the estimates of parameters a and b. As observed, the Kn values were relatively well estimated, except when samples came from the largest length range. In this case, Kn values were underestimated.

Table 6 shows the estimated means of parameters a and b for the five sampling configurations, as well as their comparisons with the true values of the parameters. The first two configurations yielded similar estimates of the true values of a and b.

In the first two sampling configurations, the ratios were close to one and the differences were almost zero. This is why the plots of Kn for the first two configurations in Figure 1 are very similar to those presented in the reference plot (upper-left). Sampling from the central and larger length segments resulted in monotonic growth in the estimates $\widehat{\text{a}}$ of a, and therefore monotonic decrease in the ratios $\text{a}/\widehat{\text{a}}$ , which ranged from 0.929 to 0.854.

A comprehensive visual analysis of the relative performance of the five sampling schemes can be observed by comparing the growth plots obtained using estimated parameters and plots using the baseline values (Fig. 2). The least biased plot resulted when samples contained the whole range of lengths (configuration 1) and with the smallest and largest length segments (configuration 2). Poor fits resulted when samples came from configurations 3, 4 and 5; the worst fit resulted when samples came only from the smallest length range.

DISCUSSION

Sampling different size segments of a stock or population will influence the estimates of individual growth, length-weight parameters, maximum theoretical individual weight and two widely used measures of growth efficiency: growth performance index ϕ’ and condition factor Kn. Using three relative measures of accuracy (bias, standard error and mean square error) and an a priori segmentation of length, we obtained results with practical applications. The standard deviation (sd) of 0.1 used here covers +/-20% of possible values around the mean; however, due to chance, some values could lay outside such boundaries. Accuracy of parameter estimates can vary when different sd values are used, and this issue merits further research. The accuracy of parameter estimates does not necessarily depend on sampling the entire length-age range possible, and the error is not the same for all parameters. Our results provide useful guidance to develop sampling schemes in the common case when time and resources are scarce. We caution that erred estimates of basic parameters, particularly von Bertalanffy individual growth, can lead to wrong values of other key parameters used in fisheries management, for example, M natural mortality rate (see Maunder et al. 2023Maunder M.N., Hamel O.S., Lee H-H., et al. 2023. A review of estimation methods for natural mortality and their performance in the context of fishery stock assessment. Fish. Res. 257: 106489. https://doi.org/10.1016/j.fishres.2022.106489 for a recent review of various methods to estimate M).

To estimate the accuracy of parameters in our simulations we used three relative indices based on length-stratified age samples. Our main purpose in this paper was to simulate and compare samples representing a balanced number of biased length samples to analyse possible effects in estimated values of parameters of general interest. Much larger sample sizes can and are often obtained; in our case, however, we struggled to provide insights for the common case of data-poor fisheries or limited resources for sampling. Other works (Xiao 1996Xiao Y. 1996. A framework for evaluating experimental designs for estimating rates of fish movement from tag recoveries. Can. J. Fish. Aqua. Sci. 53: 1272-1280. https://doi.org/10.1139/f96-050 , Perreault et al. 2020Perreault A.M.J., Zheng N., Cadigan N.G. 2020. Estimation of growth parameters based on length-stratified age samples. Can. J. Fish. Aqua. Sci. 77: 439-450. https://doi.org/10.1139/cjfas-2019-0129 ) used relative root mean squared error and RB of simulated and true values. In our case, we also used the RSE (Barbaro et al. 1981Barbaro S., Cannata G., Coppolino S., et al. 1981. Diffuse solar radiation statistics for Italy. J. S. Energy 26: 429-435. https://doi.org/10.1016/0038-092X(81)90222-X , Vølstad et al. 2011Vølstad J.H., Korsbrekke K., Nedreaas K.H., et al. 2011. Probability-based surveying using self-sampling to estimate catch and effort in Norway’s coastal tourist fishery. ICES J. Mar. Sci. 68: 1785-1791. https://doi.org/10.1093/icesjms/fsr077 ), another measure of accuracy of the parameter estimates as a function of the sampling configurations we tested. This statistic evenly distributes the deviations between sampling configurations.

Incomplete length sampling may be caused by selectivity of fishing gear and when catch is graded at sea or upon landing before samples can be taken. Natural behaviour of individuals may also cause misrepresentation of length structures in samples. For the crab Ranina ranina, for example, juveniles and adults spend considerable time buried and segregate by life history stage. Gear used for commercial fishing seldom catches juveniles, which can produce erred estimates of vB growth parameters (Kirkwood et al. 2005Kirkwood J.M., Brown I.W., Gaddes S.W., Hoyle S. 2005. Juvenile length-at-age data reveal that spanner crabs (Ranina ranina) grow slowly. Mar. Biol. 147: 331-339. https://doi.org/10.1007/s00227-005-1574-0 ).

Failure to produce robust estimates of growth parameters will inevitably curtail our ability to conduct good stock assessments to inform management (Gwinn et al. 2010Gwinn D.C., Allen M.S., Rogers M.W. 2010. Evaluation of procedures to reduce bias in fish growth parameter estimates resulting from size-selective sampling. Fish. Res. 105: 75-79. https://doi.org/10.1016/j.fishres.2010.03.005 ). For example, given that fecundity generally increases with individual size (Marshall et al. 2019Marshall D.J., Gaines S., Warner R., et al. 2019. Underestimating the benefits of marine protected areas for the replenishment of fished populations. Front. Ecol. Environ. 17: 407-413. https://doi.org/10.1002/fee.2075 ), to maximize catch from a cohort, a yield-per-recruit analysis is often performed (Beverton and Holt 1957Beverton R.J.H., Holt S.J. 1957. On the dynamics of exploited fish populations. Fish. Invest., London 19: 1-533., Die et al. 1988Die D.J., Restrepo V.R., Hoenig J.M. 1988. Utility-per-recruit modeling: A neglected concept. Trans. Am. Fish. Soc. 117: 274-281. https://doi.org/10.1577/1548-8659(1988)117<0274:UMANC>2.3.CO;2 , Zhai and Pauly 2019Zhai L., Pauly D. 2019. Yield-per-recruit, utility-per-recruit, and relative biomass of 21 exploited fish species in China’s coastal seas. Front. Mar. Sci. 6: 724. https://doi.org/10.3389/fmars.2019.00724 ), which depends on individual growth parameters.

In the present work we sought parameter values that were closest to the true values that represented correct growth trajectories (Pardo et al. 2013Pardo S.A., Cooper A.B., Dulvy N.K. 2013. Avoiding fishy growth curves. Methods Ecol. Evol. 4: 353-360. https://doi.org/10.1111/2041-210x.12020 ). Previous simulation studies (e.g. Wilson et al. 2015Wilson K.L., Matthias B.G., Barbour A.B., et al. 2015. Combining samples from multiple gears helps to avoid fishy growth curves. N. Am. J. Fish. Man. 35: 1121-1131. https://doi.org/10.1080/02755947.2015.1079573 ) recommend combining samples from fishing gears with different selectivity to improve growth parameter estimates. It has been also proposed that to maximize accuracy of individual growth parameters, a complete representation of organisms from different sizes is needed, with two conditions: 1) evenly distributed sample sizes across age/size segments, and 2) sample sizes as large as possible (Quinn and Deriso 1999Quinn T.J., Deriso R.B. 1999. Quantitative fish dynamics. Oxford University Press, New York, USA. 542 pp., Pilling et al. 2002Pilling G.M., Kirkwood G.P., Walker S.G. 2002. An improved method for estimating individual growth variability in fish, and the correlation between von Bertalanffy growth parameters. Can. J. Fish. Aqua. Sci. 59: 424-432. https://doi.org/10.1139/f02-022 , Shelton and Mangel 2012Shelton A.O., Mangel M. 2012. Estimating von Bertalanffy parameters with individual and environmental variations in growth. J. Biol. Dyn. 6, sup2: 3-30. https://doi.org/10.1080/17513758.2012.697195 ). Our simulation results indicate that more subtle characteristics underly the final estimated values of growth parameters. A key element is to a priori take into consideration the configuration of possible size segments sampled (e.g. Table 1). Using simulations, Goodyear (1995)Goodyear C.P. 1995. Mean size at age: An evaluation of sampling strategies with simulated Red Grouper data. Trans. Am. Fish. Soc. 124: 746-755. https://doi.org/10.1577/1548-8659(1995)124<0746:MSAAAE>2.3.CO;2 concluded that reliable estimates of mean size-at age require random sampling of lengths within ages, and that stratifying samples by length biased the estimates of mean length-at-age. It is not clear if samples were generated by splitting age into equal or different sizes. Goodyear (2019)Goodyear C.P. 2019. Modeling growth: consequences from selecting samples by size. Trans. Am. Fish. Soc. 148: 528-551. https://doi.org/10.1002/tafs.10152 simulated samples using two strategies relevant for the present work: samples stratified by age and by length; size-stratification produced biased estimates of length-at-age and vB parameters. In this case, age strata were of one year, and length strata were constant.

In general, for the vB growth model, it was found that underrepresented small/large individuals yield small-biased k/large-biased ${\mathrm{L}}_{\mathrm{\infty }}$ estimates, respectively (Taylor et al. 2005Taylor N.G., Walters C.J., Martell S.J.D. 2005. A new likelihood for simultaneously estimating von Bertalanffy growth parameters, gear selectivity, and natural and fishing mortality. Can. J. Fish. Aqua. Sci. 62: 215-223. https://doi.org/10.1139/f04-189 ). Because of the inverse relationship between k and ${\mathrm{L}}_{\mathrm{\infty }}$ (Gubiani et al. 2012Gubiani E.A., Gomes L.C., Agostinho A.A. 2012. Estimates of population parameters and consumption/biomass ratio for fishes in reservoirs, Paraná State, Brazil. Neotrop. Ichthyol. 10: 177-188. https://doi.org/10.1590/S1679-62252012000100017 ), biases of these two parameters vary in opposite directions. Also, stability of ϕ’ with respect to the sampling segments results from the inverse correlation between the base-10 logarithms of estimated parameters ${\mathrm{L}}_{\mathrm{\infty }}$ and k (Pauly 1998Pauly D. 1998. Beyond our original horizons: the tropicalization of Beverton and Holt. Rev. Fish Biol. Fish. 8: 307-334. https://doi.org/10.1023/A:1008863215253 ). If one parameter decreases, the other increases, so the product of terms that define ϕ’ remains stable. Useful estimates of ϕ’ requires sound sampling that, whenever possible, accounts for seasonal or annual variations in length composition of stocks (Mathews and Samuel 1990Mathews C.P., Samuel M. 1990. The relationships between maximum and asymptotic length in fishes. Fishbyte 8(2): 14-16.).

Estimated Kn depends on values of the weight-length function. In the present work, because the $b-\widehat{\mathrm{b}}$ are very close to zero and their contribution to Kn is through an exponential function, their effect is practically negligible. Hence, the values of Kn are merely a reflection of the values taken by the ratios $\text{a}/\widehat{\text{a}}$ . In practical terms, symmetric sampling configurations that include extreme length values result in practically unbiased estimates of parameters a and b (cf. Fig. 1). Samples with the largest individuals increased bias in coefficient a. This is in turn reflected in the scatter plots of Kn. For parameter b, when samples came from separate segments, RSE values were very similar to the RMSE, which means that the variability component is greater than the bias.

Our simulations showed that erroneous vB parameter estimates result when most or all individuals in the samples are of similar lengths. For practical purposes, when the vB model fits the data and sampling resources are scarce, it is convenient to actively include the smallest and largest individuals in a sample (segments 1 and 3). If researchers are interested in estimating growth performance ϕ’, it would be advisable to sample from the entire size range available using different fishing gears or sampling methods. These considerations are intended to guide sampling schemes and minimize erred estimates of growth and growth efficiency arising from lack of appropriate data.