Estimating fish growth is a fundamental step to evaluate fish populations and plays a crucial role in stock assessment models and fisheries management (Lee et al., 2024). Generally, growth can be defined as the change in length or weight over the life of an individual and directly influences fish survival, sexual maturity, reproductive success, as well as movement and migration (Peters, 1986). Accurate growth estimation enhances the understanding of fish population dynamics, productivity, and sustainability, particularly for short-lived, small pelagic species like anchovy, which are ecologically and commercially significant (Costalago et al., 2011; Gebremedhin et al., 2021). Within most stock assessment models used, growth is typically estimated through age-based approaches, relying on age determination techniques such as otolith analysis or length–frequency data (Rodríguez Mendoza, 2006). These models, when incorporating data for all the ages and lengths, help to evaluate the size composition and age structure and to estimate biological parameters accounting for growth changes produced by the effect of fishing. These advantages, among others (Lee et al., 2024), facilitate an informed decision-making process (Cope et al., 2023). However, estimating fish growth within stock assessment models may also present limitations since these models, used to handle a large amount of data, often assume that growth is uniform across the entire population, provide parameters not representative of the species biology influencing the stock status estimation (Lee et al., 2024), or require high computational time. This is the case of the integrated model Gadget (Globally applicable Area Disaggregated General Ecosystem Toolbox) (Begley and Howell, [[NoYear]]; Begley, 2004), which is an extensively used and a very flexible model that includes age–length dynamics [compared to others that only include length or age (Punt et al., 2020)]. However, Gadget could be very slow to run because it is not coded such that the gradients of the objective function with respect to the parameters are computed automatically, and it does not allow to include random effects on the parameters (Punt et al., 2020). It should be noted that other models extensively used, such as Stock synthesis, MULTIFAN, or CASAL, do not account random effects properly and that the ability to do so for several parameters is one of the essential features that the next generation of stock assessment models should have (Punt et al., 2020).

In contrast, estimating growth empirically, external or “outside” stock assessment models, offers greater flexibility in data selection, allowing researchers to incorporate representative age–length data and length–frequency distributions and to choose parameter values according to their expert knowledge on the species’ biology. In addition, estimation “outside” enables the use of more appropriate calculation methods, such as the use of random effects on the parameters through non-linear mixed models which have been proven to fit better when there is a high variability in the data (Pilling et al., 2002; Stewart et al., 2022). However, implementing these methods “inside” remains computationally intensive (Maunder and Punt, 2013). Despite these advantages, empirical growth estimation is subject to several challenges arising from uncertainties in age determination (Campana, 2001), sampling bias, and increased variability driven by environmental factors, including gear selectivity and oceanographic conditions (Pennino et al., 2020; Fernández-Corredor et al., 2021), among others (Lee et al., 2024).

European anchovy (Engraulis encrasicolus) inhabiting the Gulf of Cádiz (hereafter, GoC), ICES Subdivision 9.a, serves as a compelling case study in the context of estimating fish growth. This species is a small pelagic species (<20 cm) with a short life cycle (living up to approximately a maximum of 3 years), whose population fluctuations are governed by environmental drivers (particularly temperature, wind regimes, and the Guadalquivir river discharges), and therefore recruitment success depends highly on meteorological and oceanographic conditions during the early stages of their vital development (Ruiz et al., 2006; Ruiz et al., 2009; Rincón et al., 2016). The GoC is a highly productive coastal ecosystem that supports a significant anchovy fishery, yet growth patterns have fluctuated in recent decades, necessitating a better understanding of growth dynamics (Costalago et al., 2011). Research on anchovy growth in the GoC dates back to the 1990s, with early studies reporting life cycle characteristics such as length–frequency distributions, length–weight relationships, age structure, and reproduction (Rodríguez-Roda, 1977). Thereafter (Millán, 1999), described the reproductive biology of this species, highlighting variations in growth in response to environmental conditions. The first growth parameter estimates for GoC anchovy were derived from monthly length–frequency samples collected in three commercial ports in the area during 5 years, yielding estimates of   L ∞ =18.95 cm and k = 0.90 using ELEFAN I and   L ∞ =18.69 and k = 0.90 using Powell–Wheterall methods (Bellido et al., 2000). Currently, for the European anchovy inhabiting Atlantic waters, including the GoC, the ICES Working Group on southern horse mackerel, anchovy, and sardine (WGHANSA) provides scientific advice on fishing opportunities (ICES, 2022) using a Gadget stock assessment model since 2018. This model growth parameters estimates are as follows:   L ∞ =29.17 cm and k = 0.08 (ICES, 2018b).

