- Email: [email protected]
Use of generalized linear models to analyze catch rates having zero values: the Kuwait driftnet fishery
Fisheries Research 53 (2001) 151±168
Use of generalized linear models to analyze catch rates having zero values: the Kuwait driftnet ®shery Yimin Ye*, M. Al-Husaini, A. Al-Baz Mariculture and Fisheries Department, Kuwait Institute for Scienti®c Research, P.O. Box 1638, 22017 Salmiyah, Kuwait Received 18 May 2000; received in revised form 29 September 2000; accepted 23 October 2000
Abstract Fishery survey data often contain zero values. This paper uses a gamma-based generalized linear model for non-zero values and a Bernoulli-based model for the proportion of non-zero catches. Yearly variation in catch rates is then estimated by integrating the estimates from both models. The results for the Kuwait driftnet ®shery show that the catch rate of silver pomfret (Pampus argenteus) decreased from 115 kg per 1000 m clearances in 1984 to about 10 kg per 1000 m clearances in 1999, but its alternative target species, hilsa shad (Tenualosa ilisha), exhibited no clear pattern with great ¯uctuation over time. The proportion of non-zero catches for silver pomfret decreased with time during the study period. In contrast, the proportion for hilsa shad increased, indicating ®shermen shifted their targeting from silver pomfret to hilsa shad with the decreasing catch rate of silver pomfret. Other effect factors like month, boat type and landing port were also included in the models, and their effects on catch rates were estimated and tested. # 2001 Elsevier Science B.V. All rights reserved. Keywords: Catch rate; Generalized linear models; Driftnet ®sheries; Silver pomfret; Hilsa shad; Kuwait
1. Introduction The driftnet ®shery is one of the most important ®sheries in Kuwait. It involves two categories of ®shing boats, dhows and speedboats, and operates year-round. Although the driftnet ®shery catches more than 10 ®sh species, the majority of its catch comes from two species, hilsa shad (local name: suboor, Tenualosa ilisha) and silver pomfret (local name: zobaidy, Pampus argenteus), contributing 49.2 and 12.7%, respectively, on average, from 1984 to 1999. Although the ®shery is a traditional industry and an important source of food supply in Kuwait, the *
Corresponding author. Tel.: 965-5711295; fax: 965-5711293. E-mail address: [email protected] (Y. Ye).
driftnet ®shery remained unmanaged and no detailed ®n ®shery data were available until 1982 when a program for collecting ®n ®shery statistics and biological data started (Morgan, 1982). The Kuwait Institute for Scienti®c Research (KISR) has since carried out dockside interviews of skippers to collect data such as catch by species, ®shing effort expended, types of boat and ®shing gear, and landing port. It was expected that basic data for stock assessment such as catch and ®shing effort could be estimated from the interview survey. Catch per unit effort (CPUE) from commercial ®sheries has been used to derive indices of relative abundance or to estimate ®shing effort for many world ®sheries (Gulland, 1956; Robson, 1966; Large, 1992; Stefansson, 1996; Grif®n et al., 1997; Goni et al., 1999). However, the use of catch rates in constructing
0165-7836/01/$ ± see front matter # 2001 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 7 8 3 6 ( 0 0 ) 0 0 2 8 7 - 3
152
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
abundance indices or estimating ®shing effort requires standardization to take into account changes in the ability to catch ®sh, and ¯eet composition, and to adjust catch rate estimates for other factors that may affect the catch rates such as year, month, boat type, landing port, or abundance of other target species in the catch (Hilborn and Walters, 1992). The generalized linear modeling (GLM) technique has been used to estimate the effective extents of various factors because GLM allows identi®cation of the factors that in¯uence catch rates as well as computation of standardized catch rates, represented by the year effect factor after taking into account the effects of other factors, which are used in many stock assessment methods. The Kuwait driftnet ®shery has two main target species. When targeting one species, ®shermen usually catch none of the other species, although sometimes they may catch both. Consequently, for a single species, its interview catch records consist of many zero values. The available data show that 50% of the driftnet interview data recorded zero hilsa shad catches and 56% did not record any silver pomfret. These zero catch values create great dif®culty in estimating overall catch rates. A common way to deal with the zeros is to add a small arbitrary constant either to zero catches or to all catch records (Berry, 1987; Caputi, 1996; Fox and Starr, 1996; Ma and Hsul, 1997) when a linear model using log-transformed data is ®tted to estimate catch rates. However, if there are a reasonable number of zero catches with varying effort, this may introduce a signi®cant bias in that zero catches with low effort will appear to have ``higher'' catch rates than those with higher effort (Caputi, 1996). Different values of the arbitrary quantity result in different estimates, although the trends of the main year effect may be similar (Caputi, 1996). Of the two target species, the catch rate of either species is not represented only by the non-zero catches. Fish distribution may change from the average in a number of ways when its abundance declines. The density may stay constant in some areas, but its distributional area may diminish. This kind of change could result in the non-zero catch rates having the same mean over time, but the number of zero catches would increase (Stefansson, 1996). In the opposite case, the spatial distribution may stay the same, but the density at each point may drop, although never to zero.
This kind of change may lead to decreased non-zero catches, but the chance to catch no ®sh would stay the same. Fishermen would also switch from one target species to another, depending on their relative abundance and distribution. When a species is abundant, ®shermen may target the species and, consequently, both the probability of catching this species and its catch rate remain high. However, for the alternative target species, ®shermen may only catch them where a high-density area is located, and the probability of zero catches for this species then increases, but the non-zero catch rate may remain very high. Therefore, stock abundance or catch rate is represented by two sources of information in this case. One is the level of non-zero catch rates and the other is the probability of catching the species. A comprehensive index of abundance should integrate both kinds of information. This study is an attempt to incorporate the two sources of information in standardizing catch rates. For the Kuwait driftnet ®shery, the analysis is intended to: (1) identify factors that have signi®cant effects on the driftnet catch rates; (2) produce a time series of standardized catch rate estimates that can be used for stock assessment. 2. Materials and methods 2.1. Data Kuwait's driftnet ®shery data were derived from the ®n®shery interview program that also covers other ®sheries. The interview was carried out 6±7 days a month. A systematic sample design was used (1-in-5), but when the interview day fell on a weekend, it was deferred to the ®rst day of the following week. On each interview day, two major landing ports, Sief and Doha, were visited, and all ®shing boats harboring during the peak period between 1500 and 1700 h at Sief and between 1000 and 1300 h at Doha were contacted. The interview was made on board ®shing vessels with captains. Details of each boat's catch by species and various components of its ®shing effort were recorded, together with other registration information (Morgan, 1982). Although the interview program started in 1982, some years' data were lost during the Iraqi occupation of Kuwait in 1990. The data used in this study are
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
therefore restricted to the years 1984±1987, 1990 and 1992±1999. A total of 10 093 observations are available for the current analysis. Catch rates were calculated by catch in kilogram divided by effort in 1000 m clearances. The Kuwait driftnet ®shery captures multi-species. However, most species constitute a very small portion (<6%) of the total catch, except hilsa shad and silver pomfret (49.2 and 12.7%, respectively). Hilsa shad contributes a larger portion, but silver pomfret has a market price three times as high as hilsa shad. These two species, therefore, should be considered as the target species of the ®shery, on which the analysis of this study focuses. 2.2. The method The correct procedure in any analysis to account for both the probability of catching a species and the level of catch rates when a species is captured is to classify explicitly this dichotomy of the interview data into two categories, zero values and non-zero values. Consider a random variable (catch rate), x, with the following properties. There is a probability, p, where x is non-zero, Prfx > 0g p, and, hence, a probability of 1 p, where x is zero, Prfx 0g 1 p. Furthermore, the distribution of x conditional on x 6 0 is some well-known distribution of a positive variable. Then, the cumulative distribution function (cdf) of an individual catch x at port p in year y becomes (Aitchison, 1955; Smith, 1988) PrfXpy xg 1
ppy ppy Gpy x;
x>0
(1)
where Gpy is a continuous cdf describing the distribution of positive values and ppy the probability of interviewing a boat of non-zero catches in year y at port p. Following the approach of Stefansson (1996), the likelihood corresponding to the given cdf is given by Y Y L 1 ppy npy rpy prpypy f xpy (2) p;y
p;y:xpy >0
where npy is the number of repetitions at port p in year y, rpy the number of positive values at port p obtained in year y and f(xpy) the density function of the positive values. In Eq. (2), there are two distinct components: the probability of a non-zero value and the distribution of
153
the non-zero values. These can, therefore, be modeled and ®tted separately to obtain ®rst a ®tted probability of non-zero values and then the expected catch rate, given that some were caught. The ®tted mean value for port p in year y is then given by ppyupy (Pennigton, 1983; Stefansson, 1996). Particularly, when the distribution of f(xpy) is a member of the exponential family, both parts of Eq. (2) can easily be ®tted by GLM (McCullagh and Nelder, 1989). The frequency distribution of the silver pomfret catch rates was skewed, having a large number of zero values and a heavy tail. The histograms using log x 0:001 show that the zero values do not form a natural part of the same continuous distribution as the one producing positive values (Fig. 1). When zero values are eliminated, it is seen that the data may be close to lognormal, which implies that a lognormal or gamma distribution may be appropriate for positive values (Stefansson, 1996). A log±log plot of the variance versus average catch rate within year (Fig. 2) shows that the regression slope is close to 2 for both species, meaning that the variance is proportional to the square of the mean response and a gamma model should be used (McCullagh and Nelder, 1989). Therefore, the following gamma model was used for both species in this study f xpy
xrpy 1 e rxpy =m G rm=rr
(3)
where m is the mean and r the shape parameter (Dobson, 1990). After determining the distribution, the following generalized linear model was used to analyze the nonzero values with the extension of main effect factors to be investigated to year, month, boat type and port ln mymbp a by ym lb tp
(4)
where mymbp is the expected catch rate for year y, month m, boat type b and port p; a the catch rate obtained in January 1984 by dhow boat at Doha port; by the catch rate in year y relative to 1984; ym the abundance of month m relative to January; lb the ef®ciency of speedboat relative to dhow boat; tp the difference between Sief port and Doha port. The model was ®tted by setting the S-Plus GLM routine at a gamma distribution with a log-link function, and the ®rst-order interactions can be added easily.
154
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
Fig. 1. Frequency distribution of logged slightly shifted CPUE data for silver pomfret, log CPUE 0:001.
To model the probability of non-empty catches, the data were ®rst recorded so that, for each interview, the value 0 was assigned if no ®sh was caught and the value 1 was recorded for non-zero catches to obtain Bernoulli-type 0/1-measurements. The model for probabilities is via the logit function with binomial distribution, so that if the probability of a non-zero value is thought to depend on year, month, boat and port, then it would be appropriate to model the existence of a non-zero catch in an interview as a Bernoulli random variable with probability, p, given by pymbp ln (5) a0 b0y y0m l0b t0p 1 pymbp where the parameters have meanings similar to those in Eq. (3). The value of the probability of a non-zero catch can be obtained by solving Eq. (5). This model is a special case of generalized linear models, which are available in the statistical package of S-Plus. The software does this transformation automatically. Naturally, interaction terms can be added to Eq. (5), as in any GLMs. After obtaining ®tted values for the probability, p, of a non-zero catch rate, and for the expected catch rate, m, conditional on it being positive, the predicted
unconditional catch rate is given by pm. The present model combines two sub-GLMs, and the overall estimate of catch rate is dependent on both sub-models. If the model for the proportion (Eq. (5)) indicates that the probability of a non-zero catch rate differs between ports and the model for the mean of the non-zero catch rates (Eq. (4)) indicates that these depend on the type of boat, then the overall mean, pm, depends on both the port and the boat type (Stefansson, 1996). One of the main purposes of the GLM analysis of catch rates is to provide year effects. For a main effect model, year effects can be derived from the coef®cients by setting ``options contrasts c ``contr: treatment''; ``contr:poly'''' in S-Plus. However, if the year effects are to be extracted and used, special care needs to be taken if the model contains interaction terms with year. The common approach of extracting the year effect alone from a log-linear model can be replaced by an integral of the ®tted model over the entire scale under consideration. This approach, however, yields catch rate indices that are equivalent to the year effects when the model contains no interaction terms (Stefansson, 1996). Therefore, only the main effect model was used to estimate the variation of catch rates over year in this study.
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
155
Fig. 2. Scatterplot of log-variance versus log-mean for yearly catch rates of silver pomfret (a) and hilsa shad (b).
The present model consists of two sub-models, and even when the year effect is only present as a main effect in each of the sub-models, the overall model will contain a non-linear year effect. Thus, integration of the ®tted response surface is required to obtain the
overall catch rate. For simplicity, the present approach is to compute the average ®tted catch rates across boat types and months with each landing port, and the overall integral is taken as the direct average over all landing ports (Stefansson, 1996).
156
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
3. Results 3.1. Silver pomfret The analysis of deviance for silver pomfret shows that the difference in catch rates, expressed by CPUE, between years is signi®cant as well as the variation of catch rates with month, boat type and landing port P < 0:0001 (Table 1). Beside the main effects, the ®rst-order interaction between year and month is also signi®cant P < 0:0001, meaning the monthly variation in catch rates is not the same in all years. However, the difference in catch rates between boat types is not signi®cant and does not change with year P 0:3145, month P 0:2590 and port P 0:5546. Although the port effect is consistent with month P 0:5138 and boat type P 0:5546, it changes with year P 0:0126. The model explained 46.2% of the deviance, most of which is explained by the difference between years, 33.5% (Table 1). However, although there are signi®cant ®rst-order interactions present, the variation in ln(CPUE) explained by each effect is relatively small in comparison with the variation explained by the corresponding main effects in the model. The contribution of each of the main effects included in the model to the variation of silver pomfret catch rates accounting for other factors is shown in Fig. 3. In general, ln(CPUE) has declined over years (Fig. 3a).
There is a strong seasonal trend in catch rate, being high from May to November and low from December to April (Fig. 3b). The dhow boats have a higher catch rate than speedboats (Fig. 3c) and the catch rates of the boats that land in Sief are higher than those of the boats landing in Doha (Fig. 3d). Departures from the general patterns shown in Fig. 3 are indicated by the signi®cance of particular interactions. For example, the signi®cance of the year:month interaction, which accounts for 7.9% of the explained model deviance (Table 1), indicates that monthly variation in catch rates changes with year. If year-speci®c models were ®tted, the results would show that the estimated seasonal trends for some years follow the general pattern shown in Fig. 3b, but those for other years may not be consistent with the general pattern. Although the interaction year:port is also signi®cant (P 0:0126, Table 1), it accounts only 0.5% of the explained deviance and can be disregarded. The main effect model estimates the coef®cients expressing the difference between each level of the factors and the ®rst level. For year, the estimates represent the difference between the mean ln(CPUE) for each year and the mean ln(CPUE) for dhow boat at Doha port in January 1984 after adjusting the main effects of boat type, month and landing port. Unbiased estimates of catch rates can be obtained by the anti-logarithm of the estimates of ln(CPUE) by adding a correction factor (0:5 residual mean square
Table 1 Analysis of deviance table for the gamma-based GLM ®tted to silver pomfret catch rate data Source of variation
d.f.
Deviance
% Explained
Null Main effects Year Month Boat type Port
Res. d.f.
Res. deviance
5896
15682
F
Pr(F)
12 11 1 1
5253 324 70 163
33.5 2.1 0.5 1.0
5974 5963 5962 5961
10374 10050 9979 9815
113.5 7.6 18.4 42.5
<0.0001 <0.0001 <0.0001 <0.0001
Interactions Year:month Year:boat type Year:port Month:boat type Month:port Boat type:port
100 10 8 10 10 1
1237 45 75 42 35 1
7.9 0.3 0.5 0.3 0.2 0.0
5861 5851 5843 5833 5823 5822
8577 8532 8457 8409 8374 8372
3.2 1.2 2.4 1.2 0.9 0.3
<0.0001 0.3145 0.0126 0.2590 0.5138 0.5546
Total explained
164
7256
46.2
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
157
Fig. 3. Results of the main effect model ®tted to non-zero silver pomfret catch rates incorporating logarithmic link with gamma variance functions and main effects: (a) year, (b) month, (c) boat type and (d) port. Each plot represents the contribution of the corresponding variable to the ®tted linear predictor. The ®tted values are adjusted to average zero and the broken bars indicate two standard errors. The width of the solid bars at the base of the plots is proportional to the number of observations at each level of the factor.
158
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
Fig. 3. (Continued ).
deviation) to reduce any bias that may be introduced. However, for the comparison of yearly differences in catch rate, the intercept was replaced by the estimated annual mean ln(CPUE) in 1984 for the real combination
of boat type, month and landing port (Large, 1992), and the correction factor can be simply integrated into the intercept. The annual variation of catch rates was then plotted directly against time as year effects (Fig. 4).
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
159
Fig. 4. Trend of silver pomfret catch rates (CPUE) over the study period, as estimated by the main effect model with the non-zero values.
The percentage of the non-zero catch rates was based on the Bernoulli model (Eq. (5)). Because the data are 0/1-measurements, the deviance will have a w2-distribution, and, therefore, a w2-statistic was used to test for signi®cance in Table 2. Although the
reduction in deviance is small in relation to the total deviance, all factors are considerable in comparison with the degrees of freedom expended and, therefore, highly signi®cant, except the ®rst-order interaction between boat type and landing port (Table 2). The
Table 2 Analysis of deviance table for the Bernoulli-based model ®tted to silver pomfret catch rate data Source of variation
d.f.
Deviance
% Explained
Null Main effects Year Month Boat type Port
Res. d.f.
Res. deviance
10092
13639
Pr(Chi)
12 11 1 1
1157 525 11 636
8.5 3.8 0.1 4.7
10080 10069 10068 10067
12481 11956 11944 11308
<0.0001 <0.0001 0.0008 <0.0001
Interactions Year:month Year:boat type Year:port Month:boat type Month:port Boat type:port
112 11 11 10 11 1
1316 50 170 35 118 0
9.6 0.4 1.2 2.3 0.9 0.0
9955 9944 9933 9923 9912 9911
9992 9942 9772 9738 9619 9619
<0.0001 <0.0001 <0.0001 0.0001 <0.0001 0.9317
Total explained
181
4020
29.5
160
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
Fig. 5. Results of the main effect model ®tted to silver pomfret catch rates incorporating logistic link with Bernoulli variance functions and main effects: (a) year, (b) month, (c) boat type and (d) port. Each plot represents the contribution of the corresponding variable to the ®tted linear predictor. The ®tted values are adjusted to average zero and the broken bars indicate two standard errors. The width of the solid bars at the base of the plots is proportional to the number of observations at each level of the factor.
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
161
Fig. 5. (Continued ).
variation in percentage of positive values with month, boat type and landing port changes from year to year, and the differences between boat types and between landing ports are not the same in all the months. The model
explained 29.5% of the variation, and the main effect factors explained 17.1% (Table 2). This small percent is not unusual because the current model attempts to explain data on a very ®ne scale (Stefansson, 1996).
162
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
The contribution of each main effect factor to the variation in the percentage of success is presented in Fig. 5. A general declining trend is quite clear for year effect (Fig. 5a). The monthly trend is similar to the two-stage pattern appearing in non-zero catch rates (Fig. 3b), being high from June to December and low from January to May, and a peak is seen in June (Fig. 5b). Of the two boat types, speedboats have a higher probability to catch ®sh than dhow boats (Fig. 5c), and the boats at Doha port show a greater percent of having non-zero catches than those at Sief port (Fig. 5d). Based on the same method used in calculating catch rates, the yearly variation in the percentage of nonzero catches was estimated and shown in Fig. 6. A decreasing trend is seen, with exceptions of 1990 and 1992. In 1984, over 90% of boats recorded positive silver pomfret catches, but this percentage fell to 20% in 1999 (Fig. 6). The combined estimate of catch rates was about 115 kg per 1000 m clearances in 1984 but suddenly dropped to about 25 kg per 1000 m clearances in 1985 (Fig. 7). A decreasing trend has since remained, and the estimated catch rate in 1999 was about 10 kg per 1000 m clearances, less than one-tenth of the 1984 level (Fig. 7).
3.2. Hilsa shad Hilsa shad's analysis of deviance shows that all the main effect factors in the model are signi®cant with P < 0:0001 for year, month and landing port, and P 0:0122 for boat type (Table 3). The signi®cance of the interactions between year and month P < 0:0001 and between year and port P < 0:0001 means that the monthly pattern in catch rates as well as the difference between landing ports are not the same in all years. The monthly pattern in catch rates also varies with landing port P < 0:0001. However, the ®rst-order interactions between boat type and year and between boat type and month are not signi®cant, P 0:0698 and 0.5303, respectively, indicating that the difference in catch rates between the two boat types remains regardless of year and month. With the same method as in the analysis for silver pomfret, the yearly variation of positive catch rate values after taking into account other factors was estimated from the main effect model coef®cients and depicted in Fig. 8. The year 1984 shows an extremely high catch rate of over 200 kg per 1000 m clearances, but the 3 years of 1990, 1992 and 1993 have rather low catch rates from 10 to
Fig. 6. Trend in proportion of non-zero silver pomfret catches over the study period estimated by the Bernoulli-based main effect model.
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
163
Fig. 7. Time series of silver pomfret catch rates estimated by integrating the results from the gamma-based and Bernoulli-based models.
15 kg per 1000 m clearances. The catch rate ¯uctuates at around 50 kg per 1000 m clearances during the rest of the study period. The main effect Bernoulli model for hilsa shad explained 11.3% of the data variance (Table 4). All the four main factors are highly signi®cant, with year
and month being the main effect factors that explain great portions of the variance (Table 4). The Bernoulli model with ®rst-order interactions explained 23.2% of the total variation in the hilsa shad data. Of the six ®rst-order interactions, only the one between boat type and port is not signi®cant with P 0:4851 (Table 4).
Table 3 Analysis of deviance table for the gamma-based GLM ®tted to hilsa shad catch rate dataa Source of variation
d.f.
Deviance
% Explained
Null
Res. d.f.
Res. deviance
5329
14202
F
Pr(F)
Main effects Year Month Boat type Port
12 11 1 1
1966 600 19 684
13.8 4.2 0.1 4.8
5317 5306 5306 5304
12235 11725 11706 11022
54.4 15.4 6.3 227.4
<0.0001 <0.0001 0.0122 <0.0001
Interactions Year:month Year:boat type Year:port Month:boat type Month:port
95 8 10 9 11
1464 44 124 24 134
10.3 0.3 0.9 0.2 0.9
5209 5201 5191 5182 5171
9558 9514 9390 9366 9232
5.1 1.8 4.1 0.9 4.1
<0.0001 0.0698 <0.0001 0.5304 <0.0001
Total explained
164
7256
35.6
a
The interaction between boat type and port was not included because it would increase residual deviance.
164
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
Fig. 8. Trend of hilsa shad catch rates (CPUE) over the study period, as estimated by the main effect model with the non-zero values.
The year±month interaction explains 10.2% of the sum of the squares and is an important interaction because of its effect on the interpretation of the annual variation in catch rate. The monthly pattern in the percentage of non-zero catches changes with year P < 0:0001. Although the difference in catch rates
between boat types varies with year P 0:0019 and month P < 0:0001, it remains regardless of landing ports P 0:4851 (Table 4). The proportion of positive catches, calculated from the coef®cients of the main effect model for the hilsa shad ®shery, shows an increasing trend over time, with
Table 4 Analysis of deviance table for the Bernoulli-based model ®tted to hilsa shad catch rate data Source of variation
d.f.
Deviance
% Explained
Null Main effects Year Month Boat type Port
Res. d.f.
Res. deviance
10091
13958
Pr(Chi)
12 11 1 1
598 601 11 396
4.3 4.3 0.1 2.8
10079 10068 10067 10066
13360 12759 12748 12352
<0.0001 <0.0001 0.0010 <0.0001
Interactions Year:month Year:boat type Year:port Month:boat type Month:port Boat type:port
112 11 11 10 11 1
1422 29 99 44 91 1
10.2 0.2 0.7 0.3 0.7 0.0
9954 9943 9932 9922 9911 9910
10929 10900 10800 10756 10664 10664
<0.0001 0.0019 <0.0001 <0.0001 <0.0001 0.4851
Total explained
181
3294
23.6
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
165
Fig. 9. Trend in the proportion of non-zero hilsa shad catches estimated by the Bernoulli-based main effect model.
Fig. 10. Time series of hilsa shad catch rates estimated by integrating the results from the gamma-based and Bernoulli-based models.
166
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
exceptions of 1986 and 1987. It was only 34.8% in 1984; however, it started to increase in 1993 and reached 83.1% in 1999 (Fig. 9). The hilsa shad catch rates were ®nally estimated by integrating the gamma-based model with the Bernoulli-based one. Fig. 10 shows the year-to-year variation in the combined estimates of catch rates, exhibiting a great variation with exceptionally low values in 1990, 1992 and 1993. The highest catch rate occurred in 1986, and there is no clear pattern that can be identi®ed. 4. Discussion The results from the silver pomfret ®shery show that the annual catch rate estimated from non-zero records decreased from 1984 to 1990, but it remained stable at a low level for the remaining years (Fig. 4). The proportion of positive catch rates has been declining over time, in general (Fig. 6). The decreasing non-zero catch rate suggests the declining abundance of the target species, which may also be the reason for the increased proportion of zero catches because ®shermen would have a smaller probability of successful ®shing with the declining stock abundance. Both the proportion and level of non-zero catches are related to stock abundance. Therefore, the combined estimates represent more comprehensive indices for catch rates. In the case of the silver pomfret ®shery, a more drastic decrease in catch rate is seen, and the decreasing trend becomes clearer over the entire study period (Fig. 7), in contrast with the catch rate estimates derived by using only non-zero values that showed that the catch rate has been relatively stable since 1992 (Fig. 4). The increase in targeting an alternative target species, hilsa shad in this case, may have caused the decreasing proportion of positive catches of silver pomfret. Although market price of silver pomfret is three times higher than that of hilsa shad, if hilsa shad becomes more abundant, ®shermen may switch to target it, and, consequently, the proportion of zero catches of silver pomfret increases. Under this circumstance, however, the level of non-zero catch rates of silver pomfret should not exhibit a decrease, but rather an increase because ®shermen attracted by the high abundance of the alternative species will keep
targeting it only when an area of higher density than normal can be located. For the silver pomfret ®shery, both the proportion and level of non-zero catch rates show a decreasing trend over time (Figs. 4 and 6), and, therefore, the decrease in the proportion could not result from the unusually high abundance of the alternative target species. The catch rate of hilsa shad estimated from the non-zero records had an exceptionally high value in 1984, but rather low values from 1990 to 1993 (Fig. 8). In general, the catch rate has been ¯uctuating at a low level, except for the 3 years affected by the Gulf war. In contrast, the estimated proportion of non-zero catches shows an increasing trend after 1990 (Fig. 9), which may be caused by the declining probability of success in silver pomfret ®shing. The catch rates estimated by combining the two explicit models for zero and non-zero values differ from those derived from only the non-zero model, ¯uctuating greatly without a clear decreasing trend if the 3 years having exceptionally low levels caused by the war are excluded (Fig. 10). Both of the species, silver pomfret and hilsa shad, are the targets of the driftnet ®shery. Fishermen may target one rather than the other depending on considerations of stock abundance, ®shing location, market price, etc. Increasing the proportion of non-zero catches for one species may consequently lead to the rising of the other species' proportion of zero catches. In the annual variation, the proportion of silver pomfret's non-zero catches decreased over time (Fig. 6), but hilsa shad's proportion increased (Fig. 9), substituting the target species that has a decreasing abundance with the more abundant one. Such a substitution should also exist over a ®ner time scale like a month. The swing of targeting species is natural, and is an effective way for ®shermen to survive the drastic decline or collapse of some particular ®sh species. The more comprehensive method used in this study should provide more reliable estimates of catch rates than the traditional averaging method. The analysis of deviance shows that the difference in catch rates between dhow and speedboat is signi®cant for both species, and remains regardless of year, month and landing port (Tables 1 and 4). A dhow boat is quite different from a speedboat. Dhows have a wooden structure with a length of about 20 m. The boats are ®tted with inboard diesel engines ranging
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
from 132 to 365 hp, and have a crew of 4±6 ®shermen (Abdul-Ghaffar and Al-Ghunaim, 1994). Speedboats are smaller, normally 7±8 m long, with two engines of about 100 hp each. Only two ®shermen are involved with one speedboat. Although this study used 1000 m clearances as the unit of ®shing effort that can measure ®shing capacity of different kinds of boats relatively better, the differences in size and net structure still form the reasons for the discrepancy in catch rates. This result indicates that the overall catch rate of the driftnet ¯eet changes if the composition of ®shing boats varies, which must be taken into account when year-to-year variation in catch rates is to be estimated. However, the difference in the proportion of nonzero catches between the two boat types always varies with year, month and landing port (Tables 3 and 5). This may indicate that targeting which species in a ®shing operation is a human decision based on combined considerations of catch volume, location, market price, etc., rather than determined by speci®c time, landing place and boat type. There is no difference in making this decision between the two boat types. Results of this study indicate that using combined generalized linear models to analyze catch rates of the Kuwait driftnet ®shery is a sensible method for obtaining standardized indices. However, this ®shery is rarely documented and there is no other information available for verifying the reliability of the estimated year-to-year variation in catch rates. It would be desirable to compare the standardized catch rates with other estimates derived from resource surveys or catch
167
rate studies if they were available. Caution should be exercised when catch rates are compared with stock abundance indices because catch rates may not be closely related to abundance, particularly when some factors that in¯uence catch rates are not included in the GLM analysis. The data used in this study are highly unbalanced. The solid bars at the base of Figs. 3 and 5 indicate the relative proportion of observations at each level of the factors. In general, there are very few observations from December to February because many ®shermen stop ®shing during the winter season (Table 5). The years before 1990 had fewer interviews than the other years (Table 5). The GLM in S-Plus can take the unbalanced data into account, which can be a great advantage, but the precision would be improved by a more balanced design. The precision of the mean catch rate calculated for a certain month of a year depends on the number of observations. If no interviews were made for a certain month, then it was excluded in the GLM calculation in S-Plus, and this would affect the estimates of the coef®cients. This should be a concern in a future interview survey. The approach used in this study is based on two explicit generalized linear models, one for non-zero values and the other for the proportion of the non-zero catches in the entire data set. This provides an analytical technique where many problems usually associated with zero values are alleviated (Stefansson, 1996). This technique should result in a major improvement in the interpretation of the catch rate
Table 5 Number of observations for silver pomfret in each month Year Total
January February
March
April
May
June
July
August
September
October
November December
1984 88 1985 145 1986 139 1987 171 1990 108 1992 864 1993 329 1994 747 1995 797 1996 360 1997 4482 1998 1075 1999 787
0 3 2 11 0 0 0 7 29 0 0 2 56
20 19 16 21 47 95 21 36 51 47 213 105 58
7 18 13 15 26 206 34 92 174 44 289 210 63
6 8 18 14 17 251 26 191 126 2 113 101 63
14 11 15 11 9 97 19 108 86 46 2132 116 54
12 17 8 16 1 74 30 109 110 65 348 114 38
9 26 11 10 0 46 39 89 55 41 400 121 130
7 0 10 12 0 20 25 44 67 60 369 65 58
12 18 23 13 0 27 75 26 64 27 442 89 88
0 13 14 6 0 33 54 16 22 16 144 67 109
1 11 0 21 8 15 6 25 13 12 2 31 34
0 1 0 21 0 0 0 4 0 0 30 54 36
168
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
data compared with other traditional methods such as average catch rates. The models can accommodate temporal and spatial variability as well as the variability of other categories such as gear type, boat type and environmental factors. This potential for incorporating, estimating and testing effects of various variations is very useful for analysis of ®sh catch rates. References Abdul-Ghaffar, A.R., Al-Ghunaim, A.Y.Y., 1994. Review of Kuwait's shrimp ®sheries, their development and present status. In: Proceedings of the Technical Consultation on Shrimp Management in the Arabian Gulf, Al-Khobar, Saudi Arabia, November 6±8, 1994, pp. 1±26. Aitchison, J., 1955. On the distribution of a positive random variable having a discrete probability mass at the origin. J. Am. Statist. Assoc. 50, 901±908. Berry, D.A., 1987. Logarithmic transformations in ANOVA. Biometrics 43, 439±456. Caputi, N., 1996. Analysis of catch and effort data for red®sh, blue grenadier and orange roughy using generalized linear models. Report for Bureau of Resource Sciences. Department of Primary Industries and Energy, Australia. Dobson, A.J., 1990. An Introduction of Generalized Linear Models. Chapman & Hall, London. Fox, D.S., Starr, R.M., 1996. Comparison of commercial ®shery and research catch data. Can. J. Aquat. Sci. 53, 2681±2694.
Goni, R., Alvarez, F., Adlerstein, S., 1999. Application of generalized linear modeling to catch rate analysis of western Mediterranean ®sheries: the Castellon trawl ¯eet as a case study. Fish. Res. 42, 291±302. Grif®n, W.L., Shah, A.K., Nance, J.M., 1997. Estimation of standardized effort in the heterogeneous Gulf of Mexico shrimp ¯eet. Mar. Fish. Rev. 59 (3), 23±33. Gulland, J.A., 1956. On the ®shing effort in English demersal ®sheries. Fish. Invest., London Ser. 2 (20), 1±41. Hilborn, R., Walters, C.J., 1992. Quantitative Fisheries Stock Assessment, Choice, Dynamics and Uncertainty. Chapman & Hall, London. Large, P.A., 1992. Use of a multiplicative model to estimate relative abundance from commercial CPUE data. ICES J. Mar. Sci. 49, 253±261. Ma, C.T., Hsul, C.C., 1997. Standardization of CPUE for yellow®n tuna caught by Taiwanese longline ®shery in the Atlantic. Collect-Vol-Sci-Pap-Iccat 46 (4), 197±204. McCullagh, P., Nelder, J.A., 1989. Generalized Linear Models. Chapman & Hall, London. Morgan, G.R., 1982. Methodology for estimating catch and ®shing effort in Kuwait's ®sh ®sheries. Technical Report KISR 645. Kuwait Institute for Scienti®c Research, Kuwait. Pennigton, M., 1983. Ef®cient estimators of abundance for ®sh and plankton surveys. Biometrics 39, 281±286. Robson, D.S., 1966. Estimation of relative ®shing power of individual ships. Research Bulletin No. 3. International Commission for the Northwest Atlantic Fisheries, pp. 5±14. Smith, S.J., 1988. Evaluating the ef®ciency of the delta-distribution mean estimator. Biometrics 44, 485±493. Stefansson, G., 1996. Analysis of ground®sh survey abundance data: combining the GLM and delta approaches. ICES J. Mar. Sci. 53, 577±588.
Use of generalized linear models to analyze catch rates having zero values: the Kuwait driftnet ®shery Yimin Ye*, M. Al-Husaini, A. Al-Baz Mariculture and Fisheries Department, Kuwait Institute for Scienti®c Research, P.O. Box 1638, 22017 Salmiyah, Kuwait Received 18 May 2000; received in revised form 29 September 2000; accepted 23 October 2000
Abstract Fishery survey data often contain zero values. This paper uses a gamma-based generalized linear model for non-zero values and a Bernoulli-based model for the proportion of non-zero catches. Yearly variation in catch rates is then estimated by integrating the estimates from both models. The results for the Kuwait driftnet ®shery show that the catch rate of silver pomfret (Pampus argenteus) decreased from 115 kg per 1000 m clearances in 1984 to about 10 kg per 1000 m clearances in 1999, but its alternative target species, hilsa shad (Tenualosa ilisha), exhibited no clear pattern with great ¯uctuation over time. The proportion of non-zero catches for silver pomfret decreased with time during the study period. In contrast, the proportion for hilsa shad increased, indicating ®shermen shifted their targeting from silver pomfret to hilsa shad with the decreasing catch rate of silver pomfret. Other effect factors like month, boat type and landing port were also included in the models, and their effects on catch rates were estimated and tested. # 2001 Elsevier Science B.V. All rights reserved. Keywords: Catch rate; Generalized linear models; Driftnet ®sheries; Silver pomfret; Hilsa shad; Kuwait
1. Introduction The driftnet ®shery is one of the most important ®sheries in Kuwait. It involves two categories of ®shing boats, dhows and speedboats, and operates year-round. Although the driftnet ®shery catches more than 10 ®sh species, the majority of its catch comes from two species, hilsa shad (local name: suboor, Tenualosa ilisha) and silver pomfret (local name: zobaidy, Pampus argenteus), contributing 49.2 and 12.7%, respectively, on average, from 1984 to 1999. Although the ®shery is a traditional industry and an important source of food supply in Kuwait, the *
Corresponding author. Tel.: 965-5711295; fax: 965-5711293. E-mail address: [email protected] (Y. Ye).
driftnet ®shery remained unmanaged and no detailed ®n ®shery data were available until 1982 when a program for collecting ®n ®shery statistics and biological data started (Morgan, 1982). The Kuwait Institute for Scienti®c Research (KISR) has since carried out dockside interviews of skippers to collect data such as catch by species, ®shing effort expended, types of boat and ®shing gear, and landing port. It was expected that basic data for stock assessment such as catch and ®shing effort could be estimated from the interview survey. Catch per unit effort (CPUE) from commercial ®sheries has been used to derive indices of relative abundance or to estimate ®shing effort for many world ®sheries (Gulland, 1956; Robson, 1966; Large, 1992; Stefansson, 1996; Grif®n et al., 1997; Goni et al., 1999). However, the use of catch rates in constructing
0165-7836/01/$ ± see front matter # 2001 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 7 8 3 6 ( 0 0 ) 0 0 2 8 7 - 3
152
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
abundance indices or estimating ®shing effort requires standardization to take into account changes in the ability to catch ®sh, and ¯eet composition, and to adjust catch rate estimates for other factors that may affect the catch rates such as year, month, boat type, landing port, or abundance of other target species in the catch (Hilborn and Walters, 1992). The generalized linear modeling (GLM) technique has been used to estimate the effective extents of various factors because GLM allows identi®cation of the factors that in¯uence catch rates as well as computation of standardized catch rates, represented by the year effect factor after taking into account the effects of other factors, which are used in many stock assessment methods. The Kuwait driftnet ®shery has two main target species. When targeting one species, ®shermen usually catch none of the other species, although sometimes they may catch both. Consequently, for a single species, its interview catch records consist of many zero values. The available data show that 50% of the driftnet interview data recorded zero hilsa shad catches and 56% did not record any silver pomfret. These zero catch values create great dif®culty in estimating overall catch rates. A common way to deal with the zeros is to add a small arbitrary constant either to zero catches or to all catch records (Berry, 1987; Caputi, 1996; Fox and Starr, 1996; Ma and Hsul, 1997) when a linear model using log-transformed data is ®tted to estimate catch rates. However, if there are a reasonable number of zero catches with varying effort, this may introduce a signi®cant bias in that zero catches with low effort will appear to have ``higher'' catch rates than those with higher effort (Caputi, 1996). Different values of the arbitrary quantity result in different estimates, although the trends of the main year effect may be similar (Caputi, 1996). Of the two target species, the catch rate of either species is not represented only by the non-zero catches. Fish distribution may change from the average in a number of ways when its abundance declines. The density may stay constant in some areas, but its distributional area may diminish. This kind of change could result in the non-zero catch rates having the same mean over time, but the number of zero catches would increase (Stefansson, 1996). In the opposite case, the spatial distribution may stay the same, but the density at each point may drop, although never to zero.
This kind of change may lead to decreased non-zero catches, but the chance to catch no ®sh would stay the same. Fishermen would also switch from one target species to another, depending on their relative abundance and distribution. When a species is abundant, ®shermen may target the species and, consequently, both the probability of catching this species and its catch rate remain high. However, for the alternative target species, ®shermen may only catch them where a high-density area is located, and the probability of zero catches for this species then increases, but the non-zero catch rate may remain very high. Therefore, stock abundance or catch rate is represented by two sources of information in this case. One is the level of non-zero catch rates and the other is the probability of catching the species. A comprehensive index of abundance should integrate both kinds of information. This study is an attempt to incorporate the two sources of information in standardizing catch rates. For the Kuwait driftnet ®shery, the analysis is intended to: (1) identify factors that have signi®cant effects on the driftnet catch rates; (2) produce a time series of standardized catch rate estimates that can be used for stock assessment. 2. Materials and methods 2.1. Data Kuwait's driftnet ®shery data were derived from the ®n®shery interview program that also covers other ®sheries. The interview was carried out 6±7 days a month. A systematic sample design was used (1-in-5), but when the interview day fell on a weekend, it was deferred to the ®rst day of the following week. On each interview day, two major landing ports, Sief and Doha, were visited, and all ®shing boats harboring during the peak period between 1500 and 1700 h at Sief and between 1000 and 1300 h at Doha were contacted. The interview was made on board ®shing vessels with captains. Details of each boat's catch by species and various components of its ®shing effort were recorded, together with other registration information (Morgan, 1982). Although the interview program started in 1982, some years' data were lost during the Iraqi occupation of Kuwait in 1990. The data used in this study are
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
therefore restricted to the years 1984±1987, 1990 and 1992±1999. A total of 10 093 observations are available for the current analysis. Catch rates were calculated by catch in kilogram divided by effort in 1000 m clearances. The Kuwait driftnet ®shery captures multi-species. However, most species constitute a very small portion (<6%) of the total catch, except hilsa shad and silver pomfret (49.2 and 12.7%, respectively). Hilsa shad contributes a larger portion, but silver pomfret has a market price three times as high as hilsa shad. These two species, therefore, should be considered as the target species of the ®shery, on which the analysis of this study focuses. 2.2. The method The correct procedure in any analysis to account for both the probability of catching a species and the level of catch rates when a species is captured is to classify explicitly this dichotomy of the interview data into two categories, zero values and non-zero values. Consider a random variable (catch rate), x, with the following properties. There is a probability, p, where x is non-zero, Prfx > 0g p, and, hence, a probability of 1 p, where x is zero, Prfx 0g 1 p. Furthermore, the distribution of x conditional on x 6 0 is some well-known distribution of a positive variable. Then, the cumulative distribution function (cdf) of an individual catch x at port p in year y becomes (Aitchison, 1955; Smith, 1988) PrfXpy xg 1
ppy ppy Gpy x;
x>0
(1)
where Gpy is a continuous cdf describing the distribution of positive values and ppy the probability of interviewing a boat of non-zero catches in year y at port p. Following the approach of Stefansson (1996), the likelihood corresponding to the given cdf is given by Y Y L 1 ppy npy rpy prpypy f xpy (2) p;y
p;y:xpy >0
where npy is the number of repetitions at port p in year y, rpy the number of positive values at port p obtained in year y and f(xpy) the density function of the positive values. In Eq. (2), there are two distinct components: the probability of a non-zero value and the distribution of
153
the non-zero values. These can, therefore, be modeled and ®tted separately to obtain ®rst a ®tted probability of non-zero values and then the expected catch rate, given that some were caught. The ®tted mean value for port p in year y is then given by ppyupy (Pennigton, 1983; Stefansson, 1996). Particularly, when the distribution of f(xpy) is a member of the exponential family, both parts of Eq. (2) can easily be ®tted by GLM (McCullagh and Nelder, 1989). The frequency distribution of the silver pomfret catch rates was skewed, having a large number of zero values and a heavy tail. The histograms using log x 0:001 show that the zero values do not form a natural part of the same continuous distribution as the one producing positive values (Fig. 1). When zero values are eliminated, it is seen that the data may be close to lognormal, which implies that a lognormal or gamma distribution may be appropriate for positive values (Stefansson, 1996). A log±log plot of the variance versus average catch rate within year (Fig. 2) shows that the regression slope is close to 2 for both species, meaning that the variance is proportional to the square of the mean response and a gamma model should be used (McCullagh and Nelder, 1989). Therefore, the following gamma model was used for both species in this study f xpy
xrpy 1 e rxpy =m G rm=rr
(3)
where m is the mean and r the shape parameter (Dobson, 1990). After determining the distribution, the following generalized linear model was used to analyze the nonzero values with the extension of main effect factors to be investigated to year, month, boat type and port ln mymbp a by ym lb tp
(4)
where mymbp is the expected catch rate for year y, month m, boat type b and port p; a the catch rate obtained in January 1984 by dhow boat at Doha port; by the catch rate in year y relative to 1984; ym the abundance of month m relative to January; lb the ef®ciency of speedboat relative to dhow boat; tp the difference between Sief port and Doha port. The model was ®tted by setting the S-Plus GLM routine at a gamma distribution with a log-link function, and the ®rst-order interactions can be added easily.
154
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
Fig. 1. Frequency distribution of logged slightly shifted CPUE data for silver pomfret, log CPUE 0:001.
To model the probability of non-empty catches, the data were ®rst recorded so that, for each interview, the value 0 was assigned if no ®sh was caught and the value 1 was recorded for non-zero catches to obtain Bernoulli-type 0/1-measurements. The model for probabilities is via the logit function with binomial distribution, so that if the probability of a non-zero value is thought to depend on year, month, boat and port, then it would be appropriate to model the existence of a non-zero catch in an interview as a Bernoulli random variable with probability, p, given by pymbp ln (5) a0 b0y y0m l0b t0p 1 pymbp where the parameters have meanings similar to those in Eq. (3). The value of the probability of a non-zero catch can be obtained by solving Eq. (5). This model is a special case of generalized linear models, which are available in the statistical package of S-Plus. The software does this transformation automatically. Naturally, interaction terms can be added to Eq. (5), as in any GLMs. After obtaining ®tted values for the probability, p, of a non-zero catch rate, and for the expected catch rate, m, conditional on it being positive, the predicted
unconditional catch rate is given by pm. The present model combines two sub-GLMs, and the overall estimate of catch rate is dependent on both sub-models. If the model for the proportion (Eq. (5)) indicates that the probability of a non-zero catch rate differs between ports and the model for the mean of the non-zero catch rates (Eq. (4)) indicates that these depend on the type of boat, then the overall mean, pm, depends on both the port and the boat type (Stefansson, 1996). One of the main purposes of the GLM analysis of catch rates is to provide year effects. For a main effect model, year effects can be derived from the coef®cients by setting ``options contrasts c ``contr: treatment''; ``contr:poly'''' in S-Plus. However, if the year effects are to be extracted and used, special care needs to be taken if the model contains interaction terms with year. The common approach of extracting the year effect alone from a log-linear model can be replaced by an integral of the ®tted model over the entire scale under consideration. This approach, however, yields catch rate indices that are equivalent to the year effects when the model contains no interaction terms (Stefansson, 1996). Therefore, only the main effect model was used to estimate the variation of catch rates over year in this study.
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
155
Fig. 2. Scatterplot of log-variance versus log-mean for yearly catch rates of silver pomfret (a) and hilsa shad (b).
The present model consists of two sub-models, and even when the year effect is only present as a main effect in each of the sub-models, the overall model will contain a non-linear year effect. Thus, integration of the ®tted response surface is required to obtain the
overall catch rate. For simplicity, the present approach is to compute the average ®tted catch rates across boat types and months with each landing port, and the overall integral is taken as the direct average over all landing ports (Stefansson, 1996).
156
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
3. Results 3.1. Silver pomfret The analysis of deviance for silver pomfret shows that the difference in catch rates, expressed by CPUE, between years is signi®cant as well as the variation of catch rates with month, boat type and landing port P < 0:0001 (Table 1). Beside the main effects, the ®rst-order interaction between year and month is also signi®cant P < 0:0001, meaning the monthly variation in catch rates is not the same in all years. However, the difference in catch rates between boat types is not signi®cant and does not change with year P 0:3145, month P 0:2590 and port P 0:5546. Although the port effect is consistent with month P 0:5138 and boat type P 0:5546, it changes with year P 0:0126. The model explained 46.2% of the deviance, most of which is explained by the difference between years, 33.5% (Table 1). However, although there are signi®cant ®rst-order interactions present, the variation in ln(CPUE) explained by each effect is relatively small in comparison with the variation explained by the corresponding main effects in the model. The contribution of each of the main effects included in the model to the variation of silver pomfret catch rates accounting for other factors is shown in Fig. 3. In general, ln(CPUE) has declined over years (Fig. 3a).
There is a strong seasonal trend in catch rate, being high from May to November and low from December to April (Fig. 3b). The dhow boats have a higher catch rate than speedboats (Fig. 3c) and the catch rates of the boats that land in Sief are higher than those of the boats landing in Doha (Fig. 3d). Departures from the general patterns shown in Fig. 3 are indicated by the signi®cance of particular interactions. For example, the signi®cance of the year:month interaction, which accounts for 7.9% of the explained model deviance (Table 1), indicates that monthly variation in catch rates changes with year. If year-speci®c models were ®tted, the results would show that the estimated seasonal trends for some years follow the general pattern shown in Fig. 3b, but those for other years may not be consistent with the general pattern. Although the interaction year:port is also signi®cant (P 0:0126, Table 1), it accounts only 0.5% of the explained deviance and can be disregarded. The main effect model estimates the coef®cients expressing the difference between each level of the factors and the ®rst level. For year, the estimates represent the difference between the mean ln(CPUE) for each year and the mean ln(CPUE) for dhow boat at Doha port in January 1984 after adjusting the main effects of boat type, month and landing port. Unbiased estimates of catch rates can be obtained by the anti-logarithm of the estimates of ln(CPUE) by adding a correction factor (0:5 residual mean square
Table 1 Analysis of deviance table for the gamma-based GLM ®tted to silver pomfret catch rate data Source of variation
d.f.
Deviance
% Explained
Null Main effects Year Month Boat type Port
Res. d.f.
Res. deviance
5896
15682
F
Pr(F)
12 11 1 1
5253 324 70 163
33.5 2.1 0.5 1.0
5974 5963 5962 5961
10374 10050 9979 9815
113.5 7.6 18.4 42.5
<0.0001 <0.0001 <0.0001 <0.0001
Interactions Year:month Year:boat type Year:port Month:boat type Month:port Boat type:port
100 10 8 10 10 1
1237 45 75 42 35 1
7.9 0.3 0.5 0.3 0.2 0.0
5861 5851 5843 5833 5823 5822
8577 8532 8457 8409 8374 8372
3.2 1.2 2.4 1.2 0.9 0.3
<0.0001 0.3145 0.0126 0.2590 0.5138 0.5546
Total explained
164
7256
46.2
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
157
Fig. 3. Results of the main effect model ®tted to non-zero silver pomfret catch rates incorporating logarithmic link with gamma variance functions and main effects: (a) year, (b) month, (c) boat type and (d) port. Each plot represents the contribution of the corresponding variable to the ®tted linear predictor. The ®tted values are adjusted to average zero and the broken bars indicate two standard errors. The width of the solid bars at the base of the plots is proportional to the number of observations at each level of the factor.
158
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
Fig. 3. (Continued ).
deviation) to reduce any bias that may be introduced. However, for the comparison of yearly differences in catch rate, the intercept was replaced by the estimated annual mean ln(CPUE) in 1984 for the real combination
of boat type, month and landing port (Large, 1992), and the correction factor can be simply integrated into the intercept. The annual variation of catch rates was then plotted directly against time as year effects (Fig. 4).
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
159
Fig. 4. Trend of silver pomfret catch rates (CPUE) over the study period, as estimated by the main effect model with the non-zero values.
The percentage of the non-zero catch rates was based on the Bernoulli model (Eq. (5)). Because the data are 0/1-measurements, the deviance will have a w2-distribution, and, therefore, a w2-statistic was used to test for signi®cance in Table 2. Although the
reduction in deviance is small in relation to the total deviance, all factors are considerable in comparison with the degrees of freedom expended and, therefore, highly signi®cant, except the ®rst-order interaction between boat type and landing port (Table 2). The
Table 2 Analysis of deviance table for the Bernoulli-based model ®tted to silver pomfret catch rate data Source of variation
d.f.
Deviance
% Explained
Null Main effects Year Month Boat type Port
Res. d.f.
Res. deviance
10092
13639
Pr(Chi)
12 11 1 1
1157 525 11 636
8.5 3.8 0.1 4.7
10080 10069 10068 10067
12481 11956 11944 11308
<0.0001 <0.0001 0.0008 <0.0001
Interactions Year:month Year:boat type Year:port Month:boat type Month:port Boat type:port
112 11 11 10 11 1
1316 50 170 35 118 0
9.6 0.4 1.2 2.3 0.9 0.0
9955 9944 9933 9923 9912 9911
9992 9942 9772 9738 9619 9619
<0.0001 <0.0001 <0.0001 0.0001 <0.0001 0.9317
Total explained
181
4020
29.5
160
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
Fig. 5. Results of the main effect model ®tted to silver pomfret catch rates incorporating logistic link with Bernoulli variance functions and main effects: (a) year, (b) month, (c) boat type and (d) port. Each plot represents the contribution of the corresponding variable to the ®tted linear predictor. The ®tted values are adjusted to average zero and the broken bars indicate two standard errors. The width of the solid bars at the base of the plots is proportional to the number of observations at each level of the factor.
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
161
Fig. 5. (Continued ).
variation in percentage of positive values with month, boat type and landing port changes from year to year, and the differences between boat types and between landing ports are not the same in all the months. The model
explained 29.5% of the variation, and the main effect factors explained 17.1% (Table 2). This small percent is not unusual because the current model attempts to explain data on a very ®ne scale (Stefansson, 1996).
162
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
The contribution of each main effect factor to the variation in the percentage of success is presented in Fig. 5. A general declining trend is quite clear for year effect (Fig. 5a). The monthly trend is similar to the two-stage pattern appearing in non-zero catch rates (Fig. 3b), being high from June to December and low from January to May, and a peak is seen in June (Fig. 5b). Of the two boat types, speedboats have a higher probability to catch ®sh than dhow boats (Fig. 5c), and the boats at Doha port show a greater percent of having non-zero catches than those at Sief port (Fig. 5d). Based on the same method used in calculating catch rates, the yearly variation in the percentage of nonzero catches was estimated and shown in Fig. 6. A decreasing trend is seen, with exceptions of 1990 and 1992. In 1984, over 90% of boats recorded positive silver pomfret catches, but this percentage fell to 20% in 1999 (Fig. 6). The combined estimate of catch rates was about 115 kg per 1000 m clearances in 1984 but suddenly dropped to about 25 kg per 1000 m clearances in 1985 (Fig. 7). A decreasing trend has since remained, and the estimated catch rate in 1999 was about 10 kg per 1000 m clearances, less than one-tenth of the 1984 level (Fig. 7).
3.2. Hilsa shad Hilsa shad's analysis of deviance shows that all the main effect factors in the model are signi®cant with P < 0:0001 for year, month and landing port, and P 0:0122 for boat type (Table 3). The signi®cance of the interactions between year and month P < 0:0001 and between year and port P < 0:0001 means that the monthly pattern in catch rates as well as the difference between landing ports are not the same in all years. The monthly pattern in catch rates also varies with landing port P < 0:0001. However, the ®rst-order interactions between boat type and year and between boat type and month are not signi®cant, P 0:0698 and 0.5303, respectively, indicating that the difference in catch rates between the two boat types remains regardless of year and month. With the same method as in the analysis for silver pomfret, the yearly variation of positive catch rate values after taking into account other factors was estimated from the main effect model coef®cients and depicted in Fig. 8. The year 1984 shows an extremely high catch rate of over 200 kg per 1000 m clearances, but the 3 years of 1990, 1992 and 1993 have rather low catch rates from 10 to
Fig. 6. Trend in proportion of non-zero silver pomfret catches over the study period estimated by the Bernoulli-based main effect model.
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
163
Fig. 7. Time series of silver pomfret catch rates estimated by integrating the results from the gamma-based and Bernoulli-based models.
15 kg per 1000 m clearances. The catch rate ¯uctuates at around 50 kg per 1000 m clearances during the rest of the study period. The main effect Bernoulli model for hilsa shad explained 11.3% of the data variance (Table 4). All the four main factors are highly signi®cant, with year
and month being the main effect factors that explain great portions of the variance (Table 4). The Bernoulli model with ®rst-order interactions explained 23.2% of the total variation in the hilsa shad data. Of the six ®rst-order interactions, only the one between boat type and port is not signi®cant with P 0:4851 (Table 4).
Table 3 Analysis of deviance table for the gamma-based GLM ®tted to hilsa shad catch rate dataa Source of variation
d.f.
Deviance
% Explained
Null
Res. d.f.
Res. deviance
5329
14202
F
Pr(F)
Main effects Year Month Boat type Port
12 11 1 1
1966 600 19 684
13.8 4.2 0.1 4.8
5317 5306 5306 5304
12235 11725 11706 11022
54.4 15.4 6.3 227.4
<0.0001 <0.0001 0.0122 <0.0001
Interactions Year:month Year:boat type Year:port Month:boat type Month:port
95 8 10 9 11
1464 44 124 24 134
10.3 0.3 0.9 0.2 0.9
5209 5201 5191 5182 5171
9558 9514 9390 9366 9232
5.1 1.8 4.1 0.9 4.1
<0.0001 0.0698 <0.0001 0.5304 <0.0001
Total explained
164
7256
35.6
a
The interaction between boat type and port was not included because it would increase residual deviance.
164
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
Fig. 8. Trend of hilsa shad catch rates (CPUE) over the study period, as estimated by the main effect model with the non-zero values.
The year±month interaction explains 10.2% of the sum of the squares and is an important interaction because of its effect on the interpretation of the annual variation in catch rate. The monthly pattern in the percentage of non-zero catches changes with year P < 0:0001. Although the difference in catch rates
between boat types varies with year P 0:0019 and month P < 0:0001, it remains regardless of landing ports P 0:4851 (Table 4). The proportion of positive catches, calculated from the coef®cients of the main effect model for the hilsa shad ®shery, shows an increasing trend over time, with
Table 4 Analysis of deviance table for the Bernoulli-based model ®tted to hilsa shad catch rate data Source of variation
d.f.
Deviance
% Explained
Null Main effects Year Month Boat type Port
Res. d.f.
Res. deviance
10091
13958
Pr(Chi)
12 11 1 1
598 601 11 396
4.3 4.3 0.1 2.8
10079 10068 10067 10066
13360 12759 12748 12352
<0.0001 <0.0001 0.0010 <0.0001
Interactions Year:month Year:boat type Year:port Month:boat type Month:port Boat type:port
112 11 11 10 11 1
1422 29 99 44 91 1
10.2 0.2 0.7 0.3 0.7 0.0
9954 9943 9932 9922 9911 9910
10929 10900 10800 10756 10664 10664
<0.0001 0.0019 <0.0001 <0.0001 <0.0001 0.4851
Total explained
181
3294
23.6
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
165
Fig. 9. Trend in the proportion of non-zero hilsa shad catches estimated by the Bernoulli-based main effect model.
Fig. 10. Time series of hilsa shad catch rates estimated by integrating the results from the gamma-based and Bernoulli-based models.
166
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
exceptions of 1986 and 1987. It was only 34.8% in 1984; however, it started to increase in 1993 and reached 83.1% in 1999 (Fig. 9). The hilsa shad catch rates were ®nally estimated by integrating the gamma-based model with the Bernoulli-based one. Fig. 10 shows the year-to-year variation in the combined estimates of catch rates, exhibiting a great variation with exceptionally low values in 1990, 1992 and 1993. The highest catch rate occurred in 1986, and there is no clear pattern that can be identi®ed. 4. Discussion The results from the silver pomfret ®shery show that the annual catch rate estimated from non-zero records decreased from 1984 to 1990, but it remained stable at a low level for the remaining years (Fig. 4). The proportion of positive catch rates has been declining over time, in general (Fig. 6). The decreasing non-zero catch rate suggests the declining abundance of the target species, which may also be the reason for the increased proportion of zero catches because ®shermen would have a smaller probability of successful ®shing with the declining stock abundance. Both the proportion and level of non-zero catches are related to stock abundance. Therefore, the combined estimates represent more comprehensive indices for catch rates. In the case of the silver pomfret ®shery, a more drastic decrease in catch rate is seen, and the decreasing trend becomes clearer over the entire study period (Fig. 7), in contrast with the catch rate estimates derived by using only non-zero values that showed that the catch rate has been relatively stable since 1992 (Fig. 4). The increase in targeting an alternative target species, hilsa shad in this case, may have caused the decreasing proportion of positive catches of silver pomfret. Although market price of silver pomfret is three times higher than that of hilsa shad, if hilsa shad becomes more abundant, ®shermen may switch to target it, and, consequently, the proportion of zero catches of silver pomfret increases. Under this circumstance, however, the level of non-zero catch rates of silver pomfret should not exhibit a decrease, but rather an increase because ®shermen attracted by the high abundance of the alternative species will keep
targeting it only when an area of higher density than normal can be located. For the silver pomfret ®shery, both the proportion and level of non-zero catch rates show a decreasing trend over time (Figs. 4 and 6), and, therefore, the decrease in the proportion could not result from the unusually high abundance of the alternative target species. The catch rate of hilsa shad estimated from the non-zero records had an exceptionally high value in 1984, but rather low values from 1990 to 1993 (Fig. 8). In general, the catch rate has been ¯uctuating at a low level, except for the 3 years affected by the Gulf war. In contrast, the estimated proportion of non-zero catches shows an increasing trend after 1990 (Fig. 9), which may be caused by the declining probability of success in silver pomfret ®shing. The catch rates estimated by combining the two explicit models for zero and non-zero values differ from those derived from only the non-zero model, ¯uctuating greatly without a clear decreasing trend if the 3 years having exceptionally low levels caused by the war are excluded (Fig. 10). Both of the species, silver pomfret and hilsa shad, are the targets of the driftnet ®shery. Fishermen may target one rather than the other depending on considerations of stock abundance, ®shing location, market price, etc. Increasing the proportion of non-zero catches for one species may consequently lead to the rising of the other species' proportion of zero catches. In the annual variation, the proportion of silver pomfret's non-zero catches decreased over time (Fig. 6), but hilsa shad's proportion increased (Fig. 9), substituting the target species that has a decreasing abundance with the more abundant one. Such a substitution should also exist over a ®ner time scale like a month. The swing of targeting species is natural, and is an effective way for ®shermen to survive the drastic decline or collapse of some particular ®sh species. The more comprehensive method used in this study should provide more reliable estimates of catch rates than the traditional averaging method. The analysis of deviance shows that the difference in catch rates between dhow and speedboat is signi®cant for both species, and remains regardless of year, month and landing port (Tables 1 and 4). A dhow boat is quite different from a speedboat. Dhows have a wooden structure with a length of about 20 m. The boats are ®tted with inboard diesel engines ranging
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
from 132 to 365 hp, and have a crew of 4±6 ®shermen (Abdul-Ghaffar and Al-Ghunaim, 1994). Speedboats are smaller, normally 7±8 m long, with two engines of about 100 hp each. Only two ®shermen are involved with one speedboat. Although this study used 1000 m clearances as the unit of ®shing effort that can measure ®shing capacity of different kinds of boats relatively better, the differences in size and net structure still form the reasons for the discrepancy in catch rates. This result indicates that the overall catch rate of the driftnet ¯eet changes if the composition of ®shing boats varies, which must be taken into account when year-to-year variation in catch rates is to be estimated. However, the difference in the proportion of nonzero catches between the two boat types always varies with year, month and landing port (Tables 3 and 5). This may indicate that targeting which species in a ®shing operation is a human decision based on combined considerations of catch volume, location, market price, etc., rather than determined by speci®c time, landing place and boat type. There is no difference in making this decision between the two boat types. Results of this study indicate that using combined generalized linear models to analyze catch rates of the Kuwait driftnet ®shery is a sensible method for obtaining standardized indices. However, this ®shery is rarely documented and there is no other information available for verifying the reliability of the estimated year-to-year variation in catch rates. It would be desirable to compare the standardized catch rates with other estimates derived from resource surveys or catch
167
rate studies if they were available. Caution should be exercised when catch rates are compared with stock abundance indices because catch rates may not be closely related to abundance, particularly when some factors that in¯uence catch rates are not included in the GLM analysis. The data used in this study are highly unbalanced. The solid bars at the base of Figs. 3 and 5 indicate the relative proportion of observations at each level of the factors. In general, there are very few observations from December to February because many ®shermen stop ®shing during the winter season (Table 5). The years before 1990 had fewer interviews than the other years (Table 5). The GLM in S-Plus can take the unbalanced data into account, which can be a great advantage, but the precision would be improved by a more balanced design. The precision of the mean catch rate calculated for a certain month of a year depends on the number of observations. If no interviews were made for a certain month, then it was excluded in the GLM calculation in S-Plus, and this would affect the estimates of the coef®cients. This should be a concern in a future interview survey. The approach used in this study is based on two explicit generalized linear models, one for non-zero values and the other for the proportion of the non-zero catches in the entire data set. This provides an analytical technique where many problems usually associated with zero values are alleviated (Stefansson, 1996). This technique should result in a major improvement in the interpretation of the catch rate
Table 5 Number of observations for silver pomfret in each month Year Total
January February
March
April
May
June
July
August
September
October
November December
1984 88 1985 145 1986 139 1987 171 1990 108 1992 864 1993 329 1994 747 1995 797 1996 360 1997 4482 1998 1075 1999 787
0 3 2 11 0 0 0 7 29 0 0 2 56
20 19 16 21 47 95 21 36 51 47 213 105 58
7 18 13 15 26 206 34 92 174 44 289 210 63
6 8 18 14 17 251 26 191 126 2 113 101 63
14 11 15 11 9 97 19 108 86 46 2132 116 54
12 17 8 16 1 74 30 109 110 65 348 114 38
9 26 11 10 0 46 39 89 55 41 400 121 130
7 0 10 12 0 20 25 44 67 60 369 65 58
12 18 23 13 0 27 75 26 64 27 442 89 88
0 13 14 6 0 33 54 16 22 16 144 67 109
1 11 0 21 8 15 6 25 13 12 2 31 34
0 1 0 21 0 0 0 4 0 0 30 54 36
168
Y. Ye et al. / Fisheries Research 53 (2001) 151±168
data compared with other traditional methods such as average catch rates. The models can accommodate temporal and spatial variability as well as the variability of other categories such as gear type, boat type and environmental factors. This potential for incorporating, estimating and testing effects of various variations is very useful for analysis of ®sh catch rates. References Abdul-Ghaffar, A.R., Al-Ghunaim, A.Y.Y., 1994. Review of Kuwait's shrimp ®sheries, their development and present status. In: Proceedings of the Technical Consultation on Shrimp Management in the Arabian Gulf, Al-Khobar, Saudi Arabia, November 6±8, 1994, pp. 1±26. Aitchison, J., 1955. On the distribution of a positive random variable having a discrete probability mass at the origin. J. Am. Statist. Assoc. 50, 901±908. Berry, D.A., 1987. Logarithmic transformations in ANOVA. Biometrics 43, 439±456. Caputi, N., 1996. Analysis of catch and effort data for red®sh, blue grenadier and orange roughy using generalized linear models. Report for Bureau of Resource Sciences. Department of Primary Industries and Energy, Australia. Dobson, A.J., 1990. An Introduction of Generalized Linear Models. Chapman & Hall, London. Fox, D.S., Starr, R.M., 1996. Comparison of commercial ®shery and research catch data. Can. J. Aquat. Sci. 53, 2681±2694.
Goni, R., Alvarez, F., Adlerstein, S., 1999. Application of generalized linear modeling to catch rate analysis of western Mediterranean ®sheries: the Castellon trawl ¯eet as a case study. Fish. Res. 42, 291±302. Grif®n, W.L., Shah, A.K., Nance, J.M., 1997. Estimation of standardized effort in the heterogeneous Gulf of Mexico shrimp ¯eet. Mar. Fish. Rev. 59 (3), 23±33. Gulland, J.A., 1956. On the ®shing effort in English demersal ®sheries. Fish. Invest., London Ser. 2 (20), 1±41. Hilborn, R., Walters, C.J., 1992. Quantitative Fisheries Stock Assessment, Choice, Dynamics and Uncertainty. Chapman & Hall, London. Large, P.A., 1992. Use of a multiplicative model to estimate relative abundance from commercial CPUE data. ICES J. Mar. Sci. 49, 253±261. Ma, C.T., Hsul, C.C., 1997. Standardization of CPUE for yellow®n tuna caught by Taiwanese longline ®shery in the Atlantic. Collect-Vol-Sci-Pap-Iccat 46 (4), 197±204. McCullagh, P., Nelder, J.A., 1989. Generalized Linear Models. Chapman & Hall, London. Morgan, G.R., 1982. Methodology for estimating catch and ®shing effort in Kuwait's ®sh ®sheries. Technical Report KISR 645. Kuwait Institute for Scienti®c Research, Kuwait. Pennigton, M., 1983. Ef®cient estimators of abundance for ®sh and plankton surveys. Biometrics 39, 281±286. Robson, D.S., 1966. Estimation of relative ®shing power of individual ships. Research Bulletin No. 3. International Commission for the Northwest Atlantic Fisheries, pp. 5±14. Smith, S.J., 1988. Evaluating the ef®ciency of the delta-distribution mean estimator. Biometrics 44, 485±493. Stefansson, G., 1996. Analysis of ground®sh survey abundance data: combining the GLM and delta approaches. ICES J. Mar. Sci. 53, 577±588.