Trawlers from the Barcelona port operate on the continental shelf and slope fishing grounds (50-1000 m depth) located between 1°50’ and 2°50’ longitude east and 40°50’ and 41°30’ latitude north (Sardà et al. 1997). The fleet operates on a daily basis (with mandatory exit from port after 6 am and return to port before 6 pm) and a limited license system whereby total effort in the area has been frozen since 1986. New boats are only permitted if an existing boat is decommissioned. In addition to effort control, the only other measure of control is limiting mesh sizes (minimum 40 mm cod-end stretch mesh size) and neither in the study area nor in any other Mediterranean areas is there output control. Fish is auctioned daily at the premises of the port fish market and all transactions are recorded electronically for statistical purposes by the Barcelona Fishers’ Association.

We used the daily sale slips containing all transactions of red shrimp (A. antennatus) over the period 1994-2008 to calculate the total monthly landings (kg month^–1, lands), the effort measured as total number of trips performed monthly by each vessel (number of trips per month, trips), and the monthly average ex-vessel shrimp price (€ kg^–1, shprice). The same data set contained information on the engine power (horse power, hp) and gross registered tonnage (tons, grt), which we used as boat capacity indicators. As an indicator of fishing costs we used the average monthly fuel price (10^–3 € L^–1, fprice) from the EUROSTAT website: http://ec.europa.eu/energy/observatory/oil/bulletin_en.htm. As the environmental driver we used the NAO index, taken from the website of the Climate Analysis Group of the University of Exeter (UK):http://www1.secam.ex.ac.uk/cat/NAO. We considered: nao1, nao2 and nao3, corresponding to one, two and three years before the year of observed landings, using the significant time lags reported in Maynou (2008).

The response variable, lpue, was estimated as

$$lpue=\frac{lands}{trips}$$

hence, lpue is the monthly average landings of one vessel in each trip, corresponding to one fishing day (kg boat^–1 day^–1), which will form the basis for providing a standardized abundance index.

To assess the individual effect of each vessel, a numeric variable code was assigned to each of the 21 vessels in the fleet. Each observation was attributed to a sequential time variable from January 1994 to December 2008 (180 months) and a month variable describing the month of the year. Two more variables were derived a posteriori from code and month, after checking for statistical differences among their categories in the model, and then by performing the Tukey Honest Significant Differences test (TukeyHSD). We thus derived the new variables, group and period, for code and month, respectively. The group variable combines the 21 trawlers into three groups of increasing lpue and period is a binary variable which classifies the months in “high effect” (period1: all months excluding June and November) and “low effect” (period2: June and November). All variables are shown in Table 1.

We used generalized additive models, GAM, as described by Hastie and Tibshirani (1990)

G(y) = xβ+ s_j(x_j)+ε

where G(∙) is the link function, y is the response variable, x is the vector of linear predictors (explanatory variables), β is the corresponding vector of parameters, x_j are scalar predictors with unknown nonlinear effects represented by the functions s_j(∙), and ε is the random error.

The model building process consists of the following steps: 1) selection of the underlying distribution of the response (see Section 2.2.2 for more details); 2) selection of predictors building independent models for each covariant deleting insignificant effects in the final model; 3) selection between correlated predictors through the Pearson correlation coefficient (threshold value: ρ=|0.6|) to avoid problems of collinearity (Brauner and Shacham 1998) using the covariant with the most explanatory potential; and 4) analysis of residuals diagnostics.

All analyses were performed in R3.0.1 (mgcv-Rpackage: Wood 2006). The generalized cross validation (GCV: Craven and Wahba 1979) and the outer Newton iteration procedure were used to estimate model parameters. GCV is preferable to the UBRE/AIC method in the case of unknown smoothing parameters λ (Wood 2006). Second order P-spline as defined by Marx and Eilers (1998) was used as a smoother for nonlinear functions.

Landings per unit effort are usually modelled following Gaussian or gamma distribution functions, often without formal justifications (Stefánsson 1996). Here, we assigned a theoretical probability function to the lpue, using nonparametric techniques described by Wassermann (2005). The density function, f(y_n), was estimated using a histogram and the Gaussian kernel, where the nrd described by Scott 1992 was the rule-of-thumb used to select the bandwidth.

The empirical cumulative distribution function, F(y_n), and the 95% lower, L(y_n) and upper, U(y_n), confidence intervals were calculated as follows:

