| Title: | Reference point computation for advice rules |
|---|---|
| Description: | Blah |
| Authors: | Henning Winker [aut, cre] |
| Maintainer: | Henning Winker <[email protected]> |
| License: | EUPL |
| Version: | 1.11.13 |
| Built: | 2026-06-05 14:44:19 UTC |
| Source: | https://github.com/flr/FLRef |
ABItgt() Computes ABI for target F, e.g. ABImsy (Griffith et al. 2023)
ABItgt(stock, ftgt = NULL, thresh = 0.9, ...)ABItgt(stock, ftgt = NULL, thresh = 0.9, ...)
stock |
object of class FLStock or FLStockR |
ftgt |
target F at equilibrium, e.g. Fmsy |
thresh |
quantile ageref treshold, default 0.9 |
*FLQuant*
data(ple4) ABImsy = ABItgt(ple4,ftgt=0.22,thresh=0.9) plot(ABImsy)+ylim(0,2)+ geom_hline(yintercept = c(0.8,1),col=c(2,1),linetype=c(2,1))+ylab(expression(ABI[MSY]))data(ple4) ABImsy = ABItgt(ple4,ftgt=0.22,thresh=0.9) plot(ABImsy)+ylim(0,2)+ geom_hline(yintercept = c(0.8,1),col=c(2,1),linetype=c(2,1))+ylab(expression(ABI[MSY]))
ALK function
ALK(N_a, iALK)ALK(N_a, iALK)
N_a |
numbers at age sample for single event |
iALK |
from iALK() outout |
FLPar of ALK
generates annual ALK sample with length stratified sampling
alk.sample(lfds, alks, nbin = 20, n.sample = 1)alk.sample(lfds, alks, nbin = 20, n.sample = 1)
lfds |
length frequency *FLQuant* |
alks |
annual ALK proportions at age output form ALKs() *FLPars* |
nbin |
number of samples per length bin |
n.sample |
sample size of lfd |
FLPars of sampled ALK
annual ALK function
ALKs(object, iALK)ALKs(object, iALK)
object |
FLQuant with numbers at age |
iALK |
from iALK() outout |
FLPars of ALK
applyALK function to length to age
applyALK(lfds, alks)applyALK(lfds, alks)
alks |
*FLPars* annual ALKs |
lfd |
*FLQuant* with numbers at length |
FLQuant for numbers at age
applyAR() Function to apply ICES Advice rule for F < MSY Btrigger
applyAR(b, btrigger, fmsy, bmin = 0, fmin = 0.001)applyAR(b, btrigger, fmsy, bmin = 0, fmin = 0.001)
b |
biomass |
btrigger |
if b below btrigger F is reduced |
fmsy |
Fmsy target |
bmin |
Option for de facto fishing closure (e.g. bmin = blim) |
fmin |
F below bmin |
F advice
Converts equilibrium reference-point information from an FLBRP object
into approximate Pella-Tomlinson surplus production model parameters. The
function derives the shape parameter m from the equilibrium ratio
Bmsy / B0, and derives the intrinsic productivity parameter r
from the harvest rate at MSY.
asem2spm( object, quant = c("vb", "ssb"), fmsy = NULL, rel = FALSE, spcurve = FALSE ) asem2spm( object, quant = c("vb", "ssb"), fmsy = NULL, rel = FALSE, spcurve = FALSE )asem2spm( object, quant = c("vb", "ssb"), fmsy = NULL, rel = FALSE, spcurve = FALSE ) asem2spm( object, quant = c("vb", "ssb"), fmsy = NULL, rel = FALSE, spcurve = FALSE )
object |
An |
quant |
Character. Biomass metric used to derive the SPM approximation.
Options are |
fmsy |
Optional fishing mortality at MSY. If |
rel |
Logical. If |
spcurve |
Logical. If |
The biomass quantity used for the approximation can be either vulnerable
biomass, quant = "vb", or spawning-stock biomass, quant = "ssb".
Optionally, the function can also return the equilibrium surplus production
curve implied by the supplied FLBRP object.
The function approximates the age-structured equilibrium yield curve with a Pella-Tomlinson production function:
The shape parameter m is solved numerically from the relationship:
and r is calculated from:
where Hmsy = MSY / Bmsy. For quant = "vb", Bmsy
and B0 are based on vulnerable biomass calculated by
vbbrp(). For quant = "ssb", they are based on
ssb().
This approximation is useful for deriving SPM-compatible prior means for
r and m, or for computing SPM-style process deviations from
age-structured stochastic projections.
prior means for r and m *FLPar*
If spcurve = FALSE, an FLPar object with approximate
surplus production parameters:
rIntrinsic productivity parameter of the Pella-Tomlinson surplus production model.
mPella-Tomlinson shape parameter.
BmsyKEquilibrium biomass ratio Bmsy / B0.
HmsyHarvest rate at MSY, calculated as
MSY / Bmsy.
MSYMaximum sustainable yield.
BmsyBiomass at MSY for the selected biomass metric.
B0Unfished biomass for the selected biomass metric.
FmsyFishing mortality at MSY.
If spcurve = TRUE, a list is returned with:
parsThe FLPar object described above.
curveA data frame with columns Biomass,
SP, and quant. If rel = TRUE, Biomass
is scaled by B0 and SP by MSY.
data(ple4) ple4 <- merge_catch(ple4) sr <- srrTMB(as.FLSR(ple4,model=bevholt),spr0=mean(spr0y(ple4))) brp = brp(FLBRP(ple4,sr)) asem2spm(brp) plot_pf(brp) plot_pf(brp,rel=TRUE) ## Not run: data(ple4) ple4 <- merge_catch(ple4) sr <- srrTMB( as.FLSR(ple4, model = bevholt), spr0 = mean(spr0y(ple4)) ) brp <- brp(FLBRP(ple4, sr)) # Approximate SPM parameters using vulnerable biomass pars_vb <- asem2spm(brp, quant = "vb") pars_vb # Approximate SPM parameters using spawning biomass pars_ssb <- asem2spm(brp, quant = "ssb") pars_ssb # Return the surplus production curve in absolute units sp <- asem2spm(brp, quant = "vb", spcurve = TRUE) head(sp$curve) # Return the surplus production curve in relative units sp_rel <- asem2spm(brp, quant = "vb", spcurve = TRUE, rel = TRUE) head(sp_rel$curve) plot_pf(brp) plot_pf(brp, rel = TRUE) ## End(Not run)data(ple4) ple4 <- merge_catch(ple4) sr <- srrTMB(as.FLSR(ple4,model=bevholt),spr0=mean(spr0y(ple4))) brp = brp(FLBRP(ple4,sr)) asem2spm(brp) plot_pf(brp) plot_pf(brp,rel=TRUE) ## Not run: data(ple4) ple4 <- merge_catch(ple4) sr <- srrTMB( as.FLSR(ple4, model = bevholt), spr0 = mean(spr0y(ple4)) ) brp <- brp(FLBRP(ple4, sr)) # Approximate SPM parameters using vulnerable biomass pars_vb <- asem2spm(brp, quant = "vb") pars_vb # Approximate SPM parameters using spawning biomass pars_ssb <- asem2spm(brp, quant = "ssb") pars_ssb # Return the surplus production curve in absolute units sp <- asem2spm(brp, quant = "vb", spcurve = TRUE) head(sp$curve) # Return the surplus production curve in relative units sp_rel <- asem2spm(brp, quant = "vb", spcurve = TRUE, rel = TRUE) head(sp_rel$curve) plot_pf(brp) plot_pf(brp, rel = TRUE) ## End(Not run)
generates FLIndexBiomass with random observation error from an FLStock
bioidx.sim(object, sel = catch.sel(object), sigma = 0.2, q = 0.001, rho = 0)bioidx.sim(object, sel = catch.sel(object), sigma = 0.2, q = 0.001, rho = 0)
object |
FLStock |
sel |
FLQuant with selectivity.pattern |
sigma |
observation error for log(index) |
q |
catchability coefficient for scaling |
FLIndexBiomass
data(ple4) sel = newselex(catch.sel(ple4),FLPar(S50=1.5,S95=2.1,Smax=4.5,Dcv=1,Dmin=0.1)) ggplot(sel)+geom_line(aes(age,data))+ylab("Selectivity")+xlab("Age") object = propagate(ple4,10) sel = newselex(catch.sel(object),FLPar(S50=2.5,S95=3.2,Smax=3.5,Dcv=0.6,Dmin=0.2)) idx = bioidx.sim(object,sel=sel,q=0.0001) # Checks ggplot([email protected])+geom_line(aes(age,data))+ylab("Selectivity")+xlab("Age") ggplot(idx@index)+geom_line(aes(year,data,col=ac(iter)))+theme(legend.position = "none")+ylab("Index")data(ple4) sel = newselex(catch.sel(ple4),FLPar(S50=1.5,S95=2.1,Smax=4.5,Dcv=1,Dmin=0.1)) ggplot(sel)+geom_line(aes(age,data))+ylab("Selectivity")+xlab("Age") object = propagate(ple4,10) sel = newselex(catch.sel(object),FLPar(S50=2.5,S95=3.2,Smax=3.5,Dcv=0.6,Dmin=0.2)) idx = bioidx.sim(object,sel=sel,q=0.0001) # Checks ggplot(idx@sel.pattern)+geom_line(aes(age,data))+ylab("Selectivity")+xlab("Age") ggplot(idx@index)+geom_line(aes(year,data,col=ac(iter)))+theme(legend.position = "none")+ylab("Index")
Uses a bisection algorithm to find the value of a forecast control quantity
that gives a specified target value for a performance statistic. Typical use
cases include finding the constant fbar that gives a chosen risk level,
such as , over a projection period.
bisect( stock, sr, deviances = rec(stock) %=% 1, metrics, refpts, statistic, years, pyears = years, tune, prob, tol = 0.01, maxit = 15, verbose = TRUE )bisect( stock, sr, deviances = rec(stock) %=% 1, metrics, refpts, statistic, years, pyears = years, tune, prob, tol = 0.01, maxit = 15, verbose = TRUE )
stock |
An |
sr |
A stock-recruitment object or model passed to |
deviances |
Recruitment deviances used in the stochastic forecast.
Defaults to |
metrics |
A named list of metric functions passed to
|
refpts |
Reference points passed to |
statistic |
A named list defining the performance statistic to evaluate,
passed to |
years |
Projection years over which the forecast control is applied. |
pyears |
Years over which the performance statistic is evaluated.
Defaults to |
tune |
A named numeric vector or list giving the forecast quantity to
tune and its lower and upper bounds. For example,
|
prob |
Numeric. Target value of the performance statistic. For example,
|
tol |
Numeric. Absolute tolerance used to decide whether the target
statistic has been reached. Default is |
maxit |
Integer. Maximum number of bisection iterations. Default is
|
verbose |
Logical. If |
The function evaluates the forecast at the lower and upper bounds supplied in
tune. If the requested probability lies between the two outcomes, the
interval is repeatedly bisected until the performance statistic is within
tol of prob, or until maxit iterations are reached.
The bisection search requires that the target objective is bracketed by the
two values supplied in tune. In practice, this means that:
where is the performance statistic produced by forecasting with
control value , and is the target value supplied by
prob. If both bounds give performance statistics on the same side of
the target, the function cannot identify a unique bracketed solution.
The function is most useful for risk-based advice calculations, such as
finding the fishing mortality that gives a pre-specified probability of
falling below Blim over a future evaluation period.
If a solution is found, an FLStock object returned by
fwd() or ffwd() for the tuned control value. If the requested
target is already achieved at the lower or upper bound, the corresponding
forecast object is returned. If the supplied range does not bracket the
target statistic, a list with elements min and max is
returned, containing the forecasts at the lower and upper bounds.
Burden, R. L. and Faires, J. D. (1985). Numerical Analysis. 3rd edition. PWS Publishers. Section 2.1: The Bisection Algorithm.
## Not run: data(ple4) # Extend and propagate stock stock <- propagate(stf(ple4, end = 2118), 100) # Define a Beverton-Holt stock-recruitment relationship srr <- predictModel( model = rec ~ a * ssb * exp(-b * ssb), params = FLPar(a = 5.20, b = 1.65e-6) ) # Generate autocorrelated recruitment deviances devs <- ar1rlnorm( rho = 0.4, years = 2018:2118, iters = 100, meanlog = 0, sdlog = 0.5 ) # Define an ICES-style risk statistic statistic <- list( FP05 = list( ~ yearMeans((SB / SBlim) < 1), name = "P.05", desc = "ICES P.05" ) ) # Find fbar that gives P(SSB < Blim) = 0.05 fp05fwd <- bisect( stock = stock, sr = srr, deviances = devs, metrics = list(SB = ssb), refpts = FLPar(SBlim = 150000), statistic = statistic, years = 2018:2118, pyears = 2069:2118, tune = list(fbar = c(0.1, 1.0)), prob = 0.05 ) ## End(Not run)## Not run: data(ple4) # Extend and propagate stock stock <- propagate(stf(ple4, end = 2118), 100) # Define a Beverton-Holt stock-recruitment relationship srr <- predictModel( model = rec ~ a * ssb * exp(-b * ssb), params = FLPar(a = 5.20, b = 1.65e-6) ) # Generate autocorrelated recruitment deviances devs <- ar1rlnorm( rho = 0.4, years = 2018:2118, iters = 100, meanlog = 0, sdlog = 0.5 ) # Define an ICES-style risk statistic statistic <- list( FP05 = list( ~ yearMeans((SB / SBlim) < 1), name = "P.05", desc = "ICES P.05" ) ) # Find fbar that gives P(SSB < Blim) = 0.05 fp05fwd <- bisect( stock = stock, sr = srr, deviances = devs, metrics = list(SB = ssb), refpts = FLPar(SBlim = 150000), statistic = statistic, years = 2018:2118, pyears = 2069:2118, tune = list(fbar = c(0.1, 1.0)), prob = 0.05 ) ## End(Not run)
function to assign B[y+1] to B[y]. Warning correlation structure of B[y+1] and F[y] is meaningless
blag(mvn, verbose = TRUE)blag(mvn, verbose = TRUE)
mvn |
output list of quant posteriors and mle's
Henning Winker (GFCM)
generates catch.n with lognormal annual and multinomial age composition observation error
ca.sim( object, ess = 200, what = c("catch", "landings", "discards")[1], rescale = FALSE )ca.sim( object, ess = 200, what = c("catch", "landings", "discards")[1], rescale = FALSE )
object |
FLQuant |
ess |
effective sample size for age composition |
what |
c("catch", "landings", "discards") |
rescale |
rescale to the abudance scale of the input |
sel |
FLQuant with selectivity.pattern e.g. catch.sel() |
FLQuant with catch.n samples
data(ple4) object = propagate(catch.sel(ple4),10) ca = ca.sim(object,ess=200) # Checks ggplot(ca)+geom_line(aes(year,data,col=ac(iter)))+facet_wrap(~age)+ theme(legend.position = "none")+ylab("Index")data(ple4) object = propagate(catch.sel(ple4),10) ca = ca.sim(object,ess=200) # Checks ggplot(ca)+geom_line(aes(year,data,col=ac(iter)))+facet_wrap(~age)+ theme(legend.position = "none")+ylab("Index")
computeFbrp() Computes biological reference points corresponding to the proxy Fbrp
computeFbrp( stock, sr = "missing", proxy = NULL, x = NULL, blim = 0.1, type = c("b0", "btgt", "value"), btri = "missing", bpa = "missing", bthresh = "missing", verbose = T, fmax = 10, ... )computeFbrp( stock, sr = "missing", proxy = NULL, x = NULL, blim = 0.1, type = c("b0", "btgt", "value"), btri = "missing", bpa = "missing", bthresh = "missing", verbose = T, fmax = 10, ... )
stock |
object of class FLStock |
sr |
stock recruitment model of class FLSR |
proxy |
choice of Fmsy proxies (combinations permitted)
|
x |
basis in percent for sprx and bx, e.g. 40 for spr40 |
blim |
values < 1 are taken as fraction to B0 and blim > 1 as absolute values unless specified otherwise |
type |
type of blim input, values < 1 are
|
btri |
Btrigger can specified as absolute value |
bpa |
Bpa can specified as absolute value |
bthresh |
Bthresh (GFCM) interchangeable use with Bpa |
verbose |
|
fmax |
maximum Flim = max(Flim,fmax*Fbrp) |
brp object of class FLBRP with computed Fbrp reference points
data(ple4) srr = srrTMB(as.FLSR(ple4,model=rickerSV),spr0=spr0y(ple4)) brp = computeFbrp(stock=ple4,sr=srr,proxy=c("sprx","f0.1"),blim=0.1,type="b0") ploteq(brp,obs=TRUE,refpts="msy")data(ple4) srr = srrTMB(as.FLSR(ple4,model=rickerSV),spr0=spr0y(ple4)) brp = computeFbrp(stock=ple4,sr=srr,proxy=c("sprx","f0.1"),blim=0.1,type="b0") ploteq(brp,obs=TRUE,refpts="msy")
computeFbrps() Computes biological reference points corresponding to the proxy Fbrp
computeFbrps( stock, sr = "missing", proxy = c("sprx", "bx", "all"), fmsy = FALSE, f0.1 = TRUE, fmax = 5, verbose = T, ... )computeFbrps( stock, sr = "missing", proxy = c("sprx", "bx", "all"), fmsy = FALSE, f0.1 = TRUE, fmax = 5, verbose = T, ... )
stock |
object of class FLStock |
sr |
stock recruitment model of class FLSR |
fmsy |
if TRUE, Fmsy is computed (not suggest for segreg or geomean sr) |
f0.1 |
if TRUE, F0.1 is computed |
fmax |
maximum Flim = minfmax*Fbrp) |
verbose |
|
proxies |
choice of Fmsy proxies
|
brp object of class FLBRP with computed Fbrp reference points
Fbrp() Extract Fbrp based reference points from output of computeFbrp
Fbrp(brp)Fbrp(brp)
brp |
input of class FLBRP from ComputeFbrp |
FLPar object with computed Fbrp reference points
Fe40() Patterson estimator for Fmsy
Fe40(stock, nyears = 3)Fe40(stock, nyears = 3)
stock |
input of class FLStock |
nyears |
number of years to average |
value
flr2stars()
flr2stars(object, uncertainty = NULL, quantiles = c(0.05, 0.95))flr2stars(object, uncertainty = NULL, quantiles = c(0.05, 0.95))
object |
of class FLStockR (MLE) |
uncertainty |
of class FLStock with iters |
quantities |
default is 90CIs as c(0.05,0.95) |
STARS list with $timeseris and $refpts
Fmmy() Uses opt.bisect to derive the F at Maximum Median Yield from stochastic simulations
Fmmy( brp, sigmaR = 0.5, rho = 0, nyears = 100, iters = 250, yrs.eval = NULL, range = "missing", tol = 0.001, maxit = 15, verbose = TRUE )Fmmy( brp, sigmaR = 0.5, rho = 0, nyears = 100, iters = 250, yrs.eval = NULL, range = "missing", tol = 0.001, maxit = 15, verbose = TRUE )
brp |
output object from computeFbrp() of class FLBRP |
sigmaR |
lognormal recruitment standard deviation |
rho |
AR1 recruitment autocorrelation coefficient |
nyears |
number of simulation years |
iters |
number simulation iterations |
yrs.eval |
last years to be used evaluation period, default nyears/2 |
range |
range of Fbar value to be evaluated |
tol |
tolerance |
maxit |
number of steps |
verbose |
cat comments |
list of FLPar, FLStock and FLBRP objects
data(ple4) bh = srrTMB(as.FLSR(ple4,model=bevholtSV),spr0=spr0y(ple4)) brp = computeFbrp(ple4,bh,proxy=c("bx","msy"),x=35,blim=0.1) fmmy = Fmmy(brp,sigmaR=0.7,rho=0.3) getF(fmmy) # FMMY value plotFsim(fmmy) brpfmmy = computeFbrp(ple4,bh,proxy=getF(fmmy),blim=0.1) fsim = Fsim(brpfmmy,sigmaR=0.7,rho=0.3) plotFsim(fsim)data(ple4) bh = srrTMB(as.FLSR(ple4,model=bevholtSV),spr0=spr0y(ple4)) brp = computeFbrp(ple4,bh,proxy=c("bx","msy"),x=35,blim=0.1) fmmy = Fmmy(brp,sigmaR=0.7,rho=0.3) getF(fmmy) # FMMY value plotFsim(fmmy) brpfmmy = computeFbrp(ple4,bh,proxy=getF(fmmy),blim=0.1) fsim = Fsim(brpfmmy,sigmaR=0.7,rho=0.3) plotFsim(fsim)
Calculates the Fbar value giving a maximum probability of ssb being below Blim of 5 percent
Fp05( object, iters = "missing", range = "missing", tol = 0.001, maxit = 20, verbose = TRUE )Fp05( object, iters = "missing", range = "missing", tol = 0.001, maxit = 20, verbose = TRUE )
object |
output from Fsim() |
iters |
Number of iterations, cannot exceed input object |
range |
range of Fbar value to be evaluated |
verbose |
Should progress be shown, TRUE. |
list
data(ple4) bh = srrTMB(as.FLSR(ple4,model=bevholtSV),spr0=spr0y(ple4)) brp = computeFbrp(ple4,bh,proxy="bx",x=35,blim=0.2) # set Blim higher fsim = Fsim(brp,sigmaR=0.7,rho=0.3,iters=500) plotFsim(fsim) fp.05 = Fp05(fsim) plotFsim(fp.05,panels=c(2,4)) # black line is Fp0.05 getF(fp.05)data(ple4) bh = srrTMB(as.FLSR(ple4,model=bevholtSV),spr0=spr0y(ple4)) brp = computeFbrp(ple4,bh,proxy="bx",x=35,blim=0.2) # set Blim higher fsim = Fsim(brp,sigmaR=0.7,rho=0.3,iters=500) plotFsim(fsim) fp.05 = Fp05(fsim) plotFsim(fp.05,panels=c(2,4)) # black line is Fp0.05 getF(fp.05)
Fsim() Simulates stochastic stock dynamics under under constant Fbrp
Fsim( brp, Ftgt = NULL, sigmaR = 0.5, rho = 0, nyears = 100, iters = 250, yrs.eval = NULL, verbose = TRUE )Fsim( brp, Ftgt = NULL, sigmaR = 0.5, rho = 0, nyears = 100, iters = 250, yrs.eval = NULL, verbose = TRUE )
brp |
output object from computeFbrp() of class FLBRP |
sigmaR |
lognormal recruitment standard deviation |
rho |
AR1 recruitment autocorrelation coefficient |
nyears |
number of simulation years |
iters |
number simulation iterations |
yrs.eval |
last years to be used evaluation period, default nyears/2 |
verbose |
cat comments |
list of FLPar, FLStock and FLBRP objects
data(ple4) hs = srrTMB(as.FLSR(ple4,model=segreg),spr0=spr0y(ple4),lplim=0.05,uplim=0.25) blim = params(hs)[[2]] brp = computeFbrp(ple4,hs,proxy=c("sprx","f0.1","msy"),x=40,blim=blim) ploteq(brp) fsim = Fsim(brp,sigmaR=0.7,rho=0.3) plotFsim(fsim) plotFsim(fsim,panels=2)data(ple4) hs = srrTMB(as.FLSR(ple4,model=segreg),spr0=spr0y(ple4),lplim=0.05,uplim=0.25) blim = params(hs)[[2]] brp = computeFbrp(ple4,hs,proxy=c("sprx","f0.1","msy"),x=40,blim=blim) ploteq(brp) fsim = Fsim(brp,sigmaR=0.7,rho=0.3) plotFsim(fsim) plotFsim(fsim,panels=2)
Runs stochastic FLR forward projections with ffwd() over a grid of
fishing mortality values. By default, the grid is taken from
fbar(brp), with the first value replaced by Fmsy. For each
F value, a corresponding FLBRP, FLStock, and FLSR
object is constructed, propagated over iterations, and projected forward
using recruitment deviances generated by rlnormar1().
Fsim_grid( brp, Fgrid = NULL, sigmaR = 0.5, rho = 0, nyears = 100, iters = 10, yrs.eval = NULL, parallel = FALSE, workers = NULL, progress = TRUE, verbose = TRUE )Fsim_grid( brp, Fgrid = NULL, sigmaR = 0.5, rho = 0, nyears = 100, iters = 10, yrs.eval = NULL, parallel = FALSE, workers = NULL, progress = TRUE, verbose = TRUE )
brp |
An |
Fgrid |
Optional numeric vector of fishing mortality values. If
|
sigmaR |
Numeric. Lognormal recruitment standard deviation. |
rho |
Numeric. AR(1) autocorrelation coefficient for recruitment deviances. |
nyears |
Integer. Number of projection years. |
iters |
Integer. Number of stochastic iterations. |
yrs.eval |
Optional integer. Number of final years intended for
evaluation. If |
parallel |
Logical. If |
workers |
Optional integer. Number of parallel workers. |
verbose |
Logical. If |
keep.stock |
Logical. Currently unused; retained for compatibility. |
The function is intended for exploring the emergent biomass process error
implied by an age-structured operating model under different constant
fishing mortality levels. The output can be passed to pe_sampler()
to obtain a data frame of depletion and process-error estimates.
The function constructs one FLBRP object per F value by replacing
fbar(brp) with a constant F. These objects are then converted to
FLStock and FLSR objects before stochastic forward projection.
Recruitment deviances are generated once using rlnormar1() and then
passed to each ffwd() call. This allows the F scenarios to be compared
under the same recruitment-deviation structure.
An FLStocks object containing one projected FLStock
per F value in the grid. The stocks are named using the corresponding
F values.
## Not run: data(ple4) ple4 <- merge_catch(ple4) bh <- srrTMB( as.FLSR(ple4, model = bevholt), spr0 = spr0y(ple4), r0 = NULL ) brp <- computeFbrp(ple4, bh, proxy = c("msy")) sim <- Fsim_grid( brp, sigmaR = 0.5, rho = 0.3, nyears = 100, iters = 20, parallel = TRUE, workers= 12 ) class(sim) names(sim) ## End(Not run)## Not run: data(ple4) ple4 <- merge_catch(ple4) bh <- srrTMB( as.FLSR(ple4, model = bevholt), spr0 = spr0y(ple4), r0 = NULL ) brp <- computeFbrp(ple4, bh, proxy = c("msy")) sim <- Fsim_grid( brp, sigmaR = 0.5, rho = 0.3, nyears = 100, iters = 20, parallel = TRUE, workers= 12 ) class(sim) names(sim) ## End(Not run)
generates an up-down-constant F-pattern
fudc( object, fref = 0.2, fhi = 2.5, flo = 0.8, sigmaF = 0.2, breaks = c(0.5, 0.75) )fudc( object, fref = 0.2, fhi = 2.5, flo = 0.8, sigmaF = 0.2, breaks = c(0.5, 0.75) )
object |
An *FLStock* |
fref |
reference denominator for fbar |
fhi |
factor for high F as fhi = fbar/fref |
flo |
factor for low F as flo = fbar/fref |
sigmaF |
variation on fbar |
breaks |
relative location of directional change |
FLQuant
data(ple4) sr <- srrTMB(as.FLSR(ple4,model=bevholtSV),spr0=mean(spr0y(ple4))) brp = computeFbrp(ple4,sr,proxy="msy") fmsy = Fbrp(brp)["Fmsy"] stki = propagate(ple4,100) fy = fudc(ple4,fhi=2,flo=0.9,fref=fmsy,sigmaF=0) fyi = fudc(stki,fhi=2,flo=0.9,fref=fmsy,sigmaF=0.2) plot(fy,fyi)+ylab("F") #Forcasting om <- FLStockR(ffwd(stki,sr,fbar=fyi)) om@refpts = Fbrp(brp) plotAdvice(window(om,start=1960))data(ple4) sr <- srrTMB(as.FLSR(ple4,model=bevholtSV),spr0=mean(spr0y(ple4))) brp = computeFbrp(ple4,sr,proxy="msy") fmsy = Fbrp(brp)["Fmsy"] stki = propagate(ple4,100) fy = fudc(ple4,fhi=2,flo=0.9,fref=fmsy,sigmaF=0) fyi = fudc(stki,fhi=2,flo=0.9,fref=fmsy,sigmaF=0.2) plot(fy,fyi)+ylab("F") #Forcasting om <- FLStockR(ffwd(stki,sr,fbar=fyi)) om@refpts = Fbrp(brp) plotAdvice(window(om,start=1960))
Function to summarise forecast results
Function to summarise ICES fishing opportunities
fwd2ices(stock, uncertainty, eval.yrs = NULL, dB = NULL, refyr = NULL) fwd2ices(stock, uncertainty, eval.yrs = NULL, dB = NULL, refyr = NULL)fwd2ices(stock, uncertainty, eval.yrs = NULL, dB = NULL, refyr = NULL) fwd2ices(stock, uncertainty, eval.yrs = NULL, dB = NULL, refyr = NULL)
uncertainty |
*FLStocks* with list of *FLStockR* objects and iters |
eval.yrs |
evaluation years of forecast |
dB |
computes change in percentage biomass to refyr |
refyr |
change in percentage biomass to refyr |
object |
*FLStocks* with list of *FLStockR* objects |
rel |
if TRUE ratios B/Btgt and F/Ftgt are shown |
stk |
*FLStocks* with list of *FLStockR* objects |
data.frame
data.frame
Function to summarise forecast results
fwd2stars(object, eval.yrs = NULL, rel = NULL, dB = NULL, refyr = NULL)fwd2stars(object, eval.yrs = NULL, rel = NULL, dB = NULL, refyr = NULL)
object |
*FLStocks* with list of *FLStockR* objects |
eval.yrs |
evaluation years of forecast |
rel |
if TRUE ratios B/Btgt and F/Ftgt are shown |
dB |
computes change in percentage biomass to refyr |
data.frame
Function to summarise Type 1 GFCM fishing opportunities
fwd2type1( stock, uncertainty = NULL, eval.yrs = NULL, dB = NULL, refyr = NULL, rel = NULL, osa = FALSE, maxdec = 2 )fwd2type1( stock, uncertainty = NULL, eval.yrs = NULL, dB = NULL, refyr = NULL, rel = NULL, osa = FALSE, maxdec = 2 )
uncertainty |
*FLStocks* with list of *FLStockR* objects and iters |
eval.yrs |
evaluation years of forecast |
refyr |
change in percentage biomass to refyr |
rel |
values of biomass and F are expressed relative to the target |
osa |
one-step-ahead for the biomass evaluation year |
maxdec |
maximum decimal |
stk |
*FLStocks* with list of *FLStockR* objects |
data.frame
fwdF4B () Bisection Function to search F for a given biomass
fwdF4B( stock, sr, btgt, nfy = 3, niy = 1, ival = niy, imet = "TAC", ftune = c(0, 2), tol = 0.001, verbose = TRUE, ssbQ = 1, recQ = 1 )fwdF4B( stock, sr, btgt, nfy = 3, niy = 1, ival = niy, imet = "TAC", ftune = c(0, 2), tol = 0.001, verbose = TRUE, ssbQ = 1, recQ = 1 )
stock |
FLStock object |
sr |
stock recruitment relationship |
btgt |
target biomass of the search |
nfy |
number of forecast (default 3) |
niy |
number of intermediate years (default 1) |
ival |
intermediate year value |
imet |
intermediate year metric ("TAC","F") |
ftune |
tuning limits for F search c(0.1,2) |
tol |
precision tolerance |
verbose |
*FLStock*
Helper functio to extract F from various FLRef output
getF(x)getF(x)
x |
output object from computeFbrp() of class FLBRP |
huecol
huecol(n, alpha = 1)huecol(n, alpha = 1)
n |
number of colors |
alpha |
transluscency |
inverse ALK function with lmin added to FLCore::invALK
iALK( params, model = vonbert, age, cv = 0.1, lmin = 5, lmax = 1.2, bin = 1, max = ceiling(linf * lmax), reflen = NULL )iALK( params, model = vonbert, age, cv = 0.1, lmin = 5, lmax = 1.2, bin = 1, max = ceiling(linf * lmax), reflen = NULL )
params |
growth parameter, default FLPar(linf,k,t0) |
model |
growth model, only option currently vonbert |
age |
age vector |
cv |
of length-at-age |
lmin |
minimum length |
lmax |
maximum upper length specified lmax*linf |
bin |
length bin size, dafault 1 |
max |
maximum size value |
reflen |
evokes fixed sd for L_a at sd = cv*reflen |
timing |
t0 assumed 1st January, default seq(0,11/12,1/12), but can be single event 0.5 |
unit |
default is "cm" |
FLPar age-length matrix
generates FLIndex with lognormal annual and multinomial age composition observation error
idx.sim( object, sel = catch.sel(object), ages = NULL, years = NULL, ess = 200, sigma = 0.2, q = 0.01 )idx.sim( object, sel = catch.sel(object), ages = NULL, years = NULL, ess = 200, sigma = 0.2, q = 0.01 )
object |
FLStock |
sel |
FLQuant with selectivity.pattern |
ages |
define age range |
years |
define year range |
ess |
effective sample size for age composition sample |
sigma |
annual observation error for log(q) |
q |
catchability coefficient for scaling |
FLIndex
data(ple4) sel = newselex(catch.sel(ple4),FLPar(S50=1.5,S95=2.1,Smax=4.5,Dcv=1,Dmin=0.1)) ggplot(sel)+geom_line(aes(age,data))+ylab("Selectivity")+xlab("Age") object = propagate(ple4,10) idx = idx.sim(object,sel=sel,ess=200,sigma=0.2,q=0.01,years=1994:2017) # Checks ggplot([email protected])+geom_line(aes(age,data))+ylab("Selectivity")+xlab("Age") ggplot(idx@index)+geom_line(aes(year,data,col=ac(iter)))+facet_wrap(~age,scales="free_y")+ theme(legend.position = "none")+ylab("Index")data(ple4) sel = newselex(catch.sel(ple4),FLPar(S50=1.5,S95=2.1,Smax=4.5,Dcv=1,Dmin=0.1)) ggplot(sel)+geom_line(aes(age,data))+ylab("Selectivity")+xlab("Age") object = propagate(ple4,10) idx = idx.sim(object,sel=sel,ess=200,sigma=0.2,q=0.01,years=1994:2017) # Checks ggplot(idx@sel.pattern)+geom_line(aes(age,data))+ylab("Selectivity")+xlab("Age") ggplot(idx@index)+geom_line(aes(year,data,col=ac(iter)))+facet_wrap(~age,scales="free_y")+ theme(legend.position = "none")+ylab("Index")
jabba2FLStockR()
jabba2FLStockR()
jabba2FLStockR(jabba, blim = 0.3, bthr = 0.5, thin = 10, rel = FALSE) jabba2FLStockR(jabba, blim = 0.3, bthr = 0.5, thin = 10, rel = FALSE)jabba2FLStockR(jabba, blim = 0.3, bthr = 0.5, thin = 10, rel = FALSE) jabba2FLStockR(jabba, blim = 0.3, bthr = 0.5, thin = 10, rel = FALSE)
jabba |
fit from JABBA fit_jabba() or jabba$kbtrj |
blim |
biomass limit reference point as fraction of Bmsy |
thin |
thinnig rate of retained iters |
rel |
if TRUE ratios BBmsy and FFmsy are stored |
bpa |
biomass precautionary reference point as fraction of Bmsy |
FLStockR with refpts
FLStockR with refpts
jabba2stars()
jabba2stars(jabba, quantiles = c(0.05, 0.95), blim = 0.3, bthr = 0.5)jabba2stars(jabba, quantiles = c(0.05, 0.95), blim = 0.3, bthr = 0.5)
jabba |
fit from JABBA fit_jabba() or jabba$kbtrj |
quantiles |
default is 90CIs as c(0.05,0.95) |
blim |
biomass limit point as fraction of Bmsy, default 0.3Bmsy (ICES) |
bthr |
biomass precautionary point as fraction of Bmsy, default 0.5Bmsy (ICES) |
STARS list with $timeseris and $refpts
function to generate survey (pulse) and continuous LFDs
len.sim( N_a, params, model = vonbert, ess = 250, timing = seq(0, 11/12, 1/12), unit = "cm", scale = TRUE, reflen = NULL, bin = 1, cv = 0.1, lmin = 5, lmax = 1.2 )len.sim( N_a, params, model = vonbert, ess = 250, timing = seq(0, 11/12, 1/12), unit = "cm", scale = TRUE, reflen = NULL, bin = 1, cv = 0.1, lmin = 5, lmax = 1.2 )
N_a |
numbers at age sample |
params |
growth parameter, default FLPar(linf,k,t0) |
model |
growth model, only option currently vonbert |
ess |
effective sample size |
timing |
t0 assumed 1st January, default seq(0,11/12,1/12), but can be single event 0.5 |
unit |
default is "cm" |
scale |
if TRUE scaled to N_a input |
reflen |
evokes fixed sd for L_a at sd = cv*reflen |
bin |
length bin size, dafault 1 |
cv |
variation in L_a |
lmin |
minimum length |
lmax |
maximum upper length specified lmax*linf |
FLQuant for length
function to generate survey (pulse) and continuous LFDs
lfd.sim( object, stock, sel = catch.sel(stock), params, model = vonbert, ess = 250, timing = seq(0, 11/12, 1/12), timeref = 0.5, unit = "cm", scale = TRUE, reflen = NULL, bin = 1, cv = 0.1, lmin = 5, lmax = 1.2 )lfd.sim( object, stock, sel = catch.sel(stock), params, model = vonbert, ess = 250, timing = seq(0, 11/12, 1/12), timeref = 0.5, unit = "cm", scale = TRUE, reflen = NULL, bin = 1, cv = 0.1, lmin = 5, lmax = 1.2 )
object |
*FLQuant* numbers at age sample |
stock |
*FLStock* object |
sel |
selectivity, default catch.sel(stock) |
params |
growth parameter, default FLPar(linf,k,t0) |
model |
growth model, only option currently vonbert |
ess |
effective sample size |
timing |
default constinoues seq(0,11/12,1/12), but can be single event 0.5 |
timeref |
reference timing of the sample, default 0.5 (e.g. survey or catch.n) |
unit |
default is "cm" |
scale |
if TRUE scaled to N_a input |
reflen |
evokes fixed sd for L_a at sd = cv*reflen |
bin |
length bin size, dafault 1 |
cv |
variation in L_a |
lmin |
minimum length |
lmax |
maximum upper length specified lmax*linf |
FLQuant for length
Converts an FLStock object with total catch numbers into a
landings-only stock by assigning catch.n to landings.n and
setting discards.n to zero. Landings and discards weights are then
recomputed.
merge_catch(stock)merge_catch(stock)
stock |
An |
This is a convenience function for cases where catch should be treated as landed catch, for example before forward projection or when preparing simplified single-fleet inputs. The function modifies:
landings.n(stock): set equal to catch.n(stock)
discards.n(stock): set to zero
landings(stock): recomputed with computeLandings()
discards(stock): recomputed with computeDiscards()
An FLStock object where all catch numbers are assigned to
landings and discards are set to zero.
## Not run: stock <- merge_catch(stock) # Check that discards have been removed range(discards.n(stock)) discards(stock) ## End(Not run)## Not run: stock <- merge_catch(stock) # Check that discards have been removed range(discards.n(stock)) discards(stock) ## End(Not run)
computes Lorenzen M with scaling option
Mlorenzen(object, Mref = "missing", Aref = 2)Mlorenzen(object, Mref = "missing", Aref = 2)
object |
weight-at-age of class *FLQuant* |
Mref |
reference M for scaling |
Aref |
reference Age for scaling |
FLQuant m()
data(ple4) Ml = Mlorenzen(stock.wt(ple4)) # Scale Ms = Mlorenzen(stock.wt(ple4),Mref=0.2,Aref=2) flqs = FLQuants(Lorenzen=Ml,Scaled=Ms)data(ple4) Ml = Mlorenzen(stock.wt(ple4)) # Scale Ms = Mlorenzen(stock.wt(ple4),Mref=0.2,Aref=2) flqs = FLQuants(Lorenzen=Ml,Scaled=Ms)
generates flexible 5-paramater selex curves
newselex(object, selexpars)newselex(object, selexpars)
object |
FLQuant from catch.sel() or sel.pattern() |
selexpars |
Selectivity Parameters selexpars S50, S95, Smax, Dcv, Dmin
|
FLquant with selectivity pattern
data(ple4) sel = newselex(catch.sel(ple4),FLPar(S50=2,S95=3,Smax=4.5,Dcv=0.6,Dmin=0.3)) ggplot(sel)+geom_line(aes(age,data))+ylab("Selectivity")+xlab("Age") # Simulate harvest(ple4)[] = sel sr <- srrTMB(as.FLSR(ple4,model=bevholtSV),spr0=mean(spr0y(ple4))) brp = computeFbrp(ple4,sr,proxy="msy") fbar(brp) = FLQuant(rep(0.01,70)) stk = as(brp,"FLStock") units(stk) = standardUnits(stk) its = 100 stk <- FLStockR(propagate(stk, its)) stk@refpts= Fbrp(brp) b0=an(Fbrp(brp)["B0"]) control = FLPar(Feq=0.15,Frate=0.1,Fsigma=0.15,SB0=b0,minyear=2,maxyear=70,its=its) run <- rffwd(stk, sr=sr,control=control,deviances=ar1rlnorm(0.3, 1:70, its, 0, 0.6)) plotAdvice(run)data(ple4) sel = newselex(catch.sel(ple4),FLPar(S50=2,S95=3,Smax=4.5,Dcv=0.6,Dmin=0.3)) ggplot(sel)+geom_line(aes(age,data))+ylab("Selectivity")+xlab("Age") # Simulate harvest(ple4)[] = sel sr <- srrTMB(as.FLSR(ple4,model=bevholtSV),spr0=mean(spr0y(ple4))) brp = computeFbrp(ple4,sr,proxy="msy") fbar(brp) = FLQuant(rep(0.01,70)) stk = as(brp,"FLStock") units(stk) = standardUnits(stk) its = 100 stk <- FLStockR(propagate(stk, its)) stk@refpts= Fbrp(brp) b0=an(Fbrp(brp)["B0"]) control = FLPar(Feq=0.15,Frate=0.1,Fsigma=0.15,SB0=b0,minyear=2,maxyear=70,its=its) run <- rffwd(stk, sr=sr,control=control,deviances=ar1rlnorm(0.3, 1:70, its, 0, 0.6)) plotAdvice(run)
The plain bisection algorithm (Burden & Douglas, 1985) is employed here to find the value of a given forecast target quantity (e.g. 'fbar') for which a selected value of a performance statistic is obtained over a chosen period.
opt.bisect( stock, sr, deviances = rec(stock) %=% 1, metrics, statistic, years, pyears = years, tune, tol = 0.001, maxit = 15, log = TRUE, verbose = TRUE )opt.bisect( stock, sr, deviances = rec(stock) %=% 1, metrics, statistic, years, pyears = years, tune, tol = 0.001, maxit = 15, log = TRUE, verbose = TRUE )
stock |
object class FLStock |
sr |
object class FLSR |
metrics |
FLQuant of FLStock to be defined |
statistic |
|
years |
years to be evaluated |
tune |
range for input x |
tol |
tolerance level |
maxit |
number of optimisation steps |
log |
if TRUE, optimise on log-scale |
Credits to Iago Mosqueira
Burden, Richard L.; Faires, J. Douglas (1985), "2.1 The Bisection Algorithm", Numerical Analysis (3rd ed.), PWS Publishers, ISBN 0-87150-857-5
data(ple4) stock <- propagate(stf(ple4, end=2118), 200) srr <- predictModel(model=rec ~ ifelse(ssb <= b, a * ssb, a * b), params=FLPar(a=1.29, b=1.35e+06)) # GENERATE SRR deviances devs <- ar1rlnorm(rho=0.4, 2018:2118, iters=200, meanlog=0, sdlog=0.5) # DEFINE MMY statistic statistic <- list(MMY=list(~apply(L,1,median), name="MMY", desc="ICES Maximum Median Yield")) # CALL bisect over 100 years, Fmmy calculated over last 50. fmmy <- opt.bisect(stock, sr=srr, deviances=devs, metrics=list(L=landings), statistic=statistic, years=2018:2118, pyears=2069:2118, tune=list(fbar=c(0.01, 0.2))) # fmmy mean(fbar(fmmy)[,ac(2069:2118)])data(ple4) stock <- propagate(stf(ple4, end=2118), 200) srr <- predictModel(model=rec ~ ifelse(ssb <= b, a * ssb, a * b), params=FLPar(a=1.29, b=1.35e+06)) # GENERATE SRR deviances devs <- ar1rlnorm(rho=0.4, 2018:2118, iters=200, meanlog=0, sdlog=0.5) # DEFINE MMY statistic statistic <- list(MMY=list(~apply(L,1,median), name="MMY", desc="ICES Maximum Median Yield")) # CALL bisect over 100 years, Fmmy calculated over last 50. fmmy <- opt.bisect(stock, sr=srr, deviances=devs, metrics=list(L=landings), statistic=statistic, years=2018:2118, pyears=2069:2118, tune=list(fbar=c(0.01, 0.2))) # fmmy mean(fbar(fmmy)[,ac(2069:2118)])
Applies sample_proc_dev() to each projected stock in an
FLStocks object and returns a compact data frame of F, depletion,
and process-error estimates.
pe_sampler(stks, brp, biomass = c("vb", "ssb"), eps = 1e-12)pe_sampler(stks, brp, biomass = c("vb", "ssb"), eps = 1e-12)
stks |
An |
brp |
An |
biomass |
Character. Biomass metric to use. Either |
eps |
Numeric. Small positive value used to avoid logarithms of zero or negative predicted biomass. |
This is a convenience wrapper around sample_proc_dev(). It is useful
for producing the main diagnostic relationship:
The output can be plotted directly to visualise how the emergent biomass process error changes with stock depletion or fishing pressure.
A data.frame with one row per F scenario and columns:
Median fishing mortality in the projected stock.
Sampled biomass process-error standard deviation.
Median biomass relative to Bmsy.
Median biomass relative to B0.
## Not run: data(ple4) ple4 <- merge_catch(ple4) bh <- srrTMB( as.FLSR(ple4, model = bevholt), spr0 = spr0y(ple4), r0 = NULL ) brp <- computeFbrp(ple4, bh, proxy = c("msy")) sim <- Fsim_grid( brp, sigmaR = 0.5, rho = 0.3, nyears = 100, iters = 100, parallel = TRUE, workers = 4 ) ps <- pe_sampler(sim, brp, biomass = "vb") head(ps) plot_pe(ps) ## End(Not run)## Not run: data(ple4) ple4 <- merge_catch(ple4) bh <- srrTMB( as.FLSR(ple4, model = bevholt), spr0 = spr0y(ple4), r0 = NULL ) brp <- computeFbrp(ple4, bh, proxy = c("msy")) sim <- Fsim_grid( brp, sigmaR = 0.5, rho = 0.3, nyears = 100, iters = 100, parallel = TRUE, workers = 4 ) ps <- pe_sampler(sim, brp, biomass = "vb") head(ps) plot_pe(ps) ## End(Not run)
sets plus group on FLQuant
pgquant(object, pg)pgquant(object, pg)
object |
FLQuant |
pg |
FLQuant
Plots sampled biomass process error against biomass relative to Bmsy. Optionally both axes can be normalised relative to a reference row, typically row 1, which corresponds to the Fmsy scenario in Fsim_grid().
plot_pe(ps, rel = TRUE, ref = 1, line = TRUE, points = TRUE, reverse_x = FALSE)plot_pe(ps, rel = TRUE, ref = 1, line = TRUE, points = TRUE, reverse_x = FALSE)
ps |
data.frame returned by pe_sampler(), with columns sigma and bbmsy. |
rel |
Logical. If TRUE, both B/Bmsy and sigma are divided by their reference values. |
ref |
Integer. Reference row used for normalisation. Defaults to 1. |
line |
Logical. Add line. |
points |
Logical. Add points. |
reverse_x |
Logical. Reverse x-axis. |
A ggplot object.
plots production functions
plot_pf(object, quant = c("vb", "ssb"), fmsy = NULL, rel = FALSE)plot_pf(object, quant = c("vb", "ssb"), fmsy = NULL, rel = FALSE)
object |
An *FLBRP* |
quant |
choose between vb and ssb or both |
fmsy |
default if Fmsy |
rel |
if TRUE ratios are produced for spcurve |
ggplot
data(ple4) ple4 <- merge_catch(ple4) sr <- srrTMB(as.FLSR(ple4,model=bevholtSV),spr0=mean(spr0y(ple4))) brp = brp(FLBRP(ple4,sr)) asem2spm(brp)[1:4] plot_pf(brp) plot_pf(brp,rel=TRUE)data(ple4) ple4 <- merge_catch(ple4) sr <- srrTMB(as.FLSR(ple4,model=bevholtSV),spr0=mean(spr0y(ple4))) brp = brp(FLBRP(ple4,sr)) asem2spm(brp)[1:4] plot_pf(brp) plot_pf(brp,rel=TRUE)
plotAdvice Plots stochastic stock dynamics against refpts for constant Fsim()
plotAdvice( object, rpts = "missing", type = NULL, yield = c("catch", "landings")[1], plotrefs = TRUE, probs = c(0.05, 0.2, 0.5, 0.8, 0.95), colour = "dodgerblue", ncol = NULL, osa = FALSE, label.size = 2.5, ssbQ = 1, recQ = 1 )plotAdvice( object, rpts = "missing", type = NULL, yield = c("catch", "landings")[1], plotrefs = TRUE, probs = c(0.05, 0.2, 0.5, 0.8, 0.95), colour = "dodgerblue", ncol = NULL, osa = FALSE, label.size = 2.5, ssbQ = 1, recQ = 1 )
type |
age-structured "asm" or surplus production "spm" plotting style |
yield |
option to select "catch" (default) or "landings" |
plotrefs |
if TRUE reference points are plotted |
probs |
determine credibility intervals, default 80th, 90th percentiles #' @param ncol number of plot panel columns |
colour |
color of CIs |
osa |
if TRUE it shows one-step-ahead for SSB and Rec |
label.size |
size of refpts labels |
ssbQ |
spawning quarter for seasonal models |
recQ |
recruitment quarter seasonal models |
stock |
FLStock or FLStockR |
refpts |
as FLPar or Fbrp() if FLStockR is not provided or should be overwritten |
ggplot
data(ple4) srr = srrTMB(as.FLSR(ple4,model=rickerSV),spr0=spr0y(ple4)) brp = computeFbrp(stock=ple4,sr=srr,proxy=c("sprx","f0.1","fe40"),blim=0.1,type="b0") plotAdvice (ple4,brp)data(ple4) srr = srrTMB(as.FLSR(ple4,model=rickerSV),spr0=spr0y(ple4)) brp = computeFbrp(stock=ple4,sr=srr,proxy=c("sprx","f0.1","fe40"),blim=0.1,type="b0") plotAdvice (ple4,brp)
plotAR Plots the new proposed ICES advice rule
plotAR( pars, ftgt = 1, btrigger = "missing", bpa = "missing", bthresh = "missing", fpa = "missing", fthresh = "missing", bclose = 0, fmin = 0, obs = "missing", kobe = TRUE, alpha = 1, xmax = 1.2, ymax = 1.5, ylab = "missing", xlab = "missing", rel = FALSE, expand = TRUE, labels = TRUE, label.cex = 3.5, critical = TRUE )plotAR( pars, ftgt = 1, btrigger = "missing", bpa = "missing", bthresh = "missing", fpa = "missing", fthresh = "missing", bclose = 0, fmin = 0, obs = "missing", kobe = TRUE, alpha = 1, xmax = 1.2, ymax = 1.5, ylab = "missing", xlab = "missing", rel = FALSE, expand = TRUE, labels = TRUE, label.cex = 3.5, critical = TRUE )
pars |
FLPar object or computeFbrp() ouput
|
ftgt |
factor to adjust Fmsy or its proxy e.g. 0.8Fmsy |
btrigger |
biomass trigger below which F is linearly reduced, if > 10 value, else factor*Btgt |
bpa |
precautionary biomass threshold, if > 10 value, else factor*Blim |
fpa |
option to input Fpa value |
bclose |
biomass that invokes fishing closure |
fmin |
minimum allowable (bycatch) fishing mortality under closure |
obs |
obtion to show observation with input class 'FLStock' |
kobe |
add kobe colour-coding |
alpha |
transparency of shading |
xmax |
multiplier for upper default xlim |
ymax |
multiplier for upper default ylim |
ylab |
option customize ylab |
xlab |
option customize xlab |
rel |
option to denote x,y labs as relative B/Btgt and F/Ftgt |
expand |
option to expand the plot area to border - default TRUE |
labels |
annotate reference point labels |
critical |
option to highlight critical zone below blim |
labelslabel.cex=3.5 |
set size of labels |
ggplot
data(ple4) srr = srrTMB(as.FLSR(ple4,model=segreg),spr0=spr0y(ple4)) blim = params(srr)[[2]] brp = computeFbrp(stock=ple4,sr=srr,proxy="f0.1",blim=blim) rpt = Fbrp(brp) plotAR(rpt,btrigger=an(0.8*rpt["Btgt"])) # Use Bpa as trigger (ICES style) plotAR(rpt,obs=ple4,bpa=1.4) # Change kobe to greyscale plotAR(rpt,obs=ple4,bpa=1.4,kobe=FALSE) # add fishing closure with minimum unavoidable F and Btrigger plotAR(rpt,obs=ple4,bpa=1.4,btrigger=0.7,kobe=TRUE,bclose=1,fmin=0.01) # show a relative plotAR(rpt,obs=ple4,rel=TRUE,bpa=1.4,btrigger=0.7,kobe=TRUE,bclose=1,fmin=0.02)data(ple4) srr = srrTMB(as.FLSR(ple4,model=segreg),spr0=spr0y(ple4)) blim = params(srr)[[2]] brp = computeFbrp(stock=ple4,sr=srr,proxy="f0.1",blim=blim) rpt = Fbrp(brp) plotAR(rpt,btrigger=an(0.8*rpt["Btgt"])) # Use Bpa as trigger (ICES style) plotAR(rpt,obs=ple4,bpa=1.4) # Change kobe to greyscale plotAR(rpt,obs=ple4,bpa=1.4,kobe=FALSE) # add fishing closure with minimum unavoidable F and Btrigger plotAR(rpt,obs=ple4,bpa=1.4,btrigger=0.7,kobe=TRUE,bclose=1,fmin=0.01) # show a relative plotAR(rpt,obs=ple4,rel=TRUE,bpa=1.4,btrigger=0.7,kobe=TRUE,bclose=1,fmin=0.02)
plotbioage() Plots stock N_a, W_a, M_a and Mat_a by year
plotbioage( stk, metrics = c("Weight", "Maturity", "M", "Selectivity"), bysex = FALSE, ncol = 2 )plotbioage( stk, metrics = c("Weight", "Maturity", "M", "Selectivity"), bysex = FALSE, ncol = 2 )
stk |
stock object class FLStock |
metrics |
choose Weight, Maturity, M, Selectivity |
bysex |
plot by sex default FALSE |
ncol |
number of columns in multiplot |
ggplot
data(ple4) plotbioage(ple4)data(ple4) plotbioage(ple4)
plotbioyr() Plots stock N_a, W_a, M_a and Mat_a across years
plotbioyr(stk, ncol = 2)plotbioyr(stk, ncol = 2)
stk |
stock object class FLStock |
ncol |
number of columns in multiplot |
ggplot
data(ple4) plotbioyr(ple4)data(ple4) plotbioyr(ple4)
Produces candidate GFCM advice plot
plotCECAF( ftgt = 1, btgt = 1, blim = 0.3, btrigger = "missing", bthr = 0.8, fthr = 1.2, bclose = 0.3, fmin = 0.2, fadv = 0.8, obs = "missing", kobe = TRUE, alpha = 1, xmax = 1.3, ymax = 1.5, ylab = "missing", xlab = "missing", rel = TRUE, expand = TRUE, labels = TRUE, critical = kobe, text = TRUE )plotCECAF( ftgt = 1, btgt = 1, blim = 0.3, btrigger = "missing", bthr = 0.8, fthr = 1.2, bclose = 0.3, fmin = 0.2, fadv = 0.8, obs = "missing", kobe = TRUE, alpha = 1, xmax = 1.3, ymax = 1.5, ylab = "missing", xlab = "missing", rel = TRUE, expand = TRUE, labels = TRUE, critical = kobe, text = TRUE )
btgt |
Btarget corresponding Ftgt |
blim |
biomass limit point used as fraction of Btgt, default 0.25-0.3Btgt |
btrigger |
biomass trigger (can differ from Bthr) |
bthr |
biomass precautionary point as fraction of Bmsy, default 0.5Btgt |
fthr |
upper Ftgt threshold, default 1.2 Ftgt |
bclose |
fishing closure |
fmin |
minimum F at closure |
fadv |
advice F from HCR |
obs |
"missing" option to show observed ssb,fba |
kobe |
if TRUE option to add KOBE colors |
alpha |
option for transparency |
xmax |
multiplier of btgt for ylim |
ymax |
multiplier of ftgt for xlim |
ylab |
|
xlab |
|
rel |
illustrate relative to btgt and ftgt if TRUE |
expand |
= TRUE graphic param |
labels |
if true add labels |
critical |
show area below Blim as dark red |
text |
show text |
ftg |
Ftarget |
GFCM advice plot
plotGFCM(fadv=0.95,btrigger=0.5)plotGFCM(fadv=0.95,btrigger=0.5)
plotdyn() Plots stock trajectories at age
plotdyn(stk, ncol = 2)plotdyn(stk, ncol = 2)
stk |
stock object class FLStock |
ncol |
number of columns in multiplot |
ggplot
data(ple4) plotdyn(ple4)data(ple4) plotdyn(ple4)
ploteq() Modification of method plot('FLBRP') to plot equilibrium output of computeFbrp()
ploteq( brps, refpts = "missing", obs = FALSE, rel = FALSE, rpf = TRUE, dashed = rpf, colours = "missing", panels = NULL, ncol = 2 )ploteq( brps, refpts = "missing", obs = FALSE, rel = FALSE, rpf = TRUE, dashed = rpf, colours = "missing", panels = NULL, ncol = 2 )
brps |
output object from computeFbrp of class FLBRP |
refpts |
Reference points, defaults are computed refpts from computeFbrp()
|
obs |
Should observations be plotted? Defaults to 'FALSE'. |
rel |
option to denote x,y labs as relative B/Btgt and F/Ftgt |
rpf |
adds refpts in plots |
dashed |
plots vertical dashed lines to highlight refpts locations |
colours |
refpts colours, default is designed for computeFbrp() output |
panels |
plot panel option 1:4 |
ncol |
number of plot panel columns |
ggplot
data(ple4) srr = srrTMB(as.FLSR(ple4,model=rickerSV),spr0=spr0y(ple4)) brp = computeFbrp(stock=ple4,sr=srr,proxy=c("sprx","f0.1","msy"),blim=0.1,type="b0") ploteq(brp,obs=TRUE) ploteq(brp,obs=TRUE,refpts="msy",rel=TRUE) brp.pa = computeFbrp(stock=ple4,sr=srr,proxy=c("msy","sprx","f0.1"),blim=0.1,bpa=Fbrp(brp)["Blim"]*2,type="b0") ploteq(brp.pa,obs=TRUE,rel=TRUE)data(ple4) srr = srrTMB(as.FLSR(ple4,model=rickerSV),spr0=spr0y(ple4)) brp = computeFbrp(stock=ple4,sr=srr,proxy=c("sprx","f0.1","msy"),blim=0.1,type="b0") ploteq(brp,obs=TRUE) ploteq(brp,obs=TRUE,refpts="msy",rel=TRUE) brp.pa = computeFbrp(stock=ple4,sr=srr,proxy=c("msy","sprx","f0.1"),blim=0.1,bpa=Fbrp(brp)["Blim"]*2,type="b0") ploteq(brp.pa,obs=TRUE,rel=TRUE)
plotFsim Plots stochastic stock dynamics against refpts for constant Fsim()
plotFsim( object, worms = TRUE, thinning = 10, probs = c(0.05, 0.2, 0.5, 0.8, 0.95), plotrefs = TRUE, colour = "missing", ncol = "missing", label.size = 3, yrs.eval = NULL, panels = "missing" )plotFsim( object, worms = TRUE, thinning = 10, probs = c(0.05, 0.2, 0.5, 0.8, 0.95), plotrefs = TRUE, colour = "missing", ncol = "missing", label.size = 3, yrs.eval = NULL, panels = "missing" )
object |
output object from Fsim() |
worms |
option to show individual iterations |
thinning |
thinning rate of iterations shows, e.g. 10 shows every 10th |
probs |
determine credibility intervals, default 80th, 90th percentiles |
plotrefs |
if TRUE reference points are plotted |
colour |
color of CIs |
ncol |
number of plot panel columns |
label.size |
size of reference points |
yrs.eval |
last years to be used evaluation period, default half nyears |
ggplot
Produces candidate GFCM advice plot
plotGFCM( ftgt = 1, btgt = 1, blim = 0.25, btrigger = "missing", bthr = 0.5, fthr = 1.2, bclose = 0, fmin = 0, fadv = 0.8, obs = "missing", kobe = TRUE, alpha = 1, xmax = 1.3, ymax = 1.5, ylab = "missing", xlab = "missing", rel = TRUE, expand = TRUE, labels = TRUE, critical = kobe, text = TRUE )plotGFCM( ftgt = 1, btgt = 1, blim = 0.25, btrigger = "missing", bthr = 0.5, fthr = 1.2, bclose = 0, fmin = 0, fadv = 0.8, obs = "missing", kobe = TRUE, alpha = 1, xmax = 1.3, ymax = 1.5, ylab = "missing", xlab = "missing", rel = TRUE, expand = TRUE, labels = TRUE, critical = kobe, text = TRUE )
btgt |
Btarget corresponding Ftgt |
blim |
biomass limit point used as fraction of Btgt, default 0.25-0.3Btgt |
btrigger |
biomass trigger (can differ from Bthr) |
bthr |
biomass precautionary point as fraction of Bmsy, default 0.5Btgt |
fthr |
upper Ftgt threshold, default 1.2 Ftgt |
bclose |
fishing closure |
fmin |
minimum F at closure |
fadv |
advice F from HCR |
obs |
"missing" option to show observed ssb,fba |
kobe |
if TRUE option to add KOBE colors |
alpha |
option for transparency |
xmax |
multiplier of btgt for ylim |
ymax |
multiplier of ftgt for xlim |
ylab |
|
xlab |
|
rel |
illustrate relative to btgt and ftgt if TRUE |
expand |
= TRUE graphic param |
labels |
if true add labels |
critical |
show area below Blim as dark red |
text |
show text |
ftg |
Ftarget |
GFCM advice plot
plotGFCM(fadv=0.95,btrigger=0.5)plotGFCM(fadv=0.95,btrigger=0.5)
Plots the new proposed ICES advice rule
plotMajuro( ftgt = 1, fthresh = 1.1, btgt = 1, blim = 0.1, btrigger = 0.8 * btgt, bthresh = 0.5 * btgt, bclose = 0, fmin = 0, obs = "missing", kobe = TRUE, alpha = 1, xmax = 1.5, ymax = 1.5, ylab = "missing", xlab = "missing", rel = FALSE, expand = TRUE, labels = TRUE, critical = kobe )plotMajuro( ftgt = 1, fthresh = 1.1, btgt = 1, blim = 0.1, btrigger = 0.8 * btgt, bthresh = 0.5 * btgt, bclose = 0, fmin = 0, obs = "missing", kobe = TRUE, alpha = 1, xmax = 1.5, ymax = 1.5, ylab = "missing", xlab = "missing", rel = FALSE, expand = TRUE, labels = TRUE, critical = kobe )
ftgt |
Target F = min(Fbrp,Fp0.5) |
btgt |
Biomass target corresponding to Fbrp |
blim |
biomass limit |
btrigger |
biomass trigger below which F is linearly reduced |
bthresh |
biomass threshold beyond which biomass is classified sustainable |
bclose |
biomass that invokes fishing closure |
fmin |
minimum allowable (bycatch) fishing mortality under closure |
obs |
obtion to show observation with input class 'FLStock' |
kobe |
add kobe colour-coding |
alpha |
transparency of shading |
xmax |
multiplier for upper default xlim |
ymax |
multiplier for upper default ylim |
ylab |
option customize ylab |
xlab |
option customize xlab |
rel |
option to denote x,y labs as relative B/Btgt and F/Ftgt |
expand |
option to expand the plot area to border - default TRUE |
labels |
annotate reference point labels |
critical |
option to highlight critical zone below blim |
ggplot
plotMajuro()plotMajuro()
plotspr() Plots current vs unfished spawning biomass per recruit at age
plotspr(stk, nyears = 3)plotspr(stk, nyears = 3)
stk |
stock object class FLStock |
nyears |
number of current last years, default is 3 |
ncol |
number of columns in multiplot |
ggplot
data(ple4) plotbioage(ple4)data(ple4) plotbioage(ple4)
plotWKREF Plots the newploiomass target corresponding to Fbrp
plotWKREF( ftgt = 1, btgt = 1, blim = 0.2, btrigger = 0.9 * btgt, bthresh = 0.8 * btgt, bclose = 0, fmin = 0, obs = "missing", kobe = TRUE, alpha = 1, xmax = 1.3, ymax = 1.5, ylab = "missing", xlab = "missing", rel = FALSE, expand = TRUE, labels = TRUE, critical = kobe )plotWKREF( ftgt = 1, btgt = 1, blim = 0.2, btrigger = 0.9 * btgt, bthresh = 0.8 * btgt, bclose = 0, fmin = 0, obs = "missing", kobe = TRUE, alpha = 1, xmax = 1.3, ymax = 1.5, ylab = "missing", xlab = "missing", rel = FALSE, expand = TRUE, labels = TRUE, critical = kobe )
blim |
biomass limit |
btrigger |
biomass trigger below which F is linearly reduced |
bthresh |
biomass threshold beyond which biomass is classified sustainable |
bclose |
ratio biomass/blim that invokes fishing closure relative to blim |
fmin |
minimum allowable (bycatch) fishing mortality under closure |
obs |
obtion to show observation with input class 'FLStock' |
kobe |
add kobe colour-coding |
alpha |
transparency of shading |
xmax |
multiplier for upper default xlim |
ymax |
multiplier for upper default ylim |
ylab |
option customize ylab |
xlab |
option customize xlab |
rel |
option to denote x,y labs as relative B/Btgt and F/Ftgt |
expand |
option to expand the plot area to border - default TRUE |
labels |
annotate reference point labels |
critical |
option to highlight critical zone below blim |
ggplot
plotWKREF() # Close fishery at Blim and adjust axis labels to relative plotWKREF(blim=0.2,bclose=0.2,rel=TRUE) # Close fishery at Blim, but allow fmin (e.g. bycatch) plotWKREF(blim=0.2,bclose=0.2,fmin=0.1,rel=TRUE) # Change Btrigger above Btgt plotWKREF(blim=0.2,bclose=0.2,fmin=0.1,btrigger=0.80,rel=TRUE) # Plot stock data data(ple4) plotWKREF(ftgt=0.25,btgt=8e+05,btrigger = 0.9*8e+05, blim=2e5,bclose=3e5,fmin=0.03,obs=ple4)plotWKREF() # Close fishery at Blim and adjust axis labels to relative plotWKREF(blim=0.2,bclose=0.2,rel=TRUE) # Close fishery at Blim, but allow fmin (e.g. bycatch) plotWKREF(blim=0.2,bclose=0.2,fmin=0.1,rel=TRUE) # Change Btrigger above Btgt plotWKREF(blim=0.2,bclose=0.2,fmin=0.1,btrigger=0.80,rel=TRUE) # Plot stock data data(ple4) plotWKREF(ftgt=0.25,btgt=8e+05,btrigger = 0.9*8e+05, blim=2e5,bclose=3e5,fmin=0.03,obs=ple4)
r4sscol
rc4(n, alpha = 1)rc4(n, alpha = 1)
n |
number of colors |
alpha |
transluscency |
vector of color codes
The plain bisection algorithm (Burden & Douglas, 1985) is employed here to find the value of a given forecast target quantity (e.g. 'fbar') for which a selected value of a performance statistic is obtained over a chosen period.
ref.bisect( stock, sr, deviances = rec(stock) %=% 1, metrics, refpts, statistic, years, pyears = years, tune, target = 1, tol = 0.01, maxit = 15, verbose = TRUE )ref.bisect( stock, sr, deviances = rec(stock) %=% 1, metrics, refpts, statistic, years, pyears = years, tune, target = 1, tol = 0.01, maxit = 15, verbose = TRUE )
stock |
object class FLStock |
sr |
object class FLSR |
metrics |
FLQuant of FLStock to be defined |
statistic |
|
years |
years to be evaluated |
tune |
range for input x |
target |
target with default 1 for ratios |
tol |
tolerance level |
maxit |
number of optimisation steps |
log |
if TRUE, optimise on log-scale |
Credits to Iago Mosqueira
Burden, Richard L.; Faires, J. Douglas (1985), "2.1 The Bisection Algorithm", Numerical Analysis (3rd ed.), PWS Publishers, ISBN 0-87150-857-5
rffwd() Project forward an FLStock with evolutionary Fbar
rffwd(object, sr, fbar = control, control = fbar, deviances = "missing")rffwd(object, sr, fbar = control, control = fbar, deviances = "missing")
object |
An *FLStock* |
sr |
A stock-recruit relationship, *FLSR* or *predictModel*. |
fbar |
Yearly target for average fishing mortality, *FLQuant*. |
control |
Yearly target for average fishing mortality, *FLPar*. |
deviances |
Deviances for the strock-recruit relationsip, *FLQuant*. |
The projected *FLStock* object.
data(ple4) sr <- srrTMB(as.FLSR(ple4,model=bevholtSV),spr0=mean(spr0y(ple4))) brp = computeFbrp(ple4,sr,proxy="msy") fbar(brp) = FLQuant(rep(0.01,70)) stk = as(brp,"FLStock") units(stk) = standardUnits(stk) its = 100 stk <- FLStockR(propagate(stk, its)) stk@refpts= Fbrp(brp) b0=an(Fbrp(brp)["B0"]) control = FLPar(Feq=0.15,Frate=0.1,Fsigma=0.15,SB0=b0,minyear=2,maxyear=70,its=its) run <- rffwd(stk, sr=sr,control=control,deviances=ar1rlnorm(0.3, 1:70, its, 0, 0.6)) plotAdvice(run)data(ple4) sr <- srrTMB(as.FLSR(ple4,model=bevholtSV),spr0=mean(spr0y(ple4))) brp = computeFbrp(ple4,sr,proxy="msy") fbar(brp) = FLQuant(rep(0.01,70)) stk = as(brp,"FLStock") units(stk) = standardUnits(stk) its = 100 stk <- FLStockR(propagate(stk, its)) stk@refpts= Fbrp(brp) b0=an(Fbrp(brp)["B0"]) control = FLPar(Feq=0.15,Frate=0.1,Fsigma=0.15,SB0=b0,minyear=2,maxyear=70,its=its) run <- rffwd(stk, sr=sr,control=control,deviances=ar1rlnorm(0.3, 1:70, its, 0, 0.6)) plotAdvice(run)
Function to characterize Productivity and refpts based on r and Generation
rGclass(r = NULL, gt = NULL)rGclass(r = NULL, gt = NULL)
r |
value of the intrinsic rate of population increase |
gt |
generation time G |
list with Productivity category and suggest Fbrps
sam2FLQuant()
sam2FLQuant( object, metric = c("ssb", "fbar", "catch", "rec")[1], what = c("mle") )sam2FLQuant( object, metric = c("ssb", "fbar", "catch", "rec")[1], what = c("mle") )
object |
fit from FLSAM |
forecast |
TRUE/FALSE |
FLQuant
adopted from Laurie Kell (biodyn)
sam2FLStockR()
sam2FLStockR(res, itsCI = 1000)sam2FLStockR(res, itsCI = 1000)
res |
fit from FLSAM |
itsCI |
number of iterations to depict uncertainty in plots |
FLStockR with refpts
sam2stars()
sam2stars(sam, refpts = NULL, quantiles = c(0.05, 0.95))sam2stars(sam, refpts = NULL, quantiles = c(0.05, 0.95))
sam |
fit from FLSAM |
refpts |
optional FLPar object with reference points |
STARS list with $timeseris and $refpts
Calculates biomass process deviations from a projected FLStock
by comparing the realised biomass transition with the deterministic
surplus-production expectation implied by asem2spm().
sample_proc_dev(stk, brp, biomass = c("vb", "ssb"), eps = 1e-12)sample_proc_dev(stk, brp, biomass = c("vb", "ssb"), eps = 1e-12)
stk |
A projected |
brp |
An |
biomass |
Character. Biomass metric to use. Either |
eps |
Numeric. Small positive value used to avoid taking the logarithm of zero or negative predicted biomass. |
The process deviation is calculated as:
where , , , and are obtained from
asem2spm(brp, quant = biomass).
This corresponds to sampling the residual process variation around the SPM-equivalent biomass transition:
A list with elements:
A one-row data.frame with median F, process-error
standard deviation, median B/Bmsy, and median B/B0.
An FLQuant of log biomass process deviations.
The equivalent surplus-production parameters from
asem2spm().
Standard deviation of the sampled process deviations.
Biomass relative to Bmsy.
Median biomass relative to B0.
The biomass metric used.
## Not run: sim <- Fsim_grid( brp, sigmaR = 0.5, rho = 0.3, nyears = 100, iters = 100 ) one <- sample_proc_dev( stk = sim[[1]], brp = brp, biomass = "vb" ) one$output one$sigma ## End(Not run)## Not run: sim <- Fsim_grid( brp, sigmaR = 0.5, rho = 0.3, nyears = 100, iters = 100 ) one <- sample_proc_dev( stk = sim[[1]], brp = brp, biomass = "vb" ) one$output one$sigma ## End(Not run)
generates a Schafer surplus production model with process and observation error
schaefer.sim( k = 10000, r = 0.3, q = 0.5, pe = 0.1, oe = 0.2, bk = 0.9, years = 1980:2022, f0 = 0.2, fhi = 2.2, flo = 0.8, sigmaF = 0.15, iters = 1, blim = 0.3, bthr = 0.5, rel = FALSE )schaefer.sim( k = 10000, r = 0.3, q = 0.5, pe = 0.1, oe = 0.2, bk = 0.9, years = 1980:2022, f0 = 0.2, fhi = 2.2, flo = 0.8, sigmaF = 0.15, iters = 1, blim = 0.3, bthr = 0.5, rel = FALSE )
k |
carrying capacity |
r |
intrinsic rate of population increase |
q |
catchability coefficient |
pe |
process error |
oe |
process error |
bk |
initial fraction of b/k |
years |
time horizon |
f0 |
factor for initial year as f0 = f/fmsy |
fhi |
factor for high F as fhi = f/fmsy |
flo |
factor for low F as flo = fbar/fmsy |
sigmaF |
variation on f trajectory |
iters |
number of iterations |
rel |
if TRUE metrics B/Bmsy and F/Fmsy are produced |
FLQuants
stk = schaefer.sim(iters=100,q=0.5) plotAdvice(stk) plot(FLIndex(index=iter(stk@stock,1))) # indexstk = schaefer.sim(iters=100,q=0.5) plotAdvice(stk) plot(FLIndex(index=iter(stk@stock,1))) # index
scales catch-at-age to total catch with error (optional)
sops(object, stock, sigma = 0.1, what = c("catch", "landings", "discards")[1])sops(object, stock, sigma = 0.1, what = c("catch", "landings", "discards")[1])
object |
FLQuant catch.n, discard.n, landings.n |
stock |
FLStock |
sigma |
observation error |
what |
type c("catch", "landings", "discards") |
FLQuant
spict2FLQuant()
spict2FLQuant( x, metric = c("ssb", "fbar", "catch", "stock", "harvest")[1], osa = FALSE, forecast = F, what = c("mle") )spict2FLQuant( x, metric = c("ssb", "fbar", "catch", "stock", "harvest")[1], osa = FALSE, forecast = F, what = c("mle") )
x |
fit from SPICT |
osa |
add one-step-ahead forecast |
forecast |
TRUE/FALSE |
what |
mle or log.sd |
FLQuant
adopted from Laurie Kell (biodyn)
spict2FLStockR()
spict2FLStockR( res, blim = 0.3, bthr = 0.5, rel = FALSE, osa = FALSE, forecast = NULL, itsCI = 1, seed = 123 )spict2FLStockR( res, blim = 0.3, bthr = 0.5, rel = FALSE, osa = FALSE, forecast = NULL, itsCI = 1, seed = 123 )
res |
fit from SPICT |
blim |
biomass limit reference point as fraction of Bmsy |
bthr |
biomass precautionary reference point as fraction of Bmsy |
rel |
if TRUE ratios BBmsy and FFmsy are stored |
osa |
add one-step-ahead forecast |
forecast |
extract forecast TRUE/FALSE |
itsCI |
number of iterations to depict uncertainty in plots |
seed |
random seed for consistent sampling across scenarios |
FLStockR with refpts
spict2stars()
spict2stars(spict, blim = 0.3, bthr = 0.5, quantiles = c(0.05, 0.95))spict2stars(spict, blim = 0.3, bthr = 0.5, quantiles = c(0.05, 0.95))
spict |
fit from fit.spict() |
blim |
biomass limit point as fraction of Bmsy, default 0.3Bmsy (ICES) |
bthr |
biomass precautionary point as fraction of Bmsy, default 0.5Bmsy (ICES) |
STARS list with $timeseris and $refpts
ss2FLStockR()
ss2FLStockR(mvln, thin = 1, output = NULL, rel = FALSE)ss2FLStockR(mvln, thin = 1, output = NULL, rel = FALSE)
mvln |
output from ssmvln() |
thin |
thinnig rate of retained iters |
output |
expected outputs presented as "mle" or median of "iters" |
rel |
if TRUE ratios B/Btgt and F/Ftgt are shown |
FLStockR with refpts
ss2ices()
ss2ices(mvln, quantiles = c(0.05, 0.95))ss2ices(mvln, quantiles = c(0.05, 0.95))
mvln |
output of ssmvln() |
quantiles |
default is 95CIs as c(0.025,0.975) |
output |
choice c("iters","mle")[1] |
Fmsy |
if specified the ratio F/Fmsy is calculated (required for ensembles) |
Btrigger |
if specified the ratio SSB/Btrigger is calculated (required for ensembles) |
ICES list with $timeseris and $refpts
ss2stars()
ss2stars(mvln, output = c("iters", "mle")[2], quantiles = c(0.05, 0.95))ss2stars(mvln, output = c("iters", "mle")[2], quantiles = c(0.05, 0.95))
mvln |
output of ssmvln() |
output |
choice c("iters","mle")[1] |
quantiles |
default is 95CIs as c(0.025,0.975) |
STARS list with $timeseris and $refpts
ss3col
ss3col(n, alpha = 1)ss3col(n, alpha = 1)
n |
number of colors |
alpha |
transluscency |
vector of color codes
function to generate MVN assessment error deviations of F and SSB
ss3devs(om, vcv, Fphi = 0.423, bias.correct = TRUE, ...)ss3devs(om, vcv, Fphi = 0.423, bias.correct = TRUE, ...)
om |
*FLom* or *FLStock* object |
vcv |
covariance matrix from ss3vcv() |
Fphi |
autocorrelation of F error |
bias.correct |
lognormal bias correction if TRUE |
FLQuants with devs of F and SSB
Henning Winker (GFCM)
function to generate to extract variance-covariance matrix for F and SSB
ss3vcv(ss3rep)ss3vcv(ss3rep)
ss3rep |
from r4ss::SS_output |
covariance matrix for F and SSB end year
Henning Winker (GFCM)
function to generate uncertainty for Stock Synthesis
ssmvln( ss3rep, Fref = NULL, years = NULL, virgin = FALSE, mc = 1000, weight = 1, run = "MVLN", addprj = FALSE, ymax = NULL, xmax = NULL, legendcex = 1, verbose = TRUE, seed = 123, observed.catch = FALSE )ssmvln( ss3rep, Fref = NULL, years = NULL, virgin = FALSE, mc = 1000, weight = 1, run = "MVLN", addprj = FALSE, ymax = NULL, xmax = NULL, legendcex = 1, verbose = TRUE, seed = 123, observed.catch = FALSE )
ss3rep |
from r4ss::SS_output |
Fref |
Choice of Fratio c("MSY","Btgt","SPR","F01"), correponding to F_MSY and F_Btgt |
years |
single year or vector of years for mvln |
virgin |
if FALSE (default) the B0 base for Bratio is SSB_unfished |
mc |
number of monte-carlo simulations |
weight |
weighting option for model ensembles weight*mc |
run |
qualifier for model run |
addprj |
include forecast years |
ymax |
ylim maximum |
xmax |
xlim maximum |
verbose |
Report progress to R GUI? |
seed |
retains interannual correlation structure like MCMC |
observed.catch |
if FALSE expected catch is used |
out |
choice c("iters","mle") |
plot |
option to show plot |
legendcex=1 |
Allows to adjust legend cex |
output list of quant posteriors and mle's
Henning Winker (GFCM)
stock2ratios()
stock2ratios(object)stock2ratios(object)
object |
of class *FLStockR* |
FLStockR with ratios F/Ftgt and B/Btgt
stockMedians() Converts FLStock into simplified FLStock with Median FLQuants
stockMedians(object, FUN = iterMedians, ssbQ = 1, recQ = 1)stockMedians(object, FUN = iterMedians, ssbQ = 1, recQ = 1)
object |
of class *FLStock* or *FLStockR* or *FLStocks* |
FUN |
computes iterMedians, iterMeans |
ssbQ |
SSB quarter seasonal models |
recQ |
recruitment quarter seasonal models |
FLStockR with *FLQuants*
computes catch.n for a given harvest
updCatch.n(stk)updCatch.n(stk)
object |
An *FLStock* |
FLQuant
updates sr in brp after changing biology
updsr(object, s = 0.7, v = NULL)updsr(object, s = 0.7, v = NULL)
object |
An *FLBRP* |
s |
assumed steepness s |
v |
input option new SB0 |
FLBRP
data(ple4) sr <- srrTMB(as.FLSR(ple4,model=bevholtSV),spr0=mean(spr0y(ple4))) brp = FLBRP(ple4,sr) s = sr@SV[[1]] params(brp) # change m(brp) = Mlorenzen(stock.wt(brp),Mref=0.15) brpupd =updsr(brp,s) params(brp)data(ple4) sr <- srrTMB(as.FLSR(ple4,model=bevholtSV),spr0=mean(spr0y(ple4))) brp = FLBRP(ple4,sr) s = sr@SV[[1]] params(brp) # change m(brp) = Mlorenzen(stock.wt(brp),Mref=0.15) brpupd =updsr(brp,s) params(brp)
updstars()
updstars(star, newrefpts)updstars(star, newrefpts)
star |
output of star list with new refpoints |
newrefpts |
FLPar manually adjusted reference points |
STARS list with $timeseris and $refpts
Calculates vulnerable biomass as the stock numbers-at-age multiplied by stock weight-at-age and total selectivity. Total selectivity is defined as the sum of landings and discards selectivity, scaled to a maximum of one.
vbbrp(x)vbbrp(x)
x |
An object with |
Vulnerable biomass is calculated as:
where is total selectivity-at-age, calculated as
landings.sel(x) + discards.sel(x) and normalised to a maximum of one.
This helper is useful when deriving an SPM-equivalent biomass quantity based on vulnerable biomass rather than spawning biomass or total biomass.
An FLQuant containing vulnerable biomass by year and
iteration. Units are copied from stock(x).
## Not run: vb <- vbbrp(brp) plot(vb) ## End(Not run)## Not run: vb <- vbbrp(brp) plot(vb) ## End(Not run)