© 2005 International Council for the Exploration of the Sea
A framework to model shrimp (Pandalus borealis) stock dynamics and to quantify the risk associated with alternative management options, using Bayesian methods
Pinngortitaleriffik – Greenland Institute of Natural Resources PO Box 570, DK-3900 Nuuk, Greenland
*Correspondence to C. Hvingel: tel: +299 361200; fax: +299 361212. e-mail: hvingel{at}natur.gl.
A new integrated Bayesian framework for making quantitative assessments, predictions, and risk analyses of shrimp (Pandalus borealis) stock development is constructed. A biomass dynamic model, based on the logistic function but including an explicit term for cod predation, suggests that the quantity of shrimp consumed by cod could equal that taken by the fishery. The model proved superior to an alternative model in its ability to estimate parameters central to the assessment; the alternative model subsumed cod predation as part of an overall population growth effect without a time trend. Two series of shrimp biomass indices, catch, cod biomass estimates, cod consumption estimates, and prior distributions of model parameters provided information to the models. Process and observation errors were incorporated simultaneously using a state-space modelling framework. A Bayesian approach was used to construct posterior probability distributions of model parameters and derived variables relevant for management advice, including quantification of future risk of transgressing reference points in relation to alternative management options.
Keywords: assessment, Bayesian, management advice, model, predation effect, risk, shrimp
Received 11 June 2004; accepted 19 June 2005.
 |
Introduction |
|---|
 Top Introduction Model developmentResultsDiscussionAppendixReferences |
