R/est_growth.R

In nwfsc-assess/nwfscSurvey: Northwest Fisheries Science Center Survey

#' Estimate length-at-age using the von Bertanlaffy
#' parametization
#'
#'
#' @template dir
#' @param dat Data frame of biological data from [pull_bio()]
#' @param return_df TRUE/FALSE If set to TRUE the dat data frame is
#' returned with added columns for the estimated Lhat_low, Lhat_pred, and Lhat_high.
#' @param bySex Logical to indicate if plot by sex
#' @param Par  Data frame of starting parameters for K, Linf, L0, CV0, and CV2 based on the
#' Stock Synthesis parameterization of von Bertanlaffy growth.
#' @param estVB Logical. Estimate vonB growth to plot against predicted length. If F, it uses the paramters in \code{parStart}.
#' @param bins The bins to put ages into. If NULL then simply uses the ages as recorded.
#' @param sdFactor The number of standard deviations to include in the
#' low and high calculations. The default is 1.0.
#' @param dopng Deprecated with {nwfscSurvey} 2.1 because providing a non-NULL
#'   value to `dir` can serve the same purpose as `dopng = TRUE` without the
#'   potential for errors when `dopng = TRUE` and `dir = NULL`. Thus, users
#'   no longer have to specify `dopng` to save the plot as a png.
#'
#' @seealso [fit_vbgrowth()]
#'
#' @author Chantel Wetzel
#' @export

est_growth <- function(
    dir = NULL, dat, return_df = TRUE,
    Par = data.frame(K = 0.13, Linf = 55, L0 = 15, CV0 = 0.10, CV1 = 0.10),
    bySex = TRUE, estVB = TRUE, bins = NULL, sdFactor = 1,
    dopng = lifecycle::deprecated()) {
  if (lifecycle::is_present(dopng)) {
    lifecycle::deprecate_warn(
      when = "2.1",
      what = "nwfscSurvey::PlotMap.fn(dopng =)"
    )
  }

  #### Initialize the three new columns
  dat$Lhat_pred <- NA
  dat$Lhat_low <- NA
  dat$Lhat_high <- NA

  if (is.null(bins)) {
    dat$Age_2 <- dat$Age
  } else {
    dat$Age_2 <- findInterval(dat$Age, bins)
  }

  if (bySex) {
    use_data <- which(
      !is.na(dat$Length_cm) &
        !is.na(dat$Age)
    )
    la_data <- dat[use_data, ]

    sex_list <- list(which(
      !is.na(dat$Length_cm) &
        !is.na(dat$Age) &
        dat$Sex %in% c("F")
    ), which(
      !is.na(dat$Length_cm) &
        !is.na(dat$Age) &
        dat$Sex %in% c("M")
    ), which(
      !is.na(dat$Length_cm) &
        !is.na(dat$Age) &
        dat$Sex %in% c("U")
    ))

    la_data_list <- list(
      female = la_data[la_data$Sex == "F", ],
      male = la_data[la_data$Sex == "M", ],
      unsexed = la_data[la_data$Sex == "U", ]
    )
    sex_vec <- c("F", "M", "U")
    nn <- 3
  } else {
    use_data <- which(
      !is.na(dat$Length_cm) &
        !is.na(dat$Age)
    )
    la_data <- dat[use_data, ]
    sex_list <- list(use_data)

    la_data_list <- list(allSex = la_data)
    nn <- 1
  }

  check_dir(dir = dir)

  out <- vector(mode = "list", length = nn)
  xpar <- vector(mode = "list", length = nn)
  names(out) <- paste0(names(la_data_list), "_sd_cv")
  names(xpar) <- paste0(names(la_data_list), "_growth")

  # Loop by sex
  for (i in 1:length(la_data_list)) {
    if (estVB) {
      ests_log <- stats::optim(
        fn = fit_vbgrowth,
        par = log(Par),
        hessian = FALSE,
        par_logspace = TRUE,
        Ages = la_data_list[[i]]$Age,
        Lengths = la_data_list[[i]]$Length_cm
      )$par
      xpar[[i]] <- exp(ests_log)
      # cat("Estimated VB parameters for", names(la_data_list)[i], xpar, "\n")
    } else {
      xpar[[i]] <- Par
    }

    # Predicts from 1 to the maximum observed age by sex
    # predL <- VB.fn(1:max(la_data_list[[i]]$Age), xpar[1], xpar[2], xpar[3])
    predL <- fit_vbgrowth(
      Par = xpar[[i]],
      par_logspace = FALSE,
      Ages = la_data_list[[i]]$Age,
      Lengths = la_data_list[[i]]$Length_cm,
      ReturnType = "Pred",
      sdFactor = sdFactor
    )
    rownames(predL) <- as.character(la_data_list[[i]]$Age)
    dat[sex_list[[i]], c("Lhat_low", "Lhat_pred", "Lhat_high")] <- predL

    x <- split(la_data_list[[i]]$Length_cm, la_data_list[[i]]$Age_2)
    xsd <- unlist(lapply(x, sd))
    xcv <- xsd / predL[names(xsd), "Lhat_pred"]
    if (is.null(bins)) {
      ages <- as.numeric(names(xsd))
    } else {
      ages <- bins[as.numeric(names(xsd))]
    }
    out[[i]] <- data.frame(ages = ages, sd = xsd, cv = xcv)
  } # end sex loop

  if (!is.null(dir)) {
    save(xpar, file = file.path(dir, "growth_vonB_estimates.Rdata"))
  }

  ests <- c(xpar, out)

  if (return_df) {
    return(dat)
  } else {
    return(ests)
  }
}