R/REstd.R

In FESta: Fishing Effort Standardization

#' Relative Effort Based Standardization (REstd)
#'
#' @importFrom rlang .data

#' @importFrom stats aggregate
#' @importFrom stats as.formula
#' @importFrom stats var
#' @importFrom stats residuals


#' @importFrom stats model.matrix
#' @importFrom stats median
#' @importFrom ggplot2 aes
#' @importFrom ggplot2 margin

#' @importFrom dplyr group_by mutate ungroup
#' @param data A data frame containing the columns of year, vessel information,
#' number of boats, days of fishing and CPUE.
#' See the example dataset \link[=REstd_dataset]{REstd_dataset}.
#' @param year_col Specify the year column name (eg. "Year").
#' @param vessel_col Specify the vessel types column name (eg. "Vessel_Type").
#' @param boats_col Specify the number of boats column name (eg. "Number_of_Boats").
#' @param days_col Specify the number of fishing days column name
#' (eg. "Avg_Fishing_Days").
#' @param cpue_col Specify the CPUE column name (eg. "CPUE").
#' @param standard_vessel Specifying the reference vessel based on which
#'  relative fishing power will be calculated. If NULL, the function
#'   automatically selects a standard vessel from the available data.
#' @description
#' Standardizes fishing effort by adjusting for differences in fishing power
#' among vessel types relative to a selected standard vessel. The method
#' estimates the relative fishing power of each vessel type from observed
#' CPUE values and converts raw fishing effort into a common-efficiency effort
#' scale.
#'
#' Relative fishing power is calculated as:
#'
#' \deqn{
#' PA(i)=\frac{CPUE(i)}
#' {CPUE(standard)}
#' }
#'
#' where \eqn{CPUE(i)} is the catch-per-unit-effort of vessel type \eqn{i}
#' and \eqn{CPUE(standard)} is the CPUE of the selected standard vessel.
#'
#' Standardized effort is then computed as:
#'
#' \deqn{
#' E_{std}
#' =
#' \sum_i
#' \left[
#' PA(i)\times N(i)\times d(i)
#' \right]
#' }
#'
#' where \eqn{N(i)} is the number of boats and \eqn{d(i)} is the average
#' number of fishing days for vessel type \eqn{i}.
#' @return
#' A data frame containing year-wise standardized effort (used as standardized
#' CPUE index) along with a printed summary table and a dual-axis trend plot
#' showing standardized CPUE and mean relative fishing power by year.
#'
#' @note
#' CPUE should not contain zero values, if still zero values are present
#' then all are replaced with minimum CPUE value.
#' @references
#'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.
#'
#' Robson, D.S. (1966).
#' Estimation of the relative fishing power of individual ships.
#' ICNAF Research Bulletin, 3, 5-14.
#'
#' 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.
#'
#' @export
#'
#' @examples
#' \dontrun{
#' library(FESta)
#' data("REstd_dataset")
#' result<-REstd(data = REstd_dataset, year_col = "Year", vessel_col = "Vessel",
#'       boats_col = "Boats", days_col = "Days",
#'       cpue_col = "CPUE", standard_vessel = 'A')
#'       print(result)
#' }
REstd <- function(
    data,
    year_col,
    vessel_col,
    boats_col,
    days_col,
    cpue_col,
    standard_vessel
) {

  # -----------------------------------
  # Validation
  # -----------------------------------
  required_cols <- c(
    year_col,
    vessel_col,
    boats_col,
    days_col,
    cpue_col
  )

  missing_cols <- setdiff(required_cols, names(data))

  if (length(missing_cols) > 0) {
    stop(
      paste("Missing columns:", paste(missing_cols, collapse = ", "))
    )
  }


  # =============================================
  # ZERO HANDLING FOR CPUE
  # =============================================
  # Remove rows with CPUE <= 0 (zero or negative)
  zero_cpue_rows <- which(data[[cpue_col]] <= 0 | is.na(data[[cpue_col]]))

  if (length(zero_cpue_rows) > 0) {
    note_msg <- sprintf(
      "NOTE: Replaced %d row(s) containing CPUE <= 0 (zero or negative values)
      with minimum CPUE value.",
      length(zero_cpue_rows)
    )
    message(note_msg)
    # Find minimum positive CPUE
    min_positive <- min(data[[cpue_col]][data[[cpue_col]] > 0], na.rm = TRUE)

    # Replace only the zero rows
    data[[cpue_col]][zero_cpue_rows] <- min_positive
  }
  # -----------------------------------
  # Choose standard vessel
  # -----------------------------------
  if (is.null(standard_vessel)) {
    standard_vessel <- data[[vessel_col]][1]
  }

  cat("Standard Vessel Used:", standard_vessel, "\n\n")

  # -----------------------------------
  # Extract standard vessel CPUE
  # -----------------------------------
  standard_cpue <- data[
    data[[vessel_col]] == standard_vessel,
    cpue_col
  ][1]

  if (length(standard_cpue) == 0 || is.na(standard_cpue)) {
    stop("Standard vessel not found.")
  }

  # -----------------------------------
  # Relative Fishing Power
  # -----------------------------------
  data$Fishing_Power_PA <- data[[cpue_col]] / standard_cpue
  data$Fishing_Power_PA[data[[vessel_col]] == standard_vessel] <- 1

  # -----------------------------------
  # Raw Effort
  # -----------------------------------
  data$Raw_Effort <- data[[boats_col]] * data[[days_col]]

  # -----------------------------------
  # Standardized Effort
  # -----------------------------------
  data$Standardized_Effort <- data$Fishing_Power_PA * data$Raw_Effort

  # -----------------------------------
  # Total Standardized Effort
  # -----------------------------------
  total_standardized_effort <- sum(data$Standardized_Effort)

  # -----------------------------------
  # Summary Table
  # -----------------------------------
  summary_table <- data.frame(
    Year                 = data[[year_col]],
    Vessel_Type          = data[[vessel_col]],
    Number_of_Boats      = data[[boats_col]],
    Average_Fishing_Days = data[[days_col]],
    CPUE                 = round(data[[cpue_col]], 4),
    Relative_Fishing_Power = round(data$Fishing_Power_PA, 4),
    Raw_Effort           = round(data$Raw_Effort, 4),
    Standardized_Effort  = round(data$Standardized_Effort, 4)
  )

  # -----------------------------------
  # Year-wise aggregation for plot
  # -----------------------------------
  yearly_std_cpue <- stats::aggregate(
    Standardized_Effort ~ Year,
    data = summary_table,
    FUN  = sum
  )
  names(yearly_std_cpue) <- c("Year", "Standardized_CPUE")

  yearly_fp <- stats::aggregate(
    Relative_Fishing_Power ~ Year,
    data = summary_table,
    FUN  = mean
  )
  names(yearly_fp) <- c("Year", "Mean_Fishing_Power")

  plot_df <- merge(yearly_std_cpue, yearly_fp, by = "Year")
  plot_df  <- plot_df[order(plot_df$Year), ]

  # -----------------------------------
  # Needed summary to return
  # -----------------------------------
  needed_summary <- cbind(plot_df[, c("Year")],plot_df$Mean_Fishing_Power, plot_df[, c("Standardized_CPUE")])

  # -----------------------------------
  # Colours — matched to calculate_and_plot_cpue
  # -----------------------------------
  col_std <- "skyblue3"   # Standardized CPUE  (same as in reference function)
  col_fp  <- "#6a4c93"    # Fishing Power       (same red as Nominal CPUE line)

  # -----------------------------------
  # Year-wise trend plot
  # -----------------------------------
  plot_yearwise <- function(plot_df) {

    yr <- plot_df$Year
    sc <- plot_df$Standardized_CPUE

    all_yr <- sort(unique(as.integer(yr)))

    # --- Dynamic y-axis: breaks scale with the data, then set the
    # limit just above the last break so nothing floats in empty space ---
    data_max   <- max(sc, na.rm = TRUE)
    y_breaks   <- scales::extended_breaks(n = 6)(c(0, data_max * 1.05))
    y_breaks   <- y_breaks[y_breaks >= 0]              # drop any negative breaks
    y_max      <- max(y_breaks) * 1.03                 # tiny headroom above top break

    p <- ggplot2::ggplot(plot_df, ggplot2::aes(x = .data$Year)) +

      ggplot2::geom_line(
        ggplot2::aes(y = .data$Standardized_CPUE, colour = "Standardized CPUE"),
        linewidth = 1.1,
        linetype  = "solid"
      ) +
      ggplot2::geom_point(
        ggplot2::aes(y = .data$Standardized_CPUE, colour = "Standardized CPUE"),
        shape  = 21,
        fill   = "white",
        size   = 2.5,
        stroke = 1.3
      ) +

      ggplot2::scale_y_continuous(
        name   = "Standardized CPUE",
        limits = c(0, y_max),
        breaks = y_breaks,
        expand = ggplot2::expansion(mult = c(0, 0)),
        labels = scales::number_format(accuracy = 0.01)
      ) +

      ggplot2::scale_x_continuous(
        breaks = all_yr,
        labels = as.character(all_yr),
        limits = c(min(all_yr) - 0.5, max(all_yr) + 0.5),
        expand = ggplot2::expansion(mult = c(0, 0))
      ) +

      ggplot2::scale_colour_manual(
        name   = NULL,
        values = c("Standardized CPUE" = col_std),
        guide = ggplot2::guide_legend(
          override.aes = list(linetype = "solid", shape = 21, fill = "white")
        )
      ) +

      ggplot2::labs(title = "Standardized CPUE by Year", x = "Year") +

      ggplot2::theme_minimal(base_size = 11) +
      ggplot2::theme(
        plot.title       = ggplot2::element_text(face = "bold", hjust = 0.5, size = 12),
        axis.title.y     = ggplot2::element_text(colour = col_std, face = "bold", size = 10),
        axis.text.y      = ggplot2::element_text(colour = "black", size = 9),
        axis.text.x      = ggplot2::element_text(size = 9, angle = 45, hjust = 1),
        axis.title.x     = ggplot2::element_text(size = 10),
        legend.position  = "bottom",
        legend.key.width = grid::unit(1.4, "cm"),
        legend.text      = ggplot2::element_text(size = 9),
        panel.grid.minor = ggplot2::element_blank(),
        panel.grid.major = ggplot2::element_line(colour = "grey92", linewidth = 0.35),
        plot.margin      = ggplot2::margin(8, 8, 6, 8)
      )

    print(p)
    invisible(p)
  }

  plot_yearwise(plot_df)
  # -----------------------------------
  # Print summary table
  # -----------------------------------
  cat("\n", strrep("=", 65), "\n", sep = "")
  cat("  CPUE Summary Table\n")
  cat(strrep("=", 65), "\n\n", sep = "")
  # print(needed_summary, row.names = FALSE, digits = 4L)
  # cat("\n", strrep("=", 65), "\n\n", sep = "")
  colnames(needed_summary)<-c("Year","Mean_Fishing_Power","Standardized_CPUE")
  return(needed_summary)
}