|---|
Annual landings of northern shrimp (Pandalus borealis) in the North Atlantic increased from about 100 000 t in the early 1980s to about 400 000 t in 2002. Although shrimp are now one of the most important target species in these waters, little progress has been made in developing a standardized predictive tool for stock assessment (cf. Stefánsson et al., 1994; Cadrin and Clark, 1999; Hvingel and Kingsley, 2000; Koeller et al., 2000; ICES, 2001; NAFO, 2001).
Management advice for most shrimp stocks in the North Atlantic is basically formulated by qualitative assessment of trends in various indices of stock condition in response to the catch history (Koeller et al., 2000; ICES, 2001; NAFO, 2001). Typically, biomass estimates and length compositions from research surveys, plus series of commercial catch rates, constitute the main data source, but additional observations may also be considered, such as predator abundance and temperature. Advice is given as an annual Total Allowable Catch (TAC) or as a statement about the sustainability of current fishing practice as authorized by the assessment board. One recent extension of the qualitative assessment procedure is the "Traffic-light" approach (Koeller et al., 2000). This multi-indicator system uses a systematized list of available indicators whose values are qualitatively categorized as red (bad), yellow (neither bad nor good), or green (good). A summation of the status of all the listed indicators provides an overall status of the stock, which is then used as guidance for the next year's harvest strategy. For either of the qualitative assessment types, the method of deriving the advice is not explicitly stated, and the uncertainty associated with the process is not quantified. Such methods lack predictive rigour, including formal statements of uncertainty, and are therefore not suited to quantitative comparisons between alternative management options.
A quantitative assessment may be reached through construction of mathematical models describing stock dynamics and their links to data and ancillary information. Following the ideas of Punt and Hilborn (1997), the first step in the construction of an assessment framework is to specify alternative hypotheses. In this paper we explore logistic biomass models because on the basis of the data normally available for shrimp stock assessments, they seem suitable. Additional models can be added to the assessment framework described here, to account for model uncertainty. Process and observation errors are incorporated using a state-space modelling framework (Schnute, 1994), and Bayesian inference (Congdon, 2001) is used to estimate probability distributions of model parameters.
Age-structured population models that account in great detail for age-dependent growth, growth-dependent recruitment to the fishery, and age-specific fishing and natural mortality are now standard tools for assessing most fish stocks. However, shrimp have no permanent hard parts to accumulate annual layers, so ages cannot be measured. Age distributions have to be constructed indirectly by some type of modal analysis of length distributions (e.g. Macdonald and Pitcher, 1979), but the results are dependent on estimates of growth rate, which generally cannot be extracted from the data and have to be assumed. This is particularly true for the relatively slow-growing shrimp of northern waters and especially for the important larger size classes. The use of standard VPA software (e.g. the Lowestoft VPA suite: http://www.cefas.co.uk) to assess shrimp stocks has not been successful (Savard et al., 1991; ICES, 2001; NAFO, 2001). Length-based models are currently not used for shrimp stock assessment (Fournier et al., 1991; Quinn et al., 1998).
Models without demographic structure are an alternative. Such models are usually of a stock-production type (Pella and Tomlinson, 1969), and they describe stock dynamics simply in terms of rates of change of total biomass, rather than by details of individual, age-specific growth and mortality (see Hilborn and Walters, 1992, for a review of production models). However, natural mortality is thought to be at least as important as fishing mortality in the dynamics of shrimp stocks, and this limits the application of traditional fishery models, which assume natural mortality to be constant and/or partly a function of stock density over time. Explicit factors in natural mortality, e.g. predation and influences of the physical environment, should be considered when constructing shrimp stock assessment models.
Several species of fish and marine mammal prey on shrimp (Magnússon and Pálsson, 1991a; Pedersen and Riget, 1993; Grundwald, 1998; Nilssen et al., 2000). Atlantic cod (Gadus morhua) cooccur with northern shrimp, are bottom-feeding predators, and can attain high biomass densities, qualities that give them the potential to be a controlling predator (Lilly et al., 2000, and references therein). The rapid development of the West Greenland shrimp fisheries during the 1970s is thought to be connected to a major decline in the cod stock in that area; shrimp fisheries off the Canadian east coast appear to have, similarly, thrived after the collapse of the northern cod stocks in NAFO Divisions 2J + 3KL (Lilly et al., 2000). In the Barents Sea, cod are estimated to consume 3–4 times as much shrimp as is taken by the fishery (NAFO, 2001), and in northern Icelandic waters consumption estimates since 1983 range from 0.25 to 4.6 times the catches (H. BjÖrnsson, M. R. I. Reykjavík, pers. comm.). Therefore, for several shrimp stocks, predation by cod may constitute a dominating component of mortality and should be considered for explicit inclusion in any shrimp stock assessment model. Other predators of shrimp may also be important for particular stocks, and could be included as well if their effect generates a sufficiently strong signal, and if data are available.
The physical environment also affects shrimp populations (Anderson, 2000; Koeller, 2000) and may be the cause of the rapid changes in abundance seen in some stocks (Apollonio et al., 1986; Anderson, 2000). However, at this stage of the modelling, environmental effects were not considered as explicit variables.
This paper presents an integrated framework for assessment and management of a shrimp stock that can use all the available information on the stock dynamics, provide realistic estimates of the uncertainty associated with the assessment conclusions, and convey the information in a usable form to fishery managers.
|
Model development |
|---|
|
TopIntroductionModel developmentResultsDiscussionAppendixReferences |
|---|
Modelling framework
The models were built in a state-space framework (see, for example, Schnute, 1994; Meyer and Millar, 1999, and references therein). This time-series method offers straightforward implementation of process and observation error, and great flexibility in the mathematical construction approximating stock dynamics and data–stock relations.
The state-space conception of a stock-dynamic model regards a series of unobserved state values as defining the stock trajectory through time. States are related to one another and to the data by simultaneous, usually stochastic, equations. The models were therefore defined in terms of a set of parameters (
), the values of which defined the dynamics of the shrimp stock.
The posterior distribution for the parameters of the model, p(
|data), given a joint prior distribution, p(), and the likelihood of the data, p(data|), is determined using Bayes' (1763) theorem:
|
|
In applying Bayes' equation to the present problem, the posterior probability distribution of is derived by Monte-Carlo–Markov-Chain (MCMC) sampling methods (see Congdon, 2001). The programming framework WinBUGS v.1.3, made available by the Medical Research Council and the Imperial College of Science, Technology and Medicine in England (http://www.mrc-bsu.cam.ac.uk/bugs; Gilks et al., 1994; Spiegelhalter et al., 2000), provided a means of specifying and analysing a Bayesian model, including selection and implementation of appropriate algorithms; for numerical integration WinBUGS uses "Metropolis-Hastings within Gibbs sampling" (Gilks et al., 1993, 1996; see also www.mrc-bsu.cam.ac.uk/bugs/).
State equations
The basic equation was a generalization of the logistic model of population growth (Pella and Tomlinson, 1969):
|
| (1) |
A discrete form of this model, modified to include fishing mortality and predation by cod, and parametrization in terms of MSY (Maximum Sustainable Yield) rather than r (intrinsic growth rate) (cf. Fletcher, 1978) described the state transition from time t to time t + 1:
|
| (2) |
Absolute biomass estimates provided by population-dynamic models are liable to great uncertainty if no explicit information is available to scale the biomass indices to real stock size. Therefore, for management purposes it is desirable to work with biomass on a relative scale in order to cancel out the uncertainty of the "catchability" parameter (the scaler). This was accomplished by dividing Equation (2) by BMSY, the biomass that produces MSY. The variability (ratio of the inter-quartile range to the median) of estimated biomass ratios was about 67% lower than that of absolute estimates of B. This reparameterization also reduced autocorrelation in the chains of values sampled by the Gibbs sampler and thus hastened convergence to the posterior distribution (Meyer and Millar, 1999). Finally, a term for the process error was applied, and the state equation took the form
|
| (3) |
, is normally, independently and identically distributed with mean 0 and variance 
2. The consumption, Ot, was taken as a Holling type III functional-response function that includes prey-switching (Holling, 1959):
|
| (4) |
Under the basic stock-dynamic model, the biomass at MSY is given by
|
| (5) |
|
| (6) |
In the absence of predation, these are also the MSY available to a fishery and the corresponding standing stock. Under prey-biomass-dependent predation, the fishable yield will be less if the predator switches to this prey at prey biomass levels below BMSY. However, the standing stock giving maximum fishable yield may also be altered (reduced) if prey-switching would occur below BMSY, because it might in that case pay to fish the stock to a lower biomass, to discourage the predator from switching and thereby to reduce predation. The maximum fishable yield will, however, always be equal to at least the MSY as defined above, minus cod consumption at BMSY, and will never occur at stock levels higher than BMSY.
Observation equations
Two independent series of shrimp biomass, one of shrimp catches, one of cod biomass, and data on consumption of shrimp by cod provided the data input to the model (Table 1). The two series of shrimp biomass indices were a standardized and combined series of annual commercial vessel catch rates for the years 1976–2002, cpuet (Hvingel et al., 2000; Hvingel, 2002); and a trawl-survey biomass index for 1988–2002, survt (Carlsson et al., 2000; Kanneworff and Wieland, 2002). These indices were scaled to true biomass by catchability constants, qc and qs. Lognormal observation errors, 
and 
, were applied, giving
|
| (7) |
|
| (8) |
and , are normally, independently and identically distributed, with mean 0 and variance 2 and 2. The error for the final year, N, of the cpue index was assumed to be 1.5 times the error for the rest of the series, because this data point is an interim one based on fishery data until October (the annual assessment takes place in November). Likewise, the first year of the survey was assigned a 50% larger error than the rest of the series, to allow for the process of learning to conduct the survey.
|
Total reported catch in NAFO Subarea 1 + Division 0A 1970–2002 (Hvingel, 2002) was used as yield data (Table 1). The fishery being without major discarding problems or variable misreporting, reported catches were entered into the model as error-free.
Predation was entered into the model with two sets of data. The first was a time-series of cod biomass, codt. It was constructed as follows: absolute estimates of stock size for West and East Greenland combined were available from VPA up to 1992 (ICES, 1996). For the period 1982–1992 the West Greenland components could be estimated using the annual relative East/West survey-based biomass indices (mean and standard error: 0.6 ± 0.24), and the West Greenland series was extended back in time by multiplying the total biomass estimates for earlier years by 0.6, the 1982–1992 average West Greenland fraction. The absolute estimates for 1993 and later were obtained from the survey indices by applying the 1990–1992 mean ratio of the VPA estimate to the survey biomass, assuming a constant relationship between them.
The second series was a partial series of shrimp consumption (Table 1). The estimates of annual consumption of shrimp by cod, 
, were based on stomach samples for the years 1989–1992 (Table 26 of Grundwald, 1998). These estimates were linked to the dynamics of the shrimp stock through Equation (4) and the cod biomass series. An observation error term, 
, normally, independently and identically distributed with mean 0 and variance 2, was applied:
|
| (9) |
Priors
Bayesian philosophy considers that an observer maintains a model, perhaps mental or conceptual, of reality that is subject to being modified (updated) by observations. As a quantitative version of this, Bayesian statistics considers that quantitative observations (data) can be used to update pre-existing probability distributions of the values of parameters defining a quantitative model. In such a discrete updating process, the prior distributions pre-date and are therefore independent of the study that furnishes the data on which the updating is based. The prior distribution for a parameter should incorporate all information already available, but if none can be identified, a low-information "reference" prior (Kass and Wasserman, 1996) is used. The following paragraphs describe the prior distributions constructed for the parameters of this model and the reasoning underlying them.
There is usually little information about the pristine size of a stock, but it is inescapable that some assumption about it is made. It is commonly assumed that B1 = K, i.e. that the stock was at carrying capacity when records started (corresponding to P1 = 2 as a prior on P1 for a logistic model). This assumption is questionable under any circumstances, but especially so when, as in the present case, the stock had already been fished for some years before the start of the observations and was also subject to predation. Instead, we used a much less informative lognormally distributed prior on P1. For the model without a predation effect, we chose P1 
log N(0.44,0.2)1, a wide distribution with a mode at about 1.5; in other words, we considered that the stock might have been fished down halfway from carrying capacity to MSY level by the time record-keeping started. In the model with predation included, we chose P1 log N(–0.05,0.2). This is also a wide distribution, but it has a smaller mode at approximately 0.9, acknowledging the likely effect of a large cod stock before the start of the data series (S. A. Horsted, Pinngortitaleriffik, Greenland Institute of Natural Resources, Nuuk, pers. comm.).
The prior distributions for the error terms associated with the biomass indices were assigned inverse gamma distributions (the gamma distribution, G(r,µ), is defined by: µrxr–1e–µx/
(r); x > 0), because error standard deviations typically follow this kind of distribution. Existing analyses indicated an estimated CV of around 10% for the cpue series (Hvingel et al., 2000) and about 17% for the survey series (Carlsson et al., 2000). Their standard deviations were therefore given inverse gamma distributions with modes at 0.10 and 0.17. The error parameter for the cod consumption series, Ot, was assigned an inverse gamma prior with a mode corresponding to a CV of around 25%, as a result of the expected larger uncertainty associated with the construction of the series.
The data could not be expected to contain much information about Vmax. An informative prior was therefore constructed using estimates of growth and food conversion coefficients for cod fed to satiation in laboratory experiments (BjÖrnsson and Steinarsson, 2002). Assuming an average cod population structure for the years 1975–1989 (ICES, 1996), a level for Vmax was estimated as 3 kg-shrimp kg-cod–1 y–1. The prior was given a relatively tight distribution of the form Vmax N(3,0.1).
Low-information or reference priors were given to MSY, qc, qs, m, K, , and P50%, because we had little or no information on how their probability distributions might look. MSY was given a uniform prior between 0 and 1 000 000 t, the upper limit sufficiently high not to truncate the posterior distribution.
The catchabilities qc and qs are confounded with the carrying capacity K. A uniform distribution was therefore not uninformative, and a prior distribution uniform on a log scale was preferred as a reference prior (Gelman et al., 1995; Punt and Hilborn, 1997; McAllister and Kirkwood, 1998). A similar distribution was used as a slightly informative prior for K, to discourage sampling of very high, unrealistic, values.
A uniform prior for m would not be uninformative with respect to the location of BMSY relative to K. To make a reference prior for the ratio BMSYK–1 in terms of m, an approximate function was used for the conversion. No analytical expression for m in terms of BMSY/K exists, so m as a function of BMSYK–1 in the interval 0.37 < BMSYK–1 < 0.63 was approximated by m = a exp(b(BMSYK–1)c) with appropriate values of a (0.1817), b (5.1174), and c (1.0938). BMSYK–1 was then assigned a uniform distribution on the interval 0.37–0.63. Values of m < 1 correspond to stock trajectories with non-asymptotic lower limbs and were considered unrealistic; they equate to values of BMSYK–1 below approximately 0.37. The upper truncation point was chosen to make the prior symmetric around BMSYK–1 = 0.5.
P50% will have an important effect in setting the predation rate, V(P), because Vmax is supposed to be relatively well determined. A prior for P50% that was approximately non-informative with respect to the predation rate was constructed by using a rearranged version of the predation rate function (see Equation (4)):
|
| (10) |
Convergence diagnostics
Conclusions based on an MCMC output depend on the assumption that the chain of sampled values for each parameter can be considered an unbiased random sample from the target distribution, i.e. the posterior. A chain is usually serially correlated, and at the start is affected by the initial values assigned to parameters, which are not random samples from the posterior distributions. The influence of initial values is avoided by discarding a long stretch (several thousand iterations) at the start of the chain, retaining values only after sampling is expected to have converged to stationarity, and autocorrelation is reduced by retaining only every n(th) member of the chain, n typically being of order 102. To check that these measures were effective, the Bayesian Output Analysis (BOA) programme v.0.5.0 by Brian Smith (http://www.public-health.uiowa.edu/boa/) was used in S-Plus 2000 (http://www.mathsoft.com/splus).
A number of parallel chains with different starting points and random number seeds were analysed by the Brooks, Gelman, and Rubin convergence diagnostic (Gelman and Rubin, 1992; Brooks and Gelman, 1998) to evaluate if the samples could be considered to have arisen from the target distribution. A stationarity test (Heidelberger and Welch, 1983) was applied to individual chains. If evidence of non-stationarity is found, iterations were discarded from the beginning of the chain until the remaining chain passed the test. Raftery and Lewis (1992) tests for convergence to the stationary distribution and estimation of the run-lengths needed to accurately estimate quantiles were used, and finally the Geweke convergence diagnostic was applied (Geweke, 1992).
Model check
In order to check whether the model was a good fit to the data, different goodness-of-fit statistics were computed. First, we calculated the simple difference between each observed data point and its trial value in each MCMC sampling step. The summary statistics of the distributions of these residuals indicated, by their central tendency, whether the modelled values were biased with respect to the observations. Second, the overall posterior distribution was investigated for potential effects of model deficiencies by comparing each data point with its posterior predictive distribution (Posterior Predictive Checks; Gelman et al., 1995, 1996). Simulated sets of observed data, datarep, were drawn from the sampling distributions for the observations as output from the model. Therefore, datarep has distribution
|
| (11) |
In this expression, the term p(datarep|) represents the sampling of observations from the distributions set up in the model to define them, and the term p(dataobs) represents the sampling of the parameters of those distributions from their own posterior distributions. If the model fitted the observed data well, the observed data should, collectively, be a "likely" drawing from its simulated distribution, i.e. the observed data and the replicate data should look alike. The degree of similarity between the original and the replicate data points was summarized in a vector of p-values, calculated as the proportion of n simulations in which a sampling of the posterior distribution for an observed parameter exceeded its input value:
|
| (12) |
Third, the Conditional Predictive Ordinate (CPO; Gelfand and Dey, 1994) was calculated as a harmonic mean of the likelihood:
|
| (13) |
Derived parameters and risk calculations
Probability distributions of quantities such as mortality or other statistics relevant for the assessment/management procedure, but not readily available as model parameters, may also be generated by the MCMC sampling process if they can be derived from the existing parameters of the model. This is done simply by adding the appropriate equations to the model. The mortality caused by cod predation and the fishery, Z, is scaled to ZMSY (the combined fishing and predation mortality that yields MSY) for the same reasons as relative biomass was used instead of absolute. The equations added for generating posterior distributions of the Z-ratio were
|
| (14) |
As the MCMC sampler (when converged) is drawing samples from the probability distribution of model parameters or derived quantities, the probability of their being smaller or larger than certain reference values may readily be derived. Such reference levels may be of stock or fishery status defined as guidelines of management, e.g. to comply with the concept of the precautionary approach (FAO, 1996). As stated above, ZMSY as defined here is sustainable, and corresponds to a stock size equal to or greater than that at which the fishable yield is maximized. The Z-ratio therefore gives a conservative view of the prevailing management regime.
For example, to determine the risk of transgressing a reference limit in 10 years under a given regime of catch and cod predation, one merely adds the values pertaining to the management scenario under investigation to the data input series of catch and cod biomass. The risk is then simply the relative frequency with which the MCMC-sampled values transgress the reference points. Using biomass, B, as an example, the risk of biomass in year t being below the reference can be calculated by:
|
| (15) |
Alternative model runs
Two model versions were investigated. Model 1 had no predation by cod and was fitted to the data series on cpue, survey, and catch. This model was similar in concept to that fitted by maximum likelihood methods in Hvingel and Kingsley (2000). Model 2 corresponded to the full model with predation, as defined above. It used the additional data sets on consumption and cod stock biomass.
In each run, BUGS was set up to do 550 000 iterations, recording only every 50th to compensate for serial correlation within the chains of some of the sampled parameter values. The run-time was about 1–2 h on a 1400-MHz laptop, depending on the model and the management questions investigated. The first 1000 samples of the recorded chain were discarded for "burn-in", leaving 10 000 samples as the final result. Convergence diagnostics were calculated for the parameters listed in Table 2 to confirm that the model and sampling set-up was appropriate to ensure convergence.
|
|
Results |
|---|
|
TopIntroductionModel developmentResultsDiscussionAppendixReferences |
|---|
Model 1 (without cod predation)
Judged from the flat shape of the generated posterior distribution, the data did not contain much information about K (Figure 1). A mode was evident at around 1.5 million tonnes, but the prior used, truncated at 6 million tonnes, interfered with the posterior. Fisheries data might determine MSY reasonably well, but often cannot specify whether the yield comes from a stock that is close to or far from its carrying capacity: determining this requires a long-range extrapolation from available data on yield and stock size. We had no strong evidence for what the upper limit of the carrying capacity might be, so instead of discarding this model right away, we changed the upper truncation point of the K-prior to an arbitrarily chosen value of 15 million tonnes. This would still interfere with the posterior for K, but increasing the upper limit further did not seem to influence the posteriors of the other parameters in the model significantly.
|
Model 1 was then able to produce a reasonable simulation of the observed data (Figure 2). The probabilities of obtaining a more extreme observation than those realized in the two data series on stock size were in the range of 0.06–0.5, i.e. the observations did not lie in the tails of their posterior distributions. For the cpue series, the data point of 1987 was suggested by a low CPO to be a relatively poor fit to the model. Minor problems in capturing the extreme catch rates of 1979 and 1982 were also indicated by these points having relatively large residuals and small CPOs. The survey series was generally less well estimated; the 1991 value had the largest residual and smallest CPO. Some of the parameters showed high linear correlations (Table 2). The catchabilities (qc and qs) were negatively correlated with the parameter for carrying capacity, K, and to some extent with MSY. MSY and K were only slightly positively correlated.
|
For the parameters m (represented by BMSYK–1), K, and P1, the posterior distributions tended to approximate the input priors (Figure 1). The posterior for MSY was positively skewed and showed a mode at 90 000 t and upper and lower quartiles at 82 700 t and 187 000 t. The catchabilities, qs and qc, showed marked peaks at 0.045 and 0.00019, respectively, but had relatively wide posterior distributions. The estimated CV of the observed cpue series had a median at about 8.4%, and for the survey series at 14.9%. The process error,
, had a median of 10.3%.
Model 2 (cod predation effect included)
Entering a predation effect in the model only had minor influence on the model diagnostics. The probabilities of obtaining a more extreme observation set than the realized ones ranged from 0.07 to 0.5. For the cpue series it was still the peaks in 1979, 1982, and 1987 that had the largest log-residuals and smallest CPOs. Similarly, the greater variability of the survey series was, as before, less well captured, with the values of 1991 being the most obvious deviation. Large correlations between parameters still prevailed (Table 3).
|
However, precision of the key parameters had improved substantially. Not only did the posterior medians of the MSY differ substantially between the two models but a decline in its relative inter-quartile range (range/median), from 86% to 22%, was also observed. K now had a more peaked posterior distribution with a distinct mode around 740 000 t (Figure 1). The catchabilities, qs and qc, were larger and also more precisely estimated than in model 1, with modes at 0.52 and 0.0022, respectively. The posterior distribution for BMSYK–1, approximately uniform in the absence of predation, changed to favour values at the lower end of its range. The prior of the initial state, P1, was still not updated.
The estimated CVs of the observed biomass indices, and , were similar to those estimated by model 1. However, the mode of the posterior for the process error, , had, in line with the observations noted above, decreased from 0.094 to 0.071. The posterior for the CV of the observed cod biomass series, , markedly updated the prior (Figure 1), with a clear mode at 0.28.
The parameter principally determining the cod predation rate, P50%, was also markedly updated with a posterior showing a mode at 3.63 (Figure 1). As expected, there was no information in the data about Vmax, and the posterior copied the prior.
Model comparison
The goodness-of-fit statistics did not provide clear advice on model selection. However, there were large differences in the precision of parameters estimates (Figures 1 and 2). The process error was reduced by more than one-fourth by the inclusion of the predation effect, and the relative inter-quartile ranges of MSY and K went down by about two-thirds.
How well the stock-production curve is defined by the data may be visualized by plotting the stock size and corresponding production calculated directly from the observed biomass indices. This was done by applying the MCMC-sampled catchabilities to the actual index values, to obtain biomass estimates. Production was calculated by subtracting biomass in the current year from biomass in the next, then adding catch and if necessary estimated predation. Although only a minor difference is seen in the estimated median stock-recruitment curve, the larger variability of the model 1 estimates was apparent (Figure 3).
|
The two models produced different simulations of likely stock biomass trajectories back in time before 1976, when the cpue series starts (Figure 2). Having only the time-series of catch available as input data for that period, model 1 predicts a more or less constant stock size around the high level of the 1976 value. Model 2, using both catch and cod data, estimates the shrimp biomass during the period of high cod abundance in the 1960s to have been at about one-third of its 1976–2002 level. This scenario agrees with the general belief that the stock in those years was at a lower level (S. A. Horsted, pers. comm.). Therefore, model 2 was chosen for use in the assessment process, and subsequent comments refer to that model.
Cod consumption
The estimated annual consumption of shrimp by cod in the years 1956–2002 was in the range 0 t to about 120 000 t (Figure 2), of the same order of magnitude as the catches taken by the fishery. Estimated consumption has declined since 1960 (Figure 2), in step with a decline in cod abundance. A short-lived resurgence of the cod stock in the late 1980s caused consumption estimates to increase. The cod disappeared in the early 1990s, and estimates of consumption dropped to 0.
Stock development
The shrimp stock has been exposed since the 1950s to two different environmental regimes: one with high and the other with low cod abundance. The model indicates that the stock dynamics have responded to the difference. The trajectory of the median estimate of "biomass ratio" (Pt = BtBMSY–1) plotted against "mortality ratio" (ZtZMSY–1; Figure 4) starts in 1956 at half the optimum biomass ratio (=1), and at a mortality ratio considerably higher than 1. The stock maintained itself in this region during the years when cod were abundant. When the cod stock declined in the late 1960s and predation pressure was lifted (Figure 2), shrimp stock biomass increased and eventually began cycling in the left upper corner of the graph (Figure 4) during the current regime of scarce cod. From the early 1970s, the estimated median biomass ratio ranged from about 0.96 to 1.70, and the probability that it had been below the optimum level was small for most years, i.e. it seemed likely that the stock had been at or above its MSY level throughout the modern fishery. A steep decline in cpue was noted in the late 1980s and early 1990s following a short-lived resurgence of the cod stock, and the median estimate of biomass ratio dipped just below the optimum in 1990 and 1991 (Figure 4). The stock has since increased.
|
Production potential
The median estimate of the maximum annual production surplus available equal to the fishery and cod (MSY) was estimated to be 100 600 t. The risk function relating to the probability of exceeding MSY to the combined removal by fishery and cod predation is given as the integral of this distribution (Figure 5).
|
Prediction and risk
The probabilities of transgressing chosen references in response to different management options (optional catch series of 80, 90, 100, 110, and 120 thousand tonnes) were derived for a 10-year period assuming that the cod stock would remain at its current level. For demonstration, a biomass reference was chosen as equal to BMSY and a mortality reference equal to ZMSY. Note that the references used here are not single points, but probability distributions of possible values.
Increasing the catches will (not surprisingly) increase the risks of overfishing and of drive the biomass to levels below its most productive, and for each catch option the estimate of theserisks increases with time (Figure 6). Even for catch levels at 80 000 t, considerably below estimated median MSY, risk estimates will increase with time alone as a result of increasing uncertainty in the predictions of stock size.
|
|
Discussion |
|---|
|
TopIntroductionModel developmentResultsDiscussionAppendixReferences |
|---|
The basic ideas embedded in applied production models are the concept of carrying capacity (K) and Maximum Sustainable Yield (MSY). Resources (food, space, etc.) available to the individual in the ecosystem are limited. Density-dependent effects will become effective as the population grows, and along with the growth potential of the stock decide the form of the dome-shaped surplus-production curve of this class of model. MSY and K are not in reality constant over time, but vary with changes in a host of environmental factors. When modelled with a time-invariant value, these parameters represent averages over time of their temporary values. The dynamics of the stock under the influence of fishing and predation are also responding to variation in these basic parameters, which is therefore subsumed in the process error.
Parameter estimation
In the Bayesian framework, fundamental absence of information in the data will yield posteriors as a copy of the input priors. For the data to carry information on all the parameters of any such model, the biomass should vary widely both above and below BMSY. If the available data do not span these conditions, problems in fitting stock-production models by any method can be expected (Hilborn and Walters, 1992).
The plot of the estimated stock-production function and the underlying data points (Figure 3) suggest that the available time-series of indexed stock biomass for the West Greenland shrimp stock does not span the range from 0 to K. It looks as though the stock size has been fluctuating above BMSY between 1976 and 2002. For these reasons in particular, the shape parameter, m, could not be well determined. However, in the applied parametrization using relative biomass instead of absolute, both the starting point and the location of the maximum of the stock-production function are fixed (maximum at B/BMSY = 1), so the suggested regression line in Figure 3 may not be as speculative as it might look at first sight.
Even though the conditions for estimation of some parameters are not optimal, it may still be possible to obtain good estimates of parameters relevant for management. Fortunately, MSY is the easiest single parameter to estimate. If the range of biomass includes BMSY, good estimates of MSY can be obtained independently of other parameters. In some cases, this can be done even if the range of biomass merely approaches BMSY from one side or the other, as seen in this case. Model 2 yields estimates of production of stock sizes close to BMSY, and of a large part of the upper range of possible stock sizes (Figure 3), and the estimated posterior distribution of MSY is reasonably tight (Figure 1).
More uncertainty was associated with the estimation of the carrying capacity, K. An informative prior for q or K itself to help in scaling absolute stock size could improve conditions, but such prior knowledge was not available. When compared with model 1 results, the inclusion of predation increased the degree to which the prior probability distributions of K were updated by the Bayesian fitting process and decreased the variance of the resulting posterior distribution (Figures 1
and 2).
The posterior for the initial biomass, P1, tended to approximate the input priors whether or not these were uninformative (Figure 1). P1 is the biomass level of 1955, i.e. 20 years before data on relative stock size were available. Therefore, signals from the data about stock biomass before 1976 would be reflected in overall trends in that series rather than just defining the initial state. Altering the prior for P1 would have some effect on the first couple of years of biomass estimates, but after that, series conditioned on different priors would merge. Therefore, the first couple of years in the estimated stock trajectory is likely to be biased and should be interpreted with caution. The estimates of other parameters were only marginally affected by altering the P1-prior, even with the use of a prior uniform between 0 and 3.
We conclude that the lack of contrast in the stock dynamics in this case is not an impediment to obtaining reasonable estimates of parameters relevant to management. However, an adaptive management programme, such as investigating stock response to higher quotas, would enable us further to explore the production potential of the stock and to make firmer statements on its response to different exploitation levels. The security of adaptive management is, however, dependent on having confidence in the basic model of stock dynamics. MSY is estimated in the absence of fishing (which is usual) and predation. Fishing is assumed to be under management control and therefore capable of being set at an appropriate level of exploitation, so given adequate information about stock status and dynamics, optimal management is achievable. Cod predation on shrimp apparently needs to be allowed for in the management of the shrimp stock. A fishable yield equal to the calculated MSY minus the estimated predation at BMSY is available, and if a larger fishable yield is available, it is available from a standing stock smaller than BMSY. Managing the stock relative to the estimated BMSY therefore appears to be conservative. However, predation is not controlled by the managers of the shrimp stock, and therefore introduces an element of unpredictability to the system. Furthermore, the cod stock itself may be subject to management with its own objectives; a joint objective for the two stocks may also be possible given adequate analysis of the social and economic, as well as the biological, factors.
Predation by cod
The biological interactions between species introduced in model 2 involve predation by cod and its resulting primary effect on the mortality of shrimp. The secondary factor of predation on the growth of the predator was in this case not included. Cod are opportunistic predators with a wide food spectrum, and the dynamics of a single prey population, e.g. shrimp, is likely cancelled out by variations in the wide selection of other prey populations, unless it is part of a major directed change in the ecosystem. Also, the cod biomass data series entered would presumably already include any effect of such a feedback process.
In this study a significant improvement in the precision of estimates of principal parameters was gained by the inclusion of the cod predation effect in model 2, presumably mainly because predation provides an explanation of features of the biomass trajectory that otherwise are hard to explain: the boom in cod predation in the late 1980s coincides with the drop in cpue, the moderate levels of cod predation in the 1970s coincide with the fairly stable biomass at that time, and low predation in the 1990s coincides with increasing biomass.
The model estimated the medians of the posterior distributions of the annual consumption by cod in the years 1956–2002 to range from 200 t to about 100 000 t. This is of the same order of magnitude as the catches taken by the fishery (Figure 2). Although it is well known that cod eat shrimp, little information could be found for independent validation of the intake quantities that this model estimated. The effect of cod predation could well be large. However, the mechanics behind the shrimp–cod interaction are not fully known, so an apparent effect of cod on the shrimp stock could also perhaps be attributed to opposite, but independent, stock responses to the same shifts in environmental conditions (Lilly et al., 2000).
In the light of the investigations by Grundwald (1998), the assumed Omax at 3 kg kg–1 y–1 looks to be of the right order of magnitude. Grundwald (1998) calculated a daily ration of 0.4–0.9% body weight per day for 1989–1992. Shrimp made up 12–37% of the diet, implying that an average kg of cod ate 0.2–1.2 kg of shrimp per year. In the Barents Sea, annual point estimates of predation rate during the period 1985–2001 were 0.1–0.4 kg-shrimp kg-cod–1 (calculation based on ICES, 2002). The estimates of P50% had a median of 3.8, meaning that shrimp densities well above the estimated carrying capacity are needed before predation rate can reach even half its maximum. This suggests that, for this system, we are only dealing with the lower limb of the sigmoidal functional response, and that within a narrower range of shrimp stock sizes (maybe for most practical cases), cod predation rate could be considered linearly related to shrimp abundance.
The increased precision of model 2 parameter estimates reduces uncertainties in projections of future developments of stock biomass (Figure 2), at least during the current regime of extreme scarcity of cod where the error contribution from the estimation of predation is small. Consumption could probably be estimated with more precision, and confidence in the precision itself could be augmented if more direct estimates of cod predation rates under various scenarios of cod and shrimp abundance were available. Therefore, should the cod stock enter a rebuilding phase, it would be desirable to initiate, perhaps at some threshold stock level, a programme of direct estimation of consumption rates to improve information in this area and thus the predictive capability of the model.
Advantage of model
The model combines the flexibility of state-space modelling and Bayesian methods, and offers a quantitative alternative to current models used in northern shrimp stock assessments. Among its advantages are (i) the state-space structure permits inclusion of both process and observation errors, and allows for time-variant equations of population dynamics and of data relations to the modelled process. Therefore, the model may easily be tailored to account for particular observations in individual years such as an atypical recruitment variation, a temporary change in discarding practice, or, as used in the current model, increased observation error for the survey and cpue series in certain years. (ii) Bayesian methods can include ancillary knowledge in a model as priors and aid parameter estimation in data-scarce situations (i.e. in most marine population dynamics modelling). This potential was not explored to its full extent in the present modelling exercise, because most priors were constructed as non-informative. However, all in all, the model, and the conclusions drawn from it, showed low sensitivity to changes in prior distributions, indicating that the data were sufficiently informative. (iii) As the full distributions of model parameters or derived quantities are estimated, it is technically straightforward to quantify the risk associated with different management options and so provide a quantitative decision framework for the fishery managers.
|
Appendix |
|---|
|
TopIntroductionModel developmentResultsDiscussionAppendixReferences |
|---|
Overview of stochastic model of shrimp stock biomass dynamics incorporating predation by Atlantic cod to be solved using Bayesian inference. The t indices time-steps in years, t = (1,...,N); "
" means "distributed as"; N(µ,2) is the normal distribution with mean µ and variance 2. Correspondingly, log N is the lognormal distribution, G the gamma distribution, and U is the uniform distribution. Other parameters are explained in the text.
Observables, data
- cpuet = qcBMSYPt exp(t)
- survt = qsBMSYPt exp(t)
- codt = codt
- Ct = Ct
Unobservables,
Biomass states
- (6)
- (7)
- (7)
- (8)
- (9)
- (10)
- (11)
- (9)
- (12) log qs U(–10, 1)
- (13) log qc
U(–15, 1)
- (13) log qc
- (14) Vmax N(3, 0.1)
- (15)
- (15)
- (16)
- (17) MSY
U(35,1000)
- (18) log(K)
U(ln(100), ln(6000))
- (17) MSY
Joint prior density, p()
- (19)
Likelihood, p(data|)
- (20)
Joint posterior probability density, p(|data)
- (21) p(|data)
p()p(data|)(Bayes theorem)
|
Acknowledgements |
|---|
We thank colleagues and reviewers, especially Andre Punt and Verena Trenkel, for careful and helpful reviews and comments. These models could not have been explored as they have been without the long data series carefully assembled and checked over the years by our colleagues at Pinngortitaleriffik.
|
Footnotes |
|---|
1 The symbol "
" is to be interpreted as "is distributed as". 
|
References |
|---|
|
TopIntroductionModel developmentResultsDiscussionAppendixReferences |
|---|
-
Anderson P.J. (2000) Pandalid shrimp as indicators of ecosystem regime shift. Journal of Northwest Atlantic Fishery Science 27:1–10.
Apollonio S., Stevenson D. K., Dunton E. E. (1986) Effects of temperature on the biology of the northern shrimp, Pandalus borealis, in the Gulf of Maine. Technical Report National Oceanic and Atmospheric Administration, NMFS 42. 22 pp.
Bayes T. (1763) An essay towards solving a problem in the doctrine of chances. Philosophical Transactions of the Royal Society 330–418 Reprinted in: Biometrika: 45 (1958), 293–315.
Berenboim B.I., Dolgov A.V., Korzhev V.A., Yaragina N.A. (2000) The impact of cod on the dynamics of Barents Sea shrimp (Pandalus borealis) as determined by multispecies models. Journal of Northwest Atlantic Fishery Science 27:69–75.
BjÖrnsson B. and Steinarsson A. (2002) The food-unlimited growth rate of Atlantic cod (Gadus morhua). Canadian Journal of Fisheries and Aquatic Sciences 59:494–502.
Brooks S. and Gelman A. (1998) General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics 7:4434–455.[CrossRef][Web of Science]
Cadrin S.X. and Clark S.H. (1999) Application of catch-survey models to the northern shrimp fishery in the Gulf of Maine. North American Journal of Fisheries Management 19:551–568.[CrossRef]
Carlsson D.M., Folmer O., Kanneworff P., Kingsley M.C.S., Pennington M. (2000) Improving the West Greenland trawl survey for Pandalus borealis. Journal of Northwest Atlantic Fishery Science 27:151–160.[CrossRef]
Congdon P. (2001) Bayesian Statistical Modelling(Wiley, London) 556 pp.
FAO. (1996) Precautionary approach to fisheries. FAO Technical Paper, 350/2. 210 pp.
Fletcher R.I. (1978) Time-dependent solutions and efficient parameters for stock–production models. Fishery Bulletin US 76:377–388.
Fournier D.A., Sibert J.R., Terceiro M. (1991) Analysis of length frequency samples with relative abundance data for the Gulf of Maine northern shrimp (Pandalus borealis) by the MULTIFAN method. Canadian Journal of Fisheries and Aquatic Sciences 48:591–598.
Gelfand A.E. and Dey D.K. (1994) Bayesian model choice: asymptotic and exact calculations. Journal of the Royal Statistical Society B56:501–514.
Gelman A., Carlin J.C., Stern H., Rubin D.B. (1995) Bayesian Data Analysis(Chapman & Hall, New York) 525 pp.
Gelman A., Meng X.L., Stern H.S. (1996) Posterior predictive assessment of model fitness via realized discrepancies. Statistica Sinica 6:733–807.[Web of Science]
Gelman A. and Rubin D.B. (1992) Inference from iterative simulation using multiple sequences. Statistical Science 7:457–511.[CrossRef]
Geweke J. (1992) Evaluating the accuracy of sampling-based approaches to calculating posterior moments. In Bernado J.M., Berger J.O., Dawid A.P., Smith A.F.M. (Eds.). Bayesian Statistics(Clarendon Press, Oxford, UK) vol. 4: pp. 169–194.
Gilks W.R., Clayton D.G., Spiegelhalter D.J., Best N.G., McNeil A.J., Sharples L.D., Kirby A.J. (1993) Modelling complexity: applications of Gibbs sampling in medicine. Journal of the Royal Statistical Society B55:39–52.
In Gilks W.R., Richardson S., Spiegelhalter D.J. (Eds.). Markov Chain Monte Carlo in Practice (1996) (Chapman & Hall, London) 512 pp.
Gilks W.R., Thomas A., Spiegelhalter D.J. (1994) A language and program for complex Bayesian modelling. Statistician 43:169–178.[CrossRef][Web of Science]
Grundwald E. M. (1998) NahrungsÖkologische Untersuchungen an Fischbeständen im Seegebiet vor WestgrÖnland. PhD dissertation, Christian-Albrechts-Universität, Kiel, Germany. 208 pp.
Heidelberger P. and Welch P. (1983) Simulation run length control in the presence of an initial transient. Operations Research 31:1109–1144.
Hilborn R. and Walters C.J. (1992) Quantitative Fisheries Stock Assessment(Chapman & Hall, New York) 570 pp.
Holling C.S. (1959) Some characteristics of simple types of predation and parasitism. Canadian Entomologist 91:385–398.
Hvingel C. (2002) The fishery for northern shrimp (Pandalus borealis) off West Greenland, 1970–2002. NAFO SCR Document, 02/151, Serial No. N4780. 19 pp.
Hvingel C. and Kingsley M.C.S. (2000) The uncertainty of an assessment procedure for the West Greenland stock of northern shrimp, Pandalus borealis (Krøyer). Journal of Northwest Atlantic Fishery Science 27:183–192.
Hvingel C., Lassen H., Parsons D.G. (2000) A biomass index for northern shrimp (Pandalus borealis) in Davis Strait based on multiplicative modelling of commercial catch-per-unit-effort data (1976–1997). Journal of Northwest Atlantic Fishery Science 26:25–36.
ICES. (1996) Report of the North Western Working Group. ICES Document, CM 1996/ACFM: 20. 311 pp.
ICES. (2001) Report of the Advisory Committee on Fishery Management, 2001. ICES Cooperative Research Report, 246. 895 pp.
ICES. (2002) Report of the Arctic Fisheries Working Group. ICES Document, CM 2002/ACFM: 18. 451 pp.
Kanneworff P. and Wieland K. (2002) Stratified-random trawl survey for northern shrimp (Pandalus borealis) in NAFO Subareas 0 + 1 in 2002. NAFO SCR Document, 02/148, Serial No. N4777. 25 pp.
Kass R.E. and Wasserman L. (1996) The selection of prior distributions by formal rules. Journal of the American Statistical Association 91:1345–1370.
Koeller P.A. (2000) The relative importance of abiotic and biotic factors to the management of northern shrimp (Pandalus borealis) fishery on the Scotian shelf. Journal of Northwest Atlantic Fishery Science 27:21–33.
Koeller P.A., Savard L., Parsons D.G., Fu C. (2000) A precautionary approach to assessment and management of shrimp stocks in the northwest Atlantic. Journal of Northwest Atlantic Fishery Science 27:235–246.
Lilly G.R., Parsons D.G., Kulka D.W. (2000) Was the increase in shrimp biomass on the Northeast Newfoundland shelf a consequence of a release in predation pressure from cod? Journal of Northwest Atlantic Fishery Science 27:45–61.
McAllister M.K. and Kirkwood G.P. (1998) Bayesian stock assessment: a review and example application using the logistic model. ICES Journal of Marine Science 55:1031–1060.
Macdonald P.D.M. and Pitcher T.J. (1979) Age-groups from size–frequency data: a versatile and efficient method of analysing distribution mixtures. Journal of the Fisheries Research Board of Canada 36:987–1011.[Web of Science]
Magnússon K. J. and Pálsson ó.K. (1991a) The predatory impact of cod on shrimps in Icelandic waters. ICES Document, CM 1991/K: 31. 17 pp.
Magnússon K.J. and Pálsson ó.K. (1991) Predator–prey interactions of cod and capelin in Icelandic waters. ICES Marine Science Symposia 193:153–170.
Meyer R. and Millar R.B. (1999) Bayesian stock assessment using a state-space implementation of the delay difference model. Canadian Journal of Fisheries and Aquatic Sciences 56:37–52.
NAFO. (2001) Scientific Council Reports 2000(NAFO, Halifax, Canada) 303 pp.
Nilssen K.T., Pedersen O-P., Folkow L., Haug T. (2000) Food consumption estimates of Barents Sea harp seals. The North Atlantic Marine Mammal Commission Scientific Publications 2:9–27.
Pedersen S.A. and Riget F. (1993) Feeding habits of redfish (Sebastes spp.) and Greenland halibut (Reinhardtius hippoglossoides) in West Greenland waters. ICES Journal of Marine Science 50:445–459.
Pella J.J. and Tomlinson P.K. (1969) A generalized stock production model. Bulletin of the Inter-American Tropical Tuna Commission 13:419–496.
Punt A. and Hilborn R. (1997) Fisheries stock assessment and decision analysis: the Bayesian approach. Reviews in Fish Biology and Fisheries 7:35–63.[CrossRef][Web of Science]
Quinn T.J., Turnbull C.T., Fu C. (1998) A length-based population model for hard-to-age invertebrate populations. In Funk F. (Ed.), et al. Fishery Stock Assessment Models pp. 531–556 Alaska Sea Grant College Program Report, AK-SG-98-01.
Raftery A.L. and Lewis S. (1992) How many iterations in the Gibbs sampler? In Bernardo J.M., Berger J.O., Dawid A.P., Smith A.F.M. (Eds.). Bayesian Statistics(Calrendon Press, Oxford, UK) vol. 4: pp. 763–774.
Savard L., McQuinn I. H., Parsons D. G. (1991) Investigating the potential of sequential population analysis for northern shrimp (Pandalus borealis) in Division 0A. NAFO SCR Document, 91/71, Serial No. N1955. 12 pp.
Schaefer M.B. (1954) Some aspects of the dynamics of populations important to the management of the commercial marine fisheries. Bulletin of the Inter-American Tropical Tuna Commission 1:27–56.
Schnute J. (1994) A general framework for developing sequential fisheries models. Canadian Journal of Fisheries and Aquatic Sciences 51:1676–1688.
Spiegelhalter D.J., Thomas A., Best N. (2000) WinBUGS version 1.3 User Manual(MRC Biostatistics Unit, Institute of Public Health, Cambridge, UK).
Stefánsson G., Skúladóttir U., Pétursson G. (1994) The use of a stock production type model in evaluating the offshore Pandalus borealis stock of north Icelandic waters, including the predation of northern shrimp by cod. ICES Document, CM 1994/K: 25. 13 pp.
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  M. C. S. Kingsley, K. Wieland, B. Bergstrom, and M. Rosing Calibration of bottom trawls for northern shrimp ICES J. Mar. Sci., September 1, 2008; 65(6): 873 - 881. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||