R/eqsr_fit.R

In ices-tools-prod/msy: Estimation of Equilibrium Reference Points for Fsisheries

#' Stock recruitment fitting
#'
#' Fits one or more stock recruitment relationship to data containted in an
#' FLStock object. If more than one stock recruit relationship is provided, the
#' models are weighted based on smooth AIC weighting (See Buckland et al., 1997).
#'
#' @param stk FLStock object
#' @param nsamp Number of samples (iterations) to take from the stock recruitment
#'              fit (default is 1000).  If 0 (zero) then only the fits to the
#'              data are returned and no simulations are made.
#' @param models A character vector containing stock recruitment models to use
#'               in the model averaging. User can set any combination of
#'               "Ricker", "Segreg", "Bevholt", "Smooth_hockey".
#' @param id.sr A character vector specifying an id or name for the stock
#'              recruitment fit being run. The default is to use the slot "name"
#'              in the stk parameter is provided
#' @param remove.years A vector specifying the years to remove from the model
#'                     fitting.
#' @param rshift lag ssb by aditional years (default = 0).  As an example, for
#'   some herring stocks, age 1 (1 winter ring) fish were spawned 2 years
#'   previously, in this case, rshift = 1.
#' @return A list containing the following objects:
#' \itemize{
#'   \item `sr.sto` data.frame containing the alpha (a), beta (b), cv and model
#'         names. The number of rows correspond to the value set of `nsamp` in
#'         the function call.
#'   \item `sr.det` The parameters in the stock recruitment model corresponding
#'         to the "best fit" of any given model.
#'   \item `stk` An FLStock object, same as provided as input by the user.
#'   \item `rby` A data.frame containing the recruitment (rec), spawning stock
#'         biomass (ssb) and year used in the fitting of the data.
#'   \item `id.sr` A string containing run name (taken from the `id.sr` argument)
#' }
#'
#' @references
#'
#' Buckland, S.T., K.P. Burnham & N.H. Augustin (1997). Model selection: An integral part of inference.
#' Biometrics 53, 603-618.
#' DOI: \href{https://doi.org/10.2307/2533961}{10.2307/2533961}
#'
#' @seealso
#' \code{\link{eqsr_plot}} plots a simulation of predictive recruitment
#' from the fit, and shows a summary of the contributions of each stock
#' recruitment model to the model average fit.
#'
#' @examples
#' data(icesStocks)
#' FIT <- eqsr_fit(icesStocks$saiNS,
#'                 nsamp = 0,
#'                 models = c("Ricker", "Segreg"))
#'
#' # summary of individual fits
#' FIT$sr.det
#' eqsr_plot(FIT)
#'
#' # fit a bounded segmented regression
#' Segreg_bounded  <- function(ab, ssb) {
#'   ab$b <- min_ssb + ab$b
#'   Segreg(ab, ssb)
#' }
#' min_ssb <- min(FLCore::ssb(icesStocks$saiNS))
#'
#' FIT <- eqsr_fit(icesStocks$saiNS,
#'                 nsamp = 0,
#'                 models = c("Segreg", "Segreg_bounded"))
#'
#' # summary of individual fits
#' FIT$sr.det
#' FIT$sr.det$b[2] + min_ssb
#' eqsr_plot(FIT)
#'
#' \dontrun{
#' FIT <- eqsr_fit(icesStocks$saiNS,
#'                 nsamp = 2000,
#'                 models = c("Segreg", "Segreg_bounded"))
#'
#' # summary of individual fits
#' FIT$sr.det
#' eqsr_plot(FIT)
#' }
#'
#' @export
eqsr_fit <- function(stk, nsamp = 1000, models = c("Ricker","Segreg","Bevholt"),
                     id.sr = FLCore::name(stk), remove.years = NULL, rshift = 0)
{
  # some checks on the model argument
  if (!is.character(models)) stop("models arg should be character vector giving names of stock recruit models")
  if(any(models %in% c("ricker","segreg","bevholt"))) {
    stop("Please note that the msy stock-recruitment functions have been renamed:
   ricker -> Ricker
   bevholt -> Bevholt
   segreg ->  Segreg
   smooth_hockey -> Smooth_hockey
   This was done to resolve conflicts with same named functions in the FLCore-package.
   ERGO: use a capital in the first letter if you want to call these functions")
  }

  # get correct recruitment vector for each SSB
  # dims$min is the minimum age => recruitment age
  dms <- FLCore::dims(stk)
  rec <- c(FLCore::rec(stk))
  ssb_lag <- dms$min + rshift
  if (ssb_lag > 0)
  {
    rec <- c(rec[-seq(ssb_lag)], rep(NA, ssb_lag))
  }

  # combine all required data together
  # note year is the year that SSB had that value
  data <-
    data.frame(year = with(dms, minyear:maxyear),
               rec = rec,
               ssb = c(FLCore::ssb(stk)),
               fbar = c(FLCore::fbar(stk)),
               landings = c(FLCore::landings(stk)),
               catch = c(FLCore::catch(stk)),
               ssb_lag = ssb_lag)

  # remove years with NA recruitment
  data <- data[stats::complete.cases(data),]

  # which years to use in the fit
  data$remove.years <- FALSE
  if (!is.null(remove.years)) {
    remove.years <- data$year[data$year %in% remove.years]
    message("removing (ssb) years:\n\t",
            paste(remove.years, collapse = ", "),
            "\n  from the recruitment fitting procedure.")
    data$remove.years <- data$year %in% remove.years
  }

  # run model averaging
  srfit <- eqsr_Buckland(data[!data$remove.years,c("year", "rec", "ssb")],
                         nsamp,
                         models)

  # create output object
  out <- c(srfit, list(stk = stk, rby = data, id.sr = id.sr))
  class(out) <- c("eqsr_fit", c("list", "vector"))

  out
}