In haddonm/datalowSA: A Package to Faciliate Application of Data poor Stock Assessments

knitr::opts_chunk$set(
    echo = TRUE,
    message = FALSE,
    warning = FALSE)
options(knitr.kable.NA = "",
        knitr.table.format = "pandoc")

options("show.signif.stars"=FALSE,"stringsAsFactors"=FALSE,
        "max.print"=50000,"width"=240)

library(datalowSA)
suppressPackageStartupMessages(library(knitr))

Introduction

Which Stock Assessment?

Which stock assessment method to apply to fisheries for data-poor to data-moderate species will depend upon what fisheries and biological data are available but also, importantly, on what management objectives need to be met within the jurisdiction in question. It may be the case that the fishery for a particular species is of sufficient size and value to warrant on-going monitoring and management towards some defined goal for the stock. In such a case the assessment used should obviously be capable of generating some notion of the current state of the fishery and indicate what management actions may be required to eventually achieve the agreed management goals. But some fisheries may be so minor that trying to actively manage them would be inefficient both practically and economically. Nevertheless, to meet the requirements of the Status of key Australian Fish Stocks (SAFS) one still requires some form of defensible stock assessment capable of determining whether the current level of fishing is sustainable.

Modified Catch-MSY

The Catch-MSY method (Martell and Froese, 2013) could be termed a 'model-assisted' stock assessment method. It only requires a time-series of catches and a set of strong assumptions to conduct a stock assessment. As only a brief description of how it is considered to work is given here, it is recommended that users read the original paper to gain an understanding of what the method does and how it does it.

The underlying stock dynamics are described by the simple model used, which in the case implemented here is a Schaefer surplus production model with parameters r, the population growth rate, and K, the population carrying capacity or unfished biomass. The model uses ratios of the initial and final catches relative to the maximum catch to set up arrays of potential values for the initial and final depletion levels as well as for the potential range of r and K values (all of which are now modifiable by the user). The method sequentially steps through the years of the fishery by randomly selecting pairs of r-K values from the wide initial ranges, which defines the initial biomass, subtracting the catches, and moving the population dynamics forward each year using the predictions from the simple model. Essentially this is a stock reduction that removes catches from a known set of dynamics. However, the very many r-K pairs used (at least 20000) are combined with a fixed set of initial depletion levels (about 20 steps between the minimum and maximum initial depletion set) to generate often 100,000s of possible stock reduction trajectories. Criteria are included (e.g. no trajectory is kept if it predicted zero biomass or biomass above K) that lead to numerous potential trajectories being rejected. Those that are left after all criteria for acceptance have been completed constitute the set of trajectories deemed to be consistent with the known catches. The implications of these successful trajectories are used to produce an assessment of the possible status of the stock.

In this section we will describe how to conduct a catch-MSY analysis, how to extract the results from that analysis, as well as plot out illustrations of those results. In addition, we will examine how to project the successful trajectories under constant catch scenarios to determine what level of catches should lead to the majority of trajectories moving in a desired direction (rebuilding to a target, staying stable, or declining to a target depletion).

A standard workflow might consist of:

1. read in data and use checkdata to determine which analyses are possible
2. set the run-time parameters to conduct the desired catch-MSY analysis (see later)
3. use run_cMSY to conduct a catch-MSY analysis. One sensitivity you should run is to set sigpR (the proxy for process error) to 1e-10 (a very small number) to see deterministic trajectories. Do other sensitivities.
4. use cMSYphaseplot to determine status for SAFS, if more is wanted:
5. use summarycMSY to generate a summary object of the answer.
6. use plottrajectory to illustrate either all successful trajectories along with the successful harvest rates (use option oneplot = TRUE), or just a sample of 7 or 15 biomass trajectories.
7. use plotcMSY6 to plot up the overview of successful and failed r - K combinations.
8. use pulloutStats to obtain summary statistics regarding MSY and depletion
9. use plotconstC with $deplet from the object output from pulloutstats to plot up the successful trajectories depicting stock depletion.
10. use doconstC to conduct constant catch projections to find a level of catch that lead, on average, to the population decreasing, staying stable, or increasing (possibly to some selected target within a particular time).
11. optionally, use trendMSY to illustrate the relationship between the estimated MSY and the average MSY of the successful trajectories.

Some Formal Details

The Catch-MSY method described here can be regarded as a model-assisted data-poor method. It uses a form of stock reduction analysis where the productivity of a given stock (its unfished biomass and its reproductive rate) is characterized within the parameters of a simple mathematical model, and how that modelled stock responds to the history of known catches (a stock reduction analysis) forms the basis of the alternative methods used to characterize productivity in management useable terms.

The Catch-MSY method (Martell and Froese, 2013) uses the relatively simple Schaefer surplus production model as the basis for describing the dynamics of the stock being described.

$${B}{t+1}={B}{t}+r{B}{t}\left( 1-\frac{{B}{t}}{K} \right)-{{C}{t}}$$ where _B~t~ represents the stock biomass in year t, r represents a population growth rate that includes the balance between recruitment and natural mortality, K is the maximum population size (the carrying capacity), and C~t~ being the catch in year t. The $\left( 1-\frac{{B}{t}}{K} \right)$ represents a density dependent term that trends linearly to zero as _B~t~ tends towards K.

