R/DEstd.R

In FESta: Fishing Effort Standardization

utils::globalVariables(
  c(
    "Year",
    "Standardized_CPUE"
  )
)

#' Derived Effort Based Standardization (DEstd)

#' @importFrom utils globalVariables
#' @importFrom rlang .data
#' @importFrom stats glm
#' @importFrom stats predict
#' @importFrom stats median
#' @importFrom stats aggregate
#' @importFrom stats as.formula
#' @importFrom stats var
#' @importFrom stats residuals
#' @importFrom stats gaussian
#' @importFrom stats poisson
#' @importFrom stats binomial
#' @importFrom stats Gamma
#' @importFrom stats glm.control
#' @importFrom stats model.matrix
#' @importFrom stats optim
#' @importFrom ggplot2 aes
#' @importFrom ggplot2 margin
#' @importFrom statmod tweedie
#' @importFrom dplyr group_by mutate ungroup
#' @param data A data frame containing the columns of year, gear, catch, effort and
#' total annual catch.
#' See the example dataset \link[=DEstd_dataset]{DEstd_dataset}.
#' @param year_col Specify the year column name (eg. "Year").
#' @param gear_col Specify the gear types column name (eg. "Gears").
#' @param catch_col Specify the catch column name (eg. "Catch").
#' @param effort_col Specify the effort column name (eg. "Effort").
#'@param total_catch_col Specify the total annual catch column name
#'(eg. "Annual Catch").
#' @note The unit of the standardized CPUE will be based on the units of
#' Catch and Effort.
#' @description
#' The derived effort approach (Sparre, 1998) assumes effort is a good measure
#' when it relates linearly to catch rate. Since different gears use
#' incompatible effort units, each is converted to CPUE and then to a relative
#' CPUE so they can be combined. Dividing total yield by the yield-weighted sum
#' of these relative CPUEs gives a standardized effort series that
#' reflects relative abundance.
#'
#' The catch per unit effort for gear \eqn{i} in year \eqn{y} is:
#'
#' \deqn{
#' CPUE_i(y)=\frac{Y_i(y)}{f_i(y)}
#' }
#'
#' where \eqn{Y_i(y)} is the catch and \eqn{f_i(y)} is the corresponding effort.
#'
#' Relative CPUE is computed as:
#'
#' \deqn{
#' R_i(y)=\frac{CPUE_i(y)}
#' {\mathrm{Mean}[CPUE_i]}
#' }
#'
#' where \eqn{\mathrm{Mean}[CPUE_i]} is the average CPUE of gear \eqn{i}
#' across all years.
#'
#' Annual relative effort is estimated as:
#'
#' \deqn{
#' R(y)=\sum_{i=1}^{k}
#' \left[
#' R_i(y)\times\frac{Y_i(y)}{Y_E(y)}
#' \right]
#' }
#'
#' where \eqn{Y_E(y)} is the total catch from gears for which effort
#' information is available.
#'
#' The Standardized CPUE is then calculated as:
#'
#' \deqn{
#' E(y)=
#' \frac{Y_T(y)/R(y)}
#' {\mathrm{Mean}[Y_T/R]}
#' }
#'
#' where \eqn{Y_T(y)} is the total annual catch
#' (including gears for which effort is not known).
#' @export
#' @return
#' The output includes a summary table containing Year, Total Catch,
#' Nominal CPUE, and Standardized CPUE. In addition, two plots are
#' produced: Nominal CPUE versus Total Catch and Standardized CPUE
#' versus Total Catch.
#'@note
#'If effort column has zero values then the rows corresponding
#'to them are removed to calculate CPUE.
#'
#' @references
#' Sparre, P., and Venema, S.C. (1992).
#' Introduction to Tropical Fish Stock Assessment.
#' FAO Fisheries Technical Paper No. 306/1, 376 pp.
#'
#' \strong{Acknowledgements:}

#'The authors sincerely thank the Director,

#'ICAR–Central Marine Fisheries Research Institute (ICAR-CMFRI), Kochi,

#'for providing the necessary facilities and institutional support.

#'The authors also gratefully acknowledge the support provided by

#'the Indian Council of Agricultural Research (ICAR),

#'Department of Agricultural Research and Education (DARE),

#'Government of India, through the ICAR-National Fellow Project.
#'
#' @examples
#'
#' \dontrun{
#'
#' library(FESta)
#' data("DEstd_dataset")
#' result<-DEstd(data=DEstd_dataset,year_col = "Year",gear_col = "Gear",
#' catch_col = "Catch",total_catch_col="Total_Catch", effort_col = "Effort")
#' print(result)
#' }
#'

