Abstract: The ornate rock lobster (Panulirus ornatus) in Torres Strait supports locally important artisanal and commercial fisheries in both Australia and Papua New Guinea. In this paper we developed a new age-structured population model for the lobster fishery, cast the model into the framework of dynamic state-space modelling of modern time series methodology and used a Bayesian approach to produce posterior distributions of the model parameters. A stock recruitment relationship was incorporated into the age-structured model and its parameters were fitted together with other model parameters. The mean estimates of natural mortality rate (M) range from 0.47 to 0.77 year-1. Using the mode values of natural mortality rate and the two parameters defining the stock recruitment relationship, the maximum sustainable yield was estimated to be around 265 t (tail weight) with a corresponding fishing mortality rate (F) 0.4 year-1. The estimated fishing mortality rate exceeded this level in 1994, 1996, 1997 and 1999. To consider the natural variation in these parameters and the uncertainty involved in their estimates, Monte Carlo simulation was used to quantify the uncertainty in the relationship between sustainable yield and fishing mortality. The results show that the 50% credible interval of the maximum sustainable yield ranges from 110 to 390 t and its corresponding F ranges from 0.2 to 0.5 year-1. State-space models can take into account both measurement and process errors and the Markov chain Monte Carlo techniques produce posterior probability density functions for all the parameters of interest. All these should provide more information about the uncertainty involved in stock assessment and improve risk management in fisheries.
Introduction
The ornate rock lobster (Panulirus ornatus) is distributed across the Torres Strait between Australia and Paua New Guinea. Lobster fishing has been an ancient tradition for local islanders and has become the major source of income for them over the last two decades (Pitcher and Bishop, 1995). Along with the local islanders’ fishing, the rock lobster also supports an industrialized fleet consisting of 24 small freezer boats. In Australia, the freezer boat sector produced 59% of the total lobster catch and the local islander sector 41% in 2003/2004 (Australian Fisheries Management Authority, unpublished data). The annual total catches in Torres Strait over the last decade averaged around 300 t of tail weight, about 75% landed in Australia and 25% caught in Paua New Guinea (Ye et al., 2004a).
The rock lobsters are caught by divers fishing from dinghies. Divers use hookah compressors to fish in deep grounds (5-25 m) between reefs, whereas free divers fish on reefs in shallow waters (Pitcher et al., 1997). In 1987, a minimum legal size (100 mm tail length) was introduced to prevent growth overfishing and the number of licenses of the freezer boat sector was frozen to stop the further expansion of fishing effort in this sector (TPZJA, 1988). A two month (October and November) ban on the use of hookah diving equipment was introduced in 1993. In 2002, the minimum legal size was raised to 115 mm tail length, the hookah ban was extended to four months (October to January) and a two month (October-November) closure to free diving was introduced.
The lobster stock consists of mainly three age groups, 0, 1 and 2 year olds, with very few large 3-year-old lobsters. Larvae hatch during the Austral summer and settle in Torres Strait at an age of about 0.5 years. They recruit into the fishery around September the following year and are fished for about a year before most lobsters migrate out of Torres Strait to breed just before they reach 3 years old (Skewes et al., 1994; Pitcher et al., 1997). Lobsters from Torres Strait spawn once in their life only and experience a severe mortality after breeding (Dennis et al., 1992).
As age 2 lobsters constitute the major part of the fishable stock, the dynamics of the stock biomass exhibits strong seasonal fluctuation caused by variation in fishing pressure within a year. Production models that use only catch and effort data are usually not suitable for assessing stocks with high seasonal variation in fishing pressure. Also, production models are not suitable for assessing the consequences of introducing controls over minimum fishable size and target age groups, which are the principal regulatory measures adopted for the lobster fishery. Therefore, the population dynamics of the lobster fishery has been studied through the use of age-structured models (Pitcher et al., 1997; Pitcher et al., 2002; Ye et al., 2004a). The proportion and mean size of recruiting 1+ and fished 2 + year classes were estimated from tail size distribution using the Mix program (MacDonald and Pitcher, 1979).
Age composition is an important source of information for understanding fish population dynamics. With length, weight and age information, age-based models can follow year-classes through time to monitor population changes over time (Quinn and Deriso, 1999). However, it is difficult to take into account the two stochastic components, process error and measurement error, when fitting population dynamic models using maximum likelihood methods (Meyer and Millar, 1999).We, therefore, develop an age-structured model for the Torres Strait lobster fishery, cast the model into the framework of dynamic state-space modelling of modern time series methodology and use a Bayesian approach to produce posterior distributions of the model parameters.
Materials and Methods
The Age-structured State-space Model
Process Equation
The dynamics of a cohort of lobsters can be described by the following equation:
| (1) |
where, Ny,a is the number of lobsters at age a in year y, My is natural mortality in year y, γy,a is the ratio that gauges the extent to which the seasonal pattern in catch rate affects the population analysis (Mertz and Myers, 1996) and Cy,a is the number of age a lobsters caught in year y.
Among the three age groups in Torres Strait, the fishery targets two age groups, 1-and 2-year-olds. These two age groups are sub-adults and their natural mortality rates are not expected to be very different. However, due to highly variable environmental conditions, annual natural mortality rates are believed to be highly variable, ranging from 0.17 year-1 to 1.5 year-1 (Pitcher et al., 1997). Therefore, Eq. 1 assumes that the natural mortality rate varies between years, but remains constant over age. The model can assume that natural mortality varies with both year and age, but the number of parameters to be estimated will become too large to obtain reliable estimates.
γy,a is determined by the seasonal fishing pattern and varies with year and age (See Mertz and Myers, 1996). However, by assuming that fishing takes place in a pulse form in the middle of a fishing season that starts at time t (0 ≤ t ≤ 1) and lasts for T (0 ≤ T ≤ 1) within a year, we can simplify the calculation of γy,a to the following:
| (2) |
This model is known to be insensitive to the assumption of pulse fishing (Mertz and Myers, 1996). In fact, by specifying a pulse fishery at midyear, Eq. 2 immediately yields γy,a = γ = e-M/2 and after substitution into Eq. 1, one recovers Pope’s (1972) cohort equation.
The lobster fishery has had a minimum legal size limit in place since 1987. This regulation generally prevents age 1 lobsters from being captured before August. So, the fishing season for age 1 group is from September to December each year, meaning T1 = 0.33 and t1 = 0.67. Age 2 lobsters are subject to fishing until they migrate out of the Torres Strait fishing grounds in September; therefore T2 = 0.67 and t2 = 0.
After estimating the stock numbers at age, solving the following equation yields fishing mortality rate Fy,a:
| (3) |
where, Ta is the fraction of fishing season during the year as specified above.
Equation 1 is an extended version of the cohort equation (Pope, 1972), which is normally used to describe the dynamics of a single cohort. There are two major groups of methods in applying the cohort dynamic equation or its general catch equation to catch-at-age data. The first group of methods, commonly called Virtual Population Analysis (VPA), (Gulland, 1965; Laurec and Shepherd, 1983) or cohort analysis (Pope, 1972), are recursive algorithms that calculate stock size based on catches with no underlying statistical assumptions. The second group comprise methods that rely on formal statistical models and are normally called statistical catch-at-age methods (Doubleday, 1976; Paloheimo, 1980; Fournier and Archibald, 1982; Megrey, 1989). The recursive methods used in VPA or cohort analysis work well when there are many age groups. Their use would be doubtful for a fishery of only two major age groups. Also, VPA or cohort analyses require the assumption that the natural mortality is known. In contrast, the statistical catch-at-age methods are able to estimate natural mortality. For these reasons, statistical catch-at-age methods are used in this study. The lobster fishery has no fishing effort data available and therefore, the cohort equation (Eq. 1) was used to describe the dynamics.
Paloheimo (1980) and Doubleday (1976) allow the recruitment for each cohort to be a totally free parameter to be estimated. We, however, force the estimated recruitments fit some stock-recruitment relationship (Fournier and Archibald, 1982). This has a big advantage in reducing the parameters to be estimated for the lobster fishery. The lobster fishery stock has only two major age groups. For a 12-year catch-at-age matrix, a normal catch-at-age analysis needs to estimate 12 recruitments plus the number of age 2 lobsters in the first year. To reduce the number of parameters to be estimated, recruitment is modelled stochastically and the recruitment deviations are parameters of the model.
Spawning capacity of a stock is better measured by spawning biomass than by number of individuals because egg production of most fishes and marine invertebrates is more closely related to biomass (Beverton and Holt, 1956). With population numbers at age, spawning biomass is calculated as the average of age 2 lobsters that survived to October, November and December (the breeding months),
| (4) |
where,
| (5) |
The
| (6) |
where, TW is tail width (mm) and
where, L is carapace length calculated from the von Bertalanffy growth equation (Trendall et al., 1988)
| (7) |
The stock-recruitment relationship is defined by Ricker’s model
| (8) |
where, α is the recruits-per-spawner at low stock sizes, β is a parameter
describing how quickly the recruits-per-spawner drop as S increases and τy
is an error term assuming independent normal random variates with a zero mean
and unknown variance
Observation Equation
The lobster survey was carried out in May/June each year. Abundance indices
for each lobster age group were constructed using generalized linear models
(Ye et al., 2005). There are a variety of choices for modelling the distribution
of the observed relative abundance index Iy,a. The lognormal distribution
is generally considered a sensible choice for a relative abundance index,
| (9) |
where, NIy,a are the numbers-at-age at the time of the survey and calculated by Eq. 5, εy,a is an error term describing observation errors, which are assumed to be independently and normally distributed with a mean of zero and variance of σ2ε for each age group and q is a scale parameter that makes the survey index and stock abundance comparable and should not be confused with the catchability coefficient used to define the relationship between fishing mortality rate and effort. No explicit catchability coefficients are used in this study because we do not use fishing effort at all.
The catch data are assumed to have measurement errors, ξ, that are normally distributed with a mean of zero and unknown variance σ2ξ,
| (10) |
where, Cy,a is the estimated catch and Ĉy,a the observed catch in year y and age a.
The Bayesian Approach
In the Bayesian approach, one proceeds by defining the probability distribution
(or likelihood function) for the observation y, that is, the observations in
Equations 9 and 10. Here catches, referred
to as control variables in state-space models, are assumed to be known with
random errors. The probability model can be constructed as follows
| (11) |
where, the parameter set θ refers to all of the unknown parameters given above, except for the variance terms. These unknown parameters are further assumed to be random quantities with prior distributions π(θ|η), where η is a vector of the variance terms, generally referred to as hyperparameters in our model. The prior distributions reflect our knowledge of the state of nature prior to collecting the data. Having collected the data, we update our view of the state of nature and express this update as posterior distributions.
| (12) |
Based on the posterior distributions, we can summarize our updated knowledge about the elements of θ in terms of means, medians, etc. The integral in the denominator of Eq. 12 is a normalization factor and involves high-dimensional integration (equal to the dimension of the state μ). In the past, this was the main difficulty that hindered the application of Bayesian inference (Dowd and Meyer, 2003). For modern posterior computation using simulation techniques, this normalization constant is not required. Equation 12 shows how the prior joint probability density function (pdf) of the parameters is updated using the likelihood function p(y|θ)to yield the posterior pdf of the parameters. The posterior pdf is thus regarded as a synthesis of different sources of information.
The general solution for p(θ|y, η) relies on Markov Chain Monte Carlo (MCMC) methods. There are two major techniques devised to create Markov chains with the desired stationary distribution. They are Gibbs sampling (Gelfand and Smith, 1990; Carlin et al., 1992; Carlin and Louis, 1996; Meyer and Millar, 1999) and the Metropolis-Hasting algorithm (Metropolis et al., 1953; Hastings, 1970). The former was used here and Bayesian modelling was carried out here using the windows version of the public domain package WinBugs (Speigelhalter et al., 1995).
Prior Distributions
Prior distributions are required for all unobserved quantities in the Bayesian
implementation of this model. These priors encompass what we know about the
population from previous studies or from related or similar species. In the
absence of prior knowledge, it has been common practice to seek an appropriate
non-informative prior distribution, that is, letting the likelihood function
for the observation dominate estimation of the posterior distribution (Carlin
and Louis, 1996). Intuitively, this would suggest using a uniform prior distribution
thus giving each possible value of any parameter an equal chance of being chosen.
Unfortunately, uniform distributions are not invariant under reparameterization
and what might be non-informative on one scale may not be on another (Smith
and Lundy, 2002).
Alternatively, Box and Tao (1973) suggest choosing a prior that is diffuse enough that the data will dominate whatever information there is in the prior. The search for priors is far from clear cut. Kass and Wasserman (1996) review more than a dozen formal methods for deriving priors and find that none of them are guaranteed to produce good priors. Millar (2002) gives details on specifying default priors for several common fisheries models. This is the general approach that was followed here.
Variance parameters are often excluded from the parameter vector θ and are presumed a priori independent from other parameters under the state of prior ignorance (Box and Tao, 1973; Kass and Wasserman, 1996). Fisheries data often have lognormal errors and logσ has a flat prior (Millar, 2002). We assume therefore that the error term for the relative abundance index has the following prior, logσε~Gamma(0.001,0.001). For the variance of catch in Equation 10, the priors are set at values based on the real catch data for each age group, σξ1 ~ Gamma (0.5,0.08) for the age 1 group and σξ2 ~ Gamma (0.5,0.007) for the age 2 group, respectively.
Catchability q is a scale parameter and Jeffreys’ prior applies to it (Millar, 2002). We model its variance using the inverse gamma distribution-the recommended form of prior for scale variables (Carlin and Louis, 1996; Meyer and Millar, 1999): 1/q~Gamma (0.001, 0.001). This form of prior approximates a Jeffreys’ prior which is both noninformative and invariant to changes in scale (Congdon, 2001).
The Jeffreys’ prior for α in Equation 7 is uniform on the log scale and that for β is uniform on [0,∞] (Millar, 2002). We then set 1/logα~Gamma(0.01, 0.01) and β~uniform(0, 1) for α and β. The prior for β was limited based on the previous study on stock and recruitment relationship (Pitcher et al., 2002).
The Jeffreys’ prior of instantaneous natural mortality rate is uniform on [0,∞] (Millar, 2002) and we use M ~ lognormal (-0.32, 5) to limit its value between 0.17 and 1.15 based on the results from Pitcher et al., (1997).
The estimates of recruitment in previous years give some information about the order of magnitude and variability in recruits. We, therefore, use the abundance estimates derived from the full-scale survey in 1989 and the following year’s monitoring survey to set priors for the recruitments of the first two years and the initial numbers-at-age for the first year; the priors are, N1,1~lognormal(6.7, 10), N1,2~lognormal(6.5, 10) and N2,1~lognormal(7, 10).
Monitoring Convergence of MCMC
Convergence in the context of MCMC methods means that the sampled values
adequately characterize the posterior distribution for any particular node or
variable. It is very difficult to say conclusively that a chain has converged,
but it is much easier to diagnose when it definitely has not. Tierny (1996)
presents the theoretical conditions for convergence but there is no one omnibus
test for convergence and authors generally recommend using a variety of tools
(Carlin and Louis, 1996).
| Table 1: | Summary of posterior distributions for model parameters. The columns labelled 0.025 and 0.975 refer to the lower and upper limits of the 95% density regions for the posterior distribution of the parameter |
There are four aspects of an MCMC run that need to be considered (Gelman, 1996; Smith and Lundy, 2002). The first is to determine how many iterations it takes before the Markov chain is stationary, i.e., how long the burn-in period should be. The second is to decide on the frequency of the thinning process (sub-sampling every kth iteration) to break the autocorrelation of the samples from MCMC methods in order to obtain independent and identically distributed observations for variance estimation. The third concerns how long the sampling should be conducted, that is, how many iterations are needed once the chain has converged. Finally, to ensure robustness and consistency in the results of the Monte Carlo integration, multiple realizations of the chain are generated using different starting values for the random variables as well as different random seeds.
We used a burn-in period of 1000 interactions and sampled every 20th iteration to break the autocorrelation. Approximately 2x104 quasi-independent samples of p(θ|y,η) were drawn from the simulated Markov Chain and used to construct estimates of the posterior pdfs of the parameters and states. After the trace and history plots of two different chains with WinBUGS suggested no evidence against convergence, we employed Gelman-Rubin statistics, as modified by Brooks and Gelman (1998) and the Raftery and Lewis (1992) convergence diagnostic to further check convergence. BOA (Smith, 2004) was used for these checks. The resulting Monte Carlo error for each parameter reached less than 5% of the sample standard deviation (Table 1), a rule of thumb criterion of convergence (Spiegelhalter et al., 2003).
Results
Model Parameter Estimates
A summary of the model parameters for the age-structured state-space model
is given in Table 1. Medians are also listed and they are
quite similar to the means. However, it is worth noting that some parameters
are transformed.
The posterior pdfs of the model parameters and initial state variables, recruitments in 1989 and 1990 and the number of age 2 lobsters in 1989, have been plotted along with their associated prior distributions for easy comparison (Fig. 1).
| Fig. 1: | Posterior distributions of the estimates of the model parameters and initial state variables contrasted with their associated priors. Dashed line indicates prior and solid line the posterior |
| Fig. 2: | Prior and posterior distributions of the standard deviations in the age-structured model. Sigma for log (index) = σε, sigma for log (recruitment) = στ, sigma for log (age 1 catch) = σζ1, sigma for log (age 2 catch) = σζ2. Dashed line indicates prior and solid line the posterior |
All these estimates show significant difference between posterior and prior distributions, except the natural morality rate in 2000, which closely resembles the prior pdf, with only a slight shift in the mode of the distribution. The additional information associated with the data has allowed the refinement of the posterior pdfs of most of the parameters. The mean estimates of natural mortality rate M ranged from 0.47 year-1 in 1998 to 0.77 year-1 in 1990 (Table 1 and Fig. 1). The maximum value was almost double the lowest estimate, suggesting a large inter-year variation in natural mortality. The 95% credible intervals are also fairly wide.
The likelihood was informative for the process variance term
| Fig. 3: | Bivariate posterior distribution of α and β that defines the stock-recruitment relationship |
| Fig. 4: | Comparison of age 1 lobster population size between the model estimated medians and survey relative abundance indices |
Most parameters listed in Table 1 are not highly correlated with each other, except α and β in Eq. 8. The bivariate plot of the joint posterior distribution highlights this correlation (Fig. 3). For dependent parameters, the distribution of the one variable changes as a result of drawing a value for the other variable. This modified distribution is the conditional distribution used in Monte Carlo simulations.
The estimates of age 1 lobsters in June, when the annual fishery independent surveys were conducted, were contrasted with those estimated from the relative abundance indices and scale parameter through Eq. 9. The estimated population size captures the variation seen in relative abundance indices quite well (Fig. 4). A large discrepancy between the survey and model estimates can be seen only in 1989.
Age 2 lobsters exhibited a declining trend consistent with the estimates from relative abundance indices (Fig. 5). In contrast with the age 1 group, the model-based estimate in 1989 is much lower than the survey index. Estimates for age 1 and age 2 groups are confounded. This is caused by the integrated approach to model fit for all age groups. In general 1989 is likely to be an outlier (Ye et al., 2004c).
The estimates of fishing mortality for age 2 lobsters show a generally increasing trend from 1989 to 1996, but decreasing thereafter (Fig. 6). The age 1 group had an overall annual variation pattern similar to that of age 2 lobsters. However, the relative fishing mortality rate of age 1 lobsters was only 20% that of age 2 lobsters in 1989. It increased steadily over time and reached 60% in 1999. Unfortunately, the 2000 fishing mortality for age 1 lobsters could not be estimated because Equation 3 requires knowing the number of age 2 lobsters in 2001.
| Fig. 5: | Comparison of age 2 lobster population size between the model estimated medians and survey relative abundance indices |
| Fig. 6: | Estimates of fishing mortality rates from 1989 to 2000. The horizontal line indicates the fishing mortality level required for the maximum sustainable yield as defined in Fig. 8 |
| Fig. 7: | Model-estimated stock recruitment curve. The squares and bars indicate medians and 95% creditable intervals of the recruitment estimates |
The fishing mortality on the age 1 group has risen steadily, which may indicate the industry has increased the targeting of age 1 lobsters over recent years.
The estimated stock-recruitment curve is shown in Fig. 7. The 95% creditable intervals for recruitment seem quite wide. The high variation reflects the fact that recruitment of lobsters is highly variable, like many other fisheries. The median estimates for recruitment cannot be considered as point estimates because the parameters in Fig. 7 were estimated based on the distributions of recruitment rather than just the mean.
| Fig. 8: | Relationship between fishing mortality (F) and sustainable yield (C in tonnes of tail weight), assuming M, α and β take their median values |
In this sense, estimation of the stock-recruitment model should be more representative and reliable.
The spawning stock biomass was the highest, 230 tonnes, in 1989, but dropped down to almost one fourth of that level, 60 tonnes, in 1996 (Fig. 7). Recruitment correspondingly decreased from a high of 11 million in 1989 to 6.5 million in 1996. The stock recruitment curve indicates that spawning stock biomass should be kept at around 150 t to obtain the maximum recruitment.
Productivity of the Lobster Fishery
For management purposes, it is always interesting to know the productivity
of a fishery. Under average conditions, the population produces surplus production
each year, which can be taken as equilibrium yield without changing the population’s
size (Quinn and Deriso, 1999). The corresponding strategy for managing fish
populations is then a surplus production policy (Chapman et al., 1962;
Ricker, 1975). One calculates the surplus production in a given year and then
sets a harvest below, equal to, or above this level depending on whether one
wishes to rebuild, keep stable, or reduce the current population size (Quinn
and Deriso, 1999). As fish population size is not directly measurable, surplus
production is often expressed as a function of fishing mortality or effort for
convenience. This can be easily done with production models. However, with age-structured
models, it is rather complex to derive analytical results. Numeric methods were
used to calculate the sustainable yields under different fishing mortalities
in this study.
If the natural mortality rate, α and β assume their posterior mode values, the sustainable yield of the fishery is simply a function of fishing mortality (Fig. 8). Sustainable yield increases as fishing mortality increases up to some point, at which point yield begins to decline with further increase in fishing mortality. Unlike the curve derived from Schaefer’s production model (Schaefer, 1957), this curve is not parabolic, but has a long concave tail on the right, which differs from the convex tail of the Pella-Tomlinson model (Pella and Tomlinson, 1969). This may indicate that the spawning migration out of the fishing grounds provides some kind of protection from over-fishing and that even under very high fishing pressure the stock can still maintain a certain level of production, rather than collapsing quickly, as observed in many other fish populations.
The Bayesian approach adopted in this study enables us to go beyond point-estimates and consider the uncertainties involved in the estimates of natural mortality and the two parameters defining the stock recruitment relationship. We carried out multivariate Monte Carlo simulation based on a univariate distribution of M and a bivariate distribution for α and β to estimate the sustainable yields under various fishing mortalities (Fig. 9).
| Fig. 9: | The sustainable yields (C) under various fishing mortality rates (F) taking into account the uncertainties in the estimates of natural mortality (M) and the stock recruitment parameters (α and β in Eq. 8). The marks, *, o and x, indicate 25, 50 and 75% quantiles, respectively |
The middle line indicates the median estimates of the sustainable yields. It should be noted that the maximum yield of the middle line in Fig. 9 is slightly lower than that of Fig. 8 that was deterministically estimated.
The other two lines in Fig. 9 are the 25% and 75% quantile estimates, respectively. Most noticeable is that the fishing mortality that produces a maximum sustainable yield differs between the quantiles. The 25% quantile line peaks at F = 0.2 year-1, the mode line has a maximum at F = 0.4 year-1 and the 75% quantile line at F = 0.5 year-1. The stochastic relationship shows that there is a large uncertainty in the sustainable yield even if the control variable, fishing mortality rate, is fixed at a specific level because of the potential uncertainties in natural mortality and stock recruitment relationship.
Discussion
The Fishery
Natural mortality is a key parameter in most stock assessment models and
also one of the most difficult parameters to determine, due to the confounding
effects of recruitment, fishing mortality and natural mortality (Quinn and Deriso,
1999). Pitcher et al., (1997) estimated annual natural mortality rates
of the lobsters in Torres Strait from 0.17 year-1 to 1.15 year-1.
As their estimates were based on two points of abundance data (e.g. number of
age 2 group in year t and number of age 1 group in year t-1), the
higher variation is quite likely due to data noise or measurement errors. Our
estimates were derived from fitting the population dynamic model to the catch
and survey data and the estimation was done in an integrated way, which prevents
the loss of information contained in the year-to-year variation. This estimation
should be relatively robust to measurement errors. However, Pitcher et al.,’s
(1997) estimates are still within our posterior distributions (Fig.
1).
The highest natural mortality estimate of 0.77 year-1 was seen in 1993 and 2000 (Table 1). This coincides with periods of seagrass dieback recorded in 1992/93 and 1999/2000 (Ye et al., 2004b). However, it must be noted that the difference in natural mortality between 1993 and neighbouring years and between 2000 and 1999 is not very large (Table 1). A potential reason for this is suggested by two observations. First, the seagrass dieback was recorded mainly in the northwest area of Torres Strait, close to the freshwater discharge along the PNG coast that is implicated as the cause of the seagrass mortality. Second, the lobsters in the northwest area are normally found at a relatively low density. So, the impact of seagrass dieback on lobster mortality may be local and limited only to the northwest area. The natural mortality estimated in this study is for the entire population and therefore, the limited increase in natural mortality may reflect the restricted area of impact.
The recruitment estimates exhibit wide credible regions (Fig. 7), suggesting high uncertainties in recruitment, i.e., a relatively weak stock-recruitment relationship. However, the estimated stock-recruitment relationship is similar to that based on the estimates of age 1 and 2 lobsters directly derived from the May/June surveys from 1989 to 2002 (Ye et al., 2004c). The new age-structured model used the mean biomass of the age 2 lobsters in September, October and November as a measure of the stock’s spawning capability (Eq. 4) and the recruitment was defined as the number of age 1 lobsters at the beginning of each year. The survey-based estimates of age 2 lobsters in May/June do not directly contribute to estimates of spawning as the breeding season is from November to February (MacFarlane and Moore, 1986). The four-month period between the survey and spawning could contribute a great uncertainty to spawning- stock biomass depending on fishing intensity. Equation 4 should provide a better measure for spawning stock biomass. However, the similarity of Fig. 7 to the previously derived stock recruitment relationship of Ye et al. (2004c) suggests that the improvement in measurement of spawning stock may offer little reduction in the uncertainty of the stock recruitment relationship in this case.
The estimated fishing mortality rate of the lobster fishery exceeded 0.4 year-1 in 1994, 1996-97 and 1999 and was very close to that level in 1998 (Fig. 5). This certainly suggests that the fishery was under high fishing pressure during these years and might well have been overfished. This conclusion is supported by the estimates of low spawning stock biomass (Fig. 7) for that period. Although the abundance estimates for age 1 and 2 lobsters were relatively low in 1999 and 2000 (Fig. 4 and 5), the difference from estimates for other years is not clear cut. However, the relative abundance indices estimated from the annual surveys show very alarming signs in 2001 and 2002 (Ye et al., 2004c).
The stochastic relationship between sustainable yield and fishing mortality (Fig. 9) quantifies uncertainty in the response of three major variables (M, α and β) after a complex transfer function of sustainable yields. Recruitment models are widely used because strategic fishery management demands an assessment of the likely future response of fish populations to exploitation (Needle, 2002). Assessment of uncertainty in recruitment models can lead to probabilistic prediction of management consequences and the scientific formation of realistic regulations. It is undeniable that both natural mortality and stock-recruitment relationship vary from year to year, likely due to variation in factors related to larvae settlement and the survival of juveniles and lobsters recruited to the fishery. The natural variation in these parameters and the uncertainties in their estimates make the estimation of sustainable yield under a specific fishing mortality rate very uncertain. For example, the 50% credible interval of the maximum sustainable yield ranges from 110 t to 390 t and its corresponding fishing mortality rate from F = 0.2 year-1 to F = 0.5 year-1 (Fig. 9). The uncertainty in sustainable yield can be decomposed into the contribution of each single factor by drawing only one factor from its posterior distribution. This can help to pinpoint the factor of the greatest contribution to uncertainty. Although sustainable yield is a highly hypothetical concept, it has been used extensively in fishery management. The stochastic sustainable yield provides more information about the uncertainty involved and should improve risk management in fisheries.
Fishery management is about choosing harvest strategies to achieve the goals and objectives pre-set for a specific fishery. However, both the population dynamic model and the decision process involve stochasticity. To improve management, fishery managers need to assess the risks associated with alternative management actions. The posterior distributions of the model parameters and states produced by the Bayesian approaches can easily be used for quantitative risk analyses (Francis 1992; Francis and Shotton 1997). The results of this study should facilitate the harvest strategy evaluation of the Torres Strait lobster fishery in the near future.
State-space Models and Bayesian Approaches
In state-space models, the current state of the process summarizes all of
the information from the past that is necessary for predictions of future states
(Abraham and Ledolter, 1983). The measurement equations (Eq. 9
and 10) describe the generation of observations from the
current state with measurement error. The system equations (Eq.
1 and 7) describe how the states are supposed to evolve
through time with random shocks captured by the process error. Models without
process error imply deterministic state trajectories and time-invariant parameters
(Millar and Meyer, 2000). Deterministic models predict stock size into the infinite
future as a deterministic function of the current stock and harvesting policy.
In general, the measurement error characterizes our uncertainty with respect
to the current state while process error captures the uncertainty moving from
the current state to the next one in time (Smith and Lundy, 2002). State-space
models are more realistic and do have superior performance compared with deterministic
models. However, classical likelihood-based statistics restricts the applicability
of state-space methodology. The MCMC techniques have provided solutions to problems
that were hitherto considered computationally intractable and will almost surely
have an impact on fisheries stock assessment (Meyer and Millar, 1999).
The ability to account for both measurement and process error has made state-space models more and more popular in fisheries stock assessment over the last few years (Meyer and Millar, 1999; Millar and Meyer, 2000). This study provides an extension to the existing application by considering the simultaneous joint estimation of both the stock-recruitment relationship and other model parameters, which can also reduce the number of parameters to be estimated in the age-structured model of a stock having a low number of age groups.
Deriso et al. (1985) found it was very helpful to include a stock-recruitment curve. However, most age-structured models serve to produce estimates for recruitment and spawning stock that are then used to build a stock-recruitment model. The integrated single step method used in this study may have some advantage in preventing loss of information and allowing the stock-recruitment relationship to influence the estimation of model parameters. The traditional two-step sequential methods may help break the confounding between years or parameters. Which method performs better is likely to depend on the case. Nevertheless, we should try the two-step method and compare the resulting parameter estimates with those estimated from one-step approaches in the future.
Bayesian models require the specification of prior distributions for all the parameters to be estimated. The specification of prior information is, however, a contentious issue. The formal methods derive the prior as the square root of the determinant of the information matrix obtained from the likelihood function of the data (Jeffreys, 1961). There are limitations of using the formal methods. The prior may need to be chosen by consensus when formal methods fail. The inverse prior on scale parameters is widely accepted by Bayesian modellers as the appropriate default prior (Millar, 2002). Other priors may be derived based on real data or information of existing studies. For example, the priors for the natural mortality rates used in this study were based on previous study from Pitcher et al. (1997). However, it must be noted that there is no one way for the choices of prior assumption. Sensitivity to prior assumptions should be assessed thoroughly. We have repeated the analysis under different prior assumptions, but with the same situation in order to aid direct comparison of results.
Acknowledgements
This research was supported by AFMA Torres Strait research funds and by CSIRO. We thank Nick Ellis and Wayne Rochester for their valuable and constructive comments on early drafts of the manuscript.