R/GAMstd.R

In FESta: Fishing Effort Standardization

#library(nlme)
#' Generalized Additive Model Based Standardization (GAMstd)

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


#' @importFrom stats glm.control
#' @importFrom stats model.matrix

#' @importFrom stats median
#' @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, catch, effort,
#' fixed effects and smooth terms.
#' See the example dataset \link[=GAMstd_dataset]{GAMstd_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 smooth_terms Specify the column names of smooth terms
#' in vector format (eg. c("SST","Depth")).
#' @param k Integer value for smooth terms. If \code{NULL}, an appropriate
#' value is selected automatically.
#'@param log_transform Specify TRUE or FALSE. By default set TRUE.
#'It ensures CPUE values will be shown after required
#'log transformation of data.
#' @description
#' Generalized Additive Models (GAMs) provide a flexible approach for fisheries
#' CPUE standardization by combining linear effects of categorical variables
#' with smooth, non-parametric effects of continuous environmental covariates.
#'
#' In this implementation, categorical variables such as year, gear, area, or
#' season are incorporated as fixed effects, while continuous variables such as
#' sea surface temperature (SST), depth, salinity, or other environmental
#' covariates are modeled using spline-based smooth functions.
#'
#' The fitted GAM may be expressed as:
#'
#' \deqn{
#' g(\mu_i)
#' =
#' \beta_0
#' +
#' \sum_{k=1}^{p}\beta_kX_{ik}
#' +
#' \sum_{j=1}^{q}f_j(Z_{ij})
#' +
#' \log(E_i)
#' }
#'
#' where \eqn{g(.)} is the link function,
#' \eqn{X_{ik}} represents categorical predictors,
#' \eqn{f_j(.)} are smooth functions of continuous covariates,
#' \eqn{E_i} denotes fishing effort,
#' and \eqn{\mu_i} is the expected catch.
#'
#' The offset term adjusts the expected catch for
#' differences in fishing effort, allowing predictions to be
#' standardized to a common unit of effort resulting in standardized predictions
#' on a common unit-effort basis.
#'@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.
#'
#' @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.