F(y_n) = [rank(y_n) – 0.5] / N

L(y_n) = max{F(y_n) – ε_n , 0}

U(y_n) = min{F(y_n) + ε_n , 1},

where rank(y_n) is the ranked vector of observations and ε_n is the associated vector of errors resulting from the DKW (Dvoretzky-Kiefer-Wolfowitz) inequality

$${\epsilon }_n=\sqrt{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2n$}\right.{\log }_e\left(\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$\alpha $}\right.\right)}$$

The hypothesis tested is H₀: F(y_n)=F₀(y_n) versus the alternative hypothesis H₁: F(y_n)≠F₀(y_n), where F₀(y_n) is a theoretical function belonging to the exponential family, especially F₀(y_n)=Ga(a,b), the gamma distribution, whose density function is

p(y) ∝ y^(a–1)e ^(–by)

for y>0. Parameters a and b are derived from the estimated expectation, E(y)=a/b, and variance Var(y)= a/b². Then the resulting distribution functions were graphically compared.

Selection criteria and explained deviance

Both AIC (Akaike 1973) and percentage of explained deviance (ED) were used as selection criteria: the selected model presented both the lowest AIC and the highest ED and all term parameters significantly different from zero.

The ED for each variable was also calculated in order to determine their relative importance in the final model. The residual deviance of the full model and the deviances of reduced models (i.e. the model excluding variable x_i) were calculated to derive the proportion explained by variable x_i:

DEx_i = [D(reduced model)–D(full model)]/D(null model)

where D(∙) is the deviance for a given model and in the reduced model variable x_i is omitted.

LPUE standardization

The model used for standardization was built in order to avoid dependency on fleet variables, maintaining environmental variables, which are expected to be related to the natural abundance of the species. The standardized LPUE is then

where lpue is the “nominal” or “raw” LPUE defined in Equation 1, E(lpue) is the unconditional expectation of the LPUE and E(lpue|_{Θ=θ, Λ=λ}) is the conditional expectation of the LPUE given the vectors of parameters θ and λ, estimated using the appropriate standardization model.

Finally we compared our standardized LPUE with an alternative abundance index derived from fisheries-independent data, available in the technical report SGMED-12-11 (Cardinale et al. 2012, pp. 136-150) after normalization of both variables.

Figure 1 provides the raw lpue time series as reckoned in Equation 1 with its spline estimation (upper plot) and annual average of the NAO index (lower plot). The lpue series, initially constant, started to decline at the end of 1998 with a sharp maximum low in the period 1999-2000. Then the trend changed to inter-annual variation. Conversely, the NAO index shows the lowest records in the period 1995-1996. The lpue begins to decrease after three years from the first year of negative NAO. In Figure 2 the time series of the fuel price shows an increasing trend during the observed period (upper panel). The total number of trips per month performed by the whole fleet declines during the same period (middle panel), but it declines only at the beginning of the fuel price rise, then remaining almost constant (lower panel). Low lpue in the period 2000-2001 is also related to the peak of fprice in the same period (compare upper panels of Figs 1 and 2).

Characteristics of lpue are plotted in Figure 3. The probability density function (pdf) is positively skewed (upper panels). Data hold atypical values in the right tail (see the box-plot, left middle panel) and the distribution function of the gamma distribution lies approximately inside the 95% confidence intervals of the empirical cumulative distribution function (ecdf) of lpue (right middle panel). Finally, the QQ-plots for the gamma and the Gaussian distributions provide evidence of a better fit of data to the gamma rather than the Gaussian distribution (on the left and the right lower panels respectively).

The model building process is summarized in Table 2. All linear terms are reported before all smooth functions. Models are sorted in ascending order of total ED, with the exception of models 11 and 12, corresponding to the GLM version of the final model 7 and the standardization model respectively. In comparison with the GLM (model 11), GAM’s AIC decreases and the ED increases substantially (1.56 times), though the number of degrees of freedom increases considerably (df_GAM=27.7 versus df_GLM=9). Variables grt and time were selected as better predictors than hp and fprice, respectively, according to the Pearson correlation coefficients (grt-hp was ρ=0.61; time-fprice ρ=0.69) and larger ED of former variables. All other correlations were ρ<|0.6|.

