catchmsy {fishmethods}

Estimating MSY from catch and resilience

Description

This function estimates MSY following Martell and Froese(2012).

Usage

catchmsy(year = NULL, catch = NULL, catchCV = NULL, 
catargs = list(dist = "none", low = 0, up = Inf, unit = "MT"), 
l0 = list(low = 0, up = 1, step = 0), lt = list(low = 0, up = 1, 
refyr = NULL), 
sigv = 0, k = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0), 
r = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0), 
M = list(dist = "unif", low = 0.2, up = 0.2, mean = 0, sd = 0), 
nsims = 10000, catchout = 0, grout = 1, 
graphs = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11), 
grargs = list(lwd = 1, pch = 16, cex = 1, nclasses = 20, mains = " ", 
cex.main = 1, 
cex.axis = 1, cex.lab = 1), 
pstats = list(ol = 1, mlty = 1, mlwd = 1.5, llty = 3, llwd = 1, ulty = 3, 
ulwd = 1), 
grtif = list(zoom = 4, width = 11, height = 13, pointsize = 10))

Arguments

year

vector containing the time series of numeric year labels.

catch

vector containing the time series of catch data (in weight). Missing values are not allowed.

catchCV

vector containing the time series of coefficients of variation associated with catch if resampling of catch is desired; otherwise, catchCV = NULL.

catargs

list arguments associated with resampling of catch. dist is the specification of the resampling distribution to use ("none" = no resampling, "unif"=uniform, "norm" = normal, and "lnorm" =log-normal). If "lnorm" is selected, catch is log transformed and standard deviation on the log scale is calculated from the specified CVs using the relationship sdlog=sqrt(log(CV^2+1)). low and up are the lower and upper limit of distribution (if truncation is desired). unit is the weight unit of catch (used in graph labels; default="MT"). If "unif", the catch must be incorporated in low and up arguments. For instance, if the lower limit to sample is the value of catch, specify low=catch. If some maximum above catch will be the upper limit, specify up=50*catch. The limits for each year will be applied to catch internally.

l0

list arguments for the relative biomass in year 1. low and up are the lower and upper bounds of the starting value of relative biomass (in relation to k) in year 1. step is the step increment to examine. If step=0, then l0 is randomly selected from a uniform distribution using the lower and upper starting values. If step>0, then step increments are used (in this case, the number of simulations (nsims) are used for each increment).

lt

list arguments for the depletion level in the selected reference year (refyr). low and up are the lower and upper bounds of depletion level in refyr. refyr can range from the first year to the year after the last year of catch (t+1).

sigv

standard deviation of the log-normal random process error. signv = 0 for no process error.

k

list arguments for the carrying capacity. dist is the statistical distribution name from which to sample k. low and up are the lower and upper bounds of k in the selected distribution. mean and sd are the mean and standard deviation for selected distributions. The following are valid distributions: "none", "unif" - uniform, "norm" - normal, "lnorm" - log-normal, "gamma" - gamma, and "beta" - beta distributions. "unif" requires non-missing values for low and up. "norm", "lnorm", "gamma" and "beta", require non-missing values for low,up, mean and sd. If "lnorm" is used, mean and sd must be on the natural log scale (keep low and up on the original scale). If dist = "none", the mean is used as a fixed value.

r

list arguments for the intrinsic growth rate. dist is the statistical distribution name from which to sample r. low and up are the lower and upper bounds of r in the selected distribution. mean and sd are the mean and standard deviation for selected distributions. Valid distributions are the same as in k. If dist = "none", the mean is used as a fixed value.

M

list arguments for natural mortality. dist is the statistical distribution name from which to sample M. low and up are the lower and upper bounds of M in the selected distribution. mean and sd are the mean and standard deviation for selected distributions. Valid distributions are the same as in k. M is used to determine exploitation rate (Umsy) at MSY. If dist = "none", the mean is used as a fixed value.

nsims

number of Monte Carlos samples.

catchout

If resampling catch, save catch trajectories from each Monte Carlos simulation - 0 = No (default), 1 = Yes.

grout

numeric argument specifying whether graphs should be printed to console only (1) or to both the console and TIF graph files (2).Use setwd before running function to direct .tif files to a specific directory. Each name of each file is automatically determined.

graphs

vector specifying which graphs should be produced. 1 = line plot of observed catch versus year,2 = point plot of plausible k versus r values, 3 = histogram of plausible r values, 4 = histogram of plausible k values, 5 = histogram of M values, 6 = histogram of MSY from plausible values of l0,k,r, and Bmsy/k, 7 = histogram of Bmsy from plausible values of l0,k,r, and Bmsy/k, 8 = histogram of Fmsy from plausible values of l0,k,r, and Bmsy/k, 9 = histogram of Umsy values from Fmsy and M, 10 = histogram of overfishing limit (OFL) in last year+1 values from Umsys, and 11 = line plots of accepted and rejected biomass trajectores with median and 2.5th and 97.5th percentiles (in red). Any combinations of graphs can be selected within c(). Default is all.

