A generic statistical catch-at-age model (single fleet, single season) that uses catch, index, and catch-at-age composition
data. SCA parameterizes R0 and steepness as leading productivity parameters in the assessment model. Recruitment is estimated
as deviations from the resulting stock-recruit relationship. In 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). The standard deviation is set high
so that the recruitment is almost like free parameters. Unfished and MSY reference points are inferred afterwards from the assessment output
(SSB and recruitment estimates). SCA_Pope is a variant of SCA that fixes the expected catch to the observed
catch, and Pope's approximation is used to calculate the annual harvest rate (U).
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",
max_age = Data@MaxAge,
start = NULL,
fix_h = TRUE,
fix_F_equilibrium = TRUE,
fix_omega = 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",
max_age = Data@MaxAge,
start = NULL,
fix_h = TRUE,
fix_F_equilibrium = TRUE,
fix_omega = 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(),
...
)
SCA_Pope(
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",
max_age = Data@MaxAge,
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(),
...
)
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.
Integer, the maximum age (plus-group) in the model.
Optional list of starting values. Entries can be expressions that are evaluated in the function. 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 fishing mortality prior to the first year of the model
is estimated. If TRUE, F_equilibrium is fixed to value provided in start (if provided),
otherwise, equal to zero (assumes unfished conditions).
Logical, whether the standard deviation of the catch is fixed. If TRUE,
sigma is fixed to value provided in start (if provided), otherwise, value based on Data@CV_Cat.
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.
Numeric or character string describing the years for which recruitment deviations are estimated in SCA. By default,
equal to "comp_onegen", where 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.
If numeric, the number of years after the first year of the model for which to start estimating rec devs. Use negative numbers for years prior to the first year.
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 random effects/state-space variable). Otherwise, recruitment deviations are penalized parameters.
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.
Logical, same as `fix_F_equilibrium` for `SCA_Pope`.
SCA, SCA_Pope, and SCA_Pope: Cat, Ind, Mort, L50, L95, CAA, vbK, vbLinf, vbt0, wla, wlb, MaxAge
SCA: Rec, steep, sigmaR, CV_Ind, CV_Cat
SC2: Rec, steep, CV_Ind, CV_Cat
SCA_Pope: Rec, steep, sigmaR, CV_Ind
The basic data inputs are catch (by weight), index (by weight/biomass), and catch-at-age matrix (by numbers).
Annual F's are estimated parameters assuming continuous fishing over the year. Note: prior to version 1.2, catches were assumed
to be known perfectly with an annual harvest rate from pulse fishing in SCA. That feature has now moved to SCA_Pope.
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.
F_equilibrium Fishing mortality 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.
F A vector of F's of length nyears, length(Data@Year). If not provided, defaults to 0.1.
omega Standard deviation of catch. If not provided, the value based on Data@CV_Cat is used.
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. By default, fixed to a small value so that the dome is effectively
a three-parameter function.
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.
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 summary.Assessment retrospective profile make_MP
# NOT RUN {
res <- SCA(Data = DLMtool::SimulatedData)
res2 <- SCA2(Data = DLMtool::SimulatedData)
compare_models(res, res2)
SCA_assess <- SCA2(Data = DLMtool::Simulation_1)
# }
# NOT RUN {
plot(res)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab