High-Dimensional Propensity Score (hdPS)

In plain language

The high-dimensional propensity score (hdPS) is a method that automatically scans hundreds of diagnosis, procedure, and drug codes recorded in a patient's claims history BEFORE they start a new medication, and uses those codes as stand-ins (proxies) for health factors the data never directly measured -- like how sick a patient really is or how often they seek care. Instead of a researcher hand-picking a short list of confounders, the algorithm ranks every candidate code by how differently it appears between the two treatment groups AND how strongly it relates to the outcome, then builds a single summary score from the top-ranked codes. That score is then used to match or weight patients so the two comparison groups look more similar on all those hidden health differences -- though hdPS cannot fix confounding from factors that left no footprint at all in the claims data.

The high-dimensional propensity score (hdPS) automates the selection of confounder proxies in large healthcare databases. Investigator-specified propensity scores depend on a short list of variables a human thought to measure; in claims and EHR data the strongest confounders (frailty, disease severity, functional status, access, health-seeking behavior, socioeconomic position) are never coded directly. hdPS (Schneeweiss et al. 2009) exploits the fact that the thousands of diagnosis, procedure, drug, and utilization codes recorded before exposure are noisy proxies for those latent factors, and lets the data rank which proxies to include. The canonical algorithm has seven steps: (1) specify the cohort, exposure (typically a new-user, active-comparator contrast), outcome, and a pre-exposure lookback (commonly 365 days); (2) within each data dimension — inpatient diagnoses, outpatient diagnoses, procedures, dispensed drugs, and utilization counts — identify the most prevalent codes; (3) for each candidate code, generate up to three binary recurrence covariates (once / sporadic / frequent relative to the cohort median); (4) rank every candidate by its potential confounding impact using the Bross multiplicative bias formula, which combines the covariate's prevalence in the exposed and unexposed and its marginal association with the outcome; (5) select the top k per dimension (Schneeweiss's default selects ~200 across dimensions); (6) fit a logistic PS on the selected proxies plus force-included investigator confounders (age, sex, calendar time, indication); (7) apply the PS for matching, weighting, stratification, or SMR, then assess balance, trim non-overlap, and run the outcome model with robust variance.

Core estimand distinction

hdPS is a covariate-selection and PS-construction method, not an estimator — it does not itself define a target estimand. The estimand is fixed by how the resulting score is used: 1:1 nearest-neighbor matching on the logit-PS targets the ATT (effect in the treated); IPTW targets the ATE; overlap (Li) weights target the ATO (effect in the equipoise population, with exact finite-sample balance); stratification/SMR weighting targets an ATT-like quantity in the standard population. This must be pre-specified — switching from IPTW to matching after seeing imbalance silently changes the population the result generalizes to. Critically, hdPS only reduces confounding bias for confounders that have proxies in the data; it has no leverage on confounders with no measured footprint, and it does not address selection bias, measurement error, or model misspecification of the outcome stage.

Pros, cons, and trade-offs

- vs investigator-specified PS / multivariable regression: hdPS empirically surfaces proxies a human would never enumerate (e.g., frequency of ophthalmology visits proxying diabetic retinopathy severity), and in plasmode simulations and the small-sample work of Rassen et al. (2011) it reduces residual confounding and improves balance on measured covariates. Cost: the selection is opaque, it can pull in instruments (predict exposure, not outcome — these amplify residual bias and inflate variance) or colliders/mediators if the lookback leaks post-exposure information, and the selected set is database- and era-specific. Prefer hdPS in large claims/EHR studies where measured confounding is known to be incomplete; keep an investigator PS as a transparent reference analysis. - vs instrumental-variable (IV) methods: hdPS stays inside the measured data and needs no valid instrument (instruments are rare and nearly impossible to verify in pharmacoepi). Cost: hdPS cannot touch confounding with no proxy in the data, whereas a valid IV can — but IV targets a different (local average treatment) effect under strong, untestable assumptions. Prefer hdPS when the database is proxy-rich and no credible instrument exists. - vs disease-risk scores (DRS) / doubly-robust (TMLE, AIPW): hdPS is a single-model exposure-side approach; DRS models the outcome side and doubly-robust methods combine both for protection against one-sided misspecification. hdPS-selected covariates are routinely fed into TMLE/super-learner pipelines (Schneeweiss 2017 frames hdPS as one screening option among LASSO and high-dimensional DRS). Prefer plain hdPS for transparent, regulator-facing comparative-safety work; escalate to doubly-robust when both stages risk misspecification and the team can defend the added complexity. - vs negative controls / quantitative bias analysis (QBA): these are complements, not substitutes. hdPS reduces bias up front; negative-control outcomes and E-value/QBA detect and bound the residual bias hdPS cannot remove. A defensible analysis uses both.

