A generic statistical catch-at-age model (single fleet, single season) that uses catch, index, and catch-at-age composition data. An annual harvest rate is calculated (assuming a pulse fishery) as described in Forrest et al. (2008). There are two parameterizations for estimation of recruitment deviations, the stock-recruit relationship, and reference points (see functions section below).
SCA(x = 1, Data, SR = c("BH", "Ricker"), vulnerability = c("logistic",
"dome"), CAA_dist = c("multinomial", "lognormal"), CAA_multiplier = 50,
I_type = c("B", "VB", "SSB"), rescale = "mean1", start = NULL,
fix_h = TRUE, fix_U_equilibrium = TRUE, fix_sigma = FALSE,
fix_tau = TRUE, early_dev = c("comp_onegen", "comp", "all"),
late_dev = "comp50", integrate = FALSE, silent = TRUE,
opt_hess = FALSE, n_restart = ifelse(opt_hess, 0, 1),
control = list(iter.max = 2e+05, eval.max = 4e+05),
inner.control = list(), ...)SCA2(x = 1, Data, SR = c("BH", "Ricker"), vulnerability = c("logistic",
"dome"), CAA_dist = c("multinomial", "lognormal"), CAA_multiplier = 50,
I_type = c("B", "VB", "SSB"), rescale = "mean1", start = NULL,
fix_h = TRUE, fix_U_equilibrium = TRUE, fix_sigma = FALSE,
fix_tau = TRUE, common_dev = "comp50", integrate = FALSE,
silent = TRUE, opt_hess = FALSE, n_restart = ifelse(opt_hess, 0, 1),
control = list(iter.max = 2e+05, eval.max = 4e+05),
inner.control = list(), ...)
A position in the Data object (by default, equal to one for assessments).
An object of class Data
Stock-recruit function (either "BH" for Beverton-Holt or "Ricker").
Whether estimated vulnerability is "logistic" or "dome" (double-normal).
See details for parameterization.
Whether a multinomial or lognormal distribution is used for likelihood of the catch-at-age matrix. See details.
Numeric for data weighting of catch-at-age matrix if CAA_hist = "multinomial". Otherwise ignored. See details.
Whether the index surveys population biomass (B; this is the default in the DLMtool operating model), vulnerable biomass (VB), or spawning stock biomass (SSB).
A multiplicative factor that rescales the catch in the assessment model, which
can improve convergence. By default, "mean1" scales the catch so that time series mean is 1, otherwise a numeric.
Output is re-converted back to original units.
Optional list of starting values. See details.
Logical, whether to fix steepness to value in Data@steep in the model for SCA. This only affects
calculation of reference points for SCA2.
Logical, whether the equilibrium harvest rate prior to the first year of the model
is estimated. If TRUE, U_equilibrium is fixed to value provided in start (if provided),
otherwise, equal to zero (assumes virgin conditions).
Logical, whether the standard deviation of the index is fixed. If TRUE,
sigma is fixed to value provided in start (if provided), otherwise, value based on Data@CV_Ind.
Logical, the standard deviation of the recruitment deviations is fixed. If TRUE,
tau is fixed to value provided in start (if provided), otherwise, value based on Data@sigmaR.
Character string describing the years for which recruitment deviations are estimated in SCA. By default, "comp_onegen"
rec devs are estimated one full generation prior to the first year when catch-at-age (CAA) data are available. With "comp", rec devs are
estimated starting in the first year with CAA. With "all", rec devs start at the beginning of the model.
Typically, a numeric for the number of most recent years in which recruitment deviations will
not be estimated in SCA (recruitment in these years will be based on the mean predicted by stock-recruit relationship).
By default, "comp50" uses the number of ages (smaller than the mode)
for which the catch-at-age matrix has less than half the abundance than that at the mode.
Logical, whether the likelihood of the model integrates over the likelihood of the recruitment deviations (thus, treating it as a state-space variable).
Logical, passed to MakeADFun, whether TMB
will print trace information during optimization. Used for dignostics for model convergence.
Logical, whether the hessian function will be passed to nlminb during optimization
(this generally reduces the number of iterations to convergence, but is memory and time intensive and does not guarantee an increase
in convergence rate). Ignored if integrate = TRUE.
The number of restarts (calls to nlminb) in the optimization procedure, so long as the model
hasn't converged. The optimization continues from the parameters from the previous (re)start.
A named list of agruments for optimization to be passed to
nlminb.
A named list of arguments for optimization of the random effects, which
is passed on to newton.
Other arguments to be passed.
Typically, a numeric for the number of most recent years in which a common recruitment deviation will
be estimated (in SCA2, uninformative years will have a recruitment closer to the mean, which can be very misleading,
especially near the end of the time series). By default, "comp50" uses the number of ages (smaller than the mode)
for which the catch-at-age matrix has less than half the abundance than that at the mode.
SCA: The parameterization with R0 and steepness as leading parameters. Recruitment is estimated
as deviations from the resulting stock-recruit relationship.
SCA2: The mean recruitment in the time series is estimated and recruitment deviations around this mean are estimated
as penalized parameters (similar to Cadigan 2016). This version is generally very fast and robust. Virgin and MSY reference points
are estimated after the assessment run.
SCA: Cat, Ind, Mort, L50, L95, CAA, vbK, vbLinf, vbt0, wla, wlb, MaxAge
SCA2: Cat, Ind, Mort, L50, L95, CAA, vbK, vbLinf, vbt0, wla, wlb, MaxAge
SCA: Rec, steep, sigmaR, CV_Ind
SC2: Rec, steep, CV_Ind
For the statistical catch-at-age model, the basic data inputs are catch (by weight), index (by weight/biomass), and catch-at-age matrix (by numbers). Catches are assumed to be known perfectly (the harvest rate in a given year is the ratio of the observed catch to the vulnerable biomass at the beginning of the year). The maximum age in the model is a plus-group.
By default, steepness is fixed in the model to the value in Data@steep.
The annual sample sizes of the catch-at-age matrix is provided to the model (used in the
likelihood for catch-at-age, assuming a multinomial distribution),
and is manipulated via argument CAA_multiplier. This argument is
interpreted in two different ways depending on the value provided.
If CAA_multiplier > 1, then this value will cap the annual sample sizes
to that number. If CAA_multiplier <= 1, then all the annual samples sizes
will be re-scaled by that number. By default, sample sizes are capped at 50.
Alternatively, a lognormal distribution with inverse proportion variance can be used for the catch at age (Punt and Kennedy, 1994, as cited by Maunder 2011).
For start (optional), a named list of starting values of estimates can be provided for:
R0 Virgin recruitment, only for SCA.
h Steepness, only for SCA. If not provided, the value in Data@steep is used.
meanR Mean recruitment, only for SCA2.
U_equilibrium Harvest rate prior to the first year of model, e.g. zero means unfished conditions. Defaults to zero.
vul_par Vulnerability parameters (length 2 vector for logistic or length 4 for dome, see below). Users should provide
estimates of the parameters in normal space, e.g. vul_max between 0-1, and the function will perform the appropriate transformations for the model.
sigma Standard deviation of index. If not provided, the value based on Data@CV_Ind is used.
tau Standard deviation of recruitment deviations. If not provided, the value in Data@sigmaR is used.
Vulnerability can be specified to be either logistic or dome. If logistic, then the parameter
vector vul_par is of length 2:
vul_par[1]: a_95, the age of 95% vulnerability, via logit transformation to constrain a_95 to less than 75%
of the maximum age: a_95 = 0.75 * max_age * plogis(vul_par[1]).
vul_par[2]: a_50, the age of 50% vulnerability as an offset, i.e., a_50 = a_95 - exp(vul_par[2]).
A vague prior for vul_par[2] ~ N(0, sd = 3) is used to aid convergence, for example, when vulnerability >> 0.5 for the youngest age class.
With dome vulnerability, a double Gaussian parameterization is used, where vul_par
is an estimated vector of length 4:
vul_par[1]: a_asc, the first age of full vulnerability for the ascending limb, via logit transformation
to constrain a_95 to less than 75% of the maximum age: a_asc = 0.75 * maxage * plogis(vul_par[1]).
vul_par[2]: a_50, the age of 50% vulnerability for the ascending limb as an offset, i.e.,
a_50 = a_asc - exp(vul_par[2]).
vul_par[3]: a_des, the last age of full vulnerability (where the descending limb starts) via logit transformation
to constrain between a_asc and max_age,
i.e., a_des = (max_age - a_asc) * plogis(vul_par[3]) + a_asc.
vul_par[4]: vul_max, the vulnerability (in logit space) at the maximum age.
Vague priors of vul_par[2] ~ N(0, sd = 3) and vul_par[3] ~ N(0, 3) are used to aid convergence,
for example, when vulnerability >> 0.5 for the youngest age class.
Cadigan, N.G. 2016. A state-space stock assessment model for northern cod, including under-reported catches and variable natural mortality rates. Canadian Journal of Fisheries and Aquatic Science 72:296-308.
Forrest, R.E., Martell, S.J.D., Melnychuk, M.C., and Walters, C.J. 2008. An age-structured model with leading management parameters, incorporating age-specific selectivity and maturity. Canadian Journal of Fisheries and Aquatic Science 65:286-296.
Maunder, M.N. 2011. Review and evaluation of likelihood functions for composition data in stock-assessment models: Estimating the effective sample size. Fisheries Research 209:311-319.
Punt, A.E. and Kennedy, R.B. 1997. Population modelling of Tasmanian rock lobster, Jasus edwardsii, resources. Marine and Freshwater Research 48:967-980.
plot,Assessment,ANY-method summary,Assessment-method retrospective profile_likelihood make_MP
# NOT RUN {
res <- SCA(Data = DLMtool::SimulatedData)
res2 <- SCA2(Data = DLMtool::Simulation_1)
# }
Run the code above in your browser using DataLab