Learning Who to Treat When Treatment is Missing
(Supplementary Material)
July 15, 2026
Policy learning methods are increasingly used to inform treatment allocation under budget constraints. Most proposed methods assume complete treatment data, yet applications frequently suffer from missingness that can bias estimates and lead to suboptimal policies. We address this gap by extending efficient estimators for average treatment effect (ATE) estimation to policy value and conditional average treatment effect (CATE) estimation under missing at random (MAR) and missing completely conditionally at random (MCCAR) treatment data. Through asymptotic efficiency analysis, we prove that the MAR estimator, which leverages partially-observed units, is both valid and more efficient than the MCCAR estimator when MCCAR assumptions hold. This result provides formal justification for preferring MAR-based estimation in policy learning under both missing data settings. Our comprehensive experiments using synthetic and semi-synthetic datasets confirm that correctly specifying the missingness mechanism is crucial: misspecified estimators remain biased regardless of sample size, while our estimators achieve near-oracle performance when assumptions are satisfied. Our work provides practitioners with theoretically grounded, empirically validated tools for robust policy learning in the presence of missing treatment data.
Machine learning is increasingly used to allocate limited treatments between individuals to improve outcomes across domains such as healthcare [1], [2] and social services [3], [4]. In these settings, one potential allocation strategy is to treat those who will benefit most from treatment by estimating the conditional average treatment effect (CATE). While many CATE estimators have been proposed, nearly all methods assume complete data (e.g. [5]–[8]). This assumption is frequently violated, for example due to incomplete records, measurement errors, and data integration challenges [9], [10]. Though recent work addresses missing outcomes [11], in this work, we focus on the challenge of missing treatment data.
Missing treatment data is prevalent in observational studies. For example, [12] found that cardiovascular study subjects often had incomplete physical activity information due to missed follow-ups caused by poor health, while [13] encountered substantial missing alcohol exposure data when studying tuberculosis infection in Uganda. Similarly, [14] found substantial missingness in body mass index (BMI), but not in other covariates or outcomes, when studying the effect of maternal BMI on birthweight. These examples reflect common missingness patterns in linked administrative and self-reported survey data [9], [15]. In these settings we would not expect missingness to be completely random.
Rather, we may assume that missingness depends on other observed variables. If missingness is independent of treatment and observed outcomes given covariates, we say that treatment is Missing Completely Conditionally at Random (MCCAR). Under MCCAR assumptions, the identified CATE is a complete-case quantity, discarding observations with missing treatment. Although valid, we show this default is suboptimal. The source of this suboptimality becomes clear once we consider a less restrictive missingness mechanism. A more efficient approach is to allow missingness to also depend on observed outcomes, which we call Missing at Random (MAR). Under MAR, the identified CATE leverages partially-observed units and, since MCCAR is a special case of MAR, remains valid under MCCAR assumptions.
Several prior works have developed causal estimation methods under the MAR assumption. While early work addressed average treatment effect (ATE) estimation with semiparametric estimators [14], [16], these approaches rely on parametric nuisance function estimation. Recently, [17] developed efficient, fully nonparametric estimators for the ATE under MAR assumptions, and [13] studied longitudinal exposure effects under both missing treatments and outcomes. For CATE estimation, [18] proposed MTRNET, which uses adversarial learning to balance representations across two covariate shifts (treated-control, observed-missing) under the assumption that missingness is jointly independent of potential outcomes and treatment conditional on observed covariates. In practice, missingness may also be related to observed outcomes. Furthermore, given the difficulty of learning reliable CATE estimates from real-world data [19], flexible nonparametric CATE estimation methods that handle both MCCAR and MAR missingness are desirable. To address these challenges, our work makes the following contributions:
We extend efficient influence function methods from ATE to policy value and CATE estimation under both MAR and MCCAR assumptions.
We prove that the MAR estimator is more efficient than the MCCAR estimator whenever both are valid, even when the stronger MCCAR assumption holds. This provides a formal argument against complete-case analysis, and a theoretical basis for preferring MAR-based methods in practice.
We present comprehensive empirical validation of our theoretical results on synthetic datasets and real datasets with generated missingness.
Let \(Z_1, \dots, Z_n \sim \mathbb{P}\) be an iid sample from population distribution \(\mathbb{P} \in \mathcal{P}\) with \(Z_i = (\mathbf{X}_i, R_i, R_iA_i, Y_i)\) where \(\mathcal{P}\) is a statistical model. Here, \(\mathbf{X}\in \mathbf{\mathcal{X}} \subseteq \mathbb{R}^d\) represents an observed vector of covariates, \(A \in \{0, 1\}\) indicates treatment receipt, \(Y \in \mathcal{Y} \subseteq \mathbb{R}\) is the observed outcome, and \(R \in \{0, 1\}\) indicates whether treatment assignment is observed (\(R=1\)) or missing (\(R=0\)). Therefore, for an individual \(i\), treatment is only observed when \(R_i = 1\). When \(R_i = 0\), the product \(R_i A_i =0\) regardless of the true value of \(A_i\). We denote generic treatment decision policies by \(d: \mathbf{X} \mapsto \{0,1\}\), and we use \(\mathcal{D}\) to denote a class of decision rules. Following the potential outcomes framework, we use \(Y(d(\mathbf{X}))\) and \(Y(A)\) to denote the potential outcome observed by setting treatment according to a deterministic treatment policy \(d\) and under treatment \(A\), respectively [20]. Lastly, \(\kappa \in [0,1]\) represents the budget constraint on the proportion of the population able to be treated.
Throughout the paper, we use \(\mathbb{E}[\cdot]\) to represent expectation with respect to \(\mathbb{P}\), \(\mathbb{P}_n\) to denote sample averages, and \(\widehat{\mathbb{P}}\) for an estimate of the distribution. We place hats over quantities estimated using the observed data, i.e. \(\hat{\psi} = \psi(\widehat{\mathbb{P}})\). \(\widehat{\mathbb{E}}\) is a generic estimated conditional expectation. \(f_n = o_{\mathbb{P}}(r_n)\) denotes that \(f_n / r_n\) converges in probability to zero. Lastly, \(\| f \|^2 = \int f(z)^{2}d\mathbb{P}(z)\) is the squared \(L_2(\mathbb{P})\) norm.
We consider a setting in which a policymaker, constrained by a treatment budget \(\kappa\), seeks to design a treatment assignment policy \(d(\mathbf{X})\) maximizing utility \(V(d) = \mathbb{E} \left[Y(d)\right]\). The policymaker has access to observational data from potentially suboptimal historical treatment assignments. Formally, the policymaker’s objective may be summarized as \[\begin{align} \label{eq:optobjective} \max_{d \in \mathcal{D}} \mathbb{E}[Y(d(\mathbf{X}))] \text{ s.t. } \mathbb{E}[d(\mathbf{X})] \leq \kappa. \end{align}\tag{1}\] The optimal solution to this objective is to treat individuals with the highest treatment benefit until the budget is exhausted: \(d^*_{\kappa}(\mathbf{X}) = \mathbb{1}(\tau(\mathbf{X}) > \tau^*_{\kappa})\), where \(\tau(\mathbf{X}) = \mathbb{E} \left[ Y(1) - Y(0) \mid \mathbf{X}\right]\) is the CATE and \(\tau^*_{\kappa} = \max\{\inf{ \tau: \mathbb{P}(\tau(\mathbf{X}) \geq \tau) \leq \kappa}, 0\}\) represents the minimum benefit threshold for treatment assignment [21], [22]. When \(\tau^* = 0\), the budget constraint is non-binding and all individuals who benefit receive treatment. It is well known that the utility under policy \(d^*(\mathbf{X})\) is \[\begin{align} V(d^*) &= \mathbb{E}[Y(1)d^*(\mathbf{X}) + Y(0)(1 - d^*(\mathbf{X}))]. \end{align}\] With fully observed treatments, the standard causal assumptions are sufficient for identification of \(V(d^*)\) and \(\tau(\mathbf{X})\),
Assumption 1. Consistency: If \(A = a\), then \(Y(a) = Y\). A-Ignorability: \(Y(1), Y(0) \mathpalette{\independenT}{\perp}A \mid \mathbf{X}\). A-Positivity: \(0 < P(A = 1 \mid \mathbf{X}= \mathbf{x}) < 1\).
Assumption 1 ensures that the observed outcome is the potential outcome under the observed treatment, there is no confounding between potential outcomes and treatment assignments, and lastly all individuals have some chance at receiving treatment. Under Assumption 1 and resource constraint \(\kappa\), [21] studied the CATE and expected utility under \(d^*\): \[\begin{align} \tag{2} &\tau(\mathbf{X}) = \mathbb{E}[Y \mid \mathbf{X}, A = 1] - \mathbb{E}[Y \mid \mathbf{X}, A = 0] \\ \tag{3} &V(d^*(\mathbf{X})) = \mathbb{E}[\mathbb{E}[Y \mid \mathbf{X}, A =1]d^*(\mathbf{X}) \\ &\qquad+ \mathbb{E}[Y \mid \mathbf{X}, A=0](1 - d^*(\mathbf{X}))] \nonumber. \end{align}\] When treatment assignments are missing for some observations, it is not possible to condition on treatment assignment \(A\). Consequently, the CATE \(\tau(\mathbf{X})\) (2 ) and value function \(V(d^*(\mathbf{X}))\) (3 ) cannot be recovered from observed data and are therefore not identified.
In many applications relying on observational data, past treatment decisions may be unobserved for a portion of the population. If treatment missingness occurs completely independently of observed or unobserved covariates, we say that treatment is Missing Completely at Random (MCAR). This is a strong assumption and is unlikely to be satisfied in practice. To relax this assumption, we say that treatment is Missing Completely Conditionally at Random (MCCAR) when missingness depends only on the observed pre-treatment covariates. For example, MCCAR may occur in electronic health records when there are systematic differences in treatment observation between patients explained by observed covariates, such as if insured patients have more complete treatment records than uninsured patients. Below, we complement Assumption \(\ref{assum:standard}\) with assumptions on the missingness mechanism in the MCCAR setting:
Assumption 2. (MCCAR) R-Ignorability: \(R \mathpalette{\independenT}{\perp}(A, Y) \mid \mathbf{X}\). \(R-Positivity\): \(0 < P(R =1 | \mathbf{X}= \mathbf{x})\).
Under Assumption 2, missingness is conditionally independent of treatment given observed covariates and all individuals have some probability of treatment observation. In practice, however, missingness may also be dependent on the observed outcomes \(Y\). Outcome-dependent missingness may occur if some records are lost or fail to be matched, or when survey non-response is present. In this case, we say treatment is Missing at Random (MAR) as missingness depends on only fully observed variables [23], [24]. For example, in a study of the effect of maternal BMI on infant birthweight, the rich set of observed covariates and outcomes makes it plausible that the substantial missingness in BMI is explained by observed factors alone [14]. [16], [14], and [17] all previously studied ATE estimation in this setting under the standard Assumptions \(\ref{assum:standard}\) in addition to the MAR assumption:
Assumption 3. (MAR) R-Ignorability: \(R \mathpalette{\independenT}{\perp}A \mid (\mathbf{X}, Y)\). \(R-Positivity\): \(0 < P(R =1 | \mathbf{X}= \mathbf{x}, Y = y)\).
Assumption 3 excludes Missing Not at Random (MNAR) treatment, which would occur if missingness were dependent on unobserved factors associated with treatment status itself even after conditioning on \(\mathbf{X}\) and \(Y\). For example, MNAR missingness may occur if participants who were less likely to receive treatment for reasons not captured by recorded covariates were also less likely to have their treatment status documented. As treatment effects under MNAR are nonidentifiable without further assumptions, we restrict attention to MCCAR and MAR going forward.
As Fig. 1 shows, in the MAR setting missingness can depend on the outcome Y, whereas under MCCAR (left panel) no such path exists since missingness depends only on observed covariates. Fig. 1 also illustrates that MCCAR R-Ignorability (2) implies MAR R-Ignorability (3). Therefore, MCCAR R-Ignorability is strictly stronger than MAR R-Ignorability.
Although MAR R-Ignorability is strictly weaker than MCCAR R-Ignorability, the MAR R-Positivity condition is possibly stronger, as it requires conditioning on \(Y\) in addition to the observed covariates \(X\), as summarized in Table 1. However, under MCCAR Assumption 2, MCCAR R-Ignorability and MCCAR R-Positivity together imply MAR R-Positivity. Therefore, MAR Assumption 3 is strictly weaker than MCCAR Assumption 2.
Our two sets of assumptions lead to two identification strategies. Because MCCAR R-Ignorability implies MAR R-Ignorability, the MAR identification approach is valid in both settings. See Appendix 6.1 for identification proofs.
Under the standard Assumptions 1 and either MCCAR Assumption 2 or MAR Assumption 3, [17] identified the ATE using the following nuisance functions: \[\begin{align} \lambda_{d}(\mathbf{X}, Y) &= P(A = d(\mathbf{X}) \mid \mathbf{X}, Y, R = 1) \\ \pi(\mathbf{X}, Y) &= P(R = 1 \mid \mathbf{X}, Y) \\ \beta_{d}(\mathbf{X}) &= \mathbb{E}\left[Y\lambda_{d}(\mathbf{X}, Y) \mid \mathbf{X}\right] \\ \gamma_{d}(\mathbf{X}) &= \mathbb{E}\left[\lambda_{d}(\mathbf{X}, Y) \mid \mathbf{X}\right] = P(A = d(\mathbf{X}) \mid \mathbf{X}). \end{align}\] Here, \(\gamma_{d}\) and \(\lambda_{d}\) represent the probability the observational treatment assignment matches the recommendation under policy \(d\) with and without outcome information, \(\pi\) is the treatment observation probability, and \(\beta_{d}\) is a propensity-weighted outcome expectation. We can identify the value function under policy \(d\) as \[\begin{align} \label{eq:maroutcomeestimator} V(d) = \mathbb{E}\left[Y(d(\mathbf{X}))\right]= \mathbb{E}\left[\frac{\beta_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})}\right] = \psi_{\text{MAR}}. \end{align}\tag{4}\] Importantly, \(\psi_{\text{MAR}}\) (4 ) utilizes all data by reweighting observed outcomes by probability of treatment \(\lambda_{d}\). The CATE is then identified as \[\begin{align} \label{eq:marcateestimator} \tau(\mathbf{X}) = \frac{\beta_1(\mathbf{X})}{\gamma_1(\mathbf{X})} - \frac{\beta_0(\mathbf{X})}{\gamma_0(\mathbf{X})}. \end{align}\tag{5}\]
The MCCAR Assumption 2 permits a second identification strategy. We define \(\nu_{d}(\mathbf{X})= \mathbb{E}[Y | \mathbf{X}, A = a, R =1]\) as the outcome function among the observed data. Then, we may write the expected utility under policy \(d\) as \[\begin{align} \label{eq:mccarval} V(d) = \mathbb{E}\left[Y(d(\mathbf{X}))\right]= \mathbb{E}\left[ \nu_{d}(\mathbf{X}) \right] = \Phi_{\text{MCCAR}}. \end{align}\tag{6}\] \(\Phi_{\text{MCCAR}}\) (6 ) is a complete-case estimand as only observations without missing data are used. The same strategy yields an identified quantity for the CATE: \[\begin{align} \label{eq:cccate95estimand} \tau(\mathbf{X}) &= \nu_1(\mathbf{X}) - \nu_0(\mathbf{X}). \end{align}\tag{7}\]
Unlike Strategy One, this quantity is only valid under MCCAR Assumption 2. The first strategy, however, is valid under both MCCAR and MAR Assumption 3. If the MCCAR assumption applies, both Strategy One and Strategy Two identify the same population CATE. As we will show, meaningful differences between the two strategies emerge in estimation, particularly in terms of efficiency and robustness to nuisance function misspecification. We now proceed to deriving efficient and consistent estimators of these identified quantities.
| MCCAR | MAR | |
|---|---|---|
| R-Ignorability | Stronger | Weaker |
| R-Positivity | Weaker | Stronger\(^*\) |
| Valid IFs | \(\ifMCCAR\), \(\ifMAR\) | \(\ifMAR\) |
Our goal is to learn optimal treatment policies \(d^*\) that maximize expected outcomes subject to budget constraints. This requires efficient methods for both policy evaluation and treatment effect estimation.
To derive our estimators, we use Influence Functions (IFs) (also known as Neyman Orthogonal scores), which yield bias-corrected estimators that achieve the semiparametric efficiency bound and often possess double robustness (double ML) properties [24]–[27]. Doubly robust estimators remain consistent when some model components are misspecified, which is a valuable property in real-world settings with unknown nuisance functions, and can converge faster than their nuisance estimators. We provide technical details and all proofs in Appendix 6.2.
Remark 4. We implicitly assume \(Y\) is bounded so our estimators have finite variance. We also assume that all nuisance functions are estimated from an independent sample, a common procedure known as sample splitting, which is easily adaptable to cross-fitting. [27]–[29].
In this section, we establish nonparametric methods for evaluating the performance of a given policy \(d\). Under each missingness assumption, we construct \(\sqrt n\)-consistent, asymptotically normal estimators. We then compare their asymptotic efficiency under MCCAR Assumption \(\ref{assum:mccar}\), where both are valid.
We consider two statistical models: \(\mathcal{P}_{MAR}\), all distributions on \(Z\) satisfying Assumptions 1 and 3, and \(\mathcal{P}_{MCCAR}\), all distributions on \(Z\) satisfying Assumptions 1 and 2. Under \(\mathcal{P}_{MCCAR}\), \(V(d)\) admits uncentered IFs \(\phi_{\text{MCCAR}}\) and \(\varphi_{\text{MAR}}\); under \(\mathcal{P}_{MAR}\), only \(\varphi_{\text{MAR}}\) is valid: \[\begin{align} &\varphi_{\text{MAR}}(d) = \left(\frac{Y - \frac{\beta_{d}(\mathbf{X})}{\gamma_{d}{(\mathbf{X})}}}{\gamma_{d}{(\mathbf{X})}}\right)\biggr(\frac{R(\mathbb{1}(A = d(\mathbf{X})) - \lambda_{d}(\mathbf{X}))}{\pi(\mathbf{X}, Y)} \\ &\qquad \qquad \qquad + \lambda_{d}(\mathbf{X}, Y)\biggr) + \frac{\beta_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})} \nonumber \\ &\phi_{\text{MCCAR}}(d) = \frac{R(\mathbb{1}(A = d(\mathbf{X}))}{\eta_{d}(\mathbf{X})}\left(Y - \nu_{d}(\mathbf{X})\right) + \nu_{d}(\mathbf{X}) \end{align}\] where \(\eta_{d}= P(A = d(\mathbf{X}), R = 1 \mid \mathbf{X})\). We now establish convergence properties of estimators based on these IFs. We first present our doubly robust, nonparametric estimator for the MCCAR setting and establish its asymptotic properties:
Proposition 5. Define the bias-corrected estimator derived from \(\phi_{\text{MCCAR}}\) as \(\widehat{\Phi}_{\text{MCCAR}}= \mathbb{P}_n(\widehat{\phi}_{\text{MCCAR}})\). Assume 1, 2, in addition to:
\(\| \widehat \eta_{d}- \eta_{d}\|\, \| \widehat \nu_{d}- \nu_{d}\| = o_{\mathbb{P}}(n^{-1/2})\)
\(\| \widehat{\phi}_{\text{MCCAR}}- \phi_{\text{MCCAR}}\| = o_{\mathbb{P}}(1)\)
Then, \(\sqrt n(\widehat{\Phi}_{\text{MCCAR}}- V(d)) \rightsquigarrow \mathcal{N}(0, \textrm{Var}(\phi_{\text{MCCAR}})),\) with \(\textrm{Var}(\phi_{\text{MCCAR}}) =\) \[\mathbb{E}\left[\frac{\sigma_{d}^{2}(\mathbf{X})}{\eta_{d}(\mathbf{X})}\right] + \mathbb{E}\left[\left( \nu_{d}(\mathbf{X}) - \Phi_{\text{MCCAR}}\right)^{2}\right] ,\] where \(\sigma_{d}^{2}(\mathbf{X}) = \textrm{Var}(Y \mid \mathbf{X}, A=d(\mathbf{X}), R = 1)\).
The conditions of Proposition 5 establish asymptotic normality under a product-rate condition on the nuisance functions and consistent estimation of \(\phi_{\text{MCCAR}}\) at any rate. In particular, the latter empirical process condition rules out nuisance estimators that diverge, even when the product-rate condition is satisfied. Importantly, Proposition 5 establishes that the estimator \(\widehat{\Phi}_{\text{MCCAR}}\) is doubly robust, achieving consistency when either \(\widehat \nu_{d}\) or \(\widehat \eta_{d}\) is correctly specified.
In the MAR setting, \(\phi_{MCCAR}\) is not a valid IF. We therefore adapt the IF from [17], originally developed for ATE estimation, to policy value estimation and derive its asymptotic variance:
Proposition 6. Define the bias-corrected estimator derived from \(\varphi_{\text{MAR}}\) as \(\widehat{\psi}_{\text{MAR}}= \mathbb{P}_n(\widehat{\varphi}_{\text{MAR}})\). Assume 1 and either 3 or 2, as well as:
\(\| \widehat \pi- \pi\| \| \widehat \lambda_{d}- \lambda_{d}\| = o_{\mathbb{P}}(n^{-1/2})\)
\(\| \widehat \beta_{d}- \beta_{d}\| \| \widehat \gamma_{d}- \gamma_{d}\| = o_{\mathbb{P}}(n^{-1/2})\)
\(\| \widehat \gamma_{d}- \gamma_{d}\| = o_{\mathbb{P}}(n^{-1/4})\)
\(\| \widehat{\varphi}_{\text{MAR}}- \varphi_{\text{MAR}}\| = o_{\mathbb{P}}(1)\)
Then, \(\sqrt n(\widehat{\psi}_{\text{MAR}}- V(d)) \rightsquigarrow \mathcal{N}(0, \textrm{Var}(\varphi_{\text{MAR}}))\) and \[\begin{align} &\textrm{Var}(\varphi_{\text{MAR}}) = \mathbb{E}\biggr[\biggr(\frac{Y - \beta_{d}(\mathbf{X}) / \gamma_{d}(\mathbf{X}) }{\gamma_{d}(\mathbf{X})}\biggr)^{2} \\& \quad \times \biggr(\frac{\lambda_{d}(\mathbf{X},Y)}{\pi(\mathbf{X},Y)} - \lambda_{d}^{2}(\mathbf{X},Y)\cdot \biggr\{\frac{1 - \pi(\mathbf{X},Y)}{\pi(\mathbf{X},Y)}\biggr\}\biggr)\biggr] \\& \quad + \mathbb{E}\biggr[\biggr(\frac{\beta_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})} - \psi_{\text{MAR}}\biggr)^{2}\biggr]. \end{align}\]
Proposition 6 demonstrates that \(\widehat{\psi}_{\text{MAR}}\) achieves asymptotic normality under product-rate conditions on the nuisance functions and consistent estimation of \(\varphi_{\text{MAR}}\). As in Proposition 5, the empirical process condition requires consistent estimation of all nuisance functions at any rate. Proposition 6 also reveals that consistency is achieved when \(\widehat \gamma_{d}\) is estimated consistently, along with either \(\widehat \lambda_{d}\) or \(\widehat \pi\). Unlike standard doubly robust estimators, \(\widehat{\psi}_{\text{MAR}}\) requires \(\widehat \gamma_{d}\) to be consistently estimated, which is expected given that \(\varphi_{\text{MAR}}\) is nonlinear in \(\gamma_{d}\).
Since \(\varphi_{\text{MAR}}\) is the IF under a nonparametric model, it must be the efficient IF under Assumption \(\ref{assum:mar}\) [26]. Unlike MAR Assumption \(\ref{assum:mar}\), MCCAR Assumption \(\ref{assum:mccar}\) defines a semiparametric model, which is a subset of the nonparametric model with the additional, testable condition \(R \mathpalette{\independenT}{\perp}Y \mid X\). Under this model, both estimators are valid, motivating a comparison of their efficiency.
For this comparison, we use the Asymptotic Relative Efficiency
\(\text{ARE}(f_1, f_2) := \frac{\textrm{Var}(f_1)}{\textrm{Var}(f_2)}\) , where \(f_1\) and \(f_2\) are asymptotically normal estimators of the same quantity. When \(\textrm{ARE} > 1\), \(f_2\) is more efficient.
Theorem 7. Under Assumptions \(\ref{assum:standard}\) and \(\ref{assum:mccar}\) where \(\widehat{\psi}_{\text{MAR}}\) and \(\widehat{\Phi}_{\text{MCCAR}}\) are both asymptotically normal estimators of \(V(d)\), \(\widehat{\psi}_{\text{MAR}}\) is the more efficient estimator with asymptotic relative efficiency \[\left( 1 - \frac{ \mathbb{E}\left[\left(\frac{Y - \nu_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})}\right)^{2}\cdot \lambda_{d}(\mathbf{X}, Y)^{2} \cdot \left(\frac{1-\pi(\mathbf{X})}{\pi(\mathbf{X})}\right)\right] }{\textrm{Var}(\phi_{\text{MCCAR}})}\right)^{-1}.\]
Theorem 7 establishes that the complete-case-based MCCAR estimator is not the efficient choice when the conditions of both Propositions 5 and 6 are met, demonstrating that the MAR estimator is the preferred choice under those conditions. This efficiency gain, however, is contingent upon the conditions of Proposition 6, notably that \(\widehat \gamma_{d}\) is estimated consistently, though at nonparametric rates.
As Theorem 7 shows, the efficiency gain increases with the augmented propensity score (\(\lambda_{d}\)), missingness rate (\(1-\pi\)), and the covariance between these terms. To demonstrate these efficiency gains, we conduct a synthetic experiment where we systematically vary these factors. The experimental results in Fig. 2 strongly confirm the predictions of Theorem 7. For example, ARE increases from 1.5 to 2.0 by increasing \(\lambda_d\) from 45% to \(70\%\). Additional experimental details are provided in Appendix 7.
With policy evaluation methods in place, we now consider the complementary challenge of learning optimal treatment assignments \(d^*_\kappa(\mathbf{X})\) by estimating the CATE function \(\tau(\mathbf{X})\) using the DR-Learner [30], [31]. The DR-Learner enables direct adaptation of our IF-based value function estimators to doubly-robust CATE estimation, requiring only covariates at inference time, and permits the use of generic machine learning methods for both nuisance function and CATE estimation, thus providing flexibility that approaches built on specific methods lack (such as [18], [32]).
We first define pseudo-outcomes based on the IFs of the value functions: \(\widetilde{\varphi}_{\textrm{MAR}} = \widehat{\varphi}_{\text{MAR}}(1) - \widehat{\varphi}_{\text{MAR}}(0)\) and \(\widetilde{\phi}_{\textrm{MCCAR}} = \widehat{\phi}_{\text{MCCAR}}(1) - \widehat{\phi}_{\text{MCCAR}}(0)\). When we regress these pseudo-outcomes on covariates \(\mathbf{X}\), the resulting CATE estimator exploits the robustness properties of the pseudo-outcome construction. Under additional smoothness or sparsity conditions on the CATE function, this approach can achieve faster convergence rates than methods that estimate the difference in outcome regression functions directly [31]. We propose two DR-Learner variants corresponding to our MAR and MCCAR assumptions: \[\begin{align} \tag{8} \widehat \tau_{\textrm{MAR}}&= \widehat{\mathbb{E}} \left[\widetilde{\varphi}_{\textrm{MAR}} \mid \mathbf{X}\right] \\ \tag{9} \widehat \tau_{\textrm{MCCAR}}&= \widehat{\mathbb{E}}\left[\widetilde{\phi}_{\textrm{MCCAR}} \mid \mathbf{X}\right] \end{align}\] While \(\widehat{\tau}_{MCCAR}\) follows standard two-way sample splitting, \(\widehat{\tau}_{MAR}\) requires three-way splitting for continuous \(Y\) to avoid conditional density estimation when constructing \(\widehat{\beta}_{d}\) and \(\widehat{\gamma}_d\). Instead, we first fit \(\widehat{\lambda}_d(\mathbf{X}, Y)\). Using the second fold, we regress \(Y\widehat{\lambda}_d(\mathbf{X}, Y)\) on \(\mathbf{X}\) to obtain \(\widehat{\beta}_{d}(\mathbf{X})\), and separately regress \(\widehat{\lambda}_d(\mathbf{X}, Y)\) on \(\mathbf{X}\) to obtain \(\widehat {\gamma}_d(\mathbf{X})\). The last fold is reserved to estimate the CATE function \(\widehat \tau_{\textrm{MAR}}\). Example algorithms for both estimators are detailed in Appendix 8.
Both CATE estimates can then be used directly in the budget-constrained policy \(\widehat{d}^*_\kappa(\mathbf{X}) = \mathbb{1}(\widehat{\tau}(\mathbf{X}) > \widehat{\tau}^*_\kappa)\), where the threshold \(\widehat{\tau}^*_\kappa\) is chosen to satisfy the budget constraint \(\mathbb{E}[\hat{d}^*_\kappa(\mathbf{X})] \leq \kappa\) [21], [22].
With our CATE estimation strategy in place, we now verify our estimators’ performance on both synthetic and semi-synthetic experiments.
Remark 8. Convergence rates of \(\widehat \tau_{\textrm{MAR}}\) and \(\widehat\tau_{\textrm{MCCAR}}\) can be established following the arguments of [17], [8], and [30] using the results of Propositions 5 and 6.
Remark 9. Asymptotic normality of the estimated policy value \(\widehat V\) at the data-dependent optimal policy \(\widehat d^*_\kappa\) can be established under a margin condition on the CATE near threshold \(\tau_{\kappa}^{^*}\) [33]. See Theorems 15 and 14 in Appendix 6.3.
We evaluate our two CATE estimators in two experimental settings. First, we use synthetic data to verify their theoretical convergence properties. We then use semi-synthetic data from two randomized trials to evaluate policy performance.1.
We conduct two synthetic experiments with known potential outcomes to evaluate our estimators. First, we verify the robustness properties of our estimators by systematically misspecifying nuisance functions. Second, we evaluate their performance in a high-dimensional, sparse setting. In both experiments, treatment and outcome rates are each 50%. Throughout, we use \(\sigma(\cdot)\) for the expit function. Additional implementation details are available in Appendix 9.
We create our first dataset according to the following model, \[\begin{gather} \text{ X \sim Uni(-1, 1)} \\ \text{ Y \mid X \sim Ber(.45X^3 + .5)} \\ \text{ A \mid X, Y \sim Ber(.85Y(1 - X^{\frac{2}{5}}) + .1(1-Y)X + .12)}, \end{gather}\] which generates a smooth CATE. The smoothness of the resulting CATE and nuisance functions is favorable to our proposed approach. We design two different missingness mechanisms corresponding to the MAR and MCCAR settings, \[\begin{gather} \text{ R \mid X, Y \sim Ber(C\sigma(2YX + 2(1-X)(1-Y)) + \delta)} \\ \text{ R \mid X \sim Ber(C\sigma(X + 1.2 + \delta))} \end{gather}\] where we set \(\delta=0.01\) to ensure positivity and choose \(C\) to match the desired treatment observation rate. To test convergence properties from Propositions 5 and 6, we add controlled noise, \(\epsilon \sim \mathcal{N}(-n^{-\alpha}, n^{-2\alpha})\), at varying rates \(\alpha\) to all nuisance functions. This approach allows us to precisely control the rate of convergence to evaluate how nuisance function estimation error affects CATE performance. We compare both estimators to an oracle, using noise-free nuisance functions and estimating the CATE with smoothing splines.
Our findings confirm our theoretical predictions. As Fig. 3 reveals, \(\widehat \tau_{\textrm{MCCAR}}\) displays persistent bias in the MAR setting, confirming it is misspecified. Both estimators converge to the oracle in the MCCAR setting, as expected. At very slow convergence rates, corresponding to poor nuisance function estimation, \(\widehat \tau_{\textrm{MAR}}\) has higher error than \(\widehat \tau_{\textrm{MCCAR}}\), but converges more quickly to the oracle. While the observed treatment proportion does not affect asymptotic convergence rates, it does affect finite-sample performance through the asymptotic variance in Propositions 5 and 6. Intuitively, more missing data inflates the inverse-probability weights and thus the variance, and so the oracle RMSE is lower at higher observation rates. The MAR estimator partially mitigates this through the imputation model \(\lambda_d\).
Our second experiment evaluates performance in a sparse setting with different CATE structure. In addition to oracle models, we compare our estimators to outcome regressions based on Eqs. 5 and 7 and propensity score models appropriate for each setting. We consider a setting with \(d = 50\) covariates, where \(\mathbf{X}\sim N(0, I_d)\) and \[\begin{gather} \text{ Y \mid \mathbf{X}\sim Ber(.8\sigma(3X_1 - \mathbb{1}(X_4 < 0) + .75\mathbb{1}(X_3 >.5)) + .1)} \\ \text{ A \mid \mathbf{X}, Y \sim Ber(.7Y\mathbb{1}(X_3 < .5) + .2\mathbb{1}(-1.2 < X_{4}) + .1)}, \end{gather}\] with missingness mechanisms for the MAR and MCCAR settings respectively: \[\begin{gather} \text{ R \mid \mathbf{X}, Y \sim Ber(0.1\mathbb{1}(X_1 > 0) + 0.7 Y + 0.1)} \\ \text{ R \mid \mathbf{X}\sim Ber(0.1\mathbb{1}(X_1 > 0) + 0.7\mathbb{1}(X_3 < 0) + 0.1),} \end{gather}\] and so roughly half the units are observed, half are treated, and half have \(Y=1\). We vary \(N \in \{1000, 2000, 5000, 10,000\}\) over 50 iterations using Random Forests for nuisance function and CATE estimation [34].
Our results confirm that in the MAR setting, the MCCAR estimator can be extremely biased. In Fig. 4, \(\widehat \tau_{\textrm{MCCAR}}\) RMSE remains large even as \(N\) increases, highlighting that model misspecification cannot be overcome simply by collecting more data. In contrast, both estimators attain similar error in the MCCAR setting, with neither estimator clearly outperforming the other at large sample sizes.
Lastly, we test the robustness implications of Propositions 6 and 5 by systematically misspecifying nuisance functions using Random Forests with limited depth. Fig. 5 validates the robustness properties of our MAR estimator. The estimator remains consistent when either \(\lambda_d\) or \(\pi\) is misspecified, but not when both are misspecified simultaneously or when \(\gamma_d\) is misspecified. Similarly, the MCCAR estimator is robust to misspecification of \(\eta_d\) or \(\nu_d\), but not both together. Complete results for the MCCAR estimator are in Appendix 9.
To evaluate our estimators under known missingness mechanisms while preserving realistic covariate and outcome structure, we construct semi-synthetic datasets by generating synthetic missingness on treatment assignment using complete data from randomized controlled trials (RCTs) in policy-relevant domains. This approach allows us to assess policy performance across varying missingness levels, which is infeasible with naturally occurring missingness where ground truth is unknown, using true covariates and outcomes, and has previously been used in prior works studying treatment missingness ([18]). We provide complete preprocessing, missingness generation, and model details in Appendix 10.
We use two datasets varying in size, covariates, and outcome type. Voting: A large-scale field experiment of 230,000 individuals studying social pressure effects on voter turnout [35]. We estimate the CATE of a voting history mailer on binary voting outcomes. STAR: A class-size experiment with 7,000 students randomly assigned to small, regular, or regular-with-aide classes [36]. We estimate the CATE of small vs. regular class size on kindergarten test scores. For both datasets, we create treatment missingness under MAR and MCCAR assumptions and vary the observed treatment proportion from 10% to 90%.
We compare our proposed estimators against several baselines spanning outcome regression, propensity weighting, and deep learning approaches. We estimate the CATE with our proposed DR-MAR (Eq. 8 ) and DR-MCCAR (Eq. 9 ) estimators accounting for MAR and MCCAR missingness, respectively, OR outcome regression models based on Eq. 5 in the MAR experiment and Eq. 7 in the MCCAR experiment, PW-Learner which regresses IPW-based pseudo scores for MAR and MCCAR estimators on observed covariates (see Appendix for identification and implementation) [37], MTRNET, the adversarial learning approach of [18]. In Appendix 10.2, we compare our models to Risk models, which directly model outcomes without using treatment and are common in policy applications [19], and the standard doubly robust estimator using complete-cases only. We tune MTRNET via Bayesian hyperparameter search [38]. For all other models, we use ensembles of tuned linear models, k-Nearest Neighbors, Random Forests, and XGBoost [34], [39].
We exploit the random assignment from the original experiments. While we introduce missingness during model training, we evaluate on the complete dataset with true treatment assignments. We rank individuals by their estimated treatment benefit and assign treatment to the top \(\kappa-\) fraction under budget constraint \(\kappa\). We calculate policy values on a held-out test set, \(V(\hat{d}_{\kappa}) = \mathbb{E} \left[Y\mathbb{1}(A = \hat{d}_{\kappa}(\mathbf{X}))\right]\), under decision \(d\) and budget constraint \(\kappa\). To summarize performance across budgets, we use the Area-Under-the-Prescriptive-Effect-Curve (AUPEC) [40], which measures average policy value relative to random assignment. A higher AUPEC score indicates a policy has achieved higher utility, and any score over zero denotes gains over a random policy. We evaluate all models across three sources of randomness: treatment missingness, model seeds, and train-test splits, fixing two and varying the third over 10 seeds (STAR) and 5 seeds (Voting). We report results varying treatment missingness here and defer other results to Appendix 10.2.
MAR-based estimators demonstrate robust performance across varying treatment missingness levels and budget constraints in both the MCCAR and MAR settings, particularly when treatment data is severely limited. Even with only 10-30% of treatment data observed, the MAR-based OR and DR-MAR estimators maintain positive AUPEC clearly outperforming the random assignment baseline. On the Voting dataset, the OR outperforms the DR-MAR estimator even at low levels of treatment observation. The Voting dataset consists predominantly of binary variables that are strongly predictive of outcomes, and so the outcome and treatment assignments models may be estimated well in this easier predictive setting, diminishing the relative advantage of the DR-MAR estimator. On the other hand, the STAR dataset is smaller and more heterogeneous, making the treatment assignment model harder to estimate reliably. Since the OR estimator in the MAR setting depends on this model, estimation error here can degrade its performance relative to DR-MAR. Additionally, in the MAR setting, the DR-MCCAR estimator achieves similar performance to DR-MAR at high levels of treatment observation. While the DR-MCCAR estimator may be biased under MAR, the resulting policies can remain effective provided the sign and relative ordering of treatment effects are preserved. In the MCCAR setting, most estimators perform well, with the exception of MTRNET on the STAR dataset. We posit that this result is in part due to the smaller (\(<10k\) observations), more heterogeneous STAR dataset, which is representative of the type of setting where tree-based models typically outperform neural networks [41]. In general, the IPW-based estimators do not perform well on our datasets under our generated missingness mechanisms. Overall, our results confirm both the importance of flexible estimation methods and the value of properly modeling the missingness mechanism rather than discarding incomplete observations.
In many settings with missing data, complete-case analysis is the default practice, and yet we show it is not the efficient choice even when complete-case assumptions hold. We prove that MAR-based estimators achieve lower asymptotic variance under MCCAR assumptions. Our experiments demonstrate that estimators constructed using MAR assumptions substantially outperform alternatives when missingness depends on outcomes, while also remaining competitive under MCCAR assumptions.
Several limitations and directions for future work merit discussion. First, the MAR assumption may be violated in practice if unobserved factors simultaneously influence treatment assignment and the likelihood of treatment observation. This is the MNAR setting. Second, while our R-Positivity assumption is necessary for identification and estimation in both the MAR and MCCAR settings, it may be violated when certain subpopulations systematically lack treatment records. Third, the efficiency gain of the MAR estimator is an asymptotic result. In particular, it requires the nuisance estimators to satisfy the convergence rate conditions of Proposition 6, which may not hold with very small sample sizes. Finally, our results depend on accurate CATE estimation for the entire population, and so if treatment budgets are known in advance, methods that optimize directly for accurate ranking near the treatment threshold may prove more effective. Despite these limitations, the choice between MAR and MCCAR estimators is asymmetric. Using MCCAR estimators under MAR assumptions causes persistent bias, while using MAR estimators under MCCAR assumptions can achieve higher efficiency. This asymmetry reinforces the practical case for defaulting to MAR-based estimation when the missingness mechanism is uncertain.
Overall, our results have important implications for personalized decision-making across domains with missing data. We prove that MAR-based estimators are more efficient and require fewer identifying assumptions than MCCAR estimators, which are misspecified when missingness is also dependent on outcomes. Our efficient, robust, and nonparametric estimators make this recommendation practical, requiring only standard machine learning tools. Our work enables unbiased CATE estimation in settings where complete-case approaches fail, extending personalized recommendations to the many settings where treatment data is missing.
EK was supported by NSF CAREER Award 2047444.
In this section, we identify the causal quantity of interest by linking our estimand to the observed data. We present the identification results under both our MAR and MCCAR assumptions on the missingness mechanism.
Under the Standard Assumptions 1 and MAR Assumptions 3, [17] identified the ATE under treatment \(A\), which we can easily use to identify the value under policy \(d\): \[\begin{align} \mathbb{E}\left[Y(d(\mathbf{X}))\right]&= \int_{\mathcal{X}} \mathbb{E}\left[Y| \mathbf{X}= \mathbf{x}, A = d(\mathbf{X}) \right] d \mathbb{P}(\mathbf{x}) \\ &= \int_{\mathcal{X}} \int_{\mathcal{Y}} \frac{y \mathbb{P}(A = d(\mathbf{X}) \mid \mathbf{X}= \mathbf{x}, Y = y)}{\mathbb{P}(A = d(\mathbf{X}) \mid \mathbf{X}= \mathbf{x})}d\mathbb{P}(Y \leq y \mid \mathbf{X}=\mathbf{x})d\mathbb{P}(\mathbf{x}) \\ &= \int_{\mathcal{X}} \int_{\mathcal{Y}} \frac{y \lambda_{d}(\mathbf{x}, y)}{\int_{\mathcal{Y}} \lambda_{d}(\mathbf{x}, y) d\mathbb{P}(Y \leq y \mid \mathbf{X}=\mathbf{x})} d\mathbb{P}(Y \leq y \mid \mathbf{X}= \mathbf{x})d\mathbb{P}(\mathbf{x}) \\&=\mathbb{E}\left[ \frac{\beta_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})} \right]= \psi_{\text{MAR}}, \end{align}\] where the second equality follows by Bayes’ rule and third by our exchangeability assumptions. The identification for the CATE follows by the same argument and dropping the outer expectation: \[\label{eq:MarCateID} \tau(\mathbf{X}) = \mathbb{E} \left[ Y(1) - Y(0) \mid \mathbf{X}\right] = \frac{\beta_{1}(\mathbf{X})}{\gamma_1{(\mathbf{X})}} - \frac{\beta_0(\mathbf{X})}{\gamma_{0}(\mathbf{X})}.\tag{10}\] Notably, both \(\beta_{d}\) and \(\gamma_{d}\) can be estimated using the full sample, including observations with missing treatments.
Using the Standard Assumptions 1 and MCCAR Assumptions 2, we identify the value under policy \(d\): \[\begin{align} \mathbb{E}\left[Y(d(\mathbf{X}))\right]&= \int_{\mathcal{X}} \mathbb{E}\left[Y| \mathbf{X}= \mathbf{x}, A = d(\mathbf{X})\right] d\mathbb{P}(\mathbf{x}) \\ &= \int_{\mathcal{X}} \mathbb{E}\left[Y| \mathbf{X}= \mathbf{x}, A = d(\mathbf{X}), R = 1\right] d \mathbb{P}(\mathbf{x}) \\ &= \mathbb{E}\left [\nu_{d}(\mathbf{X}) \right]= \Phi_{\text{MCCAR}}. \end{align}\] And so, the CATE is identified as: \[\label{eq:MccarCateId} \tau(\mathbf{X}) = \mathbb{E} \left[ Y(1) - Y(0) \mid \mathbf{X}\right] = \nu_1(\mathbf{X}) - \nu_0(\mathbf{X}).\tag{11}\] We remark that under Assumption \(\ref{assum:mccar}\) the value of policy \(d\) is still identified by \(\mathbb{E}\left[ \frac{\beta_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})} \right]\). This result follows by the same proof as the MAR setting and noting that the MCCAR R-ignorability assumption, \(R \mathpalette{\independenT}{\perp}(A, Y) \mid \mathbf{X}\) implies the weaker MAR R-ignorability assumption, \(R \mathpalette{\independenT}{\perp}A \mid (\mathbf{X}, Y)\).
In this section, we prove Propositions 5 and 6. We then use these results to prove Theorem 7. Note that while Section 3 presents the uncentered IFs \(\phi_{\text{MCCAR}}\) and \(\varphi_{\text{MAR}}\) used in constructing our estimators, in the proofs in this section we mainly use the centered IFs: \[\begin{align} \mathbb{IF}_{\textrm{MCCAR}}(z, \widehat{\mathbb{P}}) &= \phi_{\text{MCCAR}}(z, \widehat{\mathbb{P}}) - \Phi_{\text{MCCAR}}\\ \mathbb{IF}_{\textrm{MAR}}(z, \widehat{\mathbb{P}}) &= \varphi_{\text{MAR}}(z, \widehat{\mathbb{P}}) - \psi_{\text{MAR}}. \end{align}\] For each proposition, we verify our candidate IFs by the von Mises Expansion (also known as the distributional Taylor Expansion): \[\label{eq:vmexpansion} \psi(P) - \psi(\overline{P}) = \int \mathbb{IF}(z, \overline{P})d(\bar P - P)(z) + R_2(\overline{P}, P),\tag{12}\] where \(\psi\) is our estimand of interest for population \(P\), \(\overline{P}\) is a perturbed distribution on \(P\), influence function \(\mathbb{IF}\) is a mean-zero function with finite variance, and \(R_2(\overline{P}, P)\) is a second-order error term. As a result of this second-order remainder, functions \(\mathbb{IF}\) which satisfy this expansion achieve faster convergence to the true estimand \(\psi\) and under nonparametric statistical models also achieve the semiparametric efficiency bound [26]. Using the results of the von Mises Expansion, we can construct doubly robust estimators and analyze the asymptotic properties of these estimators:
Lemma 10. Let \(\mathbb{IF}(z, \widehat{\mathbb{P}}) = \varphi(z, \widehat{\mathbb{P}}) - \mathbb{P}_n(\psi(\widehat{\mathbb{P}}))\) satisfy Eq. 12 for estimand \(\psi\) under estimated distribution \(\widehat{\mathbb{P}}\). Then, we may define the bias-corrected estimator as: \[\begin{align} \widehat \psi = \mathbb{P}_n(\mathbb{IF}(Z, \widehat{\mathbb{P}})) + \psi(\widehat{\mathbb{P}}) = \mathbb{P}_n(\varphi(Z, \widehat{\mathbb{P}})). \end{align}\] Assume that all nuisance functions are estimated from an independent sample and that \(\|\mathbb{IF}(z, \widehat{\mathbb{P}}) - \mathbb{IF}(z, \mathbb{P})\| = o_{\mathbb{P}}(1)\). Then, by Proposition 1 from [26], we have: \[\begin{align} \widehat \psi - \psi &= \mathbb{P}_n(\mathbb{IF}(Z, \widehat{\mathbb{P}})) + \psi(\widehat{\mathbb{P}}) - \psi(\mathbb{P}) \\ &= \mathbb{P}_n(\mathbb{IF}(Z, \mathbb{P})) - \mathbb{E} \left[\mathbb{IF}(Z; \mathbb{P})\right] + R_2(\widehat{\mathbb{P}}, \mathbb{P}) +o_{\mathbb{P}}\left(1/\sqrt{n}\right). \end{align}\] The first two terms are the difference of a sample average and the true population mean. If \(R_2(\widehat{\mathbb{P}}, \mathbb{P}) = o_{\mathbb{P}}\left(1 /\sqrt n\right)\), then by the Central Limit Theorem and Slutsky’s Theorem, \[\begin{align} \sqrt n(\widehat \psi - \psi) \rightsquigarrow N\left(0, \textrm{Var}(\mathbb{IF}(Z; \mathbb{P}))\right). \end{align}\]
Therefore, to establish \(\sqrt{n}\)-consistency and asymptotic normality under Lemma 10, it suffices to provide conditions under which \(R_2(\widehat{\mathbb{P}}, \mathbb{P}) = o_{\mathbb{P}}(1/\sqrt{n})\) for an IF satisfying the von Mises Expansion (Eq. 12 ).
Recall our bias-corrected estimator \(\widehat{\Phi}_{\text{MCCAR}}= \mathbb{P}_n(\mathbb{IF}_{\textrm{MCCAR}}(Z, \widehat{\mathbb{P}})) + \Phi_{\text{MCCAR}}(\widehat{\mathbb{P}}) = \mathbb{P}_n(\phi_{\text{MCCAR}}(Z, \widehat{\mathbb{P}}))\). We first verify that \(\mathbb{IF}_{\textrm{MCCAR}}\) satisfies the von Mises Expansion (Eq. 12 ) and derive the conditions under which the remainder term satisfies \(R_2(\widehat{\mathbb{P}}, \mathbb{P}) = o_{\mathbb{P}}(1/\sqrt{n})\). We conclude by deriving the asymptotic variance.
Recall our candidate influence function for the value function in the MCCAR setting: \[\begin{align} \label{eq:MccarIf} \mathbb{IF}_{\textrm{MCCAR}}= \frac{R\mathbb{1}(A = d(\mathbf{X}))(Y - \nu_{d}(\mathbf{X}))}{\eta_{d}(\mathbf{X})} + \nu_{d}(\mathbf{X}) - \Phi_{\text{MCCAR}} . \end{align}\tag{13}\] In this section, we will prove that the candidate IF (Eq. 13 ) admits a von Mises Expansion (Eq. 12 ) with second-order remainder bounded by \[R_2(\widehat{\mathbb{P}}, \mathbb{P}) \lesssim \|\eta_{d}- \widehat \eta_{d}\|\|\nu_{d}- \widehat \nu_{d}\|,\] where we take the estimated distribution \(\widehat{\mathbb{P}}\) to be the perturbed distribution \(\overline{P}\). We begin by showing that \(\mathbb{IF}_{\textrm{MCCAR}}\) (Eq. 13 ) is mean zero: \[\begin{align} \mathbb{E}\left[\mathbb{IF}_{\textrm{MCCAR}}\right] &= \mathbb{E}\left[\frac{R\mathbb{1}(A = d(\mathbf{X}))(Y - \nu_{d}(\mathbf{X}))}{\eta_{d}(\mathbf{X})} + \nu_{d}(X)\right] - \mathbb{E}\left[\nu_{d}(\mathbf{X})\right] \\ &= \mathbb{E}\left[\frac{\eta_{d}(\mathbf{X})}{\eta_{d}(\mathbf{X})}\biggr( \mathbb{E}\left[Y\mid \mathbf{X}= \mathbf{x}, A = d(\mathbf{X}), R =1 \right] - \nu_{d}(\mathbf{X}) \biggr)\right] = 0, \end{align}\] where the second equality follows by iterated expectation over \(\mathbf{X}\) and the fact that \(R, \mathbb{1}(A = d(\mathbf{X}))\) are binary. Using this result, we may arrange the von Mises Expansion (Eq. 12 ) in terms of the remainder \[\label{app:eqmccarremainder} R_2(\widehat{\mathbb{P}}, \mathbb{P}) = \Phi_{\text{MCCAR}}(\widehat{\mathbb{P}}) - \Phi_{\text{MCCAR}}(\mathbb{P}) + \int \mathbb{IF}_{\textrm{MCCAR}}(z; \widehat{\mathbb{P}})d \mathbb{P}(z).\tag{14}\] Expanding the final term of the remainder (14 ), \[\begin{align} \int \mathbb{IF}_{\textrm{MCCAR}}(z; \widehat{\mathbb{P}})d\mathbb{P}(z) &= \int \biggr\{ \frac{R \mathbb{1}(A = d(\mathbf{x}))(Y - \widehat \nu_{d}(\mathbf{x}))}{\widehat \eta_{d}(\mathbf{x})} + \widehat \nu_{d}(\mathbf{x}) - \Phi_{\text{MCCAR}}(\widehat{\mathbb{P}})\biggr\}d\mathbb{P}(z) \\ &= \int \biggr\{\frac{\eta_{d}(\mathbf{x})}{\widehat \eta_{d}(\mathbf{x})} \cdot (\nu_{d}(\mathbf{x}) - \widehat \nu_{d}(\mathbf{x})) + \widehat \nu_{d}(\mathbf{x})\biggr\}d\mathbb{P}(z) - \Phi_{\text{MCCAR}}(\widehat{\mathbb{P}}) \\ &= \int \biggr\{\frac{\eta_{d}(\mathbf{x})}{\widehat \eta_{d}(\mathbf{x})} \cdot (\nu_{d}(\mathbf{x}) - \widehat \nu_{d}(\mathbf{x})) - (\nu_{d}(\mathbf{x}) - \widehat \nu_{d}(\mathbf{x}))\biggr\}d\mathbb{P}(z) + \Phi_{\text{MCCAR}}(\mathbb{P}) - \Phi_{\text{MCCAR}}(\widehat{\mathbb{P}}) \\ &= \int \left(\frac{1}{\widehat \eta_{d}(\mathbf{x})} - \frac{1}{\eta_{d}(\mathbf{x})}\right)(\nu_{d}(\mathbf{x}) - \widehat \nu_{d}(\mathbf{x}))(\eta_{d}(\mathbf{x}))d\mathbb{P}(z) + \Phi_{\text{MCCAR}}(\mathbb{P}) - \Phi_{\text{MCCAR}}(\widehat{\mathbb{P}}) \end{align}\] where the second equality follows by iterated expectation over \(\mathbf{X}\) and the fact that \(R, \mathbb{1}(A = d(\mathbf{X}))\) are binary. And thus, the remainder (14 ) is equivalently \[R_2(\widehat{\mathbb{P}}, \mathbb{P}) = \int \left(\frac{1}{\widehat \eta_{d}(\mathbf{x})} - \frac{1}{\eta_{d}(\mathbf{x})}\right)(\nu_{d}(\mathbf{x}) - \widehat \nu_{d}(\mathbf{x}))(\eta_{d}(\mathbf{x}))d\mathbb{P}(\mathbf{x}).\] Under \(A-\) and \(R-\)positivity (Assumptions 1 and 2), we have \[\begin{align} |R_2(\widehat{\mathbb{P}}, \mathbb{P})| \leq \frac{1}{\epsilon} \int | \eta_{d}(\mathbf{x}) - \widehat \eta_{d}(\mathbf{x})| |\nu_{d}(\mathbf{x}) - \widehat \nu_{d}(\mathbf{x})|d\mathbb{P}(\mathbf{x}) \lesssim \| \eta_{d}(\mathbf{x}) - \widehat \eta_{d}(\mathbf{x}) \| \| \nu_{d}(\mathbf{x}) - \widehat \nu_{d}(\mathbf{x})\|, \end{align}\] where the final inequality follows by Cauchy-Schwarz. As a result, the remainder term is second-order and we conclude that the von Mises Expansion (Eq. 12 ) holds.
And so, Lemma 10 applies when \(\| \eta_{d}(\mathbf{x}) - \widehat \eta_{d}(\mathbf{x}) \| \| \nu_{d}(\mathbf{x}) - \widehat \nu_{d}(\mathbf{x})\| = o_{\mathbb{P}}(1/\sqrt n)\), establishing the convergence properties of Proposition 5.
To complete the proof of Proposition 5, we derive the asymptotic variance.
Using the fact that \(\mathbb{E}\left[\mathbb{IF}_{\textrm{MCCAR}}\right] = 0\) (shown above), \[\begin{align} \label{app:eq:varMCCARlong} \textrm{Var}(\phi_{\text{MCCAR}}) &= \mathbb{E}\left[\left( \frac{R\mathbb{1}(A = d(\mathbf{X}))(Y - \nu_{d}(\mathbf{X}))}{\eta_{d}(\mathbf{X})} + \nu_{d}(\mathbf{X}) - \Phi_{\text{MCCAR}} \right)^{2}\right] \nonumber \\&= \mathbb{E}\left[\left(\frac{R\mathbb{1}(A = d(\mathbf{X}))(Y-\nu_{d}(\mathbf{X}))}{\eta_{d}(\mathbf{X})}\right)^{2}\right] + \mathbb{E}\left[\left( \nu_{d}(\mathbf{X}) - \Phi_{\text{MCCAR}}\right)^{2}\right] \nonumber \\&\quad + \mathbb{E}\left[ \left(\frac{R\mathbb{1}(A = d(\mathbf{X}))(Y-\nu_{d}(\mathbf{X}))}{\eta_{d}(\mathbf{X})}\right) \cdot \left( \nu_{d}(\mathbf{X}) - \Phi_{\text{MCCAR}}\right)\right]. \end{align}\tag{15}\] We begin by showing that the final term of Eq. 15 is zero. By iterated expectation, \[\begin{align} &\mathbb{E}\left[\mathbb{E}\left[\left(\frac{R\mathbb{1}(A = d(\mathbf{X}))(Y - \nu_{d}(\mathbf{X}))}{\eta_{d}(\mathbf{X})} \right) \left(\nu_{d}(\mathbf{X}) - \Phi_{\text{MCCAR}}\right) \mid \mathbf{X}\right]\right] \\ &= \mathbb{E}\left[\frac{\eta_{d}(\mathbf{X})}{\eta_{d}(\mathbf{X})} \left(\mathbb{E}\left[Y|\mathbf{X}, R = 1, A = d(\mathbf{X}) \right] - \nu_{d}(\mathbf{X})\right)(\nu_{d}(\mathbf{X}) - \Phi_{\text{MCCAR}})\right] = 0, \end{align}\] which holds because \(R, \mathbb{1}(A = d(\mathbf{X}))\) are binary. Then, rewriting the first term of the Eq. 15 , \[\begin{align} \label{app:eq:varstep} \mathbb{E}\left[\left(\frac{R\mathbb{1}(A = d(\mathbf{X}))(Y-\nu_{d}(\mathbf{X}))}{\eta_{d}(\mathbf{X})}\right)^{2}\right] &= \mathbb{E}\left[\frac{R\mathbb{1}(A = d(\mathbf{X}))(Y - \nu_{d}(\mathbf{X}))^{2}}{\eta_{d}(\mathbf{X})^{2}}\right] \equiv \mathbb{E}\left[\frac{\sigma_{d}^{2}(\mathbf{X})}{\eta_{d}(\mathbf{X})}\right] \end{align}\tag{16}\] where the second equality follows by iterated expectation. It follows that \[\label{app:eq:varMCCAR} \textrm{Var}(\phi_{\text{MCCAR}}) = \mathbb{E}\left[\frac{\sigma_{d}^{2}(\mathbf{X})}{\eta_{d}(\mathbf{X})}\right] + \mathbb{E}\left[\left( \nu_{d}(\mathbf{X}) - \Phi_{\text{MCCAR}}\right)^{2}\right] .\tag{17}\]
For completeness, we include our derivation of the candidate IF. We follow the "tricks" from [26], namely (1) treating the data as discrete, (2) treating the IF as a derivative, and (3) utilizing known IFs, such as the IF of a conditional expectation or density. First, following [26], we define \(\text{InfFun} : \psi \to L_2(\mathbb{P})\) as an operator that maps a function \(\psi\) to its influence function under a nonparametric model. Then, \[\begin{align} \text{InfFun}(\Phi_{\text{MCCAR}}) &= \text{InfFun}\left(\sum_{x} \nu_d(x)p(x)\right) \\ &= \sum_{x}\big(\text{InfFun}(\nu_d(\mathbf{x}))p(\mathbf{x}) + \nu_d(\mathbf{x})\text{InfFun}(p(\mathbf{x}))\big) \\ &= \sum_{x}\left( \frac{\mathbb{1}(\mathbf{X}= \mathbf{x}, R = 1, A = d(\mathbf{x}))}{p(\mathbf{X}= \mathbf{x}, R =1 , A = d(\mathbf{x}))}\left(Y - \nu_d(\mathbf{x})\right)p(\mathbf{x}) + \nu_d(\mathbf{x})(\mathbb{1}(\mathbf{X}= \mathbf{x}) - p(\mathbf{x})) \right) \\ &= \frac{R\mathbb{1}(A = d(\mathbf{X}))}{\eta_d(\mathbf{X})}(Y - \eta_d(\mathbf{x})) + \nu_d(\mathbf{X}) - \Phi_{\text{MCCAR}}, \end{align}\] which is exactly the centered influence function in (13 ).
Recall our bias-corrected estimator \(\widehat{\psi}_{\text{MAR}}= \mathbb{P}_n(\mathbb{IF}_{\textrm{MAR}}(Z, \widehat{\mathbb{P}})) + \psi_{\text{MAR}}(\widehat{\mathbb{P}}) = \mathbb{P}_n(\varphi_{\text{MAR}}(Z, \widehat{\mathbb{P}}))\). We now verify that \(\varphi_{\text{MAR}}\) satisfies the von Mises Expansion (Eq. 12 ) and provide sufficient conditions on the nuisance functions for \(\sqrt{n}\)-consistency and asymptotic normality by analyzing the remainder term and applying Lemma 10. We then derive the asymptotic variance.
[17] verified the influence function for the average effect of treatment \(a\) under MAR assumptions. We first note that the effect under policy \(d\) can be seen as a generalization as the average effect of treatment \(a\), as this can be represented by a policy that gives all individuals treatment \(a\). Therefore, using this fact, the remainder term Lemma 1 of [17] may be written as: \[\begin{align} R_2(\widehat{\mathbb{P}}, \mathbb{P}) &= \int \left(\frac{y - \frac{\widehat \beta_{d}(\mathbf{x})}{\widehat \gamma_{d}(\mathbf{x})}}{\widehat \gamma_{d}(\mathbf{x})}\right)\left(\frac{\pi(\mathbf{x}, y) - \widehat \pi(\mathbf{x}, y)}{\widehat \pi(\mathbf{x}, y)}\right)\left(\lambda_{d}(\mathbf{x}, y) - \widehat \lambda_{d}(\mathbf{x}, y)\right) \\ &\qquad + \left(\frac{\beta_{d}(\mathbf{x})}{\gamma_{d}(\mathbf{x})} - \frac{\widehat \beta_{d}(\mathbf{x})}{\widehat \gamma_{d}(\mathbf{x})}\right)\left(\frac{\gamma_{d}(\mathbf{x}) - \widehat \gamma_{d}(\mathbf{x})}{\widehat \gamma_{d}(\mathbf{x})}\right)d\mathbb{P}(z) \end{align}\]
And so, the remainder term is second-order and we conclude that the von Mises Expansion (Eq. 12 ) holds. As [17] noted, the conditions of Lemma 10 are met when \(\|\pi- \widehat \pi\|\|\lambda_{d}-\widehat \lambda_{d}\| + \|\beta_{d}- \widehat \beta_{d}\|\|\gamma_{d}-\widehat \gamma_{d}\| + \|\gamma_{d}- \widehat \gamma_{d}\|^{2} = o_{\mathbb{P}}(1/\sqrt n)\).
To complete the proof of Proposition 6, we derive the asymptotic variance.
Because \(\mathbb{E}\left[\mathbb{IF}_{\textrm{MAR}}\right] = 0\), \(\textrm{Var}(\varphi_{\text{MAR}}) =\) \[\begin{align} \label{app:eq:MARvarexpanded} &\mathbb{E}\left[\left( \left(\frac{Y - \beta_{d}(\mathbf{X}) / \gamma_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})}\right) \cdot \left(\frac{R \left(\mathbb{1}(A = d(\mathbf{X})) - \lambda_{d}(\mathbf{X},Y) \right)}{\pi(\mathbf{X},Y)} + \lambda_{d}(\mathbf{X},Y) \right) + \frac{\beta_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})} - \psi_{\text{MAR}} \right)^{2}\right] \nonumber \\ &= \mathbb{E}\left[\biggr\{\left(\frac{Y - \beta_{d}(\mathbf{X}) / \gamma_{d}(\mathbf{X}) }{\gamma_{d}(\mathbf{X})}\right) \left(\frac{R \left(\mathbb{1}(A = d(\mathbf{X})) - \lambda_{d}(\mathbf{X},Y) \right)}{\pi(\mathbf{X},Y)} + \lambda_{d}(\mathbf{X},Y) \right)\biggr\}^{2}\right] \nonumber \\&\quad+ \mathbb{E}\left[\left(\frac{\beta_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})} - \psi_{\text{MAR}}\right)^{2}\right] \nonumber \\&\quad + \mathbb{E}\left[\left(\frac{Y - \beta_{d}(\mathbf{X}) / \gamma_{d}(\mathbf{X}) }{\gamma_{d}(\mathbf{X})}\right) \left(\frac{R \left(\mathbb{1}(A = d(\mathbf{X})) - \lambda_{d}(\mathbf{X},Y) \right)}{\pi(\mathbf{X},Y)} + \lambda_{d}(\mathbf{X},Y) \right)\left(\frac{\beta_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})} - \psi_{\text{MAR}}\right)\right]. \end{align}\tag{18}\] We first show that the cross-product term in Eq. 18 is zero. By iterated expectation over \(\mathbf{X}\) and \(Y\), \[\begin{align} &\mathbb{E}\left[ \mathbb{E}\left[\left(\frac{Y - \beta_{d}(\mathbf{X}) / \gamma_{d}(\mathbf{X}) }{\gamma_{d}(\mathbf{X})}\right) \left(\frac{R \left(\mathbb{1}(A = d(\mathbf{X})) - \lambda_{d}(\mathbf{X},Y) \right)}{\pi(\mathbf{X},Y)} + \lambda_{d}(\mathbf{X},Y) \right)\left(\frac{\beta_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})} - \psi_{\text{MAR}}\right) \mid \mathbf{X}, Y \right] \right] \\ &\quad= \mathbb{E}\left[\left(\frac{\beta_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})} - \psi_{\text{MAR}}\right) \left(\frac{Y - \beta_{d}(\mathbf{X}) / \gamma_{d}(\mathbf{X}) }{\gamma_{d}(\mathbf{X})}\right) \left(\frac{\pi(\mathbf{X}, Y) \left(\lambda_{d}(\mathbf{X}, Y) - \lambda_{d}(\mathbf{X},Y) \right)}{\pi(\mathbf{X},Y)} + \lambda_{d}(\mathbf{X},Y) \right)\right] \\ &\quad= \mathbb{E}\left[\left(\frac{\beta_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})} - \psi_{\text{MAR}}\right) \left(\frac{\mathbb{E}[Y\lambda_{d}(\mathbf{X},Y) \mid \mathbf{X}] - \beta_{d}(\mathbf{X}) / \gamma_{d}(\mathbf{X}) \mathbb{E}[\lambda_{d}(\mathbf{X},Y)\mid \mathbf{X}] }{\gamma_{d}(\mathbf{X})}\right)\right] = 0. \end{align}\] Rearranging the first-term in Eq. 18 , \[\begin{align} &\mathbb{E}\left[\biggr\{\left(\frac{Y - \beta_{d}(\mathbf{X}) / \gamma_{d}(\mathbf{X}) }{\gamma_{d}(\mathbf{X})}\right) \cdot \left(\frac{R \left(\mathbb{1}(A = d(\mathbf{X})) - \lambda_{d}(\mathbf{X},Y) \right)}{\pi(\mathbf{X},Y)} + \lambda_{d}(\mathbf{X},Y) \right)\biggr\}^{2}\right] \\&=\mathbb{E}\left[\left(\frac{Y - \beta_{d}(\mathbf{X}) / \gamma_{d}(\mathbf{X}) }{\gamma_{d}(\mathbf{X})}\right)^{2} \cdot \biggr\{\left(\frac{R \left(\mathbb{1}(A = d(\mathbf{X})) - \lambda_{d}(\mathbf{X},Y) \right)}{\pi(\mathbf{X},Y)}\right)^{2} + \lambda_{d}^{2} (\mathbf{X},Y)\biggr\}\right] \\&=\mathbb{E}\left[\left(\frac{Y - \beta_{d}(\mathbf{X}) / \gamma_{d}(\mathbf{X}) }{\gamma_{d}(\mathbf{X})}\right)^{2} \cdot \biggr\{\left(\frac{R \left(\mathbb{1}(A = d(\mathbf{X})) - 2\mathbb{1}(A = d(\mathbf{X}))\lambda_{d}(\mathbf{X},Y) + \lambda_{d}^{2}(\mathbf{X},Y) \right)}{\pi^{2}(\mathbf{X},Y)}\right) + \lambda_{d}^{2} (\mathbf{X},Y)\biggr\}\right] \\&= \mathbb{E}\left[\left(\frac{Y - \beta_{d}(\mathbf{X}) / \gamma_{d}(\mathbf{X}) }{\gamma_{d}(\mathbf{X})}\right)^{2} \cdot \left(\frac{\lambda_{d}(\mathbf{X},Y)- \lambda_{d}^{2}(\mathbf{X},Y)}{\pi(\mathbf{X},Y)} + \lambda_{d}^{2} (\mathbf{X},Y)\right)\right] \\&= \mathbb{E}\left[\left(\frac{Y - \beta_{d}(\mathbf{X}) / \gamma_{d}(\mathbf{X}) }{\gamma_{d}(\mathbf{X})}\right)^{2} \cdot \left(\frac{\lambda_{d}(\mathbf{X},Y)}{\pi(\mathbf{X},Y)} - \lambda_{d}^{2}(\mathbf{X},Y)\cdot \biggr\{\frac{1 - \pi(\mathbf{X},Y)}{\pi(\mathbf{X},Y)}\biggr\}\right)\right] \\ &= \mathbb{E}\left[\left(\frac{Y - \beta_{d}(\mathbf{X}) / \gamma_{d}(\mathbf{X}) }{\gamma_{d}(\mathbf{X})}\right)^{2} \cdot \left(\frac{\lambda_{d}(\mathbf{X},Y)}{\pi(\mathbf{X},Y)} - \lambda_{d}^{2}(\mathbf{X},Y)\cdot \biggr\{\frac{1 - \pi(\mathbf{X},Y)}{\pi(\mathbf{X},Y)}\biggr\}\right)\right] \end{align}\] where the first equality holds because the cross term is eliminated by iterated expectation over \(\mathbf{X}, Y\) (\(\mathbb{E}[\mathbb{1}(A = d(\mathbf{X})) |\mathbf{X}, Y)] - \lambda_{d}(\mathbf{X}, Y) = 0)\), and the third equality by again taking iterated expectation over \(\mathbf{X}\) and \(Y\). Therefore, we have shown \[\begin{align} \label{eq:varMAR} &\textrm{Var}(\varphi_{\text{MAR}})= \nonumber \\& \mathbb{E}\left[\left(\frac{Y - \beta_{d}(\mathbf{X}) / \gamma_{d}(\mathbf{X}) }{\gamma_{d}(\mathbf{X})}\right)^{2} \cdot \left(\frac{\lambda_{d}(\mathbf{X},Y)}{\pi(\mathbf{X},Y)} - \lambda_{d}^{2}(\mathbf{X},Y)\cdot \biggr\{\frac{1 - \pi(\mathbf{X},Y)}{\pi(\mathbf{X},Y)}\biggr\}\right)\right] + \mathbb{E}\left[\left(\frac{\beta_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})} - \psi_{\text{MAR}}\right)^{2}\right] . \end{align}\tag{19}\]
Remark 11. We note the variance of an estimand with contrasts (such as the ATE) may have an additional term from the covariance. Using the ATE as an example, the covariance between the residual terms is: \[\begin{align} &-2\mathbb{E}\biggr[\biggr\{\left(\frac{Y - \beta_1(\mathbf{X})/\gamma_1(\mathbf{X})}{\gamma_1(\mathbf{X})}\right) \left(\frac{R(\mathbb{1}(A = 1) - \lambda_1(\mathbf{X}, Y))}{\pi(\mathbf{X}, Y)} + \lambda_1(\mathbf{X}, Y)\right) \biggr\} \\& \qquad \times \left(\frac{Y - \beta_0(X)/\gamma_0(\mathbf{X})}{\gamma_0(\mathbf{X})}\right) \left(\frac{R(\mathbb{1}(A = 0) - \lambda_0(\mathbf{X}, Y))}{\pi(\mathbf{X}, Y)} + \lambda_0(\mathbf{X}, Y)\right) \biggr\}\biggr] \\ &= -2\mathbb{E}\biggr[ \biggr\{\left(\frac{Y - \beta_1(\mathbf{X})/\gamma_1(\mathbf{X})}{\gamma_1(\mathbf{X})}\right) \left(\frac{Y - \beta_0(\mathbf{X})/\gamma_0(\mathbf{X})}{\gamma_0(\mathbf{X})}\right) \biggr\} \biggr\{ \frac{-\lambda_1(\mathbf{X}, Y) \lambda_{0}(\mathbf{X}, Y)}{\pi(\mathbf{X}, Y)} + \lambda_{1}(\mathbf{X}, Y) \lambda_{0}(\mathbf{X}, Y)\biggr\} \biggr] \\&= 2\mathbb{E}\left[\biggr\{\frac{1 - \pi(\mathbf{X}, Y)}{\pi(\mathbf{X}, Y)}\biggr\} \biggr\{\frac{(Y - \beta_1(\mathbf{X})/\gamma_1(\mathbf{X}))}{\gamma_1(\mathbf{X})} \biggr\} \biggr\{\frac{(Y - \beta_0(\mathbf{X})/\gamma_0(\mathbf{X})) }{\gamma_0(\mathbf{X})} \biggr\}\biggr\{\lambda_{1}(\mathbf{X}, Y)\lambda_{0}(\mathbf{X}, Y)\biggr\}\right]. \end{align}\] Under complete data, the covariance of the residuals, i.e. \(A(Y - \nu_1)(1-A)(Y - \nu_0)\), are zero, which as we see above does not happen with \(\varphi_{\text{MAR}}\). Because of this, we may lose efficiency if the residuals have the same sign.
As \(\varphi_{\text{MAR}}\) is also a valid IF under MCCAR Assumptions \(\ref{assum:mccar}\), we conclude by deriving an expression for the asymptotic relative efficiency of estimators built using \(\varphi_{\text{MAR}}\) compared to \(\phi_{\text{MCCAR}}\).
We begin by rewriting \(\textrm{Var}(\varphi_{\text{MAR}})\) in terms of \(\textrm{Var}(\phi_{\text{MCCAR}})\). First, we note the last terms of both expressions are equivalent: \[\begin{align} \label{app:eq:outcomevar} \mathbb{E}\biggr[\biggr(\frac{\beta_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})} - \psi_{\text{MAR}}\biggr)^{2}\biggr] = \mathbb{E}\left[\left( \nu_{d}(\mathbf{X}) - \Phi_{\text{MCCAR}}\right)^{2}\right], \end{align}\tag{20}\] which follows because \(\psi_{\text{MAR}}= \mathbb{E}\biggr[\frac{\beta_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})}\biggr] = V(d) = \mathbb{E} \left[\nu_{d}(\mathbf{X})\right] = \Phi_{\text{MCCAR}}\). It is also possible to verify this equality by reversing the steps of the MAR Identification strategy (see 6.1). Now, considering the first term of \(\textrm{Var}(\varphi_{\text{MAR}})\) (Eq. 19 ), \[\begin{align} \label{app:eq:vargain} \mathbb{E}\biggr[&\biggr\{\frac{Y - \beta_{d}(\mathbf{X})/\gamma_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})}\biggr\}^{2}\cdot \biggr\{\frac{\lambda_{d}(\mathbf{X}, Y)}{\pi(\mathbf{X}, Y)} - \lambda_{d}^2(\mathbf{X},Y)\left(\frac{1-\pi(\mathbf{X},Y)}{\pi(\mathbf{X}, Y)}\right)\biggr\}\biggr] \nonumber \\ &= \mathbb{E}\biggr[\biggr\{\frac{(Y - \nu_{d}(\mathbf{X}))^{2}}{\gamma_{d}(\mathbf{X})^{2}}\biggr\}\cdot \biggr\{\frac{\mathbb{1}(A = d(\mathbf{X})}{\pi(\mathbf{X},Y)} - \lambda_{d}^2(\mathbf{X},Y)\left(\frac{1-\pi(\mathbf{X},Y)}{\pi(\mathbf{X}, Y)}\right)\biggr\}\biggr] \nonumber \\ &= \mathbb{E}\left[\frac{\sigma_{d}^{2}(\mathbf{X})}{\eta_{d}(\mathbf{X})}\right] - \mathbb{E}\left[\left(\frac{Y - \nu_{d}(\mathbf{X})}{\gamma_{d}(\mathbf{X})}\right)^{2}\cdot\left(\frac{\lambda_{d}(\mathbf{X}, Y)}{\pi(\mathbf{X})} \right) \cdot \lambda_{d}(\mathbf{X}, Y)(1 - \pi(\mathbf{X}))\right] \nonumber \\ &\equiv \mathbb{E}\left[\frac{\sigma_{d}^{2}(\mathbf{X})}{\eta_{d}(\mathbf{X})}\right] - g(d, \mathbb{P}) \end{align}\tag{21}\] where the first equality follows by iterated expectation over \(\mathbf{X}\) and \(Y\) and the fact that \(A\) is binary (\(\mathbb{E}\left[\mathbb{E}\left[\mathbb{1}(A | \mathbf{X}, Y) \mid \mathbf{X}, Y\right]\right]=A)\) as well as R-Ignorability from Assumption 2. The third equality also follows from R-Ignorability of Assumption 2. Therefore, by Eq. 20 and Eq. 21 , we may rewrite \(\textrm{Var}(\varphi_{\text{MAR}})\) as \[\begin{align} \textrm{Var}(\varphi_{\text{MAR}}) &= \mathbb{E}\left[\left( \nu_{d}(\mathbf{X}) - \Phi_{\text{MCCAR}}\right)^{2}\right] + \mathbb{E}\left[\frac{\sigma_{d}^{2}(\mathbf{X})}{\eta_{d}(\mathbf{X})}\right] - g(d, \mathbb{P}) \\&= \textrm{Var}(\phi_{\text{MCCAR}}) - g(d, \mathbb{P}). \end{align}\] Using this result, and omitting the argument for the decision policy \(d\) in \(g\), we can derive an expression for the ARE, \[\begin{align} \textrm{ARE}(\textrm{Var}(\phi_{\text{MCCAR}}), \textrm{Var}(\varphi_{\text{MAR}})) &= \frac{\textrm{Var}(\phi_{\text{MCCAR}})}{\textrm{Var}(\varphi_{\text{MAR}})} \\ &= \frac{\textrm{Var}(\phi_{\text{MCCAR}})}{\textrm{Var}(\phi_{\text{MCCAR}}) - g(\mathbb{P})} \\ &= \biggr(1 - \frac{g(\mathbb{P})}{\textrm{Var}(\phi_{\text{MCCAR}})}\biggr)^{-1}. \end{align}\] Hence, as \(g(\mathbb{P}) \geq 0\), we conclude that estimators based on \(\varphi_{\text{MAR}}\) are more efficient.
Remark 12. Recall that the asymptotic variance of contrasts of \(\varphi_{\text{MAR}}\) will have an additional covariance term (Remark 11). We now show that despite this penalty \(\varphi_{MAR}\) is still more efficient. Notice \[\begin{align} \textrm{Var}(\phi_{\text{MCCAR}}(1) - \phi_{\text{MCCAR}}(0)) &= \mathbb{E}\left[\frac{\sigma_{1}^{2}(\mathbf{X})}{\eta_{1}(\mathbf{X})}\right] + \left[\frac{\sigma_{0}^{2}(\mathbf{X})}{\eta_{0}(\mathbf{X})}\right] + \mathbb{E} \left[(\tau(\mathbf{X}) - \tau)^{2}\right] \end{align}\] And then using the argument above and covariance term from Remark 11, \[\begin{align} &\textrm{Var}(\varphi_{\text{MAR}}(1) - \varphi_{\text{MAR}}(0)) = \textrm{Var}(\phi_{\text{MCCAR}}(1) - \phi_{\text{MCCAR}}(0)) - g(1,\mathbb{P}) - g(0, \mathbb{P}) \\&\quad +2\mathbb{E}\left[\biggr\{\frac{1 - \pi(\mathbf{X})}{\pi(\mathbf{X})}\biggr\} \biggr\{\frac{(Y - \nu_1)}{\gamma_1(\mathbf{X})} \biggr\} \biggr\{\frac{(Y - \nu_0) }{\gamma_0(\mathbf{X})} \biggr\}\biggr\{\lambda_{1}(\mathbf{X}, Y)\lambda_{0}(\mathbf{X}, Y)\biggr\}\right] \end{align}\] If the covariance term (last term) is negative, which occurs when the product of the expectation of the residuals is negative, then \(\varphi_{\text{MAR}}\) is still more efficient. If the covariance term is positive, we only achieve gains if \(g(1, \mathbb{P}) + g(0, \mathbb{P})\) is greater than the covariance term. Dropping arguments for notational convenience, \[\begin{gather} g(1, \mathbb{P}) + g(0, \mathbb{P}) \geq 2\mathbb{E} \left[\left(\frac{1-\pi}{\pi}\right)\left( \frac{Y - \nu_1}{\gamma_1}\right)\left( \frac{Y - \nu_0}{\gamma_0}\right)\lambda_1\lambda0\right] \\ \mathbb{E} \left[\left(\frac{1-\pi}{\pi}\right)\biggr\{\left(\frac{(Y - \nu_1)\lambda_1}{\gamma_1}\right)^{2}-2\left(\frac{(Y - \nu_1)\lambda_1}{\gamma_1}\right)\left( \frac{(Y - \nu_0)\lambda_0}{\gamma_0}\right)+\left(\frac{(Y - \nu_0)\lambda_0}{\gamma_0}\right)^{2} \biggr\}\right] \geq 0 \end{gather}\] which clearly holds. And so, despite the covariance penalty in the ATE case we achieve efficiency gains.
We establish the asymptotic normality of the policy value under a fixed policy \(d\) in the main text. In this section, we extend these results to policy value estimation at the data-dependent optimal policy \(\widehat d^*_{\kappa}\). In order to obtain asymptotic normality in this case, we need to invoke an additional assumption:
Assumption 13 (Margin Condition). There exists constants \(C > 0\) and \(\alpha \geq 0\) such that \[\mathbb{P}(0 < |\tau(\mathbf{x}) - \tau^*_{\kappa} | \leq t) \leq Ct^\alpha \text{ for all } t > 0.\]
The margin condition, originating from [33], limits the probability mass concentrating near the optimal decision boundary \(\tau^{*}_{\kappa}\). It is satisfied whenever the CATE distribution has a continuous density near \(\tau^*_{\kappa}\), and therefore is mild in practice.
We can now use Assumption [ass:margincondition] to extend asymptotic normality to policy value estimation an estimated policy \(\widehat d^{*}_{\kappa}\) in both the MAR and MCCAR settings:
Theorem 14. Define the bias-corrected estimator derived from \(\phi_{\text{MCCAR}}\) at the estimated policy \(\widehat d_{\kappa}^{*}\) as \(\widehat{\Phi}_{\text{MCCAR}}(\widehat d_{\kappa}^{*}) = \mathbb{P}_n(\widehat{\phi}_{\text{MCCAR}}(\widehat d_{\kappa}^{*}))\). Assume the conditions of Proposition 5, in addition to Assumption [ass:margincondition] and \(\|(\widehat \tau - \widehat \tau^*_{\kappa}) - (\tau - \tau^*_{\kappa}) \|^{1 + \alpha}_{\infty} = o_{\mathbb{P}}(n^{-1/2} )\). Then, \(\sqrt{n}(\widehat{\Phi}_{\text{MCCAR}}(\widehat d_{\kappa}^{*}) - V(d^*_{\kappa})) \rightsquigarrow \mathcal{N}(0, \textrm{Var}(\phi_{\text{MCCAR}}(d^*_{\kappa})))\).
Theorem 15. Define the bias-corrected estimator derived from \(\varphi_{\text{MAR}}\) at the estimated policy \(\widehat d_{\kappa}^{*}\) as \(\widehat{\psi}_{\text{MAR}}(\widehat d_{\kappa}^{*}) = \mathbb{P}_n(\widehat{\varphi}_{\text{MAR}}(\widehat d_{\kappa}^{*}))\). Assume the conditions of Proposition 6, in addition to Assumption [ass:margincondition] and \(\|(\widehat \tau - \widehat \tau^*_{\kappa}) - (\tau - \tau^*_{\kappa}) \|^{1 + \alpha}_{\infty} = o_{\mathbb{P}}(n^{-1/2} )\). Then, \(\sqrt{n}(\widehat{\psi}_{\text{MAR}}(\widehat d_{\kappa}^{*}) - V(d^*_{\kappa})) \rightsquigarrow \mathcal{N}(0, \textrm{Var}(\varphi_{\text{MAR}}(d^*_{\kappa})))\).
We state the proof for Theorem 14; the proof for Theorem 15 is analogous. In this section, we let \(\mathbb{P}(f) \equiv \mathbb{E} \left[f(\mathbf{X})\right]\) for a generic function \(f\). First, we rewrite the identified value function in the MCCAR setting as \[\begin{align} V(d^*_{\kappa}) &= \mathbb{E} \left[Y(d^*_{\kappa})\right] = \mathbb{E} \left[\tau(\mathbf{X}) \mathbb{1}(\tau(\mathbf{X}) > \tau^*_{\kappa})] + \mathbb{E}[\nu_0(\mathbf{X})\right]. \end{align}\] Then, recall that we can rewrite the influence function as \[\begin{align} \phi_{\text{MCCAR}}(d) &= \frac{R\mathbb{1}(A = d)}{\eta_{d}(\mathbf{X})}\left(Y - \nu_{d}(\mathbf{X})\right) + \nu_{d}(\mathbf{X}) \nonumber \\ &= \left( \left\{\frac{RA}{\eta_1(\mathbf{X})} - \frac{R(1 - A)}{\eta_0(\mathbf{X})}\right\}(Y - \nu_a(\mathbf{X})) + \nu_1(\mathbf{X}) - \nu_0(\mathbf{X})\right)\big(\mathbb{1} (\tau > \tau^*_{\kappa}) \big) \\ & \qquad \qquad + \left(\frac{R(1-A)}{\eta_0(\mathbf{X})}\right)(Y - \nu_0(\mathbf{X})) + \nu_0(\mathbf{X}) \\ &= \xi(\mathbf{X}) d^{*}_{\kappa} + \left(\frac{R(1-A)}{\eta_0(\mathbf{X})}\right)(Y - \nu_0(\mathbf{X})) + \nu_0(\mathbf{X}), \end{align}\] where \(\xi(\mathbf{X}) \equiv \left\{\frac{RA}{\eta_1(\mathbf{X})} - \frac{R(1 - A)}{\eta_0(\mathbf{X})}\right\}(Y - \nu_a(\mathbf{X})) + \nu_1(\mathbf{X}) - \nu_0(\mathbf{X})\).
Dropping arguments for notational convenience, we decompose the difference \[\begin{align} \widehat{\Phi}_{\text{MCCAR}}(\widehat d_{\kappa}^{*}) - V(d^*_{\kappa}) &= \mathbb{P}_n\{\widehat{\phi}_{\text{MCCAR}}(\widehat d_{\kappa}^{*})\} - \mathbb{E} \{Y(d^*_{\kappa})\} \nonumber \\ &= \mathbb{P}_n(\widehat \xi\widehat d^{*}_{\kappa}) - \mathbb{P} (\xi d_{\kappa}^{*}) \tag{22} \\ &\quad + \mathbb{P}_n\Big(\Big\{\frac{R(1-A)}{\widehat \eta_0}\Big\}(Y - \widehat \nu_0) + \widehat \nu_0\Big) - \mathbb{P}(v_0), \tag{23} \end{align}\] where (22 ) follows by iterated expectation. (23 ) can be controlled following standard arguments. Similar to [21], we may expand (22 ) as \[\begin{align} \mathbb{P}_n(\widehat \xi\widehat d^{*}_{\kappa}) - \mathbb{P} (\xi d_{\kappa}^{*}) &= \mathbb{P}_n(\widehat \xi\widehat d^{*}_{\kappa} - \tau^*_{\kappa}(\widehat d^*_{\kappa} - \kappa)) - \mathbb{P} (\xi d_{\kappa}^{*} - \tau^*_{\kappa}(d^*_{\kappa} - \kappa)) \\ &\quad + \mathbb{P}_n(\tau^*_{\kappa}(\widehat d^*_{\kappa} - \kappa)) - \mathbb{P}(\tau^*_{\kappa}(d^{*}_{\kappa} - \kappa)). \end{align}\]
By the definition of \(d^{*}_{\kappa}\), \(\mathbb{P}(d^{*}_{\kappa}) = \kappa\), rendering the last population term zero. To simplify the handling of the third term, we assume that there are no ties in the estimated scores \(\widehat \tau(\mathbf{X})\) at the boundary \(\widehat \tau^{*}_{\kappa}\), and so \(\mathbb{P}_n(\widehat d^{*}_{\kappa}) = \kappa\) holds exactly by construction thus eliminating the third term. If this condition is violated, \(\widehat d^{*}_{\kappa}\) may instead be implemented as a stochastic policy that randomly treats a subset of tied individuals. We refer to [21] for a further discussion on how to modify the proof in this case.
We then further decompose (22 ) as \[\begin{align} \mathbb{P}_n(\widehat \xi \widehat d^{*}_{\kappa} - \tau^{*}_{\kappa}(\widehat d^{*}_{\kappa} - \kappa)) - \mathbb{P} (\xi d^{*}_{\kappa} -\tau^{*}_{\kappa}(d^*_{\kappa} - \kappa)) &= (\mathbb{P}_n- \mathbb{P})((\widehat \xi \widehat d^{*}_{\kappa} - \xi d^{*}_{\kappa}) - \tau^{*}_{\kappa}(\widehat d^{*}_{\kappa} - d^{*}_{\kappa})) \tag{24} \\ &\quad + (\mathbb{P}_n- \mathbb{P})((\xi - \tau^{*}_{\kappa})d_{\kappa}^{*}) \tag{25} \\ & \quad + \mathbb{P}((\widehat \xi - \tau^{*}_{\kappa}) \widehat d^{*}_{\kappa} - (\xi - \tau^{*}_{\kappa})d^{*}_{\kappa}). \tag{26} \end{align}\]
The empirical process term (24 ) can be handled by standard sample splitting and Donsker class arguments. The second term (25 ) will converge by the Central Limit Theorem. We now focus on the final term of the decomposition (26 ), which may be rewritten as \[\begin{align} \mathbb{P}((\widehat \xi - \tau^{*}_{\kappa}) \widehat d^{*}_{\kappa} - (\xi - \tau^{*}_{\kappa})d^{*}_{\kappa}) = \mathbb{P}((\widehat \xi - \xi)\widehat d^{*}_{\kappa}) + \mathbb{P}((\xi - \tau^{*}_{\kappa})(\widehat d^{*}_{\kappa} - d^{*}_{\kappa})). \end{align}\]
After applying iterated expectation, the first term will be controlled by exactly the same arguments as the derivation of the remainder term (14 ) from the proof of Proposition 5, resulting in the second-order term \[\begin{align} \mathbb{P}((\widehat \xi - \xi)\widehat d^{*}_{\kappa}) = \mathbb{E} \left[ \mathbb{E}[\widehat \xi(\mathbf{X}) - \xi(\mathbf{X}) \mid \mathbf{X}] \widehat d^{*} _{\kappa}(\mathbf{X})\right] \lesssim \sum_{a \in \{0, 1\}} \| \eta_a(\mathbf{x}) - \widehat \eta_a(\mathbf{x})\|\|\nu_a(\mathbf{x}) - \widehat \nu_a(\mathbf{x})\|. \end{align}\] Now, focusing on the second term, \[\begin{align} \mathbb{P}((\xi - \tau_{\kappa}^{*})(\widehat d^{*}_{\kappa} - d^{*}_{\kappa})) &= \int (\tau(\mathbf{x}) - \tau_{\kappa}^{*})(\widehat d^{*}_{\kappa} - d^{*}_{\kappa}) d \mathbb{P}(\mathbf{x}) \\ & \leq \int | \tau (\mathbf{x}) - \tau^{*}_{\kappa}| | \widehat d^{*}_{\kappa} - d^{*}_{\kappa}| d \mathbb{P}(\mathbf{x}) \\ & \leq \int | \tau(\mathbf{x}) - \tau_{\kappa}^{*}||\mathbb{1}(|\tau(\mathbf{x}) -\tau^*_{\kappa}| \leq | (\widehat \tau(\mathbf{x}) - \tau(\mathbf{x})) - (\widehat \tau^*_{\kappa} - \tau^*_{\kappa})|) d\mathbb{P}(\mathbf{x}) \\ & \leq \int | (\widehat \tau(\mathbf{x}) - \tau(\mathbf{x})) - (\widehat \tau^*_{\kappa} - \tau^*_{\kappa})||\mathbb{1}(|\tau(\mathbf{x}) -\tau^*_{\kappa}| \leq | (\widehat \tau(\mathbf{x}) - \tau(\mathbf{x})) - (\widehat \tau^*_{\kappa} - \tau^*_{\kappa})|) d\mathbb{P}(\mathbf{x}) \\ & \leq \| (\widehat \tau(\mathbf{x}) - \tau(\mathbf{x})) - (\widehat \tau^*_{\kappa} - \tau^*_{\kappa}) \|_{\infty} \mathbb{P}(|\tau(\mathbf{x}) -\tau^*_{\kappa}| \leq \| (\widehat \tau(\mathbf{x}) - \tau(\mathbf{x})) - (\widehat \tau^*_{\kappa} - \tau^*_{\kappa})\|_{\infty}) \\ & \lesssim \| (\widehat \tau(\mathbf{x}) - \tau(\mathbf{x})) - (\widehat \tau^*_{\kappa} - \tau^*_{\kappa}) \|_{\infty}^{1 + \alpha}, \end{align}\] where the first line holds by iterated expectation, the last inequality by Assumption [ass:margincondition], and the second inequality by the fact that \[\begin{align} |\mathbb{1}\{\widehat \tau(\mathbf{x}) > \widehat \tau^*_{\kappa}\} - \mathbb{1}\{ \tau(\mathbf{x}) > \tau^*_{\kappa}\} | \leq \mathbb{1}(|\tau(\mathbf{x}) -\tau^*_{\kappa}| \leq | (\widehat \tau(\mathbf{x}) - \tau(\mathbf{x})) - (\widehat \tau^*_{\kappa} - \tau^*_{\kappa})|). \end{align}\]
For our efficiency simulations, we use the following base settings, \[\begin{align} X &\sim \text{Uniform}(-1, 1) \\ \pi(X) &= .4\sin(2\pi X) + .5\\ R \mid X &\sim \text{Bernoulli}(\pi(X)) \\ \gamma(X) &= .4\sin(2\pi X - .72\pi) + .5\\ \lambda(X, Y) &= \gamma(X) \\ A \mid X &\sim \text{Bernoulli}(\gamma(X)) \\ Y \mid A &\sim \text{Bernoulli}(0.5). \end{align}\] We do not use a treatment effect in order to have a constant residual term. Hence, \(\lambda(X, Y) = P(A|X,Y) = P(A|X) = \gamma(X)\). All default constants were chosen so \(\pi, \gamma, \lambda \approx 0.5\). We run 100 iterations for each configuration and report the average. For each iteration, we sample \(n = 10,000\) observations. In total, all three experiments complete in 1-2 minutes on a personal laptop. All code to reproduce these experiments is available in the supplementary materials.
Let \(l^* \in \{0.1, 0.15, \dots, 0.8, 0.85\}\) be the desired level for \(\lambda\). We then set \(\lambda(X, Y) = .09\sin(2\pi X - .72\pi) + l^*\) while keeping all other variables constant.
Let \(p^* \in \{0.1, 0.15, \dots, 0.8, 0.85\}\) be the desired level for \(\pi\). We then set \(\pi(X, Y) = .09\sin(2\pi X) + p^*\) while keeping all other variables constant.
Let \(c^* \in \{0.0, 0.05, \dots, 0.9, 0.95\}\) be the correlation. We then set \(\lambda(X, Y) = .4\sin(2\pi X - c^*\pi) + .5\) while keeping all other variables constant.
Recall that \(\widehat \beta_{d}(\mathbf{X}) = \mathbb{E}[Y\widehat \lambda_{d}(\mathbf{X}, Y) \mid \mathbf{X}] = \int y \widehat \lambda_{d}(\mathbf{X}, y) d \mathbb{P}(Y \leq y \mid X)\) and \(\widehat \gamma_{d} = \mathbb{E}[\widehat \lambda_{d}(\mathbf{X}, Y) \mid \mathbf{X}] = \int \widehat \lambda_{d}(\mathbf{X}, y) d \mathbb{P}(Y \leq y \mid X)\). We may thus estimate \(\widehat \beta_{d}\) and \(\widehat \gamma_{d}\) through either regression or conditional density estimation. Due to the difficulties of the latter, we propose estimating the density only for discrete outcomes and using regression and a three-way sample splitting scheme when \(Y\) is continuous. Algorithm 8 handles this first case for binary \(Y\), but can be easily modified to accommodate categorical outcomes. Below, we let \(\mu(\mathbf{X}) = P(Y = y \mid \mathbf{X})\) be the conditional density.
Algorithm 8 is not valid for continuous outcomes. To avoid conditional density estimation, we instead estimate \(\widehat \beta\) and \(\widehat \gamma\) using a nested cross-fitting procedure. There are many ways to split the data this way. For example, one could split the data into thirds and estimate six different CATE models (one for each permutation of which fold estimates nuisances, which estimates \(\widehat \beta, \widehat \gamma\) and which estimates the CATE). To reduce computational complexity, Algorithm 9 maintains two main splits for two different CATE models \(\widehat \tau_{1}, \widehat \tau_{2}\), but further splits each nuisance fold in half to separately estimate \(\widehat \lambda\) versus \(\widehat \beta\), \(\widehat \gamma\) in nested cross-fitting loop.
All code to reproduce the following experiments is available at https://github.com/jesund/learningwithmissingtreatments.
We generate \(n = 1,000\) samples using the data-generating processes described in Section 4.1. We calibrated the constant \(C\) in the
missingness functions \(\pi\) empirically to match target observation rates. We report our results with observation rates of \(25\%\) and \(75\%\) in the
main text, and below report our results with the other observation rates we tested in Fig. 10. For each observation rate in the MAR and MCCAR setting, we ran 500 iterations. All nuisance functions used their true values
and were not estimated. We then added randomly sampled noise \(\epsilon \sim N(-n^{-\alpha}, n^{-2\alpha})\) to each nuisance function, where \(\alpha\) is the rate of convergence,
i.e. setting \(\alpha = \frac{1}{2}\) corresponds to \(\sqrt n-\)rates of convergence. The noisy nuisance functions were then used to compute pseudo-outcomes \(\widetilde{\varphi}_{\textrm{MAR}}\) and \(\widetilde{\varphi}_{\textrm{MCCAR}}\). Pseudo-outcomes for oracle models were computed using the noise-free nuisance functions. We used a
smooth.spline from base R with the default parameters to estimate the CATE for all models, including the oracle, by regressing the pseudo-outcomes on \(X\). We then calculated the RMSE on a test set of size
\(1,000\) and report the median over the 500 iterations. All experiments were run on a personal laptop with 12 CPU cores and took 5-10 minutes to complete.


