RCT Rejection Sampling for Causal Estimation Evaluation

A causal model of a directed acyclic graph (causal DAG) ${\mathcal{G}}(V)$ can be viewed as the set of distributions induced by a system of structural equations: For each variable $V_{i}\in V$ there exists a structural equation $V_{i}\leftarrow f_{i}(\operatorname{pa}_{i},\epsilon_{i})$ (Pearl, 2009). This function maps the variable’s parents’ values—$\operatorname{pa}_{i}$ of $V_{i}$ in ${\mathcal{G}}(V)$—and an exogenous noise term³³3Typically these noise terms are assumed to be mutually independent, but this assumption is not strictly necessary (Richardson & Robins, 2013). Our results are non-parametric in the sense that we do not require any distributional assumptions on the noise terms or the specific form (e.g., linearity) of the structural equations $f_{i}(\cdot)$., $\epsilon_{i}$, to values of $V_{i}$. The system of equations induces a joint distribution $P(V)$ that is Markov relative to the DAG ${\mathcal{G}}(V)$, i.e., $P(V=v)=\prod_{V_{i}\in V}P(V_{i}=v_{i}\mid\operatorname{Pa}_{i}=\operatorname{pa}_{i})$. Independences in the distribution can be read off from ${\mathcal{G}}$ via the well-known d-separation criterion (Pearl, 2009). Interventions in the model are typically formalized using the do-operator (Pearl, 2009), where $Y|\operatorname{do}(T=t)$ denotes the value of an outcome $Y$ under an intervention that sets the treatment $T$ to value $t$.

Here, our causal estimand of interest is the average treatment effect (ATE) defined as,

where $t$ and $t^{\prime}$ denote distinct values of $T$. A causal parameter is said to be identified if it can be expressed as a function of the observed data $P(V)$. Given a set of variables $Z\subset V$ that satisfy the backdoor criterion w.r.t $T$ and $Y$⁴⁴4The set $Z$ satisfies the backdoor criterion if no variable in $Z$ is a causal descendant of $T$, and $Z$ blocks all backdoor paths between $T$ and $Y$, i.e., all paths of the form $T\leftarrow\cdots\rightarrow Y$., the ATE is identified via the well-known backdoor adjustment functional (Pearl, 1995).

We now describe the specific setup and objectives of RCT subsampling.⁵⁵5Here our target of interest is the ATE, though similar principles apply to other causal parameters like conditional average treatment effects. We start with a dataset $D_{\text{RCT}}$ consisting of $n$ iid draws from an RCT of pre-treatment covariates $C=\{C_{1},\dots,C_{k}\}$, treatment $T$, and outcome $Y$. Since this data is assumed to come from an RCT, the observed data distribution $P(C,T,Y)$ is Markov relative to the causal DAG shown in Fig. 1(a) where $T\perp\!\!\!\perp C$. The goal of RCT subsampling is to construct an observational dataset $D_{\text{OBS}}$ that consists of $m\leq n$ iid draws from $D_{\text{RCT}}$ that satisfies the following conditions which enable appropriate evaluation of causal estimation methods:

1. (I)

   Dependence induced. $D_{\text{OBS}}$ consists of samples drawn according to a new distribution $P^{*}(C,T,Y)$ that satisfies the dependence relation $T\not\perp\!\!\!\perp C$.

2. (II)

   ATE identified. There exists a functional $g$ of the RCT distribution and a functional $h$ of the subsampled data distribution such that $\text{ATE}=g(P(C,T,Y))=h(P^{*}(C,T,Y))$.

Here, (II) is an important identification pre-condition that ensures that it is possible, at least in theory, to compute estimates of the ATE from $D_{\text{OBS}}$ to match the ATE from $D_{\text{RCT}}$, the latter of which is treated as ground-truth in evaluation. From Fig. 1(a) it is clear that two sets of variables satisfy the backdoor criterion: the set $C$ and the empty set. Thus, the ATE is identified from the RCT distribution $P(C,T,Y)$ via the following two backdoor adjustment functionals,

Thus, a subsampling algorithm satisfies (II) if there is a functional $h(P^{*}(C,T,Y))$ that is equal to equation 2 or equation 3.

