Chapter 15: Outcome Regression and Propensity Scores
This chapter covers outcome regression and propensity score methods (stratification, standardization, and matching on the propensity score): the most commonly used parametric methods for causal inference. They work well for treatments fixed at a single time, but they are not designed for the complexities of time-varying treatments, which is why the book presents them only after the g-methods.
This chapter is based on Hernán and Robins (2020, chap. 15, pp. 199-208).
Key theme: outcome regression and propensity score methods need the same identifying conditions as every other method (exchangeability, positivity, and consistency / well-defined interventions), plus correctly specified models for the quantities they estimate along the way. Unlike the g-methods of Chapters 12-14, they have limited applicability to complex longitudinal data (Hernán and Robins 2020, 199).
1 15.1 Outcome Regression (pp. 199-201)
1.1 A Structural Model with Parameters for \(L\)
Structural nested models (Chapter 14) contain parameters for \(A\) and for \(A \times L\) product terms only, so g-estimation is agnostic about the \(L\)-\(Y\) relation. If instead we are willing to specify the \(L\)-\(Y\) relation within levels of \(A\), we can use a structural model that also has parameters for \(L\).
In Chapter 12 this model was called a faux marginal structural model: it has the form of a marginal structural model, but because it conditions on the entire vector \(L\) (not a subset \(V\)), the stabilized IP weights \(SW^A(L)\) all equal 1 and no weighting is needed (Hernán and Robins 2020, 199).
\(\beta_3\) is often called the “main effect” of \(L\), but it need not have a causal interpretation (the \(L\)-\(Y\) relation may itself be confounded); it only describes how the mean of \(Y^{a=0,c=0}\) varies with \(L\).
1.2 Conditional Effects from the Structural Model
Proof. Evaluate the model at \(a = 1\) and \(a = 0\) and subtract:
\[ \begin{aligned} &\operatorname{E}\mathopen{}\left[Y^{a=1,c=0} \mid L = l\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{a=0,c=0} \mid L = l\right]\mathclose{} \\ &\quad = (\beta_0 + \beta_1 \cdot 1 + 1 \cdot {\beta_2}^{\top} l + {\beta_3}^{\top} l) - (\beta_0 + \beta_1 \cdot 0 + 0 \cdot {\beta_2}^{\top} l + {\beta_3}^{\top} l) \\ &\quad = (\beta_0 + \beta_1 + {\beta_2}^{\top} l + {\beta_3}^{\top} l) - (\beta_0 + {\beta_3}^{\top} l) \\ &\quad = \beta_1 + {\beta_2}^{\top} l. \end{aligned} \]
So the stratum-specific effects depend only on \((\beta_1, \beta_2)\); the counterfactual means under no treatment depend on \((\beta_0, \beta_3)\).
1.3 Estimating the Structural Model by Outcome Regression
Proof. Part 1:
\[ \begin{aligned} \operatorname{E}\mathopen{}\left[Y^{a,c=0} \mid L = l\right]\mathclose{} &= \operatorname{E}\mathopen{}\left[Y^{a,c=0} \mid A = a, C = 0, L = l\right]\mathclose{} && \text{(conditional exchangeability; positivity makes the conditioning event non-null)} \\ &= \operatorname{E}\mathopen{}\left[Y \mid A = a, C = 0, L = l\right]\mathclose{} && \text{(consistency)}. \end{aligned} \]
Part 2: by part 1 and the structural model, \(\operatorname{E}\mathopen{}\left[Y \mid A = a, C = 0, L = l\right]\mathclose{} = \beta_0 + \beta_1 a + a {\beta_2}^{\top} l + {\beta_3}^{\top} l\) for every \((a, l)\) with positive probability among the uncensored, so the linear model with \(\alpha = \beta\) equals the conditional mean and is correctly specified; in vector form, \(\operatorname{E}\mathopen{}\left[Y \mid A, L, C = 0\right]\mathclose{} = {X}^{\top} \beta\). The moment conditions make every expectation in the normal equations finite. The population least-squares coefficients \(\alpha\) are the solutions of the normal equations \(\operatorname{E}\mathopen{}\left[X (Y - {X}^{\top} \alpha) \mid C = 0\right]\mathclose{} = 0\). The vector \(\beta\) solves them:
\[ \begin{aligned} \operatorname{E}\mathopen{}\left[X (Y - {X}^{\top} \beta) \mid C = 0\right]\mathclose{} &= \operatorname{E}\mathopen{}\left[\operatorname{E}\mathopen{}\left[X (Y - {X}^{\top} \beta) \mid A, L, C = 0\right]\mathclose{} \mid C = 0\right]\mathclose{} && \text{(iterated expectations)} \\ &= \operatorname{E}\mathopen{}\left[X (\operatorname{E}\mathopen{}\left[Y \mid A, L, C = 0\right]\mathclose{} - {X}^{\top} \beta) \mid C = 0\right]\mathclose{} && \text{(} X \text{ is a function of } (A, L)\text{)} \\ &= \operatorname{E}\mathopen{}\left[X ({X}^{\top} \beta - {X}^{\top} \beta) \mid C = 0\right]\mathclose{} && \text{(correct specification)} \\ &= 0. \end{aligned} \]
Because \(\operatorname{E}\mathopen{}\left[X {X}^{\top} \mid C = 0\right]\mathclose{}\) is invertible, the normal equations have a unique solution, so \(\alpha = \beta\).
So the structural parameters can be estimated by ordinary least squares, fitting the outcome regression model of Proposition 2.
Like stratification (Chapter 3), outcome regression adjusts for confounding by estimating the effect within strata of \(L\). In Section 13.2, the outcome regression was an intermediate step toward a standardized (marginal) mean; here the regression is the end of the procedure: we compare conditional means instead of standardizing them (Hernán and Robins 2020, 200).
The same approach works for discrete outcomes, e.g., a logistic model for \(\Pr[Y = 1 \mid A = a, C = 0, L]\) when \(Y\) is dichotomous.
Two methods for the same effect rely on different nuisance parameters (Definition 2):
- Outcome regression for \(\beta_1, \beta_2\): consistent (in general) only if \(\beta_0 + {\beta_3}^{\top} L\) correctly models how \(\operatorname{E}\mathopen{}\left[Y^{a=0,c=0} \mid L\right]\mathclose{}\) depends on \(L\). The nuisance parameters are \(\beta_0\) and \(\beta_3\).
- G-estimation of the structural nested model \(\operatorname{E}\mathopen{}\left[Y^{a,c=0} - Y^{a=0,c=0} \mid L\right]\mathclose{} = \beta_1 a + a {\beta_2}^{\top} L\): consistent (in general) only if the treatment model for \(\Pr[A = 1 \mid L]\) is correct. The nuisance parameters are those of the treatment model, e.g., \(\gamma_0, \gamma_1\) in \(\operatorname{logit}\Pr[A = 1 \mid L] = \gamma_0 + {\gamma_1}^{\top} L\).
Example: if \(L\) should enter as \(\beta_3 L + \beta_4 L^2\) but is modeled as \(\beta_3 L\) only, outcome regression is biased; g-estimation is unaffected by that error but is biased if the \(L\)-\(A\) relation is misspecified. For a time-fixed treatment, choosing between the two amounts to deciding which nuisance model we can specify more accurately. When possible, a better option is a doubly robust method (Fine Point 13.2) (Hernán and Robins 2020, 200).
2 15.2 Propensity Scores (pp. 201-202)
Quitters having a higher average estimated probability of quitting is exactly what we expect when \(L\) predicts quitting.
2.1 The Balancing Property
Proof. The book argues graphically (Technical Point 15.1); the same result follows from iterated expectations. Since \(\pi(L)\) is a function of \(L\), conditioning on \((L, \pi(L))\) is the same as conditioning on \(L\):
\[ \Pr[A = 1 \mid L, \pi(L)] = \Pr[A = 1 \mid L] = \pi(L). \]
Conditioning on \(\pi(L)\) alone, by the law of iterated expectations:
\[ \begin{aligned} \Pr[A = 1 \mid \pi(L)] &= \operatorname{E}\mathopen{}\left[\Pr[A = 1 \mid L] \mid \pi(L)\right]\mathclose{} \\ &= \operatorname{E}\mathopen{}\left[\pi(L) \mid \pi(L)\right]\mathclose{} \\ &= \pi(L). \end{aligned} \]
The two probabilities are equal, so \(A\) does not depend on \(L\) once \(\pi(L)\) is fixed: \(A \perp\!\!\!\perp L \mid \pi(L)\).
In words: individuals with the same \(\pi(L)\) may differ in \(L\) (e.g., smoking intensity and exercise), but among all individuals with a given value of \(\pi(L)\) in the super-population, the distribution of \(L\) is the same in the treated and the untreated. A function \(b(L)\) is a balancing score exactly when \(\pi(L)\) can be written as a function of \(b(L)\) (Rosenbaum and Rubin 1983, Theorem 2), so the propensity score is the coarsest balancing score.
If the distribution of \(\pi(L)\) were the same in the treated and the untreated, there would be no association between \(L\) and \(A\): then \(A \perp\!\!\!\perp\pi(L)\), and combining this with \(A \perp\!\!\!\perp L \mid \pi(L)\) from Proposition 3 gives \(A \perp\!\!\!\perp L\).
Two caveats from the book (Hernán and Robins 2020, 201–2):
- The propensity score balances only the measured covariates \(L\); unlike randomization, it does nothing about unmeasured confounders.
- In a finite study population the true propensity score balances \(L\) only approximately, because of sampling variability. An estimated propensity score from a correct model generally gives better balance; in a randomized experiment, adjusting for the estimated \(\pi(L)\) corrects both systematic and chance imbalances, whereas the true \(\pi(L)\) ignores the chance imbalances.
Rosenbaum and Rubin proved that exchangeability and positivity given \(L\) imply exchangeability and positivity given any balancing score (Rosenbaum and Rubin 1983, Theorem 3; Hernán and Robins 2020, 202). Graphically, \(\pi(L)\) sits between \(L\) and \(A\), with a deterministic arrow \(L \to \pi(L)\) (book Figure 15.2); conditioning on \(\pi(L)\) blocks the backdoor paths from \(A\) to \(Y\) through \(L\) (\(A \leftarrow \pi(L) \leftarrow L \to Y\)), so adjusting for it suffices.
Methods based on prognostic scores need stronger assumptions and do not extend readily to time-varying treatments; see Hansen (2008) and Abadie et al. (2013), as cited in Hernán and Robins (2020, 202).
Proof (Proof of part 1). \[ \begin{aligned} \Pr[A = 1 \mid Y^a, \pi(L)] &= \operatorname{E}\mathopen{}\left[\Pr[A = 1 \mid Y^a, L] \mid Y^a, \pi(L)\right]\mathclose{} && \text{(iterated expectations; } \pi(L) \text{ is a function of } L\text{)} \\ &= \operatorname{E}\mathopen{}\left[\Pr[A = 1 \mid L] \mid Y^a, \pi(L)\right]\mathclose{} && (Y^a \perp\!\!\!\perp A \mid L) \\ &= \operatorname{E}\mathopen{}\left[\pi(L) \mid Y^a, \pi(L)\right]\mathclose{} && \text{(definition of } \pi(L)\text{)} \\ &= \pi(L), \end{aligned} \]
which does not depend on \(Y^a\); hence \(Y^a \perp\!\!\!\perp A \mid \pi(L)\).
Proof (Proof of part 2). The proof of Proposition 3 showed \(\Pr[A = 1 \mid \pi(L)] = \pi(L)\), so \(\Pr[A = 1 \mid \pi(L) = s] = s\) and \(\Pr[A = 0 \mid \pi(L) = s] = 1 - s\) for every \(s\) with \(\Pr[\pi(L) = s] > 0\). Because \(L\) is discrete, \(\Pr[\pi(L) = s] = \sum_{l : \pi(l) = s} \Pr[L = l]\), so the values \(s\) with \(\Pr[\pi(L) = s] > 0\) are exactly the values \(\pi(l)\) with \(\Pr[L = l] > 0\). Positivity within levels of \(\pi(L)\) therefore holds if and only if \(0 < \pi(l) < 1\) for every \(l\) with \(\Pr[L = l] > 0\). Positivity within levels of \(L\) requires \(\Pr[A = 1 \mid L = l] = \pi(l)\) and \(\Pr[A = 0 \mid L = l] = 1 - \pi(l)\) to be positive for every \(l\) with \(\Pr[L = l] > 0\), which is the same condition.
Hence, if \(L\) is sufficient to adjust for confounding, so is \(\pi(L)\) (Hernán and Robins 2020, 202); part 2 of Theorem 1 states the positivity condition given \(\pi(L)\).
2.2 Using the Propensity Score
3 15.3 Propensity Stratification and Standardization (pp. 202-204)
Proof. For each \(a \in \{0, 1\}\),
\[ \begin{aligned} \operatorname{E}\mathopen{}\left[Y^{a,c=0} \mid \pi(L) = s\right]\mathclose{} &= \operatorname{E}\mathopen{}\left[Y^{a,c=0} \mid A = a, C = 0, \pi(L) = s\right]\mathclose{} && \text{(exchangeability given } \pi(L)\text{; positivity)} \\ &= \operatorname{E}\mathopen{}\left[Y \mid A = a, C = 0, \pi(L) = s\right]\mathclose{} && \text{(consistency)}. \end{aligned} \]
Subtracting the \(a = 0\) line from the \(a = 1\) line gives the result.
For treatment alone, without censoring, Theorem 1 shows that exchangeability and positivity given \(L\) carry over to \(\pi(L)\). When no one is censored, the joint \((A, C)\) conditions of Proposition 4 reduce to these treatment-only conditions, so they hold whenever exchangeability and positivity given \(L\) hold (with \(L\) discrete, as in Proposition 4). With censoring, \(\pi(L)\) alone need not suffice; one would instead condition on a score \(b(L)\) with \((A, C) \perp\!\!\!\perp L \mid b(L)\), for example the vector of joint probabilities \(\Pr[A = a, C = c \mid L]\) over all \((a, c)\).
In practice \(\pi(L)\) is generally continuous, so almost no two individuals share a value \(s\) (e.g., only one NHEFS individual had an estimated \(\pi(L)\) of 0.6563).
3.1 Stratifying on Deciles
With wide confidence intervals, differences between deciles may reflect chance.
3.2 A Problem with Categorized Propensity Scores
Within a decile, the distribution of the continuous \(\pi(L)\) may still differ between treated and untreated (e.g., higher average \(\pi(L)\) in the treated), so the groups may not be exchangeable within that decile.
Remedy: use the estimated \(\pi(L)\) as a continuous covariate in the outcome model \(\operatorname{E}\mathopen{}\left[Y \mid A, C = 0, \pi(L)\right]\mathclose{}\).
IP weighting and g-estimation never had this problem, because they used the numerical value of the estimated probability rather than a categorized version (Hernán and Robins 2020, 203).
For a dichotomous treatment the IP weight denominator is not \(\pi(L)\) itself but a function of it: \(\pi(L)\) for the treated (\(A = 1\)) and \(1 - \pi(L)\) for the untreated (\(A = 0\)).
The outcome model now requires a correct specification of the \(\pi(L)\)-\(Y\) relation. Because \(\pi(L)\) is one-dimensional, this is easy to protect with flexible terms such as cubic splines; IP weighting and g-estimation were agnostic about that relation. Note, though, that \(\pi(L)\) itself must still be estimated from a model of \(A\) on the high-dimensional \(L\), as in IP weighting and g-estimation.
3.3 Standardizing over the Propensity Score
The outcome model on \(\pi(L)\) estimates effects within levels \(s\) of \(\pi(L)\).
4 15.4 Propensity Matching (pp. 204-205)
Matching on \(\pi(L)\) is analogous to matching on a single continuous variable (Chapter 4).
- E.g., pair each treated individual with one (or more) untreated individuals with the same, or a close, propensity score.
- Under exchangeability and positivity given \(\pi(L)\), exact matching on the true \(\pi(L)\) makes the treated and the untreated in the matched population exchangeable, so, under consistency, association measures in the matched population equal the effect measures in the matched population.
- In practice we match on an estimate \(\hat\pi(L)\). Equal estimated scores need not mean equal true scores, so even with a correctly specified propensity model this exchangeability holds only approximately, with the approximation improving as the sample grows.
- Caliper matching adds a second approximation: some bias may remain, and it typically shrinks as the caliper narrows.
- With censoring, outcomes are observed only for the uncensored members of the matched population. One way to account for that is to weight the uncensored by the inverse of \(\Pr[C = 0 \mid L, A]\) (Chapter 12), which needs a correct censoring model, \(Y^{a,c=0} \perp\!\!\!\perp C \mid A, L\), and \(\Pr[C = 0 \mid L, A] > 0\).
The matched population can be given the \(\pi(L)\) distribution of the treated, of the untreated, or any other distribution (Hernán and Robins 2020, 204).
Variance estimation for matched estimators used to be an open problem; Abadie and Imbens (2006) solved it (Hernán and Robins 2020, 204).
- Matched and unmatched estimates may differ because of effect modification: when untreated far outnumber treated, matching keeps nearly all treated and drops many untreated, so the matched population resembles the treated and the estimate approaches the effect in the treated (Technical Point 4.1 gives IP weighting and standardization alternatives for that effect).
- Effect modification across propensity strata may suggest that decision makers treat those likely to benefit (Kurth et al. 2006, as cited in Hernán and Robins (2020, 206)), but statements like “beneficial for propensity scores between 0.11 and 0.93” have little policy relevance because they are not expressed in terms of measured variables \(L\).
- Other reasons matched and overall estimates may differ: positivity violations among the unmatched, or an unmeasured confounder that is more (or less) prevalent, or better (or worse) measured, in the matched population; apparent effect modification may reflect different residual confounding across propensity strata.
4.1 Defining Closeness
Exact matches on a continuous \(\pi(L)\) are rare, so we match on a close value, e.g., untreated individuals within \(\pm 0.05\) of the treated individual’s estimated \(\pi(L)\) (e.g., a treated individual with estimated \(\pi(L) = 0.6563\) could be matched to an untreated individual with 0.6579).
Bias-variance trade-off:
- criteria too loose: the \(\pi(L)\) distribution differs between matched treated and untreated, so exchangeability fails;
- criteria too tight: many individuals are excluded; approximate exchangeability, but wider confidence intervals.
4.2 Matching and Positivity
Matching does not distinguish random from structural nonpositivity.
The book also states this count as a percentage of the study population (0.01%); that percentage is not reproduced here because it is inconsistent with the counts: 2 individuals are about 0.1% of a study population of roughly 1,600.
4.3 Who Is in the Matched Population?
Restricting to, say, estimated \(\pi(L) < 0.67\) yields a population that is hard to describe: “individuals do not come with a propensity score tattooed on their forehead” (Hernán and Robins 2020, 205), so transportability to other populations is hard to judge.
4.4 Restrict on Real-World Variables Instead
Using propensity scores to find the region of overlap can be useful, but restricting the study population to that region is a “lazy” way to secure positivity: the automatic positivity of matching must be weighed against the difficulty of assessing transportability (Hernán and Robins 2020, 205). Even if everyone’s propensity score were visible, the population could remain ill-characterized, because the same propensity score value can mean different things in different settings.
5 15.5 Propensity Models, Structural Models, Predictive Models (pp. 205-208)
Part II used two kinds of models for causal inference.
Propensity models have been used for matching and stratification (this chapter), IP weighting (Chapter 12), and g-estimation (Chapter 14). Their parameters lack a causal interpretation because \(L\) and \(A\) can be associated for many reasons: the association reflects the effect of \(L\) on \(A\) under the book’s Figure 7.1, but not under Figures 7.2 and 7.3 (Hernán and Robins 2020, 205).
5.1 Two Classes of Structural Models
See Fine Point 14.1 of the book for the relation between structural nested models and faux semiparametric marginal structural models (Hernán and Robins 2020, 206).
5.2 Causal Versus Predictive Uses of Outcome Regression
5.3 Variable Selection: Prediction Tools Misapplied
Forward selection, backward elimination, stepwise selection, and machine learning are built to improve prediction. Applied to propensity or structural models, they can add superfluous or harmful covariates, causing
- inflated variances, and
- greater bias.
It is tempting to judge a propensity model by how accurately it predicts \(A\). It need not; it only needs to include the \(L\) that guarantee exchangeability.
Reporting predictive measures such as Mallows’s \(C_p\) for propensity models is common, but their relevance for causal inference is questionable (Hernán and Robins 2020, 207). Covariates strongly associated with \(A\) but not needed for exchangeability do not reduce bias; including them can yield very large variances or even amplify existing bias.
In the limit of perfect prediction, every treated individual has \(\pi(L) = 1\) and every untreated individual has \(\pi(L) = 0\): there is no overlap and the analysis is impossible (Hernán and Robins 2020, 207).
5.4 Self-Inflicted Bias and Model Specification
- Including colliders can cause systematic bias, even though colliders may be good predictors.
- Including instruments (Chapter 16) can amplify bias from unmeasured confounders.
- Chapter 18 returns to these issues.
6 Summary
| Method | Model used | Additional modeling requirement |
|---|---|---|
| Outcome regression (faux MSM) | \(\operatorname{E}\mathopen{}\left[Y \mid A, C = 0, L\right]\mathclose{}\) | Correct outcome model (nuisance parameters \(\beta_0, \beta_3\); Definition 2) |
| G-estimation (SNM) | \(\Pr[A = 1 \mid L]\); with censoring, also \(\Pr[C = 0 \mid L, A]\) for the censoring weights | Correct treatment model and correct structural nested model; with censoring, also a correct censoring model |
| Propensity stratification / outcome regression on \(\pi(L)\) | \(\pi(L)\) and \(\operatorname{E}\mathopen{}\left[Y \mid A, C = 0, \pi(L)\right]\mathclose{}\); with censoring, the joint score \(b(L)\) in place of \(\pi(L)\) | Correct propensity model and correct score-\(Y\) relation; with censoring, also a correct censoring model |
| Propensity standardization | as above | as above; targets the population effect |
| Propensity matching | \(\pi(L)\); with censoring, also \(\Pr[C = 0 \mid L, A]\) | Correct propensity model; with censoring, also a correct censoring model; estimate refers to the matched population |
Takeaways (Hernán and Robins 2020, 199–208):
- Outcome regression and propensity score methods need the same identifying conditions as the g-methods, plus correct models for their nuisance parameters (Definition 2).
- The propensity score balances measured covariates only.
- Matching guarantees positivity automatically but can produce a target population that is hard to characterize; restricting on real-world variables is preferable.
- Propensity models should include the variables needed for exchangeability, not the best predictors of treatment.
- Unlike the g-methods, these methods do not extend readily to time-varying treatments (Part III).