© 2003 by ICES/CIEM International Council for the Exploration of the Sea/Conseil International pour l'Exploration de la Mer
Modelling stochastic fish stock dynamics using Markov Chain Monte Carlo
Danish Institute for Fisheries Research, Charlottenlund Slot 2920 Charlottenlund, Denmark
*Correspondence to P. Lewy. e-mail: pl{at}dfu.min.dk.
A new age-structured stock dynamics approach including stochastic survival and recruitment processes is developed and implemented. The model is able to analyse detailed sources of information used in standard age-based fish stock assessment such as catch-at-age and effort data from commercial fleets and research surveys. The stock numbers are treated as unobserved variables subject to process errors while the catches are observed variables subject to both sampling and process errors. Results obtained for North Sea plaice using Markov Chain Monte Carlo methods indicate that the process error by far accounts for most of the variation compared to sampling error. Comparison with results from a simpler separable model indicates that the new model provides more precise estimates with fewer parameters.
Keywords: stochastic survival, stock dynamics, process error, sampling error, age-structured assessment, Markov Chain Monte Carlo estimation, North Sea plaice
Received 13 August 2002; accepted 26 March 2003.
 |
Introduction |
|---|
 Top Introduction Population dynamics modelsEstimation methodsSimulation experimentsMaterials and software usedResultsResults of simulation...DiscussionAppendix AReferences |