For our purposes, we add the condition (I) so that estimation in the observational data does not reduce to equation 3. That is, we aim to produce samples according to a distribution $P^{*}$ such that some adjustment is in fact necessary to produce unbiased ATE estimates. We note that (I) by itself is not sufficient to guarantee this; RCT subsampling procedures also require that there exists at least one pre-treatment covariate correlated with the outcome, i.e., $\exists\ C_{i}\in C$ such that $C_{i}\not\perp\!\!\!\perp Y$ in $P(C,T,Y)$ (Gentzel et al., 2021). However, this condition is easily testable, and we implement these checks in our synthetic experiments and real-world proof of concept (§4.2).

We now show a theoretical gap in existing approaches to subsampling RCTs, and propose a new algorithm that is theoretically guaranteed to satisfy conditions (I) and (II).

We claim that prior work that proposes RCT subsampling can result in observational samples from which the causal effect is not identified non-parametrically unless additional constraints are placed on the subsampling process. We consider Algorithm 2 in Gentzel et al. (2021) which does not explicitly impose such constraints and can be summarized as follows. Let $S$ be a binary variable indicating selection into the observational data from $D_{\text{RCT}}$. A structural equation $S\leftarrow\mathbbm{1}(T=\text{Bernoulli}(f(C)))$ is used to generate the selection variable, where $f$ is a function defined by the researcher and $\mathbbm{1}$ corresponds to the indicator function. $D_{\text{{OBS}}}$ is created by retaining only samples from $D_{\text{RCT}}$ where $S=1$. This results in $P^{*}(C,T,Y)=P(C,T,Y\mid S=1)$ which is Markov relative to the causal DAG in Fig. 1(b). From this DAG, it is easy to check via d-separation that condition (I) is satisfied as $T\not\perp\!\!\!\perp C\mid S=1$. However, the following proposition shows that condition (II) is not satisfied.

We provide a proof in Appendix A.

The intuition behind the proof of Proposition 3.1 is as follows. Identification of the ATE relies on two pieces: the conditional mean of the outcome given treatment and covariates and the marginal distribution of covariates. From Fig. 1(b), we have $\mathbb{E}[Y|T,C]=\mathbb{E}[Y|T,C,S=1]$, but $P(C)\not=P(C|S=1)$. Indeed this marginal distribution cannot be identified via any non-parametric functional of the subsampled distribution $P^{*}(C,T,Y)$ (Bareinboim & Tian, 2015). However, this non-identification result holds assuming that there is no additional knowledge/constraints on how $P^{*}$ is generated; in the next section we modify the sampling to place constraints on the generated distribution $P^{*}$ that mitigate this issue.

We propose Algorithm 1, which uses a rejection sampling procedure to subsample RCTs. Rejection sampling is useful when the target distribution is difficult to sample from but there exists a proposal distribution which is easier to sample from and the proposal distribution (times a constant) forms an “upper envelope” for the target distribution (Murphy, 2012, Chapter 23.2). Similar ideas on resampling data based on ratios of propensity scores appear in Thams et al. (2023) and Bhattacharya & Nabi (2022) in the context of testing independence constraints in post-intervention distributions. Though the rejection sampler also selects samples based on a function of $T$ and $C$, as in Fig. 1(b), we prove that additional constraints placed by the sampling strategy ensure identification holds in the new observed data distribution.

The intuition behind our algorithm is as follows. Sufficient constraints for maintaining identifiability of the ATE in $P^{*}(C,T,Y)$ via the functional in equation 2 are to ensure that $P^{*}(C)=P(C)$ and $P^{*}(Y\mid T,C)=P(Y\mid T,C)$.⁶⁶6One could also consider maintaining equality of just the conditional mean of $Y$ rather than the full conditional density. When this holds, it follows that equation 2 is equivalent to the adjustment functional $h(P^{*}(C,T,Y))=\sum_{c}P^{*}(c)\times(\mathbb{E}^{*}[Y\mid T=t,c]-\mathbb{E}^{*}[Y\mid T=t^{\prime},c])$, where $\mathbb{E}^{*}$ denotes the expectation taken w.r.t $P^{*}(Y\mid T,C)$. To also satisfy (I), we propose resampling with weights that modify $P(T)$ to a new conditional distribution $P^{*}(T\mid C)$.