Importantly, for our purposes, one of the properties of the discrete Schaefer surplus production model is that MSY can be estimated very simply from the parameter estimates:

$${MSY}=\frac{rK}{4}$$ which reflects the symmetric production function implied by the model dynamics. A relatively simple future possible development would be to include the option of using Fox model dynamics instead of the Schaefer.

Such surplus-production (or biomass dynamic) models usually require both a time-series of total catches (landings plus discards) and a time-series of an index of relative abundance (Haddon, 2011). In Australia the index of relative abundance is most often a time-series of CPUE (ideally standardized CPUE). For the truly data-poor fisheries only a catch series may be available and hence the nature of this stock reduction method entails exploring the parameter space for plausible trajectories rather than fitting a model to data to find the most statistically likely trajectory.

Empirical Harvest Strategies

In the Southern and Eastern Scalefish and Shark Fishery (SESSF), rather than using surplus production models or other simple approaches that attempt to model the underlying population dynamics of a stock, empirical harvest strategies have been developed that use such time-series in empirical relationships that give rise directly to management related advice on catch levels (Little et al., 2011; Haddon, 2014). Such empirical harvest strategies can provide the needed management advice but do not determine stock status unless the reference period, often used in such approaches, is assumed to be a proxy for the target reference point (and associated limit reference point) for sustainability. A weight-of-evidence, often qualitative argument would need to be made to support the use of such a proxy. In the SESSF, this so-called Tier 4 harvest strategy is used to determine whether a stock is over-fished or not but currently cannot be used to determine whether over-fishing is occurring. In addition, there is the strong assumption made that commercial catch-rates are a direct reflection of the stock biomass. There are, however, some species, for example mirror dory (Zenopsis nebulosa) where catch rates increase when catches increase and then decline once catches begin to decline. They appear to be fisheries based on availability rather than the fishery being the major influence on the stock biomass, and other aspects of the environment of the species appear to be driving its dynamics. The use of CPUE may thus be misleading in such cases or at best lead to simply reactive management decisions (cpue goes up so can catches, cpue goes down so must catches).

More widely than the SESSF there are many fisheries within Australia that may only have a time-series of catches with only limited information related to a usable index of relative abundance. In addition, such catch time-series may not be available from the beginning of the fishery, which means that methods such as Depletion-Based Stock Reduction Analysis (Dick and MacCall, 2011) cannot be validly applied (although, as shown in Haddon et al, 2015, if sufficient years of catches are present (perhaps >25) then the method can still provide approximate estimates of management related parameters). Under such data-limited situations other catch-only based assessment methods can provide the required estimates of management interest.

Stock Reduction Analyses

As with many of the more capable catch-only data-poor approaches the Catch-MSY method evolved from the stock reduction analyses of Kimura and Tagart (1982), Kimura et al. (1984), and eventually Walters et al. (2006). It uses a discrete version of the Schaefer surplus production model (Schaefer, 1954, 1957; Haddon, 2011) to describe the stock dynamics in each case. The Catch-MSY requires a time-series of total removals, prior ranges for the r and K parameters of the Schaefer model, and possible ranges of the relative stock size (depletion levels) in the first and last years of the time-series. As described by Martell and Froese (2013) the range of initial depletion levels can be divided into a set of initial values, and a stock reduction using the known total removals, applied to each of these multiple initial depletion levels combined with pairs of r-K parameters randomly drawn from uniform distributions across the prior ranges of those parameters. Each of these parameter pairs plus each of the initial depletion levels are projected using the total catch trajectory leading to a stock biomass trajectory which is either accepted or rejected depending on whether the stock collapses or exceeds the carrying capacity, and whether the final depletion level falls within the assumed final range.

The initial and final depletion ranges can be relatively broad. Other criteria can be included to further constrain the biomass trajectories if extra evidence is available. Such additional constraints are still under development. For example, in some of the examples you will notice that the annual harvest rates for some accepted trajectories can be very high (> 0.5), which for many (though not all) Australian species can be considered to be implausible. Now it is possible to conduct a sensitivity analysis where trajectories implying some pre-defined harvest rate will also be rejected. These high fishing mortality trajectories are only possible for the more productive parameter combinations so removing such trajectories will likely reduce the predicted MSY (maximum productivity).

We can use the invert data set of catches to exemplify the process of applying the Catch-MSY method.

The data-file format

As can be seen in the R-code above, the data for the catch_MSY analysis is read into an R object containing fish, glb, and others if present. Using the data command places the data object, in the example it is called invert, into the global environment. Using the readdata command with a "named.csv" file also produces one large object. The format of the standard data files is given in three example appendices at the end of this document. The actual data requirements of the different methods available differs greatly. The catch-MSY method really only requires two columns in fish one being "year" the other being "catch". The surplus production modelling also requires a "cpue" column, and the age-structured surplus production model also requires various biological properties relating to weight-at-age, maturity-at-age, and selectivity-at-age. These are tested by the checkdata function

Table 1. The format of the fish object. Not all fields are required for all analyses.

data(invert)
fish <- invert$fish   # contains available fishery data
head(fish,15)

#library(datalowSA)

glb <- invert$glb     # contains available biological data
checkdata(invert)
              # normally one would run at least 20000 iterations, preferably more
