Iss20

NAMESPACE

@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(print,eval_data)
+S3method(print,k_marginals)
 export("%>%")
 export(add_discriminator_auc)
 export(add_pmse)
@@ -26,6 +27,7 @@ export(prep_discrete_eval_data)
 export(util_ci_overlap)
 export(util_co_occurrence)
 export(util_corr_fit)
+export(util_k_marginals)
 export(util_ks_distance)
 export(util_moments)
 export(util_percentiles)

NEWS.md

@@ -1,8 +1,9 @@
-# syntheval (development version)
+# syntheval 0.1.0 (development)
 
-* Add optional parallelization to the bootstrapped null pMSEs in `add_pmse_ratio()` via `furrr`; enable it with `future::plan()`. Results for a given seed change relative to earlier versions because bootstrap iterations now use parallel-safe L'Ecuyer-CMRG streams.
+* Add `util_k_marginals()` to calculate the k-marginals metric for 1-, 2-, and 3-way marginals, with worst-marginal and worst-cell output, marginal sampling with priority variables, sample weights, a `synth_vars` flag to restrict marginals to synthesized variables, `na.rm` handling of missing values, `group_by` stratification with per-stratum scores, and discretization of numeric variables. (#20)
+* Add optional parallelization to the bootstrapped null pMSEs in `add_pmse_ratio()` via `furrr`; enable it with `future::plan()`. Results for a given seed change relative to earlier versions because bootstrap iterations now use parallel-safe L'Ecuyer-CMRG streams. (#113)
 
-# syntheval 0.0.5
+# syntheval 0.0.5 (release)
 
 * Remove `util_tails()`
 * Update deprecated dplyr code. 

R/compute_marginal_cells.R

@@ -0,0 +1,105 @@
+#' @title Compute per-cell proportion differences for the k-marginals metric
+#'
+#' @description For each supplied variable combination, computes the marginal
+#' cell proportions in the synthetic and confidential data and their absolute
+#' differences. Cells absent from one dataset count as proportion zero.
+#'
+#' @param synth_data A tibble with synthetic data.
+#' @param conf_data A tibble with confidential data.
+#' @param combos A character matrix with one row per variable combination.
+#' @param weight_var Optional character name of a numeric sample-weight
+#' column; when set, cell proportions are weight shares instead of row
+#' shares. Defaults to `NULL` (unweighted).
+#' @param na.rm A logical for dropping rows with a missing value from each
+#' marginal that uses the affected variable.
+#' @param allow_empty_synth A logical for permitting the synthetic data to
+#' contribute no rows to a marginal (a stratum the synthesis never produced),
+#' in which case its proportions are zero. A confidential marginal with no
+#' rows always errors, since there is nothing to score against.
+#'
+#' @return A tibble with one row per cell: `variables`, `cell`,
+#' `prop_synth`, `prop_conf`, and `abs_diff`.
+#'
+.compute_marginal_cells <- function(
+  synth_data,
+  conf_data,
+  combos,
+  weight_var = NULL,
+  na.rm = FALSE,
+  allow_empty_synth = FALSE
+) {
+  # cell proportions for one dataset over one set of variables; weighted
+  # proportions are weight shares instead of row shares
+  process_data <- function(data, vars, prop_name, allow_empty) {
+    if (na.rm) {
+      data <- dplyr::filter(
+        data,
+        !dplyr::if_any(.cols = dplyr::all_of(vars), .fns = is.na)
+      )
+
+      if (nrow(data) == 0 && !allow_empty) {
+        stop(
+          "no rows remain for the marginal over ",
+          paste(vars, collapse = ", "),
+          " after removing missing values"
+        )
+      }
+    }
+
+    if (is.null(weight_var)) {
+      counts <- dplyr::count(data, dplyr::across(dplyr::all_of(vars)))
+    } else {
+      counts <- dplyr::count(
+        data,
+        dplyr::across(dplyr::all_of(vars)),
+        wt = .data[[weight_var]]
+      )
+    }
+
+    props <- counts |>
+      dplyr::mutate("{prop_name}" := .data$n / sum(.data$n)) |>
+      dplyr::select(-"n")
+
+    return(props)
+  }
+
+  # per-cell differences for one set of variables; cells absent from one
+  # dataset count as 0
+  marginal_cells <- function(vars) {
+    # only the synthetic side may be empty (a stratum the synthesis never
+    # produced); a confidential marginal with no rows has nothing to score
+    # against and errors inside process_data
+    cells <- dplyr::full_join(
+      process_data(
+        data = synth_data, vars = vars, prop_name = "prop_synth",
+        allow_empty = allow_empty_synth
+      ),
+      process_data(
+        data = conf_data, vars = vars, prop_name = "prop_conf",
+        allow_empty = FALSE
+      ),
+      by = vars
+    ) |>
+      tidyr::replace_na(replace = list(prop_synth = 0, prop_conf = 0)) |>
+      tidyr::unite(col = "cell", dplyr::all_of(vars), sep = ", ") |>
+      dplyr::mutate(
+        variables = paste(vars, collapse = ", "),
+        abs_diff = abs(.data$prop_synth - .data$prop_conf)
+      ) |>
+      dplyr::select(
+        "variables", "cell", "prop_synth", "prop_conf", "abs_diff"
+      )
+    # variables disambiguates cells across combinations and drives the
+    # per-combination summary; the prop columns show the direction of the
+    # discrepancy, not just its size
+    return(cells)
+  }
+
+  cells <- purrr::map(
+    .x = seq_len(nrow(combos)),
+    .f = \(i) marginal_cells(vars = combos[i, ])
+  ) |>
+    purrr::list_rbind()
+
+  return(cells)
+}

R/discretize_k_marginal_vars.R

@@ -0,0 +1,94 @@
+#' @title Discretize numeric variables for the k-marginals metric
+#'
+#' @description Discretizes every numeric variable among `vars` into `bins`
+#' bins with interior cut points derived from the observed (non-missing)
+#' confidential values and applied to both datasets. The outer bins extend to
+#' +/-Inf so synthetic values outside the confidential range land in edge
+#' bins.
+#'
+#' @param synth_data A tibble with synthetic data.
+#' @param conf_data A tibble with confidential data.
+#' @param vars Character vector of candidate variables; only numeric ones are
+#' discretized.
+#' @param bins Single integer >= 2 giving the number of bins. Fewer bins are
+#' produced, with a warning, if tied cut points collapse.
+#' @param discretize_method Method used to place bin breaks: "width" for
+#' fixed binwidths, "ntile" for quantile bins, or "cluster" for univariate
+#' k-means clustering (set a seed before calling for reproducible clusters).
+#'
+#' @return A list with the discretized `synth_data` and `conf_data`.
+#'
+.discretize_k_marginal_vars <- function(
+  synth_data,
+  conf_data,
+  vars,
+  bins,
+  discretize_method
+) {
+  if (!(is.numeric(bins) && length(bins) == 1 && !is.na(bins) &&
+    bins >= 2 && bins == floor(bins))) {
+    stop("`bins` must be a single integer >= 2")
+  }
+
+  numeric_vars <- vars[
+    purrr::map_lgl(.x = vars, .f = \(v) is.numeric(conf_data[[v]]))
+  ]
+
+  for (var in numeric_vars) {
+    # cut points derive from observed values; missing values follow na.rm
+    # like every other variable, becoming an NA bin or dropped rows
+    conf_values <- conf_data[[var]][!is.na(conf_data[[var]])]
+
+    if (!all(is.finite(conf_values)) || length(conf_values) == 0) {
+      stop(
+        "observed numeric values must be finite to discretize; `", var,
+        "` is not"
+      )
+    }
+
+    if (dplyr::n_distinct(conf_values) < bins) {
+      stop(
+        "`", var, "` has fewer distinct confidential values than `bins`"
+      )
+    }
+
+    # interior cut points always derive from the confidential data; each
+    # method yields bins - 1 of them
+    cut_points <- switch(
+      EXPR = discretize_method,
+      width = seq(
+        from = min(conf_values),
+        to = max(conf_values),
+        length.out = bins + 1
+      )[2:bins],
+      ntile = stats::quantile(
+        x = conf_values,
+        probs = seq(from = 0, to = 1, length.out = bins + 1),
+        names = FALSE
+      )[2:bins],
+      cluster = {
+        centers <- sort(
+          stats::kmeans(x = conf_values, centers = bins)$centers[, 1]
+        )
+
+        (centers[-1] + centers[-length(centers)]) / 2
+      }
+    )
+
+    # outer bins extend to +/-Inf so out-of-range synthetic values land in
+    # edge bins; unique() collapses ties from skewed quantiles
+    breaks <- unique(c(-Inf, cut_points, Inf))
+
+    if (length(breaks) - 1 < bins) {
+      warning(
+        "`", var, "` was discretized into ", length(breaks) - 1,
+        " bins instead of ", bins, " because of tied cut points"
+      )
+    }
+
+    synth_data[[var]] <- cut(x = synth_data[[var]], breaks = breaks)
+    conf_data[[var]] <- cut(x = conf_data[[var]], breaks = breaks)
+  }
+
+  return(list(synth_data = synth_data, conf_data = conf_data))
+}

R/select_k_marginal_combos.R

@@ -0,0 +1,47 @@
+#' @title Select the variable combinations for the k-marginals metric
+#'
+#' @description Enumerates all unique k-combinations of the supplied
+#' variables and, when the total exceeds `n_marginals`, samples the
+#' combinations down to the cap while always retaining combinations that
+#' contain a priority variable.
+#'
+#' @param shared_vars Character vector of variables available for marginals.
+#' @param k Scalar order of the k-marginal.
+#' @param n_marginals Single integer target maximum for the number of
+#' variable combinations, or `Inf` for no cap. All priority combinations are
+#' retained even when they alone exceed the target.
+#' @param priority_vars Optional character vector of variable names whose
+#' combinations always survive sampling. Defaults to `NULL`.
+#'
+#' @return A character matrix with one row per selected combination and `k`
+#' columns.
+#'
+.select_k_marginal_combos <- function(
+  shared_vars,
+  k,
+  n_marginals,
+  priority_vars = NULL
+) {
+  kmarginals_vars <- t(utils::combn(x = shared_vars, m = k))
+
+  # sample combinations down to n_marginals, always keeping combinations that
+  # contain a priority variable
+  if (nrow(kmarginals_vars) > n_marginals) {
+    is_priority <- apply(
+      X = kmarginals_vars,
+      MARGIN = 1,
+      FUN = \(vars) any(vars %in% priority_vars)
+    )
+
+    n_sampled <- min(max(n_marginals - sum(is_priority), 0), sum(!is_priority))
+
+    sampled_rows <- sample(x = which(!is_priority), size = n_sampled)
+
+    kmarginals_vars <- kmarginals_vars[
+      sort(c(which(is_priority), sampled_rows)), ,
+      drop = FALSE
+    ]
+  }
+
+  return(kmarginals_vars)
+}

R/stratify_k_marginals.R

@@ -0,0 +1,105 @@
+#' @title Compute stratified k-marginals results
+#'
+#' @description Splits the synthetic and confidential data into the strata
+#' observed in the confidential data, computes per-cell differences and
+#' per-combination MabsDDs within each stratum, and rolls the per-stratum
+#' scores up weighted by each stratum's confidential share.
+#'
+#' @param synth_data A tibble with synthetic data.
+#' @param conf_data A tibble with confidential data.
+#' @param combos A character matrix with one row per variable combination.
+#' @param group_by Character vector of grouping variable names.
+#' @param weight_var Optional character name of a numeric sample-weight
+#' column; when set, cell proportions and stratum shares are weight shares
+#' instead of row shares. Defaults to `NULL` (unweighted).
+#' @param na.rm A logical for dropping rows with a missing value from each
+#' marginal that uses the affected variable.
+#'
+#' @return A list with `score` (the share-weighted mean of per-stratum
+#' scores), `marginals` and `cells` (stacked per-stratum tables with the
+#' grouping columns, worst first), and `group_scores` (one row per stratum
+#' with its share and score, worst first).
+#'
+.stratify_k_marginals <- function(
+  synth_data,
+  conf_data,
+  combos,
+  group_by,
+  weight_var = NULL,
+  na.rm = FALSE
+) {
+  # strata are defined by the confidential data; a stratum with no synthetic
+  # rows scores against all-zero synthetic proportions. Shares are the
+  # confidential row (or weight) share of each stratum, computed once
+  conf_totals <- if (is.null(weight_var)) {
+    rep(1, nrow(conf_data))
+  } else {
+    conf_data[[weight_var]]
+  }
+
+  strata <- conf_data |>
+    dplyr::mutate(.stratum_total = conf_totals) |>
+    dplyr::summarize(
+      .share = sum(.data$.stratum_total),
+      .by = dplyr::all_of(group_by)
+    ) |>
+    dplyr::mutate(.share = .data$.share / sum(.data$.share))
+
+  per_stratum <- purrr::map(
+    .x = seq_len(nrow(strata)),
+    .f = \(i) {
+      stratum <- strata[i, group_by, drop = FALSE]
+
+      synth_g <- dplyr::semi_join(synth_data, stratum, by = group_by)
+      conf_g <- dplyr::semi_join(conf_data, stratum, by = group_by)
+
+      cells_g <- .compute_marginal_cells(
+        synth_data = synth_g,
+        conf_data = conf_g,
+        combos = combos,
+        weight_var = weight_var,
+        na.rm = na.rm,
+        allow_empty_synth = TRUE
+      )
+
+      marginals_g <- cells_g |>
+        dplyr::summarize(madd = mean(.data$abs_diff), .by = "variables")
+
+      list(
+        cells = dplyr::bind_cols(stratum, cells_g),
+        marginals = dplyr::bind_cols(stratum, marginals_g),
+        group_scores = dplyr::bind_cols(
+          stratum,
+          tibble::tibble(
+            share = strata$.share[i],
+            score = (1 - mean(marginals_g$madd)) * 1000
+          )
+        )
+      )
+    }
+  )
+
+  cells <- purrr::list_rbind(purrr::map(per_stratum, "cells")) |>
+    dplyr::arrange(dplyr::desc(.data$abs_diff))
+
+  marginals <- purrr::list_rbind(purrr::map(per_stratum, "marginals")) |>
+    dplyr::arrange(dplyr::desc(.data$madd))
+
+  group_scores <- purrr::list_rbind(
+    purrr::map(per_stratum, "group_scores")
+  ) |>
+    dplyr::arrange(.data$score)
+
+  # per-stratum scores roll up weighted by confidential shares, so small
+  # strata surface in group_scores without dominating the headline
+  score <- sum(group_scores$share * group_scores$score)
+
+  return(
+    list(
+      score = score,
+      marginals = marginals,
+      cells = cells,
+      group_scores = group_scores
+    )
+  )
+}