The considerations listed in the prior paragraph inform our choice of an acceptance probability of $\frac{1}{M}\times\frac{P^{*}(T\mid C)}{P(T)}$ in the rejection sampler, where $M$ is the usual upper bound on the likelihood ratio used in the rejection sampler, which in our case is $\frac{P^{*}(T|C)}{P(T)}$.⁷⁷7In practice, we approximate $M$ from $D_{\text{RCT}}$ as $\frac{\max_{i}\in\{1,\dots,n\}P^{*}(T=t_{i}|C_{i})}{\min_{i\in\{1,\dots,n\}}\widehat{P}(T=t_{i})}$. Here, $P^{*}(T\mid C)$ is a function specified by the evaluation designer that satisfies positivity ($\forall c,0<P^{*}(T\mid C=c)<1$ almost surely), and is a non-trivial function of $C$ in the sense that $P^{*}(T\mid C)\not=P^{*}(T)$ for at least some values of $T$ and $C$.

Since $P^{*}$ satisfies $T\not\perp\!\!\!\perp C$ and yields identification via the usual adjustment functional obtained in a conditionally ignorable causal model, Algorithm 1 can be thought of as producing samples exhibiting confounding bias similar to the causal DAG in Fig. 1(c), despite the selection mechanism. A longer argument for this qualitative claim is in Appendix B.

We conclude this section by noting that similar to prior works on RCT subsampling algorithms, the subsampling strategy in Algorithm 1 only requires researchers to specify a single function, $P^{*}(T\mid C)$. Hence, our procedure satisfies our original desideratum of limited researcher degrees of freedom, while providing stronger theoretical guarantees for downstream empirical evaluation. However, specification of $P^{*}(T\mid C)$ may still be challenging when $C$ is high-dimensional. In Section 4.4, we discuss this finite data consideration and we use a proxy strategy for our proof of concept in which we have a low-dimensional confounding set $C$ along with a set of high-dimensional covariates $X$ that serve as proxies of this confounding set.

Using synthetic DGPs for $D_{\text{RCT}}$, we produce $D_{\text{OBS}}$ using Algorithm 2 from Gentzel et al. (2021) and separately via our RCT rejection sampler. We then compute ATE estimates using equation 2 for $D_{\text{OBS}}$ and compare it to the ground-truth estimates using equation 3 in $D_{\text{RCT}}$. Section C in the appendix gives the full details of the data-generating processes (DGPs) for three settings. Briefly, the DGP in Setting 1 has a single confounding covariate $C$, sets $P(T=1)=0.3$, and has an interaction term $TC$ in the structural equation for $Y$. Setting 2 is the same as Setting 1 except we set $P(T=1)=0.5$. Setting 3 is a non-linear DGP with five covariates, $C_{1},\dots,C_{5}$. All methods are provided with the true adjustment set and functional form for the outcome regression, i.e., our experiments here use oracle estimators to validate the identification theory proposed in the previous subsection.

We construct 95% confidence intervals via bootstrapping and the percentile method Wasserman (2004) and report confidence interval coverage. After obtaining a sample $D_{\text{OBS}}$ from the RCT via either RCT rejection sampling or Algorithm 2 from Gentzel et al., we resample $D_{\text{OBS}}$ with replacement and calculate the ATE for that bootstrap sample. We repeat this for 1000 bootstrap samples and obtain a 95% confidence interval by taking the 2.5% and 97.5% points of the bootstrap distribution as the endpoints of the confidence interval. Across our 1000 random seeds of the synthetic DGPs, we measure the proportion of these confidence intervals that contain the true (RCT) ATE and report this metric as confidence interval coverage. See Appendix Section H for additional confidence interval plots.