#'
#' @references
#' Hastie, T., & Tibshirani, R. (1986). Generalized additive models.
#' Statistical science, 1(3), 297-310.
#'
#' 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("GAMstd_dataset")
#' result<-GAMstd(data=GAMstd_dataset,year_col='Year',catch_col='Catch',
#' effort_col='Effort',fixed_effects = c("Year", "Gear"),
#' smooth_terms = c("SST", "Depth"), k=10,log_transform = TRUE)
#' print(result)
#' }
GAMstd <- function(
    data,
    year_col,
    catch_col,
    effort_col,
    fixed_effects,
    smooth_terms,
    k = NULL,
    log_transform = TRUE

) {
  k_fixed=k
  if(is.null(k_fixed) == TRUE){
    k_adaptive = TRUE}else{
      k_adaptive = FALSE
    }

  show_table = TRUE
  draw_plot = TRUE

  index_variable=year_col
  # ===========================================================
  # 1. VALIDATION
  # ===========================================================
  required_cols <- c(
    year_col,
    catch_col,
    effort_col,
    fixed_effects,
    smooth_terms,
    index_variable
  )

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

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


  # ===========================================================
  # 2. DATA PREPARATION
  # ===========================================================
  df <- data

  # Convert fixed effects to factors
  for (v in fixed_effects) {
    df[[v]] <- as.factor(df[[v]])
  }

  # Effort must be strictly positive — offset(log(Effort)) and the CPUE ratio
  # itself are undefined at zero/negative effort, so those rows are genuinely
  # unusable and must be dropped.
  n_removed <- sum(df[[effort_col]] <= 0, na.rm = TRUE)
  if (n_removed > 0) {
    message("[GAMstd] Removing ", n_removed, " rows with non-positive effort.")
    df <- df[df[[effort_col]] > 0, ]
  }

  # Catch = 0 (or negative, e.g. a data-entry artifact) is a real observation,
  # not an invalid one — dropping it discards information. Instead, substitute
  # the smallest positive observed catch in that column, which keeps the row
  # usable under log_transform = TRUE (log(Catch)) or Gamma-log-link.
  n_zero_catch <- sum(df[[catch_col]] <= 0, na.rm = TRUE)
  if (n_zero_catch > 0) {
    min_positive_catch <- min(df[[catch_col]][df[[catch_col]] > 0], na.rm = TRUE)
    message("[GAMstd] Replacing ", n_zero_catch,
            " non-positive catch value(s) with minimum positive catch (",
            round(min_positive_catch, 4), ").")
    df[[catch_col]][df[[catch_col]] <= 0] <- min_positive_catch
  }

  if (nrow(df) == 0) {
    stop("No valid observations remain after removing non-positive effort.")
  }
  # ===========================================================
  # 3. ADAPTIVE k SELECTION (per smooth term)
  # ===========================================================
  n_obs    <- nrow(df)
  n_smooth <- length(smooth_terms)

  # Global upper bound from data size
  if (k_adaptive) {
    global_k <- max(3L, min(10L, floor(n_obs / (n_smooth * 3L))))
    message("[GAMstd] Data-size-based k upper bound = ", global_k)
  } else {
    global_k <- k_fixed
    #message("[GAMstd] Fixed k = ", global_k)
  }

  # Per-term k: must not exceed (unique covariate values - 1)
  k_per_term <- vapply(smooth_terms, function(v) {
    n_unique <- length(unique(df[[v]]))
    k_cap    <- n_unique - 1L
    k_final  <- min(global_k, k_cap)
    if (k_final < global_k) {
      message("[GAMstd] k for '", v, "' capped at ", k_final,
              " (only ", n_unique, " unique values)")
    }
    as.integer(max(3L, k_final))   # GAM needs at least 3 basis functions
  }, integer(1L))

  # ===========================================================
  # 4. BUILD FORMULA (per-term k)
  # ===========================================================
  fixed_formula  <- paste(fixed_effects, collapse = " + ")

  smooth_formula <- paste0(
    "s(", smooth_terms, ", k = ", k_per_term, ")",
    collapse = " + "
  )

  # The offset is ALWAYS on the log scale — this is what turns a catch
  # model into a CPUE (rate) model, regardless of family/link below.
  # log_transform only decides whether the RESPONSE is logged, and
  # therefore which family/link is statistically appropriate:
  #   log_transform = TRUE  -> response = log(Catch), gaussian(), identity link
  #   log_transform = FALSE -> response = Catch,       Gamma(),   log link
  response_term <- if (log_transform) paste0("log(", catch_col, ")") else catch_col

  formula_text <- paste(
    response_term,
    "~",
    fixed_formula,
    "+",
    smooth_formula,
    "+ offset(log(",
    effort_col,
    "))"
  )

  #message("[GAMstd] Formula: ", formula_text)

  # ===========================================================
  # 5. FIT GAM
  # ===========================================================
  gam_family <- if (log_transform) gaussian() else Gamma(link = "log")

  gam_model <- mgcv::gam(
    formula = as.formula(formula_text),
    data = df,
    family = gam_family,
    method = "REML"  # Better for smoothness selection
  )
  # ---- Model selection criteria ----
  # cat("  Model Selection Criteria\n")
  # cat(strrep("-", 65), "\n\n", sep = "")
  model_aic <- stats::AIC(gam_model)
  model_bic <- stats::BIC(gam_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")


  # ===========================================================
  # 6. PREDICTION DATA
  # ===========================================================
  # Fixed effects other than index_variable are held at a single reference
  # level (their first factor level) so the predicted series isolates the
  # index_variable's effect. If index_variable is itself a fixed effect,
  # every level of it is kept (one prediction row per year).
  prediction_values <- list()
  for (v in fixed_effects) {
    if (v == index_variable) {
      prediction_values[[v]] <- levels(df[[v]])
    } else {
      prediction_values[[v]] <- levels(df[[v]])[1]
    }
  }
  predict_df <- expand.grid(prediction_values)

  # If index_variable isn't a fixed effect, it must be a smooth term
  # (continuous). Build one row per unique observed value of it, crossed
  # with the reference levels of the fixed effects above.
  if (!index_variable %in% fixed_effects) {
    index_values <- sort(unique(df[[index_variable]]))
    predict_df <- merge(predict_df, data.frame(x = index_values), by = NULL)
    names(predict_df)[names(predict_df) == "x"] <- index_variable
  }

  # ===========================================================
  # 7. SMOOTH TERMS AT MEAN VALUES
  # ===========================================================
  # index_variable already has its real values assigned in section 6
  # (whichever branch applied); only the *other* smooth terms get held
  # at their mean.
  for (v in smooth_terms) {
    if (v == index_variable) next
    predict_df[[v]] <- mean(df[[v]], na.rm = TRUE)
  }
  # ===========================================================
  # 8. STANDARD EFFORT
  # ===========================================================
  predict_df[[effort_col]] <- 1

  # ===========================================================
  # 9. PREDICT STANDARDIZED CPUE
  # ===========================================================
  # No back-transformation — predictions are returned on whatever scale
  # the model was fit on:
  #   log_transform = TRUE  -> prediction IS log(CPUE) directly
  #                             (since offset(log(Effort)) = log(1) = 0
  #                             and response was log(Catch))
  #   log_transform = FALSE -> prediction is CPUE directly
  #                             (Gamma/log-link already gives response scale)
  standardized_cpue <- predict(
    gam_model,
    newdata = predict_df,
    type = "response"
  )
  if (log_transform) {
    standardized_cpue <- exp(standardized_cpue)
  }

  # ===========================================================
  # 10. CREATE RESULTS
  # ===========================================================
  results <- data.frame(
    Prediction_Group = seq_len(nrow(predict_df)),
    predict_df,
    Standardized_CPUE = round(as.numeric(standardized_cpue), 4)
  )

  # Remove effort column from results (optional)
  results[[effort_col]] <- NULL

  row.names(results) <- NULL

  # ===========================================================
  # 11. AGGREGATE CATCH FOR PLOTTING
  # ===========================================================
  # For plotting, we need catch by the index variable (usually Year)
  # But if index_variable is not Year, we need to handle it properly
  calculate_and_plot_cpue(year = data[[year_col]],
                          std_cpue =standardized_cpue,effort = data[[effort_col]], total_catch = NULL,catch = data[[catch_col]],log_transform= log_transform,aic = NULL, bic = NULL)


}