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

```
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))
```

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.

dataframe containing the initial values for each explored parameter.

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

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)

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

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

.

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.

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.

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

```
# NOT RUN {
# }
# 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)
# }
```

Run the code above in your browser using DataLab