Semi-supervised learning and the question of true versus estimated propensity scores

Throughout the paper we use the following conventions:

• •

  Random variable: italicized, uppercase letter (e.g. $X$, $A$, …)

• •

  Scalar realization of random variable: italicized, lowercase letter (e.g. $x$, $a$, …)

• •

  Vector realization of random variable: lowercase letter (e.g. $\mathrm{x}$, $\mathrm{a}$, …)

• •

  Matrix realization of random variable: bold, uppercase letter (e.g. $\mathbf{X}$, $\mathbf{A}$, …)

• •

  Metric space of random variable’s support: math calligraphy letter (e.g. $\mathcal{X}$, $\mathcal{A}$, …)

We first draw a distinction between observational and experimental studies for treatment effect estimation. In experimental settings, the treatment of interest is deliberately randomized according to a pre-specified design, typically with the goal of enabling straightforward identification of the treatment effect. In observational settings, data on treatment, outcome and covariates are observed without any direct manipulation. Researchers wishing to infer treatment effects with observational data must rely on a series of assumptions, which we now introduce.

Let $Y$ refer to the outcome of interest, $Z$ denote a binary treatment assignment, and $X$ be a vector of covariates drawn from covariate space $\mathcal{X}$. We use the potential outcomes notation $Y^{1}$ and $Y^{0}$ to refer to the counterfactual values of $Y$ when $Z=1$ and $Z=0$ (Hernan and Robins (2020)). Our estimand of interest, is the average treatment effect: $\mbox{E}[Y^{1}-Y^{0}]$.

Since the potential outcomes $(Y^{1},Y^{0})$ are never observed simultaneously, we rely on a number of identifying assumptions common to the causal inference literature (Rubin (1980), Rosenbaum and Rubin (1983), Hernan and Robins (2020)):

1. 1.

   Stable unit treatment value assumption (SUTVA): for any sample of size $n$ with $Y\in\mathcal{Y}$ and $Z\in\mathcal{Z}$, $(Y_{i}^{1},Y_{i}^{0})\perp Z_{j}$ for all $i,j\in\{1,...,n\}$ with $j\neq i$

2. 2.

   Positivity: $0<\mbox{P}(Z=1\mid X=\mathbf{X})<1$ for all $\mathbf{X}\in\mathcal{X}$

3. 3.

   Conditional exchangeability: $(Y^{1},Y^{0})\perp Z\mid X$

Taken together, these assumptions enable identification of average treatment effects using observational data after adjusting for the effect of $X$ on $Z$ and $Y$. Thus, we can rewrite our estimand as $\textrm{ATE}=\mbox{E}[Y^{1}-Y^{0}]=\mbox{E}_{X}[\mbox{E}[Y\mid X,Z=1]-\mbox{E}[Y\mid X,Z=0]]$. There are many ways to adjust for $X$ in estimating the ATE, but this paper focuses on methods that use the propensity score. Rosenbaum and Rubin (1983) define the propensity score as $p(X)=\mbox{P}(Z=1\mid X)$ and show that $p(X)$ has several desirable properties:

1. 1.

   $X\perp Z\mid p(X)$

2. 2.

   $(Y^{1},Y^{0})\perp Z\mid p(X)$

In short, the propensity score can be used to ensure covariate balance across treatment groups and the resulting covariate balance deconfounds the treatment effect estimate of $Z$ on $Y$. The deconfounding property of $p(X)$ implies that the treatment effect estimand can be rewritten as $\textrm{ATE}=\mbox{E}[Y^{1}-Y^{0}]=\mbox{E}_{p(X)}[\mbox{E}[Y\mid p(X),Z=1]-\mbox{E}[Y\mid p(X),Z=0]]$.

Above, we let $X$ refer to the set of covariates of $Y$ and $Z$, but for the purposes of this discussion, we introduce a taxonomy of covariate types. First, observe that, using a similar framing as Hahn et al. (2020), the relationship between $Y$, $Z$ and $X$ can be expressed as:

where $X_{\mu}$ refers to the set of prognostic variables, $X_{\tau}$ refers to the set of moderator variables, and $X_{\pi}$ refers to the set of propensity variables. We consider that any variable in $X$ can be classified as follows:

1. 1.

   Confounder: Any variable which appears in $X_{\pi}$ and either of $X_{\mu}$ or $X_{\tau}$

2. 2.

   Pure prognostic / moderator variable: Any variable which appears in either $X_{\mu}$, $X_{\tau}$, or both, but not in $X_{\pi}$

3. 3.

   Instrument: Any variable which appears in $X_{\pi}$ and not in $X_{\mu}$ nor $X_{\tau}$ (this is a “pure propensity” variable)

4. 4.

   Extraneous variable: Any variable in $X$ which appears in none of $X_{\tau}$, $X_{\mu}$, or $X_{\pi}$)

Hirano, Imbens and Ridder (2003) discuss the inverse-propensity weighted estimator of average treatment effect, which estimates the ATE as follows

They show that even when the true propensities are known, using them directly in the IPW estimator is asymptotically inefficient. We consider the finite sample variance properties of this estimator on a simple data generating process, which we define below:

• •

  Covariates: $X_{1},X_{2}\sim\mbox{Bernoulli}(1/2)$, letting $X=(X_{1},X_{2})$

• •

  Propensity function: $p(X)=1/\left(1+\exp\{-(\alpha_{1}+\alpha_{2}X_{1})\}\right)$

• •

  Treatment: $Z\sim\mbox{Bernoulli}(p(X))$

• •

  Outcome: $Y\sim\mathcal{N}(\gamma_{1}X_{1}+\gamma_{2}X_{2}+\tau Z,1)$

We let the random variable $N$ denote the overall sample size and define subset-specific sample sizes as follows:

• •

  $N_{x}$: the number of observations with $X=x$

• •

  $N_{x,z}$: the number of observations with $X=x$ and $Z=z$

First, observe in this case that $p(X)$ takes on two distinct values: $1/\left(1+\exp\{-\alpha_{1}\}\right)$ if $X_{1}=0$, and $1/\left(1+\exp\{-\alpha_{1}-\alpha_{2}\}\right)$ if $X_{1}=1$. We can split in the sum into eight parts, stratifying on the unique values of ($X_{1}$, $X_{2}$, $Z$). For a given stratum, the weights in the denominator, as well as the $Z_{i}$ and $1-Z_{i}$ indicators, are all the same and can be factored out. We refer to “true propensities” estimator as $\mbox{IPW}_{p(X)}$ to indicate that the true propensity function, $p(X)$ is the basis for the weights in the denominator. Let $p_{x}=\mbox{P}(Z=1\mid X=x)=\mbox{P}(Z=1\mid X_{1}=x_{1})$. The portion of the “true propensities” estimator represented by a given $\{X=x\}$ (i.e. $\{X_{1}=x_{1},X_{2}=x_{2}\}$) stratum is given by

Where $\mbox{IPW}_{p(X)}=\mbox{IPW}_{p(X),X=(0,0)}+\mbox{IPW}_{p(X),X=(1,0)}+\mbox{IPW}_{p(X),X=(0,1)}+\mbox{IPW}_{p(X),X=(1,1)}$ and $\hat{p}_{x}=\left(N_{x,\mathrm{z}=1}/N_{x}\right)$. We see that replacing the true propensity weights with these empirical weights results in the following estimator:

We can compare the variances of these two estimators, by first observing that the data are i.i.d., so that conditional on the number of observations with $Z=1$ and $X=x$, the covariances of the strata means are 0. We derive the variances of the two estimators in Appendix A, but the key to understanding why $\mbox{V}\left[\mbox{IPW}_{\hat{p}(X),x}\right]\leq\mbox{V}\left[\mbox{IPW}_{p(X),x}\right]$ is to notice that the ratios $\left(\hat{p}_{x_{1}}/p_{x_{1}}\right)$ and $\left(1-\hat{p}_{x}\right)/\left(1-p_{x}\right)$ are random variables which converge to $1$ as $n\rightarrow\infty$ but are noisy in finite samples. These finite sample imbalances are an ancillary statistic for $\tau$, and the $\hat{p}(X)$ estimator conditions on this ancillary statistic for more efficient inference. For a detailed discussion of the role of ancillary statistics in causal inference, readers are referred to Chapter 10 of Hernan and Robins (2020) and the references therein.