When to use

Large administrative claims, EHR, or linked databases with a new-user, active-comparator design where (a) the outcome is reasonably common, (b) thousands of pre-exposure codes are available as proxies, and (c) key clinical confounders are unmeasured or poorly captured. It is the de facto standard confounding-control layer in FDA Sentinel and distributed-network safety studies, and is most valuable precisely when investigator covariate lists are thin.

When NOT to use — and when it is actively misleading or dangerous

- Small cohorts or rare outcomes. With few exposed or few events, the Bross ranking is unstable, PS separation/positivity violations appear, and selecting hundreds of sparse covariates overfits. Below ~25–50 events the algorithm can worsen bias; fall back to a parsimonious investigator PS or exact matching (Rassen 2011 documents the small-sample failure mode and proposed remedies). - The lookback can leak post-exposure information. If candidate codes are drawn from a window that includes or straddles the index date, hdPS will happily select mediators or colliders (e.g., early treatment toxicities), inducing collider-stratification bias and over-adjustment — the bias gets built into the score. Hard-stop candidate extraction strictly before time zero. - Strong instruments are abundant. In settings dominated by exposure-predictive-but-outcome-irrelevant codes (formulary, provider, or plan markers), unconstrained selection amplifies residual confounding and variance. Screen for, and exclude, near-instruments. - Cross-database transport. A proxy set tuned to one payer/era does not transport; re-running the algorithm per database is mandatory, not optional. Treating a frozen code list as portable produces silently biased estimates in the new source.

Data-source operational depth

- Claims (FFS): The natural habitat — fee-for-service Medicare and commercial claims pay per service, so diagnosis/procedure/NDC capture is dense and the proxies are rich. Build dimensions explicitly from claim type (facility/inpatient dx and procedures, professional/outpatient dx and procedures, Part D / pharmacy NDCs, DME). Use `service_date` to enforce a strictly pre-index window. Failure mode: codes on rule-out claims inflate apparent comorbidity; the once/sporadic/frequent recurrence coding partially mitigates this. - Claims (Medicare Advantage): MA-only person-time lacks complete FFS claims (capitated encounters are under-reported for non-risk services), so candidate prevalence is distorted and "absence of a code" is missingness, not a true negative. Conversely, HCC risk-adjustment coding is enriched under MA, so a handful of HCC-driving diagnoses can dominate selection. Exclude MA-only spans, or treat MA as a separate stratum/dimension and never pool raw candidate prevalences across FFS and MA. - EHR: Add clinical dimensions — labs (LOINC + abnormal-flag), vitals, problem-list entries, NLP-derived concepts — which often out-proxy billing codes for severity. But EHR is visit-driven: encounter frequency is simultaneously a powerful proxy and a potential collider if care-seeking is itself affected by the (pre-index) disease, and a patient who leaves the system is differentially unobserved. Link to claims for complete pharmacy/procedure history before trusting the proxy set. - Linked claims–EHR–registry: The richest substrate (registry staging/biomarkers + EHR labs + claims completeness), but linkage selects the linkable subset and introduces order/fill/service date discrepancies that must be reconciled before the pre-index window is fixed. Differential competing risks matter in elderly claims cohorts: if one arm has higher background mortality, naive hdPS-adjusted cause-specific analyses can mislead — pair the score with an appropriate competing-risks outcome model.

