Error Distribution Model to Standardize LPUE, CPUE and Survey-Derived Catch Rates of Target and Non-Target Species

Abstract: Indices of abundance are usually a key input parameter used for fitting a stock assessment model, as they provide abundance estimates representative of the fraction of the stock that is vulnerable to fishing. These indices can be estimated from catches derived from fishery-dependent sources, such as catch per unit effort (CPUE) and landings per unit effort (LPUE), or from scientific survey data (e.g., relative population number—RPN). However, fluctuations in many factors (e.g., vessel size, period, area, gear) may affect the catch rates, bringing the need to evaluate the appropriateness of the statistical models for the standardization process. In this research, we analyzed different generalized linear models to select the best technique to standardize catch rates of target and non-target species from fishery dependent (CPUE and LPUE) and independent (RPN) data. The examined error distribution models were gamma, lognormal, tweedie, and hurdle models. For hurdle, positive observations were analyzed assuming a lognormal (hurdle–lognormal) or gamma (hurdle–gamma) error distribution. Based on deviance table analyses and diagnostic checks, the hurdle–lognormal was the statistical model that best satisfied the underlying characteristics of the different data sets. Finally, catch rates (CPUE, LPUE and RPN) of the thornback ray Raja clavata, blackbelly rosefish Helicolenus dactylopterus, and common mora Mora moro from the NE Atlantic (Azores region) were standardized. The analyses confirmed the spatial and temporal nature of their distribution.


Introduction
Fisheries research is a subject of high interest with both economic and ecological relevance, mainly focused on guaranteeing the sustainability of the resource and the economic performance of the fishery. The integration of information on species' life history, fisheries monitoring, and resource surveys for assessing the stock size and harvest rate relative to sustainable reference points is known as stock assessment [1]. Its results often become scientific advice that is delivered to governmental agencies to be applied as law enforcement and fishing directives [2].
The most common method for stock assessment is the use of mathematical models that fit the available data to provide simplified representations of population and fishery dynamics [3]. Models used for stock assessment are usually based on several parameters, including mortality rates, reproductive aspects, size composition, and indices of relative abundance, to estimate current population status [1]. The type of model and analysis to be used for stock assessment relies on the species' available information and the data quality [4,5].
Indices of relative abundance are usually a key input parameter used for fitting a stock assessment model [1], and they can be estimated from catches derived from fishery-dependent sources, such as catch per unit effort (CPUE) and landings per unit effort (LPUE), or from
The LPUE (kg per landing per vessel) series was compiled from daily landing reports obtained from the auction service of the Azores-Lotaçor S.A. database during the period 1985-2017. Data recorded for each fishing operation included vessel identification, métier, and catch in kg by species. The CPUE (catch per day at sea) data came from the database of the Department of Oceanography and Fisheries, University of the Azores (DOP/UAz), and was collected between 1990 and 2017 as part of the European Union's Data Collection Framework (DCF; [27]). Standardized fishing inquiries (n = 31,616) were performed by clerks to the vessels' captains during the landings. Each inquiry included the vessel identification, the number of days at sea, and detailed information on the fishing operation, such as the gear type, mean depth of fishing, and catch in kg by species. Survey-derived abundance indices (RPN; individuals per 1000 hooks) were calculated for the period 1996-2019 and came from the Azorean spring bottom longline survey [10]. Surveys followed a stratified random design and covered the main islands and major seamounts. Each set included the area, moon phase, depth stratum, soak time, and catch number by species.
Factors considered in the analyses varied among the different data sets and are detailed in Table 1. Thornback ray Raja clavata (RJC), blackbelly rosefish Helicolenus dactylopterus (BRF), and common mora Mora moro (RIB) were the fish species selected as case studies. In addition to being species that typify the depth-aligned demersal fish assemblage structure in the Azores [28][29][30][31], both the thornback ray and common mora are usually caught as by-catch of the Azorean hook-and-line fisheries targeting demersal fishes such as blackbelly rosefish H. dactylopterus, blackspot seabream Pagellus bogaraveo, and alfonsinos Beryx spp. [29,30,32,33].
Records with structural zeros (see [34,35]) and missing effort data were excluded from the analyses.