To test model effectiveness in capturing the growth dynamics for short-lived small pelagic species, we propose a comparative analysis following recommended practices outlined by (Lee et al., 2024) for incorporating age–length data into a Gadget-integrated stock assessment model for the specific case of anchovy (Engraulis encrasicolus) stock in the Gulf of Cadiz. Our approach involves a three-step process: initially, conducting an exploratory analysis of the length-at-age datasets available to identify potential variability and biases; secondly, fitting a von Bertalanffy growth function in a random effects’ framework aligned with the data properties previously identified; and finally, comparing the outcomes of the existing model with two alternative Gadget implementations ( Figure 1 ). In one approach, we fix the growth parameters based on values derived from the previous step, while the other permits the model to autonomously calculate the parameters. For these alternative implementations, goodness of fit (likelihood scores), computational time, and their ability to accurately reflect the biology of the species were considered as performance indicators.

The overarching goal of this analysis is to determine the optimal model configuration in terms of both efficiency and accuracy, exploring the possibility to include random effects in growth estimation “outside” the integrated model by fixing the resulting parameters. This study seeks to determine the best approach for incorporating length–age data into integrated models when datasets exhibit high variability and limited age classes, particularly when random effects cannot be directly implemented “inside” the model. This could offer a valuable solution to improve the precision and effectiveness of stock assessments and, therefore, to inform adaptive fisheries management strategies.

The methods used here include the identification of potential variability and biases in the age–length datasets available, an exploratory analysis of non-linear models with fixed and random effects to fit a von Bertalanffy growth function “outside” the model, the selection of the most suitable methodology considering the data properties, the conditioning of two different Gadget models with the data available: one fixing the values for growth parameters using those estimated in the previous step and another allowing the model to estimate them, and lastly, the comparison of the outcomes of the model used to provide scientific advice in 2023 (ICES, 2023) with the two alternative Gadget implementations described in the previous step in terms of computational time and total likelihood score.

In this study, the guidelines proposed by (Lee et al., 2024) were applied to incorporate length-at-age data into a specific integrated fisheries stock assessment model (Maunder and Punt, 2013), demonstrating that the efficacy of the random effects approach for estimating growth in short-lived species was also characterized by high variability in age–length data, like the anchovy in the GoC, where existing growth parameter estimates were inaccurate and outdated.

Notably, for the anchovy in the Gulf of Cádiz, a comprehensive data and methodological revision was necessary due to the limitations of existing growth parameter estimates. The initial estimation for the growth parameters of this species, published in 2000 (Bellido et al., 2000), had been widely referenced. Since 2018, these parameters are also used as initial values in the current model for providing scientific advice on TACs (ICES, 2022). However, the estimated growth parameters in this model fail to accurately capture the species’ biology, highlighting a need for improvement. This limitation led to a re-examination of the approach to integrating length-age data in models, resulting in a new research question: how to effectively incorporate length-age data in integrated models when dealing with highly variable datasets and limited age classes.

In order to address this question, the guidelines proposed by (Lee et al., 2024) were applied to the anchovy in the GoC. A thorough analysis of the anchovy data available until 2023 was conducted, resulting in the inclusion of a new dataset focused on the recruited population. The data analysis revealed a great variability, and also with length classes overlapping for ages 0 and 1, where the minimum length is the same for both ages and the maximum length is only 20.5 cm higher for age 1 (see Table 2 and Figure 3 ), which could be one of the reasons explaining the model miss estimates for L ∞   and k when they are estimated by the model ( Table 4 ).

The observed variability motivated a dedicated methodology analysis, aiming to identify alternative approaches to address the issue. Another possible reason for the misestimates could be related to the fact that growth parameter calculation inside the model does not account for the variability and overlapping in the different age–length datasets, a problem documented by (Sainsbury, 2011). To address this, non-linear mixed-effects regression techniques were considered, as they have been used for observational studies that are replicated across sites or times. In this case, it was assumed that repeated measures taken over time represent lengths at various ages for the same individuals [this is analogous to the methodology proposed by (Helser and Lai, 2004)]. This allows for the incorporation of individual variability in a robust way by explicitly assuming that the growth parameters for each individual represent samples from a normal multivariate population of growth parameters. This approach has been successfully applied in various fields, including human growth (Rogol et al., 2000). As (Pilling et al., 2002) explains, the key advantage of this methodology is that it estimates the mean vector and the covariance matrix of the population of growth parameters, which are precisely the parameters needed to assess individual growth variability. Several authors (Cope and Punt, 2007; Pilling et al., 2002; Castillo-Jordán et al., 2010; Weisberg et al., 2010; Stewart et al., 2022) have documented the benefits of non-linear mixed-effects models, which facilitate the incorporation of different sources of growth variability by including fixed effects and random effects. These models are recommended when a great variation around the mean growth pattern is observed in most length-at-age plots ( Figure 3 ), and they have proven to overcome many of the shortcomings of the fixed-effects models approach (Wang and Ellis, 1998; Weisberg et al., 2010).