Worked claims example

Question: incident diabetic ketoacidosis with SGLT2 inhibitors vs DPP-4 inhibitors among adults with type 2 diabetes in a 100% Medicare FFS (Parts A/B/D) + commercial database. (1) Cohort: age ≥18, ≥2 T2DM diagnoses, 365 days continuous FFS-observable enrollment before the first study fill (`fill_date`), new-user washout = no prior SGLT2i/DPP-4 NDC in that window; exclude MA-only person-time. Index date = first qualifying fill; arm = `ndc` dispensed that day. (2) Force-include investigator confounders: age, sex, calendar year, baseline HbA1c-proxy, prior insulin, diabetes-complication count. (3) Candidate generation, strictly in `[index_date − 365, index_date)`: in each dimension take the 200 most prevalent codes; for each, create once/sporadic/frequent recurrence covariates → ~3,000 candidates. (4) Bross ranking: for a candidate such as "frequent (>median) ophthalmology E/M visits" (CPT 9920x/9921x), with prevalence ~0.18 among DPP-4 initiators vs ~0.11 among SGLT2i initiators and an outcome relative risk ~1.6, the multiplicative bias term is large, so this retinopathy/severity proxy ranks near the top; a near-instrument like "mail-order pharmacy flag" (predicts arm, not DKA) ranks low and is dropped. (5) Select top k=50 per dimension (~250 covariates), fit the PS by logistic regression on selected + forced terms. (6) Apply 1:1 logit-PS matching with a 0.2-SD caliper (estimand = ATT), confirm standardized mean differences <0.1 on selected and forced covariates, trim non-overlap. (7) Outcome model: Cox or pooled logistic on the matched set with robust variance; censor at disenrollment, death, end of data, and (as-treated) last `days_supply` end + grace period. (8) Sensitivity: vary k (100/250/500), the dimension split, the recurrence cut, an investigator-PS reference, and a negative-control outcome (e.g., influenza vaccination) to detect residual confounding.

Interpreting the output

Consider the SGLT2i versus DPP-4i DKA study above. After incorporating the three selected proxy codes — frequent office visits, diabetic retinopathy, and insulin fills — alongside investigator-specified covariates in the hdPS, suppose the analysis reports HR = 1.38 (95% CI 1.08–1.75) for diabetic ketoacidosis.

Formal interpretation: In the hdPS-reweighted pseudo-population, the instantaneous rate of DKA was 38% higher among SGLT2i initiators than DPP-4i initiators (HR 1.38, 95% CI 1.08–1.75). This estimates the ATT (or ATE, depending on weight type) in the weighted pseudo-population. The central identification assumption is no unmeasured confounding given all included covariates — both investigator-specified variables and the empirically selected proxy codes. The proxy codes are data-driven surrogates for disease severity and care-seeking behavior not directly coded; they reduce but cannot eliminate residual confounding from factors leaving no trace in pre-treatment claims patterns. Positivity must also hold: every patient must have had a plausible probability of initiating either drug class.

Practical interpretation: After accounting for measurable differences in diabetes burden — using both clinical covariates and the claims-derived proxies that hdPS surfaced — SGLT2i users developed DKA at a 38% higher rate. The proxy codes captured dimensions of disease burden invisible to a standard propensity model. However, proxies are associative signals, not the unmeasured confounders themselves; bias from unmeasured indication factors may persist. A sensitivity analysis varying k (100, 250, 500 proxy codes) and a negative-control outcome test should accompany the primary estimate to assess residual confounding.

Worked example

Scenario

A researcher is comparing SGLT2 inhibitors vs DPP-4 inhibitors for risk of diabetic ketoacidosis in Medicare claims. She cannot directly measure disease severity or how actively a patient manages their diabetes -- those factors are never coded. She runs hdPS on the 365-day pre-treatment claims window to let the data surface proxy codes for those unmeasured factors. The table below shows five candidate codes the algorithm evaluates, their prevalence in each treatment group, and whether they get selected.

Dataset