Table 2 incorporates the NAO terms, though their explanatory potential in the full model is too weak to be significant. Adding NAO terms to model 7 (Table 2: models 8, 9, 10) did not improve model fit in terms of ED and AIC. Nevertheless, we found significant NAO effect in less complex models.

The final model for A. antennatus LPUE is

where lpue~Ga(a,b), log_e is the link function, θ is the intercept, α is the parameter associated with the linear effect of grt, and β and γ_k are the effects associated with the categorical variables period and group. s_j(∙) are the smooth functions associated with trips, time and shprice, while ε represents the random errors of the model.

Figure 4 shows the diagnostic plots for the selected model. These plots show the results for the best model (corresponding to Equation 3 in the text and model 7 in Table 2). Residual quantiles lie on the straight line of the theoretical quantiles, although slightly heavy-tailed; in the histogram, residuals are consistent with normality and the relationship between response and fitted values is linear and positive. Residuals versus the linear predictor (that is, the sum of all partial effects) show a faint heteroscedasticity.

Table 3 shows results related to (a) the linear part, (b) the smooth functions and (c) the global parameters of the final model (Eq. 3), with a total explained deviance of 43.00%. The predictor with the highest explanatory potential was time (ED=12.40%), which captured the intra-annual fluctuations in red shrimp, lpue. The second predictor in terms of explained deviance was trips (ED=11.38%). Red shrimp price was the third most important predictor (ED=9.30%). Other variables such as grt, group and period had less impact. The model returned all significant parameters (p<<0.001). All partial effects are reported in Figure 5. The partial effect for time shows a substantial difference before and after the period 1999-2000, when a clear drop is exhibited in the shape. Before this threshold the shape is almost constant, whereas after the abrupt decay increasing variation is observed over recent years. Predictor trips represents a positive and monotonic relationship, reaching a plateau for trips>15; shprice reached its maximum value around 30 € kg^–1; lpue was significantly lower for period2, representing June and November, than for the rest of the year (period1). There were significant differences between the three groups of vessels. Variable grt showed a positive linear effect, meaning that larger vessels had higher lpue.

The model used to standardize A. antennatus LPUE is

where lpue~Ga (a,b), θ is the intercept, a the parameter for grt and γ_j the effects associated with code. s(trips) is the smooth function for trips and ε the random errors. With the aim of standardizing, this model comprises all fleet-dependent variables, whose effect must be removed from the nominal value. The diagnostic plots show reasonably good outputs (not shown).

Finally, LPUE index and the SGMED in 2002-2008 with their confidence intervals after normalization of variables are plotted in Figure 6. The normalized Barcelona fishery’s LPUE, calculated using Equations 4 and 2, has narrower confidence intervals than the normalized SGMED index. Only in 2004 were the two indices statistically different and in 2005 index estimates were at the limit of the other CI. Thus, five out of seven estimates can be considered equivalent.

We incorporated effort, temporal, economic and environmental variables into a global regression model to evaluate their relative importance. Model 7 captures LPUE variability with a total deviance of 43% explained by six predictors. In order to quantify the different sources of LPUE variability, we found that the set of fishery-related variables (trips, grt and groups) was the most important source, with an ED of 20.58%, followed by temporal (time and period, ED=13.12%) and finally economic variables (shprice, ED=9.30%). Among the variables taken from fishery data, trips has provided the greatest impact to date. A low number of trips per month could be associated with generalist trawlers operating usually on the continental shelf and with less knowledge about red shrimp fishing grounds than trawlers specialized in this fishery (Maynou et al. 2003), suggesting that the higher activity of non-specialized trawlers in deep-water harbours yields lower LPUE values (see Fig. 5: partial effect for trips). Conversely, boats performing a higher number of trips per month are expert in deep-water fishing grounds and their skippers are more likely to find high-concentration shoals through a process of trial and error (as hypothesized by Sardà and Maynou 1998). Boat characteristics (grt) also influence LPUE, as was to be expected. The greater the gross registered tonnage, the higher the LPUE that is observed. The variable group captures other capabilities of fishermen and technical characteristics of vessels, such as the type of engine, net shape and skipper’s expertise, which have been shown to positively bias the LPUE (Maynou et al. 2003, Marriott et al. 2011) and are expected to be important in many fisheries (Maunder and Punt 2004). The importance of these predictors implies that their influence should be eliminated during standardization.