This result seems to imply a paradox. Even if a propensity function is known exactly, is the analyst better off ignoring that information and estimating the propensity scores? While it is true that directly weighting by the true propensities yields higher variances, we show that the true propensity scores can still be used to achieve efficiency. Consider in this case that $p(X)=1/\left(1+\exp\{-\alpha_{1}-\alpha_{2}X_{1}\}\right)$, so $Z$ can be modeled via logistic regression on $X_{1}$, $X_{2}$, and an intercept term. Fitting this model on a sample $\mathbf{X}=(\mathrm{x}_{1},\mathrm{x}_{2})$ and $\mathrm{z}$ would yield coefficient estimates $\hat{\beta}_{1}$, $\hat{\beta}_{2},$ and $\hat{\beta}_{3}$ such that $\hat{p}(\mathbf{X})=\frac{1}{1+\exp\{-\hat{\beta}_{1}-\hat{\beta}_{2}\mathrm{x}_{1}-\hat{\beta}_{3}\mathrm{x}_{2}\}}$, where the true values of the parameters are $\beta_{1}=\alpha_{1}$, $\beta_{2}=\alpha_{2}$, and $\beta_{3}=0$. Thus, if we know the true value of $p(\mathbf{X})$, we can calculate the true $\beta_{1}+\beta_{2}\mathrm{x}_{1}+\beta_{3}\mathrm{x}_{2}$ by logit transformation

Since both $\beta_{1}+\beta_{2}\mathrm{x}_{1}+\beta_{3}\mathrm{x}_{2}$ and $\hat{\beta}_{1}+\hat{\beta}_{2}\mathrm{x}_{1}+\hat{\beta}_{3}\mathrm{x}_{2}$ are one-dimensional linear functions of $\mathbf{X}$, we could express $\hat{\beta}_{1}-\hat{\beta}_{2}\mathrm{x}_{1}-\hat{\beta}_{3}\mathrm{x}_{2}=\hat{\eta}_{0}-\log\left(\frac{1}{p(\mathbf{X})}-1\right)\hat{\eta}_{1}$ for some $\hat{\eta}_{0},\hat{\eta}_{1}\in\mathbb{R}$ and observe that this is equivalent to a logistic regression of $Z$ on a 1-dimensional vector $-\log\left(\frac{1}{p(\mathbf{X})}-1\right)$ with an intercept term.

This 1-dimensional regression adjustment is similar to the sample correction discussed in the context of the $\hat{p}(X)$ estimator, with one key difference: the differences between true and empirical treatment probabilities are adjusted across the support of $p(X)$ rather than $X$. In our two-covariate example, the difference is that $p(X)$ is only a function of $X_{1}$, so has two unique values, rather than the four distinct values of $X$.

We introduce the following shorthand to refer to the three estimators:

1. 1.

   $p(X)$: IPW estimator using the true propensity function, $p(X)$ as weights

2. 2.

   $\hat{p}(X)$: IPW estimator using the covariate-estimated propensity function, $\hat{p}(X)$, as weights

3. 3.

   $\hat{p}(p(X))$: IPW estimator using the true propensity function, $p(X)$, as a 1-dimensional variable for estimating sample propensity weights

Hirano, Imbens and Ridder (2003) and the example above demonstrate that the $\hat{p}(X)$ estimator is lower variance than the $p(X)$ estimator. However, we see in Appendix B that the $\hat{p}(p(X))$ estimator has lower variance than the $\hat{p}(X)$ estimator as long as the following two conditions are satisfied:

• •

  $X_{2}$ is not a pure prognostic / moderator variable: $\mbox{E}\left[\bar{Y}_{X_{1}=x_{1},X_{2}=1,Z=z}\right]=\mbox{E}\left[\bar{Y}_{X_{1}=x_{1},X_{2}=0,Z=z}\right]$ for all $x_{1},z$. In our simple data generating process, this is true if $\gamma_{2}=0$.

• •

  The outcome variance does not differ along the support of $X_{2}$: $\mbox{V}\left[\bar{Y}_{X_{1}=x_{1},X_{2}=1,Z=z}\right]=\mbox{V}\left[\bar{Y}_{X_{1}=x_{1},X_{2}=0,Z=z}\right]$ for all $x_{1},z$. In our simple data generating process, this is true as the variance of $Y$ is constant with respect to $X$.

