R/GLMstd.R

In FESta: Fishing Effort Standardization

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

# =========================================================
# GLMstd  Multi family GLM CPUE Standardization
# =========================================================

#' Generalized Linear Model Based Standardization (GLMstd)
#'
#' @importFrom stats glm predict aggregate as.formula var sd residuals dlnorm
#' @importFrom stats gaussian poisson Gamma glm.control
#' @importFrom stats model.matrix optim median setNames
#' @importFrom ggplot2 ggplot aes geom_area geom_line geom_point
#' @importFrom ggplot2 scale_y_continuous scale_x_continuous scale_colour_manual
#' @importFrom ggplot2 expansion sec_axis guide_legend labs theme_minimal theme
#' @importFrom ggplot2 element_text element_blank element_line margin ggplotGrob
#' @importFrom scales comma number_format
#' @importFrom grid unit textGrob gpar grid.newpage grid.draw
#' @importFrom statmod tweedie
#' @importFrom MASS glm.nb
#'@importFrom tweedie tweedie_AIC
#' @param data A data frame containing the columns of year, catch, effort,
#' and fixed effects.
#' See the example dataset \link[=GLMstd_dataset]{GLMstd_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 family_type Select one or more distributions from the list:
#' \code{"gamma"},
#' \code{"tweedie"},
#' \code{"gaussian"},
#' \code{"lognormal"},
#' \code{"poisson"},
#' \code{"nbinom"}.
#'
#' @param link_function Select any of the given link functions:
#' "log", "identity", "inverse", "logit", "probit", "cloglog". If  \code{NULL}
#' then by default it will take "log".
#'
#' @param maxit Maximum number of iterations. Default \code{100}.
#'@description
#'The \code{GLMstd()} function performs catch per unit effort (CPUE)
#' standardization using generalized linear models (GLMs). The function
#' allows the user to fit one or multiple probability distributions to
#' the catch data while incorporating fishing effort as an offset term.
#'
#' Several commonly used distributions in fisheries standardization are
#' supported, including Gamma, Tweedie, Gaussian, Lognormal, Poisson,
#' and Negative Binomial distributions. The function estimates
#' standardized CPUE indices by accounting for the effects of year and
#' other explanatory variables specified as fixed effects.
#'
#' @return
#' When a single distribution is selected, the function returns a
#' summary table containing Year, Total Catch, Nominal CPUE, and
#' Standardized CPUE, together with plots of Nominal and
#' Standardized CPUE versus Total Catch.
#'
#' When multiple distributions are selected, the function provides
#' a comparative summary table and corresponding graphical displays
#' for all selected distributions.
#'
#' For each of the cases AIC and SBC/BIC values will be provided.
#'@note
#'If catch column value of has zero value(s) then they will be replaced
#'by minimum value of the catch column (in gamma and lognormal case) and further
#'if effort column has zero values then the corresponding rows will be
#'removed for doing the CPUE calculation.
#'
#' @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.
#'
#' Maunder, M.N. and Punt, A.E. (2004).
#' Standardizing catch and effort data:
#' a review of recent approaches.
#' Fisheries Research, 70, 141-159.
#'
#' Nelder, J.A. and Wedderburn, R.W.M. (1972).
#' Generalised linear models. J. R. Statist. Soc. A 137, 370-384.
#'
#' \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('GLMstd_dataset')
#' # Single family
#' result1 <- GLMstd(
#'   data          = GLMstd_dataset,
#'   year_col="Year",
#'   catch_col     = "Catch",
#'   effort_col    = "Effort",
#'   fixed_effects = c("Year", "Gear"),
#'   family_type   = "gamma"
#' )
#' print(result1)
#'
#' # Multiple families - comparison table + tiled plots
#' library(FESta)
#' data('GLMstd_dataset')
#' result2 <- GLMstd(
#'   data          = GLMstd_dataset,
#'   year_col      = "Year",
#'   catch_col     = "Catch",
#'   effort_col    = "Effort",
#'   fixed_effects = c("Year", "Gear"),
#'   family_type   = c("gamma", "lognormal", "tweedie", "poisson", "nbinom","gaussian")
#' )
#' print(result2)
#' }
GLMstd <- function(
    data,
    year_col,
    catch_col,
    effort_col ,
    fixed_effects,
    family_type,
    link_function  = NULL,
    maxit          = 100
) {

  index_variable=year_col

  # ---- FIX: consolidated logging instead of scattered message() calls.
  # Every family's data-prep/fitting notes get appended here as they happen,
  # then printed ONCE as a single grouped summary at the end, instead of
  # interleaving with R's own console output mid-loop.
  .fam_log <- new.env(parent = emptyenv())
  .fam_log$store <- list()
  .log_fam <- function(fam, msg) {
    .fam_log$store[[fam]] <- c(.fam_log$store[[fam]], msg)
  }

  # ===========================================================
  # SECTION A: Internal panel builder (reused for both plot sets)
  .build_panel <- function(year,
                           Total_Catch,
                           CPUE,
                           panel_title,
                           show_y_left  = FALSE,
                           show_y_right = FALSE,
                           show_legend  = FALSE,
                           line_colour  = "skyblue3",
                           cpue_label   = "Standardized CPUE") {

    plot_df <- data.frame(
      Year        = as.numeric(year),
      Total_Catch = Total_Catch,
      CPUE        = CPUE
    )

    # ── Dual-axis scaling ─────────────────────────────────────────────────
    # CPUE is on LEFT axis (primary).
    # Total Catch bars are scaled onto the CPUE range for display,
    # then back-transformed on the RIGHT axis.
    cpue_max  <- max(CPUE,        na.rm = TRUE)
    catch_max <- max(Total_Catch, na.rm = TRUE)

    # ---- Compute clean breaks FIRST, then derive headroom from them ----
    cpue_breaks  <- pretty(c(0, cpue_max), n = 6)
    cpue_breaks  <- cpue_breaks[cpue_breaks >= 0]
    cpue_range   <- max(cpue_breaks) * 1.15      # headroom above the top break

    catch_breaks <- pretty(c(0, catch_max), n = 6)
    catch_breaks <- catch_breaks[catch_breaks >= 0]
    catch_limit  <- max(catch_breaks) * 1.15

    catch_scale_f <- cpue_range / catch_limit
    plot_df$catch_scaled <- Total_Catch * catch_scale_f

    # ── Integer x-axis breaks ─────────────────────────────────────────────
    all_yr  <- sort(unique(as.integer(year)))
    #step    <- max(1L, ceiling(length(all_yr) / 8L))
    yr_brks <- all_yr[seq(1L, length(all_yr), by = 2)]
    #yr_brks<-all_yr
    # ── Axis title strings ────────────────────────────────────────────────
    ytitle_left  <- if (show_y_left)  "CPUE"          else ""
    ytitle_right <- if (show_y_right) "Total Catch"   else ""

    # ── Colours — matched to calculate_and_plot_cpue ──────────────────────
    col_catch <- "orange"    # Total Catch bars
    col_cpue  <- line_colour    # CPUE line colour (overridable per panel)

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

      # Total Catch bars (scaled to CPUE range, alpha 0.25)
      ggplot2::geom_bar(
        ggplot2::aes(y = catch_scaled, fill = "Total Catch"),
        stat  = "identity",
        width = 0.6,
        alpha = 0.6
      ) +

      # CPUE line + hollow circle points
      ggplot2::geom_line(
        ggplot2::aes(y = CPUE, colour = cpue_label),
        linewidth = 0.9,
        linetype  = "solid"
      ) +
      ggplot2::geom_point(
        ggplot2::aes(y = CPUE, colour = cpue_label),
        shape  = 21,
        fill   = "white",
        size   = 2.2,
        stroke = 1.2
      ) +

      # LEFT axis: CPUE  |  RIGHT axis: Total Catch (back-transform)
      ggplot2::scale_y_continuous(
        name     = ytitle_left,
        labels   = scales::number_format(accuracy = 0.01),
        expand   = ggplot2::expansion(mult = c(0, 0)),
        limits   = c(0, cpue_range),
        breaks   = cpue_breaks,                          # explicit, fixed breaks
        sec.axis = ggplot2::sec_axis(
          ~ . / catch_scale_f,
          name   = ytitle_right,
          breaks = catch_breaks,                         # in catch units, NOT multiplied
          labels = scales::comma(catch_breaks)
        )
      )  +

      # x-axis integer breaks
      ggplot2::scale_x_continuous(
        breaks = yr_brks,
        labels = as.character(yr_brks)
      ) +

      # Colour scale for CPUE line
      ggplot2::scale_colour_manual(
        name   = NULL,
        values = stats::setNames(col_cpue, cpue_label),
        guide  = if (show_legend)
          ggplot2::guide_legend(
            order = 2,
            override.aes = list(
              linetype = "solid",
              shape    = 21,
              fill     = "white"
            )
          )
        else "none"
      ) +

      # Fill scale for catch bars
      ggplot2::scale_fill_manual(
        name   = NULL,
        values = c("Total Catch" = col_catch),
        guide  = if (show_legend)
          ggplot2::guide_legend(order = 1)
        else "none"
      ) +

      ggplot2::labs(title = panel_title, x = NULL) +

      ggplot2::theme_minimal(base_size = 10) +
      ggplot2::theme(
        plot.title         = ggplot2::element_text(face  = "bold",
                                                   hjust = 0.5,
                                                   size  = 11),
        # LEFT axis — CPUE colour
        axis.title.y       = ggplot2::element_text(colour = col_cpue,
                                                   face   = "bold",
                                                   size   = 9),
        axis.text.y        = ggplot2::element_text(colour = "black",
                                                   size   = 8),
        # RIGHT axis — catch colour
        axis.title.y.right = ggplot2::element_text(colour = col_catch,
                                                   face   = "bold",
                                                   size   = 9,
                                                   angle=90),
        axis.text.y.right  = ggplot2::element_text(colour = "black",
                                                   size   = 8),
        axis.text.x        = ggplot2::element_text(size  = 8,
                                                   angle = 45,
                                                   hjust = 1),
        axis.title.x       = ggplot2::element_blank(),
        legend.position    = "bottom",
        legend.key.width   = grid::unit(1.2, "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(6, 6, 4, 6)
      )

    p
  }
  # ===========================================================
  # SECTION B: Tile n panels into one figure
  # ===========================================================
  .tile_panels <- function(panel_data_list, plot_title) {

    n     <- length(panel_data_list)
    ncols <- ceiling(sqrt(n))
    nrows <- ceiling(n / ncols)

    is_left_col  <- function(i) ((i - 1L) %% ncols) == 0L
    is_right_col <- function(i) (i %% ncols == 0L) || (i == n)

    panels <- lapply(seq_len(n), function(i) {
      pd <- panel_data_list[[i]]
      .build_panel(
        year         = pd$year,
        Total_Catch  = pd$Total_Catch,
        CPUE         = pd$CPUE,
        panel_title  = pd$title,
        show_y_left  = is_left_col(i),
        show_y_right = is_right_col(i),
        show_legend  = (i == n)
      )
    })

    if (requireNamespace("patchwork", quietly = TRUE)) {
      pw <- patchwork::wrap_plots(panels, ncol = ncols, nrow = nrows) +
        patchwork::plot_annotation(
          title   = plot_title,
          # caption = paste0(
          #   " CPUE (right axis)   ",
          #   " Species Catch (left axis)"
          # ),
          theme = ggplot2::theme(
            plot.title   = ggplot2::element_text(face  = "bold",
                                                 hjust = 0.5,
                                                 size  = 14),
            plot.caption = ggplot2::element_text(hjust  = 0.5,
                                                 size   = 9,
                                                 colour = "grey45")
          )
        )
      print(pw)
      return(invisible(pw))
    }

    if (requireNamespace("gridExtra", quietly = TRUE)) {
      grobs <- lapply(panels, ggplot2::ggplotGrob)
      ttl   <- grid::textGrob(
        plot_title,
        gp = grid::gpar(fontface = "bold", fontsize = 12)
      )
      arr <- gridExtra::arrangeGrob(
        grobs = grobs, ncol = ncols, nrow = nrows, top = ttl
      )
      grid::grid.newpage()
      grid::grid.draw(arr)
      return(invisible(arr))
    }

    op <- graphics::par(mfrow = c(nrows, ncols))
    on.exit(graphics::par(op), add = TRUE)
    for (p in panels) print(p)
    invisible(NULL)
  }

  #########New section added
  .tile_panels_grouped <- function(nominal_data, family_data_list, plot_title) {

    n_fam <- length(family_data_list)
    ncols <- ceiling(sqrt(n_fam))
    nrows <- ceiling(n_fam / ncols)

    is_left_col  <- function(i) ((i - 1L) %% ncols) == 0L
    is_right_col <- function(i) (i %% ncols == 0L) || (i == n_fam)

    nominal_plot <- .build_panel(
      year         = nominal_data$year,
      Total_Catch  = nominal_data$Total_Catch,
      CPUE         = nominal_data$CPUE,
      panel_title  = nominal_data$title,
      show_y_left  = TRUE,
      show_y_right = TRUE,
      show_legend  = TRUE,
      line_colour  = "#6a4c93",
      cpue_label   = "Nominal CPUE"
    )

    family_panels <- lapply(seq_len(n_fam), function(i) {
      pd <- family_data_list[[i]]
      .build_panel(
        year         = pd$year,
        Total_Catch  = pd$Total_Catch,
        CPUE         = pd$CPUE,
        panel_title  = pd$title,
        show_y_left  = is_left_col(i),
        show_y_right = is_right_col(i),
        show_legend  = (i == n_fam),
        cpue_label   = "Standardized CPUE"
      )
    })

    if (requireNamespace("patchwork", quietly = TRUE)) {

      # Nominal occupies ONLY the first cell of row 1 (same size as any other panel)
      design <- patchwork::area(t = 1, l = 1, b = 1, r = 1)

      idx <- 1L
      for (r in seq_len(nrows)) {
        for (c in seq_len(ncols)) {
          if (idx <= n_fam) {
            design <- c(design, patchwork::area(t = r + 1, l = c, b = r + 1, r = c))
            idx <- idx + 1L
          }
        }
      }

      pw <- patchwork::wrap_plots(c(list(nominal_plot), family_panels), design = design) +
        patchwork::plot_annotation(
          title = plot_title,
          theme = ggplot2::theme(
            plot.title = ggplot2::element_text(face = "bold", hjust = 0.5, size = 14)
          )
        )
      print(pw)
      return(invisible(pw))
    }

    # Fallback: nominal printed alone, then family grid (equal size, separate layout)
    print(nominal_plot)
    .tile_panels(family_data_list, plot_title)
  }
  # ===========================================================
  # SECTION C: Compute nominal CPUE from raw data (year-wise)
  # Formula: mean of (catch / yearly_total_catch) * yearly_total_effort
  # ===========================================================
  .compute_nominal_cpue <- function() {

    df_nom <- data.frame(
      Year   = as.integer(as.character(data[[index_variable]])),
      Catch  = as.numeric(data[[catch_col]]),
      Effort = as.numeric(data[[effort_col]])
    )

    # Yearly totals
    yr_catch  <- stats::aggregate(Catch  ~ Year, data = df_nom, FUN = sum)
    yr_effort <- stats::aggregate(Effort ~ Year, data = df_nom, FUN = sum)
    names(yr_catch)  <- c("Year", "Total_Catch")
    names(yr_effort) <- c("Year", "Total_Effort")

    # One row per year: ratio-of-totals nominal CPUE
    out <- merge(yr_catch, yr_effort, by = "Year")
    out$Nominal_CPUE <- out$Total_Catch / out$Total_Effort
    out <- out[order(out$Year), ]
    row.names(out) <- NULL
    out
  }


  # ===========================================================
  # SECTION D: Plot set 1 — single Nominal CPUE panel
  # ===========================================================
  .plot_nominal <- function(nom_df) {

    p <- .build_panel(
      year         = nom_df$Year,
      Total_Catch  = nom_df$Total_Catch,
      CPUE         = nom_df$Nominal_CPUE,
      panel_title  = " ",
      show_y_left  = TRUE,
      show_y_right = TRUE,
      show_legend  = TRUE
    )

    if (requireNamespace("patchwork", quietly = TRUE)) {
      pw <- patchwork::wrap_plots(list(p), ncol = 1L) +
        patchwork::plot_annotation(
          title   = "Nominal CPUE vs Total Catch",
          #caption = " CPUE (right axis)   Species Catch (left axis)",
          theme   = ggplot2::theme(
            plot.title   = ggplot2::element_text(face  = "bold",
                                                 hjust = 0.5,
                                                 size  = 14),
            plot.caption = ggplot2::element_text(hjust  = 0.5,
                                                 size   = 9,
                                                 colour = "grey45")
          )
        )
      print(pw)
      return(invisible(pw))
    }

    print(p)
    invisible(p)
  }

  # ===========================================================
  # Helper: Jacobian-corrected AIC/BIC for the lognormal model,
  # so it's comparable to families fit directly on catch_col.
  # lm_fit was fit to log(catch); logLik(lm_fit) is on the log(catch)
  # scale. Converting back to the catch scale via change-of-variables:
  #   logLik(catch-scale) = logLik(log(catch)-scale) - sum(log(catch))
  # This is the standard correction for comparing AIC/BIC of a
  # log-transformed model against models fit on the raw response.
  # ===========================================================
  .lognormal_ic <- function(lm_fit, y) {
    ll_log_scale <- as.numeric(stats::logLik(lm_fit))
    k            <- attr(stats::logLik(lm_fit), "df")   # params incl. sigma
    n            <- stats::nobs(lm_fit)

    ll_orig_scale <- ll_log_scale - sum(log(y))

    list(
      aic = -2 * ll_orig_scale + 2 * k,
      bic = -2 * ll_orig_scale + log(n) * k
    )
  }

  # ===========================================================
  # Helper: flag families whose Standardized CPUE is nearly
  # identical, so it's obvious this isn't a coding bug -- it
  # usually means the extra flexibility one family offers over
  # another (e.g. NB's dispersion parameter vs Poisson's fixed
  # variance=mean, or Tweedie's var.power converging near 2,
  # i.e. Gamma-like) isn't doing anything on this dataset.
  # ===========================================================
  .check_similar_fits <- function(cpue_table, family_type, all_results) {
    if (length(family_type) < 2L) return(invisible(NULL))

    pairs <- utils::combn(family_type, 2L, simplify = FALSE)
    for (p in pairs) {
      f1 <- p[1]; f2 <- p[2]
      v1 <- cpue_table[[f1]]; v2 <- cpue_table[[f2]]

      # Columns are already rounded to 4 decimals, so exact equality here
      # means "identical up to 4 decimal places".
      if (isTRUE(all.equal(v1, v2, tolerance = 0))) {
        cat(sprintf(
          "  NOTE: %s and %s produced identical Standardized CPUE (matched to 4 decimal places).\n",
          toupper(f1), toupper(f2)))
        cat("        Likely cause: these families share the same linear predictor and log link here, and differ only in the variance function used to weight observations during fitting. When the data is reasonably balanced (or that weighting has little effect), different families can converge to the same fitted coefficients even though their distributional assumptions -- and AIC/BIC -- differ.\n\n")
      }
    }
  }

  # ===========================================================
  # SECTION E: Fit ONE family, return result + plot_data
  # ===========================================================
  .fit_one <- function(fam) {

    optim_method   <- "Nelder-Mead"

    canonical_link <- list(
      gamma     = "log",
      tweedie   = "log",
      gaussian  = "log",
      lognormal = "log",
      poisson   = "log",
      nbinom    = "log"
    )
    zero_tolerant <- c("tweedie", "poisson", "nbinom", "gaussian")

    df <- data
    model_effects <- union(fixed_effects, index_variable)
    model_effects <- c(index_variable, sort(setdiff(model_effects, index_variable)))
    for (v in model_effects) df[[v]] <- as.factor(df[[v]])

    lnk <- if (is.null(link_function)) canonical_link[[fam]] else link_function

    valid_links <- c("log", "identity", "inverse", "logit", "probit", "cloglog")
    if (!lnk %in% valid_links)
      stop(sprintf("link_function must be one of: %s",
                   paste(valid_links, collapse = ", ")))
    if (fam == "nbinom" && lnk != "log")
      stop("'nbinom' supports only link_function = 'log'.")
    if (fam %in% c("gamma", "tweedie", "poisson", "lognormal") &&
        lnk %in% c("logit", "probit", "cloglog"))
      stop(sprintf("link '%s' is invalid for family '%s'.", lnk, fam))
    if (fam == "gaussian" && lnk %in% c("logit", "probit", "cloglog"))
      stop(sprintf("link '%s' is invalid for 'gaussian'.", lnk))

    .log_fam(fam, sprintf("link = '%s'", lnk))

    # ---- FIX: symmetric zero-handling messages across ALL families.
    # Previously, zero_tolerant families (tweedie/poisson/nbinom/gaussian)
    # printed a message() unconditionally (even "0 zero catch values"),
    # while gamma/lognormal only printed when n_zeros > 0 -- meaning gamma
    # and lognormal never printed anything when there were no zeros. Now
    # every family reports its zero-handling status exactly once per run,
    # regardless of whether zeros were present.
    n_zeros <- sum(df[[catch_col]] == 0, na.rm = TRUE)
    if (fam %in% zero_tolerant) {
      if (n_zeros > 0) {
        .log_fam(fam, sprintf("%d zero(s) in '%s' handled natively", n_zeros, catch_col))
        message(sprintf("[GLMstd] %s: %d zero catch values handled natively (no replacement needed)",
                        toupper(fam), n_zeros))
      } else {
        #message(sprintf("[GLMstd] %s: No zero catch values present", toupper(fam)))
      }
    } else {
      if (n_zeros > 0) {
        # ---- FIX: data-driven substitute instead of a fixed 1e-6 ----
        # A fixed constant can be arbitrarily far from the scale of a given
        # dataset (too large relative to tiny catches, or too small/unstable
        # relative to large ones). Using half the smallest observed non-zero
        # value keeps the substitute anchored to the actual data scale, and
        # strictly below every real observation so it doesn't distort the
        # ordering of catches.
        min_nonzero <- suppressWarnings(min(df[[catch_col]][df[[catch_col]] > 0], na.rm = TRUE))
        if (!is.finite(min_nonzero)) {
          stop(sprintf(
            "'%s' has no non-zero, finite values in '%s'; cannot substitute zeros.",
            fam, catch_col))
        }
        zero_sub <- min_nonzero
        df[[catch_col]][df[[catch_col]] == 0] <- zero_sub
        .log_fam(fam, sprintf(
          "%d zero(s) in '%s' replaced with %.6g (half the minimum observed non-zero value)",
          n_zeros, catch_col, zero_sub))
        message(sprintf("[GLMstd] %s: Replaced %d zero catch values with minimum positive catch (%.6g)",
                        toupper(fam), n_zeros, zero_sub))
      } else {
        #message(sprintf("[GLMstd] %s: No zero catch values present", toupper(fam)))
      }
    }

    n_bad <- sum(df[[effort_col]] <= 0, na.rm = TRUE)
    if (n_bad > 0) {
      df <- df[df[[effort_col]] > 0, ]
      .log_fam(fam, sprintf("%d row(s) with effort <= 0 removed", n_bad))

      # --- ADD THIS NOTE ---
      if (n_bad > 0) {
        message(sprintf("[GLMstd] %s: Removed %d row(s) with zero/negative effort",
                        toupper(fam), n_bad))
      }
    }
    df <- droplevels(df)
    if (nrow(df) == 0L) stop("No valid rows remain after data preparation.")
    # ---- NEW: integer coercion for count-based families ----
    if (fam %in% c("poisson", "nbinom")) {
      non_int <- sum(df[[catch_col]] != floor(df[[catch_col]]), na.rm = TRUE)
      if (non_int > 0) {
        df[[catch_col]] <- round(df[[catch_col]])
        .log_fam(fam, "catch data was non-integer; rounded to nearest integer for fitting")
      }
    }
    ff   <- paste(model_effects, collapse = " + ")
    ftxt <- if (fam == "lognormal") {
      paste("log(", catch_col, ") ~", ff, "+ offset(log(", effort_col, "))")
    } else {
      paste(catch_col, "~", ff, "+ offset(log(", effort_col, "))")
    }

    glm_model <- NULL
    var_power_used <- NA_real_

    if (fam == "lognormal") {
      lm_formula <- stats::as.formula(
        paste("log(", catch_col, ") ~", ff, "+ offset(log(", effort_col, "))")
      )
      lm_fit <- stats::lm(lm_formula, data = df)

      glm_model <- list(
        type      = "lognormal_lm",
        beta_hat  = stats::coef(lm_fit),
        X_formula = paste("~", ff),
        sigma     = stats::sigma(lm_fit),   # residual standard error, exact MLE-consistent estimate
        lm_fit    = lm_fit
      )

    } else if (fam == "gamma") {
      suppressWarnings(
        glm_model <- stats::glm(
          stats::as.formula(ftxt), data = df,
          family  = stats::Gamma(link = lnk),
          control = stats::glm.control(maxit = maxit),
          x = FALSE, y = FALSE, model = FALSE)
      )

    } else if (fam == "tweedie") {

      # ---- FIX: estimate var.power from the data instead of a fixed 1.5 ----
      # A hardcoded var.power assumes every dataset this function is ever run
      # on shares the same mean-variance relationship. That's not robust:
      # different species/gear/effort combinations can have genuinely
      # different variance structure. tweedie::tweedie.profile() finds the
      # power (between Poisson-like 1 and Gamma-like 2) best supported by
      # THIS data's likelihood, so the function adapts automatically.
      var_power_used <- tryCatch({
        if (!requireNamespace("tweedie", quietly = TRUE)) {
          .log_fam(fam, "'tweedie' package not installed; used fallback var.power = 1.5 (install.packages('tweedie') for data-driven estimation)")
          1.5
        } else {
          # capture.output swallows tweedie.profile()'s own cat()/print()
          # progress chatter (e.g. "Smooth perhaps inaccurate...") so it
          # doesn't leak into the console mid-loop.
          prof <- NULL
          invisible(utils::capture.output(
            prof <- suppressWarnings(suppressMessages(
              tweedie::tweedie.profile(
                stats::as.formula(ftxt),
                data      = df,
                p.vec     = seq(1.1, 1.9, by = 0.1),
                link.power= 0,
                do.plot   = FALSE,
                method    = "series",
                verbose   = FALSE
              )
            ))
          ))
          est <- prof$p.max
          if (is.null(est) || !is.finite(est)) {
            .log_fam(fam, "tweedie.profile() did not return a usable estimate; used fallback var.power = 1.5")
            1.5
          } else {
            .log_fam(fam, sprintf("var.power estimated from data = %.3f (via tweedie.profile)", est))
            est
          }
        }
      }, error = function(e) {
        .log_fam(fam, sprintf("tweedie.profile() failed (%s); used fallback var.power = 1.5", conditionMessage(e)))
        1.5
      })

      suppressWarnings(
        glm_model <- stats::glm(
          stats::as.formula(ftxt), data = df,
          family  = statmod::tweedie(link.power = 0, var.power = var_power_used),
          control = stats::glm.control(maxit = maxit),
          x = FALSE, y = TRUE, model = FALSE)
      )

    } else if (fam == "gaussian") {
      mustart_vals <- pmax(df[[catch_col]], 1e-3)  # initialization only; y itself stays untouched
      glm_model <- stats::glm(
        stats::as.formula(ftxt), data = df,
        family   = stats::gaussian(link = lnk),
        mustart  = mustart_vals,
        control  = stats::glm.control(maxit = maxit),
        x = FALSE, y = FALSE, model = FALSE)

    } else if (fam == "poisson") {
      #if (any(df[[catch_col]] != floor(df[[catch_col]]), na.rm = TRUE))
      #message("[GLMstd] NOTE: Non-integer catch detected. Poisson fitted as quasi-count approximation.")
      suppressWarnings(
        glm_model <- stats::glm(
          stats::as.formula(ftxt), data = df,
          family  = stats::poisson(link = lnk),
          control = stats::glm.control(maxit = maxit),
          x = FALSE, y = FALSE, model = FALSE)
      )

    } else if (fam == "nbinom") {
      suppressWarnings(
        glm_model <- MASS::glm.nb(
          stats::as.formula(ftxt), data = df,
          control = stats::glm.control(maxit = maxit),
          x = FALSE, y = FALSE, model = FALSE)
      )
      #message(sprintf("[GLMstd] NegBin theta = %.4f.", glm_model$theta))
    }

    # ---- Model selection criteria ----
    if (fam == "lognormal") {
      # model_aic <- stats::AIC(glm_model$lm_fit)
      # model_bic <- stats::BIC(glm_model$lm_fit)  #Replaced this because AIC/BIC was not comparable
      ic        <- .lognormal_ic(glm_model$lm_fit, df[[catch_col]])
      model_aic <- ic$aic
      model_bic <- ic$bic

    } else if (fam == "tweedie") {
      if (requireNamespace("tweedie", quietly = TRUE)) {
        model_aic <- tweedie::tweedie_AIC(glm_model, k = 2, verbose = FALSE)
        model_bic <- tweedie::tweedie_AIC(glm_model, k = log(stats::nobs(glm_model)), verbose = FALSE)
      } else {
        model_aic <- NA_real_
        model_bic <- NA_real_
        # message("[GLMstd] NOTE: Install 'tweedie' package (install.packages('tweedie')) ",
        #         "to compute AIC/BIC for the tweedie family.")
      }

      # } else if (fam == "poisson" && any(df[[catch_col]] != floor(df[[catch_col]]), na.rm = TRUE)) {
      #   model_aic <- NA_real_
      #   model_bic <- NA_real_
      #   message("[GLMstd] NOTE: AIC/BIC for 'poisson' set to NA (non-integer catch makes likelihood undefined).")

    } else {
      model_aic <- stats::AIC(glm_model)
      model_bic <- stats::BIC(glm_model)
    }
    pred_lv  <- levels(df[[index_variable]])

    # Reference levels only for OTHER modeling effects -- index_variable is
    # the one that varies across the prediction loop below, so it's excluded
    # here rather than held fixed.
    other_effects <- setdiff(model_effects, index_variable)
    ref_vals <- stats::setNames(
      lapply(other_effects, function(v) levels(df[[v]])[1]),
      other_effects
    )
    results <- data.frame()

    for (lv in pred_lv) {
      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(df[[v]]))
      nd[[effort_col]] <- 1

      cpue <- if (fam == "lognormal") {
        Xn <- stats::model.matrix(
          stats::as.formula(glm_model$X_formula), data = nd)
        exp(as.numeric(Xn %*% glm_model$beta_hat) + glm_model$sigma^2 / 2)
      } else {
        suppressWarnings(
          as.numeric(predict(glm_model, newdata = nd, type = "response")[1])
        )
      }

      results <- rbind(results,
                       data.frame(Index_Level       = as.character(lv),
                                  Standardized_CPUE = round(cpue, 4L)))
    }
    row.names(results) <- NULL

    out <- list(
      model              = glm_model,
      family_used        = fam,
      link_used          = lnk,
      formula_used       = ftxt,
      fixed_effects      = fixed_effects,
      index_variable     = index_variable,
      standardized_index = results
    )
    out$aic <- model_aic
    out$bic <- model_bic

    if (fam == "tweedie")   out$estimated_var_power <- var_power_used
    if (fam == "lognormal") { out$coefficients <- glm_model$beta_hat
    out$sigma        <- glm_model$sigma }
    if (fam == "nbinom")    out$theta          <- glm_model$theta

    ci <- stats::aggregate(df[[catch_col]] ~ df[[index_variable]], FUN = sum)
    names(ci) <- c("Index_Level", "Total_Catch")
    ci$Index_Level <- as.character(ci$Index_Level)

    pd <- merge(
      transform(results, Index_Level = as.character(Index_Level)),
      ci, by = "Index_Level", all.x = TRUE
    )
    pd$Total_Catch[is.na(pd$Total_Catch)] <- 0
    pd$Year <- as.numeric(pd$Index_Level)

    plot_data <- list(
      year        = pd$Year,
      Total_Catch = pd$Total_Catch,
      CPUE        = pd$Standardized_CPUE,
      title       = sprintf("%s Dist.", tools::toTitleCase(fam))
    )

    list(result = out, plot_data = plot_data)
  }


  # ===========================================================
  # SECTION F: Main body — validate, loop, compare, plot
  # ===========================================================

  # F-1  Validate families
  valid_fam   <- c("gamma", "tweedie", "gaussian",
                   "lognormal", "poisson", "nbinom")
  family_type <- unique(tolower(trimws(as.character(family_type))))
  bad_fam     <- setdiff(family_type, valid_fam)
  if (length(bad_fam) > 0L)
    stop(sprintf(
      "Unknown family: %s. Choose from: %s",
      paste(bad_fam,  collapse = ", "),
      paste(valid_fam, collapse = ", ")
    ))

  # F-2  Package checks
  if ("nbinom" %in% family_type && !requireNamespace("MASS", quietly = TRUE))
    stop("Install 'MASS': install.packages('MASS')")

  # F-3  Column validation
  req  <- c(catch_col, effort_col, fixed_effects, year_col)
  miss <- setdiff(req, names(data))
  if (length(miss) > 0L)
    stop(paste("Missing columns:", paste(miss, collapse = ", ")))

  # F-4  Fit all families
  # ---- FIX: one family's convergence failure no longer aborts the whole run.
  # Families like gamma/nbinom can fail on datasets with quasi-separation
  # (e.g. a gear with mostly-zero catch mixed with sparse extreme years)
  # while other families fit fine. We catch failures per-family, warn with
  # the reason, and continue with whatever converged. If NONE converge, we
  # stop with a clear message rather than returning an empty result.
  orig_family_type <- family_type   # preserve original order/set for later reference, even after filtering

  all_results  <- vector("list", length(family_type))
  all_plotdata <- vector("list", length(family_type))
  names(all_results) <- names(all_plotdata) <- family_type
  failed_fam <- character(0)

  for (fam in family_type) {
    fit <- tryCatch(
      .fit_one(fam),
      error = function(e) {
        .log_fam(fam, sprintf("FAILED - %s", conditionMessage(e)))
        NULL
      }
    )
    if (is.null(fit)) {
      failed_fam <- c(failed_fam, fam)
    } else {
      all_results[[fam]]  <- fit$result
      all_plotdata[[fam]] <- fit$plot_data
    }
  }

  if (length(failed_fam) > 0L) {
    all_results[failed_fam]  <- NULL
    all_plotdata[failed_fam] <- NULL
    family_type <- setdiff(family_type, failed_fam)
  }

  # ---- FIX: nothing is printed here anymore. The per-family fitting notes
  # (zero handling, link used, failure reasons, etc.) are kept in
  # .fam_log$store for programmatic access (attached as an attribute on the
  # final return value below) but are no longer echoed to the console, so
  # the first thing the user sees is the Model Selection Criteria table.
  fitting_notes <- .fam_log$store

  if (length(family_type) == 0L) {
    cat("  NOTE: None of the selected distribution(s) converged for this dataset.\n\n")

    return(invisible(NULL))
  }

  # F-5  Compute nominal CPUE (same for all families)
  nom_df <- .compute_nominal_cpue()

  # F-6  Single family path
  # F-6  Single family path
  if (length(family_type) == 1L) {

    single_result <- all_results[[1L]]
    si            <- single_result$standardized_index

    ci <- stats::aggregate(data[[catch_col]] ~ data[[index_variable]], FUN = sum)
    names(ci) <- c("Index_Level", "Total_Catch")
    ci$Index_Level <- as.character(ci$Index_Level)

    single_table <- merge(
      transform(si, Index_Level = as.character(Index_Level)),
      ci, by = "Index_Level", all.x = TRUE
    )
    single_table$Total_Catch[is.na(single_table$Total_Catch)] <- 0

    # ---- merge in Nominal CPUE ----
    nom_merge <- nom_df[, c("Year", "Nominal_CPUE")]
    names(nom_merge) <- c("Index_Level", "Nominal_CPUE")
    nom_merge$Index_Level <- as.character(nom_merge$Index_Level)
    single_table <- merge(single_table, nom_merge, by = "Index_Level", all.x = TRUE)

    names(single_table)[names(single_table) == "Index_Level"] <- index_variable

    numeric_idx <- suppressWarnings(as.numeric(single_table[[index_variable]]))
    if (all(!is.na(numeric_idx))) {
      single_table[[index_variable]] <- numeric_idx
      single_table <- single_table[order(single_table[[index_variable]]), ]
    }

    single_table <- single_table[, c(index_variable, "Total_Catch",
                                     "Nominal_CPUE", "Standardized_CPUE")]
    single_table$Nominal_CPUE      <- round(single_table$Nominal_CPUE, 4L)
    single_table$Standardized_CPUE <- round(single_table$Standardized_CPUE, 4L)
    row.names(single_table) <- NULL

    aic_bic_table <- data.frame(
      Distribution = toupper(family_type),
      AIC = round(single_result$aic, 4L),
      BIC = round(single_result$bic, 4L)
    )
    cat("  Model Selection Criteria\n")
    cat(strrep("-", 35), "\n\n", sep = "")
    print(aic_bic_table, row.names = FALSE)
    cat("\n")
    ##To show the NOTE final
    if (length(failed_fam) > 0L)
      cat(sprintf("  NOTE: Skipped due to non-convergence: %s.\n\n",
                  paste(failed_fam, collapse = ", ")))
    cat("\n", strrep("=", 35), "\n", sep = "")
    cat("  Standardized CPUE Table [", toupper(family_type), "]\n", sep = "")
    cat(strrep("=", 35), "\n\n", sep = "")
    #print(single_table, row.names = FALSE, digits = 4L)
    #cat("\n", strrep("=", 65), "\n\n", sep = "")

    # Combine Nominal + Standardized panel into a single tiled frame
    nominal_panel_data <- list(
      year        = nom_df$Year,
      Total_Catch = nom_df$Total_Catch,
      CPUE        = nom_df$Nominal_CPUE,
      title       = "Nominal CPUE"
    )
    .tile_panels_grouped(nominal_panel_data, all_plotdata, "CPUE vs Total Catch")
    attr(single_table, "fitting_notes") <- fitting_notes
    attr(single_table, "skipped_families") <- failed_fam
    return(single_table)
  }
  # F-7  Multi-family: comparison table
  cpue_table <- data.frame(
    Index = as.numeric(as.character(
      all_results[[1L]]$standardized_index$Index_Level
    ))
  )
  names(cpue_table)[1L] <- index_variable

  for (fam in family_type)
    cpue_table[[fam]] <- round(
      all_results[[fam]]$standardized_index$Standardized_CPUE, 4L)

  ci <- stats::aggregate(data[[catch_col]] ~ data[[index_variable]], FUN = sum)
  names(ci) <- c(index_variable, "Total_Catch")
  ci[[index_variable]] <- as.numeric(as.character(ci[[index_variable]]))

  cpue_table <- merge(cpue_table, ci, by = index_variable, all.x = TRUE)
  cpue_table$Total_Catch[is.na(cpue_table$Total_Catch)] <- 0

  # ---- merge in Nominal CPUE ----
  nom_merge <- nom_df[, c("Year", "Nominal_CPUE")]
  names(nom_merge)[1L] <- index_variable
  nom_merge$Nominal_CPUE <- round(nom_merge$Nominal_CPUE, 4L)
  cpue_table <- merge(cpue_table, nom_merge, by = index_variable, all.x = TRUE)

  cpue_table <- cpue_table[order(cpue_table[[index_variable]]), ]
  cpue_table <- cpue_table[, c(index_variable, "Total_Catch",
                               "Nominal_CPUE", family_type)]
  row.names(cpue_table) <- NULL
  ##################
  aic_bic_table <- data.frame(
    Distribution = toupper(family_type),
    AIC = round(vapply(family_type, function(f) all_results[[f]]$aic, numeric(1L)), 4L),
    BIC = round(vapply(family_type, function(f) all_results[[f]]$bic, numeric(1L)), 4L)
  )
  cat("  Model Selection Criteria (AIC / BIC)\n")
  cat(strrep("-", 35), "\n\n", sep = "")
  print(aic_bic_table, row.names = FALSE)
  cat("\n")
  if (length(failed_fam) > 0L)
    cat(sprintf("  NOTE: Skipped due to non-convergence: %s.\n\n",
                paste(failed_fam, collapse = ", ")))

  cat("\n", strrep("=", 35), "\n", sep = "")
  cat("  Standardized CPUE Comparison Table\n")
  cat(strrep("=", 35), "\n\n", sep = "")

  # ---- Diagnostic: flag any pair of families with near-identical
  # Standardized CPUE, so it's clear this reflects the data (e.g. no
  # meaningful overdispersion) rather than a bug in the fitting code.
  .check_similar_fits(cpue_table, family_type, all_results)

  # F-8/F-9  Combine Nominal + all family panels into one tiled frame
  nominal_panel_data <- list(
    year        = nom_df$Year,
    Total_Catch = nom_df$Total_Catch,
    CPUE        = nom_df$Nominal_CPUE,
    title       = "Nominal CPUE"
  )
  .tile_panels_grouped(nominal_panel_data, all_plotdata, "CPUE vs Total Catch")
  # F-10 Return
  attr(cpue_table, "fitting_notes") <- fitting_notes
  attr(cpue_table, "skipped_families") <- failed_fam
  return(cpue_table)
}


# ===========================================================
# Print method for GLMstd_multi
# ===========================================================
#' @export
print.GLMstd_multi <- function(x, ...) {
  families <- setdiff(names(x), "comparison_table")
  cat("\nGLMstd multi-family result.\n")
  cat("Families fitted:", paste(families, collapse = ", "), "\n")
  if (!is.null(x$comparison_table)) {
    cat("\nStandardized CPUE Comparison Table:\n")
    print(x$comparison_table, row.names = FALSE, digits = 4L)
  }
  #invisible(x)
  return(x)
}