Inter-annual variable (time) is much more important than intra-annual variable (period) (Table 3). The former is more strongly determined by a range of sources, such as environmental and economic drivers, than the latter. The same order of importance was found by Maynou et al. (2003). The LPUE was almost constant before 1999-2000, when a sharp decline was observed. We hypothesized a relationship with low NAO in the previous three years (see Fig. 1 and the introduction) that confirms the findings of Maynou (2008) and a possible relationship with the increase in fuel prices beginning in 2000. The following years were characterized by high inter-annual variability, when the price of fuel increased and showed greater variations. We believe that the trend could be explained by economic factors (Fig. 2), especially since very high fuel-related costs are incurred in the fishing of red shrimp as a result of it being performed in deep-water (Sardà et al. 1997).

The partial effect of ex-vessel prices, shprice, shows a parabolic-like shape and significant explanatory potential (ED=9.30%). Low selling prices do not induce fishermen to practice this deep-water fishery, because they could not offset the high associated costs, and trawlers would rather switch to continental shelf fisheries, with lower costs and lower risk. When there are profits to be made, probably owing to the low availability of the product, fishermen practice this deep-water fishery more intensely and landings per unit effort also increase. At the higher sales price bracket (i.e. more than about 30 € kg^–1) decreasing landing rates mean higher sale prices. Here, an alternation of cause and effect between the two variables probably comes into play. As mentioned, fprice is also an important explanatory variable, although it was not inserted in the final model because of its correlation with time. Fuel price has a significant effect on LPUE and can reach approximately 1.6% of explained deviance (as can be deduced from the model building process in Table 2).

Obtaining information on deep-sea species population dynamics is notoriously difficult, but our analysis suggests that the peculiarity of red shrimp fishery makes it possible to use fishery-dependent data to accurately describe the relative abundance of this resource. There are no discards for this fishery and the by-catch fraction, represented for example by Merluccius merluccius and Micromesistius poutassou, is small. These characteristics enable landings to be considered equivalent to catches and interchange LPUE and CPUE as indices (Hilborn and Walters 1992, Denis et al. 2002).

In turn, the definition of the relative importance of explanatory variables enables their impact on LPUE to be understood and makes intervention on the relevant variables possible from a management perspective. Fishery-related variables tend to have a significant effect on LPUE (ED=21% in our case), and management measures aiming to reduce fishing mortality in this heavily harvested stock (Cardinale et al. 2012) could be based on limiting the size of the trawlers. Furthermore, the number of trips permitted in deep-water fishing grounds in this fishery could be limited, for example, by defining a threshold when the number of trips does not significantly increase the partial effect for LPUE (see Fig. 5).

In addition, to evaluate the impact of predictors on the LPUE, the regression analysis could be the basis to provide a standardized index for assessing species stocks. Standardization of landings data allows an index of the real species abundance to be developed, assuming that the explanatory variables available remove (or explain) most of the variation in the data that is not attributable to natural changes (Maunder and Punt 2004). We advocate the selection of effort predictors for standardization.

Trends in CPUE (and LPUE) are usually assumed to reflect changes in the abundance of marine stocks (Maunder and Punt 2004), but the raw index is often not proportional to abundance (Maunder et al. 2006). The raw LPUE is in fact dependent on many human factors that could be avoided. Time variables have a strong relationship with the abundance and environmental factors, which in turn are related to abundance, so they cannot be used in the model during standardization (e.g. see Maunder and Punt 2004). The economic source of variability should be considered, but shrimp price and LPUE realistically have a cause-effect relationship, so they could be not properly used. Conversely, fishery-related variables seem to be the most reliable for this purpose.

During the study period, the fleet was practically constant, making monthly trips a good indicator of fishing effort and landing ability, and remained almost constant despite potential technological creep (Marriott et al. 2011) because no significant changes in fishing technology have occurred in the area in the last 20 years.

Studies of deep-water systems, where harsh conditions limit methods for evaluating fisheries, often suffer from a lack of data in order to assess stock status. Although the goal of fisheries managers is to promote sustainable production of fish stocks through formal stock assessment, it is often impractical to collect fishery-independent data in isolated or harsh environments. In these cases the information collected by a fishery is the main (or only) source of abundance data available (Maunder et al. 2006).