reps <- 5000  # read the help for run_cMSY to understand each input parameter
              # try changing the value of sigpr and makeplots to TRUE                 
answer <- run_cMSY(fish,glb,n=reps,sigpR=0.025,maximumH=1.0)
str(answer,max.level=1)

The use of checkdata indicates that both catch-MSY and spm analyses are possible with this data. This does not mean such analyses will always be valid with the given data only that the required data to conduct these analyses are present.

The answer object contains all the results from the catch-MSY analysis and is used by other R functions to generate summaries and plots of those results.

SAFS Status

A summary and illustration of the stock status can be obtained by extracting the average stock biomass trend, along with the average fishery harvest rate trend from the successful trajectories.

Using the average estimate of the r - K pairs it is also possible to generate an estimate of the production curve, from which it is possible to derive estimates of B~msy~, 0.2K or 0.2B0, and H~targ~ and H~lim~, with which the phase plot can be subdivided to provide a visual representation of the how the history of catches from the fishery are reflected in predicted changes in the biomass and harvest rate. Whether over-fishing is occurring (leading to a status of 'depleting') is determined by whether the current point lies above the H~targ~ or above the H~lim~, which in turn is decided by the management objectives adopted in each jurisdiction. A status of 'depleted' or sustainable currently corresponds to whether the current year's point is to the left or right of the 0.2B~0~ line.

out <- cMSYphaseplot(answer,fish)

Figure 1. A phase plot of the man predicted biomass and harvest rates through the years observed. The first year of data is a green point and the last a red point. The axes for the bottom plot are identified through the colour of the axis titles

The results are synthesized within the contents of out.

str(out)

Other Potential Outputs from catch-MSY

A useful summary of the catch-MSY analysis containing the primary results can be obtained using the summarycMSY function.

summcMSY <- summarycMSY(answer,fish,final=TRUE)
str(summcMSY,max.level=1)  # try max.level = 2

As can be seen from the series of large objects in summcMSY, such as the r, K, and msy components, only a portion of the n=5000 trials succeeded in meeting the constraints that define an acceptable trajectory. This has removed all the parameter combinations that were not productive enough or those that were too productive and retained only those that were at least realistic/plausible enough to be consistent with what is known about the fishery. If you run the last three lines of R code a few times and examine the structure of the summcMSY object you will notice that each time the length of the r, K, and msy objects usually differs. This merely exemplifies the fact that the randomly selected parameter combinations lead to different numbers of successful trajectories each time through. To obtain a visual representation of a selection of the successful trajectories you can use the function plottrajectory (see its help function for a description of all the options).

It should be noted that the catch-MSY method uses a two stage strategy. Unless set otherwise, it first sets the initial K values between the maximum catch and 60 times the maximum catch. This invariably leads to very many successful trajectories that suggest a very large initial biomass combined with a very low population growth rate. While mathematically this may match the productivity of the stock, as suggested by the time-series of catches, it is biologically less plausible than lower K values associated with higher r values (these two parameters are negatively correlated). To avoid the less plausible combinations in their original code Martell & Froese (2013) search for the smallest K value that will still give rise the overall average MSY value across the successful trajectories found in the first run through. Once a more restricted initial K range has been found then the n replicates are repeated and the final result put into the R1 object inside answer. The results from the first run through the replicates are put inside the Rfirst object.

out <- plottrajectory(answer$R1,fish$year,fish$catch,answer$parbound,
                      oneplot=FALSE,Bmax=25000,
                      scalar=1.0,plotout=TRUE,plotall=15)

Figure 2. 15 examples of different r-K combinations only illustrating those trajectories that were plausible, the number at the top of each plot is the predicted MSY for the given parameter pair. The r and K values are the y-axis labels in each case. The final plot is the catch history with the predicted average MSY as the blue line.

\newline

Once the summcMSY object has been generated the results can be summarized by using the plotcMSY6 function. This generates a plot of the combinations of r and K and whether they succeeded or not. In addition it displays the distribution of the successful r, K, and MSY values.

plotcMSY6(summcMSY,fish[,"catch"],label=glb$spsname)

Figure 3. A plot of the catch trajectory with the MSY and its 90% percentiles. The two colour plots are a plot of the K vs r combinations with the red dots depicting failure and then the increasing colours depicting more combinations of initial depletion that succeeded for each r-K pair. The right-hand plot is the log-transformed version of the left-hand plot. The histograms describe the distributions of the successful r-K pairs and the resulting MSY. The red lines are the median and the 90th percentile confidence intervals.

Interpretation of the plotcMSY6 plot

The histograms along the bottom row of Figure 3 illustrate the distribution of the successful values of each parameter and the resulting MSY from the parameter pairs. Because of the truncation of the K values in the two phase operation of the analysis it should not be surprising if the 'Successful K' histogram is flat up against the right hand bound. However, if either of the r or left-hand side of the K plots is hard up against the bounds, it would be worthwhile exploring whether different initial values of the r and K values are less constrained. This can be done using a different resilience value or by explicitly altering the start_r and the start_K values. Around the jurisdictions examples were seen where parts of the r or K histograms were flat topped, which suggests there is insufficient information in the catches to adequately discriminate between potentially successful trajectories and failures.