|---|
The precautionary approach has become a basic concept in fish stock management (Anon., 1995). The concept implies that uncertainties have to be taken into account in the assessment of the fishery e.g. by estimating the risk that the stock biomass falls below a certain critical limit. The quantification of these uncertainties has emphasised the need for developing stochastic assessment approaches.
Numerous stochastic assessment methods including frequentist, state space, time series and Bootstrap models have been suggested (e.g. Doubleday, 1976; Fournier and Archibald, 1982; Deriso et al., 1985; Gavaris, 1988; Lewy, 1988; Methot, 1990; Powers and Restrepo, 1993; Gudmundsson, 1994; Schnute, 1994; Patterson and Melvin, 1996). In the 1990s Bayesian methods have been used in connection with biomass dynamics models (Kinas, 1996; McAllister and Kirkwood, 1998; Millar and Meyer, 2000), with models that bridge the gap between biomass and age-structured models (McAllister et al., 1994; Meyer and Millar, 1999a) and with fully age-structured models (Ianelli and Fournier, 1998; Virtala et al., 1998; Patterson, 1999).
However, only a few authors have considered the stock dynamics as a stochastic process or the stock size of a cohort as an unknown stochastic variable. Virtala et al. (1998) have formulated a consistent model, where the number of survivors in a cohort, the fish caught and the number dead from natural causes are assumed to be multinomially distributed. As noted by the authors the limitation of this model is that for stocks with millions of fish the model in practice becomes nearly deterministic. Sullivan (1992) and Gudmundsson (1994) applied Kalman filter approaches to length and age-structured state space models, respectively, and thus considered the survival in a cohort as a stochastic process. Schnute and Richards (1995) formulated an age-structured state space model including both process and measurement errors. The properties of estimates based on a simpler model that only include measurement error were evaluated using data generated from the general model.
In this paper an age-structured assessment model with structural relations between variables and parameters is developed, where stock numbers are treated as unknown stochastic variables subject to process error and the catch variables subject to both sampling and other process errors. Estimates of parameters including process variances and predicted stock numbers have been obtained using likelihood-based Markov Chain Monte Carlo (MCMC). The assessment model enables the inclusion of detailed sources of information used in standard age-structured assessment such as catch-at-age and effort data from commercial fleets and research surveys. Catch data without effort information is combined into one fleet and all catch data by fleet are treated as stochastic variables subject to sampling and process errors. The usual problem of weighting the different sources of information (catch-at-age observations by fleet, the survival and recruitment processes) is solved by estimating the associated variances. Using data for North Sea plaice the importance of the process error is investigated and the properties of the estimated biomass and mortality rates are compared with the results of a simple separable model as well as with results from extended survivors analysis (XSA, Shepherd, 1999).
|
Population dynamics models |
|---|
|
TopIntroductionPopulation dynamics modelsEstimation methodsSimulation experimentsMaterials and software usedResultsResults of simulation...DiscussionAppendix AReferences |
|---|
The stochastic age-based stock assessment models considered include the type of data used in many fish stock assessments in the North Atlantic, where the following information is available for a range of years:
- Catch-at-age in numbers and effort data for commercial fleets (Cf,a,y and ef,y)
- Catch-at-age in numbers without effort data for the remaining part of total international catches (residual catches, Cres,a,y)
- CPUE by age for research surveys (Is,a,y)
All catch-at-age observations are assumed to be stochastic variables. The model applies residual catches, Cres,a,y, as well as commercial catches and effort data by fleet as observations instead of total international catches used in deterministic VPA approaches or partly stochastic approaches, which assume that total catches in numbers are known without error. The application of residual catches ensures that the catch observations for different fleets are independent variables. Effort data were incorrectly treated as covariates assumed to be known without errors.
The numbers of survivors in a cohort are considered as unobserved variables, which are subject to stochastic variations caused by fishing and natural mortality processes. Similarly, the observed numbers of fish caught are subject to stochastic variation due to the fishing process as well. Further, the catches also are subject to sampling errors. The application of a survivor and catch dynamic model that includes process variation is natural, because even with perfect knowledge of the state of the system we would not be able to accurately predict tomorrow's survivors or catch. The following lognormal distributions have been used to describe the stochastic models used for describing the dynamics of stock and catch:
|
| (1) |
|
| (2) |
|
| (3) |
|
| (4) |
survival, res, f and s the standard deviations for the survival and fishing processes, q the catchability, e the effort, T the day of year when the survey takes place and the 
s the standardised normal distribution.
Theoretically, Equation (1) does not prevent stock numbers at a given time from exceeding stock numbers at an earlier time. To ensure this condition an alternative stock dynamics model was formulated:
|
| (1a) |
Equations (1)–(4) imply that the median of the four types of variables equals the standard deterministic stock and catch expressions 
and 
. It is assumed that the observed catches conditioned on stock numbers are mutually independent and that catchability remains unchanged over time for both commercial fleets and surveys. To simplify the model it has been assumed that effort is known without error and that the oldest age group is assumed to be a "true" age group and not a plus-group.
To reduce the numbers of parameters, fishing mortality for the residual fleet defined in Equation (2) is assumed to be multiplicative:
|
|
Recruitment to the stock is also considered a stochastic process and a Ricker stock-recruitment model has been used to model the relation between recruitment and spawning stock biomass, SSB. Recruitment is assumed to be lognormally distributed:
|
| (5) |
recruit the standard deviation of the recruitment process and where the parameters, 
, ß and recruit will be estimated simultaneously with other parameters.
The number of recruits in the first year is also assumed to be lognormally distributed:
|
| (6) |
initial are the standard deviation and variance of the stock number of the first age-class in the first year, respectively.
Finally the remaining initial stock size-at-age in the first year is modelled assuming that the expected stock size is in equilibrium in the first year:
|
| (7) |
The variance parameters, recruit, and survival, are allowed to differ because they are related to different processes of stock dynamics while initial2 is different from the two process variances because it relates to the assumption of equilibrium of the initial stock.
The 55 parameters considered are:
|
|
In Appendix A the models specified in Equations (1)–(4) are shown to be partly supported by stochastic extension of the standard deterministic stock and catch differential equations. The lognormal distribution of the survivors used, Equation (1), is shown to follow directly from the stochastic formulation in Appendix A. The variance of the survival process, survial2, is shown to be the sum of the variances associated with the mortalities due to fishing and natural causes defined in Appendix A. The corresponding probability distributions of the catch observations cannot be analytically derived in support of the lognormal catch distributions assumed in Equations (2)–(4). Nevertheless it can be shown that the expected value of the catches equals the standard deterministic expressions, which approximately equals the assumption made in Equations (2)–(4).
The stock numbers defined by Equation (1), are considered as unobserved random variables predicted from the estimates of the parameters involved. This is in contrast to most stock assessment models where the initial stock sizes are treated as parameters and estimated directly.
The catch and survey observations, Cres,a,y, Cf,a,y and Is,a,y are implicitly sampling estimates used as estimates of the true values. Using the formula for conditional variances, 
and assuming that log sampling estimates are unbiased, the variances for each of the three types of log catch observations are shown to be the sum of two components: a fishing process contribution and a sampling contribution, i.e.:
|
| (8) |
|
| (9) |
|
| (10) |
and 
etc. It should be noted that if the sampling error is known externally to the model then the process variances can be estimated by the residuals, 
and 
("
" indicates estimated values). For the case-study of North Sea plaice the sampling variance was externally known. Here the importance of the process variances, process,res2, process,f2 and process,s2, is quantified by the proportions 
.
To enable comparisons of the output from our model with other approaches, we considered a simple conventional method, which does not include the stochastic survival of the fish in the sea. The stochastic catch models in this method are exactly the same as described above, except that the initial stock numbers in the sea are treated as parameters and the process models, Equations (1) and (5)–(7), are replaced by the standard deterministic model:
|
|
The initial stock sizes, that is the recruits, 
, and stock-size-at-age in the first year, 
, are selected as parameters. The remaining stock numbers by age and year are treated as functions of the initial stock size parameters and of the total mortality. The deterministic relationship between stock size, initial stock size and total mortality is:
|
|
|
|
|
Estimation methods |
|---|
|
TopIntroductionPopulation dynamics modelsEstimation methodsSimulation experimentsMaterials and software usedResultsResults of simulation...DiscussionAppendix AReferences |
|---|
For the stochastic survival model it was only possible to formulate the likelihood for given values of the stock numbers. For the unconditional likelihood, however, it was not possible to analytically derive a closed form expression for the likelihood function and thereby obtain the ML estimates. Thus the parameters have been estimated using MCMC (Gilks et al., 1996) to simulate the (unconditional) likelihood function.
For complex models with structural relationships between variables and parameters, such as the stochastic survival model considered, the so-called single component Metropolis–Hastings or Gibbs sampling (Metropolis et al., 1953; Hastings, 1970; Gilks, 1996) is an MCMC method especially suitable for simulating the likelihood function. For each of the parameters the method sequentially simulates a chain of values given the remaining parameters and variables, i.e. a Markov chain. This means that for each step the combined set of parameter values has the simultaneous distribution, corresponding to the normalised likelihood function or the posterior distribution of the parameters, as a stationary distribution. For each parameter the mean of these simulated values is used as an estimate.
For the deterministic model the calculation of the likelihood function is possible and straightforward. However, for this model the MCMC was also used to estimate parameters since it is preferable to use the same type of estimator when comparing different models.
The quantiles of the simulated distribution will be used as confidence intervals for the estimated parameters. For the deterministic survival model this is supported by simulations showing that these quantiles are in fact reasonable estimates for the confidence limits and are actually better estimates than the usual inverse Hessian matrix used for the ML estimates (Nielsen and Lewy, 2001). The method is a strong tool for simulating the likelihood especially in cases where it is not possible to analytically derive the likelihood function. The difference between the MLE and this estimator lies in the MLE being the maximum of the (normalised) likelihood function while the new estimator being the mean.
When implementing the parameter estimation method it has been necessary to restrict parameter space to finite intervals. The limits have been chosen such that the intervals are sufficiently wide not to affect the sampled parameter distribution. For the stochastic survival model the restricted intervals are as follows: Fres,y, y > 1, Fres,a and qf,a should lie in the interval (0–10), qs,a, res2, f2, s2, initial2, recruit2 and survival2 in (0–2), ln(µstart) and in (0–100) and ß in (0–1).
Due to numerical calculations, effort data have been normalised around the mean over time to avoid excessively small values for the catchability parameters (the parameters were found to lie between 0.01 and 1). For the deterministic model the chosen interval for initial stock size, NminA,y and Na,1, a>minA, is (0–1010).
The MCMC estimates obtained may also be interpreted as Bayesian mean posterior estimates considering the uniform distributions over the restricted intervals as prior distributions of the parameters. The confidence intervals may correspondingly be treated as credibility intervals.
A basic problem with MCMC is that one has to determine a simulation step for which the simulated chain is effectively at equilibrium. The values simulated before this step, the so-called "burn-in" period, then have to be discarded. The convergence of the chains was examined partly by visual inspection and partly by formal procedures.
For the stochastic model the Gelman–Rubin convergence diagnostic (Gelman and Rubin, 1992) was computed in order to determine the "burn-in" steps. The idea of the diagnostic is to simulate a number of chains started on over dispersed values and see when they become indistinguishable. Formally, the convergence is evaluated by considering the ratio of the between and within sequence variances of the different chains. This ratio converges to one when the number of steps tends to infinity. Gelman and Rubin suggest that the burn-in period ends when the diagnostic is less than 1.1 or 1.2.
|
Simulation experiments |
|---|
|
TopIntroductionPopulation dynamics modelsEstimation methodsSimulation experimentsMaterials and software usedResultsResults of simulation...DiscussionAppendix AReferences |
|---|
Simulation experiments were carried out to investigate the properties of the estimators. This was done by simulating the catch observations based on the models described and using fixed known parameters obtained from parameter estimation. The catch observations were simulated in the following way:
- The model used was the same as described by Equations (1)–(7).
- The parameters,
, used were the values estimated applying data described in the next section.
- Fres,a,y and Ff,a,y were calculated, the latter using effort data and catchabilities, qf,a.
- NminA,1 was predicted by randomly drawing from the lognormal distribution, (Equation (6)).
- Na,1, a=2,...,A were predicted by randomly drawing from the lognormal distribution, (Equation (7)).
- SSB1was calculated.
- For y=2 recruitment NminA,y was randomly drawn from Equation (5).
- For a=2,...,A Na,y was randomly drawn from Equation (1a) and SSBy calculated.
- Steps 7 and 8 were repeated as long as y<Y.
- The catch observations, Cres,a,y, Cf,a,y and Is,a,y, were generated from the lognormal distributions (Equations (2)–(4)).
One hundred replications with 40 000 chains were generated and for each set of observations the parameters were estimated and stock numbers and SSB were predicted. The empirical mean and variance were calculated for the parameters and compared to the true values by calculating the relative bias, [(estimated–true)/true]x100. The relative bias was also calculated for predicted SSB and recruitment. As these quantities differ for each replication the relative bias was calculated separately for each. Finally, the mean and variance of the relative bias were calculated.
|
Materials and software used |
|---|
|
TopIntroductionPopulation dynamics modelsEstimation methodsSimulation experimentsMaterials and software usedResultsResults of simulation...DiscussionAppendix AReferences |
|---|
The models have been applied to a set of data for North Sea plaice for the period 1988–1997. These data consist of:
- Catch-at-age and effort data for the Dutch (ages 2–9, years 1989–1997) and English (ages 4–10, years 1988–1997) commercial beam trawl fleets.
- Catch-at-age data for the combined fleet without effort data (the residual fleet) (age groups 1–12, years 1988–1997).
- Survey indices for the Dutch beam trawl (ages 1–7, years 1988–1997) and the Sole Net Survey (ages 1–3, years 1988–1997)
These data and the mean weight-at-age used are the same as used by the ICES Assessment Working Group (ICES, 1999) except that the working group includes age groups 1–14 in their XSA analysis. The oldest age group considered here, age group 12, is treated as a real age group and not as a plus-group. However, as the catches of age groups 13 and 14 were very small the results should not be affected significantly. The XSA results presented in this paper are working group estimates (ICES, 1999).
Estimates of the sampling error for plaice in the North Sea were provided by the EU project, EMAS (EMAS, 2001). Average CVs by age group for the Dutch, English and Danish fleets in the period 1991–1998 were calculated by bootstrapping the samples (Tables 3.16, 3.28 and 3.45 of EMAS, 2001) and were used as estimates of sampling,Netherland2, sampling,England2 and sampling,res2, respectively. The Danish figures were used as estimates of sampling,res2, because the Danish catch constitutes the largest national proportion of the residual catch (about 40% in 1997).
The software package, WinBUGS 1.4 (Bayesian inference using Gibbs sampling, Spiegelhalter et al., 2000) was used to simulate the posterior distributions of the parameters. The package is designed to sample from complex models with structural relationships between variables and parameters. The package applies a single component Metropolis–Hastings algorithm to simulate the full conditional distributions (Gilks et al., 1996). An older version of the program has been described in detail for a biomass dynamics model (Meyer and Millar, 1999b).
|
Results |
|---|
|
TopIntroductionPopulation dynamics modelsEstimation methodsSimulation experimentsMaterials and software usedResultsResults of simulation...DiscussionAppendix AReferences |
|---|
For each of the two stochastic models (1) and (1a) two independent chains with 50 000 steps were generated each with over-dispersed starting points. For each of the years 1988–1997, the Gelman–Rubin diagnostics for the spawning biomass chains for the first 5000 steps were found to lie between 1.03 and 1.2. Based on this and on visual inspection of both plots of parameters versus step number and of parameter distribution, each of these initial 5000 steps was discarded as "burn-in" period and only the remaining 90 000 steps used in the results. Comparisons of the two models showed that the estimated parameters were almost identical indicating that the truncation made in model (1a) had no practical influence on the results.
For the deterministic model 11 000 steps were generated. Based on visual inspections of the same type of plots, the first 1000 steps were discarded as a "burn in" period. The results were also compared with other runs with up to 50 000 steps, which gave almost identical results.
The fit of the model was examined by inspection of standardised residual plots (Figure 1). No systematic deviations from the model were found. However, the plots should be interpreted with caution, as log catches were only assumed normal for given stock numbers implying that residuals depending on the same stock number were dependent variables. The assumption of lognormal distribution was investigated by comparing the histogram of residuals with the density of the fitted normal distribution (Figure 2A) and by a Q–Q plot (Figure 2B). No serious indications of deviations from the normal distribution were found. Even though the Q–Q plot showed that residuals were a bit too "heavy tailed", this may, however, be due to residual dependencies.
|
|
For all but one year the stochastic survival method provides lower SSB and higher
estimates than the deterministic method and XSA, which are quite similar (Figure 3). Although it is not possible to test whether these differences are significant, the SSB1997 from the stochastic approach seems to be a more conservative estimate, as it appears to be about 19% lower than that from the deterministic approach. Correspondingly, for 1997 the 
from the stochastic approach is about 15% higher than the deterministic one. The recruitment estimates indicate that there are substantial differences for the last two years, for which the XSA estimates are larger than the estimates for the two other methods.
|
The importance of the process variances,
process,Netherland2, process,England2 and process,res2 illustrated by the proportions 
etc., indicates that in the case of North Sea plaice process error is by far the most dominant (Table 1). These results might also imply that the model is inconsistent for Danish catches of age groups 1, 11 and 12, for which the sampling variance is greater than total variance estimated.
|
The CVs connected to the survival and recruitment processes, initial stock size and catch observations by fleet are given in Table 2. A comparison of the CVs of catch observations from the two models, for the two most important fleets (the Dutch and the residual), indicates that CVs for the stochastic model are lower than for the deterministic even though it is not possible to test the significance of the difference. Furthermore, the stochastic model apparently also provides more precise estimates of SSB and fishing mortality than the deterministic one, for all years (Figure 4). As the stochastic model has fewer parameters (55) than the deterministic model (70) one can conclude that the fit of the stochastic model is at least as good as that of the deterministic and apparently provides more precise estimates with fewer parameters.
|
|
Figure 4 also shows that CVs are larger for the first and last years, especially for the average F, but also for SSB. For the first year this may mainly be due to catch and effort data for the Dutch beam trawl fleet not being available. For the last year(s) the large CVs may be due to few age groups being included in the cohorts in question.
The survival process has a rather low CV of 0.12 compared to a recruitment CV of 0.35 (Table 2) indicating that the survival process model provides a better fit to generated stock numbers compared to the ability of the Ricker stock-recruitment process to model recruitment. This is not surprising, as no well-defined stock-recruitment relationship exist for North Sea plaice.
Our results also indicate that the CVs associated with the Dutch beam trawlers and the residual fleet (accounting for the main parts of the catch) were lower than the CVs of the English beam trawlers and especially lower than the two surveys (Table 2). Hence, the exclusion of these surveys from the analysis would probably only result in insignificant changes.
|
Results of simulation experiments |
|---|
|
TopIntroductionPopulation dynamics modelsEstimation methodsSimulation experimentsMaterials and software usedResultsResults of simulation...DiscussionAppendix AReferences |
|---|
The MCMC estimator of the SSB slightly underestimates the true SSB by about 8% while the estimators of
overestimate the true values by about 15% (Table 3). Recruitment estimates and standard deviations of log catches for all fleets, res, f and s, are almost unbiased while the survival process and initial standard deviations are slightly overestimated (16–18%). Stock-recruitment parameters, however, are significantly overestimated (: 339%, ß: 101% and recruit: 63%). The correlation between and ß is 0.91. For all parameters it can be shown that the true values lie within the 95% intervals of the empirical distributions.
|
|
Discussion |
|---|
|
TopIntroductionPopulation dynamics modelsEstimation methodsSimulation experimentsMaterials and software usedResultsResults of simulation...DiscussionAppendix AReferences |
|---|
Our model of fish stock assessment which includes stochastic survival and recruitment, seemed to be more flexible than the deterministic separable model as the fit of the former is at least as good as that of the latter and appears to provide more precise estimates of biomass and mortality rates with fewer parameters. Furthermore, the biomass estimates of the stochastic model seem to be more conservative compared to those of the deterministic model, which is in agreement with retrospective analyses of the spawning stock biomass (ICES, 1999), indicating that the SSB in the last year was overestimated.
The stochastic survival model estimates of SSB and recruitment were found to be lower than those of the deterministic model, and the opposite was the case for average F. Theoretically, this could be explained in that the stochastic model—contrary to deterministic model—includes a stock-recruitment model and that the average stock numbers in the first year are assumed to be in equilibrium. To investigate this, a third model was considered, which corresponds to the deterministic model but including the same stock-recruitment relationship and stock restrictions in the first year as the stochastic model. Again, the parameters for third model were estimated using WinBUGS and MCMC. Both recruitment estimates and stock size in the first year were found to be close to the estimates of the deterministic model. This implied that the differences between the results from the stochastic and deterministic survival models were caused by the survival model and not by initial stock constrains.
Simulation experiments showed that apart from the stock-recruitment parameters, the MCMC estimates of the remaining parameters—including spawning biomass and average fishing mortality—were only slightly biased. The recruitment parameters, , ß and recruit, were all significantly biased and the estimates of and ß were highly correlated. However, despite being unable to obtain reliable recruitment parameter estimates the resulting recruitment estimates obtained using these parameters and the stock-recruitment function were almost unbiased (with a large variation).
Errors associated with the catch-at-age by fleet used in stock assessment consist of sampling error and other errors denoted as process errors (including both fishing process error and model error). For the stochastic model, results for North Sea plaice indicate that process error is by far the most important factor while sampling error only plays a minor role. For assessment purposes this implies that even if information were available on the sampling errors of the catches for different fleet components, these should not be used directly to account for different variances associated with the various sources of information. The total variance of the catch by fleet still has to be estimated in the model. However, in cases where the sampling uncertainty is high relative to the process error, the total model variance may be estimated by the sampling variance. It should be noted that the concept of process error always refers to a specified model used.
The survival process of fish in a cohort has been assumed to be lognormally distributed. It has been shown that this assumption can be derived from a stochastic reformulation of the standard deterministic differential equations that describe the stock dynamics. For the formulation chosen the log variance of the survival process has been shown simply to be the sum of the process errors due to fishing and natural mortality.
The MCMC methodology, in particular the single component Metropolis–Hastings and graphical models, has proven to be a powerful tool for making inference in complex fish stock assessment models including structural relationships between variables and parameters. It is easy to implement such complex model in the WinBUGS program.
In both the stochastic and deterministic models the catchability is assumed to be constant over years. A random walk extension of the deterministic model was formulated and implemented for which catchability was allowed to vary over years: qf,a,y|qf,a,y–1
LN(ln(qf,a,y–1),2) where 2 arbitrarily was fixed at 0.0252. This model was used for the catchability of the Dutch and English commercial beam trawl fleets. Markov chains with length of 60 000 were generated for the parameters and appeared unstable, which could indicate an over-parameterisation of the model. Therefore, if residuals for a fleet indicate a trend over time, as is the case for English beam trawlers, it is probably better to change the assumption of constant catchability to a parametric catchability model allowing for temporal trends.
|
Appendix A |
|---|
|
TopIntroductionPopulation dynamics modelsEstimation methodsSimulation experimentsMaterials and software usedResultsResults of simulation...DiscussionAppendix AReferences |
|---|
Fish stock dynamics are usually based on the standard differential equations:
, 
and 
, where t is the time, Nt is the number of fish in the sea, 
is the accumulated number of fish caught to time t, 
is the cumulative number of fish dead from natural causes at time t, Ft is fishing and Mt natural mortality. These deterministic models can be extended to be stochastic differential equations:
|
| (A1) |
|
| (A2) |
|
| (A3) |
fishing and natural are standard deviations associated with the two processes. For a given time, t0, Wt0fishing and Wt0natural are normal distributions, N(0,1). The two death processes are assumed to be mutually independent.
The models chosen imply that 
, 
and dNt in a short time interval dt, is normally distributed with the mean equal to the deterministic values and the variances equal to fishing2Nt2dt, natural2Nt2dt and survival2Nt2dt, respectively, where survival2=fishing2+natural2.
Using Ito's formula (see for example Gard, 1988) one can show that the following equation for given initial stock, N0, is the solution to the stochastic differential Equation (A3):
|
|
|
|
Assuming that Zt is constant over a time step, for instance a year, one can show that this corresponds to the number of survivors at the end of the year given the number at the beginning of the year following a lognormal distribution:
|
| (A4) |
survival2.
It is not possible correspondingly to derive the probability distributions for 
and 
(or the observed catches during a period, 
). However, it can be shown that the expected value of the catches during a period equals the standard value of 
. The coefficient of variation of the catches can also be derived analytically as a rather complicated function of F, M and fishing and natural.
It should be noted that only the variance associated with the survival, 
has been estimated in the model. In principle the variances, fishing2and natural2, could be estimated separately as well as by using the complicated relationship mentioned and the equation, 
.
When estimating parameters using Gibbs sampling it is assumed that catch observations conditioned the survivors and the mortality parameters are independent. It is not evident if this assumption is fulfilled for the model based on Equations (A1)–(A3).
The WinBUGS code used in this study is available from the authors.
|
Acknowledgements |
|---|
The present work is part of a project funded by the Danish Ministry of Food, Agriculture and Fisheries entitled "Development of improved models of fisheries impact on marine fish stocks and ecosystem".
We further thank Lasse Heje Pedersen, Stern School of Business, New York University, Søren Feodor Nielsen, the Department of Statistics and Operations Research, University of Copenhagen, for giving inspiration to develop the stochastic survival model, and Uffe Høgsbro Thygesen, DIFRES, for valuable discussions. Finally, we are grateful to two anonymous referees whose comments and suggestions considerably improved the manuscript.
|
References |
|---|
|
TopIntroductionPopulation dynamics modelsEstimation methodsSimulation experimentsMaterials and software usedResultsResults of simulation...DiscussionAppendix AReferences |
|---|
-
Anon. (1995) Guidelines on the precautionary approach to capture fisheries and species introduction. FAO Fisheries Technical Paper No. 350 (1): 52 pp.
Deriso R.B., Quinn P.R. Jr., Neal T.J. (1985) Catch-at-age analysis with auxiliary information. Canadian Journal of Fisheries and Aquatic Sciences 42:815–824.
Doubleday W.G. (1976) A least square approach to analysing catch at age data. International North Pacific Fisheries Commission Bulletin 12:69–81.
EMAS. (2001) Evaluation of market sampling strategies for a number of commercially exploited stocks in the North Sea and development of procedures for consistent data storage and retrieval. Final Report to the EU Commission, CFP Study Project 98/075. European Commission, Directorate-General XIV Fisheries, Rue de la Loi 200, 1049 Brussels, Belgium.
Fournier D. and Archibald C.P. (1982) A general theory for analyzing catch at age data. Canadian Journal of Fisheries and Aquatic Sciences 39:1195–1207.
Gard T.C. (1988) Introduction to Stochastic Differential Equations. (Marcel Dekker, New York, NY) 234 pp.
Gavaris S. (1988) An adaptive framework for estimation of population size. Canadian Atlantic Fisheries Scientific Advisory Commission, 88/29. 12 pp.
Gelman A. and Rubin D.B. (1992) Inference from iterative simulation using multiple sequences (with discussion). Statistical Science 7:457–511.
Gilks W.R. (1996) Full conditional distributions. In Gilks W.R., Richardson S., Spiegelhalter D.J. (Eds.). Markov Chain Monte Carlo in Practise(Chapman & Hall, London, UK) pp. 75–86.
Gilks W.R., Richardson S., Spiegelhalter D.J. (1996) Markov chain Monte Carlo in practise(Chapman & Hall, London, UK) 486 pp.
Gudmundsson G. (1994) Time series analysis of catch-at-age observations. Applied Statistics 43:1117–1126.
Hastings W.K. (1970) Monte Carlo sampling methods using Markov chains and their applications. Biometrica 57:97–109.[CrossRef]
Ianelli J. N. and Fournier D. A. (1998) Alternative age-structured analyses of the NRC simulated stock assessment data. Analyses of simulated data sets in support of the NRC study on stock assessment methods pp. 81–96 NOAA technical Memorandum, NMFS-F/SPO-30, U.S. Departments of Commerce.
ICES. (1999) Report of the working group on the assessment of demersal stocks in the North Sea and Skagerrak. ICES CM 1999/ACFM: 8. 676 pp.
Kinas P.G. (1996) Bayesian fishery stock assessment and decision making using adaptive importance sampling. Canadian Journal of Fisheries and Aquatic Sciences 53:414–423.
Lewy P. (1988) Integrated stochastic virtual population analysis: estimates and their precision of fishing mortalities and stock sizes for the North Sea whiting stock. Journal du Conseil International pour l'Exploration de la Mer 44:217–228.
McAllister M.K. and Kirkwood G.P. (1998) Using Bayesian decision analysis to help achieve a precautionary approach for managing developing fisheries. Canadian Journal of Fisheries and Aquatic Sciences 55:2642–2661.
McAllister M.K., Pikitch E.K., Punt A.E., Hilborn R. (1994) A Bayesian approach to stock assessment and harvest decisions using the sampling/importance resampling algorithm. Canadian Journal of Fisheries and Aquatic Sciences 51:2673–2687.
Methot R.D. (1990) Synthesis model: an adaptable framework for analysis of diverse stock assessment data. International North Pacific Fisheries Commission Bulletin 50:259–277.
Metropolis N., Rosenbluth A.V., Rosenbluth M.N., Teller A.H., Teller E. (1953) Equations of state calculations by fast computing machine. Journal of Chemical Physics 21:1087–1091.[CrossRef]
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.
Meyer R. and Millar R.B. (1999) BUGS in Bayesian stock assessments. Canadian Journal of Fisheries and Aquatic Sciences 56:1078–1086.
Millar R.B. and Meyer R. (2000) Non-linear state space modelling of fisheries biomass dynamics by using Metropolis–Hastings within-Gibbs sampling. Applied Statistics 49:Part 3327–342.
Nielsen A. and Lewy P. (2001) Comparison of the frequentist properties of the Bayes and maximum likelihood estimators in a fish stock assessment model. Canadian Journal of Fisheries and Aquatic Sciences 59:136–143.
Patterson K.R. (1999) Evaluating uncertainty in harvest control law catches using Bayesian Markov chain Monte Carlo virtual population analysis with adaptive rejection sampling and including structural uncertainty. Canadian Journal of Fisheries and Aquatic Sciences 56:208–221.
Patterson K. and Melvin G.D. (1996) Integrated catch age analysis, version1.2. Scottish Fisheries Research Report 58:1–60.
Powers J.E. and Restrepo V.R. (1993) Evaluation of stock assessment research for Gulf of Mexico King Mackerel: benefits and costs to management. North American Journal of Fisheries Management 13:15–26.[CrossRef]
Schnute J.T. (1994) A General framework for developing sequential fisheries models. Canadian Journal of Fisheries and Aquatic Sciences 51:1676–1688.
Schnute J.T. and Richards L.J. (1995) The influence of error on population estimates from catch-at-age models. Canadian Journal of Fisheries and Aquatic Sciences 52:2063–2077.
Shepherd J.G. (1999) Extended survivors analysis: an improved method for the analysis of catch-at-age data and abundance indices. ICES Journal of Marine Science 56:584–591.
Spiegelhalter D.J., Thomas A., Best N.G. (2000) WinBUGS Version 1.3 User Manual(Medical Research Council Biostatistics Unit, Institute of Public Health, Cambridge, UK) http://www.mrc-bsu.cam.ac.uk/bugs.
Sullivan P.J. (1992) A Kalman filter approach to catch-at length analysis. Biometrics 48:237–257.[CrossRef][Web of Science]
Virtala M., Kuikka S., Arjas E. (1998) Stochastic virtual population analysis. ICES Journal of Marine Science 55:892–904.
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  G. Aarts and J. J. Poos Comprehensive discard reconstruction and abundance estimation using flexible selectivity functions ICES J. Mar. Sci., May 1, 2009; 66(4): 763 - 771. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
J. A. Bogaards, S. B. M. Kraak, and A. D. Rijnsdorp Bayesian survey-based assessment of North Sea plaice (Pleuronectes platessa): extracting integrated signals from multiple surveys ICES J. Mar. Sci., May 1, 2009; 66(4): 665 - 679. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||