R/percent_overlap.R

In phontrast: Contrast and Separation Metrics for Phonological Categories

#' Proportional overlap between two distributions via KDE
#'
#' Computes the proportional overlap (shared area) between two categories
#' in an n-dimensional acoustic space using multivariate kernel density
#' estimation. Despite the historical function name, the return value is a
#' 0--1 proportion: 0 = no overlap, 1 = identical.
#'
#' @param data Data frame.
#' @param features Character vector of numeric feature columns.
#' @param category_col String; exactly two categories.
#' @param bw Bandwidth selection method. Uses the same options as
#'   \code{jsd_kde_nd()}: \code{"Hpi"}, \code{"Hscv"}, \code{"Hpi.diag"},
#'   or \code{"scott.diag"}.
#' @param eval_on Where to evaluate the KDEs. Uses the same options as
#'   \code{jsd_kde_nd()}: \code{"pooled"}, \code{"group1"},
#'   \code{"group2"}, or \code{"pooled_sample"}.
#' @param eval_n Optional positive integer giving the maximum number of
#'   evaluation points to use.
#' @param eval_seed Optional integer seed used only when \code{eval_n} causes
#'   evaluation-point subsampling.
#' @param engine KDE evaluation engine. Uses the same options as
#'   \code{jsd_kde_nd()}: \code{"ks"}, \code{"fast_diag"}, or
#'   \code{"fast_diagonal"}.
#' @param chunk_size Positive integer controlling the number of evaluation
#'   points processed per chunk by \code{engine = "fast_diag"}.
#' @param method Estimator: \code{"mc"} (default) for the Monte-Carlo plug-in
#'   estimate of the overlapping coefficient, or \code{"legacy"} for the
#'   pre-1.2.0 self-normalized sample-point estimate. \code{eval_on} applies to
#'   \code{"legacy"} only. Ignored when \code{density = "mvnorm"}.
#' @param density Density model behind the estimate: \code{"kde"} (default)
#'   estimates each category's density by kernel density estimation;
#'   \code{"mvnorm"} fits one multivariate normal per category and estimates the
#'   overlapping coefficient between the two Gaussians by Monte-Carlo. Under
#'   \code{"mvnorm"} the KDE-specific arguments (\code{bw}, \code{engine},
#'   \code{eval_on}, \code{chunk_size}, \code{method}, \code{eval_n}) do not
#'   apply; the Monte-Carlo sample size is set by \code{mc_n} and
#'   \code{eval_seed} makes the draw reproducible.
#' @param mc_n Positive integer; number of Monte-Carlo samples drawn from each
#'   fitted Gaussian when \code{density = "mvnorm"} (default \code{10000}). The
#'   estimator draws \code{mc_n} fresh points from each category's fitted
#'   Gaussian to estimate the overlapping coefficient between the two Gaussians.
#'   Larger values reduce Monte-Carlo variance. Ignored when
#'   \code{density = "kde"}.
#' @param bw_scale Positive number multiplying the selected kernel bandwidth on
#'   the standard-deviation scale (univariate bandwidths by \code{bw_scale},
#'   bandwidth matrices by \code{bw_scale^2}); see \code{jsd_kde_nd()}.
#'   Ignored when \code{density = "mvnorm"}.
#' @param ... Reserved for future extensions; currently unused.
#'
#' @return Numeric scalar proportion in \code{[0, 1]}.
#' @export
percent_overlap_kde <- function(data,
                                features,
                                category_col,
                                bw = c("Hpi", "Hscv", "Hpi.diag", "scott.diag"),
                                eval_on = c("pooled", "group1", "group2", "pooled_sample"),
                                eval_n = NULL,
                                eval_seed = NULL,
                                engine = c("ks", "fast_diag", "fast_diagonal"),
                                chunk_size = 1000L,
                                method = c("mc", "legacy"),
                                density = c("kde", "mvnorm"),
                                mc_n = 10000L,
                                bw_scale = 1,
                                ...) {

  .validate_metric_inputs(data, features, category_col)
  method <- match.arg(method)
  density <- match.arg(density)
  .check_bw_scale(bw_scale)

  if (identical(density, "mvnorm")) {
    mc <- .mvnorm_mc_pair(
      data = data,
      features = features,
      category_col = category_col,
      mc_n = mc_n,
      eval_seed = eval_seed,
      metric = "percent_overlap_kde()"
    )
    return(.overlap_mc(mc))
  }

  if (identical(method, "mc")) {
    mc <- .kde_mc_pair(
      data = data,
      features = features,
      category_col = category_col,
      bw = bw,
      eval_n = eval_n,
      eval_seed = eval_seed,
      engine = engine,
      chunk_size = chunk_size,
      metric = "percent_overlap_kde()",
      bw_scale = bw_scale
    )
    return(.overlap_mc(mc))
  }

  dens <- .kde_density_pair(
    data = data,
    features = features,
    category_col = category_col,
    bw = bw,
    eval_on = eval_on,
    eval_n = eval_n,
    eval_seed = eval_seed,
    engine = engine,
    chunk_size = chunk_size,
    metric = "percent_overlap_kde()",
    bw_scale = bw_scale
  )

  # Normalize to discrete probability masses on a shared grid
  p <- dens$p / sum(dens$p)
  q <- dens$q / sum(dens$q)

  # Overlap is the shared area on the grid: sum(min(p, q))
  overlap <- sum(pmin(p, q))

  # Bound numerically
  overlap <- max(min(overlap, 1), 0)

  overlap
}

