Simulating HDPS Under Sparse Confounding

Why Simulate HDPS Under Sparse Confounding?

High-dimensional propensity score (HDPS) methods are intended for settings with many measured candidate features but comparatively few important confounders. Their performance is not determined by the number of features alone. A rare confounder may be strongly associated with treatment but poorly estimated; a common feature may be weakly associated with treatment yet strongly prognostic for the outcome. Increasing the number of irrelevant features can also make selection and propensity estimation unstable. A simulation lets us vary these conditions separately while keeping the true treatment effect known.

The goal is not to construct a simulation that proves one strategy always wins. It is to specify a data-generating process (DGP), estimand, analysis strategies, and performance measures in advance, then see how results change across a planned grid of scenarios. Results are conditional on those choices: a simulation can reveal vulnerabilities in the designed settings, but it does not establish performance in every observational data set.

A Sparse-Confounding Data-Generating Process

One transparent starting point is to generate independent binary features \(X_j\sim\operatorname{Bernoulli}(q_j)\), then generate treatment according to

Choose only a few \(j\)'s to have both \(\beta_j\ne 0\) and a nonzero outcome coefficient \(\gamma_j\). Those are the true confounders in this particular DGP. Features with \(\beta_j=0\) but \(\gamma_j\ne 0\) are prognostic but do not affect treatment; features with \(\beta_j\ne 0\) but \(\gamma_j=0\) affect treatment but not the outcome conditional on the other features. The remaining features are null in both models.

For the main simulation, draw \(\varepsilon\sim N(0,\sigma^2)\) independently of \(X\) and of the random variable used to generate \(A\). Then the potential outcomes are independent of \(A\) conditional on \(X\), so full conditional exchangeability holds in this DGP. The true ATE is \(\tau\), since \(Y(1)-Y(0)=\tau\) for every unit. This explicit independence condition matters: merely imposing \(E[\varepsilon\mid A,X]=0\) gives a conditional-mean statement, not necessarily independence of the potential outcomes and treatment given \(X\). In this additive outcome model, the mean condition can support mean-based identification, but it should not be described as full exchangeability.

For every finite linear predictor, the logistic treatment probability is between zero and one. That mathematical positivity does not guarantee useful empirical overlap: some units can still have probabilities near zero or one, yielding extreme weights. Choose \(\alpha_0\) to target a prespecified marginal treatment proportion, such as one half, and record the resulting propensity distribution in each scenario. Recalibrate the intercept when changing other DGP parameters if holding the treatment proportion fixed is part of the design.

Specify the HDPS Strategies Before Running the Simulation

“HDPS” is not one uniquely defined ranking rule. To make a simulation reproducible, state exactly how features are screened and ranked. For example, require a minimum count of exposed and unexposed individuals with each binary feature present. Among eligible features, rank them by absolute treatment association, such as the absolute log odds ratio for \(A\) comparing \(X_j=1\) with \(X_j=0\). Separately rank them by the absolute standardized difference in outcome means between \(X_j=1\) and \(X_j=0\), measured among untreated individuals. Selecting the outcome association among controls avoids using the treatment effect itself as the feature’s outcome association in this example, though it remains an outcome-informed selection procedure.

Define the joint HDPS rule as retaining the \(k\) eligible features with the smallest sum of their two ranks, where rank 1 is strongest. Also define the treatment-only rule as retaining the \(k\) features with the strongest treatment association. These are operational definitions for a simulation, not claims that all implementations use these exact rules. Compare them with an oracle strategy that includes only the true confounders and a random-\(k\) feature strategy as a reference. Fit the same prespecified treatment model and use the same effect estimator for every strategy. This isolates differences due to feature choice more clearly than changing selection, model form, and estimator simultaneously.

Vary the factors in a planned grid: prevalence \(q_j\) of the true confounders, magnitudes and signs of their \(\beta_j\)'s, magnitudes of their \(\gamma_j\)'s, and total dimension \(p\). Also consider prognostic-only and treatment-only features, and vary the number of null features. State whether coefficients are held fixed or recalibrated as prevalence changes: the same log odds ratio can imply different marginal treatment imbalance when a feature is rare. Avoid changing several factors without recording them, since an apparent “dimension effect” could otherwise be caused by altered overlap or confounding strength.