Statistical Models
Standardized abundances were estimated by a Generalized Linear Modeling (GLM) approach. Several error distributions were examined, namely:

1.
Gamma. A positive constant was added to the nominal catch rates (RPN, CPUE and LPUE) in order to deal with zero catches. Two different values were tested: 1 and c = 10% of the mean catch rate [9,36]. The logarithm was used as the link function.

2.
Lognormal. This distribution assumed that the logarithm of the catch rate was normally distributed, using the identity as the link function. This error distribution also required the addition of a positive constant (1 or c) to deal with zero catches.

3.
Tweedie. This distribution is part of the exponential family of distributions and is defined by a mean (µ) and variance (ϕµ p ), in which ϕ is the dispersion parameter and p is an index parameter [22]. Given that this distribution can handle a certain proportion of zeros, the nominal catches were used directly. The power parameter (p) of the variance function was calculated by maximum likelihood estimation. 4.
Hurdle models. The catch estimates involved fitting two sub-models to the data [19,35].
The first sub-model modeled the probability that a positive observation (non-zero catch) occurred, assuming a binomial error distribution and logit link function. Positive observations were analyzed using a second sub-model assuming (1) a lognormal (hurdle-lognormal) error distribution with an identity link function and logtransformed catch rates, and (2) a gamma (hurdle-gamma) error distribution with a log link function.
The general formulations used in the present study were expressed by the following equations: Model-checking and fit diagnostics [16] were used to select the most appropriate error-model as follows: 1.
Pearson residuals were plotted against the fitted values as a check of the assumed variance function; 2.
Standardized deviance residuals were plotted against the estimated linear predictor (η) to check for systematic deviations from the assumptions underlying the error distribution; and 3.
The dependent variable was plotted against the estimated linear predictor (η) as a check of the assumed link function.
Plots were not examined for the model of positive observations' probability (binomial error distribution) because they do not provide particularly useful or interesting information [16].
Any observation for which Cook's distance was greater than 1 was considered highly influential and therefore removed from the analysis.

Standardization Procedure
Deviance tables were used to select the explanatory factors and interactions that explained most of the variability in the data [21]. The effect of each explanatory factor/interaction was evaluated according to (1) the percent of deviance explained by the addition of a specific variable to the model, and (2) the significance of an additional variable to the total deviance explained (Chi-squared test, α = 0.05). The factors and interactions selected for the final model were chosen when their insertion accounted for equal or greater than 5% of the total deviance, and the significance of their choice was confirmed by the Chi-squared test.
After selecting the set of explanatory factors and interactions, all interactions that included the year factor were treated as random [37]. This process converted the basic models from GLM into generalized linear mixed models (GLMMs). The significance of the random interactions was evaluated using the likelihood ratio test [38], the Akaike Information Criteria (AIC), and Bayesian Information Criteria (BIC), where lower values indicated better model fitting. Once a final model was identified, diagnostic plots were revised to identify potential departures from the GLMM assumptions or observations with a large influence on the model results.
Standardized abundance indices were estimated as the least-squares means (LSmeans) of the year factor for the selected model. In the case of the selected model being a hurdle model, indices of abundance were estimated as the product of the LSmeans from each of the two analyses, after back-transforming to the response scale. Then, the variance estimation of the standardized index was calculated following Walter and Ortiz [39] for two-stage variance estimates.
The standardization procedure was then repeated for each species using the alternative error distributions considered at the beginning of the analysis. The standardized indices were summarized by the percentage differences from those from the most appropriate error-model.

Catch Trend Comparison between Datasets
To compare the general LPUE, CPUE, and RPN trends, linear trend models were applied to the annual standardized catch rates by species to find the intercept (a) and slope (b) that give the best average fit. Analysis of covariance was used to test all the regression lines to see which ones have significantly different (p < 0.05) a and b values. Differences in a were interpreted as differences in magnitude, while differences in b were interpreted as differences in the rate of change [40,41].

