Propensity Score Matching

The code that produced every figure in this report, exactly as it ran. Fingerprint ad3c3c0e5f3d50c1. The same code on the same data gives the same report.

`standard_propensity_matching_v2` <- function(pf) {
  `%||%` <- function(a, b) if (!is.null(a)) a else b
  #' Readable figures (LAT-3180, LAT-3181): whole numbers from a thousand up, one decimal from a hundred, two from
  #' one, three significant figures below one. A cell carries what the value needs, not what R prints.
  tidy <- function(x) {
    x <- as.numeric(x)
    ifelse(is.na(x), NA_real_,
      ifelse(abs(x) >= 1000, round(x, 0),
        ifelse(abs(x) >= 100, round(x, 1),
          ifelse(abs(x) >= 1, round(x, 2), signif(x, 3)))))
  }
  p_text <- function(p) if (is.na(p)) "" else if (p < 1e-4) "<0.0001" else format(signif(p, 3))
  inputs <- pf$taskList$inputs
  params <- inputs$module_parameters %||% list()
  question <- (inputs$userContext %||% list())$objective %||%
    "Did the treated group do better because of the treatment, or were they different to begin with?"

  #' ## Column mapping
  #' `outcome` (numeric), `treatment` (two values) and a `covariate_N` series (at least one, numeric or categorical).
  #' Semantic names inside; the customer's own headers live in `col_map`.
  col_map <- inputs$column_mapping %||% list()
  df <- renderObject.taskFunction.init(inputs, col_map)   # df has SEMANTIC names
  human <- function(sem) {
    v <- col_map[[sem]]
    if (is.null(v) || !nzchar(as.character(v))) sem else as.character(v)
  }
  for (sem in c("outcome", "treatment"))
    if (!(sem %in% names(df))) stop(sprintf("column_mapping must map '%s' (%s was not found).", sem, human(sem)))
  outcome_h <- human("outcome"); treat_h <- human("treatment")
  cov_cols <- grep("^covariate_[0-9]+$", names(df), value = TRUE)
  cov_cols <- cov_cols[order(as.integer(sub("^covariate_", "", cov_cols)))]
  if (length(cov_cols) < 1) stop("column_mapping must map at least one covariate (covariate_1, covariate_2, ...): the pre-existing differences to adjust for.")
  n_in <- nrow(df)

  #' ## The columns this tool did not look at, read from the raw rows (init narrows `df` to the mapped columns)
  raw_names <- local({
    ds <- inputs$dataset %||% inputs$df
    if (is.data.frame(ds)) return(names(ds))
    if (is.list(ds) && length(ds) > 0) {
      rows <- ds[seq_len(min(length(ds), 50))]
      nm <- unique(unlist(lapply(rows, function(r) if (is.list(r)) names(r) else NULL)))
      if (length(nm)) return(nm)
      if (!is.null(names(ds)) && all(nzchar(names(ds)))) return(names(ds))
    }
    character(0)
  })
  mapped_actual <- unique(as.character(unlist(col_map)))
  ignored_cols <- setdiff(raw_names, unique(c(mapped_actual, make.names(mapped_actual))))

  #' ## Parameters
  caliper_sd <- suppressWarnings(as.numeric(params$caliper %||% 0.2))
  if (is.na(caliper_sd) || caliper_sd < 0.05 || caliper_sd > 1) stop("module_parameters$caliper must be a number from 0.05 to 1 (standard deviations of the logit propensity score)")

  #' ## Data preparation
  #' The treatment must take exactly two values; the treated label is 1/yes/true/treated/member/exposed when one value
  #' matches, else the second value in sorted order. Covariates: numeric when 95% of their values parse and they take
  #' more than two values, else categorical; a covariate with one value, or one that names each row, is dropped and
  #' named. A row missing the outcome, the treatment or any kept covariate is excluded and counted. At least 40 rows
  #' and 10 in each group are required.
  y <- suppressWarnings(as.numeric(as.character(df$outcome)))
  tr <- trimws(as.character(df$treatment)); tr[is.na(tr) | tr == "" | tolower(tr) %in% c("na", "n/a", "null", "nan")] <- NA
  lv <- sort(unique(stats::na.omit(tr)))
  if (length(lv) != 2) stop(sprintf("%s takes %d distinct value%s (%s); propensity matching compares exactly two groups, such as member/non-member, yes/no or 1/0.",
                                    treat_h, length(lv), if (length(lv) == 1) "" else "s", paste(utils::head(lv, 5), collapse = ", ")))
  hits <- lv[grepl("^(1|yes|true|treated|member|exposed|y|t)$", tolower(lv))]
  treated_label <- if (length(hits) == 1) hits else lv[2]; control_label <- setdiff(lv, treated_label)[1]
  t01 <- as.integer(tr == treated_label)
  dat <- data.frame(y = y, t = t01)
  kept <- character(0); cov_kind <- character(0); dropped <- character(0); why <- character(0)
  for (cc in cov_cols) {
    v <- df[[cc]]; num <- suppressWarnings(as.numeric(as.character(v))); ok <- !is.na(v) & nzchar(trimws(as.character(v)))
    if (sum(ok) == 0) { dropped <- c(dropped, human(cc)); why <- c(why, "empty"); next }
    if (sum(!is.na(num)) / sum(ok) >= 0.95 && length(unique(num[!is.na(num)])) > 2) {
      dat[[cc]] <- num; cov_kind <- c(cov_kind, "numeric")
    } else {
      f <- trimws(as.character(v)); f[is.na(f) | f == ""] <- NA
      u <- length(unique(stats::na.omit(f)))
      if (u < 2) { dropped <- c(dropped, human(cc)); why <- c(why, "one value"); next }
      if (u > 0.5 * sum(!is.na(f)) && u > 20) { dropped <- c(dropped, human(cc)); why <- c(why, "names each row"); next }
      dat[[cc]] <- factor(f); cov_kind <- c(cov_kind, "categorical")
    }
    kept <- c(kept, cc)
  }
  if (length(kept) < 1) stop(sprintf("No usable covariate remained (%s).", paste(paste(dropped, why, sep = ": "), collapse = "; ")))
  complete <- stats::complete.cases(dat); n_missing <- sum(!complete); dat <- dat[complete, , drop = FALSE]; n <- nrow(dat)
  n_t <- sum(dat$t == 1); n_c <- sum(dat$t == 0)
  if (n < 40) stop(sprintf("Only %d complete rows; a matched comparison needs at least 40.", n))
  if (n_t < 10 || n_c < 10) stop(sprintf("%s has %d %s rows and %d %s rows; each group needs at least 10.", treat_h, n_t, treated_label, n_c, control_label))

  #' ## The propensity score: logistic regression of the treatment on the covariates
  glm_warn <- character(0)
  ps_model <- withCallingHandlers(stats::glm(stats::reformulate(kept, response = "t"), data = dat, family = stats::binomial()),
                                  warning = function(w) { glm_warn <<- c(glm_warn, conditionMessage(w)); invokeRestart("muffleWarning") })
  ps <- as.numeric(stats::fitted(ps_model)); lps <- stats::qlogis(pmin(pmax(ps, 1e-6), 1 - 1e-6))
  sep_frac <- mean(ps < 1e-3 | ps > 1 - 1e-3)
  nonconv <- any(grepl("did not converge|fitted probabilities numerically 0 or 1", glm_warn))
  t_idx <- which(dat$t == 1); c_idx <- which(dat$t == 0)
  support_share <- mean(ps[t_idx] >= min(ps[c_idx]) & ps[t_idx] <= max(ps[c_idx]))

  #' ## Matching: one nearest control per treated row on the logit score, within the caliper, without replacement,
  #' treated rows visited from the highest score down (the v1's rule, kept so the two agree pair for pair)
  sd_lps <- stats::sd(lps); caliper <- if (is.finite(sd_lps)) caliper_sd * sd_lps else 0
  avail <- rep(TRUE, length(c_idx)); pt <- integer(0); pc <- integer(0)
  for (ti in t_idx[order(-lps[t_idx])]) {
    d <- abs(lps[c_idx] - lps[ti]); d[!avail] <- Inf; j <- which.min(d)
    if (length(j) == 1 && is.finite(d[j]) && d[j] <= caliper) { pt <- c(pt, ti); pc <- c(pc, c_idx[j]); avail[j] <- FALSE }
  }
  n_pairs <- length(pt); n_unmatched <- n_t - n_pairs
  if (n_pairs < 10) stop(sprintf("Only %d matched pair%s could be formed within the caliper; the %s and %s groups are too different on the mapped covariates for a reliable match.", n_pairs, if (n_pairs == 1) "" else "s", treated_label, control_label))

  #' ## The two estimates and the balance
  yv <- dat$y
  naive <- mean(yv[t_idx]) - mean(yv[c_idx]); nt <- tryCatch(stats::t.test(yv[t_idx], yv[c_idx]), error = function(e) NULL)
  naive_ci <- if (!is.null(nt)) as.numeric(nt$conf.int) else c(NA, NA); naive_p <- if (!is.null(nt)) nt$p.value else NA_real_
  diffs <- yv[pt] - yv[pc]; att <- mean(diffs); mt <- tryCatch(stats::t.test(diffs), error = function(e) NULL)
  att_ci <- if (!is.null(mt)) as.numeric(mt$conf.int) else c(NA, NA); att_p <- if (!is.null(mt)) mt$p.value else NA_real_
  smd_num <- function(a, b) { s <- sqrt((stats::var(a) + stats::var(b)) / 2); if (!is.finite(s) || s == 0) { if (isTRUE(all.equal(mean(a), mean(b)))) 0 else NA_real_ } else (mean(a) - mean(b)) / s }
  smd_cat <- function(a, b) { l <- union(unique(as.character(a)), unique(as.character(b))); if (!length(l)) return(NA_real_); max(abs(vapply(l, function(x) mean(as.character(a) == x) - mean(as.character(b) == x), numeric(1)))) }
  hn <- vapply(kept, human, character(1))
  bal <- do.call(rbind, lapply(seq_along(kept), function(i) {
    x <- dat[[kept[i]]]
    if (is.numeric(x)) { b <- smd_num(x[t_idx], x[c_idx]); a <- smd_num(x[pt], x[pc]) } else { b <- smd_cat(x[t_idx], x[c_idx]); a <- smd_cat(x[pt], x[pc]) }
    data.frame(covariate = hn[i], smd_before = round(b, 3), smd_after = round(a, 3), balanced = ifelse(is.na(a), "unknown", ifelse(abs(a) < 0.10, "yes", "no")), stringsAsFactors = FALSE)
  }))
  rownames(bal) <- NULL
  worst_after <- if (any(is.finite(bal$smd_after))) max(abs(bal$smd_after), na.rm = TRUE) else NA_real_
  n_balanced <- sum(bal$balanced == "yes")

  #' ## The frames for the places
  int_df <- data.frame(estimate = c("Naive difference", "Matched difference"), effect = tidy(c(naive, att)), ci_low = tidy(c(naive_ci[1], att_ci[1])), ci_high = tidy(c(naive_ci[2], att_ci[2])), stringsAsFactors = FALSE)
  bba <- data.frame(covariate = rep(bal$covariate, 2), abs_smd = round(abs(c(bal$smd_before, bal$smd_after)), 3), stage = rep(c("Before matching", "After matching"), each = nrow(bal)), stringsAsFactors = FALSE)
  ov <- data.frame(propensity = round(ps, 4), group = ifelse(dat$t == 1, treated_label, control_label), stringsAsFactors = FALSE)
  if (nrow(ov) > 1000) { set.seed(20260920); ov <- ov[sort(sample.int(nrow(ov), 1000)), , drop = FALSE] }
  pairs_df <- data.frame(control_outcome = tidy(yv[pc]), treated_outcome = tidy(yv[pt]), stringsAsFactors = FALSE)
  if (nrow(pairs_df) > 1000) { set.seed(20260920); pairs_df <- pairs_df[sort(sample.int(nrow(pairs_df), 1000)), , drop = FALSE] }
  eff_df <- data.frame(estimate = c("Naive difference", "Matched difference (ATT)"), effect = tidy(c(naive, att)), ci_low = tidy(c(naive_ci[1], att_ci[1])), ci_high = tidy(c(naive_ci[2], att_ci[2])),
                       p_value = c(p_text(naive_p), p_text(att_p)), rows = c(n_t + n_c, 2L * n_pairs), stringsAsFactors = FALSE)

  #' ## Assumption checks (LAT-3138)
  matched_share <- n_pairs / n_t
  checks_df <- data.frame(
    check = c("Enough matched pairs", "Treated rows matched", "Common support", "Balance after matching", "Propensity model fit", "Outcome kept out of the score"),
    statistic = c(sprintf("%d pairs", n_pairs), sprintf("%d of %d %s rows matched (%s%%)", n_pairs, n_t, treated_label, format(round(100 * matched_share, 1))),
                  sprintf("%s%% of %s rows inside the %s score range", format(round(100 * support_share, 1)), treated_label, control_label),
                  sprintf("worst absolute SMD after = %s; %d of %d covariates under 0.10", format(round(worst_after, 3)), n_balanced, nrow(bal)),
                  if (nonconv) "convergence warning" else sprintf("%s%% of scores at 0 or 1", format(round(100 * sep_frac, 1))),
                  sprintf("%d covariate%s in the score; the outcome is not one of them", length(kept), if (length(kept) == 1) "" else "s")),
    p_value = c("", "", "", "", "", ""),
    verdict = c(if (n_pairs >= 50) "holds" else if (n_pairs >= 20) "strained" else "violated",
                if (matched_share >= 0.9) "holds" else if (matched_share >= 0.6) "strained" else "violated",
                if (support_share >= 0.9) "holds" else if (support_share >= 0.7) "strained" else "violated",
                if (is.na(worst_after)) "strained" else if (worst_after < 0.10) "holds" else if (worst_after < 0.25) "strained" else "violated",
                if (nonconv || sep_frac > 0.5) "violated" else if (sep_frac > 0.1) "strained" else "holds", "holds"),
    note = c("few pairs leave the matched interval wide", "the matched estimate is for the treated rows that found a match; the unmatched are the least comparable",
             "treated rows outside the controls' range have no comparable control", "a covariate still above 0.10 after matching leaves part of the raw gap unexplained",
             "scores at 0 or 1 mean the covariates predict the treatment almost perfectly, so few pairs are comparable",
             "a score that used the outcome would match on the very thing being compared"),
    stringsAsFactors = FALSE)

  #' ## Method, assumptions, answer
  change_word <- if (abs(att) < 0.9 * abs(naive)) "shrinks" else if (abs(att) > 1.1 * abs(naive)) "grows" else "changes little"
  excluded <- c(if (n_missing > 0) sprintf("%d row%s missing the outcome, the treatment or a covariate", n_missing, if (n_missing == 1) "" else "s"),
                if (length(dropped)) sprintf("%d covariate%s dropped (%s)", length(dropped), if (length(dropped) == 1) "" else "s", paste(paste(dropped, why, sep = ": "), collapse = "; ")),
                if (n_unmatched > 0) sprintf("%d %s row%s unmatched within the caliper", n_unmatched, treated_label, if (n_unmatched == 1) "" else "s"))
  method <- paste0(
    "Propensity score matching of ", outcome_h, " between ", treat_h, " = ", treated_label, " (treated, ", n_t, " rows) and ", control_label, " (control, ", n_c, " rows) over ", n, " complete rows of ", n_in,
    "; the propensity score from a logistic regression of the treatment on ", paste(hn, collapse = ", "),
    "; one nearest control per treated row on the logit score within a caliper of ", caliper_sd, " SD (", format(round(caliper, 4)), "), without replacement, treated rows visited from the highest score down; ",
    n_pairs, " pairs formed; the matched difference (the effect on the treated) as the mean of the pair differences with a paired-t 95% interval; the naive difference with a Welch interval; balance as standardized mean differences (numeric) or the largest category-share difference (categorical); the gap ", change_word, " after matching",
    if (length(excluded)) paste0("; excluded: ", paste(excluded, collapse = "; ")) else "",
    if (length(ignored_cols)) paste0("; not used: ", paste(utils::head(ignored_cols, 12), collapse = ", "), " (not mapped)") else "", ".")
  assumptions <- list(
    "Every factor that made a unit more likely to be treated AND affected the outcome is among the mapped covariates; an unmeasured one biases the matched difference.",
    "The covariates were measured before the treatment; a covariate the treatment changed would remove part of its own effect.",
    "The matched difference is the effect on the TREATED units that found a match, not on everyone; the unmatched treated were the least comparable.",
    sprintf("The treated label was read as '%s'; if that is the wrong side, the sign of every difference flips.", treated_label),
    "Matching balances what was measured; it is not randomisation, and the estimate is causal only if the first assumption holds.")
  answer <- list(att = tidy(att), att_low = tidy(att_ci[1]), att_high = tidy(att_ci[2]), p_value = p_text(att_p), naive_difference = tidy(naive),
                 naive_low = tidy(naive_ci[1]), naive_high = tidy(naive_ci[2]), pairs = n_pairs, treated = n_t, control = n_c,
                 treated_matched_pct = round(100 * matched_share, 1), worst_smd_after = round(worst_after, 3), covariates_balanced = n_balanced,
                 covariates = length(kept), treated_label = treated_label, control_label = control_label, gap = change_word)

  results <- list()
  #' The verdict and the headline are NOT places of a library tool (LAT-3130): the last mile writes them.
  results$summary_metrics <- place_metric(list(att = tidy(att), att_low = tidy(att_ci[1]), att_high = tidy(att_ci[2]), naive_difference = tidy(naive),
    pairs = n_pairs, treated_matched_pct = round(100 * matched_share, 1), worst_smd_after = round(worst_after, 3)), lead = "att", place = "summary_metrics")
  results$effect_interval <- place_interval(int_df, term = "estimate", value = "effect", low = "ci_low", high = "ci_high", place = "effect_interval")
  results$balance_before_after <- place_comparison(bba, category = "covariate", value = "abs_smd", series = "stage", place = "balance_before_after")
  results$propensity_overlap <- place_distribution(ov, x = "propensity", series = "group", place = "propensity_overlap")
  results$balance_table <- place_table(bal, place = "balance_table")
  if (nrow(pairs_df) >= 10) {
    results$matched_pairs <- place_relationship(pairs_df, x = "control_outcome", y = "treated_outcome", place = "matched_pairs")
  } else {
    results$matched_pairs <- place_dropped("fewer than 10 matched pairs, too few to draw", place = "matched_pairs")
  }
  results$effect_table <- place_table(eff_df, place = "effect_table")
  results$assumption_checks <- place_table(checks_df, place = "assumption_checks")
  results$matching_method <- place_method(method = method, n_in = n_in, n_used = 2L * n_pairs, assumptions = assumptions, excluded = as.list(excluded), place = "matching_method")

  objects <- list()   # filled by the object layer, not here
  list(answer = answer, method = method, n = n, results = results, objects = objects,
       json_output = list(answer = answer, method = method, n = n))
}