Nevertheless, the choice of methodology for fitting these models presents several challenges. A key issue is the need to calculate an integral that typically cannot be solved analytically. The accuracy of the methodology is highly dependent on how this integral is approximated or solved. For the purposes of this work, we implemented the approach proposed by (Lindstrom and Bates, 1990), which linearizes the non-linear model around the expected value of the random effects. This linearization simplifies the integral to a form that can be solved analytically, thus eliminating the need for numerical integration. However, a notable limitation of this method is that it may result in less accurate estimates, particularly when the random effects deviate from normality. Although the results of normality tests were inconclusive, model selection diagnostics—such as AIC, BIC, and standardized residuals ( Figure 4 )—indicate that mixed-effects models outperform non-linear approaches in this case. These diagnostics, which assess overfitting and underfitting, provide strong evidence supporting the superiority of mixed-effects models for this particular analysis.

Despite this, some misspecifications were observed for sizes exceeding 163 mm ( Figure 4 ). These inaccuracies could potentially be addressed by adopting more advanced methods, such as the Laplace approximation and Bayesian approaches, which would allow for the specification of non-normal distributions for the random effects. However, implementing these techniques would require further customization of the model to estimate random effects based on the sampling month, demanding additional expertise and incurring higher computational costs. A more comprehensive discussion of computational approaches for estimating random effects is available in (Thorson and Minto, 2015).

Additionally, it is important to note that these misspecifications could also stem from the limited sample size for larger age classes and the absence of measurement error in the aging process. As emphasized by (Cope and Punt, 2007), aging error should be taken into account when fitting growth curves. This error primarily arises from the inconsistency between otolith formation, which is subject to environmental influences (Izzo et al., 2018), and the interpretation of age (Campana and Thorrold, 2011; Morales-Nin, 2000; Campana, 2001). Another potential source of error is the decision regarding which parameters should include random effects, as this choice may be influenced by the normality assumption. Alternative approaches could yield different results. In this case, the AIC and BIC values underscore the importance of incorporating random effects for the t ₀ parameter. Notably, the two best models—those that fix two of the three parameters—as well as the best model that fixes only one parameter, all include random effects for t ₀ (see Table 3 ). This finding is likely linked to the high variability observed in the mean size of younger age groups ( Table 2 ). Accurate estimation of the t ₀ parameter is crucial, as misestimating it can introduce a significant bias into the estimates of k and L∞. Specifically, a negative covariance was observed between k and L∞, indicating that an overestimation of k (which suggests a faster growth rate) generally leads to an underestimation of L∞, implying a smaller potential maximum size. Conversely, underestimating k results in an overestimation of L∞ (Archontoulis and Miguez, 2015; ICES, 2019; Lentin and Elvarsson, 2021; Stewart et al., 2022). These relationships have important implications for the use of these parameters as proxies in life-history estimation and stock assessment models, highlighting the need for careful consideration of model structure and underlying assumptions.

In light of these findings, it is interesting to examine the parameter values estimated. The k value shows the most significant difference compared to previous attempts to estimate growth “outside” the stock assessment model. Most of the non-linear and mixed-effects models tested ( Table 3 ) provide values below 0.5. However, when t ₀ is fixed in the non-linear regression, the highest estimated value is obtained (k = 1.34). This value is closer to the k estimated by (Bellido et al., 2000), which was around 1 for most of the years studied (ranging from 0.70 to 1.07). This is particularly interesting because they also assumed t ₀ as fixed. Nevertheless, there are significant differences in the L ∞   value between both approaches. In this case, an unrealistically low L ∞   of 14.9 cm was obtained, compared to the L ∞   values between 18.3 and 19.69 cm obtained previously. Furthermore, in the non-linear framework, allowing t ₀ to be estimated results in a small k value of 0.1 but an unrealistically large L ∞     of 42.7 cm. These improbable estimates are obtained regardless of the choice of t₀ , and as mentioned before this is likely associated with the high variability in lengths for younger ages. This provides additional evidence that the non-linear approach without random effects is inadequate for this dataset, which predominantly consists of individuals aged 0 and 1 ( Table 2 ).

