R/tracer_pull_analysis.R

In OceanTracer: Marine Source Unmixing Toolkit

#' @title Analyze the Marginal Effect (Pull) of Individual Tracers
#'
#' @description Evaluates the directional bias introduced by individual tracers within 
#' Consistent Tracer Selection (CTS) subset models. This function unpacks unconstrained 
#' source apportionments to reveal if specific tracers are systematically pulling 
#' source contributions into physically impossible spaces (e.g., >1 or <0).
#'
#' @details 
#' Global unmixing models (using all tracers simultaneously) can suffer from "tracer compensation," 
#' where extreme mathematical contradictions between non-conservative tracers cancel each other out, 
#' resulting in a misleadingly plausible central solution. 
#' 
#' This function exposes that illusion. By taking the output of a CTS run (e.g., \code{\link{CTS_n_sources}}), 
#' it isolates every debate a specific tracer participated in. It then calculates the marginal effect 
#' (median, mean, and SD) that the tracer has on the fractional weight of each source.
#' 
#' Visualizations are automatically generated and stored in the returned list:
#' \itemize{
#'   \item \strong{boxplot:} Displays the distribution of calculated weights for each source, grouped by tracer. 
#'   Red dashed lines denote the physical boundaries of 0 and 1.
#'   \item \strong{heatmap:} Displays the median weight pull. Deep red indicates a strong positive pull 
#'   (overestimation), while deep blue indicates a negative pull (pushing the source into negative mass space).
#' }
#' 
#' @param cts_result A data frame containing the output from a CTS run. It must contain an 
#'   \code{id} column with space-separated tracer names, and weight columns named \code{w1}, \code{w2}, etc.
#'
#' @return A named list containing four elements:
#' \describe{
#'   \item{\code{summary}}{A \code{\link[dplyr]{tibble}} containing the Mean, Median, and SD of the weights forced by each tracer on each source.}
#'   \item{\code{boxplot}}{A \code{\link[ggplot2]{ggplot}} object showing the distribution of source weights grouped by tracer.}
#'   \item{\code{heatmap}}{A \code{\link[ggplot2]{ggplot}} object visualizing the median tracer bias across all sources.}
#'   \item{\code{raw_expanded}}{A \code{\link[dplyr]{tibble}} of the raw, pivoted data used for calculations and plotting.}
#' }
#'
#' @importFrom dplyr select group_by summarise rename all_of
#' @importFrom tidyr separate_rows pivot_longer
#' @importFrom ggplot2 ggplot aes geom_boxplot geom_hline facet_wrap theme_bw theme element_text element_rect labs geom_tile scale_fill_gradient2 theme_minimal
#'
#' @examples
#' # Create a toy dataset mimicking a CTS model output
#' cts_out <- data.frame(
#'   id = c("TracerA TracerB", "TracerA TracerC", "TracerB TracerC"),
#'   w1 = c(0.3, 0.4, 0.2),
#'   w2 = c(0.7, 0.6, 0.8)
#' )
#' 
#' # Analyze tracer pull
#' pull_results <- tracer_pull_analysis(cts_out)
#' 
#' # View summary statistics and display boxplot
#' print(pull_results$summary)
#' print(pull_results$boxplot)
#' 
#' @export
tracer_pull_analysis <- function(cts_result) {
  
  # 1. Identify weight columns dynamically (w1, w2, w3, etc.)
  w_cols <- grep("^w[0-9]+$", colnames(cts_result), value = TRUE)
  if (length(w_cols) == 0L) {
    stop("No weight columns matching pattern '^w[0-9]+$' (e.g., 'w1', 'w2') were found in 'cts_result'.")
  }
  
  # 2. Reshape data: Unpack tracer combinations and pivot weight columns
  df_long <- cts_result |>
    dplyr::select(dplyr::all_of(c("id", w_cols))) |>
    tidyr::separate_rows(.data$id, sep = " ") |>
    dplyr::rename(Tracer = .data$id) |>
    tidyr::pivot_longer(
      cols      = dplyr::all_of(w_cols),
      names_to  = "Source",
      values_to = "Weight"
    )
  
  # 3. Calculate central tendency and variance of the tracer "pull"
  summary_pull <- df_long |>
    dplyr::group_by(.data$Tracer, .data$Source) |>
    dplyr::summarise(
      Mean_Weight   = mean(.data$Weight, na.rm = TRUE),
      Median_Weight = stats::median(.data$Weight, na.rm = TRUE),
      SD_Weight     = stats::sd(.data$Weight, na.rm = TRUE),
      .groups       = "drop"
    )
  
  # 4. Plot 1: Boxplots of Tracer Pull (Highlights mathematical extreme compensations)
  p_box <- ggplot2::ggplot(df_long, ggplot2::aes(x = .data$Tracer, y = .data$Weight, fill = .data$Tracer)) +
    ggplot2::geom_boxplot(alpha = 0.7, outlier.size = 0.5, outlier.alpha = 0.3) +
    ggplot2::geom_hline(yintercept = 0, linetype = "dashed", color = "red", linewidth = 0.8) +
    ggplot2::geom_hline(yintercept = 1, linetype = "dashed", color = "red", linewidth = 0.8) +
    ggplot2::facet_wrap(~ .data$Source, scales = "free_y") +
    ggplot2::theme_bw() +
    ggplot2::theme(
      axis.text.x      = ggplot2::element_text(angle = 45, hjust = 1),
      legend.position  = "none",
      strip.background = ggplot2::element_rect(fill = "grey90")
    ) +
    ggplot2::labs(
      title    = "Tracer Pull: Distribution of Source Contributions",
      subtitle = "Values outside 0-1 (red lines) indicate the tracer forces impossible physical math.",
      x        = "Tracer Included in Model",
      y        = "Calculated Source Weight"
    )
  
  # 5. Plot 2: Heatmap of Median Pull (Highlights directional bias towards specific sources)
  p_heat <- ggplot2::ggplot(summary_pull, ggplot2::aes(x = .data$Source, y = .data$Tracer, fill = .data$Median_Weight)) +
    ggplot2::geom_tile(color = "white", linewidth = 0.5) +
    ggplot2::scale_fill_gradient2(
      low      = "dodgerblue",
      mid      = "white",
      high     = "firebrick",
      midpoint = 0.5,
      name     = "Median\nWeight"
    ) +
    ggplot2::theme_minimal() +
    ggplot2::labs(
      title    = "Heatmap of Tracer Bias",
      subtitle = "Red = Pulls strongly towards source. Blue = Pushes into negative mass space.",
      x        = "Sediment Source",
      y        = "Tracer"
    )
  
  return(list(
    summary      = summary_pull,
    boxplot      = p_box,
    heatmap      = p_heat,
    raw_expanded = df_long
  ))
}