#' Estimate proportional overlap globally or by group
#'
#' Unified front-end for KDE-based proportional overlap between two categories.
#' The returned \code{overlap} column is a 0--1 proportion, not a 0--100
#' percentage.
#'
#' @inheritParams estimate_jsd
#' @param group_col Optional character vector of one or more grouping columns.
#'   If provided, returns per-group overlap. Multiple grouping columns are
#'   combined into a labeled \code{group} value such as
#'   \code{"Sex=F | Style=read"}.
#' @param bw Bandwidth selection method passed to \code{percent_overlap_kde()}.
#' @param eval_on KDE evaluation points passed to \code{percent_overlap_kde()}.
#' @param eval_n Optional maximum number of KDE evaluation points.
#' @param eval_seed Optional integer seed for KDE evaluation-point subsampling.
#' @param engine KDE evaluation engine passed to \code{percent_overlap_kde()}.
#'   \code{"fast_diagonal"} is accepted as an alias for \code{"fast_diag"}.
#' @param chunk_size Chunk size for \code{engine = "fast_diag"}.
#' @param method Estimator passed to \code{percent_overlap_kde()}: \code{"mc"}
#'   (default) or \code{"legacy"} (pre-1.2.0 self-normalized estimate). Ignored
#'   when \code{density = "mvnorm"}.
#' @param density Density model passed to \code{percent_overlap_kde()}:
#'   \code{"kde"} (default) or \code{"mvnorm"} (fit one multivariate normal per
#'   category and estimate the overlapping coefficient between the two Gaussians
#'   by Monte-Carlo).
#' @param mc_n Positive integer; number of Monte-Carlo samples drawn from each
#'   fitted Gaussian when \code{density = "mvnorm"} (default \code{10000}).
#'   Ignored when \code{density = "kde"}.
#' @param bw_scale Positive bandwidth multiplier passed to
#'   \code{percent_overlap_kde()} (default \code{1}); see \code{jsd_kde_nd()}.
#' @param ... Additional arguments passed to \code{percent_overlap_kde()}.
#'
#' @return A tibble (global = one row; grouped = one per group) with
#'   \code{overlap} as a 0--1 proportion.
#' @export
estimate_overlap <- function(data,
                             features,
                             category_col,
                             group_col  = NULL,
                             min_tokens = 20,
                             bw = c("Hpi", "Hscv", "Hpi.diag", "scott.diag"),
                             eval_on = c("pooled", "group1", "group2", "pooled_sample"),
                             eval_n = NULL,
                             eval_seed = NULL,
                             engine = c("ks", "fast_diag", "fast_diagonal"),
                             chunk_size = 1000L,
                             method = c("mc", "legacy"),
                             density = c("kde", "mvnorm"),
                             mc_n = 10000L,
                             bw_scale = 1,
                             ...) {

  bw <- match.arg(bw)
  eval_on <- match.arg(eval_on)
  engine <- .match_kde_engine(engine)
  method <- match.arg(method)
  density <- match.arg(density)
  .check_bw_scale(bw_scale)
  .check_positive_count(min_tokens, "min_tokens")
  .validate_metric_inputs(data, features, category_col, group_col)

  if (is.null(group_col)) {
    # ---- Global ----
    keep_cols <- c(category_col, features)
    df <- .metric_data(data, keep_cols)

    n <- nrow(df)
    if (n < min_tokens)
      stop("Not enough tokens for global percent overlap.")

    ov <- percent_overlap_kde(
      data         = df,
      features     = features,
      category_col = category_col,
      bw           = bw,
      eval_on      = eval_on,
      eval_n       = eval_n,
      eval_seed    = eval_seed,
      engine       = engine,
      chunk_size   = chunk_size,
      method       = method,
      density      = density,
      mc_n         = mc_n,
      bw_scale     = bw_scale,
      ...
    )

    return(tibble::tibble(
      scope        = "global",
      n_tokens     = n,
      overlap      = ov
    ))
  }

  # ---- Grouped ----
  group_col <- .check_group_cols(group_col)
  .check_columns(data, c(group_col, category_col, features))
  data <- .metric_data(data, c(group_col, category_col, features))

  groups <- .split_groups(data, group_col)

  out <- lapply(groups, function(df_g) {
    n_tok <- nrow(df_g)
    if (n_tok < min_tokens ||
        .observed_n_categories(df_g[[category_col]]) != 2L)
      return(NULL)

    ov <- tryCatch(
      percent_overlap_kde(
        data         = df_g,
        features     = features,
        category_col = category_col,
        bw           = bw,
        eval_on      = eval_on,
        eval_n       = eval_n,
        eval_seed    = eval_seed,
        engine       = engine,
        chunk_size   = chunk_size,
        method       = method,
        density      = density,
        mc_n         = mc_n,
        bw_scale     = bw_scale,
        ...
      ),
      error = function(e) NA_real_
    )

    data.frame(
      scope    = "group",
      group    = .group_label(df_g, group_col),
      n_tokens = n_tok,
      overlap  = ov,
      stringsAsFactors = FALSE
    )
  })

  out <- do.call(rbind, out)
  if (is.null(out)) {
    out <- data.frame(
      scope = character(),
      group = character(),
      n_tokens = integer(),
      overlap = numeric(),
      stringsAsFactors = FALSE
    )
  }
  rownames(out) <- NULL
  tibble::as_tibble(.warn_failed_groups(out, "overlap", "estimate_overlap()"))
}