Table 2 shows that our proposed RCT rejection sampler results in a reduction of absolute bias compared to Algorithm 2 from  Gentzel et al. (2021) by a factor of over 24 for Setting 1 (0.22/0.009) and a factor of 21 in Setting 3 (0.252/0.012). For Setting 3, Gentzel et al.’s procedure results in almost a 100% increase in bias relative to the gold RCT ATE of $-0.26$. In Setting 2 where $P(T=1)=0.5$, the differences in absolute bias between the two algorithms is less pronounced.⁸⁸8We note Gentzel et al. (2021) primarily focus on the setting for which $P(T=1)=P(T=0)=0.5$. However, their approach does not seem to generalize well outside of this setting, theoretically and empirically. In Settings 1 and 2, the confidence interval coverage for Gentzel et al.’s procedure is 0 whereas our RCT rejection sampling algorithm results in coverage of 0.97 and 0.98 for Settings 1 and 2 respectively, both slightly above the nominal 0.95 coverage. The results of the simulation are consistent with our theoretical findings that our algorithm permits identifiability under more general settings than prior work.

A necessary component for RCT subsampling is obtaining a real-world RCT dataset. This ensures a more realistic data generating processes compared to synthetic or semi-synthetic approaches (see Table 1). As Gentzel et al. (2021) note, there are many RCT repositories from a variety of disciplines from which evaluation designers could gather data. However, Gentzel et al. find many of these existing datasets only have one or two covariates that satisfy $C\not\perp\!\!\!\perp Y$ (see Consideration #1 below).

As we briefly mentioned in Section 1, for settings with just a few covariates one can often use simple estimation strategies with theoretical guarantees—e.g., parametric models or contingency tables—and empirical evaluation may not be particularly informative in this setting. Along these lines, we recommend that evaluation designers first ask themselves, Is empirical evaluation of causal estimators appropriate and necessary for this setting? Not all settings are in need of empirical evaluation.

A potentially untapped resource for RCT rejection sampling data is A/B tests from large online platforms. Other work, e.g., Eckles & Bakshy (2021), have used these types of experiments for constructed observational studies and we use such a dataset in our proof of concept. The large scale of these experiments can be advantageous since RCT subsampling reduces the number of units in the observational dataset by roughly half. Further, they often contain rich metadata and many covariates, which can be used to induce confounding in a way that emulates a high-dimensional setting.

As we mentioned in Section 3, a necessary precondition for RCT subsampling in general is the existence of a causal edge between $C$ and $Y$, implying $C\not\perp\!\!\!\perp Y$. The relationship between $C$ and $Y$ is naturally occurring (not modified by evaluation designers) and the amount of confounding induced by sampling is, in part, contingent on this relationship (Gentzel et al., 2021). One can empirically check $C\not\perp\!\!\!\perp Y$ via independence tests, e.g., evaluating the odds ratio when both $C$ and $Y$ are binary variables. If there do not exist covariates, $C$ such that $C\not\perp\!\!\!\perp Y$, one cannot move forward in the evaluation pipeline using RCT subsampling.

Evaluation designers using RCT rejection sampling have one degree of freedom: specification of $P^{*}(T|C)$. We describe one specific approach in our proof of concept to choose $P^{*}(T|C)$, but we anticipate evaluation designers using RCT rejection sampling to create large empirical benchmarks may want to include many different parameterizations of $P^{*}(T|C)$ to evaluate empirical performance of methods under numerous settings. Consideration #3 describes approaches to diagnosing the choice of $P^{*}(T|C)$ for a specific finite dataset.

We recommend running diagnostics on the sampling procedure because, although Theorem 3.2 proves that our RCT rejection sampling permits identification, this is an asymptotic statement, and the sampler can produce finite samples where confounding bias is induced, but difficult to correct for.

The primary goal of this work is to create clear steps to follow during the evaluation design phase. Although this stage precedes a thorough modeling effort, we recommend that one runs baseline models to check for potential issues.

The broader purpose of this line of empirical evaluation is to understand the real-world settings for which certain estimation approaches are successful or unsuccessful. We believe bridging the gap between theory that holds asymptotically and finite data is important for drawing valid causal inference, but evaluation designers might encounter unforeseen challenges particular to finite data that need to be examined carefully.