run_cMSY: run_cMSY - sets up and runs the catch-MSY simulations

View source: R/cMSY_funs.R

run_cMSYR Documentation

run_cMSY - sets up and runs the catch-MSY simulations

Description

run_cMSY - sets up and runs the catch-MSY simulations The input data is merely a matrix with, as a minimum, a column of years, and a column of catches, the input object 'glob' needs to contain the resilience. Strictly it only needs the resilience. run_cMSY calls sraMSY, which, in turn calls oneSRA. Eventially it calls plottrajectory to plot out results. This code derives from example code produced by Martell and Froese (2013) A simple method for estimating MSY from catch and resilience Fish and Fisheries 14:504-514.

Usage

run_cMSY(
  indat,
  glob,
  n = 10000,
  incB = 0.025,
  sigpR = 0.025,
  multK = 1,
  finaldepl = NA,
  start_k = NA,
  start_r = NA,
  initdepl = NA,
  maximumH = 1,
  Hyear = NA
)

Arguments

indat

- a data.frame, with at least a 'catch' column containing catch at time t, a 'year' column for year. The 'fish' object from the standard data file or readdata.

glob

the globals list 'glb', from readdata or from one of the included data sets. It contains at least a 'resilience' object for resilience, which is either 'verylow', 'low', 'medium", or 'high', and finally a 'spsname' object, which, not surprisingly, is the name of the species concerned.

n

- defaults to 10000; the number of random selections of r and K. Defines the number of replicate searches within the parameter space.

incB

the increments between the bounds of the initial biomass depletion levels; defaults to 0.025, but have used 0.05 previously

sigpR

the measure of process error added to the dynamics. If set to a very small value, 1e-06, the model will act as deterministic.

multK

a multiplier for K to allow for stock biomasses to rise above K, defaults to 1.0, ie. K is the upper limit.

finaldepl

this allows the option of externally setting the final depletion where there have been major reductions in catch that have not been due to a reduction in the stock; defaults to NA, which sets the finaldepl to the pre-defined priors

start_k

allows an option to alter the starting K values; for example in Orange Roughy, gigantic initial catches possibly make up a significant proportion of the initial biomass so multiplying by 60 or 100 will lead to ridiculous initial K values. A vector of two numbers is required. The default is NA, which means it will be c(maxcatch,60*maxcatch)

start_r

allows for altering the default starting r values; for example in a species though to be of resilience verylow there may be uncertainty over how low and one might want a range from say 0.01 - 0.3 instead of 0.015 - 0.125. The default is NA, which implies the schedule of values in to code will be used. A vector of two numbers is required.

initdepl

this allows the option of externally setting the initial depletion. This may be useful where there is evidence that the stock really has been unfished and is expected to be much closer to 100 than the default setting of c(0.5, 0.975); defaults to NA, which sets the initdepl to the pre-defined priors

maximumH

upper limit of harvest rate included in the constraints; defaults to 1.0, which implies no upper limit.

Hyear

is the index to the year in which a rnage of harvest rates is to be used to constrina the acceptable trajectories.

Value

plots up nine graphs summarizing catches, r, K, and MSY. returns a list containing the vector of relative counts of successful trajectories for different combinations of r and K.

Examples

## Not run: 
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.04)
summcMSY <- summarycMSY(answer,fish,final=TRUE)
str(summcMSY,max.level=1)
out <- plottrajectory(answer$R1,fish$year,fish$catch,answer$parbound,     
                      oneplot=FALSE,scalar=1.0,plotout=TRUE,plotall=7)
plotcMSY6(summcMSY,fish[,"catch"],label=glb$spsname)
ans <- pulloutStats(answer$R1,probabs=c(0.025,0.05,0.5,0.95,0.975))
out <- plotconstC(ans$deplet,endyear=2017,constC=0,console=FALSE,intensity=NA)
outC <- doconstC(answer$R1,constCatch=50,lastyear=2017,console=FALSE,intensity=NA)

## End(Not run)  

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