Five candidate proxy codes evaluated by hdPS from the pre-treatment lookback window (365 days before first fill). Prevalence = share of patients in that arm who had the code at least once. Outcome RR = how much more common diabetic ketoacidosis is among patients with vs without that code.

Steps

• For each candidate code, count how many patients in each arm had it at least once, sporadically (at or above the group median count), or frequently (above the median) during the 365-day lookback -- these are the three recurrence covariates per code.

• Apply the Bross formula: a code scores high when (a) its prevalence differs noticeably between the SGLT2i and DPP-4i groups AND (b) patients who have the code get diabetic ketoacidosis at a meaningfully higher rate than those who do not.

• Frequent outpatient visits (CPT 99213), diabetic retinopathy diagnosis (E11.319), and insulin use all have both unequal prevalence between arms AND a raised outcome RR -- so their Bross scores are high and they are selected as proxy confounders.

• The mail-order pharmacy flag has unequal prevalence (it predicts which drug was chosen) but an outcome RR near 1.0 -- it is a near-instrument that would amplify bias if included, so it is dropped.

• The counseling code has equal prevalence and no outcome link -- it is uninformative and dropped.

• The top-ranked codes across all dimensions (up to ~50 per dimension) are combined with investigator-specified confounders (age, sex, calendar year, diabetes severity) to fit a logistic regression predicting treatment arm.

• Each patient receives a propensity score from that model; patients are then matched 1-to-1 on those scores so each SGLT2i user is paired with a similarly-scored DPP-4i user, balancing both the selected proxies and the directly measured covariates.

Result

Three proxy codes (frequent office visits, retinopathy diagnosis, insulin use) are selected because they satisfy both Bross criteria; two codes are dropped (one near-instrument, one uninformative). The final hdPS model incorporates these proxies alongside forced investigator confounders, producing a propensity score used for 1:1 matching to estimate the ATT (the treatment effect in patients similar to those who actually initiated SGLT2i).

Runnable example

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression

LOOKBACK_DAYS = 365
TOP_N_PREVALENT = 200     # most prevalent codes screened per dimension (Schneeweiss step 2)
TOP_K_PER_DIM   = 50      # covariates retained per dimension after Bross ranking (step 5)

def _recurrence_covariates(codes, cohort):
    # Keep only pre-index codes inside the lookback window.
    c = codes.merge(cohort[["person_id", "index_date"]], on="person_id")
    c = c[(c["code_date"] < c["index_date"]) &
          (c["code_date"] >= c["index_date"] - pd.Timedelta(days=LOOKBACK_DAYS))]
    counts = c.groupby(["dimension", "person_id", "code"]).size().rename("n").reset_index()

    feats = []
    for dim, dim_df in counts.groupby("dimension"):
        prev = dim_df.groupby("code")["person_id"].nunique().sort_values(ascending=False)
        for code in prev.head(TOP_N_PREVALENT).index:
            sub = dim_df[dim_df["code"] == code]
            med = sub["n"].median()
            # once / sporadic(>=median) / frequent(>median): the three hdPS recurrence covariates.
            for label, mask in (("once", sub["n"] >= 1),
                                ("sporadic", sub["n"] >= max(med, 1)),
                                ("frequent", sub["n"] > med)):
                ids = set(sub.loc[mask, "person_id"])
                feats.append((f"{dim}|{code}|{label}", ids))
    return feats

def _bross_bias(flag, exposed, outcome):
    # Multiplicative confounding-bias multiplier (Bross/Schneeweiss): prevalence of the covariate
    # in exposed (P_C1) vs unexposed (P_C0) and the covariate-outcome relative risk (RR_CD).
    p_c1 = flag[exposed == 1].mean()
    p_c0 = flag[exposed == 0].mean()
    o1, o0 = outcome[flag == 1].mean(), outcome[flag == 0].mean()
    rr_cd = (o1 / o0) if o0 > 0 else 1.0
    bias = ((p_c1 * (rr_cd - 1) + 1) / (p_c0 * (rr_cd - 1) + 1))
    return abs(np.log(bias)) if bias > 0 else 0.0