DEstd <- function(
    data,
    year_col,
    gear_col,
    catch_col,
    effort_col,
    total_catch_col
) {

  # ---------------------------
  # Validation
  # ---------------------------
  required_cols <- c(year_col, gear_col, catch_col, effort_col, total_catch_col)
  missing_cols  <- setdiff(required_cols, names(data))
  if (length(missing_cols) > 0)
    stop("Missing columns: ", paste(missing_cols, collapse = ", "))

  # =============================================
  # HANDLE ZERO EFFORT ONLY
  # =============================================
  # Remove rows with zero or missing effort
  zero_effort_rows <- which(data[[effort_col]] == 0 | is.na(data[[effort_col]]))
  if (length(zero_effort_rows) > 0) {
    message(paste("NOTE: Removed", length(zero_effort_rows),
                  "rows with zero or missing effort values"))
    data <- data[-zero_effort_rows, ]
  }

  # Check if any data remains
  if (nrow(data) == 0) {
    stop("No valid observations remaining after removing zero effort rows.")
  }

  # ---------------------------
  # CPUE per row
  # ---------------------------
  data$CPUE <- data[[catch_col]] / data[[effort_col]]
  sampled_yield <- aggregate(
    data[[catch_col]],
    by = list(Year = data[[year_col]]),
    FUN = sum
  )

  names(sampled_yield) <- c("Year", "YS")
  data <- merge(
    data,
    sampled_yield,
    by.x = year_col,
    by.y = "Year"
  )
  # ---------------------------
  # Mean CPUE per gear
  # ---------------------------
  mean_cpue <- stats::aggregate(
    data$CPUE,
    by  = list(Gear = data[[gear_col]]),
    FUN = mean,
    na.rm = TRUE
  )
  names(mean_cpue)[2] <- "Mean_CPUE"

  data <- merge(data, mean_cpue,
                by.x = gear_col, by.y = "Gear",
                all.x = TRUE)

  # ---------------------------
  # Relative CPUE and Weighted Relative CPUE
  # ---------------------------
  data$Relative_CPUE_Ri <- data$CPUE / data$Mean_CPUE

  data$Weighted_Relative_CPUE <-
    data$Relative_CPUE_Ri *
    (data[[catch_col]] / data$YS)
  # ---------------------------
  # Yearly Relative Effort
  # ---------------------------
  yearly_relative_effort <- stats::aggregate(
    data$Weighted_Relative_CPUE,
    by  = list(Year = data[[year_col]]),
    FUN = sum,
    na.rm = TRUE
  )
  names(yearly_relative_effort)[2] <- "Relative_Effort_Ry"

  # ---------------------------
  # Bring in yearly total catch
  # ---------------------------
  yearly_totals_unique <- stats::aggregate(
    data[[total_catch_col]],
    by  = list(Year = data[[year_col]]),
    FUN = function(x) x[!is.na(x)][1]
  )
  names(yearly_totals_unique) <- c("Year", "Total_Catch_YT")


  yearly_data <- merge(yearly_relative_effort, yearly_totals_unique,
                       by = "Year")

  # ---------------------------
  # Standardized effort
  # ---------------------------
  yearly_data$YT_by_R <- yearly_data$Total_Catch_YT / yearly_data$Relative_Effort_Ry
  mean_YT_by_R        <- mean(yearly_data$YT_by_R, na.rm = TRUE)
  yearly_data$Standardized_CPUE <- yearly_data$YT_by_R / mean_YT_by_R

  yearly_data <- yearly_data[order(yearly_data$Year), ]

  # ---------------------------
  # Output summaries
  # ---------------------------
  gear_summary <- data.frame(
    Year                   = data[[year_col]],
    Gear                   = data[[gear_col]],
    CPUE                   = round(data$CPUE,                   4),
    Weighted_Relative_CPUE = round(data$Weighted_Relative_CPUE, 4)
  )

  # Keep all columns needed for plotting
  yearly_summary <- data.frame(
    Year                   = yearly_data$Year,
    Total_Catch_YT         = yearly_data$Total_Catch_YT,
    Relative_Effort_Ry     = yearly_data$Relative_Effort_Ry,
    Standardized_CPUE = yearly_data$Standardized_CPUE
  )

  calculate_and_plot_cpue(year = data[[year_col]],
                          std_cpue =yearly_data$Standardized_CPUE,effort = data[[effort_col]], total_catch = NULL,catch = data[[catch_col]] )

}