The most interesting takeaway from this result is that the cases in which the $\hat{p}(X)$ estimator outperforms the $\hat{p}(p(X))$ estimator do not pertain to the intended use of the propensity score, to deconfound treatment assignment and enable identification of the ATE. The $\hat{p}(X)$ estimator can attain a lower variance, but it does so by stratifying on a pure prognostic / moderator variable, $X_{2}$. This points to an unintended use of weighting / stratification estimators, rather than a fundamental problem with using true propensities. Adjusting for pure prognostic / moderator variables can be achieved by direct regression adjustment, and Hahn et al. (2020) provide a nonparametric Bayesian approach that shows promising results on a number of difficult data generating processes.

To demonstrate the difference in variances of the three estimators in finite samples, we conduct a simulation study. We use the same data generating process discussed in the finite sample variance calculations above, in which $Y$ is normal, $X$ and $Z$ are Bernoulli, and $X$ has $p=2$ dimensions. In the simulations presented below, we generate 10,000 datasets for each value of $N$ and estimate the average treatment effect using the IPW estimator, comparing three possible weighting options:

1. a)

   $p(X)$: True propensities ($\{0.3,0.7\}$)

2. b)

   $\hat{p}(X)$: Estimated propensities (logistic regression of $Z$ on $X$)

3. c)

   $\hat{p}(p(X))$: Sample adjustment using true propensities (logistic regression of $Z$ on $\mbox{logit}\left(p(X)\right)$)

Our simulations confirm the points made in Hirano, Imbens and Ridder (2003) and Rosenbaum (1987) that using the true propensities directly leads to higher variances in finite samples. However, we also see that the true propensities can be used in a one-dimensional model of $Z$ to attain even lower variances than estimators using the fitted propensities.

In practice, knowledge of the exact true propensity function is unlikely outside of randomized experiments. However, it is certainly possible to imagine investigators having access to auxiliary data that enable more accurate estimates of the true propensities. The received wisdom of Hirano, Imbens and Ridder (2003), Lunceford and Davidian (2004), and Rosenbaum (1987) would seem to imply that such accuracy gains impede efficient estimation of the ATE. The primary challenge in using propensity scores estimated from a larger auxiliary dataset stems from the utility of conditioning on finite sample treatment probabilities in weighting estimators. However, we have shown in Section 2.3 that the same variance reduction can be achieved by conditioning on finite sample probabilities formed by $p(X)$ rather than $X$, and this result extends to an estimate of $p(X)$ derived from arbitrarily large datasets.

We can formalize the notion of using a large auxiliary dataset in estimating the ATE by framing it as a semi-supervised learning problem. Semi-supervised learning refers to a broad class of techniques that combine labeled and unlabeled data in statistical learning problems (Zhu and Goldberg (2009), Belkin, Niyogi and Sindhwani (2006), Liang, Mukherjee and West (2007)). Consider a dataset with $N$ observations. As above, we let $Y$, $Z$, and $X$ denote outcome, treatment, and covariates, respectively. We refer to the data points with observed outcomes as “labeled data,” and let $N_{o}$ refer to the number of such labeled observations and $N_{m}=N-N_{o}$ refer to the number of data points with missing outcomes. We index an outcome’s status as missing or observed by $M$.

Following the notation of Liang, Mukherjee and West (2007), we denote “unlabeled,” or marginal, data as

• •

  $X^{m}=\{X_{i};i=N_{o}+1,...,N_{o}+N_{m}\}$, and

• •

  $Z^{m}=\{Z_{i};i=N_{o}+1,...,N_{o}+N_{m}\}$

Similarly, we denote “labeled” data as

• •

  $Y^{o}=\{Y_{i};i=1,...,N_{o}\}$,

• •

  $X^{o}=\{X_{i};i=1,...,N_{o}\}$, and

• •

  $Z^{o}=\{Z_{i};i=1,...,N_{o}\}$,

We observe that $M=1$ where $X,Z\in(X^{m},Z^{m})$ and $M=0$ elsewhere. The discussion in Section 2.3 shows that the set of marginal $(X^{m},Z^{m})$ data can be put to good use in estimating the average treatment effect as follows:

1. 1.

   Obtain an estimate, $p_{e}(X)$ of the true propensities using the labeled and unlabeled data $(X,Z)$

2. 2.

   Adjust the estimated propensities to the labeled sample by modeling $Z^{o}\sim\mbox{Bernoulli}\left(p_{e}(X^{o})\right)$, and denote these adjusted estimates by $p^{*}(X^{o})$