Reading Martell and Froese (2013) one can see that the method relies on the notion that a fishery is expected to begin small, build up catches, and then those catches lead to depletion so that they cannot be maintained and the catches drop. Of course, in Australia there can be many reasons for catches to drop other than the stock becoming depleted. There may be management interventions that lead to reduced catches (imposed catch limits, large marine closures, gear restrictions, and other management steps). In such situations where catches are reduced before stock depletion leads to reduced catches the catch-MSY seems likely to generate conservative estimates of sustainable production. Other issues that can arise are that catches have only ever increased or stayed stable. Without an on-going decline in catches at some period in the fishery, the model used inside the catch-MSY will have no may in which to characterize maximum productivity and it may give simply an approximate estimate of the average catches. Obviously care is needed when interpreting the outputs from the catch-MSY analysis.

The Catch-MSY method requires a catch history and some prior notion of the relative resilience or expected productivity of the species being fished. The basic idea is that the method begins with a given range of initial depletion and final depletion levels along with initial ranges of the two parameters describing the Schaefer model; these are the r and K parameters representing the un-restricted population growth rate and the population carrying capacity respectively. The catch-MSY method is based on a review of very many catch histories and relies on an expectation of catches increasing as a fishery develops and then decreasing as the catches impact the stock. The statistics can begin at the start of a fishery or after it has already developed, hence the initial depletion is set by comparing the initial catches with the maximum catches. The initial ranges for the model parameters depend on the assumed productivity or resilience of the species of which four options exist - "verylow", "low", "medium", and "high", with associated r ranges of (0.015 - 0.125), (0.1 - 0.6), (0.3 - 0.8), and (0.6 - 1.5). In the original code associated with the Martell and Froese (2013) paper they only used three resilience categories, omitting the "medium". The initial range of the unfished biomass is also very broad with a minimum set at the maximum catch and a maximum set at 60 times the maximum catch.

In case these particular combinations do not suit what is known or suspected of a particular species, the run_cMSY function now includes the option to set your own limits on both the initial values for r and K. A call to either formals(run_CMSY) or to ?run_cMSY will illustrate the syntax and the names of start_r and start_K. In each case they require a short vector of two numbers (e.g. start_r=c(0.015,0.3))

Results from the Catch-MSY Analysis

The plots above illustrate the results and provide some information but it is also helpful to obtain a tabulation of the outcome of the analysis. For that we use the pulloutStats function. It generates another object containing a matrix of r, K, MSY, and current depletion percentiles and mean values. In addition, it contains the successful biomass trajectories and those same trajectories translated into depletion levels.

results <- pulloutStats(answer$R1)
print(str(results))

Table 2. The results$output statistics from the Catch-MSY analysis. The 2.5Perc and 97.5Perc are the respective percentile describing the spread of the trajectories (i.e. not confidence intervals around the mean). The % columns are the quantiles so that the 50% columns represents the medians.

round(results$output,2)

These results are fine as far as they go but in order to obtain some notion of stock status it is necessary to trace the successful trajectories in terms of how their depletion has changed through time. For this we can use the same plottrajectory function as used previously only this time changing the oneplot parameter to TRUE, which overrides the plotall parameter. The Bmax parameter can be adjusted to obtain an acceptable spread of the results.

out <- plottrajectory(answer$R1,fish$year,fish$catch,answer$parbound,
                      oneplot=TRUE,scalar=1.0,plotout=TRUE,plotall=7)

Figure 4. The top plot is of the successful biomass trajectories and the red line is the mean in each year. The bottom plot is of the annual harvest rate; note a few seemingly successful trajectories lead to harvest rates > 0.5. Associated with these high harvest rates note also that the final biomass for many trajectories is very low, which may also be deemed unrealistic.

Given each biomass trajectory has an associated K value it is possible to translate the biomass trajectories into depletion trajectories. There are two ways to do this, the simplest way is to use the pulloutStats function (see example above) which outputs both the biomass and depletion matrices as a byproduct to estimating the mean MSY and current depletion.

results <- pulloutStats(answer$R1)
# Note the use of constC=0, we are not doing any projections yet so no constant catches
effectC <- plotconstC(results$deplet,endyear=2017,constC=0,limit=0.2,target=0.4,
                      console=TRUE,intensity=NA,contours=TRUE) 
abline(h=c(0.2,0.3,0.4),col=3,lwd=2)

Figure 5. A plot of the successful depletion trajectories with the mean and median annual depletion marked. The lower red line is the default 0.2B0 limit reference point, while the upper is the input target reference point. The green line denotes the end of the final year in which data are available. If the abbreviated table under the plot is not wanted (all the information will be in effectC) then the console=TRUE should be set to console=FALSE. The red dashed lines are the inner 80th and 90th percentile bounds but if not wanted set contours=FALSE.

The PltLim% column in the console output is the proportion of trajectories that are < limit, whereas the PgtTarg% is the proportion of trajectories that are greater than the target. Given the great uncertainty in these analyses it should not be surprising that both of these increase, which suggests the total spread of the outcomes is increasing. The question that needs answering is whether the average is increasing or decreasing here both the mean and median are increasing though only about 4 - 5 % across the five years prior to the final year of data. As the average stock size increases the rate of increase would be expected to first increase and then decline as the stock moved on average above the biomass that produces the maximum productivity (B~MSY~).

