BUGS in Bayesian Stock Assessments

This paper illustrates the ease with which Bayesian nonlinear state-space models can now be used for practical sheries stock assessment. Sampling from the joint posterior density is accomplished using Gibbs sampling via BUGS, a freely available software package. By taking advantage of the model representation as a directed acyclic graph, BUGS automates the hitherto tedious calculation of the full conditional posterior distributions. Moreover, the output from BUGS can be read directly into the software CODA for convergence diagnostics and statistical summary. We illustrate the BUGS implementation of a nonlinear non-normal state-space model using a Schaefer surplus production model as a basic example. This approach extends to other assessment methodologies, including delaydi erence and age-structured models. Introduction State-space models are among the most powerful tools for dynamic modeling and forecasting (Fahrmeir and Tutz 1994). They have started to enjoy an increasing popularity in sheries stock assessment (Sullivan 1992, Pella 1993, Gudmundsson 1994, Schnute 1994, Freeman and Kirkwood 1995, Reed and Simons 1996, Kinas 1996, Meyer and Millar , Millar and Meyer 1998a) as they can realistically account for both measurement and process error. However, unrealistic assumptions such as linearity of state transitions and Gaussian error distributions that are imperative for maximum likeli2 hood estimation via Kalman ltering have limited the number of sheries models that can be t within the classical/frequentist paradigm. Meyer and Millar (1998) explained how delay di erence and surplus production models can be cast into the framework of state-space modeling. They demonstrate a fully Bayesian approach using Gibbs sampling for posterior computation following Carlin et al. (1992). In contrast to the classical approach, the Bayesian approach can easily handle realistic distributional assumptions as well as nonlinearities in state and observation equations. The papers by Millar and Meyer (1998a) and Meyer and Millar (1998) provide the underlying theory of tting Bayesian state-space surplus production and delay di erence models, respectively. Here we report on signi cant progress made in facilitating the routine implementation which may have a revolutionary e ect on Bayesian stock assessment in everyday practice. This is achieved through BUGS (Bayesian Inference Using Gibbs Sampling), a recently developed software package (Spiegelhalter et al. 1996) by the Medical Research Council Biostatistics Unit, Institute of Public Health, Cambridge, England. BUGS samples from the joint posterior distribution by using the Gibbs sampler (Gilks et al. 1996), i.e. by cyclically sampling from each of the full conditionals. For reviews on BUGS the reader is referred to Thomas et al. 1992, Gilks et al. 1994, and Gentleman 1997. BUGS is available free of charge from http://www.mrc-bsu.cam.ac.uk/bugs/Welcome.html for the operating systems UNIX, LINUX, and WINDOWS, among others. It comes with a complete documentation as well as two example volumes. These examples 3 demonstrate the variety of complex models such as random e ects, generalized linear, proportional hazards, latent variable, and frailty models, amenable to a Bayesian analysis via BUGS. We will show that nonlinear non-Gaussian state-space models can be added to this list. In nonlinear non-Gaussian state-space models the full conditional distributions required for Gibbs sampling are typically not logconcave. Both Meyer and Millar (1998) and Millar and Meyer (1998a) implemented the Gibbs sampler in C-code, using the C-subroutines ARS (Gilks and Wild 1992) and ARMS (Gilks et al. 1995) to sample from univariate logconcave and non-logconcave full conditional distributions, respectively. This, however, required the explicit derivation \by hand" of the full conditional distribution of each parameter in the model a nontrivial, substantial, and tedious task which may deter the practical stock assessment scientist from tting these models. We explain how BUGS alleviates this chore by making use of conditional independence assumptions in the model which are graphically represented by a directed acyclic graph (DAG). Moreover, the new version 0.6 of BUGS contains a method to sample from any not necessarily log-concave full conditional density. This makes it now possible to t nonlinear non-Gaussian state-space models. The paper is organized as follows: In the rst section, we brie y describe a biomass dynamics model, the Schaefer surplus production model, that is the most commonly used non-age-structured model in practical stock assessments. This simple, parsimonious model was chosen for illustrative purposes. We use the same catch-e ort dataset on South Atlantic albacore tuna as in Millar and Meyer (1998a) and in Polacheck et 4 al. (1993). In section 2 we will put this surplus production model into the context of nonlinear, non-Gaussian state-space methodology. This model is then represented as a DAG in section 3. In the subsequent sections we explain how this DAG assists the model's implementation in BUGS and how to perform convergence diagnostics using CODA (Best et al. 1995). 1 The Problem The data available for stock assessment purposes quite often consist of a time series of annual catches Ct; t = 1; : : : ; N , and relative abundance indices It; t = 1; : : : ; N , such as research survey catch rates or catch-per-unit-e ort (CPUE) indices from commercial sheries. For example, Table 1 gives an historical dataset of catch-e ort data of South Atlantic albacore tuna (Thunnus alalunga) from 1967 to 1989. Age-composition data Table 1 near here are not available for this stock. This dataset has previously been analysed by Polacheck et al. (1993) and Yeh et al. (1991). Objectives include the estimation of the size of the stock at the end of 1989, and management parameters such as the maximum surplus production (MSP) and the optimal e ort, the level of commercial shing e ort required to harvest MSP. When only catch-e ort data are available, biomass dynamics models are the primary assessment tools for many sheries (Hilborn and Walters 1992). They relate the current biomass to previous biomass plus terms for growth and recruitment, minus terms for natural mortality and catch. Surplus production models aggregate terms 5 for recruitment, growth, and natural mortality into one term for \surplus production" so that the biomass dynamics equations can be written in the form (Polacheck et al. 1993): Bt = Bt 1 + g(Bt 1) Ct 1 (1) where Bt, Ct, and g(Bt) denote biomass at the start of year t, catch during year t, and the surplus production function, respectively. The surplus production function is usually assumed to be non-negative with g(0) = g(K) = 0, where K is the carrying capacity (corresponding to the level of the stock biomass at equilibrium prior to commencement of the shery). The Schaefer (1954) form of the surplus production function is g(Bt 1) = rBt 1 1 Bt 1 K : (2) Substituting (2) in (1) gives a parsimonious model describing the annual biomass dynamics transitions with just the two parameters r, the intrinsic growth rate, and K: Bt = Bt 1 + rBt 1 1 Bt 1 K Ct 1: (3) Note that the annual catch is treated as a xed constant. A common, though simplifying assumption is that the relative abundance index is directly proportional to the biomass, i.e. It = qBt (4) with catchability parameter q. For a more detailed discussion on surplus production models, the interested reader is referred to Hilborn and Walters (1992, Chapter 8), Polacheck (1993), Millar and Meyer (1998a), and references therein. 6 2 Bayesian Nonlinear State-Space Model Polacheck et al. (1993) compare three commonly used statistical techniques for tting the model de ned by equations (3) and (4), process error models, observation error models, and equilibriummodels. None of these is capable of incorporating uncertainty present in both equations: natural variability underlying the annual biomass dynamics transitions (process error) and uncertainty in the observed abundance indices due to measurement and sampling error (observation error). This is possible, however, using a state-space model, as shown in Meyer and Millar (1998) and Millar and Meyer (1998a). Equations (3) and (4) are the deterministic versions of the stochastic state and observation equations. We assumed log-normal error structures and used a reparametrization (Pt = Bt=K) by expressing the annual biomass as a proportion of carrying capacity as in Millar and Meyer (1998) to speed mixing (i.e. sampling over the support of the posterior distribution) of the Gibbs sampler. The state equations are rewritten as: P1j 2 = eu1 ; PtjPt 1;K; r; 2 = (Pt 1 + rPt 1(1 Pt 1) Ct 1=K) eut ; t = 2; : : : ; N (5) and the observation equations are: ItjPt; q; 2 = qKPt evt; t = 1; : : : ; N; (6) where ut are iid normal with mean 0 and variance 2, and vt are iid normal with mean 0 and variance 2. Here, XjY denotes the conditional distribution of X given Y . A full Bayesian model consists of the joint prior distribution of all unobservables, here the ve parameters, K; r; q; 2; 2, and the unknown states, P1; : : : ; PN , and the 7 joint distribution of the observables, here the relative abundance indices I1; : : : ; IN . Bayesian inference is then based on the posterior distribution of the unobservables given the data. In the sequel, we will denote the probability density function of a parameter by p( ). We assume that the parameters K; r; q; 2; 2 are independent a priori. By a successive application of Bayes theorem and conditional independence of subsequent states, the joint prior density is given by p(K; r; q; 2; 2; P1; : : : ; PN ) = p(K)p(r)p(q)p( 2)p( 2) N Y i=1 p(PtjPt 1;K; r; 2): (7) A noninformative prior is chosen for q. Prior distributions for K; r; 2; 2 are speci ed using biological knowledge and inferences from related species and stocks as discussed in Millar and Meyer (1998): K lognormal( K = 5:04; K = 0:5162); r lognormal( r = 1:15; r = 0:8983); q / 1=q; 2 inverse-gamma(3:79; 0:0102); 2 inverse-gamma(1:71; 0:0086): For general guidelines on the choice of prior distributions for parameters in stock assessment models the reader is referred to Punt and Hilborn (1997). Because of the conditional independence assumption of the relative abundance indices given the unobserved states, the sampling distribution is p(I1; : : : ; IN jK; r; q; 2; 2; P1; : : : ; PN ) = N Y t=1 p(ItjPt; q; 2): (8) 8 Then, by Bayes theorem, the joint posterior distribution of the unobservables given the data, p(K; r; q; 2; 2; P1; : : : ; PN jI1; : : : ; IN), is proportional to the joint posterior distribution of all unobservables and observables, i.e. p(K; r; q; 2; 2; P1; : : : ; PN ; I1; : : : ; IN) = p(K)p(r)p(q)p( 2)p( 2) N Y i=1 p(PtjPt 1;K; r; 2) N Y t=1 p(ItjPt; q; 2) (9) i.e. the product of (7) and (8), the product of prior and sampling distribution. Thus, this joint distribution (9) is the crucial one for subsequent Bayesian inference. 3 Model Representation as DAG Now let us step back and look at a graphical representation of the full Bayesian model. This abstraction has the advantage that we can concentrate on the essential model structure without getting bogged down by all the details about densities. For any year t, let us represent all unobservables (K; r; q; 2; 2; Pt) and observables (It) as ellipses, and constants (Ct) as rectangles. A way to express the conditional independence assumptions is by drawing solid arrows between nodes (see Figure 1). Hollow arrows Figure 1 near here go to deterministic nodes, which are logical functions of other nodes. Pmed[t] is an example of a deterministic node as it is a function of the nodes K; r;Ct 1, and Pt 1. This renders a model representation as a directed acyclic graph (DAG) as all edges in the graph are directed and there are no cycles because of the conditional independence assumptions. Let V denote the set of all nodes in the graph. Direct predecessors of 9 a node v 2 V are called \parents", direct o spring the \children". The solid arrows indicate that given its parent nodes, each node v is independent of all other nodes except descendants of v. For instance, if we are in year t and know the biomass in year t 1 and the values of the parameters r, K, and 2, then our belief in Pt is independent of the biomass in previous years 1 to t 2, and the data of all other years except the current relative abundance index It. It is then easy to construct the joint probability distribution of all stochastic nodes using the graphical description of the conditional independence assumptions: p(V ) = Y v2V p(vjparents(v)): (10) For our speci c surplus production model (10) is the graph-theoretical version of (9). In this way, the DAG (Figure 1) assists in constructing the full Bayesian model. For further reading on conditional independence graphs and graphical chain models the interested reader is referred to Wermuth and Lauritzen (1990). 4 Bayesian Inference Using BUGS Let Vu denote the subset of unobservable nodes, and Vo the subset of observable nodes. Once p(V ) has been obtained from (10), a general technical di culty encountered in any application of Bayesian inference is calculating the high-dimensional integral necessary to nd the normalization constant in the posterior distribution of the unobservables given the data: p(VujVo) = p(Vu; Vo) p(Vo) = p(V ) R p(Vu; Vo)dVu : (11) 10 In our speci c example this would require an (N+5)-dimensional integration as we have to integrate over the unobservables K; r; q; 2; 2; P1; : : : ; PN . Calculating the marginal posterior distribution of any variable would require a subsequent (N + 4)-dimensional integration. High-dimensional integration problems can be solved via Markov chain Monte Carlo as reviewed in Gilks et al. (1996). The Gibbs sampler, a special MCMC algorithm, generates a sample from the posterior (11) by iteratively sampling from each of the univariate full conditional posterior distributions as explained in Meyer and Millar (1998). These univariate full conditional posterior distributions p(vjV nv), for v 2 Vu, can be easily constructed from the joint posterior distribution p(V ) in (10) by picking out those terms that depend on v: p(vjV nv) / p(vjparents(v)) Y v2parents(w) p(wjparents(w)): (12) This is facilitated by the graphical representation (Figure 1) as the full conditional posterior distribution of any node v depends only on its parents, children, and coparents. For instance, if v = Pt, then the full conditional posterior distribution of Pt, p(PtjK; r; q; 2 2; P1; : : : ; Pt 1; Pt+1; PN ; I1; : : : ; IN), is proportional to p(PtjPt 1;K; r; 2) p(Pt+1jPt;K; r; 2) p(ItjPt; q; 2): Here, the dependence of the deterministic nodes Pmed[t] and Imed[t] as logical functions of Pt; r;K, and q has been resolved. In this way, BUGS exploits the representation of the model as a DAG for constructing these full conditional posterior distributions for all unobservable nodes. Once this is accomplished, it uses certain sophisticated sampling methods to sample from these univariate densities. BUGS contains a small 11 expert system for choosing the best sampling method. The rst choice is to identify conjugacy, where the full conditional reduces analytically to a well-known distribution, and sample accordingly. If the density is not conjugate but turns out to be log-concave, it employs the adaptive rejection sampling (ARS) algorithm (Gilks and Wild 1992). If the density is not log-concave, BUGS uses a Metropolis-Hastings (MH) step. The MH algorithms di er across the various BUGS versions and platforms. The current UNIX version 0.6 uses the Griddy Gibbs sampler as developed by Ritter and Tanner (1992). More e cient MH implementations based on slice sampling (Neal 1998) if the variable has a restricted range and adaptive techniques (Gilks et al. 1998) for a variable with unrestricted range are currently under development; a rst version has been released under WinBUGS, the BUGS version for the WINDOWS operating system. 5 Model Implementation in BUGS For a typical BUGS run, four di erent les have to be speci ed: 1. The data are entered in a le with extension .dat, here called tuna.dat. 2. Initial values for all unobservables needed to start the Gibbs sampler are in a .in le, here surplus.in. 3. The le that contains the model descriptions has the .bug extension, here surplus.bug. 12 4. The commands for running BUGS that control the number of sampled values, which parameters to monitor and so on, go into the surplus.cmd le. InWinBUGS, these commands are options in various menus. These four les are listed in the appendix. The data le can be either in rectangular format like tuna.dat or in SPLUS format as in tuna S.dat. The format of the initial value le follows that of the data le, here for instance surplus.in in SPLUS format. If initial values for parameters are not given in the initial value le, then values will be generated by forward sampling from the prior distributions speci ed in the model. The DAG representation of our model not only serves BUGS-internal purposes but can assist us in specifying our model in the BUGS language. The le containing the model speci cation, surplus.bug consists of two sections: the declarations and the model description. The declaration part speci es the name of the model, which nodes are constants (rectangles) and which are stochastic (ellipses), and gives the name of the les containing the data and initial values. The model description forms a declarative representation of the full Bayesian model. It is enclosed in curly brackets f: : :g. This is a translation of the graphical model into BUGS syntax. Each statement consists of a relation of two kinds: which means \is distributed as" and substitutes the solid arrows, < which means \is to be replaced by" and substitutes the hollow arrows. 13 Quantities on the left of a are stochastic, those on the left of < are deterministic. In general, each quantity should appear once and only once on the left-hand-side of a statement. The order of the expressions within a pair of braces is irrelevant. WinBUGS even has an option Doodle that allows the user to specify the model graphically by drawing a DAG. It uses a hyper-diagram approach to add extra information to the graph to give a complete model speci cation. DoodleBUGS then writes the corresponding model in BUGS syntax to a le. Note that BUGS uses a nonstandard parametrization of distributions in terms of the precision (1=variance) instead of the variance. The density function of a dlnorm( ; ) distribution, for instance, is f(xj ; ) = q 2 1 xe 1 2 (logx )2 . Thus we need the nodes isigma2 and itau2. As we can not specify a function of certain variables as a parameter of a distribution, we need the deterministic nodes Pmed[t] and Imed[t]. BUGS only allows to specify proper prior distributions. The noninformative distribution for q is improper, though. Therefore, we use a gamma(0.001,0.001) prior for the inverse of q which is practically equivalent to q / 1q . Furthermore, we have to restrict the ranges of those nodes with non-logconcave full conditional distributions by specifying lower and upper bounds using the I(lower,upper) function. This is not necessary in the WinBUGS version because of the adaptive MH algorithm. Finally, the le surplus.cmd compiles the BUGS commands in surplus.bug, generates an initial 1000 iterations (the so called \burn-in" period), monitors every speci ed parameter for the next 22500 iterations, stores every 25th value, and calculates summary statistics of the sampled values for each speci ed parameter. Our preference is 14 to submit these commands as a batch job by using the command backbugs "surplus.cmd" with session output automatically directed to the bugs.log le. Alternatively, BUGS can be run interactively by using the command bugs. BUGS generates 5 les after completion: 1. The le bugs.log contains the BUGS code (surplus.bug) that was used, possibly error messages, the running time, and the requested summary statistics of the marginal posterior distribution of each parameter. The posterior summaries from this le are listed in the appendix. Note that the results from tting this surplus production model are presented and discussed in Millar and Meyer (1998a). 2. The le bugs.out contains 2 columns. The rst column gives the iteration number, the second column the corresponding sampled value. 3. The le bugs.ind tells you which line of the bugs.out le corresponds to which monitored variable. Here, lines 1 to 9000 of bugs.out are samples from variable K, lines 9001 to 18000 are samples from variable r and so on. 4. The le bugs1.out contains the results of the stats command in a rectangular format suitable for reading into statistical packages for producing graphs, tables etc. The 10 columns contain the summary statistics: mean, sd, 2.5%, observed lower percentile, 97.5%, observed upper percentile, median, number of iterations, start iteration, nish iteration. 15 5. The le bugs1.ind contains the node names for the variables listed, and the corresponding row number in the bugs1.out le. Whenever one employsMCMCmethods to sample from the posterior distribution, a question of extreme importance is \Has the chain converged to its target distribution?". Convergence diagnostics is still an area of active research. A recent review on methods used for establishing whether an MCMC algorithm has converged, i.e. whether its output can be regarded as samples from the target distribution of the Markov chain, has been given by Cowles and Carlin (1995). Some of these methods are implemented in the software CODA (Best et al. 1995). CODA is a menu-driven collection of SPLUS functions for analysing the output obtained by BUGS. Besides trace plots and the usual tests for convergence, CODA calculates statistical summaries of the posterior and kernel density estimates. CODA is being maintained and distributed by the same research group responsible for BUGS. 6 Discussion BUGS provides a means for a Bayesian analysis of the most complex sheries models such as fully age-structured models (Millar and Meyer 1998b), fully Bayesian hierarchical random e ects models for stock-recruitment data (Liermann and Hilborn 1997), and for tackling the important Bayesian hierarchic meta-analysis problems identi ed in the \future directions" of Hilborn and Liermann (1998). BUGS will have considerable implications for the day-to-day practice of the sheries stock assessment scientist. 16 She can concentrate on the essential realistic modeling of data and of prior knowledge, handing over the mathematical intricacies of model tting to a software package. Because prior and sampling distributions can be easily modi ed without having to recalculate full conditional distributions, the scientist's time is freed up to experiment with di erent scenarios a hitherto far too time-consuming procedure. The recent evolution in Bayesian computation and software within the last few years parallels and outperforms in some ways the development seen over the last decades of frequentist software packages based on maximizing the likelihood function. Once stable numerical routines for nonlinear optimization (like the Newton-Raphson algorithm) had been developed, these were implemented as modules in software packages (like SAS, GLIM, or SPLUS) for maximizing the likelihood function. Fitting a complicated model became possible for the statistical layperson just by running a certain DATA and PROC step in SAS for instance. Knowledge about mathematical di erentiation and optimization was no longer required. The inevitable di culty in calculating the high-dimensional integrals necessary for the determination of posterior probability distributions has hindered routine Bayesian data analysis. For scientists applying Bayesian inference a good grasp of numerical and Monte Carlo integration techniques has been essential. Until recently, all applications have involved writing one-o computer code in lowor intermediate-level languages such as C or Fortran. Even with reliable Markov chain Monte Carlo subroutines like the Metropolis-Hastings algorithm and the Gibbs sampler to overcome the multidimensional integrations, writing and debugging a speci c program would have taken 17 anything from a few days to several weeks. And a subsequent modi cation of the pro-gram, as e.g. an extension of the model or an application to a di erent dataset, mightwell have taken several hours. So Bayesian computations have clearly proved to beinferior when compared with model tting in SAS or SPLUS for which model speci-cation may take just a few minutes, and modi cations may take just a few seconds.Therefore, the Bayesian software package BUGS has been developed with the Gibbssampler and Metropolis-Hastings algorithms as building blocks. Instead of having tochoose the right PROC in SAS for the problem at hand, with di erent commands andoptions for each procedure, the user of BUGS only has to specify prior and samplingdistributions for his/her model. In this respect, the Bayesian software package BUGSis considerably more exible, general, and user-friendly. Moreover, as demonstrated inthe context of state-space models, the Bayesian approach in conjunction with the Gibbssampler is capable of handling far more complex models than the classical approach.The progress made in Bayesian computation via MCMC methods has alreadyrevolutionized many scienti c disciplines and provided solutions to problems whichwere hitherto considered computationally intractable. The routine implementation ofBayesian inference that is now possible will \almost surely" have an impact on sheriesstock assessment.18 ReferencesBest, N.G., Cowles, M.K., and Vines, S.K. 1995. CODA manual version 0.30. Cam-bridge, UK: MRC Biostatistics Unit.Carlin, B.P., Polson, N.G., and Sto er, D.S. 1992. A Monte Carlo approach to non-normal and nonlinear state-space modeling. J. Amer. Statist. Assoc. 87: 493{500.Cowles, M.K., and Carlin, B.P. 1995. Markov chain Monte Carlo convergence diagnos-tics: a comparative review. Technical Report, University of Minnesota.Fahrmeir, L., and Tutz, G. 1994. Multivariate Statistical Modelling Based on Gener-alized Linear Models. Springer, New York.Freeman, S.N., and Kirkwood, G.P. 1995. On a structural time series method for esti-mating stock biomass and recruitment from catch and e ort data. Fish. Res. 22,77-98.Gentleman, R. 1997. A Review of BUGS: Bayesian Inference Using Gibbs Sampling.Chance 10, 48-51.Gilks, W.R., and Wild, P. 1992. Adaptive rejection sampling for Gibbs sampling. Ap-plied Statistics 41: 337{48.Gilks, W.R., Thomas, A., and Spiegelhalter, D.J. 1994. A language and program forcomplex Bayesian modelling. The Statistician 43: 169{178.19 Gilks, W.R., Best, N.G., and Chan, K.K.C. 1995. Adaptive rejection Metropolis sam-pling within Gibbs sampling. Appl. Statist. 44: 455{472.Gilks, W.R., Richardson, S., and Spiegelhalter, D.J. 1996. Markov Chain Monte Carloin Practice. Chapman & Hall, London.Gilks, W.R., Roberts, G.O., and Sahu, S.K. 1998. AdaptiveMarkov Chain Monte Carlothrough regeneration. MCMC preprint server.Gudmundsson, G. 1994. Time series analysis of catch-at-age observations. Appl. Statist.43, 117-126.Hilborn, R., and Walters, C.J. 1992. Quantitative sheries stock assessment: choice,dynamics and uncertainty. Chapman & Hall, New York.Hilborn, R., and Liermann, M. 1998. Standing on the shoulders of giants: learningfrom experience in sheries. Rev. Fish. Biol. Fisheries 8, 1-11.Kinas, P.G. 1996. Bayesian shery stock assessment and decision making using adap-tive importance sampling. Can. J. Fish. Aquat. Sci. 53: 414{423.Liermann, M., and Hilborn, R. 1997. Depensation in sh stocks: a hierarchic Bayesianmeta-analysis. Can. J. Fish. Aquat. Sci. 54: 1976-1983.Meyer, R., and Millar, R.B. 1998. Bayesian stock assessment using a state-space im-plementation of the delay di erence model. Can. J. Fish. Aquat. Sci. xx: xxx{xxx.20 Millar, R.B., and Meyer, R. 1998a. Nonlinear state-space modeling of sheries biomassdynamics using the Gibbs sampler. Technical Report STAT9802, Department ofStatistics, The University of Auckland.Millar, R.B., and Meyer, R. 1998b. Future opportunity and challenge for sheries mod-elers. Technical Report STAT9804, Department of Statistics, The University ofAuckland.Neal, R.M. 1997. Markov Chain Monte Carlo methods based on `slicing' the densityfunction. Technical Report No. 9722. Department of Statistics, University ofToronto.Pella, J.J. 1993. Utility fo structural time series models and the Kalman lter for pre-dicting consequences of shery actions. In Proceedings of the international sym-posium on management strategies for exploited sh populations. Edited by G.Kruse, D.M. Eggers, R.J. Marasco, C. Pautzke, and T.J. Quinn II. Alaska SeaGrant College Program Report 93-02, Univ. of Alaska, Fairbanks. pp. 571{593.Polacheck, T., Hilborn, R., and Punt, A.E. (1993). Fitting surplus production models:comparing methods and measuring uncertainty. Can. J. Fish. Aquat. Sci. 50,2597-2607.Punt, A.E., and Hilborn, R. 1997. Fisheries stock assessment and decision analysis:the Bayesian approach. Rev. Fish. Bio. Fish. 7: 35{63.Reed, W.J., and Simons, C.M. 1996. Analyzing catch-e ort data by means of the Kalman21 lter. Can. J. Fish. Aquat. Sci. 53: 2157{2166.Ritter, C., and Tanner, M.A. 1992. Facilitating the Gibbs sampler: the Gibbs stopperand the griddy-Gibbs sampler. Journal of the Royal Statistical Society, Series B,59: 291{317.Schaefer, M.B. 1954. Some aspects of the dynamics of populations important to themanagement of commercial marine sheries. Inter-Am. Trop. Tuna Comm. Bull.1: 27{56.Schnute, J.T. 1994. A general framework for developing sequential sheries models.Can. J. Fish. Aquat. Sci. 51: 1676{1688.Sullivan, P.J. 1992. A Kalman lter approach to catch-at-length analysis. Biometrics48: 237{257.Spiegelhalter, D.J., Thomas, A., Best, N., and Gilks, W.R. 1996. BUGS 0.5, Bayesianinference using Gibbs sampling. Manual (version ii) Cambridge, UK: MRC Bio-statistics Unit.Thomas, A., Spiegelhalter, D.J., and Gilks, W.R. 1992. BUGS: A program to performBayesian inference using Gibbs sampling. In Bayesian Statistics 4. Edited by J.M.Bernardo, J.O. Berger, A.P. Dawid, and A.F.M. Smith. Oxford University Press,Oxford. pp. 837{42.Wermuth, N., and Lauritzen, S.L. 1990. On substantive research hypothesis, conditionalindependence graphs and graphical chain models (with discussion). J. Roy.22 Statist. Soc. B, 52, 21-72.Yeh, S.-Y., Tsou, T.-S., and Liu, H.-C. (1991). Assessment of the south Atlantic al-bacore resource by using surplus production models, 1976-1988. Colln. Vol. Sci.Pap. ICCAT 39, 166-170.23 Appendixtuna.dat15.9 61.8925.7 78.9828.5 55.5923.7 44.6125.0 56.8933.3 38.2728.2 33.8419.7 36.1317.5 41.9519.3 36.6321.6 36.3323.1 38.8222.5 34.3222.5 37.6423.6 34.0129.1 32.1614.4 26.8813.2 36.6128.4 30.0734.6 30.7537.5 23.3625.9 22.3625.3 21.91tuna S.datlist(C = c(15.9,25,7,...,25.3), I = c(61.89,78.98,...,21.91))surplus.inlist(P=c(0.99,0.98,0.96,0.94,0.92,0.90,0.88,0.86,0.84,0.82,0.80,0.78,0.76,0.74,0.72,0.70,0.68,0.66,0.64,0.62,0.60,0.58,0.56),r=0.8, K=200, iq=5, isigma2=100, itau2=100)24 surplus.bug:model surplusproduction;const N=23;var C[N], I[N], Imed[N], P[N], Pmed[N],r, K, q, iq, sigma2, isigma2, tau2, itau2;data C, I in "tuna.dat";inits in "surplus.in";{# prior distribution of K : lognormal with 10% and 90% quantile at 80 and 300K ~ dlnorm(5.042905,3.7603664)I(1,2000);# priordistribution of r: lognormal with 10% and 90% quantile at 0.1 and 1.0r ~ dlnorm(-1.151293,1.239084233)I(0.001,2.0);# prior distribution of q: instead of improper (prop. to 1/q) use just proper IGiq ~ dgamma(0.001,0.001)I(0.5,200);q <1/iq;# prior distribution of sigma2: inv. gamma with 10% and 90% qu. at 0.04 and 0.08isigma2 ~ dgamma(3.785518,0.010223);sigma2 <1/isigma2;# prior distribution of tau2: inv. gamma with 10% and 90% qu. at 0.05 and 0.15itau2 ~ dgamma(1.708603,0.008613854);tau2 <1/itau2;# (conditional) prior distribution of Ps (from state equations):Pmed[1] <0;P[1] ~ dlnorm(Pmed[1],isigma2) I(0.001,2.0)for (t in 2:N) { Pmed[t] P[t]~ dlnorm(Pmed[t],isigma2)I(0.001,2.0)}#sampling distribution:for (t in 1:N) { Imed[t] I[t] ~ dlnorm(Imed[t],itau2);}}25 surplus.cmdcompile("surplus.bug")update(1000)monitor(K,25)monitor(r,25)monitor(q,25)monitor(sigma2,25)monitor(tau2,25)monitor(P[],25)update(22500)stats(K)stats(r)stats(q)stats(sigma2)stats(tau2)stats(P[])q()bugs.logBugs>stats(K)meansd2.5% : 97.5% CI mediansample2.619E+2 6.903E+1 1.452E+2 4.193E+2 2.556E+29000Bugs>stats(r)meansd2.5% : 97.5% CI mediansample3.287E-1 1.160E-1 1.436E-1 6.184E-1 3.121E-19000Bugs>stats(q)meansd2.5% : 97.5% CI mediansample2.622E-1 7.772E-2 1.455E-1 4.708E-1 2.471E-19000Bugs>stats(sigma2)meansd2.5% : 97.5% CI mediansample3.204E-3 1.988E-3 1.137E-3 8.458E-3 2.689E-39000Bugs>stats(tau2) meansd2.5% : 97.5% CI mediansample1.240E-2 4.667E-3 5.668E-3 2.393E-2 1.160E-29000Bugs>stats(P[]) meansd2.5% : 97.5% CI mediansample[1]1.019E+0 5.475E-2 9.183E-1 1.137E+0 1.016E+09000[2]9.945E-1 7.636E-2 8.686E-1 1.173E+0 9.847E-19000[3]8.735E-1 6.702E-2 7.522E-1 1.018E+0 8.693E-19000...[23]3.253E-1 4.022E-2 2.530E-1 4.118E-1 3.234E-1900026 Table 1 : Catch (in thousand tonnes) and CPUE (in kg per 100 hooks) datafor South Atlantic albacore tuna from Polacheck (1993)Year Catch CPUE196715.9 61.89196825.7 78.98196928.5 55.59197023.7 44.61197125.0 56.89197233.3 38.27197328.2 33.84197419.7 36.13197517.5 41.95197619.3 36.63197721.6 36.33197823.1 38.82197922.5 34.32198022.5 37.64198123.6 34.01198229.1 32.16198314.4 26.88198413.2 36.61198528.4 30.07198634.6 30.75198737.5 23.36198825.9 22.36198925.3 21.9127 Figure 1 : Representation of the surplus production model as a directedacyclic graph (DAG).