Estimating Effects and Measuring Performance

For a fitted propensity score \(\hat e_i\), one common ATE estimator is the normalized inverse-probability weighted difference:

The weighting rationale follows from conditional exchangeability and a correctly specified propensity score. For example, \(E[AY/e(X)\mid X]=E[Y\mid A=1,X]\), because \(P(A=1\mid X)=e(X)\). Averaging over \(X\), exchangeability and consistency give \(E[AY/e(X)]=E[Y(1)]\). Similarly, \(E[(1-A)Y/\{1-e(X)\}]=E[Y(0)]\). The displayed normalized estimator is a sample analogue using weighted group means; it is not exactly unbiased in every finite sample. Its consistency also relies on adequate overlap and a suitable propensity estimation strategy. If selection omits important confounders or the fitted model is inadequate, weighting need not remove confounding bias.

Across \(B\) independent simulated data sets, record the estimate \(\hat\tau_b\) in each replicate. Useful summaries include empirical bias, variance, and root mean squared error (RMSE):

Report balance as well as effect performance. For each feature, a standardized mean difference (SMD) compares its treated and control means; after weighting, use the weighted means and a clearly stated scale, such as the unweighted pooled standard deviation in the original sample. Summarize balance among true confounders separately from balance across all measured features. A small average SMD can hide a severely imbalanced rare confounder, so report a maximum or upper percentile too. Balance is a diagnostic, not proof that unmeasured confounding is absent or that the causal assumptions hold.

Diagnostics, Limitations, and Common Mistakes

For each scenario, inspect treatment counts, feature counts in both treatment groups, fitted propensity distributions, maximum weights, and effective sample size, in addition to balance and effect summaries. Record the frequency of failed model fits or separation rather than silently dropping those replicates. If the treatment model is unstable, distinguish failure to estimate a propensity from poor causal performance among the replicates that succeeded. Report the number of successful replicates and assess whether excluding failures changes the conclusions. When comparing dimension, hold sample size fixed or deliberately vary it as another factor, and state which choice was made.

The number \(B\) controls Monte Carlo precision, not the sampling variability within a data set. For example, the Monte Carlo standard error of the estimated mean effect is approximately the replicate standard deviation divided by \(\sqrt{B}\). If that value is large relative to differences between strategies, increase \(B\) before interpreting a small ranking difference as meaningful. Report both Monte Carlo uncertainty and the replicate-to-replicate spread.

• Wrong: “The simulation generated an additive outcome with \(E[\varepsilon\mid A,X]=0\), so full exchangeability must hold.” Correction: A conditional-mean restriction alone does not imply equality of the full conditional distributions of potential outcomes. To claim full exchangeability in the stated DGP, generate \(\varepsilon\) independently of \((A,X)\), or separately justify the full conditional-independence condition.
• Wrong: “The maximum SMD is small, so the estimated treatment effect is unbiased.” Correction: Balance on measured features is a diagnostic, not a guarantee against unmeasured confounding, selection error, or model misspecification. Check true-confounder balance in the simulation and calculate effect bias directly against the known \(\tau\).
• Wrong: “More selected features must reduce bias.” Correction: Additional features may help control confounding or improve precision, but they may also worsen overlap, increase weight variability, or destabilize estimation. Compare bias and variance separately as well as through RMSE.
• Wrong: “A rare feature with a large treatment association should always be retained.” Correction: Its importance also depends on outcome association and the information available to estimate its effect. Report screening and fit failures, and vary prevalence rather than relying on one sparse-feature scenario.

The oracle is a useful benchmark, not a feasible observational strategy: it assumes the analyst knows which features truly confound. A treatment-only rule avoids outcome-informed ranking but can omit strongly prognostic features. Joint ranking can identify features relevant to both models, yet its performance depends on how associations are measured and on the selected-feature count. These comparisons assess specified algorithms under specified DGPs; they do not establish that any strategy is universally preferable.