3. 3.

   Estimate the ATE using only the labeled data ($Y^{o}$, $Z^{o}$, and $p^{*}(X^{o})$)

As with other constructions in the semi-supervised learning literature, this framing allows us to profitably use information about the statistical patterns in the unlabeled data (in this case, $p(Z=1\mid X)$) to better estimate an effect among the labeled data ($\mbox{E}[Y^{1}-Y^{0}]$).

We employ three estimators which use propensity scores in our simulation study:

• •

  Inverse Propensity Weighted (IPW) estimator (Hirano, Imbens and Ridder (2003)), introduced in Section 2.3 and defined as

  $$\textrm{ATE}_{\textrm{IPW}}=\frac{1}{N}\sum_{i=1}^{N}\left(\frac{Y_{i}Z_{i}}{p(X_{i})}-\frac{Y_{i}(1-Z_{i})}{(1-p(X_{i}))}\right)$$

• •

  Targeted Maximum Likelihood Estimator (TMLE) (van der Laan (2010a), van der Laan (2010b), Gruber and van der Laan (2009))

  $$\textrm{ATE}_{\textrm{TMLE}}=\frac{1}{N}\sum_{i=1}^{N}\left(Q^{*}\left(Z_{i}=1,X_{i}\right)-Q^{*}\left(Z_{i}=0,X_{i}\right)\right)$$

  where $Q^{*}\left(Z_{i},X_{i}\right)$ represents a semi-parametric model of the outcome $Y$ which incorporates the propensity score.

• •

  Bayesian Causal Forests (BCF) (Hahn et al. (2020))

  $$\textrm{ATE}_{\textrm{BCF}}=\frac{1}{N}\sum_{i=1}^{N}\left(\bar{f}(X_{i},Z_{i}=1)-\bar{f}(X_{i},Z_{i}=0)\right)$$

  where $\bar{f}(X_{i},Z_{i})$ is an average of posterior simulations (across all simulations) of two BART models, defined as $f(X_{i},Z_{i})=\mu(X_{i},p(X_{i}))+\tau(X_{i})Z_{i}$.

Given the two step nature of each of the above estimators, in which a propensity model is first constructed and its estimates then used to calculate an average treatment effect, we consider three approaches to treatment effect estimation:

• •

  Complete case analysis: We discard all unlabeled data samples, estimate $\hat{p}(\mathbf{X}^{o})$, and then compute average treatment effects on only the labeled observations.

• •

  Semi-supervised: We use the labeled and unlabeled $\mathbf{X}$ and $\mathrm{z}$ values to estimate $\hat{p}(\mathbf{X})$ and then use the $\hat{p}(\mathbf{X}^{o})$ predictions to compute average treatment effects on the labeled observations.

• •

  True propensities: Simulation studies provide us with actual probabilities of receiving treatment, so we can use these true $p(\mathbf{X}^{o})$ to compute average treatment effects on the labeled observations. Of course, in most real world applications, the exact true propensities will not be available, and we consider this scenario as a limiting case in which an arbitrarily large amount of unlabeled data were available for modeling $\hat{p}(\mathbf{X})$.

We construct confidence intervals for our simulations as follows:

• •

  IPW (logit) and IPW (BART): We compute the asymptotic standard error assuming a logistic propensity model as in Cerulli (2014)

• •

  TMLE: We use the 95% confidence interval produced by the tmle R package

• •

  BCF: We construct a 95% credible interval using posterior samples of the treatment effect

We evaluate the results of our simulations using the following metrics.

• •

  RMSE: root mean squared error of estimated treatment effects

• •

  Bias: difference between true (simulated) ATE and average estimated ATE

• •

  Coverage: fraction of 95% confidence / credible intervals that contain the true ATE

For simulated covariates, outcomes, and treatment effects, we use a subset of the data generating processes tested in Hahn et al. (2020).

There are $p=5$ covariates, the first three of which are independent standard normal variables, the fourth of which is binary, and the fifth is an unordered categorical variable with values $1,2,3$. Treatment effects are homogeneous ($\tau=3$), and the outcome is determined by a combination of treatment effects and a piecewise interaction function of three of the covariates ($\mu(X)=1+g(X_{5})+X_{1}X_{3}$, where $g(X)=2$ if $X=1$, $g(X)=-1$ if $X=2$, and $g(X)=-4$ if $x=3$). $\mu(X)$ is often referred to as a prognostic effect and can be conceptualized as the expected value of the outcome for individuals who do not receive the treatment.