With the most suitable dataset and methodology for estimating growth parameters identified, a comparative analysis was conducted to test their incorporation into the Gadget integrated model. Specifically, the parameters obtained by the mixed-effects approach were first fixed into the model, and then the model was allowed to autonomously calculate the parameters. The efficiency of both Gadget implementations was evaluated in terms of computational time and goodness of fit. For both indicators, it was determined that the first implementation was the more efficient, reducing computational time in 1 h and improving likelihood score values ( Table 4 ). A notable time reduction was observed, as the number of parameters to estimate had been reduced, but 1 h is a considerable reduction for this stock assessment model, giving some idea of the importance of these parameters in the rest of the calculations and having in mind that this model is able to handle both age and length data at the same time. This implementation is also more accurate for the t ₀ parameter, as the estimated value (-1.29) is the closest to the mean of the values obtained using the mixed-effects approach (-1.14) compared to the other Gadget implementations. Since t ₀​ depends on the estimated recruitment length, this length should be correspondingly more accurate. Furthermore, the Gadget-estimated parameter values are close to the unrealistic ones obtained by the model used for the 2023 assessment, indicating that high variability was already present in the earlier dataset.

Moreover, the dataset used in this study (including ECOCADIZ-RECLUTAS and ECOCADIZ age 0) has important implications for the management of the anchovy fishery and currently is not the same dataset used to provide TAC advice. For a short-lived, small pelagic fish like anchovy, its sustainability relies on annual recruitment to ensure a healthy adult biomass for the following year. The assessment and management of these species are particularly challenging due to their short life expectancy, characteristic aggregative behavior, rapid response to climate and environmental signals, and large, variable natural mortality (Barange et al., 2009). Therefore, management is challenged by large recruitment variability (Thorson et al., 2014), which accounts for the largest source of population interannual variability (Siple et al., 2021). Hence, an early indication of recruitment greatly benefits their management (Uriarte et al., 2023). Advantages of incorporating recruitment indicators, like juvenile surveys, in stock assessment and management have been demonstrated for the Bay of Biscay anchovy stock (Boyra et al., 2013; Uriarte et al., 2023), where the advice for TAC recommendation is provided once the results of the recruitment survey are available. In particular, the Gulf of Cádiz anchovy fishery targets a highly variable resource with a short lifespan and strong dependence on yearly recruitment (Ruiz et al., 2006; Ruiz et al., 2009; ICES, 2023). Consequently, including the ECOCADIZ-RECLUTAS autumn acoustic survey dataset as an indicator of early recruitment in the assessment model will reduce the uncertainty in estimating the strength of the next cohort. This early recruitment information, combined with an accurate growth function, can serve as a tool to reduce uncertainty by enabling precise projections of the spawning stock biomass (SSB) for the following year. The SSB, in turn, directly informs recommendations for total allowable catch (TAC), facilitating more proactive and adaptive management measures. For example, fishery managers can implement precautionary adjustments to TACs or fishing effort levels before recruitment failures are fully manifested in the fishery. Such proactive measures help minimize both ecological and economic risks, ensuring more sustainable and resilient fisheries management.

In summary, since a choice needs to be made regarding whether growth should be estimated outside or inside the stock assessment model (Maunder et al., 2016), sufficient evidence was presented here for the GoC anchovy stock to decide to estimate it outside the model or what is called an “empirical approach” according to (Lee et al., 2024), who recommend to use it when there is a high biological variability. This approach has been shown to be effective for this stock, and its adoption holds the potential to be broadened to encompass other short-lived small pelagic species that exhibit similar data variability, offering a valuable solution for improving the accuracy and robustness of stock assessments to better understand population dynamics. Moreover, our approach addresses a significant gap in existing integrated models used to provide TAC scientific advice, which often do not have the “random effects on parameters” feature; by reducing uncertainty in growth parameters, this approach contributes to more stable and scientifically justified TAC recommendations, minimizing the risk of over- or underestimating the sustainable harvest levels, thereby providing a more reliable and robust framework for informing short-lived small pelagic fisheries management decisions.