Figure 10: Median RMSE of CATE estimates comparing \(\widehat \tau_{\textrm{MAR}}\) (blue), \(\widehat \tau_{\textrm{MCCAR}}\) (orange), and \(\widehat \tau_{\textrm{Oracle}}\) (green) in the MAR (a) and MCCAR (b) settings over 500 iterations against the rate of noise added to all nuisance functions. In the MAR setting (top), \(\widehat \tau_{\textrm{MCCAR}}\) remains biased even after more treatments are observed, while in the MCCAR setting (bottom) both estimators are unbiased and \(\widehat \tau_{\textrm{MAR}}\) achieves slightly faster convergence. This pattern holds across all observation rates including \(25\%\) and \(75\%\)..
Our two high-dimensional experiments generate data following the MAR and MCCAR data-generating processes given in Section 4.1.2. We run 50 iterations for each sample size \(n \in \{1,000,
2,000, 5,000, 10,000\}\) under both settings. For each iteration, we estimate the CATE following Algorithm 8. We use RandomForest for all nuisance function estimation. All propensity scores are
clipped to be in \([0.01, 0.99]\). We found that using a shallower RandomForest for the CATE model worked best based on preliminary experiments, and so we set max_depth = 3. We compare our
estimates to an oracle CATE that uses the true values for all nuisance functions in the input to the RandomForest. All pseudo-outcomes are clipped to be within the 1st and 99th quantiles. We generate a separate test sample of the same size for
each iteration. We report the average RMSE over each iteration for each sample size and setting. All experiments were run on a personal laptop with 12 CPU cores and took 10-20 minutes to complete.
To test the robustness properties of our proposed estimators, we run a variant of the high-dimensional experiment. To introduce model misspecification, we systematically train nuisance functions with \(\verb|RandomForest|\) with \(\verb|max_depth=1|\). The shallow trees are unable to capture the complexity of the data-generating process in Section 4.1.2, as each nuisance includes multiple variables and some have interaction terms. All other hyperparameters and experimental details remain the same as described above.
Table 2 shows the results of the experiment in the MCCAR setting. Following Proposition 5, if \(\nu_d\) or \(\eta_d\) is misspecified, the MCCAR estimator still converges. On the other hand, the MCCAR estimator is not robust to misspecification of both \(\nu_d\) and \(\eta_d\).
| \(N\) | \(\widehat \tau_{MCCAR}\) | \(\widetilde \nu_d\) | \(\widetilde \eta_d\) | \(\widetilde \nu_d, \widetilde \eta_d\) |
|---|---|---|---|---|
| 1000 | 0.095 | 0.078 | 0.155 | 0.130 |
| 5000 | 0.021 | 0.019 | 0.040 | 0.048 |
| 10000 | 0.014 | 0.014 | 0.031 | 0.044 |
All instructions for downloading the data and code to run the semi-synthetic experiments are provided at https://github.com/jesund/learningwithmissingtreatments. All experiments were run on a computing cluster. MTRNET was tuned and trained on an NVIDIA L40S GPU, while other models were trained using nodes with 24 CPUs. For a single replication, hyperparameter tuning and training takes 1-2 hours for MTRNET versus 5-10 minutes for the other methods.
This dataset comes from a 2008 field experiment studying the effects of social pressure on voter turnout [35]. We focused on two experimental arms: the “Neighbors” treatment (\(A = 1\)), where potential voters received mailings showing their neighbors’ voting records, and the control group (\(A = 0\)), which received no mailings. The outcome \(Y\) is a binary indicator for whether the individual voted in the August 2006 primary election. We preprocessed the data by converting the year of birth to age, removing the household ID and cluster ID, and converting categorical variables to dummy variables using one-hot encoding, dropping the first category of each variable to avoid multicollinearity.
This dataset comes from a four-year longitudinal study examining the effect of class size on student achievement [36]. We focused on kindergarten students, comparing two experimental conditions: small class (13-17 students) and regular class (22-25 students), excluding the regular-with-aide condition to maintain binary treatment. We preprocessed the data by filtering for kindergarten students only and excluded students in the regular-with-aide class type. We created a binary treatment indicator where regular class was coded as 0 and small class as 1. We removed students with missing test scores and created a composite outcome variable \(Y\) by summing the kindergarten test scores across four domains: reading, math, listening, and word skills.
We retained key covariates including gender, race, birth month, school, teacher demographics (gender, race, years of experience, education level, career ladder), and student characteristics (attendance, free lunch eligibility, grade repetition, special education status). Categorical variables were converted to dummy variables using one-hot encoding, and missing values for attendance variables were imputed with the median values.
We artificially induce treatment missingness by generating observation scores \(s_i\) for each individual \(i\) and then sampling binary indicators \(R_i\) based on these scores. We base the missingness mechanism on the most predictive covariate \(X^\star\) for each dataset, identified by feature importance from a RandomForest:
years of teaching experience (STAR) and age (Voting).
For the continuous outcome dataset (STAR), we generate observation scores using both outcomes and \(X^\star\) for each individual \(i\): \[\begin{align} q &= F_{Y}^{-1}(r) \\ s_i &= \mathbb{1}(Y_i < q) (X_{i}^{\star} + 1) + 0.1 \end{align}\] where \(F_Y^{-1}(r)\) is the \(r\)-th quantile of the outcome distribution, \(r\) is the observation rate, and we add one to \(X_{i}^{\star}\) to set the minimum above zero. For the binary outcome dataset (Voting), we use: \[\begin{align} s_i &= r Y_i (X_{i}^{\star} - \min(X^{\star}) + 1) + 0.1 \end{align}\] where we transform \(X^{\star}\) so that the minimum value is one, as in the continuous setting. The baseline probability of 0.1 under both mechanisms ensures that all treatment assignments have some chance of being observed.
We generate missingness dependent only on \(X^{\star}\) by setting observation scores based on thresholds of \(X^{\star}\): \[\begin{align} q &= F_{X^{\star}}^{-1}(r) \\ s_i &= 0.8\mathbb{1}(X_{i}^{\star} < q) + 0.1 \end{align}\] where \(F_{X^\star}^{-1}(r)\) is the \(r\)-th quantile of \(X^{\star}\) and \(r\) is again the observation rate.
After computing observation propensity scores \(s_i\) for each individual, we employ an incremental sampling algorithm to select which treatment assignments are observed. This ensures that as the observation rate \(r\) increases, the set of observed treatments expands monotonically rather than changing completely. We first normalize scores to sampling probabilities: \(p_i = \frac{s_i}{\sum_{j=1}^{n}
s_j}\). We then sample a random permutation of all \(n\) indices without replacement using numpy.random.choice with probabilities \(p_i\), producing an ordered sequence
representing the sampling priority. To achieve an observation rate of \(r\) (e.g., \(r=0.3\) for 30% observed), we set \(R_i = 1\) for the first \(\lfloor nr \rfloor\) observations in the sequence and \(R_i = 0\) for the remainder. This nested structure ensures that increasing \(r\) adds new observations
while retaining all previously observed ones.
We compare the results between the following models in the main text:
DR-MAR Algorithm 8 for the Voting dataset and Algorithm 9 for the STAR dataset.
DR-MCCAR Algorithm 7.
MTRNET Official implementation [18] (https://github.com/mkuzma96/MTRNet).
OR Outcome regression models fit on the nuisance functions using Eq. 7 in the MCCAR setting and Eq. 5 in the MAR setting. We do not use a second-stage regression.
PW-Learner which computes pseudo-scores using only propensity scores [37] and then employs a second stage regression model to estimate the CATE. Below we give the full estimators we use in each setting.
And include the following additional models for comparison in this Appendix:
DR-CC Complete-case doubly robust learner using only observations with \(R_i = 1\).
Risk Model Outcome regression ignoring treatment, predicting \(\mathbb{E}[Y|X]\).
In the MAR setting, we compute the following pseudo-score: \[\begin{align} \frac{AYR}{\widehat \gamma_1(\mathbf{X})\widehat \pi(\mathbf{X},Y)} - \frac{(1-A)YR}{\widehat \gamma_0(\mathbf{X})\widehat \pi(\mathbf{X},Y)}, \end{align}\] where \(\widehat \gamma_1\) is estimated by first estimating \(\widehat \lambda_1\), and then regression \(\widehat \lambda_1\) on \(\mathbf{x}\) (as is done for the DR-MAR estimator). We use the same sample splitting scheme as the DR-MAR estimator. Under MAR Assumptions 3 for a specific \(\mathbf{x}\), it follows that \[\begin{align} \mathbb{E} \left[\frac{AYR}{\gamma_1(\mathbf{x})\pi(\mathbf{x}, Y)} \mid \mathbf{X} = x\right] = \mathbb{E} \left[\frac{Y\mathbb{E} \left[AR \mid \mathbf{X}= \mathbf{x}, Y\right]}{\gamma_1(\mathbf{x})\pi(\mathbf{x}, Y)} \mid \mathbf{X} = \mathbf{x}\right] = \frac{\mathbb{E} \left[Y \lambda_1(\mathbf{x}, Y)\mid \mathbf{X}= \mathbf{x}\right]}{\gamma_1(\mathbf{x})} = \frac{\beta_1(\mathbf{x})}{\gamma_{1}(\mathbf{x})}, \end{align}\] which is exactly the first term of the identified cate in Eq. 5 . The second term follows a similar proof.
In the MCCAR setting under Assumption 2, we use the following pseudo-score: \[\begin{align} \frac{AYR}{\widehat \eta_1(\mathbf{x})} - \frac{(1-A)YR}{\widehat \eta_0(\mathbf{x})}. \end{align}\] In practice, \(\widehat \eta_d\) may be estimated using separate models for \(P(A = d\mid \mathbf{X}= \mathbf{x}, R = 1)\) and \(P(R = 1 \mid \mathbf{X}= \mathbf{x})\), which we find performs slightly better than estimating the joint probability of both treatment and missingness. We follow the same sample splitting scheme as the DR-MCCAR estimator. It follows that \[\begin{align} \mathbb{E} \left[\frac{AYR}{\eta_1(\mathbf{x})} \mid \mathbf{X}= \mathbf{x}\right] = \frac{\mathbb{E} \left[Y \mid A = 1, \mathbf{X}= \mathbf{x}\right]\eta_1(\mathbf{x})}{\eta_1(\mathbf{x})} = \mathbb{E} \left[Y \mid A = 1, \mathbf{X}= \mathbf{x}\right], \end{align}\] where we used the fact that \(\eta_1(\mathbf{x}) = P(R = 1, A = 1 \mid \mathbf{X}= \mathbf{x})\). The second term follows a similar proof.
We use an ensemble for all nuisance function estimation and pseudo-outcome regression of five base learners: ridge/logistic regression, elastic net/lasso, random forest, XGBoost, and K-nearest neighbors. Propensity scores for all PW-Learner and DR-Learner variants are clipped to [0.1, 0.9]. After computing pseudo-outcomes, we winsorize at the 5th and 95th percentiles for the STAR dataset to reduce the influence of extreme values. For Voting, no winsorization is applied due to the larger sample size providing more stable estimates.
| Model | Hyperparameter | Values |
|---|---|---|
| Ridge / Logistic | \(\alpha\) / \(C\) | \(\{0.001, 0.01, 0.1, 1, 10, 100\}\) |
| Penalty | L2 (ridge penalty) | |
| Elastic Net / Lasso | \(\alpha\) / \(C\) | \(\{0.01, 0.1, 1, 10\}\) |
| \(L_1\) ratio | \(\{0.1, 0.3, 0.5, 0.7, 0.9\}\) (regression) / 1.0 (classification) | |
| Random Forest | n_estimators | 500 |
| max_depth | \(\{3, 5, 10, \text{None}\}^*\) | |
| max_features | \(\{\text{sqrt}, \text{None}\}\) | |
| min_samples_split | \(\{2, 5\}^*\) | |
| XGBoost | n_estimators | \(\{100, 300\}\) |
| max_depth | \(\{3, 5\}\) | |
| learning_rate | \(\{0.01, 0.1\}\) | |
| K-Nearest Neighbors | \(k\) | \(\{3, 5, 7, 11, 15\}\) |
| weights | \(\{\text{uniform}, \text{distance}\}\) |
\(^*\)For Voting: max_depth \(\in \{3, 5, 10, 20\}\), min_samples_split \(\in \{10, 20\}\)
We use 3-fold cross-validation for hyperparameter tuning. Classification models are optimized for negative log-loss, while regression models are optimized for negative mean squared error. We standardize all continuous features for ridge/logistic regression, elastic net/lasso, and KNN models.
Final ensemble predictions are weighted averages of the five base learner predictions. For STAR, we use uniform ensemble weights. For Voting, we learn optimal weights using non-negative least squares on 2-fold cross-validated predictions, normalized to sum to one.
Hyperparameters are selected using Optuna [38]with the Tree-structured Parzen Estimator (TPE) sampler. We use 100 trials for STAR and 50 trials for Voting (due to computational cost) using the search space in Table 4. We hold out 20% of the training data as a validation set and optimize for AUPEC on the validation split. All features are standardized to zero mean and unit variance before training. We then retrain the tuned MTRNet on the full training set using the selected hyperparameters.
| Hyperparameter | Search Space |
|---|---|
| Representation layer size | \(\{50, 100, 200\}\) |
| Hypothesis layer size | \(\{50, 100, 200\}\) |
| Learning rate | \(\{0.0001, 0.0005, 0.001, 0.005, 0.01\}\) |
| Dropout rate | \(\{0.1, 0.2, 0.3\}\) |
| Training iterations | \(\{100, 200, 300\}\) |
| Batch size | \(\{100, 300, 500, 1000, 1500\}^*\) |
| \(\alpha\) (IPW weight) | \(\{10^{k/2} : k \in \{-4, -3.5, \ldots, 2\}\}\) |
| \(\beta\) (representation penalty) | \(\{10^{k/2} : k \in \{-4, -3.5, \ldots, 2\}\}\) |
| \(\lambda\) (regularization) | \(\{0.00005, 0.0001, 0.0005\}\) |
\(^*\)For Voting: batch size \(\in \{1000, 2500, 5000\}\)
We split the data 80-20 for training and testing, and evaluate all methods on the held-out test set.
For a given treatment budget \(\kappa \in [0,1]\) , we compute the policy value as follows. First, we treat the decision rule generated by \(\widehat \tau\) as fixed because it was trained on a separate sample. We then let \(\widehat \tau^*_{\kappa} = \max\{F_{\widehat{\tau}}^{-1}(\kappa), 0\}\) be the treatment threshold. We assign treatment to individuals with \(\widehat{\tau}(\mathbf{X}_i) > \tau^*_{\kappa}\), creating treatment decisions \(D_i \in \{0,1\}\). We estimate the policy value by \(V_{\kappa}(\widehat \tau)= \mathbb{E} \left[Y(1) \mid D = 1\right]P(D = 1) + \mathbb{E}[Y(0)\mid D = 0] P(D = 0)\), which we compute empirically using the observed outcomes under random treatment assignment. The AUPEC summarizes performance across all budgets: \(\text{AUPEC}(\widehat{\tau}) = \int_0^1 \left[ V_{\kappa}(\widehat{\tau}) - V_{\kappa}^{\text{random}} \right] d\kappa\), where \(V_\kappa^{\text{random}} = \kappa \mathbb{E}[Y(1)] + (1-\kappa) \mathbb{E}[Y(0)]\) is the expected policy value under random assignment.
For each dataset and missingness mechanism, we vary the observation rate \(r \in \{0.1, 0.2, \ldots, 0.9\}\) to simulate different levels of treatment missingness. For each combination of dataset, mechanism, and observation rate, we run 10 replications (STAR) or 5 replications (Voting) with different random seeds. We conduct three separate robustness experiments, each varying:
Missingness generation: varying the random permutation used to sample
Train-test split: varying the random seed for data splitting
Model initialization: varying random seeds for model training
We report mean AUPEC for each configuration, allowing us to robustly assess both average performance and variability across different sources of randomness.
We compare the results for the standard doubly-robust DR-Learner using only complete observations, and the Risk Models which predict the risk of an adverse outcome. We tested two rules based on the risk model: (1) treat those with the highest risk of the adverse outcome, and (2) treat those in the middle of the risk distribution. These rules reflect different heuristics policy-makers use in practice when assigning treatment. We compare our results to the best performing rule for each dataset.
Figure 12: AUPEC with 95% confidence intervals, where zero indicates no gain over random policy. Results show average performance across treatment observation rates for STAR (10 iterations) and Voting (5 iterations) datasets under MAR and MCCAR missingness. Left: varying model initialization seed. Right: varying train-test split seed. Train-test split variation introduces substantially more uncertainty than model seed variation. The results support our main conclusions: in the MAR setting, DR-MAR outperforms all other estimators, while in the MCCAR setting all models are competitive. We omit MTRNET from this comparison due to computational costs.. a — Changing model seed, b — Changing train-test split seed
| DR-CC | DR-MAR | DR-MCCAR | IPW | MTRNET | OR | Random | Risk | ||
| Obs. Rate | \(\kappa\) | ||||||||
| 20% | 10% | 0.303 | 0.306 | 0.301 | 0.306 | 0.304 | 0.307 | 0.305 | 0.305 |
| 25% | 0.312 | 0.320 | 0.309 | 0.319 | 0.318 | 0.323 | 0.317 | 0.319 | |
| 50% | 0.332 | 0.341 | 0.309 | 0.340 | 0.340 | 0.346 | 0.337 | 0.337 | |
| 50% | 10% | 0.303 | 0.309 | 0.304 | 0.306 | 0.304 | 0.307 | 0.305 | 0.305 |
| 25% | 0.315 | 0.321 | 0.318 | 0.318 | 0.317 | 0.322 | 0.317 | 0.319 | |
| 50% | 0.336 | 0.343 | 0.340 | 0.340 | 0.340 | 0.345 | 0.337 | 0.337 | |
| 80% | 10% | 0.306 | 0.308 | 0.309 | 0.305 | 0.305 | 0.307 | 0.305 | 0.305 |
| 25% | 0.320 | 0.320 | 0.320 | 0.317 | 0.319 | 0.322 | 0.317 | 0.319 | |
| 50% | 0.340 | 0.343 | 0.342 | 0.340 | 0.343 | 0.345 | 0.337 | 0.337 |
| DR-CC | DR-MAR | DR-MCCAR | IPW | MTRNET | OR | Random | Risk | ||
| Obs. Rate | \(\kappa\) | ||||||||
| 20% | 10% | 0.308 | 0.307 | 0.306 | 0.306 | 0.306 | 0.307 | 0.305 | 0.305 |
| 25% | 0.321 | 0.321 | 0.320 | 0.318 | 0.318 | 0.318 | 0.317 | 0.319 | |
| 50% | 0.341 | 0.339 | 0.340 | 0.341 | 0.341 | 0.340 | 0.337 | 0.337 | |
| 50% | 10% | 0.308 | 0.308 | 0.306 | 0.307 | 0.306 | 0.305 | 0.305 | 0.305 |
| 25% | 0.321 | 0.322 | 0.319 | 0.319 | 0.319 | 0.318 | 0.317 | 0.319 | |
| 50% | 0.341 | 0.341 | 0.342 | 0.338 | 0.343 | 0.341 | 0.337 | 0.337 | |
| 80% | 10% | 0.307 | 0.308 | 0.309 | 0.307 | 0.306 | 0.305 | 0.305 | 0.305 |
| 25% | 0.320 | 0.321 | 0.321 | 0.320 | 0.318 | 0.319 | 0.317 | 0.319 | |
| 50% | 0.342 | 0.341 | 0.342 | 0.340 | 0.344 | 0.342 | 0.337 | 0.337 |
| DR-CC | DR-MAR | DR-MCCAR | IPW | MTRNET | OR | Random | Risk | ||
| Obs. Rate | \(\kappa\) | ||||||||
| 20% | 10% | 1889.8 | 1899.1 | 1888.6 | 1885.9 | 1886.7 | 1891.4 | 1888.8 | 1890.8 |
| 25% | 1894.4 | 1905.1 | 1895.5 | 1886.3 | 1889.8 | 1900.2 | 1891.1 | 1895.1 | |
| 50% | 1899.3 | 1909.7 | 1900.2 | 1883.0 | 1897.8 | 1903.2 | 1895.0 | 1898.3 | |
| 50% | 10% | 1897.6 | 1896.2 | 1893.7 | 1884.3 | 1890.2 | 1894.0 | 1888.8 | 1890.8 |
| 25% | 1902.4 | 1902.1 | 1902.4 | 1884.2 | 1895.8 | 1901.3 | 1891.1 | 1895.1 | |
| 50% | 1906.3 | 1919.2 | 1907.7 | 1888.8 | 1900.3 | 1908.2 | 1895.0 | 1898.3 | |
| 80% | 10% | 1898.3 | 1899.1 | 1894.8 | 1888.7 | 1890.8 | 1888.8 | 1888.8 | 1890.8 |
| 25% | 1903.2 | 1908.3 | 1903.4 | 1890.7 | 1894.3 | 1888.6 | 1891.1 | 1895.1 | |
| 50% | 1909.0 | 1915.9 | 1911.0 | 1895.2 | 1899.9 | 1892.8 | 1895.0 | 1898.3 |
| DR-CC | DR-MAR | DR-MCCAR | IPW | MTRNET | OR | Random | Risk | ||
| Obs. Rate | \(\kappa\) | ||||||||
| 20% | 10% | 1891.5 | 1892.8 | 1891.2 | 1884.4 | 1887.2 | 1892.4 | 1888.8 | 1890.8 |
| 25% | 1899.0 | 1896.6 | 1898.6 | 1885.7 | 1890.4 | 1897.2 | 1891.1 | 1895.1 | |
| 50% | 1905.6 | 1905.5 | 1907.2 | 1885.9 | 1891.8 | 1906.7 | 1895.0 | 1898.3 | |
| 50% | 10% | 1895.5 | 1900.0 | 1898.2 | 1890.0 | 1890.5 | 1895.7 | 1888.8 | 1890.8 |
| 25% | 1904.0 | 1906.0 | 1907.0 | 1893.4 | 1899.6 | 1904.2 | 1891.1 | 1895.1 | |
| 50% | 1916.3 | 1915.2 | 1914.4 | 1896.4 | 1907.6 | 1917.3 | 1895.0 | 1898.3 | |
| 80% | 10% | 1898.8 | 1900.3 | 1900.0 | 1890.8 | 1892.0 | 1899.0 | 1888.8 | 1890.8 |
| 25% | 1905.4 | 1906.4 | 1905.7 | 1893.7 | 1897.6 | 1906.4 | 1891.1 | 1895.1 | |
| 50% | 1916.9 | 1917.4 | 1917.0 | 1896.4 | 1906.0 | 1919.8 | 1895.0 | 1898.3 |
Code for both experiments is available on Github at https://github.com/jesund/learningwithmissingtreatments↩︎