The default plot of the depletion in Figure 5 has grey lines of equal density but if one wants to attempt to generate a plot with the density of colour matched to the density of trajectories then it is possible to include a number in the intensity parameter. The value input determines the density of the trajectories required to obtain full colour intensity and this will undoubtedly vary by fishery and so will need to be used interactively to find a value that generates a satisfactory plot. If the plot is to be printed then be sure to save the plot as a bit map or .png file so that the varying density and transparency is retained in the plot.

effectC <- plotconstC(results$deplet,endyear=2017,constC=0,limit=0.2,target=0.4,
                      console=FALSE,intensity=30) 
abline(h=c(0.2,0.3,0.4),col=3,lwd=2)

Figure 6. A plot of the successful depletion trajectories with the mean and median annual depletion marked with the density of trajectories represented by different intensity of colour. The lower red line is the default 0.2B0 limit reference point, while the upper is the input target reference point. The green line denotes the end of the final year in which data are available. The intensity value may need to be adjusted empirically to obtain a satisfactory plot.

In Figure 6 different intensity of colours are used to denote the paths followed by most trajectories. In this case between 1995 - about 2001 most lines are towards the top of the overall path but are moving downwards until between 2002 - 2008 the median line is biased low with it being closer to the bottom of all paths than the top. After that most paths are moving upwards with the median once again shifting closer to the upper margin than the lower.

Constant Catch Projections

From the table and the plot it is clear that the average depletion is quite a way below the selected target. We need to remember that this is a data-poor assessment and that currently there are no agreed harvest strategies or harvest control rules. Nevertheless, if we project the currently potentially successful trajectories forward under different constant catch scenarios we will be able to determine what levels of catch will lead to the stock increasing (on average) and what will lead to decreases (on average). Depending on where the assessment indicates a species is laying will then determine what management action to advise to manage a stock in a desired direction (whether that be to increase or decrease its current state).

output <- doconstC(answer$R1,projn=10,constCatch=250,lastyear=2017,limit=0.2,
                   target=0.48)

Figure 7. A plot of the successful depletion trajectories followed by five years of projection. The limit and target reference points are depicted by the two fine red lines. The green line denotes the end of the final year in which data are available, with projections to the right. The mean depletion varies away from the median because an array of trajectories go extinct on projection

The projections suggest that a total catch of about 250 t would be expected to lead, on average, a small increase in stock levels over the five years. By exploring the outcomes with smaller and larger constant catches the implications can be made clear and appropriate decisions about catch levels or triggers could then be made in an attempt to meet whatever objective is desired for the fishery. What seems to be implied by Figure 7 is that the stock in question was over-fished between about 2003 - 2010 but as the plot of catches vs MSY from the plotcMSY6 analysis illustrate, after 2003 catches were reduced significantly and that led to the stock rebuilding until in 2017 it was estimated (with great uncertainty) to be about 30%B0. In 2017, given the stock is rebuilding and is already above the limit reference point it can be claimed that it is sustainable and that over-fishing is not occurring.

Robustness of MSY Estimate

While it is true that the ranges of the successful r and K combinations is quite wide, because there is a strong negative correlation between these two parameters the spread of the MSY estimates is not so large in relative terms. This can be illustrated by using the trendMSY function which estimates the MSY implied by a series of slices down the K axis from the scatter of r - K combinations. Given that MSY is estimated from the r - K combinations it is possible to map the estimated MSY onto the plot of r - K combinations, along with the estimated trend from the successful r - K combinations (Figure 6).

r <- summcMSY$r
K <- summcMSY$K
meanmsy <- trendMSY(summcMSY$r,summcMSY$K,inc=300)
centr <- central(r)
centK <- central(K)
means <- central(summcMSY$msy)
avMSY <- means["Geometric","Mean"]
# plotprep(width=7,height=4.0)
par(mfrow=c(1,1),mai=c(0.45,0.45,0.05,0.05),oma=c(0.0,0,0.0,0.0)) 
par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)
plot(r,K,type="p",pch=16,col=rgb(1,0,0,1/3),panel.first=grid())
rval <- seq(0.1,0.6,length=100)
kval <- (4 * avMSY)/rval  # calculate the K value that would generate avMSY
lines(rval,kval,col=4,lwd=2)
pickvalid <- which(meanmsy[,"N"] > 0)
lines(meanmsy[pickvalid,"rcenter"],meanmsy[pickvalid,"Kcenter"],lwd=2,col=3)
points(centr[2,1],centK[2,1],pch=16,cex=2,col=1)

Figure 8. A plot of the successful r - K combinations with the estimated mean MSY over-plotted in blue and the mean MSY for horizontal slices of the cloud of successful r - K combinations derived from trendMSY. Note that the overall mean MSY (the black dot) coincides approximately with the point at which the blue and the green lines begin to diverge.

One can imagine contours of MSY running across the r - K parameter space and the catch-MSY analysis is discovering which r - K combinations are consistent with the possibilities (Figure 7).

