R/ORQGLMstd.R

In FESta: Fishing Effort Standardization

# =========================================================
# Global variable declarations for R CMD check
# =========================================================
utils::globalVariables(c("orq_cpue", "Total_Catch",
                         "Standardized_CPUE", "Year", "CPUE"))


#' Ordered Quantile Transformation GLM Based Standardization (ORQGLMstd)
#'
#'
#' @importFrom stats glm predict aggregate as.formula residuals gaussian
#' @importFrom stats glm.control model.matrix shapiro.test

#' @importFrom bestNormalize orderNorm
#'
#' @param data A data frame containing the columns of year, catch, effort
#'  and fixed effects.
#'  See the example dataset \link[=ORQGLMstd_dataset]{ORQGLMstd_dataset}.
#'@param year_col Specify the year column name (eg. "Year").
#' @param catch_col Specify the catch column name (eg. "Catch").
#' @param effort_col Specify the effort column name (eg. "Effort").
#' @param fixed_effects Specify the column names of the fixed effects in vector
#' format (eg. c("Year","Gear")).

#' @param maxit Maximum number of iterations allowed during GLM fitting.
#'   Default \code{100}.
#'
#' @description
#' Fisheries catch and CPUE data are often highly skewed, heavy-tailed, or
#' non-normal, which can violate the assumptions of conventional Gaussian
#' models. The Ordered Quantile (ORQ) transformation addresses this by
#' mapping observed CPUE values to an approximately standard normal
#' distribution while preserving rank order.
#'
#' After transformation, a Gaussian GLM is fitted using the specified
#' explanatory variables. Standardized predictions are generated for each
#' level of the index variable and back-transformed to the original CPUE
#' scale using the inverse ORQ transformation.
#'
#' @return
#' The output includes AIC and SBC/BIC values, 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 corresponding rows will be
#'removed for doing the CPUE calculation.
#' @references
#' Peterson, R.A. (2021).
#' Finding Optimal Normalizing Transformations via bestNormalize.
#' The R Journal, 13(1), 310-329.
#'
#' 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.
#'
#' @export
#'
#' @examples
#' \dontrun{
#' library(FESta)
#' data("ORQGLMstd_dataset")
#' result<-ORQGLMstd(
#'   data  = ORQGLMstd_dataset,
#'   year_col="Year",
#'   catch_col      = "Catch",
#'   effort_col     = "Effort",
#'   fixed_effects  = c("Year", "Gear"),
#'   maxit          = 100
#' )
#'print(result)
#' }
ORQGLMstd <- function(
    data,
    year_col,
    catch_col,
    effort_col,
    fixed_effects,

    maxit = 100
) {
  index_variable = year_col
  # -------------------------------------------------------
  # 1. Validate columns
  # -------------------------------------------------------
  required_cols <- c(catch_col, effort_col, fixed_effects,year_col)
  missing_cols  <- setdiff(required_cols, names(data))
  model_effects <- union(fixed_effects, index_variable)
  if (length(missing_cols) > 0)
    stop(paste("Missing columns:", paste(missing_cols, collapse = ", ")))

  if (!requireNamespace("bestNormalize", quietly = TRUE))
    stop("Package 'bestNormalize' is required. Install with: install.packages('bestNormalize')")

  if (!(index_variable %in% model_effects))
    stop(sprintf("'index_variable' (\"%s\") must be included in 'model_effects'.", index_variable))
  # -------------------------------------------------------
  # 2. Prepare data
  # -------------------------------------------------------
  dff <- data

  # Remove rows with invalid effort
  n_bad <- sum(dff[[effort_col]] <= 0, na.rm = TRUE)
  if (n_bad > 0) {
    dff <- dff[dff[[effort_col]] > 0, ]
    message(sprintf("[ORQGLMstd] %d row(s) with effort <= 0 removed.", n_bad))
  }
  if (nrow(dff) == 0L)
    stop("No valid rows remain after removing zero/negative effort.")

  # FIX 1: Always coerce fixed effects to factor, regardless of input type.
  # The original code only converted character columns, leaving integer/numeric
  # Year columns as continuous — which treats Year as a slope, not levels.
  for (v in model_effects)
    dff[[v]] <- as.factor(dff[[v]])

  dff <- droplevels(dff)

  # -------------------------------------------------------
  # 3. ORQ transformation of CPUE
  # -------------------------------------------------------
  raw_cpue <- dff[[catch_col]] / dff[[effort_col]]

  orq_obj      <- suppressWarnings(bestNormalize::orderNorm(raw_cpue))
  dff$orq_cpue <- predict(orq_obj)   # z-scores, approximately N(0,1)

  #message("[ORQGLMstd] ORQ transformation applied to CPUE.")

  # -------------------------------------------------------
  # 4. Gaussian GLM on transformed CPUE
  # No effort offset needed: response is already a rate (CPUE),
  # not raw catch.
  # -------------------------------------------------------
  #smooth_terms  <- paste(sprintf("s(%s)", smooth_variable), collapse = " + ")
  fixed_formula <- paste(model_effects, collapse = " + ")
  formula_text  <- paste("orq_cpue ~", fixed_formula)

  glm_model <- stats::glm(
    stats::as.formula(formula_text),
    data    = dff,
    family  = stats::gaussian(link = "identity"),
    control = stats::glm.control(maxit = maxit),
    x = FALSE, y = FALSE, model = FALSE
  )

  #message("[ORQGLMstd] Gaussian GLM fitted on ORQ-transformed CPUE.")
  # cat("  Model Selection Criteria\n")
  # cat(strrep("-", 65), "\n\n", sep = "")
  model_aic <- stats::AIC(glm_model)
  model_bic <- stats::BIC(glm_model)
  # lm=list(AIC=model_aic,BIC=model_bic)
  # print(lm)
  # cat("\n")
  aic_bic_table <- data.frame(
    AIC = round(model_aic, 4L),
    BIC = round(model_bic, 4L)
  )
  cat("  Model Selection Criteria\n")
  cat(strrep("-", 35), "\n\n", sep = "")
  print(aic_bic_table, row.names = FALSE)
  cat("\n")
  # -------------------------------------------------------
  # 5. Prediction scaffold
  # -------------------------------------------------------
  pred_lv  <- levels(dff[[index_variable]])
  ref_vals <- setNames(
    lapply(model_effects, function(v) levels(dff[[v]])[1]),
    model_effects
  )

  orq_min <- min(dff$orq_cpue, na.rm = TRUE)
  orq_max <- max(dff$orq_cpue, na.rm = TRUE)

  results <- data.frame()

  for (lv in pred_lv) {

    # Build single-row newdata at reference levels
    nd <- as.data.frame(matrix(nrow = 1L, ncol = 0L))
    for (v in model_effects)
      nd[[v]] <- if (v == index_variable) lv else ref_vals[[v]]
    for (v in model_effects)
      nd[[v]] <- factor(nd[[v]], levels = levels(dff[[v]]))

    # Predicted ORQ z-score for this index level
    pred_orq <- suppressWarnings(
      as.numeric(stats::predict(glm_model, newdata = nd, type = "response")[1L])
    )

    # Clamp to observed ORQ range to prevent extrapolation instability
    pred_orq <- max(orq_min, min(orq_max, pred_orq))

    # Back-transform z-score → original CPUE scale
    pred_cpue <- suppressWarnings(
      stats::predict(orq_obj, newdata = pred_orq, inverse = TRUE)
    )
    if (is.na(pred_cpue)) pred_cpue <- NA_real_

    results <- rbind(results,
                     data.frame(Index_Level = as.character(lv),
                                CPUE_backtransformed = round(pred_cpue, 6)))
  }
  row.names(results) <- NULL

  # -------------------------------------------------------
  # FIX 2: Compute Standardized CPUE as index relative to mean.
  # Previously Standardized_CPUE was identical to Predicted_CPUE —
  # the normalization step was missing entirely.
  # Dividing by the mean gives an interpretable index:
  #   > 1  →  above-average abundance year
  #   < 1  →  below-average abundance year
  #   = 1  →  average year
  # -------------------------------------------------------
  mean_cpue            <- mean(results$CPUE_backtransformed, na.rm = TRUE)
  results$Standardized_CPUE <- round(results$CPUE_backtransformed / mean_cpue, 4)

  # -------------------------------------------------------
  # 6. Residual normality check (Shapiro-Wilk)
  # -------------------------------------------------------
  shapiro_result <- tryCatch(
    stats::shapiro.test(stats::residuals(glm_model)),
    error = function(e) NULL
  )

  # if (!is.null(shapiro_result)) {
  #   cat("\n", strrep("-", 50), "\n", sep = "")
  #   cat("  Shapiro-Wilk Normality Test on GLM Residuals\n")
  #   cat(strrep("-", 50), "\n", sep = "")
  #   cat(sprintf("  W = %.4f,  p-value = %.4f\n",
  #               shapiro_result$statistic, shapiro_result$p.value))
  #   if (shapiro_result$p.value >= 0.05)
  #     cat("  Residuals appear approximately normal (p >= 0.05).\n")
  #   else
  #     cat("  Residuals deviate from normality (p < 0.05).\n")
  #   cat(strrep("-", 50), "\n\n", sep = "")
  #}
  calculate_and_plot_cpue(year = data[[year_col]],
                           std_cpue =results$Standardized_CPUE,effort = data[[effort_col]], total_catch = NULL,log_transform = NULL,catch = data[[catch_col]],aic=NULL,bic=NULL )

}