R/MGMSstd.R

In FESta: Fishing Effort Standardization

utils::globalVariables(
  c(
    'Species', 'Value', 'Mean_CPUE','Total_standardized_CPUE'
  )
)
#' Multi-Gear Mean Standardization (MGMSstd)
#'
#'
#' @param data A data frame containing year, gear, species, and CPUE columns.
#' @param year_col Specify the year column name (eg. "Year").
#' @param gear_col Specify the gear types column name (eg. "Gears").
#' @param species_col Specify the species column name (eg. "Species").
#' @param cpue_col Specify the CPUE column name (eg. "CPUE_value").
#'
#'@description
#' The Multigear Mean Standardization (MGMS) method was proposed by
#' Gibson-Reinemer et al. (2017) to combine CPUE data collected using
#' different sampling gears into a common relative scale suitable for
#' community analyses.
#' Initially, CPUE values are expressed as relative abundance:
#'
#' \deqn{
#' RA_{ij}=
#' \frac{c_{ij}/e}
#' {TC_j/e}
#' }
#'
#' where
#'
#' \itemize{
#'   \item \eqn{c_{ij}}{cij} is the catch of species \eqn{i} in year (or sample) \eqn{j}.
#'   \item \eqn{e} is the sampling effort.
#'   \item \eqn{TC_j} is the total catch of all species in year (or sample) \eqn{j}.
#' }
#'
#' To preserve both within-sample and among-sample abundance patterns,
#' the total CPUE of each sample is standardized by the mean total
#' CPUE across all samples:
#'
#' \deqn{
#' MSC_{ij}
#' =
#' \frac{c_{ij}/e}{TC_j/e}
#' \times
#' \frac{TC_j/e}{\overline{TC}/e}
#' }
#'
#' which simplifies to
#'
#' \deqn{
#' MSC_{ij}
#' =
#' \frac{c_{ij}/e}
#' {\overline{TC}/e}
#' }
#'
#' where
#'
#' \deqn{
#' \overline{TC}/e
#' }
#'
#' is the mean total catch per unit effort.
#'
#'
#' @return
#' The function produces:
#' \itemize{
#'   \item Gear-specific standardized CPUE tables.
#'   \item A multiline plot of total standardized CPUE by gear.
#' }
#' @note
#' CPUE should not contain zero values, if still zero values are present
#' then all are replaced with minimum CPUE value.
#' @references
#' Gibson-Reinemer, D. K., Ickes, B. S., & Chick, J. H. (2017).
#' Development and assessment of a new method for combining
#' catch per unit effort data from different fish sampling gears:
#' multigear mean standardization (MGMS).
#' Canadian Journal of Fisheries and Aquatic Sciences,
#' 74(1), 8--14.
#' https://doi.org/10.1139/cjfas-2016-0003
#'
#' 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.
#'
#' \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("MGMSstd_dataset")
#'
#' result <- MGMSstd(
#'   data = MGMSstd_dataset,
#'   year_col = "Year",
#'   gear_col = "Gear",
#'   species_col = "Species",
#'   cpue_col = "CPUE"
#' )
#'print(result)
#' }
#'
#' @export
#'
MGMSstd <- function(data, year_col, gear_col, species_col, cpue_col){
  value_col <- cpue_col
  # =============================================
  # ZERO HANDLING FOR CPUE (Same logic as REstd)
  # =============================================
  # Find rows with CPUE <= 0 or NA
  zero_cpue_rows <- which(data[[cpue_col]] == 0 | is.na(data[[cpue_col]]))

  if (length(zero_cpue_rows) > 0) {
    # Find minimum POSITIVE CPUE value (excluding zeros)
    min_cpue <- min(data[[cpue_col]][data[[cpue_col]] != 0], na.rm = TRUE)

    if (is.finite(min_cpue)) {
      # Replace only the zero rows with minimum positive value
      data[[cpue_col]][zero_cpue_rows] <- min_cpue

      note_msg <- sprintf(
        "NOTE: Replaced %d CPUE column values with minimum CPUE (%.6g)",
        length(zero_cpue_rows), min_cpue
      )
      message(note_msg)
    }
  }
  # =============================================
  species_list <- unique(data[[species_col]])
  result_list <- lapply(species_list, function(sp) {
    df_sp <- data[data[[species_col]] == sp, c(year_col, gear_col, value_col)]
    names(df_sp) <- c("Year", "Gear", "Value")
    # Year x Gear matrix for this species
    wide_df <- tidyr::pivot_wider(
      df_sp,
      names_from  = Gear,
      values_from = Value,
      values_fill = 0
    )
    wide_df <- wide_df[order(wide_df$Year), ]
    yrs <- wide_df$Year
    mat <- as.matrix(wide_df[, -1, drop = FALSE])
    rownames(mat) <- yrs
    # Row sums -> mean of row sums (single scalar for this species)
    row_sums   <- rowSums(mat, na.rm = TRUE)
    mean_value <- mean(row_sums, na.rm = TRUE)
    # Normalize the whole table by that mean
    norm_mat <- mat / mean_value
    norm_mat
  })
  names(result_list) <- species_list
  ###################PLOT
  # -----------------------------------------
  # Multiline plot of total CPUE by species
  # -----------------------------------------
  plot_data <- do.call(
    rbind,
    lapply(names(result_list), function(sp) {
      mat <- result_list[[sp]]
      data.frame(
        Year = as.numeric(rownames(mat)),
        Species = sp,
        Total_standardized_CPUE = rowSums(mat, na.rm = TRUE)
      )
    })
  )

  # -----------------------------------------
  # ROBUST Y-AXIS: pretty() breaks first, then
  # force the top break strictly above the max
  # data value so no point ever sits above the
  # last labelled gridline
  # -----------------------------------------
  cpue_max <- max(plot_data$Total_standardized_CPUE, na.rm = TRUE)
  cpue_min <- min(plot_data$Total_standardized_CPUE, na.rm = TRUE)

  cpue_breaks <- pretty(c(cpue_min, cpue_max), n = 6)
  step <- diff(cpue_breaks)[1]

  y_bot <- min(cpue_breaks)
  y_top <- max(cpue_breaks)
  if (y_top <= cpue_max) y_top <- y_top + step   # guarantee headroom

  cpue_breaks <- cpue_breaks[cpue_breaks >= y_bot & cpue_breaks <= y_top]
  if (max(cpue_breaks) < y_top) cpue_breaks <- c(cpue_breaks, y_top)

  p <- ggplot2::ggplot(
    plot_data,
    ggplot2::aes(
      x = Year,
      y = Total_standardized_CPUE,
      colour = Species,
      group = Species
    )
  ) +
    ggplot2::geom_line(linewidth = 1.1) +
    ggplot2::geom_point(size = 2.5) +
    ggplot2::scale_x_continuous(
      breaks = sort(unique(plot_data$Year))
    ) +
    ggplot2::scale_y_continuous(
      limits = c(y_bot, y_top),
      breaks = cpue_breaks,
      expand = ggplot2::expansion(mult = c(0, 0))
    ) +
    ggplot2::labs(
      title = "Total Standardized CPUE by Species",
      x = "Year",
      y = "Total Standardized CPUE",
      colour = "Species"
    ) +
    ggplot2::theme_minimal(base_size = 11) +
    ggplot2::theme(
      plot.title = ggplot2::element_text(
        hjust = 0.5,
        face = "bold"
      ),
      legend.position = "bottom"
    )
  print(p)
  #########################################
  # ============================================
  # COMBINE ALL SPECIES RESULTS INTO ONE DATA FRAME
  # ============================================

  result <- do.call(
    rbind,
    lapply(names(result_list), function(sp) {

      mat <- result_list[[sp]]

      df <- as.data.frame(mat)

      # Add Year and Species
      df$Year <- as.numeric(rownames(mat))
      df$Species <- sp

      # Total standardized CPUE
      df$Total_standardized_CPUE <- rowSums(mat, na.rm = TRUE)

      # Put columns in desired order
      df <- df[, c(
        "Species",
        "Year",
        setdiff(
          names(df),
          c("Species", "Year", "Total_standardized_CPUE")
        ),
        "Total_standardized_CPUE"
      )]

      df
    })
  )

  # Reset row names
  rownames(result) <- NULL

  return(result)
}
# ---- Custom print method: only fires on auto-print (not on assignment) ----
#'@export
print.cpue_species_tables <- function(x, ...) {
  for (sp in names(x)) {
    cat("\n", strrep("=", 60), "\n", sep = "")
    cat("  Species:", sp, "\n")
    cat(strrep("=", 60), "\n\n", sep = "")
    tbl <- round(x[[sp]], 4)
    total_standardized_CPUE <- round(rowSums(x[[sp]], na.rm = TRUE), 4)
    out <- data.frame(
      Year = rownames(tbl),
      tbl,
      Total_standardized_CPUE = total_standardized_CPUE,
      row.names = NULL,
      check.names = FALSE
    )
    print(out, ...)
  }
  invisible(x)
}