r <- summcMSY$r
K <- summcMSY$K
avMSY <- seq(100,600,100)  
nMSY <- length(avMSY)
# plotprep(width=7,height=4.0)
par(mfrow=c(1,1),mai=c(0.45,0.45,0.05,0.05),oma=c(0.0,0,0.0,0.0)) 
par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)
plot(r,K,type="p",pch=16,col=rgb(1,0,0,1/3),panel.first=grid(),ylim=c(2500,9000))
rval <- seq(0.1,0.6,length=100)
for (i in 1:nMSY) {
   kval <- (4 * avMSY[i])/rval  # calculate the K value that would generate avMSY
   lines(rval,kval,col=4,lwd=2)
}
x <- seq(0.15,0.375,length=nMSY)
yval <- 4 * (avMSY/x)
for (i in 1:nMSY) text(x[i],yval[i],avMSY[i],cex=1,font=7,pos=4)

Figure 9. A plot of the successful r - K combinations with an array of implied MSY values as implied contours across the r - K parameter space. Some contours would under-estimate the implied productivity while others would over-estimate what was consistent with available observations.

Significant Unknown Recreational Catches

It is often the case that obtaining accurate estimates of recreational catch is rarely done even when it is known that they can rival commercial catches. In States where there are significant recreational catches of commercial species there remains a need to assess the stock status of such species and likely also provide management advice as to what would constitute sustainable commercial catches. As long as there is some notion of the proportion of recreational catches it is still possible to apply the catch-MSY assessment method to a time-series of commercial catches. To demonstrate the effect of having significant amounts of recreational catch on top of a commercial catch we can conduct an analysis using one of the data sets internal to datalowSA, then that can be compared to the same analysis after the catches have been proportionally reduced in some manner as if the total catch was made up of known commercial catches and unknown recreational catches. In this preliminary exploration we will examine a deterministic reduction keeping just 80%, and 60% of known commercial catches while repeating the analysis. To include some variation we can reduce the original catches by a randomly varying proportion around some selected mean, finally we can conduct the analysis on a deterministically reducing commercial proportion down from say 80% to 50% over the time period of data available.

We can use the invert data set. First we obtain the expected outcome by analyzing the full data set.

library(datalowSA)
data(invert)
fish <- invert$fish
glb <- invert$glb

reps <- 5000  # one would run at least 20000, preferably more
answer <- run_cMSY(fish,glb,n=reps,sigpR=0.025)
summcMSY <- summarycMSY(answer,fish,final=TRUE)
ans <- pulloutStats(answer$R1)
round(ans$output[,1:3],3)

So the mean current depletion is estimated to be about 0.31, with an MSY of about 370 tonnes. Now we can answer what would happen if the catch data was only 80% of the original (with variation about that 80% each year; this implies there would be unknown recreational catches of about 20% plus or minus some variable amount each year) while the cpue data remains unchanged. The deterministic reductions gave very precise results so this exposition begins with the somewhat more realistic variation around a given average proportion of recreational catches

commprop <- 0.8
propcom <- rnorm(31,mean=commprop,sd=0.025) # 31 equals the number of years
fishC <- fish
fishC[,"catch"] <- fishC[,"catch"]*propcom
answerC <- run_cMSY(fishC,glb,n=5000,sigpR=0.025)
ansC <- pulloutStats(answerC$R1)
#str(ans10)
round(ansC$output[,1:3],3)
round(ans$output[,1:3],3)
msy <- ansC$output["MSY","Mean"]/ans$output["MSY","Mean"]
depl <- ansC$output["CurrDepl","Mean"]/ans$output["CurrDepl","Mean"]
cat("reduced catch MSY/full catch MSY = ", msy, mean(propcom),"\n")
cat("Proportion of Current Depetion = ",depl,"\n")

You could try different values for the commprop value to see the effect of the recreational proportion increasing. But it appears that the impact of a 20% decrease in known catches is to reduce the sustainable production available to commercial operators by about 20% (and similarly a 40% reduction leads to about a 40% reduction in production).

If instead of the mean recreational proportion staying approximately stable through time it is also possible to examine the effect of a trend in the proportion of commercial catches. Here we inspect the effect of the commercial proportion changing from 80% to 50% (unknown recreational catches change from 20% to 50%) over a 31 year period. Normally one would use at least 20000 replicates rather than 5000.

propcom <- seq(0.8,0.5,length=31) # 31 equals the number of years
fishC <- fish
fishC[,"catch"] <- fishC[,"catch"]*propcom
answerC <- run_cMSY(fishC,glb,n=5000,sigpR=0.025)
ansC <- pulloutStats(answerC$R1)
round(ansC$output[,1:3],3)
round(ans$output[,1:3],3)
msy <- ansC$output["MSY","Mean"]/ans$output["MSY","Mean"]
depl <- ansC$output["CurrDepl","Mean"]/ans$output["CurrDepl","Mean"]
cat("reduced catch MSY/full catch MSY = ", msy, mean(propcom),"\n")
cat("Proportion of Current Depletion = ",depl,"\n")

If this code chunk is run numerous times the outcomes suggest that the estimate of sustainable commercial catch may be biased slightly high but, on the other hand the depletion level is invariably very close to that obtained when using all the catch data.

