R/SVstd.R

In FESta: Fishing Effort Standardization

#' Standard Vessel Based Standardization (SVstd)
#'
#'
#' @importFrom stats aggregate
#' @importFrom stats as.formula
#' @importFrom stats var
#' @importFrom stats residuals
#' @importFrom stats model.matrix
#' @importFrom ggplot2 aes
#' @importFrom ggplot2 margin
#' @importFrom stats median
#' @importFrom dplyr ungroup
#' @param data A data frame containing the columns of year, vessel,
#' catch and effort.
#' See the example dataset \link[=SVstd_dataset]{SVstd_dataset}.
#'
#' @param year_col Specify the year column name (eg. "Year").
#' @param vessel_col Specify the vessel types column name (eg. "Vessel_Type").
#' @param catch_col Specify the catch column name (eg. "Catch").
#' @param effort_col Specify the effort column name (eg. "Effort").
#' @param standard_vessel Specify the reference (standard)
#'   vessel name (eg. "Vessel1") used to estimate relative fishing power. If \code{NULL}, the
#'   function automatically selects a standard vessel based on the
#'   available data.
#' @description
#' This method selects a reference (standard) vessel and estimates the
#' relative fishing power of all other vessels based on periods when both
#' the standard and comparison vessels operated simultaneously.
#'
#' The relative fishing power (RFP) for vessel \eqn{i} is calculated as:
#'
#' \deqn{
#' RFP_i = \frac{C_i/E_i}{C_s/E_s}
#' }
#'
#' where \eqn{C_i} and \eqn{E_i} are the total catch and effort of vessel
#' \eqn{i}, respectively, and \eqn{C_s} and \eqn{E_s} are the corresponding
#' catch and effort values for the selected standard vessel during the
#' same period.
#'
#' The standardized annual CPUE index for year \eqn{t} is then computed as:
#'
#' \deqn{
#' I_t = \frac{\sum_i C_{t,i}}
#'            {\sum_i RFP_i E_{t,i}}
#' }
#'
#' where \eqn{C_{t,i}} is the catch and \eqn{E_{t,i}} is the effort of
#' vessel \eqn{i} in year \eqn{t}.
#'
#' @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 catch column value of has zero value(s) then they will be replaced
#'by minimum value of the catch column and further
#'if effort column has zero values then the corresponding rows will be
#'removed for doing the CPUE calculation.
#'
#' @export
#' @references
#' Beverton, R.J.H., and Holt, S.J. (1957).
#' On the Dynamics of Exploited Fish Populations.
#' Fishery Investigations Series II, Volume XIX.
#'
#'Varghese, E., Jayasankar, J., Sathianandan, T.V., Kuriakose, S., Mini, K.G.,
#'Gills, R., Muktha, M., Sreepriya, V. and Gopalakrishnan, A. (2023).
#'A note on different methods for standardization of fishing efforts.
#'Marine Fisheries Information Service,
#'Technical and Extension Series, (257), 7-17.
#'
#'
#' Maunder, M.N., and Punt, A.E. (2004).
#' Standardizing catch and effort data:
#' a review of recent approaches.
#' Fisheries Research, 70, 141-159.
#'
#' \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("SVstd_dataset")
#'result<-SVstd(data=SVstd_dataset,year_col ="Year",  vessel_col = "Vessel",
#'catch_col = "Catch",effort_col = "Effort",standard_vessel = "V006")
#'print(result)
#'}
SVstd <- function(data,
                  year_col,
                  vessel_col,
                  catch_col,
                  effort_col,
                  standard_vessel = NULL) {

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


  # -------------------------------
  # Zero Handling
  # -------------------------------
  # 1. Remove rows with zero effort
  zero_effort_rows <- which(data[[effort_col]] == 0 | is.na(data[[effort_col]]))
  if (length(zero_effort_rows) > 0) {
    message(paste("Removed", length(zero_effort_rows),
                  "rows with zero or missing effort values"))
    data <- data[-zero_effort_rows, ]
  }

  # 2. Handle zero catch (replace with small positive value)
  zero_catch_rows <- which(data[[catch_col]] == 0 & !is.na(data[[catch_col]]))
  if (length(zero_catch_rows) > 0) {
    min_catch <- min(data[[catch_col]][data[[catch_col]] != 0], na.rm = TRUE)
    if (is.finite(min_catch)) {
      data[[catch_col]][zero_catch_rows] <- min_catch
      message(paste("Replaced", length(zero_catch_rows),
                    "zero catch values with minimum catch value ", round(min_catch, 4)))
    } else {
      # If all catches are zero, remove rows
      data <- data[-zero_catch_rows, ]
      message(paste("Removed", length(zero_catch_rows),
                    "rows with zero catch (no non zero catch found)"))
    }
  }

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


  # -------------------------------
  # Choose standard gear
  # -------------------------------
  if (is.null(standard_gear))
    standard_gear <- unique(data[[gear_col]])[1]
  cat("Standard Gear:", standard_gear, "\n\n")

  # -------------------------------
  # Overall CPUE per gear
  # -------------------------------
  gear_summary <- stats::aggregate(
    cbind(Catch = data[[catch_col]], Effort = data[[effort_col]]),
    by  = list(Gear = data[[gear_col]]),
    FUN = sum
  )
  gear_summary$CPUE <- gear_summary$Catch / gear_summary$Effort

  # -------------------------------
  # Standard gear CPUE
  # -------------------------------
  std_cpue <- gear_summary$CPUE[gear_summary$Gear == standard_gear]
  if (length(std_cpue) == 0)
    stop("Standard gear not found in data.")

  # -------------------------------
  # Relative Fishing Power (RFP)
  # -------------------------------
  gear_summary$RFP <- gear_summary$CPUE / std_cpue

  # -------------------------------
  # Merge RFP back to original data
  # -------------------------------
  data2 <- merge(
    data,
    gear_summary[, c("Gear", "RFP")],
    by.x = gear_col,
    by.y = "Gear"
  )

  # -------------------------------
  # Yearly standardized CPUE
  # -------------------------------
  years   <- sort(unique(data2[[year_col]]))
  results <- data.frame()

  for (y in years) {
    temp                <- data2[data2[[year_col]] == y, ]
    total_catch         <- sum(temp[[catch_col]])
    standardized_effort <- sum(temp$RFP * temp[[effort_col]])
    It                  <- total_catch / standardized_effort
    results <- rbind(
      results,
      data.frame(
        Year                = y,
        Total_Catch         = total_catch,
        Standardized_Effort = standardized_effort,
        Standardized_CPUE   = It
      )
    )
  }


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

}