For treatment assignment, we consider two cases:

1. 1.

   $P(Z=1\mid X)=0.5$

2. 2.

   $P(Z=1\mid X)=0.8\Phi(3\mu(X)/s-0.5X_{1})+0.05+u/10$

   • •

     $\Phi(\cdot)$ is the standard normal CDF

   • •

     $s$ is the sample standard deviation of simulated observations of $\mu(X)$

   • •

     $u\sim\mathrm{Uniform}(0,1)$

The first treatment assignment mechanism mirrors a simple balanced randomized experiment. This presents relatively straightforward inference, as there is no confounding, but it allows us to simulate the phenomenon of including non-confounders in the propensity model. The second treatment assignment mechanism is described as “Targeted Selection” in Hahn et al. (2020), and refers to the phenomenon of assigning treatment based on the expected value of the outcome for those who don’t receive treatment.

In order to simulate the process by which outcomes are unlabeled, we assume a fixed overall sample size of 5,000 and randomly label 50, 100, or 500 observations so that $M$ is not correlated with $Y$, $X$, or $Z$. This parallels the notion of “missing completely at random” (MCAR) in the missing data literature (Little and Rubin (2002)), though our problem is slightly different. We treat our unlabeled $X^{m}$ and $Z^{m}$ as auxiliary data that may be helpful in estimating the ATE rather than treating our unobserved $Y^{m}$ as missing data to be imputed or otherwise estimated. Before proceeding to a more detailed discussion of the simulation results, we briefly summarize the key takeaways:

• •

  Using unlabeled data will harm inferences in the case of a mis-specified propensity model, as confidence intervals narrow around a biased estimate

• •

  BCF and IPW/BART work best in the case of targeted selection

• •

  While randomized treatment assignment makes ATE estimation much easier, unlabeled data can still be useful in this case for variance reduction purposes

We first examine the targeted selection case, in which treatment assignment is correlated with the prognostic effect.

We now turn to the randomized treatment assignment scenario.

Section 2.3 reviews the paradox of using estimated propensities to reduce variance in weighting estimators, even if true propensities are known (cited in Hirano, Imbens and Ridder (2003), Rosenbaum (1987), and Lunceford and Davidian (2004)). We believe this issue deserves a detailed discussion, as it contains a number of subtleties and its central claim (i.e. that estimated propensities are “better” than true propensities from a variance perspective) persists in the causal inference community to this day.

Though there may be benefits to using estimated propensity scores, in practice, calculating propensity scores introduces a new set of risks. The simulations in Section 3 show that a misspecified propensity model can bias estimates, sometimes drastically. It may seem natural, then, to treat a propensity model as any other supervised learning problem, with variable selection, model validation, and diagnostics to ensure a proper fit. Indeed, McCaffrey, Ridgeway and Morral (2004) use boosting to estimate propensities, in part for the automatic variable selection and flexible model specification that tree ensembles provide. Their results show that using boosting for a propensity model leads to better covariate balance and lower standard errors than a logistic propensity model.

However, Hernan and Robins (2020) caution that traditional model selection techniques may be of dubious utility in constructing propensity models. They note that the purpose of the propensity score is to control for confounding variables, not predict $Z$ from $X$ arbitrarily well. They are careful to point out that a purely data-driven approach to model selection runs the risk of including variables, such as colliders, mediators, or post-treatment variables, that jeopardize identification of treatment effects. They also advise that, even when treatment effects are identified, including instruments in a propensity model can increase the variance of treatment effect estimates by pushing estimated propensities closer to 0 and 1.

At first glance, this point would seem to argue against our recommendation to use semi-supervised machine learning in ATE estimation. Indeed, naively overfitting a propensity model is a surefire way to get noisy (and perhaps even biased) treatment effect estimates. However, we believe that a number of practical steps can be taken to mitigate these concerns.