These two sets of scenarios suggest that if we consider an unknown recreational catch to reflect one of these two scenarios (randomly varying around some relatively constant proportion of the commercial catch or approximately following some trajectory from one proportion to a higher proportion of recreational catch) then it should still be possible to get an estimate of the stock's current depletion level while, at the same time generating potential management advice concerning the commercial fishery. These explorations need to be more fully investigated but, at least with the Catch-MSY approach, the conclusions appear to be acceptable. It can even provide a phase plot of predicted biomass against harvest rate, but obviously only for the commercial component. These results are consistent with those of Rudd and Branch (2017) who find very similar results using fitted stock assessments in a management strategy evaluation framework rather than catch-only methods.

cMSYphaseplot(answerC,fishC)

which can then be compared with the analysis that derived from the full data-set of catches:

cMSYphaseplot(answer,fish)

There are clear differences in that the analysis of partial catch suggests that there was a time that the stock went below the limit reference point but overall, the trajectory obtained is approximately the same as that for the full catch data set.

It is not surprising that differences occur between those analyses that capture all catches and those that include a proportion of unknown catches because the productivity now appears to be different. Nevertheless, the same general pattern observed with all data remains, even with the truncated data. In reality, with a fishery having a significant recreational proportion, we would only have the 'truncated' data but these explorations indicate that it should still be possible to obtain both a status and management advice (at least for the commercial proportion of the fishery) under the assumption that the relative proportion would not change further. Of course, if there is a trend in the proportion of recreational take relative to the commercial take, then only advice in the short term could be deemed useful.

Plausible Levels of Harvest Rate

Even though the overall productivity of the diverse Australian fisheries may be relatively high, the productivity of many individual Australian fisheries tends not to be as great as equivalent fisheries in the northern hemisphere. Many commercial Australian fish species seem to live longer than northern hemisphere equivalents and live less densely on the sea bed, and hence naturally cannot be as productive. It seems unlikely, therefore, that the fishing mortalities observed in places such as the North Sea, where F values of 1.0 (equivalent to harvest rates of between 0.6 and 0.7 each year) have been common, are not plausible and certainly not sustainable here. The catch-MSY method only uses the available time-series of catches and its random selection of r-K combinations can select highly productive combinations which enable the stock to grow well beyond plausible bounds (rejected), there are other combinations whose implied productivity (from the Schaefer model used within the catch-MSY to reflect the stock dynamics) is far too low so that any trajectories go extinct (rejected). Then there are the combinations that imply a productivity that allows for the catches taken but stays within the other constraints implied by the resilience and starting catch levels for the species (accepted trajectories). However, within those accepted trajectories there will be some having implied productivities that only barely survive within the constraints. It is also implied that some of these trajectories reflect relatively high annual harvest rates. For example, using the invert data set we can obtain some trajectories with harvest rates greater than 0.5 (50% of exploitable biomass taken every year).

# library(datalowSA)
data(invert)
fish <- invert$fish
glb <- invert$glb 
# normally one would run 20000+ replicates, but for speed here we use 5000
answer <- run_cMSY(fish,glb,n=5000,sigpR=0.025,finaldepl=c(0.05,0.5),maximumH=1.0)
ans <- pulloutStats(answer$R1)
round(ans$output[,1:3],3)
out <- plottrajectory(answer$R1,fish$year,fish$catch,answer$parbound,
                      oneplot=TRUE,scalar=1.0,plotout=TRUE,plotall=7)

Figure 10. A plot of the implied successful biomass trajectories and their implied annual harvest rates generated by the function plottrajectory with the option oneplot set to TRUE. The red lines are the median values, the green line is at a harvest rate of 0.5, for reference.

In Figure 10. some of the harvest rates in the final years approach 0.5 but between 2000 and about 2007 some trajectories fully breach the 0.5 line and often stay close to it. Part of this reflects the default final depletion range of c(0.05,0.5) given if catches in the final year are less than half the maximum catch (a selection which ignores the fact that management intervention or marketing issues may have controlled catches rather than an inability to catch). However, even if we change the lower final depletion value to, for example, 0.15, while that reduces most of the final year high harvest rates that survive, the period from 2000 still appears to be exceptional.

In order to conduct a sensitivity on the implications of there being an upper limit of harvest rate a new maximumH parameter has been added to the run_cMSY function. Its default value is set to 1.0, so that all harvest rates are possible (obviously one could not take more than 100% of what is present!). If we alter that value to 0.5 the outcome changes.

# library(datalowSA)
# normally one would run 20000+ replicates, but for speed we use 5000
answer <- run_cMSY(fish,glb,n=5000,sigpR=0.025,finaldepl=c(0.05,0.5),maximumH=0.5)
ans <- pulloutStats(answer$R1)
round(ans$output[,1:3],3)
out <- plottrajectory(answer$R1,fish$year,fish$catch,answer$parbound,
                      oneplot=TRUE,scalar=1.0,plotout=TRUE,plotall=7)

Figure 11. A plot of the implied successful biomass trajectories and their implied annual harvest rates generated by the function plottrajectory with the option oneplot set to TRUE and the maximumH parameter set to 0.5. The red lines are the median values, the green line is at a harvest rate of 0.5, for reference.