def fit_hdps(cohort, codes, forced=None):
    cohort = cohort.set_index("person_id")
    exposed = cohort["exposed"]; outcome = cohort["outcome"]
    ranked = []  # (dimension, name, |log bias|, indicator series)
    for name, ids in _recurrence_covariates(codes, cohort.reset_index()):
        flag = pd.Series(cohort.index.isin(ids).astype(int), index=cohort.index)
        if flag.nunique() < 2:
            continue
        ranked.append((name.split("|")[0], name,
                       _bross_bias(flag, exposed, outcome), flag))

    selected = []
    for dim in {r[0] for r in ranked}:
        dim_cands = sorted((r for r in ranked if r[0] == dim), key=lambda r: r[2], reverse=True)
        selected.extend(dim_cands[:TOP_K_PER_DIM])

    X = pd.concat({name: flag for _, name, _, flag in selected}, axis=1)
    if forced is not None:
        X = X.join(forced)                      # age, sex, calendar year, indication, severity
    ps = LogisticRegression(max_iter=1000, C=1.0).fit(X, exposed).predict_proba(X)[:, 1]
    cohort["ps"] = ps
    return cohort[["exposed", "outcome", "ps"]], X

library(data.table)
LOOKBACK_DAYS  <- 365L
TOP_N_PREVALENT <- 200L
TOP_K_PER_DIM   <- 50L

bross_bias <- function(flag, exposed, outcome) {
  p_c1 <- mean(flag[exposed == 1]); p_c0 <- mean(flag[exposed == 0])
  o1 <- mean(outcome[flag == 1]);   o0 <- mean(outcome[flag == 0])
  rr_cd <- if (o0 > 0) o1 / o0 else 1
  bias  <- (p_c1 * (rr_cd - 1) + 1) / (p_c0 * (rr_cd - 1) + 1)
  if (bias > 0) abs(log(bias)) else 0
}

fit_hdps <- function(cohort, codes, forced = NULL) {
  setDT(cohort); setDT(codes)
  c <- codes[cohort[, .(person_id, index_date, exposed, outcome)], on = "person_id"]
  c <- c[code_date < index_date & code_date >= index_date - LOOKBACK_DAYS]
  cnt <- c[, .(n = .N), by = .(dimension, person_id, code, exposed, outcome)]

  # Screen the most prevalent codes per dimension, then build recurrence covariates.
  feats <- list()
  for (d in unique(cnt$dimension)) {
    dd <- cnt[dimension == d]
    prev <- dd[, .(np = uniqueN(person_id)), by = code][order(-np)][seq_len(min(.N, TOP_N_PREVALENT))]
    for (cd in prev$code) {
      sub <- dd[code == cd]; med <- median(sub$n)
      for (lab in c("once", "sporadic", "frequent")) {
        ids <- switch(lab,
          once     = sub[n >= 1, person_id],
          sporadic = sub[n >= max(med, 1), person_id],
          frequent = sub[n > med, person_id])
        feats[[paste(d, cd, lab, sep = "|")]] <- list(dim = d, ids = unique(ids))
      }
    }
  }

  ex <- cohort$exposed; ou <- cohort$outcome; pid <- cohort$person_id
  ranked <- rbindlist(lapply(names(feats), function(nm) {
    flag <- as.integer(pid %in% feats[[nm]]$ids)
    if (length(unique(flag)) < 2) return(NULL)
    data.table(dim = feats[[nm]]$dim, name = nm,
               score = bross_bias(flag, ex, ou))
  }))
  keep <- ranked[order(-score), head(.SD, TOP_K_PER_DIM), by = dim]$name

  X <- as.data.table(lapply(keep, function(nm) as.integer(pid %in% feats[[nm]]$ids)))
  setnames(X, keep)
  if (!is.null(forced)) X <- cbind(X, forced)   # age, sex, calendar year, indication, severity
  fit <- glm(ex ~ ., data = cbind(ex = ex, X), family = binomial())
  cohort[, ps := predict(fit, type = "response")]
  cohort[, .(person_id, exposed, outcome, ps)]
}