Nothing
#' Surplus production model with FMSY and MSY as leading parameters
#'
#' A surplus production model that uses only a time-series of catches and a relative abundance index
#' and coded in TMB. The base model, `SP`, is conditioned on catch and estimates a predicted index.
#' Continuous surplus production and fishing is modeled with sub-annual time steps which should approximate
#' the behavior of ASPIC (Prager 1994). The Fox model, `SP_Fox`, fixes BMSY/K = 0.37 (1/e).
#' The state-space version, `SP_SS` estimates annual deviates in biomass. An option allows for setting a
#' prior for the intrinsic rate of increase.
#' The function for the `spict` model (Pedersen and Berg, 2016) is available in [MSEextra][MSEtool::MSEextra].
#'
#' @param x An index for the objects in `Data` when running in [runMSE][MSEtool::runMSE].
#' Otherwise, equals to 1 When running an assessment interactively.
#' @param Data An object of class Data.
#' @param AddInd A vector of integers or character strings indicating the indices to be used in the model. Integers assign the index to
#' the corresponding index in Data@@AddInd, "B" (or 0) represents total biomass in Data@@Ind, "VB" represents vulnerable biomass in
#' Data@@VInd, and "SSB" represents spawning stock biomass in Data@@SpInd.
#' @param rescale 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.
#' @param start Optional list of starting values. Entries can be expressions that are evaluated in the function. See details.
#' @param prior A named list for the parameters of any priors to be added to the model. See details.
#' @param fix_dep Logical, whether to fix the initial depletion (ratio of biomass to carrying capacity in the
#' first year of the model). If `TRUE`, uses the value in `start`, otherwise equal to 1
#' (unfished conditions).
#' @param fix_n Logical, whether to fix the exponent of the production function. If `TRUE`,
#' uses the value in `start`, otherwise equal to `n = 2`, where the biomass at MSY
#' is half of carrying capacity.
#' @param fix_sigma 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`.
#' @param fix_tau Logical, the standard deviation of the biomass deviations is fixed. If `TRUE`,
#' tau is fixed to value provided in `start` (if provided), otherwise, equal to 0.1.
#' @param early_dev Character string describing the years for which biomass deviations are estimated in `SP_SS`.
#' By default, deviations are estimated in each year of the model (`"all"`), while deviations could also be estimated
#' once index data are available (`"index"`).
#' @param LWT A vector of likelihood weights for each survey.
#' @param n_seas Integer, the number of seasons in the model for calculating continuous surplus production.
#' @param n_itF Integer, the number of iterations to solve F conditional on the observed catch given multiple seasons within an annual time step.
#' Ignored if `n_seas` = 1.
#' @param integrate Logical, whether the likelihood of the model integrates over the likelihood
#' of the biomass deviations (thus, treating it as a state-space variable).
#' @param Euler_Lotka Integer. If greater than zero, the function will calculate a prior for the intrinsic rate of increase to use in the estimation model
#' (in lieu of an explicit prior in argument `prior`). The value of this argument specifies the number of stochastic samples used to calculate the prior SD.
#' See section on priors below.
#' @param SR_type If `use_r_prior = TRUE`, the stock-recruit relationship used to calculate the stock-recruit alpha parameter from
#' steepness and unfished spawners-per-recruit. Used to develop the r prior.
#' @param silent Logical, passed to [TMB::MakeADFun()], whether TMB
#' will print trace information during optimization. Used for diagnostics for model convergence.
#' @param opt_hess Logical, whether the hessian function will be passed to [stats::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`.
#' @param n_restart The number of restarts (calls to [stats::nlminb()]) in the optimization procedure, so long as the model
#' hasn't converged. The optimization continues from the parameters from the previous (re)start.
#' @param control A named list of parameters regarding optimization to be passed to
#' [stats::nlminb()].
#' @param inner.control A named list of arguments for optimization of the random effects, which
#' is passed on to [newton][TMB::newton] via [TMB::MakeADFun()].
#' @param ... For `SP_Fox`, additional arguments to pass to `SP`.
#' @details
#' For `start` (optional), a named list of starting values of estimates can be provided for:
#' \itemize{
#' \item `MSY` Maximum sustainable yield.. Otherwise, 300% of mean catch by default.
#' \item `FMSY` Steepness. Otherwise, `Data@@Mort[x]` or 0.2 is used.
#' \item `dep` Initial depletion (B/B0) in the first year of the model. By default, 1.
#' \item `n` The production function exponent that determines BMSY/B0. By default, 2 so that BMSY/B0 = 0.5.
#' \item `sigma` Lognormal SD of the index (observation error). By default, 0.05. Not
#' used with multiple indices.
#' \item `tau` Lognormal SD of the biomass deviations (process error) in `SP_SS`. By default, 0.1.
#' }
#'
#' Multiple indices are supported in the model.
#'
#' Tip: to create the Fox model (Fox 1970), just fix n = 1. See example.
#' @section Priors:
#' The following priors can be added as a named list, e.g., prior = list(r = c(0.25, 0.15), MSY = c(50, 0.1). For each parameter below, provide a vector of values as described:
#'
#' \itemize{
#' \item `r` - A vector of length 2 for the lognormal prior mean (normal space) and SD (lognormal space).
#' \item `MSY` - A vector of length 2 for the lognormal prior mean (normal space) and SD (lognormal space).
#' }
#'
#' In lieu of an explicit r prior provided by the user, set argument `Euler_Lotka = TRUE` to calculate the prior mean and SD using
#' the Euler-Lotka method (Equation 15a of McAllister et al. 2001).
#' The Euler-Lotka method is modified to multiply the left-hand side of equation 15a by the alpha parameter of the
#' stock-recruit relationship (Stanley et al. 2009). Natural mortality and steepness are sampled in order to generate
#' a prior distribution for r. See `vignette("Surplus_production")` for more details.
#' @return An object of [Assessment-class] containing objects and output from TMB.
#' @note The model uses the Fletcher (1978) formulation and is parameterized with FMSY and MSY as
#' leading parameters. The default conditions assume unfished conditions in the first year of the time series
#' and a symmetric production function (n = 2).
#'
#' @section Online Documentation:
#' Model description and equations are available on the openMSE
#' [website](https://openmse.com/features-assessment-models/3-sp/).
#'
#' @author Q. Huynh
#' @references
#' Fletcher, R. I. 1978. On the restructuring of the Pella-Tomlinson system. Fishery Bulletin 76:515:521.
#'
#' Fox, W.W. 1970. An exponential surplus-yield model for optimizing exploited fish populations. Transactions of the American Fisheries Society 99:80-88.
#'
#' McAllister, M.K., Pikitch, E.K., and Babcock, E.A. 2001. Using demographic methods to construct Bayesian priors
#' for the intrinsic rate of increase in the Schaefer model and implications for stock rebuilding. Can. J. Fish.
#' Aquat. Sci. 58: 1871-1890.
#'
#' Pedersen, M. W. and Berg, C. W. 2017. A stochastic surplus production model in continuous time. Fish and Fisheries. 18:226-243.
#'
#' Pella, J. J. and Tomlinson, P. K. 1969. A generalized stock production model. Inter-Am. Trop. Tuna Comm., Bull. 13:419-496.
#'
#' Prager, M. H. 1994. A suite of extensions to a nonequilibrium surplus-production model. Fishery Bulletin 92:374-389.
#'
#' Stanley, R.D., M. McAllister, P. Starr and N. Olsen. 2009. Stock assessment for bocaccio (Sebastes
#' paucispinis) in British Columbia waters. DFO Can. Sci. Advis. Sec. Res. Doc. 2009/055. xiv + 200 p.
#' @section Required Data:
#' \itemize{
#' \item `SP`: Cat, Ind
#' \item `SP_SS`: Cat, Ind
#' }
#' @section Optional Data:
#' `SP_SS`: CV_Ind
#' @examples
#' data(swordfish)
#'
#' #### Observation-error surplus production model
#' res <- SP(Data = swordfish)
#'
#' # Provide starting values, assume B/K = 0.875 in first year of model
#' # and symmetrical production curve (n = 2)
#' start <- list(dep = 0.875, n = 2)
#' res <- SP(Data = swordfish, start = start)
#'
#' \donttest{
#' plot(res)
#' profile(res, FMSY = seq(0.1, 0.4, 0.01))
#' retrospective(res)
#' }
#'
#' #### State-space version
#' res_SS <- SP_SS(Data = swordfish, start = list(dep = 0.875, sigma = 0.1, tau = 0.1))
#'
#' \donttest{
#' plot(res_SS)
#' }
#'
#' #### Fox model
#' res_Fox <- SP(Data = swordfish, start = list(n = 1), fix_n = TRUE)
#' res_Fox2 <- SP_Fox(Data = swordfish)
#'
#' #### SP with r prior calculated internally (100 stochastic samples to get prior SD)
#' res_prior <- SP(Data = SimulatedData, Euler_Lotka = 100)
#'
#' #### Pass an r prior to the model with mean = 0.35, lognormal sd = 0.10
#' res_prior2 <- SP(Data = SimulatedData, prior = list(r = c(0.35, 0.10)))
#'
#' #### Pass MSY prior to the model with mean = 1500, lognormal sd = 0.05
#' res_prior3 <- SP(Data = SimulatedData, prior = list(MSY = c(1500, 0.05)))
#' @seealso [SP_production] [plot.Assessment] [summary.Assessment] [retrospective] [profile] [make_MP]
#' @export
SP <- function(x = 1, Data, AddInd = "B", rescale = "mean1", start = NULL, prior = list(),
fix_dep = TRUE, fix_n = TRUE, LWT = NULL,
n_seas = 4L, n_itF = 3L, Euler_Lotka = 0L, SR_type = c("BH", "Ricker"),
silent = TRUE, opt_hess = FALSE, n_restart = ifelse(opt_hess, 0, 1),
control = list(iter.max = 5e3, eval.max = 1e4), ...) {
SP_(x = x, Data = Data, AddInd = AddInd, state_space = FALSE, rescale = rescale, start = start, prior = prior,
fix_dep = fix_dep, fix_n = fix_n, fix_sigma = TRUE,
fix_tau = TRUE, LWT = LWT, n_seas = n_seas, n_itF = n_itF,
Euler_Lotka = Euler_Lotka, SR_type = SR_type, integrate = FALSE,
silent = silent, opt_hess = opt_hess, n_restart = n_restart, control = control, inner.control = list(), ...)
}
class(SP) <- "Assess"
#' @rdname SP
#' @export
SP_SS <- function(x = 1, Data, AddInd = "B", rescale = "mean1", start = NULL, prior = list(),
fix_dep = TRUE, fix_n = TRUE, fix_sigma = TRUE,
fix_tau = TRUE, LWT = NULL, early_dev = c("all", "index"), n_seas = 4L, n_itF = 3L,
Euler_Lotka = 0L, SR_type = c("BH", "Ricker"), integrate = FALSE,
silent = TRUE, opt_hess = FALSE, n_restart = ifelse(opt_hess, 0, 1),
control = list(iter.max = 5e3, eval.max = 1e4), inner.control = list(), ...) {
SP_(x = x, Data = Data, AddInd = AddInd, state_space = TRUE, rescale = rescale, start = start, prior = prior,
fix_dep = fix_dep, fix_n = fix_n, fix_sigma = fix_sigma,
fix_tau = fix_tau, early_dev = early_dev, LWT = LWT, n_seas = n_seas, n_itF = n_itF,
Euler_Lotka = Euler_Lotka,
SR_type = SR_type, integrate = integrate, silent = silent, opt_hess = opt_hess, n_restart = n_restart,
control = control, inner.control = inner.control, ...)
}
class(SP_SS) <- "Assess"
#' @rdname SP
#' @export
SP_Fox <- function(x = 1, Data, ...) {
SP_args <- c(x = x, Data = Data, list(...))
SP_args$start$n <- 1
SP_args$fix_n <- TRUE
do.call2("SP", SP_args)
}
class(SP_Fox) <- "Assess"
#' @useDynLib SAMtool
SP_ <- function(x = 1, Data, AddInd = "B", state_space = FALSE, rescale = "mean1", start = NULL, prior = list(),
fix_dep = TRUE, fix_n = TRUE, fix_sigma = TRUE,
fix_tau = TRUE, early_dev = c("all", "index"), LWT = NULL, n_seas = 4L, n_itF = 3L,
Euler_Lotka = 0L, SR_type = c("BH", "Ricker"), integrate = FALSE,
silent = TRUE, opt_hess = FALSE, n_restart = ifelse(opt_hess, 0, 1),
control = list(iter.max = 5e3, eval.max = 1e4), inner.control = list(), ...) {
dependencies = "Data@Cat, Data@Ind"
dots <- list(...)
start <- lapply(start, eval, envir = environment())
if (all(c("use_r_prior", "r_reps") %in% names(dots))) { # Backward compatibility
if (dots$use_r_prior) Euler_Lotka <- dots$r_reps
if ("r_prior" %in% names(start)) prior$r <- start$r_prior
}
early_dev <- match.arg(early_dev)
if (any(names(dots) == "yind")) {
yind <- eval(dots$yind)
} else {
ystart <- which(!is.na(Data@Cat[x, ]))[1]
yind <- ystart:length(Data@Cat[x, ])
}
Year <- Data@Year[yind]
C_hist <- Data@Cat[x, yind]
if (any(is.na(C_hist))) stop('Model is conditioned on complete catch time series, but there is missing catch.')
ny <- length(C_hist)
if (rescale == "mean1") rescale <- 1/mean(C_hist)
Ind <- lapply(AddInd, Assess_I_hist, Data = Data, x = x, yind = yind)
I_hist <- vapply(Ind, getElement, numeric(ny), "I_hist")
if (is.null(I_hist) || all(is.na(I_hist))) stop("No indices found.", call. = FALSE)
I_sd <- vapply(Ind, getElement, numeric(ny), "I_sd")
nsurvey <- ncol(I_hist)
if (state_space) {
if (early_dev == "all") est_B_dev <- rep(1, ny)
if (early_dev == "index") {
first_year_index <- which(apply(I_hist, 1, function(x) any(!is.na(x))))[1]
est_B_dev <- ifelse(1:ny < first_year_index, 0, 1)
}
} else {
if (nsurvey == 1 && all(AddInd == 0 | AddInd == "B")) {
fix_sigma <- FALSE # Override: estimate sigma if there's a single survey
}
est_B_dev <- rep(0, ny)
}
if (is.null(prior$r) && Euler_Lotka > 0L) {
r_samps <- r_prior_fn(x, Data, r_reps = Euler_Lotka, SR_type = SR_type)
prior$r <- c(mean(r_samps), sd(log(r_samps)))
}
prior <- make_prior_SP(prior)
if (is.null(LWT)) LWT <- rep(1, nsurvey)
if (length(LWT) != nsurvey) stop("LWT needs to be a vector of length ", nsurvey)
data <- list(model = "SP", C_hist = C_hist, rescale = rescale, I_hist = I_hist, I_sd = I_sd, I_lambda = LWT,
fix_sigma = as.integer(fix_sigma), nsurvey = nsurvey, ny = ny,
est_B_dev = est_B_dev, nstep = n_seas, dt = 1/n_seas, n_itF = n_itF,
use_prior = prior$use_prior, prior_dist = prior$pr_matrix,
sim_process_error = 0L)
params <- list()
if (!is.null(start)) {
if (!is.null(start$FMSY) && is.numeric(start$FMSY)) params$log_FMSY <- log(start$FMSY[1])
if (!is.null(start$MSY) && is.numeric(start$MSY)) params$MSYx <- log(start$MSY[1])
if (!is.null(start$dep) && is.numeric(start$dep)) params$log_dep <- log(start$dep[1])
if (!is.null(start$n) && is.numeric(start$n)) params$log_n <- log(start$n[1])
if (!is.null(start$sigma) && is.numeric(start$sigma)) params$log_sigma <- log(start$sigma)
if (!is.null(start$tau) && is.numeric(start$tau)) params$log_tau <- log(start$tau[1])
}
if (is.null(params$log_FMSY)) params$log_FMSY <- ifelse(is.na(Data@Mort[x]), 0.2, 0.5 * Data@Mort[x]) %>% log()
if (is.null(params$MSYx)) params$MSYx <- mean(3 * C_hist * rescale) %>% log()
if (is.null(params$log_dep)) params$log_dep <- log(1)
if (is.null(params$log_n)) params$log_n <- log(2)
if (is.null(params$log_sigma)) params$log_sigma <- rep(log(0.05), nsurvey)
if (is.null(params$log_tau)) params$log_tau <- log(0.1)
params$log_B_dev <- rep(0, ny)
map <- list()
if (fix_dep) map$log_dep <- factor(NA)
if (fix_n) map$log_n <- factor(NA)
if (fix_sigma) map$log_sigma <- factor(rep(NA, nsurvey))
if (fix_tau) map$log_tau <- factor(NA)
if (any(!est_B_dev)) map$log_B_dev <- factor(ifelse(est_B_dev, 1:sum(est_B_dev), NA))
random <- NULL
if (integrate) random <- "log_B_dev"
info <- list(Year = Year, data = data, params = params, control = control, inner.control = inner.control)
obj <- MakeADFun(data = info$data, parameters = info$params, hessian = TRUE,
map = map, random = random, DLL = "SAMtool", silent = silent)
high_F <- try(obj$report(c(obj$par, obj$env$last.par[obj$env$random]))$penalty > 0 ||
any(is.na(obj$report(c(obj$par, obj$env$last.par[obj$env$random]))$F)), silent = TRUE)
if (!is.character(high_F) && !is.na(high_F) && high_F) {
for(ii in 1:10) {
obj$par["MSYx"] <- 0.5 + obj$par["MSYx"]
if (all(!is.na(obj$report(obj$par)$F)) &&
obj$report(c(obj$par, obj$env$last.par[obj$env$random]))$penalty == 0) break
}
}
mod <- optimize_TMB_model(obj, control, opt_hess, n_restart)
opt <- mod[[1]]
SD <- mod[[2]]
report <- obj$report(obj$env$last.par.best)
Yearplusone <- c(Year, max(Year) + 1)
nll_report <- ifelse(is.character(opt), ifelse(integrate, NA, report$nll), opt$objective)
report$dynamic_SSB0 <- SP_dynamic_SSB0(obj) %>%
structure(names = Yearplusone)
Assessment <- new("Assessment", Model = ifelse(state_space, "SP_SS", "SP"),
Name = Data@Name, conv = SD$pdHess,
FMSY = report$FMSY, MSY = report$MSY,
BMSY = report$BMSY, VBMSY = report$BMSY, SSBMSY = report$BMSY,
B0 = report$K, VB0 = report$K, SSB0 = report$K,
FMort = structure(report$F, names = Year),
F_FMSY = structure(report$F/report$FMSY, names = Year),
B = structure(report$B, names = Yearplusone),
B_BMSY = structure(report$B/report$BMSY, names = Yearplusone),
B_B0 = structure(report$B/report$K, names = Yearplusone),
VB = structure(report$B, names = Yearplusone),
VB_VBMSY = structure(report$B/report$BMSY, names = Yearplusone),
VB_VB0 = structure(report$B/report$K, names = Yearplusone),
SSB = structure(report$B, names = Yearplusone),
SSB_SSBMSY = structure(report$B/report$BMSY, names = Yearplusone),
SSB_SSB0 = structure(report$B/report$K, names = Yearplusone),
Obs_Catch = structure(C_hist, names = Year),
Obs_Index = structure(I_hist, dimnames = list(Year, paste0("Index_", 1:nsurvey))),
Catch = structure(report$Cpred, names = Year),
Index = structure(report$Ipred, dimnames = list(Year, paste0("Index_", 1:nsurvey))),
NLL = structure(c(nll_report, report$nll_comp, report$penalty, report$prior),
names = c("Total", paste0("Index_", 1:nsurvey), "Dev", "Penalty", "Prior")),
info = info, obj = obj, opt = opt, SD = SD, TMB_report = report,
dependencies = dependencies)
if (state_space) {
Assessment@Dev <- structure(report$log_B_dev, names = Year)
Assessment@Dev_type <- "log-Biomass deviations"
Assessment@NLL <- structure(c(nll_report, report$nll_comp, report$penalty, report$prior),
names = c("Total", paste0("Index_", 1:nsurvey), "Dev", "Penalty", "Prior"))
} else {
Assessment@NLL <- structure(c(nll_report, report$nll_comp[1:nsurvey], report$penalty, report$prior),
names = c("Total", paste0("Index_", 1:nsurvey), "Penalty", "Prior"))
}
if (Assessment@conv) {
if (state_space) {
SE_Dev <- as.list(SD, "Std. Error")$log_B_dev
Assessment@SE_Dev <- structure(ifelse(is.na(SE_Dev), 0, SE_Dev), names = Year)
}
Assessment@SE_FMSY <- SD$sd[names(SD$value) == "FMSY"]
Assessment@SE_MSY <- SD$sd[names(SD$value) == "MSY"]
Assessment@SE_F_FMSY <- SD$sd[names(SD$value) == "F_FMSY_final"] %>% structure(names = max(Year))
Assessment@SE_B_BMSY <- Assessment@SE_SSB_SSBMSY <- Assessment@SE_VB_VBMSY <-
SD$sd[names(SD$value) == "B_BMSY_final"] %>% structure(names = max(Year))
Assessment@SE_B_B0 <- Assessment@SE_SSB_SSB0 <- Assessment@SE_VB_VB0 <-
SD$sd[names(SD$value) == "B_K_final"] %>% structure(names = max(Year))
catch_eq <- function(Ftarget) {
projection_SP(Assessment, Ftarget = Ftarget, p_years = 1, p_sim = 1, obs_error = list(matrix(1, 1, 1), matrix(1, 1, 1)),
process_error = matrix(1, 1, 1)) %>% slot("Catch") %>% as.vector()
}
Assessment@forecast <- list(catch_eq = catch_eq)
}
return(Assessment)
}
r_prior_fn <- function(x = 1, Data, r_reps = 1e2, SR_type = c("BH", "Ricker"), seed = x) {
SR_type <- match.arg(SR_type)
set.seed(x)
M <- trlnorm(r_reps, Data@Mort[x], Data@CV_Mort[x])
steep <- sample_steepness3(r_reps, Data@steep[x], Data@CV_steep[x], SR_type)
max_age <- Data@MaxAge
a <- Data@wla[x]
b <- Data@wlb[x]
Linf <- Data@vbLinf[x]
K <- Data@vbK[x]
t0 <- Data@vbt0[x]
La <- Linf * (1 - exp(-K * (c(1:max_age) - t0)))
Wa <- a * La ^ b
A50 <- min(0.5 * max_age, iVB(t0, K, Linf, Data@L50[x]))
A95 <- max(A50+0.5, iVB(t0, K, Linf, Data@L95[x]))
mat_age <- 1/(1 + exp(-log(19) * (c(1:max_age) - A50)/(A95 - A50)))
mat_age <- mat_age/max(mat_age)
log_r <- vapply(1:r_reps, function(y) uniroot(Euler_Lotka_fn, c(-6, 2), M = M[y], h = steep[y], weight = Wa,
mat = mat_age, maxage = max_age, SR_type = SR_type)$root, numeric(1))
return(exp(log_r))
}
Euler_Lotka_fn <- function(log_r, M, h, weight, mat, maxage, SR_type) {
M <- rep(M, maxage)
NPR <- calc_NPR(exp(-M), maxage)
SBPR <- sum(NPR * weight * mat)
CR <- ifelse(SR_type == "BH", 4*h/(1-h), (5*h)^1.25)
alpha <- CR/SBPR
EL <- alpha * sum(NPR * weight * mat * exp(-exp(log_r) * c(1:maxage)))
return(EL - 1)
}
SP_dynamic_SSB0 <- function(obj, par = obj$env$last.par.best, ...) {
dots <- list(...)
newdata <- obj$env$data
newdata$C_hist <- rep(1e-8, newdata$ny)
newparams <- clean_tmb_parameters(obj)
newparams$log_dep <- log(1)
obj2 <- MakeADFun(data = newdata, parameters = newparams, map = obj$env$map,
random = obj$env$random, DLL = "SAMtool", silent = TRUE)
obj2$report(par)$B
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.