Propensity Scores When Treatment Isn’t Binary
Background
You know how to work with propensity scores when treatment is binary. Treatment is on or off, you model the probability of being treated given covariates, and then you match, weight, or stratify on that single number to recover a treated-versus-control contrast. I have written about the zoo of methods that build on it, and it is one of the most reliable workhorses in the applied causal toolkit.
But many of the treatments we actually care about are not switches, but dials. Hours of job training, milligrams of a drug, dollars of ad spend, days of hospital stay. The reflex is to dichotomize: call the top half “high dose,” the bottom half “low dose,” and run the familiar binary machinery. This is a mistake I have made myself, and it quietly throws away the most interesting part of the problem. The question with a dial is not “does treatment help?” but what does the whole response curve look like? — does it rise and then plateau, is there a sweet spot, are there diminishing returns? Collapsing the dose into two buckets answers none of that.
So how do we think about propensity scores when treatment goes from a switch to a dial? The reassuring answer is that the logic of the propensity score survives almost entirely intact. Let’s take a closer look.
Notation
Let \(T \in \mathcal{T}\) be the treatment, now taking values in an interval \([t_0, t_1]\) rather than \(\{0, 1\}\). For each unit we posit a unit-level dose–response function \(Y_i(t)\), the potential outcome that would be observed at dose \(t\). The target is the average dose–response function
\[ \mu(t) = \mathbb{E}[Y_i(t)], \qquad t \in \mathcal{T}, \]
an entire curve rather than a single number like \(\mathbb{E}[Y(1)] - \mathbb{E}[Y(0)]\). We observe covariates \(X_i\), the realized dose \(T_i\), and \(Y_i = Y_i(T_i)\).
The identifying assumption generalizes unconfoundedness in a subtle way. Hirano and Imbens call it weak unconfoundedness:
\[ Y(t) \perp T \mid X \quad \text{for all } t \in \mathcal{T}. \]
The word “weak” matters. We do not require the whole family \(\{Y(t)\}_{t \in [t_0,t_1]}\) to be jointly independent of \(T\) given \(X\) — only that the independence holds separately at each dose level. That is all the machinery needs, and it is a genuinely weaker requirement.
A Closer Look
The scalar score becomes a density
With binary treatment, the propensity score \[e(X) = P(T = 1 \mid X)\] is a single number in \([0,1]\). Its magic is dimension reduction: Rosenbaum and Rubin showed that conditioning on that one scalar is enough to remove covariate bias, no matter how many covariates went into it.
When \(T\) is continuous, “the probability of receiving treatment” is no longer a probability — it is a density. The natural object is the conditional density of the treatment given covariates,
\[ r(t, x) = f_{T \mid X}(t \mid x), \]
and the generalized propensity score (GPS) is this density evaluated at the dose a unit actually received:
\[ R = r(T, X). \]
This is the one change that ripples through everything. The score is still a scalar random variable — but it is a scalar you get by plugging into a density, not by reading off a probability, and it now depends on which dose you are asking about.
What the GPS balances
The binary score’s defining property is balance: within strata of equal \(e(X)\), treatment is independent of \(X\). The GPS inherits a version of this, but with a twist worth internalizing:
\[ X \perp \mathbf{1}\{T = t\} \mid r(t, X). \]
Read carefully: balance holds at each dose level \(t\) separately, conditioning on the score evaluated at that same \(t\).You use as many propensity scores as there are levels of the treatment, yet you never use more than one at a time. This is the mental model shift. There is no single number that balances covariates across all doses simultaneously; there is a family of scores, one per dose, each balancing locally. Notably, this balancing property is mechanical — it falls out of the definition and does not itself require unconfoundedness. Combine it with weak unconfoundedness and you get the payoff: assignment is unconfounded given the GPS.
Recovering the curve
The estimation strategy is where the GPS earns its keep, and it is pleasingly concrete.
- Model the treatment. Estimate the conditional density \(r(t, x)\) — e.g. fit by regression, and form the fitted score \(\hat R_i = \hat r(T_i, X_i)\) for each unit.
- Regress outcome on dose and score. Estimate \(\beta(t, r) = \mathbb{E}[Y \mid T = t, R = r]\), treating the two scalars \(T\) and \(\hat R\) as the only regressors (say, a flexible polynomial in both).
- Average at each dose. For a target dose \(t\), average the fitted surface over the score evaluated at that dose: \(\hat\mu(t) = \frac{1}{N}\sum_i \hat\beta\big(t, \hat r(t, X_i)\big)\). Sweep \(t\) across its range to trace out the whole curve.
The step that trips people up is the third one. You do not average over \(R = r(T, X)\), the score at each unit’s own realized dose. You average over \(r(t, X)\), the score evaluated at the target dose \(t\) for everyone. Getting this wrong silently biases the curve. Note also that \(\beta(t, r)\) has no causal meaning on its own — the regression coefficients are not effects. Only after the averaging step does \(\hat\mu(t)\) become interpretable, which is a familiar theme from interpreting regression adjustments causally.
Multiple implementations
The GPS defined above is Hirano and Imbens’ route, and it is the cleanest place to build intuition. But the same core idea (reducing the covariates to a low-dimensional score summarizing treatment assignment) shows up in a few different shapes.
Imai and van Dyk arrive at the general case from a different angle. Rather than a density, they work with the propensity function: the entire conditional distribution of \(T \mid X\). When that distribution belongs to a parametric family, the propensity function collapses to its parameter (e.g. the mean of a normal), and you subclassify or adjust on that parameter. This unifies the ordinal, categorical, and continuous cases under one theory and connects back to Imbens’ original multi-valued treatment argument.
Then there is the weighting route. Instead of adjusting for the score in a regression, you can weight each observation by the inverse of its GPS — the continuous analog of inverse-probability weighting, which underpins marginal structural models. This is tempting because it mirrors the weighting logic so central to causal inference, but I would sound a caution: with a continuous treatment the weights are ratios of densities, and they blow up wherever the estimated density is small. Extreme weights are already the Achilles’ heel of ordinary IPW; with a continuum they can be worse, and stabilized weights become close to mandatory. A third route, stratification on the GPS, has been developed more recently as a way to sidestep the instability of weighting while keeping the nonparametric spirit.
The deeper pitfall cuts across all three implementations: overlap. In the binary world, positivity means every unit has some chance of being treated and untreated. With a dial, you need positive density at every dose for every covariate profile — a far more demanding requirement, and one that is easy to violate without noticing. And because the score rides on a fitted density model, misspecifying \(T \mid X\) contaminates everything downstream. Check balance dose-by-dose; do not trust a single global diagnostic.
Bottom Line
- The logic of the propensity score survives the jump from binary to continuous treatment; what changes is that the score becomes a conditional density evaluated at the realized dose.
- The estimand is a whole dose–response curve \(\mu(t)\), not a single contrast — which is exactly what dichotomizing a dose throws away.
- Balance holds one dose at a time: there are as many scores as dose levels, but you only ever use one at once.
- Estimation is a clean two-step recipe, but you must average the fitted surface over the score at the target dose, not each unit’s own dose.
- Overlap is the real bottleneck — positivity at every dose is demanding — and inverse-GPS weighting is fragile, so treat extreme weights with suspicion.
Where to Learn More
The foundational trio is short and readable. Imbens (2000) introduces the generalized propensity score for multi-valued treatments in a five-page Biometrika note. Hirano and Imbens (2004) extend it to the continuous case and lay out the two-step estimator and balancing theory — the single best starting point. Imai and van Dyk (2004) give the complementary “propensity function” formulation that unifies the whole family. For an applied illustration, Kluve et al. (2012) estimate the dose–response of training duration on employment, and Brown et al. (2021) develop the stratification variant.
References
Brown, D. W., Greene, T. J., Swartz, M. D., Wilkinson, A. V., & DeSantis, S. M. (2021). Propensity Score Stratification Methods for Continuous Treatments. Statistics in Medicine, 40(5), 1189–1203.
Hirano, K., & Imbens, G. W. (2004). The Propensity Score with Continuous Treatments. In A. Gelman & X.-L. Meng (Eds.), Applied Bayesian Modeling and Causal Inference from Incomplete-Data Perspectives (pp. 73–84). Wiley.
Imai, K., & van Dyk, D. A. (2004). Causal Inference with General Treatment Regimes: Generalizing the Propensity Score. Journal of the American Statistical Association, 99(467), 854–866.
Imbens, G. W. (2000). The Role of the Propensity Score in Estimating Dose-Response Functions. Biometrika, 87(3), 706–710.
Kluve, J., Schneider, H., Uhlendorff, A., & Zhao, Z. (2012). Evaluating Continuous Training Programmes by Using the Generalized Propensity Score. Journal of the Royal Statistical Society: Series A, 175(2), 587–617.
Robins, J. M., Hernán, M. Á., & Brumback, B. (2000). Marginal Structural Models and Causal Inference in Epidemiology. Epidemiology, 11(5), 550–560.
Rosenbaum, P. R., & Rubin, D. B. (1983). The Central Role of the Propensity Score in Observational Studies for Causal Effects. Biometrika, 70(1), 41–55.