High-Dimensional Propensity Score (hdPS)

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 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.

Decision diagram

flowchart TD
  U[Unmeasured confounder<br/>severity / frailty / access] -->|causes| E[Exposure: drug A vs B]
  U -->|causes| Y[Outcome]
  U -.->|leaves a footprint in| P[Pre-index code proxies<br/>dx / proc / rx / utilization]
  P -->|hdPS adjusts on these| Adj((blocks back-door<br/>U -> E ... U -> Y))
  I[Instrument<br/>formulary / plan marker] -->|predicts only| E
  I -.->|if wrongly selected| Amp[Bias amplification + variance inflation]
  C[Post-index code<br/>mediator / collider] -.->|if lookback leaks| Over[Collider / over-adjustment bias]

flowchart TD
  Start[Comparative-safety question in a large healthcare database] --> Q1{Rich pre-index codes<br/>as confounder proxies?}
  Q1 -->|No| Inv["Use investigator-specified PS / DRS"]
  Q1 -->|Yes| Q2{"Enough exposed and events?<br/>(>~50 events/arm)"}
  Q2 -->|No| Small["hdPS unstable -> parsimonious PS or exact matching<br/>Rassen 2011"]
  Q2 -->|Yes| Q3{Candidate window strictly pre-index?}
  Q3 -->|No| Fix["Re-extract: hard-stop before time zero<br/>avoid mediators / colliders"]
  Q3 -->|Yes| Q4{Screen out near-instruments?}
  Q4 -->|No| Drop[Remove exposure-only predictors<br/>avoid bias amplification]
  Q4 -->|Yes| Run["Run dimension-balanced hdPS<br/>force key confounders -> PS"]
  Run --> Use{Estimand?}
  Use -->|ATT| M[1:1 logit-PS matching]
  Use -->|ATE| W[IPTW]
  Use -->|ATO| O[Overlap weights]
  M --> Diag["Balance SMD<0.1 + negative controls + E-value"]
  W --> Diag
  O --> Diag

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

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).

Trade-offs

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)]
}

%let lookback = 365;
%let topk = 50;   /* covariates kept per dimension after Bross ranking */

/* Recurrence covariates from strictly pre-index codes: one binary per code (>=1 occurrence shown). */
proc sql;
  create table cand as
  select c.dimension, c.code, c.person_id, k.exposed, k.outcome
  from work.codes c
       inner join work.cohort k on c.person_id = k.person_id
  where c.code_date < k.index_date
    and c.code_date >= k.index_date - &lookback
  group by c.dimension, c.code, c.person_id, k.exposed, k.outcome;
quit;

/* Bross multiplicative bias per candidate: prevalence in exposed vs unexposed x covariate-outcome RR. */
proc sql;
  create table bross as
  select dimension, code,
         mean(exposed) as p_any,
         /* P(covariate) within exposure strata */
         sum(exposed)/ (select sum(exposed) from work.cohort)            as p_c1,
         (count(*)-sum(exposed))/
            ((select count(*)-sum(exposed) from work.cohort))            as p_c0,
         mean(outcome)                                                   as o_c1
  from cand
  group by dimension, code
  having calculated p_c1 > 0 and calculated p_c0 > 0;
quit;

data ranked;
  set bross;
  rr_cd = max(o_c1, 0.001) / 0.05;                 /* outcome RR vs cohort baseline risk (supply baseline) */
  bias  = (p_c1*(rr_cd-1)+1) / (p_c0*(rr_cd-1)+1);
  logbias = abs(log(max(bias, 1e-6)));
run;

/* Top-k per dimension -> selected covariate list. */
proc rank data=ranked out=rk descending ties=low groups=10000;
  by dimension; var logbias; ranks rnk;
run;
proc sql;
  create table selected as select dimension, code from rk where rnk < &topk;
quit;

/* Pivot selected codes to one indicator column per code, join investigator-forced covariates,
   then estimate the PS and 1:1 match (ATT). Macro-build the wide design from `selected` upstream. */
proc logistic data=work.analytic descending;     /* analytic = cohort + selected indicators + work.forced */
  class exposed (ref='0');
  model exposed = <selected_indicators> age sex cal_year indication severity;
  output out=ps_out p=ps;
run;

proc psmatch data=ps_out region=allobs;
  class exposed;
  psmodel exposed(treated='1');                   /* uses ps from the LOGISTIC step above */
  match method=greedy(k=1) distance=lps caliper=0.2;
  assess lps var=(age severity) / plots=(boxplot stddiff);  /* require post-match SMD < 0.1 */
  output out(obs=match)=matched matchid=mid;
run;

Citations

1. [1]Schneeweiss S, Rassen JA, Glynn RJ, Avorn J, Mogun H, Brookhart MA. High-dimensional propensity score adjustment in studies of treatment effects using health care claims data. Epidemiology. 2009;20(4):512-522.
2. [2]Schneeweiss S. Automated data-adaptive analytics for electronic healthcare data to study causal treatment effects. Variable selection for confounding adjustment in high-dimensional covariate spaces when analyzing healthcare databases. Epidemiology. 2017;28(2):237-248.
3. [3]Rassen JA, Blin P, Kloss S, Neugebauer RS, Platt RW, Pottegård A, Schneeweiss S, Toh S. High-dimensional propensity scores for empirical covariate selection in secondary database studies: planning, implementation, and reporting. Pharmacoepidemiol Drug Saf. 2023;32(2):93-106.

1. [4]Rassen JA, Glynn RJ, Brookhart MA, Schneeweiss S. Covariate selection in high-dimensional propensity score analyses of treatment effects in small samples. Am J Epidemiol. 2011;173(12):1404-1413.