With this extra constraint the successful biomass trajectories are much more restricted and the harvest rates far lower. This change only improves the depletion level by about 3%, but the productivity, as measured by the MSY is reduced from about 370 t down to 333 t, a reduction of 10%. It may appear counter-intuitive that removing some of the lower trajectories (note the increase in the lower 95th percentile) led to a lower productivity, but the important parts are reflected in the mean values of r, the population growth rate. By comparing the outputs related to Figures 10 and 11 the reduction in r should be clear. It is this reduction that has led to the decrease in productivity.

Such a sensitivity on this parameter (maximumH) is very dependent upon local knowledge of the history of any fishery. Productivity is also partly determined by the resilience attributed to a species. One could also run sensitivities on what resilience was given to a species. The importance of that is that the resilience determines the implied bounds on r used in the search for successful trajectories. Of course, it is now possible to modify these if it is felt necessary by directly entering values for start_r as a vector of two numbers.

Other Sensitivities

We have seen it is possible to examine the implications of changing the initdepl and finaldepl, along with changing the start_r and start_K. A further sensitivity that ought to be conducted relates to the process error term sigpR. This term aims to attempt to capture un-accounted for variation in the stock dynamics between the years. The simple model is so simple it fails to capture many sources of natural variation. It is worthwhile examining alternative values for sigpR. In particular, it is worthwhile effectively turning this source of variation off to see the effect of almost deterministic dynamics. This can be achieved by setting sigpR = 0.

# library(datalowSA)
# normally one would run 20000+ replicates, but for speed we use 5000
answer <- run_cMSY(fish,glb,n=5000,sigpR=0,finaldepl=c(0.05,0.5),maximumH=1.0)
ans <- pulloutStats(answer$R1)
round(ans$output[,1:3],3)
out <- plotconstC(ans$deplet,endyear=2017,constC=0)

Figure 12. A plot of the implied successful depletion trajectories when the sigpR parameter is set to 0. The 23 different starting depletion levels between the default 0.15 - 0.7 can now be clearly seen. The lowest values are the most sparse in terms of successful trajectories.

The assumption that a final depletion level of 5% (0.05) is plausible should be questioned for many Australian fisheries when the catches remain about 1/3 of the maximum and the CPUE barely changes. As always, such data-poor stock assessments should only form part of a weight-of-evidence argument supporting a claim for a set of management advice and stock status determinations.

plotfish(fish,glb)

Figure 13. A plot of the reported catch history and the standardized CPUE for the example data set used in the vignette.

References

Dick, E.J. and A.D. MacCall (2011) Depletion-based stock reduction analysis: a catch-based method for determining sustainable yields for data-poor fish stocks. Fisheries Research 110(2): 331-341

Haddon, M. (2014) Tier 4 analyses in the SESSF, including deep water species. Data from 1986 – 2012. Pp 352 – 461 in Tuck, G.N. (ed) (2014) Stock Assessment for the Southern and Eastern Scalefish and Shark Fishery 2013. Part 2. Australian Fisheries Management Authority and CSIRO Marine and Atmospheric Research, Hobart. 313p.

Haddon, M., Klaer, N., Wayte, S., and G. Tuck (2015) Options for Tier 5 approaches in the SESSF and identification of when data support for harvest strategies are inappro-priate. CSIRO. FRDC Final Report 2013/200. Hobart. 115p.

Kimura, D.K. and J.V. Tagart (1982) Stock Reduction Analysis, another solution to the catch equations. Canadian Journal of Fisheries and Aquatic Sciences 39: 1467 - 1472.

Kimura, D.K., Balsiger, J.W., and Ito, D.H. 1984. Generalized stock reduction analysis. Canadian Journal of Fisheries and Aquatic Sciences 41: 1325–1333.

Little, L.R., Wayte, S.E., Tuck, G.N., Smith, A.D.M., Klaer, N., Haddon, M., Punt, A.E., Thomson, R., Day, J. and M. Fuller (2011) Development and evaluation of a cpue-based harvest control rule for the southern and eastern scalefish and shark fishery of Australia. ICES Journal of Marine Science 68(8): 1699-1705.

Martell, S. and R. Froese (2013) A simple method for estimating MSY from catch and resilience. Fish and Fisheries 14: 504-514

Punt, A.E., Butterworth, D.S. and A.J. Penney (1995) Stock assessment and risk analysis for the South Atlantic population of albacore Thunnus alalunga using an age-structured production model South African Journal of Marine Science 16: 287-310. http://dx.doi.org/10.2989/025776195784156476

R Core Team (2017). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. see also https://cran.r-project.org/

RStudio (2016) www.rstudio.com

Rudd, M.B. and T.A. Branch (2017) Does unreported catch lead to over-fishing? Fish and Fisheries 18: 313-323.

Schaefer, M.B. (1954) Some aspects of the dynamics of populations important to the management of the commercial marine fisheries. Bulletin, Inter-American Tropical Tuna Commission, 1: 25-56.

Schaefer, M.B. (1957) A study of the dynamics of the fishery for yellowfin tuna in the Eastern Tropical Pacific Ocean. Bulletin, Inter-American Tropical Tuna Commission, 2: 247-285

Walters, C.J., Martell, S.J.D. and J. Korman (2006) A stochastic approach to stock reduction analysis. Canadian Journal of Fisheries and Aquatic Sciences 63: 212 - 223.

haddonm/datalowSA documentation built on Nov. 5, 2023, 6:40 p.m.