R/GLMMstd.R

In FESta: Fishing Effort Standardization

#' Generalized Linear Mixed Model Based Standardization (GLMMstd)
#'
#' @importFrom stats predict
#' @importFrom stats aggregate
#' @importFrom stats as.formula

#' @importFrom lme4 glmer

#' @description
#' Fisheries catch and effort data frequently exhibit dependency structures
#' arising from repeated observations of vessels, trips, areas, observers, or
#' other sampling units. Generalized Linear Mixed Models (GLMMs) accommodate
#' these dependencies by incorporating random effects in addition to fixed
#' explanatory variables.
#'
#' In this implementation, CPUE is calculated as:
#'
#' \deqn{
#' CPUE_i=\frac{Catch_i}{Effort_i}
#' }
#'
#' and transformed using a logarithmic transformation:
#'
#' \deqn{
#' log(CPUE_i)
#' }
#'
#' The model fitted is:
#' \deqn{g(\mu_i) = X_i\beta + Z_i u}
#' where \eqn{g} is the log link, \eqn{\mu_i = E(CPUE_i)},
#' and CPUE follows a Gamma distribution.
#'
#' Fixed effects typically include factors such as year, gear, season, or
#' fishing area, while random effects may represent vessels, trips, observers,
#' ports, or other grouping variables.
#'
#' Standardized CPUE indices are obtained by predicting CPUE for each level of
#' the selected index variable while averaging over random-effect variation.
#'
#'
#' @param data A data frame containing the columns of year, catch, effort,
#' fixed effects and random effects.
#' See the example dataset \link[=GLMMstd_dataset]{GLMMstd_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 random_effects Specify the column names of the random effects
#' in vector format (eg. c("Vessel","Area")). Should not contain any fixed
#' effect column name.
#'@param log_transform Specify TRUE or FALSE. By default set TRUE.
#'It ensures CPUE values will be shown after required
#'log transformation of data.
#'
#'@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 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.
#'
#' @references
#' Pinheiro, J.C., and Bates, D.M. (2000).
#' Mixed-Effects Models in S and S-PLUS.
#' Springer-Verlag, New York.
#'
#' 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('GLMMstd_dataset')
#' result<-GLMMstd(
#'   data=GLMMstd_dataset,
#'   year_col="Year",
#'   catch_col = "Catch",
#'   effort_col = "Effort",
#'   fixed_effects = c("Year"),
#'   random_effects = c("Vessel"),
#'   log_transform = TRUE
#' )
#' print(result)
#' }
#'
#' @export
GLMMstd <- function(
    data,
    year_col,
    catch_col,
    effort_col,
    fixed_effects,
    random_effects,
    log_transform = TRUE

) {
  index_variable=year_col
  # ---- Validation ----
  required_cols <- c(year_col,catch_col, effort_col, fixed_effects, random_effects)
  missing_cols <- setdiff(required_cols, names(data))
  if (length(missing_cols) > 0) {
    stop("Missing columns: ", paste(missing_cols, collapse = ", "))
  }

  df <- data

  # ---- Ensure index_variable is in fixed_effects ----
  if (!(index_variable %in% fixed_effects)) {
    fixed_effects <- unique(c(fixed_effects, index_variable))
    message("Added '", index_variable, "' to fixed_effects.")
  }

  # =============================================
  # # ---- Handle zeros in catch and effort ----
  # =============================================

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

  # 2. Replace zero catch with minimum positive catch value
  zero_catch_rows <- which(df[[catch_col]] == 0 & !is.na(df[[catch_col]]))
  if (length(zero_catch_rows) > 0) {
    min_positive_catch <- min(df[[catch_col]][df[[catch_col]] > 0], na.rm = TRUE)
    df[[catch_col]][zero_catch_rows] <- min_positive_catch
    message(paste("Replaced", length(zero_catch_rows),
                  "zero catch values with minimum catch value:",
                  round(min_positive_catch, 4)))
  }

  # 3. Final check for any remaining non-positive values (safety check)
  if (any(df[[catch_col]] <= 0 | df[[effort_col]] <= 0, na.rm = TRUE)) {
    # Remove any remaining problematic rows
    df <- df[df[[catch_col]] > 0 & df[[effort_col]] > 0, ]
    warning("Removed additional rows with non-positive values after handling zeros")
  }

  if (nrow(df) == 0) stop("No positive observations remaining after zero handling.")
  # =============================================

  # ---- Calculate CPUE ----
  df$CPUE    <- df[[catch_col]] / df[[effort_col]]
  df$logCPUE <- log(df$CPUE)

  # ---- Convert all factors ----
  for (v in c(fixed_effects, random_effects)) {
    df[[v]] <- as.factor(df[[v]])
  }
  df <- droplevels(df)

  # ---- Build formula ----
  # log_transform = TRUE  -> response = logCPUE, gaussian() identity link
  # log_transform = FALSE -> response = CPUE,    Gamma(link="log")
  response_col   <- if (log_transform) "logCPUE" else "CPUE"
  fixed_formula  <- paste(fixed_effects, collapse = " + ")
  random_formula <- paste0("(1|", random_effects, ")", collapse = " + ")
  formula_text   <- paste(response_col, "~", fixed_formula, "+", random_formula)

  glmm_family <- if (log_transform) stats::gaussian() else stats::Gamma(link = "log")

  # ---- Try nlminbwrap first, fall back to allFit if needed ----
  #message("Optimizing model fit...")

  if (log_transform) {
    # Gaussian/identity response (logCPUE) -> use lmer(), not glmer().
    # glmer() with family=gaussian() internally redirects to lmer-style
    # fitting but is incompatible with glmerControl()'s GLMM-specific
    # arguments (tolPwrss, compDev, nAGQ0initStep, and the "nlminbwrap"
    # optimizer), causing an "unused arguments" error.
    glmm_model <- suppressWarnings(
      lme4::lmer(
        stats::as.formula(formula_text),
        data    = df,
        control = lme4::lmerControl(
          optimizer = "bobyqa",
          optCtrl   = list(maxfun = 2e5)
        )
      )
    )
  } else {
    glmm_model <- suppressWarnings(
      lme4::glmer(
        stats::as.formula(formula_text),
        data    = df,
        family  = glmm_family,
        control = lme4::glmerControl(
          optimizer = "nlminbwrap",
          optCtrl   = list(maxfun = 2e5)
        )
      )
    )
  }

  # Check if the chosen optimizer converged cleanly
  converged <- is.null(glmm_model@optinfo$conv$lme4$messages)

  if (!converged) {
    message("Initial optimizer did not converge cleanly, trying additional optimizers...")
    all_fits <- suppressWarnings(
      lme4::allFit(glmm_model, verbose = FALSE)
    )
    valid_optimizers <- if (log_transform) {
      c("bobyqa", "Nelder_Mead", "nloptwrap")   # valid for lmer()
    } else {
      c("bobyqa", "nlminbwrap", "Nelder_Mead")  # valid for glmer()
    }
    all_fits       <- all_fits[names(all_fits) %in% valid_optimizers]
    conv_status    <- sapply(all_fits, function(fit) is.null(fit@optinfo$conv$lme4$messages))
    converged_fits <- all_fits[conv_status]

    if (length(converged_fits) == 0) {
      warning("No optimizer converged cleanly; using initial fit.")
    } else {
      logliks    <- sapply(converged_fits, stats::logLik)
      best_name  <- names(which.max(logliks))
      message("Best optimizer: ", best_name)
      glmm_model <- converged_fits[[best_name]]
    }
  } else {
    #message("Model converged successfully with nlminbwrap.")
  }
  #glmm_model@optinfo$conv$lme4$messages  # exact complaint
  # ---- Model selection criteria ----
  # cat("  Model Selection Criteria\n")
  # cat(strrep("-", 65), "\n\n", sep = "")
   model_aic <- stats::AIC(glmm_model)
   model_bic <- stats::BIC(glmm_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")

  # ---- Standardized predictions ----
  # No back-transformation: prediction is returned on whatever scale the
  # model was fit on (logCPUE when log_transform=TRUE, CPUE otherwise),
  # consistent with GAMstd's approach.
  prediction_levels <- levels(df[[index_variable]])
  results <- data.frame()

  for (lev in prediction_levels) {
    # Build newdata with reference values for all fixed effects except index_variable
    newdata <- data.frame(row.names = 1)
    for (v in fixed_effects) {
      if (v == index_variable) {
        newdata[[v]] <- factor(lev, levels = levels(df[[v]]))
      } else {
        newdata[[v]] <- factor(levels(df[[v]])[1], levels = levels(df[[v]]))
      }
    }
    # Predict population-level (re.form = NA)
    pred <- predict(glmm_model, newdata = newdata, re.form = NA, type = "response",
                    allow.new.levels = TRUE)
    if (log_transform==TRUE) {
      pred <- exp(pred)
    }
    results <- rbind(results, data.frame(
      Index_Level = lev,
      Standardized_CPUE = pred
    ))
  }
  row.names(results) <- NULL

  # ---- Aggregate total catch by index variable ----
  catch_by_index <- stats::aggregate(
    df[[catch_col]] ~ df[[index_variable]],
    FUN = sum
  )
  names(catch_by_index) <- c("Index_Level", "Total_Catch")
  catch_by_index$Index_Level <- as.character(catch_by_index$Index_Level)

  # ---- Prepare plot data ----
  plot_df <- results
  plot_df$Index_Level <- as.character(plot_df$Index_Level)
  plot_df <- merge(plot_df, catch_by_index, by = "Index_Level", all.x = TRUE)
  if (any(is.na(plot_df$Total_Catch))) {
    warning("Some index levels have no catch data; setting Total_Catch to 0.")
    plot_df$Total_Catch[is.na(plot_df$Total_Catch)] <- 0
  }

  # Convert Index_Level to numeric if possible (e.g., Year), else use position
  numeric_levels <- suppressWarnings(as.numeric(plot_df$Index_Level))

  if (all(!is.na(numeric_levels))) {
    plot_df$Year <- numeric_levels
    x_labels     <- NULL
  } else {
    plot_df$Year <- seq_len(nrow(plot_df))
    x_labels     <- plot_df$Index_Level
  }

  # ---- Plot ----
  # plot_cpue_index(
  #   year        = plot_df$Year,
  #   Total_Catch = plot_df$Total_Catch,
  #   CPUE        = plot_df$Standardized_CPUE
  #   # x_labels    = x_labels
  # )
  calculate_and_plot_cpue(year = data[[year_col]],
                          std_cpue =plot_df$Standardized_CPUE,effort = data[[effort_col]], total_catch = NULL,catch = data[[catch_col]],
                          log_transform = log_transform ,aic=NULL,bic=NULL)

  # ---- Return ----
  # return(list(
  #   model              = glmm_model,
  #   formula_used       = formula_text,
  #   standardized_index = results
  # ))
}