Nominal Catch Data
Catch distributions were found to be highly skewed, with many zero or low catches, mostly for non-target species, as expected ( Figure 1). The LPUE data comprised 4952 landing reports, and such a high proportion of zero catches was observed throughout the studied years for R. clavata, H. dactylopterus, and M. moro at 0.68, 0.62, and 0.88, respectively. For the CPUE database, a total of 31,616 inquiries were analyzed for the three species, with a relative number of zero catches of 0.86 (R.

Nominal Catch Data
Catch distributions were found to be highly skewed, with many zero or low catches, mostly for non-target species, as expected (

Error-Model Selection (Application)
The Pearson residuals were plotted against the fitted values ( Figure S1) as a check of the adequacy of the assumed variance function. The null (expected) pattern of this plot is a distribution of residuals with no trend [16]. The positive-lognormal error distribution was the best model where there was no trend for the residuals of the LPUE, CPUE, and

Error-Model Selection (Application)
The Pearson residuals were plotted against the fitted values ( Figure S1) as a check of the adequacy of the assumed variance function. The null (expected) pattern of this plot is a distribution of residuals with no trend [16]. The positive-lognormal error distribution was the best model where there was no trend for the residuals of the LPUE, CPUE, and RPN data for all species ( Figure S1). The lognormal (1 and c), gamma (1 and c), tweedie, and positive-gamma models showed discrepancies with residuals concentrated at specific plot regions or a few points far from the rest, which indicated the variance function was not appropriate. A positive trend (slope > 0) implies that the assumed variance function is rising too slowly concerning the mean, whereas a negative trend (slope < 0) implies the opposite [16]. Therefore, these error distributions would either underestimate or overstate catch rates for the observed CPUE, LPUE, and RPN.
Diagnostic plots of standardized deviance residuals against the estimated linear predictor are shown in Figure S2. The null pattern is the distribution of the residuals with a mean of zero (dashed line) and a constant range [16]. For all three databases, the gamma (1 and c) and lognormal (1 and c) error distributions showed declining trends in the mean, and the tweedie model showed a systematic change in range with the linear predictor. Curvature can be caused by a variety of factors, including incorrect link function selection, incorrect scale selection for one or more covariates, or the absence of a quadratic component in a covariate [16]. The positive-gamma and positive-lognormal models were consistent with the expected pattern of a mean close to zero.
The link function selected for each error distribution was checked through the plot of the dependent variable (LPUE, CPUE, or RPN) against the linear predictor ( Figure S3). The null pattern is a straight line [16], which was observed for all species, databases, and error distributions ( Figure S3). This means that the link function was properly chosen for the error distributions analyzed.
To summarize, the model diagnostic plots indicated that the positive-lognormal model was the most appropriate error distribution for fitting the LPUE, CPUE, and RPN data for R. clavata, H. dactylopterus, and M. moro. Consequently, the hurdle-lognormal model was selected as the most suitable error distribution for standardizing such fishery-dependent and independent data for target and non-target species.

Standardization Procedure
Deviance tables produced from the models and used for selecting relevant explanatory factors and interactions that better explained the data variability for LPUE, CPUE, and RPN data are shown in Tables S1-S3, respectively. The interactions that included the year factor were treated as random interactions, and the final best-fitted model was selected according to the lowest log-likelihood, AIC, and BIC values of binomial and positive-lognormal models of each species (Table S4).
Diagnostic plots of each final model selected from the positive-lognormal error distribution ( Figure S4) displayed no strong indication of departure from the null pattern, suggesting that the model selected was well fitted. Subsequently, the standardization procedure was performed using the hurdle-lognormal method (Figure 2). Nominal and standardized abundance indices estimated per year for each species and database are shown in Table S5.
The null pattern is a straight line [16], which was observed for all species, databases, and error distributions ( Figure S3). This means that the link function was properly chosen for the error distributions analyzed.
To summarize, the model diagnostic plots indicated that the positive-lognormal model was the most appropriate error distribution for fitting the LPUE, CPUE, and RPN data for R. clavata, H. dactylopterus, and M. moro. Consequently, the hurdle-lognormal model was selected as the most suitable error distribution for standardizing such fisherydependent and independent data for target and non-target species.

Standardization Procedure
Deviance tables produced from the models and used for selecting relevant explanatory factors and interactions that better explained the data variability for LPUE, CPUE, and RPN data are shown in Tables S1-S3, respectively. The interactions that included the year factor were treated as random interactions, and the final best-fitted model was selected according to the lowest log-likelihood, AIC, and BIC values of binomial and positive-lognormal models of each species (Table S4).
Diagnostic plots of each final model selected from the positive-lognormal error distribution ( Figure S4) displayed no strong indication of departure from the null pattern, suggesting that the model selected was well fitted. Subsequently, the standardization procedure was performed using the hurdle-lognormal method (Figure 2). Nominal and standardized abundance indices estimated per year for each species and database are shown in Table S5.

Consequences of Choosing a Wrong Error-Model
When knowledge of the true status of each species is unknown, an alternative approach for assessing the implications of basing standardization on the incorrect error distribution, in terms of the magnitude of error in the standardized catch rates, is to assume that the results from the selected model provide the most accurate representation of true stock status [21]. Therefore, although the hurdle-lognormal was selected as the best model for standardizing catch rates of all species and databases analyzed, the other explored error distributions were processed with the same method for the model selection, the best model formulation (Table S6), and the diagnostic plots ( Figure S5). Standardized catch rates (Table S7) were then summarized as the percent difference from the hurdle-lognormal model.
For the LPUE database (Figure 3), the hurdle-gamma model showed lower percentage differences from the chosen error distribution for the R. clavata and M. moro, with respective average differences of 6% and −7%, while for the H. dactylopterus, the tweedie model showed the lowest ones (−9% on average). The other error-models for R. clavata showed greater discrepancies from the hurdle-lognormal, varying from 23% (tweedie) to 59% (gamma + 1) on average. For H. dactylopterus, the average percentage differences from those from the hurdle-lognormal model varied from −14% (hurdle-gamma) to 70% (lognormal + 1). In the case of M. moro, these average differences in magnitude were greater than 130 times.
with respective average differences of 6% and −7%, while for the H. dactylopterus, the tweedie model showed the lowest ones (−9% on average). The other error-models for R. clavata showed greater discrepancies from the hurdle-lognormal, varying from 23% (tweedie) to 59% (gamma + 1) on average. For H. dactylopterus, the average percentage differences from those from the hurdle-lognormal model varied from −14% (hurdlegamma) to 70% (lognormal + 1). In the case of M. moro, these average differences in magnitude were greater than 130 times.  (Figure 3) for R. clavata from the gamma (c and 1), lognormal (c and 1), and hurdlegamma models were similar, with percentage differences ranging from −1% (hurdle-gamma) to 2% (gamma + c). The estimates from the tweedie error distribution showed greater variations, with an average difference of −6%. For H. dactylopterus, percentage differences from those estimates of standardized CPUE for the hurdle-lognormal model were less than 4% on average. However, these differences could reach about 127% in some years. For the M. moro, the hurdlegamma model showed the lower percentage difference from the hurdle-lognormal, with an average of 6%. The other error-models had greater discrepancies, ranging from 32% (lognormal + c) to 35% (tweedie and gamma + c).  Figure 3) for R. clavata from the gamma (c and 1), lognormal (c and 1), and hurdle-gamma models were similar, with percentage differences ranging from −1% (hurdle-gamma) to 2% (gamma + c). The estimates from the tweedie error distribution showed greater variations, with an average difference of −6%. For H. dactylopterus, percentage differences from those estimates of standardized CPUE for the hurdle-lognormal model were less than 4% on average. However, these differences could reach about 127% in some years. For the M. moro, the hurdle-gamma model showed the lower percentage difference from the hurdle-lognormal, with an average of 6%. The other error-models had greater discrepancies, ranging from 32% (lognormal + c) to 35% (tweedie and gamma + c).
The results of the RPN estimates from the hurdle-gamma models were very similar to those from the hurdle-lognormal for all species, with average differences ranging from −1% (R. clavata and H. dactylopterus) to 0% (M. moro; Figure 3). The other errormodels for R. clavata showed percentage differences varying from −2% (lognormal + 1) to 28% (gamma + 1) on average. For H. dactylopterus, the RPN estimates from the other models were similar in trend among them but different in magnitude from those from the hurdle-lognormal, with average percentage differences ranging from 22% (tweedie) to 47% (gamma + 1). For M. moro, these percentage differences were greater than 80%.

Catch Trend Comparison between Datasets
Standardized catch rates over the studied years showed varying trends according to each species and across different databases (Figure 2). For R. clavata, both the LPUE and RPN demonstrated an increasing trend (Figure 2), with no significant differences between their intercepts and slopes (Table 2). These trends, in turn, were significantly different from that of the CPUE (Table 2), which was greater in magnitude (intercept) but showed a decreasing tendency (Figure 2). For H. dactylopterus, LPUE, CPUE and RPN were decreasing in trend (Figure 2), with no significant difference between their intercepts and slopes (Table 2). Trend analysis for M. moro indicated significant differences between the intercepts and slopes from the CPUE and RPN compared to that from the LPUE (Table 2), which was lower in magnitude and showed a positive rate of change (i.e., an increasing trend; Figure 2). Table 2. Results of analysis of covariance (ANCOVA) for differences between regression lines from annual standardized LPUE (kg landing −1 vessel −1 ), CPUE (kg days at sea −1 vessel −1 ), and RPN (ind. 10 −3 hooks) catch rates by species (RJC: thornback ray Raja clavata, BRF: blackbelly rosefish Helicolenus dactylopterus, RIB: common mora Mora moro). Significant differences at the 0.05 significance level are in bold.

Final Considerations
This study updated and expanded the protocol developed by Ortiz and Arocha [21] for comparing error distribution models to standardize catch rates, aiming to select the best error-model and evaluate model assumptions given the data set characteristics. The use of generalized linear mixed models, particularly the hurdle-lognormal model, resulted in more agreement between model assumptions and observed catch rates. The results confirmed the previously reported consistency of these models with zero-inflated data sets [20,21,48], the importance of model checking and validation procedures for the standardization of nominal catch rates, and the implications of basing the standardization on the incorrect error distribution model [21].
The analyses of catch rates from fishery-dependent and independent data indicated that a few past or recent years could be responsible for the observed differences in the species' catch trends. For R. clavata, the standardized CPUE values observed in the last three years (2015-2017) were the main responsible for pulling the trend down. This species is primarily found in shallower waters (down to 250 m) around the islands [33,49] and those lower CPUE values are assumed to be related to the observed reduction in both the catches of the species [30] and the fishing effort (days at sea) of the Azorean fleet in coastal areas during the last few years [50] due to the fishing area restrictions [30]. At the same time, as catches declined, a rise in LPUE during the 2015-2017 period might suggest that fewer vessels landed this species and that each of these landing occurrences recorded large captures. This could be explained by the fact that the vessels that continue to operate legally with hook-and-line fishing gear in areas where R. clavata is more likely to be found are typically small (less than 12 m in length), with limited storage capacity and autonomy. This was in part confirmed by this study as métier, gear, vessel size and area, and the interaction between the last three factors and year were significant explanatory variables when modeling catch rates.
For H. dactylopterus, the LPUE, CPUE and RPN were similar in trends. The Azorean spring bottom longline survey is designed to target demersal species such as H. dactylopterus by using a similar fishing strategy to that of the commercial fleet [10]. Consequently, it was expected that after the standardization process, both the fishery-dependent and independent catch rates of H. dactylopterus would show similar trends. This species is usually caught by vessels operating with hook-and-line gears [29], and the harvested population is mainly found in seamount areas at 350-800 m depth [31]. This was validated by this study, which found that métier, gear, depth, and area were important explanatory factors. Finally, the observed declining trends in abundance indices, combined with the sedentary nature of this species, confirmed the already reported vulnerability of this stock to fishing activity [31,51]. As a result, further assessment analyses and effective management strategies are required to establish a balance between stock exploitation and conservation.
In the case of the M. moro species, the observed variations in catch patterns were due to lower LPUE values until the early 1990s and the greatest one in 2017. The first could be associated with the under-exploited phase of its fishery. It is a deep-water species, mainly found at depths ranging from 300 to 2500 m on the outer continental shelf and slope [50]. The exploitation phase of this resource began in the early 2000s; before that, it was captured as a by-catch of the demersal fishery occurring in depths between 200 and 600 m [50]. With the reorientation of the larger vessels to offshore deep waters due to the implementation of fishing area restrictions, the M. moro came to dominate the bottom longline catches in depths between 600 and 1200 m. This was consistent with the analysis of catch rates in the present study because métier, vessel size, target effect, depth zone, and area were significant explanatory factors. Finally, the highest LPUE value in 2017 can represent an error in nominal data and needs to be further investigated.
The primary objective of standardizing catch rates is to generate a relative abundance index. This index is thought to be proportional to population size (e.g., [1,4,5,13,14]). As a result, the abundance index may be utilized as a basis for management measures (e.g., definition of total allowable catches-TACs) either directly, by analyzing the population trend (i.e., rise, decrease, stability) over time, or indirectly, as an input into stock assessment models (e.g., [1][2][3][4][5][6][7][8][9]). Therefore, the definition of the best suitable statistical GLM technique to standardize catch rates made by the present study constitutes a useful outcome for fisheries research and management. However, it is critical to understand the constraints that current data can impose on model application and to plan for the future by selecting what data should be gathered in the future. For example, one of the most important factors in standardizing catch data from commercial fisheries is where the fish were caught. This was evident in this study as area was a significant factor in modeling catch rates from scientific surveys for the three studied species. Therefore, regardless of the method used to standardize the catch rate, it will not perform well unless the appropriate explanatory variables are available.

Conclusions
The use of the hurdle-lognormal model outperformed other methods in the different data sources and, therefore, it is recommended as an effective method for standardizing catch rates when it is not possible to perform a comparison with other error-models. Further simulation-based analyses are advised to explore, for example, the effect of other potential sources of bias in the estimation of abundance indices (e.g., fishing area; market demand expressed as reference price, i.e., price per kg; oceanographic factors, e.g., oceanic and atmospheric phenomena, etc.) or the use of clustering-based area stratification for catch and effort data (e.g., [52]). Furthermore, while GLMs are the most commonly used models in fisheries stock assessment and abundance standardization, there have also been analyses that demonstrate the efficacy of using generalized additive models (GAMs; [53]) and, therefore, they should be explored in future studies.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/modelling3010001/s1, Figure S1. Diagnostic plot to detect the adequacy of the assumed variance function in generalized linear models: Pearson residuals plotted versus fitted values for thornback ray Raja clavata (RJC), blackbelly rosefish Helicolenus dactylopterus (BRF), and common mora Mora moro (RIB). The null pattern is a no trend in the residuals. The red line is the loess smoother through the plotted values. A: LPUE (kg landing −1 vessel −1 ); B: CPUE (kg days at sea −1 vessel −1 ); C: RPN (ind. 10 −3 hooks). Figure S2. Diagnostic plot to check the assumed error distribution in generalized linear models: standardized deviance residuals plotted versus the estimated linear predictor (η) for thornback ray Raja clavata (RJC), blackbelly rosefish Helicolenus dactylopterus (BRF), and common mora Mora moro (RIB). The null pattern is a distribution of residuals with mean zero and constant variance. The red line is the loess smoother through the plotted values. A: LPUE (kg landing −1 vessel −1 ); B: CPUE (kg days at sea −1 vessel −1 ); C: RPN (ind. 10 −3 hooks). Figure S3. Diagnostic plot to check the link function selection in generalized linear models: dependent variable plotted versus the estimated linear predictor (η) for thornback ray Raja clavata (RJC), blackbelly rosefish Helicolenus dactylopterus (BRF), and common mora Mora moro (RIB). The null pattern is a straight line. The red line is the loess smoother through the plotted values. A: LPUE (kg landing −1 vessel −1 ); B: CPUE (kg days at sea −1 vessel −1 ); C: RPN (ind. 10 −3 hooks). Figure S4. Diagnostic plots for positive catch rates of thornback ray Raja clavata (RJC), blackbelly rosefish Helicolenus dactylopterus (BRF), and common mora Mora moro (RIB) to check (I) the adequacy of the assumed variance function, (II) the assumed error distribution, and (III) the link function selection in the best selected models. The null pattern is a no trend in the residuals (I), a distribution of residuals with mean zero and constant variance (II), and a straight line ( Institutional Review Board Statement: No ethical approval was required from the Portuguese official veterinary department, and this study was performed in accordance with relevant institutional and national guidelines and regulations.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data underlying this article will be shared upon reasonable request to the corresponding author.