First, estimated propensity scores should always be visualized or otherwise inspected to ensure that they are not numerically close to 0 and 1. Second, regularization methods (such as $\ell^{2}$ normalization or certain types of ensembles), can be used instead of logistic regression to minimize the risk of estimating propensity scores close to 0 and 1. Finally, subject matter expertise should always play a role in building propensity models. In many applied problems, especially those with a high-dimensional covariate space, the challenge of constructing an accurate causal graph is highly non-trivial. Insofar as subject matter expertise can identify variables such as colliders, post-treatment variables, or non-confounders, such variables should of course be excluded from any propensity model.

Researchers who are concerned about an underlying causal graph invalidating their propensity model can also turn to the causal discovery literature (i.e. Peters, Janzing and Schölkopf (2017)). Specifically, Entner, Hoyer and Spirtes (2013) present an approach to the problem of identifying confounders for adjustment in a causal model. In practice, the challenge of specifying (or estimating from the data) a causal graph can present a number of pitfalls beyond the scope of this paper. We suffice to say that there is no easy replacement for domain expertise in this process, but that this point applies to propensity models that use fully labelled data as well that those that use unlabeled data.

In order to better understand the pitfalls of including non-confounders in a propensity model, we test several possibilities in a simulation study. Letting $X_{j}$ be a variable in a set of $p$ covariates, we consider four cases:

1. a)

   Case 1: $X_{j}$ is a confounder and is included in a propensity model

2. b)

   Case 2: $X_{j}$ is not a confounder and is included in a propensity model

3. c)

   Case 3: $X_{j}$ is a confounder and is not included in a propensity model

4. d)

   Case 4: $X_{j}$ is not a confounder and is not included in a propensity model

Cases 1 and 4 are ideal, while cases 2 and 3 both fail to properly account for the true underlying causal model. The question we are concerned with is whether erring on the side of including non-confounders in a propensity model (case 2) is as harmful to inference as ignoring true confounders (case 3). To do this, we simulate 100 datasets of $n=$1,000 from the following data generating process:

For all variables except $X_{j}$, the $\beta$ and $\gamma$ coefficients for the outcome and propensity models are drawn uniformly at random from a range of $(-0.5,0.5)$ and $(-0.3,0.3)$, respectively. For the scenarios in which $X_{j}$ is a confounder (that is, $X_{j}$ impacts both $Y$ and $Z$), we set $\beta_{X_{j}}=0.5$. Each of our simulations assume some relationship between $\mathrm{x}_{j}$ and $\mathrm{z}$, and we vary $\gamma_{X_{j}}$ between $\{0,1,10\}$ to test the extent to which conditioning on $X_{j}$ in a propensity model separates the treatment and control groups.

Below we see the simulated bias, RMSE, and 95% interval coverage of each of the four cases, using BCF with a BART propensity model and $\gamma_{X_{j}}=1$

We see little difference in the results across each of the cases, because the influence of $X_{j}$ on $Y$ is modest in the confounded case. Below is the same set of outputs for $\gamma_{X_{j}}=10$.

We see that cases 1 and 4 exhibit the best performance (since they properly incorporate the true causal graph), but even when $X_{j}$ and $Z$ are strongly related, including $X_{j}$ in a propensity model when $X_{j}$ is not a confounder is less harmful than excluding $X_{j}$ when it is a confounder. To understand why, we briefly review the discussion of “regularization-induced confounding” (RIC) covered in Hahn et al. (2020).

As we have seen in our simulations, classic parametric models of a treatment assignment are vulnerable to misspecification. While this shortcoming suggests the use of flexible machine learning methods to estimate propensity scores, such methods run the risk of overfitting the data and estimating propensities close to 0 or 1. In addition to numeric instability of weighting estimators, such overfitting can also violate the positivity or conditional exchangeability assumption. Researchers who want to fit flexible propensity models while avoiding overfitting can use regularized machine learning methods, such as BART. Hahn et al. (2020), building on Hahn et al. (2018), show that a naive regularized model can bias treatment effect estimates. Their proposed solution, BCF, controls this bias by incorporating the estimated propensities as a feature in a regression model predicting $\mbox{E}(Y\mid Z=0,\hat{p}(X))$.

In our case, since the simulation studies above are conducted using BCF with a BART propensity model, we avoid some of the common problems of overfit propensity models while also recovering unbiased estimates of the average treatment effect. This explains why including a non-confounder, $X_{j}$, which is strongly predictive of $Z$ in the propensity model doesn’t harm inference of the ATE by nearly as much as excluding $X_{j}$ when it is a confounder.