grargs

list control arguments for plotting functions. lwd is the line width for graph = 1 and 11, pch and cex are the symbol character and character expansion value used in graph = 2, nclasses is the nclass argument for the histogram plots (graphs 3-11), mains and cex.main are the titles and character expansion values for the graphs, cex.axis is the character expansion value(s) for the x and y-axis tick labels and cex.lab is the character expansion value(s) for the x and y-axis labels. Single values of nclasses,mains, cex.main,cex.axis, cex.lab are applied to all graphs. To change arguments for specific graphs, enclose arguments within c() in order of the number specified in graphs.

pstats

list control arguments for plotting the mean and 95 and management quantities on respective graphs. ol = 0, do not overlay values on plots, 1 = overlay values on plots. mlty and mlwd are the line type and line width of the mean value; llty and llwd are the line type and line wdith of the 2.5 ulwd are the line type and line width of the 97.5

grtif

list arguments for the .TIF graph files. See tiff help file in R.

Details

The method of Martell and Froese (2012) is used to produce estimates of MSY where only catch and information on resilience is known.

The Schaefer production model is

B[t+1]<-B[t]+r*B[t]*(1-B[t]/k)-catch[t]

where B is biomass in year t, r is the instrince rate of increase, k is the carrying capacity and catch is the catch in year t. Biomass is the first year is calculated by B[1]=k*l0. For sigv>0, the production equation is multiplied by exp(rnorm(1,0,sigv)) if process error is desired. The maximum sustainable yield (MSY) is calculated as

MSY=r*k/4

Biomass at MSY is calculated as

Bmsy=k/2

Fishing mortality at MSY is calculated as

Fmsy=r/2

Exploitation rate at MSY is calculated as

Umsy=(Fmsy/(Fmsy+M))*(1-exp(-Fmsy-M))

The overfishing limit in last year+1 is calculated as

OFL=B[last year +1]*Umsy

length(year)+1 biomass estimates are made for each run.

If using the R Gui (not Rstudio), run

graphics.off() windows(width=10, height=12,record=TRUE) .SavedPlots <- NULL

before running the catchmsy function to recall plots.

Value

Initial

dataframe containing the initial values for each explored parameter.

Parameters

dataframe containing the mean, median, 2.5th and 97.5 plausible (likelihood=1) parameters.

Estimates

dataframe containing the mean, median, 2.5th and 97.5 of the management quantities (i.e., MSY, Bmsy, etc.) for the plausible parameters (likelihood=1)

Values

dataframe containing the values of l0, k, r, Bmsy/k, M and associated management quantities for all (likelihood=0 and likelihood=1) random draws.

end1yr

value of the last year of catch data + 1 for use in function dlproj.

type

designates the output object as a catchmsy object for use in function dlproj.

The biomass estimates from each simulation are not stored in memory but are automatically saved to a .csv file named "Biotraj-cmsy.csv". Yearly values for each simulation are stored across columns. The first column holds the likelihood values for each simulation (1= accepted, 0 = rejected). The number of rows equals the number of simulations (nsims). This file is loaded to plot graph 11 and it must be present in the default or setwd() directory.

When catchout=1, catch values randomly selected are saved to a .csv file named "Catchtraj-cmsy.csv". Yearly values for each simulation are stored across columns. The first column holds the likelihood values (1= accepted, 0 = rejected). The number of rows equals the number of simulations (nsims).

Use setwd() before running the function to change the directory where .csv files are stored.

Note

The random distribution function was adapted from Nadarajah, S and S. Kotz. 2006. R programs for computing truncated distributions. Journal of Statistical Software 16, code snippet 2.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [e-mail]

References

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

See Also

dbsra dlproj

Examples

  ## Not run: 
   data(lingcod)
   outpt<-catchmsy(year=lingcod$year,
     catch=lingcod$catch,catchCV=NULL,
     catargs=list(dist="none",low=0,up=Inf,unit="MT"),
    l0=list(low=0.8,up=0.8,step=0),
    lt=list(low=0.01,up=0.25,refyr=2002),sigv=0,
    k=list(dist="unif",low=4333,up=433300,mean=0,sd=0),
    r=list(dist="unif",low=0.015,up=0.1,mean=0,sd=0),
    M=list(dist="unif",low=0.18,up=0.18,mean=0.00,sd=0.00),
    nsims=30000)
 
## End(Not run)