Keywords: Continuous treatment effects, surrogate outcomes, double robustness.

1 Introduction↩︎

In many causal inference applications, the primary outcomes are missing for a non-trivial number of observations. For instance, in studies on long-term health effects of medical interventions, some measurements require expensive testing and a loss to follow-up is common [1]. In evaluating commercial online ad effectiveness, some individuals may drop out from the panel because they use multiple devices [2], leading to missing revenue measures. In many of these studies, however, there often exist short-term outcomes that are easier and faster to measure, e.g., short-term health measures or an online ad’s click-through rate, that are observed for a greater share of the sample. These outcomes, which are typically informative about the primary outcomes themselves, are referred to as surrogate outcomes or surrogates.

There is a rich causal inference literature addressing missing outcome data. Simply restricting to data with observed primary outcomes may induce strong bias [3]. Ignoring unlabeled data also reduces the effective sample size for estimating the treatment effects and inflates the variance. [4] considered the missing completely at random (MCAR) setting and showed that incorporating unlabeled data reduces variance. [5] generalized the results to missing-at-random (MAR) settings where the unlabeled data has a much larger size than the labeled data. [6] further examined the role of surrogates in datasets with limited primary outcomes and showed efficiency gains after including surrogates and unlabeled data in the analysis. [7] proposed a generalized kernel ridge regression framework, which incorporates information on surrogate outcomes in treatment effect estimation. See also [8]–[10] for relevant discussions.

Continuous treatments appear in many applications; e.g., waiting time before follow-up, percent of discount, and drug dosage. Existing estimation procedures include outcome modeling [11] and treatment process modeling [12]. [13] proposed doubly robust methods that model both the outcome and the treatment process, and enjoy appealing robustness properties. [14] further examined high-order estimators to achieve faster rates.

In this paper, we consider the estimation of dose-response functions with limited primary outcome data. We propose novel doubly robust methods using both labeled and unlabeled data, with the help of surrogates. Our approach avoids the potential selection bias caused by restricting only to labeled data and provably reduces the variance. Importantly, we study the theoretical properties of the estimator proposed, and establish its asymptotic normality to facilitate statistical inference. Our work serves as a counterpart for [6] in continuous treatment effects settings and enriches the existing dose-response estimation literature [13], [14].

The rest of this paper is organized as follows: In Section 2 we introduce the problem setup and notation. Section 3 provides assumptions to identify the continuous treatment effect. Novel methodology and theoretical guarantees are discussed in Section 4. Our simulation study in Section 5 demonstrates the performance of our method. Finally we conclude with a discussion in Section 6. All proofs, additional simulation results and a real data example are included in supplementary materials.

2 Setup and Notation↩︎

In this section we introduce the problem of estimating continuous treatment effects with surrogates. We formalize the problem using potential outcomes framework [15], [16] and introduce notation to present our results concisely.

2.1 Data Structure↩︎

Suppose we have access to two datasets \(\mathcal{L}\) and \(\mathcal{U}\) from a randomized experiment/observational study. The labeled dataset is \(\mathcal{L} = \{ \mathbf{Z}_i = (\mathbf{V}_i, \mathbf{S}_i, A_i,Y_i,R_i=1), 1\leq i \leq n_1 \}\), where \(\mathbf{V}\) is the set of covariates, \(\mathbf{S}\) is the vector of surrogates, \(A\) is the continuous treatment, and \(Y\) is the primary outcome. Here \(R\) is an indicator with \(R=1\) if a sample is from the labeled data \(\mathcal{L}\) and \(R=0\) otherwise. The unlabeled dataset consists of samples without primary outcomes \(\mathcal{U} = \{\mathbf{Z}_i = (\mathbf{V}_i, \mathbf{S}_i, A_i, R_i = 0), n_1+1 \leq i \leq n_1 + n_2 \}\). Hence the total sample size is \(n=n_1 + n_2\) and both \(\mathcal{L}\) and \(\mathcal{U}\) contain information on treatment \(A\), surrogates \(\mathbf{S}\) and covariates \(\mathbf{V}\), but the primary outcome is missing in the unlabeled dataset \(\mathcal{U}\) due to, for example, loss to follow-up. Let \(\mathbf{X}= (\mathbf{V}, \mathbf{S})\) be the union of covariates and surrogates. Note that we will present the results for \(\mathbf{S}\neq \emptyset\), but \(\mathbf{S}= \emptyset\) can be viewed as a special setting where our methodology still applies.

2.2 Estimand and Nuisance Functions↩︎

Now we define the continuous treatment effects of interest. We use the random variable \(Y^a\) to denote the potential (counterfactual) outcome we would have observed had a subject received treatment \(A = a\), which may be contrary to the observed \(Y\). The continuous treatment effect (or dose-response function) is defined as \[\label{eq:dose-response} \theta(a) = \mathbb{E}[Y^a].\tag{1}\] Without additional causal assumptions introduced in Section 3, \(\theta(a)\) involves counterfactual outcome \(Y^a\) and cannot be identified by observed data. Note that since \(\mathbf{S}\) contains post-treatment surrogate outcomes and may be affected by the treatment, we also use potential outcomes \(\mathbf{S}^a\) to denote the surrogate under treatment \(A=a\).

We further introduce some nuisance functions that are not of primary interest, but which our estimation method depends on. Let \(\mu(A, \mathbf{X}, 1) = \mathbb{E}[Y|A,\mathbf{X},R=1]\) be the regression function of the primary outcome in the labeled population \(R=1\). Denote the function obtained by further regressing \(\mu(A,\mathbf{X}, 1)\) on \((A,\mathbf{V})\) as \(\tau(A,\mathbf{V}) = \mathbb{E}[\mu(A,\mathbf{X}, 1)|A,\mathbf{V}]\) , where the expectation is over the conditional distribution of \(\mathbf{S}\) given \(A,\mathbf{V}\). Denote the conditional density of \(A\) given \(\mathbf{V}\) as \(\pi(a|\mathbf{V})\) and the marginal density of \(A\) as \(f(a) = \int \pi(a|\mathbf{v}) d\mathbb{P}(\mathbf{v})\). The ratio of the marginal density and conditional density is denoted as \(w(a, \mathbf{v}) = f(a)/\pi(a|\mathbf{v})\). \(w(a,\mathbf{v})\) is known as a stabilized weight in the literature [17], [18]. The propensity score for \(R\) (i.e. the conditional probability that the primary outcome is observed) is denoted as \(\rho(A,\mathbf{X}) = \mathbb{P}(R=1|A,\mathbf{X})\).

Let \((\mathcal{V}, \mathcal{S}, \mathcal{A})\) denote the support of \((\mathbf{V}, \mathbf{S}, A)\). For a (possibly random) function \(f\) on variables \(\mathbf{Z}\) we use \(\mathbb{P}_n [f(\mathbf{Z})]\) or \(\mathbb{P}_n [f]\) to denote the sample average \(\frac{1}{n}\sum_{i=1}^n f(\mathbf{Z}_i)\) on a sample of size \(n\). The sample over which averages are taken should be clear from context. We use \(\mathbb{P}[f] = \int f(\mathbf{z}) d \mathbb{P}(\mathbf{z})\) to denote the expectation of \(f(\mathbf{Z})\) where only randomness of \(\mathbf{Z}\) is considered and \(f\) is conditioned on. Finally we use \(\|f\|_{\infty} = \sup_{\mathbf{z}\in \mathcal{Z}} |f(\mathbf{z})|\) to denote the uniform norm, \(\|f\|_a = \left\{\int f^2(\mathbf{z}) d \mathbb{P}(\mathbf{z}|A=a)\right\}^{1/2}\) to denote the \(L_2\)-norm with respect to the conditional distribution \(\mathbf{Z}|A=a\) and \(\|f\|_2 = \left\{\int f^2(\mathbf{z}) d \mathbb{P}(\mathbf{z})\right\}^{1/2}\) to denote the usual \(L_2\)-norm.

3 Identification↩︎

In this section we discuss sufficient conditions to identify the dose-response function 1 , summarized as follows:

Assumption 1 is also known as the stable unit treatment value assumption (SUTVA), and requires an absence of interference between different individuals in the study. Assumption 2 is commonly used to identify treatment effects and holds in a randomized experiment or observational study with all confounders measured. Assumption 3 ensures whether the primary outcome is observed only depends on covariates \(\mathbf{V}\), surrogates \(\mathbf{S}\), and treatment \(A\) so that the distributions of labeled and unlabeled data are comparable after conditioning on \(\mathbf{X},A\). Assumption 4 says every subject has some chance to receive treatment \(A=a\) and has the primary outcome observed. An illustrative causal graph is shown in Figure 1. The readers are referred to [19] for detailed discussions on identifying continuous treatment effects and [6] for the role of surrogates in identifying treatment effects. With all these assumptions, the treatment effects of interest can be identified using the observable distribution as summarized in Theorem 1.

The identification formula ?? suggests the following plug-in style estimator of \(\theta(a)\): first regress \(Y\) on \(A, \mathbf{X}\) in the labeled dataset \(\mathcal{L}\) and obtain an estimator of \(\mu\) as \(\widehat{\mu}\), which is further regressed on \(A, \mathbf{V}\) to get an estimator \(\widehat{\tau}\). Finally we take an average over all the samples to get an estimator of \(\theta(a)\) as \[\label{eq:plugin} \widehat{\theta}(a) = \mathbb{P}_n [\widehat{\tau}(a,\mathbf{V})].\tag{2}\] The performance of such a plug-in style estimator highly depends on the estimation error of \(\widehat{\tau}\). To see this, assume the nuisance estimator \(\widehat{\tau}\) is independent of the samples that we average over, then the conditional bias of \(\widehat{\theta}(a)\) given \(\widehat{\tau}\) is \[\mathbb{E}\left[\widehat{\theta}(a)\right]-\theta(a) = \mathbb{E}\left[\widehat{\tau}(a,\mathbf{V}) - \tau(a,\mathbf{V})\right],\] which solely depends on the estimation error of \(\tau\). When \(\tau\) is hard to estimate (e.g., non-smooth/sparse in high-dimensional problems), the plug-in style estimator may inherit the slow convergence rate of \(\widehat{\tau}\) and have sub-optimal performance.

4 Doubly Robust Estimation↩︎

In this section we present the main results of this paper. We begin with an alternative characterization of the dose-response function 1 in Section 4.1, based on which we propose a doubly robust estimator in Section 4.2. Finally, theoretical guarantees of the proposed method are provided in Section 4.3.

4.1 Doubly Robust Characterization↩︎

Since treatment \(A\) is continuous, the function \(\theta(a)\) in ?? is not pathwise differentiable [20], [21] and we need a novel way to apply semiparametric efficiency theory. The idea is to find a pseudo-outcome \(\varphi(\mathbf{Z}):=\varphi(\mathbf{Z};\mu, \pi, \rho)\) depending on nuisance functions such that \(\mathbb{E}[\varphi(\mathbf{Z};{\mu}, {\pi}, {\rho})| A=a] = \theta(a)\), i.e., regressing \(\varphi(\mathbf{Z};\mu, \pi, \rho)\) on \(A\) yields the target dose-response function ideally with second-order dependence on nuisance estimation error. Following [13], [22], we consider the functional \(\psi = \mathbb{E}[\theta(A)]\), which is pathwise differentiable and admits an efficient influence function. Then the pseudo-outcome \(\varphi(\mathbf{Z};\mu, \pi, \rho)\) is a component in the influence function of \(\psi\). We omit the derivation of the influence function and only present the form of pseudo-outcome. Let \((\bar{\mu}, \bar{\pi}, \bar{\rho})\) be nuisance functions that may not necessarily equal the true \(({\mu}, {\pi}, {\rho})\), and define \[\varphi(\mathbf{Z}; \bar{\mu}, \bar{\pi}, \bar{\rho}) =\left[\frac{R(Y-\bar{\mu}(A,\mathbf{X},1))}{\bar{\rho}(A,\mathbf{X})} + \bar{\mu}(A,\mathbf{X}, 1 ) - \bar{\tau}(A,\mathbf{V}) \right]\frac{\int_{\mathcal{V}} \bar{\pi}(A|\mathbf{v}) d \mathbb{P}(\mathbf{v})}{\bar{\pi}(A|\mathbf{V})} + \int_{\mathcal{V}} \bar{\tau}(A,\mathbf{v}) d \mathbb{P}(\mathbf{v})\] where \(\bar{\tau}(A,\mathbf{V}) = \mathbb{E}[\bar{\mu}(A,\mathbf{X},1)|A,\mathbf{V}]\). The following theorem shows an alternative characterization of the dose-response function \(\theta(a)\) through \(\varphi(\mathbf{Z}; \bar{\mu}, \bar{\pi}, \bar{\rho})\).

Proposition 1 gives the first interpretation of the double robustness of our methods. There are two chances to obtain the dose-response function \(\theta(a)\) by regressing the pseudo-outcome \(\varphi(\mathbf{Z};\bar{\mu}, \bar{\pi}, \bar{\rho})\) on \(A\): we correctly specify either the outcome regression model \(\mu\) or both the propensity score of \(R\) and conditional density of \(A\) given \(\mathbf{V}\) (although see Proposition 2 for a slightly different parameterization).

4.2 Estimation Procedure↩︎

The doubly robust characterization in Proposition 1 motivates a two-stage procedure to estimate \(\theta(a)\): In the first step we model the nuisance functions that appear in the pseudo-outcome \(\varphi(\mathbf{Z})\) with flexible nonparametric or machine learning methods. We then construct estimated pseudo-outcomes \(\widehat{\varphi}(\mathbf{Z})\) and regress these on \(A\) to obtain an estimator of \(\theta(a)\). We will use the stability framework developed in [23] to analyze such an estimator, which regresses imputed outcomes \(\widehat{\varphi}(\mathbf{Z})\) on treatment. The formal procedure is summarized in Algorithm 2.

Sample splitting is used in Algorithm 2 to avoid complicated empirical process assumptions that are difficult to justify in practice and simplify our theoretical analysis [24]–[28]. In Step 1, one can use any appropriate regression/classification algorithms to estimate \(\mu, \tau, \rho\). However, there are fewer results on estimating stabilized weight \(w(a,\mathbf{v}) = f(a)/\pi(a|\mathbf{v})\). [18] proposed a method that directly estimates \(w\) with entropy maximization. Alternatively, one can estimate the conditional density \(\pi\) [29], and then the marginal density can be estimated by \[\widehat{f}(a) = \frac{1}{n} \sum_{i \in D_2^n} \widehat{\pi}(a| \mathbf{V}_i)\] and use their ratio to estimate \(w\).

Our analysis in Section 4.3 applies to a wide class of nuisance estimators if the product of convergence rates is fast enough. In step 2 we focus on linear smoother-based estimators since they are relatively straightforward to analyze. We believe similar results hold for a wider class of regression estimators under the stability framework in [23] and leave the theoretical analysis of applying general regression algorithms for future investigation. In applications, researchers can choose suitable parametric methods based on their domain knowledge or flexible nonparametric machine learning methods to avoid model misspecification (or ensembles thereof).

4.3 Theoretical Results↩︎

We first present estimation error guarantees of Algorithm 2, for general linear smoothers. Then we focus on a specific type of linear smoother, namely local linear regression, and study the asymptotic distribution of the estimator in Section 4.3.2.

4.3.1 Oracle Estimation Theory↩︎

In Algorithm 2 we obtain the estimator \(\widehat{\theta}(a)\) by regressing the imputed pseudo-outcome \(\widehat{\varphi}(\mathbf{Z})\) on \(A\). It is natural to compare \(\widehat{\theta}(a)\) with the “oracle" estimator \(\widetilde{\theta}(a)\) defined as \[\widetilde{\theta}(a) = \sum_{i \in T^n}W_i(a;\mathbf{A}^n)\varphi(\mathbf{Z}_i),\] which regresses the ground-truth pseudo-outcome \(\varphi(\mathbf{Z})\) on \(A\) using the same linear smoother. Intuitively it is hard for \(\widehat{\theta}\) to have a faster convergence rate than \(\widetilde{\theta}\). In the following theorem, we summarize the conditions under which \(\widehat{\theta}\) enjoys the same rate as \(\widetilde{\theta}\) and hence is”oracle efficient".

The conditions in Theorem 2 are mild: \(\varphi(\mathbf{Z})\) is not a function of \(A\) and hence its conditional variance given \(A\) should be positive; We do not impose conditions on the convergence rate of \(\widehat{\varphi}\) but only require its consistency. The key condition for \(\widehat{\theta}\) to be oracle efficient is \(\widehat{\mathbb{E}}_n\left[\widehat{b}(A)|A=a\right] = o_{\mathbb{P}}(R_n(a))\) and hence we need a bound on the conditional bias \(\widehat{b}(a)\), as summarized in the following proposition.

We note that the condition on linear smoothers will be satisfied by Nadaraya–Watson estimators [30], [31] and local polynomial estimators [32], [33] when the kernel function used has compact support, e.g., the uniform and Epanechnikov kernel. In the bound for \(\widehat{b}(a)\), the first two terms depend on the product of the convergence rates of nuisance functions. This phenomenon is commonly observed in influence functions-based doubly robust approaches [25], [34], [35] and gives the second interpretation of double robustness: Introducing an extra term to the plug-in style estimator in \(\varphi\) can correct for the first-order bias, making the remainder second-order and “doubly small". In common examples of ATE estimation, conditional bias only involves estimation error of outcome model and propensity scores. The additional dependence on the convergence rate of \(\widehat{\tau}\) (compared with Proposition 1) appears in the bias since the formula ?? is more complicated compared with ATE-style functional and we use an agnostic estimator of \(\tau\) in Algorithm 2. In most settings we expect the nuisance error rate \(\|\widehat{\mu}-\mu\|_a\) (with respect to measure \(d \mathbb{P}(\mathbf{z}|A=a)\)) has the same order as the more common conventional rate \(\|\widehat{\mu}-\mu\|\) (with respect to measure \(d \mathbb{P}(\mathbf{z})\)), which is \(n^{-\alpha/(2\alpha+d+1)}\) if \(\mu(a,\mathbf{x})\) belongs to a   class of order \(\alpha\) and \(\mathbf{X}\) is \(d\)-dimensional. Alternatively one can always upper bound \(\|\widehat{\mu}-\mu\|_{a}\) with \(\|\widehat{\mu}-\mu\|_{\infty}\) at the cost of log factors [36]. The empirical process term \(\frac{1}{n} \sum_{i \in D_2^n} \widehat{\tau}(a,\mathbf{V}_i) - \mathbb{P}[\widehat{\tau}(a, \mathbf{V})]\) would be \(O_{\mathbb{P}}(1/\sqrt{n})\) provided that \(\mathbb{E}[\widehat{\tau}^2(a,\mathbf{V})|D^n]\) is bounded. The bound on \(\widehat{\mathbb{E}}_n\left[\widehat{b}(A)|A=a\right]\) in Proposition 2 involves \[\sup_{t \in N(a)} \left | \frac{1}{n}\sum_{i \in D_2^n} \widehat{\tau}(t,\mathbf{V}_i) - \mathbb{P}[\widehat{\tau}(t,\mathbf{V})] \right|\] as a coarse bound and for specific estimators, a tighter bound can be derived. For instance, as we will see in local linear estimation this empirical process term is \(o_{\mathbb{P}}\left(1/\sqrt{nh} \right)\) and asymptotically negligible. Hence we can focus on the first two second-order terms in the conditional bias. Theorem 2 together with Proposition 2 gives conditions under which the estimator \(\widehat{\theta}(a)\) has the same rate as the oracle estimator \(\widetilde{\theta}(a)\). For instance, assume \(\mathbf{V}= \mathbf{X}\) (no surrogates) and \(|\mathbf{X}| = d\), suppose the dose-response function \(\theta(a)\) is \(\alpha\)-smooth (i.e. belongs to a   class of order \(\alpha\)) and \(w\) and \(\mu\) are \(s\)-smooth, then the rate condition for \(\widehat{\theta}(a)\) to be oracle efficient is \[n^{-\frac{2s}{2s+d+1}} \leq n^{-\frac{\alpha}{2\alpha+1}}\] or equivalently \[s \geq \frac{\alpha(d+1)}{2(\alpha+1)}.\]

4.3.2 Asymptotic Normality↩︎

In the following discussions, we will analyze the estimator \(\widehat{ \theta}(a)\) based on a particular linear smoother (i.e. the local linear regression estimator). For a scalar bandwidth parapmeter \(h\), let \(\mathbf{g}_{ha}(t) = [1, (t-a)/h]^{\top}\) be the local linear basis, \(K_{ha}(t) = h^{-1}K((t-a)/h)\) with \(K\) being a probability density. The local linear regression solves the following weighted least square problem: \[\min_{\boldsymbol{\beta}\in \mathbb{R}^2} \sum_{i \in T^n} K_{ha}(A_i) \left[ \widehat{\varphi}(\mathbf{Z}_i) -\mathbf{g}_{ha}(A_i)^{\top}\boldsymbol{\beta}\right]^2\] which gives the closed-form solution: \[\widehat{\boldsymbol{\beta}}_h(a) = \widehat{\mathbf{D}}_{ha}^{-1} \mathbb{P}_n [\mathbf{g}_{ha}(A)K_{ha}(A)\widehat{\varphi}(\mathbf{Z})],\] where \(\widehat{\mathbf{D}}_{ha} = \mathbb{P}_n[\mathbf{g}_{ha}(A)K_{ha}(A)\mathbf{g}_{ha}^{\top}(A)]\) and \(\mathbb{P}_n\) is the sample average over \(T^n\). Then the local linear estimator of \(\theta(a)\) is \(\mathbf{e}_1^{\top}\widehat{\boldsymbol{\beta}}_h(a)\), i.e. the first component of \(\widehat{\boldsymbol{\beta}}_h(a)\). We summarize the asymptotic properties of this local linear estimator in Theorem 3, which is the key contribution of our paper.

Assumptions 1-4 in Theorem 3 are standard in the nonparametric kernel regression literature. Assumption 5 guarantees that the contribution of nuisance estimation error is asymptotically negligible compared with the smoothing error. One can also use a symmetric kernel \(K\) supported on \(\mathbb{R}\) that is square-integrable and has a finite second-order moment. Then the rate condition would be \(\sup_{t \in \mathcal{A}} \|\widehat{\rho} - \rho\|_t \|\widehat{\mu} - \mu\|_t = o_{\mathbb{P}} \left(\frac{1}{\sqrt{nh}} \right)\). Similar to our discussions in Section 4.3.1, the rate conditions in assumption 5 are imposed on the product of nuisance estimation error since we use a doubly robust estimator and the conditional bias is “second-order small". The theoretically optimal bandwidth to estimate a twice continuously differentiable function is \(h \asymp n^{-1/5}\) and yields a root mean square error of order \(n^{-2/5}\). With such a choice of \(h\) the requirement on the convergence rate becomes \(\sup_{t \in \mathcal{A}} \|\widehat{\rho} - \rho\|_t \|\widehat{\mu} - \mu\|_t = o_{\mathbb{P}} \left(n^{-2/5} \right)\), which, for example, can be satisfied when \(\rho\) is consistent and \(\mu\) is estimated with rate \(O_{\mathbb{P}}(n^{-2/5})\). In applications, one can select the bandwidth using leave-one-out cross-validation [37] due to its computational ease. Specifically, after we obtain the estimated pseudo-outcome, we treat them as known and select \(h\) by \[\label{eq:loocv} \widehat{h}_{opt} = \mathop{\mathrm{argmin}}_{h \in \mathcal{H}} \sum_{i } \left\{ \frac{\widehat{\varphi}(\mathbf{Z}_i)-\widehat{\theta}_h(A_i)}{1-\widehat{W}_h(A_i)} \right\}^2.\tag{3}\] where \(\widehat{W}_h(A_i) = \mathbf{e}_1^{\top} \widehat{\mathbf{D}}_{hA_i}^{-1}\mathbf{e}_1 h^{-1} K(0)\). In the setting of Algorithm \(\ref{alg:DR-surrogates}\) we can select the bandwidth on \(D_2^n\) to avoid overfitting on \(T^n\).

Similar to most nonparametric inference methods, Theorem 3 shows the estimator \(\widehat{\theta}\) is centered around \(\bar{\theta}(a) = \theta(a) - h^2 \theta''(a) \int u^2 K(u)du /2\) instead of \(\theta(a)\) under optimal smoothing, which is known as the “bias problem" in the literature [38]. There are several methods to overcome the bias problem and each of them has its own consideration and trade-offs. For instance, one can estimate the second-order derivative and debias the estimator [39], [40] but this requires extra smoothness conditions. Another method is to undersmooth [41] and make the bias decrease asymptotically relative to the variance. Unfortunately, there does not seem to be a simple, practical rule for choosing just the right amount of undersmoothing. In this paper we choose to live with the bias and report uncertainty quantification for \(\bar{\theta}\). Theoretically, the bias shrinks to 0 as \(n \rightarrow \infty\) and the proposed estimator \(\widehat{\theta}(a)\) is still consistent for \(\theta(a)\). To construct confidence intervals for \(\bar{\theta}(a)\) one needs to estimate the variance. Define a localized functional \(\theta_h(a) = \mathbf{e}_1^{\top} \mathbf{D}_{ha}^{-1}\mathbb{E}[\mathbf{g}_{ha}(A)K_{ha}(A)\theta(A)]\) (which can be viewed as population version of local linear estimator \(\mathbf{e}_1^{\top}\widehat{\boldsymbol{\beta}}_h(a)\)) with efficient influence function \[\begin{align} \phi_{ha}(\mathbf{Z}) = &\,\mathbf{e}_1^{\top}\mathbf{D}_{ha}^{-1}\mathbf{g}_{ha}(A)K_{ha}(A)\left(\varphi(\mathbf{Z})-\mathbf{g}_{ha}^{\top}(A)\mathbf{D}_{ha}^{-1}\mathbb{E}[\mathbf{g}_{ha}(A)K_{ha}(A)\theta(A)]\right) \\ &\,+ \mathbf{e}_1^{\top}\mathbf{D}_{ha}^{-1}\int \mathbf{g}_{ha}(t)K_{ha}(t)\tau(t,\mathbf{V})f(t) dt -\theta_h(a). \end{align}\] Following [13], [40], one can show the variance of \(\widehat{\theta}(a)\) can be approximated by \(\frac{1}{n} \mathbb{P}_n \left[\left(\widehat{\phi}_{ha}(\mathbf{Z})\right)^2\right]\) for \[\widehat{\phi}_{ha}(\mathbf{Z})=\mathbf{e}_1^{\top}\widehat{\mathbf{D}}_{ha}^{-1}\mathbf{g}_{ha}(A)K_{ha}(A)\left(\widehat{\varphi}(\mathbf{Z})-\mathbf{g}_{ha}^{\top}(A)\widehat{\beta}_h(a)\right)+ \mathbf{e}_1^{\top}\widehat{\mathbf{D}}_{ha}^{-1}\int \mathbf{g}_{ha}(t)K_{ha}(t)\widehat{\tau}(t,\mathbf{V}) d \mathbb{P}_n(t) -\widehat{\theta}(a).\]

Finally, we compare the asymptotic variance in Theorem 3 with the asymptotic variance in [13], where the unlabeled dataset \(\mathcal{U}\) and surrogates \(\mathbf{S}\) are unavailable. Consider the MCAR setting where \(R\) is independent of all other variables so that \(\rho(a,\mathbf{x}) = \rho \in (0,1)\), and for simplicity assume \(n_1 = n\rho\) to show how the surrogates and unlabeled data help to reduce the variance in this special setting. Note that since \(R\) is independent of all other variables, we have \(\operatorname{Var}(Y|A=a,\mathbf{X},R=1) = \operatorname{Var}(Y|A=a,\mathbf{X})\) and \(\mu(a,\mathbf{X},1) = \mathbb{E}[Y|A=a,\mathbf{X},R=1] = \mathbb{E}[Y|A=a,\mathbf{X}] = \mu(a,\mathbf{X})\). The asymptotic variance of \(\widehat{\theta}(a)\) in our setting (i.e. in Theorem 3) is reduced to \[\label{eq:surrogates-var} \frac{1}{nh} \mathbb{E}\left\{ \left[ \frac{\operatorname{Var}(Y|A=a,\mathbf{X})}{\rho}+\operatorname{Var}(\mu(a,\mathbf{X})|A=a, \mathbf{V}) \right]w^2(a,\mathbf{V}) \bigg | A=a \right\}\tag{4}\] under the MCAR assumption (note the additional factor \(\int K^2(u) du /f(a)\) is omitted since it appears in both settings). In the setting where the unlabeled data is unavailable, the asymptotic variance is shown to be \[\label{eq:supervised-var1} \frac{1}{n_1 h}\mathbb{E}\left\{ \operatorname{Var}(Y|A=a,\mathbf{V})w^2(a,\mathbf{V})|A=a \right \}\tag{5}\] in [13]. By the property of conditional variance \[\operatorname{Var}(Y|A=a,\mathbf{V}) = \mathbb{E}[\operatorname{Var}(Y|A=a,\mathbf{X})|A=a, \mathbf{V}]+\operatorname{Var}(\mathbb{E}[Y|A=a,\mathbf{X}]|A=a,\mathbf{V}),\] 5 can be re-written as \[\label{eq:supervised-var2} \frac{1}{n_1 h}\mathbb{E}\left\{ (\operatorname{Var}(Y|A=a,\mathbf{X})+\operatorname{Var}(\mu(a,\mathbf{X})|A=a,\mathbf{V}))w^2(a,\mathbf{V})|A=a \right \} .\tag{6}\] Comparing 4 with 6 we see the first term is the same since \(n\rho = n_1\). However the second term in 6 is improved by a factor of \(\rho\) in 4 . This shows how the variance of the estimator is smaller after introducing unlabeled data and surrogate outcomes. The amount of improvement depends on the missing rate \(1-\rho\) and \(\operatorname{Var}(\mu(a,\mathbf{X})|A=a,\mathbf{V})\), which measures the variation of \(\mu(A,\mathbf{V}, \mathbf{S})\) that cannot be explained by \((A,\mathbf{V})\).

5 Simulation Study↩︎

In this section we use simulations to evaluate the performance of the proposed methods. We will illustrate the advantage of doubly robust estimation over naive plug-in style estimators. Consider the following data-generating process: The covariates \(\mathbf{V}\) have a multi-variate Gaussian distribution \[\mathbf{V}= (V_1, V_2, V_3, V_4)^{\top} \sim N(\mathbf{0}, \mathbf{I}_4) .\] Conditioning on \(\mathbf{V}\), the continuous treatment \(A\) has normal distribution \(N(\lambda(\mathbf{V}),1)\) with \[\lambda(\mathbf{V}) = 1+0.2V_1+0.2V_2 - 0.2V_3 + 0.3V_4.\] The surrogates \(\mathbf{S}\) has a normal distribution \[\mathbf{S}= (S_1,S_2)^{\top} \sim N(\mathbf{0}, \mathbf{I}_2).\] The indicator of whether the outcome is observed or not \(R\) is Bernoulli\((0.5)\) (so we assume a missing completely at random mechanism and \(\rho = 0.5\)). Finally the outcome \(Y\) conditioning on \(A,\mathbf{X},R=1\) has a normal distribution \(N(\mu(A,\mathbf{X},1), 1)\) with \[\mu(A,\mathbf{X},1) =1 + (0.1,-0.1)^{\top} \mathbf{S}+ (0.2,0.2,0.3,-0.1)^{\top}\mathbf{V}+ A(1-0.1V_1 + 0.1V_3) - A^2.\] By direct calculations we have \[\tau(A,\mathbf{V}) =1+ (0.2,0.2,0.3,-0.1)^{\top}\mathbf{V}+ A(1-0.1V_1 + 0.3V_3) - A^2,\] The dose-response function of interest is \[\theta(A) = 1 + A-A^2.\] To illustrate the performance of two estimators with different nuisance estimation errors we will manually set the estimation error, which is applicable for simulation purposes. For a fixed \(\alpha\) we let \(\epsilon_1, \dots, \epsilon_4 \sim N(n^{-\alpha}, n^{-2\alpha})\) and set \(\widehat{\lambda}(\mathbf{V})= \lambda(\mathbf{V}) + \epsilon_1\), the estimated conditional density of \(A\) is \(N\left(\widehat{\lambda}(\mathbf{V}),1 \right), \widehat{\mu}(A,\mathbf{X},1) = \mu(A,\mathbf{X},1) + \epsilon_2, \widehat{\tau}(A,\mathbf{V}) = \tau(A,\mathbf{V}) + \epsilon_3, \text{logit} \left(\widehat{\rho} \right) = \text{logit}(\rho) + \epsilon_4\). Such estimates guarantee the nuisance estimation error is of order \(n^{-\alpha}\). After we generate a sample of size \(n\), the plug-in style estimator is defined as \(\widehat{\theta}(a) = \frac{1}{n} \sum_{i=1}^n \widehat{\tau}(a,\mathbf{V}_i)\). To implement the doubly robust estimator we split the sample into two parts \(D,T\) (since the nuisance estimators are given there is no need to split the sample into three parts as in Algorithm 2). We use \(D\) to construct estimator of the initial estimator \(\widehat{\theta}_0(a) = \frac{1}{|D|}\sum_{i \in D}\widehat{\tau}(a,\mathbf{V}_i)\), marginal density \(\widehat{f}(a) = \frac{1}{|D|}\sum_{i \in D}\widehat{\pi}(a|\mathbf{V}_i)\) and select the optimal bandwidth \(h^*\) (i.e. construct pseudo-outcomes and run LOOCV as in equation 3 on \(D\)). We then construct pseudo-outcomes on \(T\) and perform local linear estimation using the optimal bandwidth \(h^*\). Finally, the roles of \(D\) and \(T\) are exchanged to obtain another estimator and we average the two estimates as the final doubly robust estimator. For sample size \(n \in \{500, 2000\}\) and convergence rate \(\alpha \in \{ 0.1, 0.13, \dots, 0.4\}\), we repeat the data generation and estimation process \(M=500\) times. We will aim at estimating \(\theta(1)\) and compare the RMSE\(=\left \{\frac{1}{M} \sum_{m=1}^M \left[\widehat{\theta}_m(1)-\theta(1)\right]^2 \right \}^{1/2}\) of plug-in estimator and doubly robust estimator, where \(\widehat{\theta}_m(1)\) is the estimate from \(m\)-th repetition. The results are summarized in Figure 3.

As we see in Figure 3, if the nuisance estimation error is large (\(\alpha\) is small), the doubly robust estimator outperforms the naive plug-in estimator. This can be explained by the second-order bias term of the doubly robust estimator in Proposition 2, i.e., the conditional bias is the product of nuisance estimation errors and is “doubly small". On the other hand, the plug-in style estimator inherits the slow convergence rate of \(\widehat{\tau}\). As \(\alpha\) increases, the estimators of nuisance functions are more accurate, and plug-in style estimators finally outperform the doubly robust estimator because the doubly robust estimator may suffer from accumulating error in constructing pseudo-outcomes, bandwidth selection, and local linear regression, which dominates the conditional bias when nuisance estimation is accurate enough. In real applications, there are typically many covariates and parametric models for nuisance functions may not be correct. The convergence rate of nonparametric nuisance estimation can be slow when the dimension of covariates is large so the doubly robust estimator with a smaller bias is recommended for use. Additional simulation results and a real data example are provided in the supplementary materials.

6 Discussion↩︎

In this work, we study the estimation of continuous treatment effects when there is limited access to the primary outcome of interest but auxiliary information on surrogate outcomes is available. We propose a doubly robust estimator that is less sensitive to nuisance estimation error and hence incorporates flexible nonparametric machine learning methods. Although nonparametric machine learning methods usually suffer from slow convergence rates, they are widely used in nuisance function estimation, especially when practitioners do not have sufficient domain knowledge to justify parametric models. Our doubly robust estimator facilitates the application of nonparametric methods and enjoys optimal estimation rates under mild conditions. Asymptotic normality is further established, which enables researchers to construct confidence intervals and perform statistical inference. We also show how incorporating information on surrogate outcomes improves the variance, compared with methods solely based on labeled data, as in [13]. In summary, our methodology provides a robust and efficient approach to leverage surrogate outcomes in continuous treatment effect estimation. However, our method could not deal with the case where only covariate is available in the unlabeled dataset and continuous treatment is also missing (corresponds to the generalizability problem as in [42]). It is also interesting to extend our results to studies with multiple outcomes [43], [44]. We leave these problems for future investigation.

7 Background on Efficiency Theory↩︎

As discussed in Section 3, the plug-in estimator suffers from first-order error and is sensitive to the estimation accuracy of nuisance functions. To address this problem, one can derive the efficient influence function of the target functional (\(\mathbb{E}[\theta(A)]\) in our problem), based on which a one-step estimator can be obtained to reduce the bias. The efficient influence function is critical in non-parametric efficiency theory [20], [45]–[48]. Mathematically, the influence function is the derivative of the target statistical functional in a Von Mises expansion (i.e., distributional Taylor expansion). In the discrete case, it coincides with the Gateaux derivative of the functional when the contamination distribution is a point mass. Influence functions are important in different respects. First, the variance of the influence function is equal to the efficiency bound of the target statistical functional, which characterizes the inherent estimation difficulty of the target functional and provides a benchmark to compare against when we construct estimators. Moreover, it allows us to correct for first-order bias in the plug-in estimator and obtain doubly robust-style estimators, which enjoy appealing statistical properties even if non-parametric methods with relatively slow rates are used in nuisance estimation.

Suppose the statistical functional of interest \(\psi=\psi(\mathbb{P})\) admits the first-order Von Mises expansion. Mathematically, we have \[\label{eq:von-mises} \psi(\widehat{\mathbb{P}})-\psi(\mathbb{P})=-\int \phi_1(\mathbf{z},\widehat{\mathbb{P}}) d \mathbb{P}(\mathbf{z})+R_2(\widehat{\mathbb{P}}, \mathbb{P}),\tag{7}\] where \(\phi_1(\mathbf{z},\mathbb{P})\) is the influence function of \(\psi(\mathbb{P})\), \(\psi(\widehat{\mathbb{P}})\) is the plug-in estimator and \(R_2(\widehat{\mathbb{P}}, \mathbb{P})\) is the second-order reminder. Under regularity conditions, Von Mises expansion implies the pathwise differentiability \[\left.\frac{\partial}{\partial \epsilon} \psi\left(\mathbb{P}_\epsilon\right)\right|_{\epsilon=0}=\int \phi_1(\mathbf{z}; \mathbb{P}) s_\epsilon(\mathbf{z}) d \mathbb{P}(z)\] where \(s_\epsilon(\mathbf{z})=\left.\frac{\partial}{\partial \epsilon} \log p_\epsilon(\mathbf{z})\right|_{\epsilon=0}\) is the submodel score function. 7 suggests that the plug-in estimator has first-order bias \(-\int \phi_1(\mathbf{z},\widehat{\mathbb{P}}) d \mathbb{P}(\mathbf{z})\). Equivalently, we can write \[\psi(\mathbb{P})=\psi(\widehat{\mathbb{P}})+\int \phi_1(\mathbf{z},\widehat{\mathbb{P}}) d \mathbb{P}(\mathbf{z}) +R_2(\widehat{\mathbb{P}}, \mathbb{P}),\] which motivates us to correct for the first-order bias and arrive at the doubly robust estimator \[\widehat{\psi}^{dr} = \psi(\widehat{\mathbb{P}})+\mathbb{P}_n\{\phi_1(\mathbf{Z},\widehat{\mathbb{P}})\}.\] Under further empirical process assumptions or sample splitting assumptions, we can show the dominating term in the conditional bias of \(\widehat{\psi}^{dr}\) is \(R_2(\widehat{\mathbb{P}}, \mathbb{P})\), which is usually a second-order error term and depends on the product of convergence rates of nuisance functions. We refer the readers to [34], [48] for a more complete review. In our problem, the estimand \(\theta(a)\) in ?? is not pathwise differentiable and the ideas above do not apply directly. However, the pseudo-outcome \(\varphi(\mathbf{Z})\) is obtained by deriving the influence function of \(\mathbb{E}[\theta(A)]\). The readers are referred to Section 3.1 of [13] for more discussion on the connection between the pseudo-outcome of \(\theta(a)\) and the influence function of \(\mathbb{E}[\theta(A)]\).

8 Proofs↩︎

8.1 Proof of Theorem 1↩︎

8.2 Proof of Proposition 1↩︎

8.3 Proof of Theorem 2↩︎

8.4 Proof of Proposition 2↩︎

8.5 Proof of Theorem 3↩︎

9 Additional Simulation Results↩︎

9.1 Nuisance Functions Estimated by Parametric Methods↩︎

We evaluate the performance of plug-in-style estimator and doubly robust estimator when nuisance functions are estimated using parametric models under the same setting as Section 5. We will fit linear regression models for \(\mu, \tau, \lambda\) and separately consider correctly specifying the outcome model \(\mu, \tau\) or not, where a misspecified model left out the quadratic term in \(a\) but keeps all other main effects and interactions. The conditional density of \(A\) given \(\mathbf{V}\) is obtained by first estimating the conditional mean of \(\mathbb{E}[A|\mathbf{V}]\) and plug-in the normal density (i.e., we assume the model for conditional density is always correct). For the plug-in estimator we randomly separate the sample into two parts \(D,T\). The first part \(D\) is used to fit the regression models for all nuisance functions and we take the average on the second part as \(\widehat{\theta}(a) = \frac{1}{|T|}\sum_{i \in T} \widehat{\tau}(a, \mathbf{V}_i)\). The roles of \(D,T\) are then exchanged to obtain another estimate and the final estimator is the average of two estimates. The doubly robust estimator is implemented according to Algorithm 2, where the bandwidth is selected on \(D_2\). We generate samples with sample size \(n \in \{ 10^{2.6}, 10^{2.8}, \dots, 10^{4.6}\}\), apply two estimators to estimate \({\theta}(1)\) and repeat the process \(500\) times. The results are summarized in Figure 4.

When the outcome model is misspecified, the plug-in estimator (solely based on outcome modeling) is no longer consistent, as shown in Figure 4. The doubly robust estimator, however, models both the outcome and treatment process and has a smaller estimation error when the outcome model is misspecified. Estimation with correctly specified outcome model corresponds to \(\alpha = 0.5\) in Figure 3, where doubly robust estimator with outcome model and propensity score both correctly specified has larger error compared with the plug-in estimator since slow rate of local linear smoothing \(O(n^{-2/5})\) dominates the nuisance estimation error \(O(n^{-1/2})\).

9.2 Nuisance Functions Estimated by Nonparametric Methods↩︎

We further evaluate the performance of the plug-in style estimator and doubly robust estimator when nuisance functions are estimated using nonparametric models under the same setting as Section 5. We will fit nuisance functions \(\mu, \tau, \rho\) by superlearner combining generalized linear models and random forests. The conditional density of \(A\) given \(\mathbf{V}\) is estimated by kernel density estimation. The plug-in estimator is implemented the same way as in Appendix 9.1 and the doubly robust estimator is implemented according to Algorithm 2. We generate samples with sample size \(n \in \{ 10^{2.6}, 10^{2.8}, \dots, 10^{4}\}\), apply two estimators to estimate \({\theta}(1)\) and repeat the process \(500\) times. The results are summarized in Figure 4.

As shown in Figure 5, the plug-in-style estimator has smaller estimation error when all the nuisance functions are fitted by nonparametric methods. This could be explained by the diffculty in nonparametrically estimating conditional density \(\pi\). Due to the curse of dimensionality, a large sample size is required for the kernel density estimator to estimate \(\pi\) well. In our simulations, we find \(\widehat{\pi}\) could be small and hence violate the positivity assumption, yielding variation in the construction of pseudo-outcome \(\widehat{\varphi}\) and larger estimation error of dose-response compared with the plug-in estimator. In applications, prior knowledge on the conditional density, including a reliable parametric model on \(\pi\) from domain knowledge or information on the lower bound on \(\pi\) due to study design might help us reduce the variation of \(\widehat{\pi}\) and estimate the conditional density better, which could improve the performance of doubly robust estimator proposed.

9.3 Surrogates Dependent on Treatment and Covariates↩︎

In this section we provide additional simulation results on the setting where surrogates may depend on the treatment and covariates. In real applications, the surrogates \(\mathbf{S}\) are post-treatment short-term outcomes and are likely to depend on pre-treatment covariates and the treatment. We follow the same setting and estimation procedures as in Section 5 (so the nuisance estimation error is set manually) but modify the conditional distribution of \(\mathbf{S}\) given \((A,\mathbf{V})\) as \[\mathbf{S}\sim N\left((V_1+A,V_2-A)^{\top}, \boldsymbol{I}_2\right).\] The RMSE is estimated from \(M=500\) replications of the data-generating process and estimation procedures. The results of DR-learner and plug-in estimators are summarized in Figure 6.

Figure 6 shows when surrogate outcomes \(\mathbf{S}\) depend on treatment and covariates, the doubly robust estimator enjoys better estimation accuracy compared with the plug-in estimator when \(\alpha\) is small (and hence nuisance estimation error is large). As \(\alpha\) increases and the nuisance estimation error becomes smaller, the plug-in estimator gradually outperforms the doubly robust estimator. The results are similar to those in Section 5 and readers are referred to Section 5 for more discussion and explanation.

9.4 Comparison Between Supervised and Semi-supervised Methods↩︎

In this section we compare the supervised estimator, which implements the method in [13] on labeled data \(\mathcal{L}\), with our semi-supervised estimator incorporating the unlabeled data \(\mathcal{U}\) with surrogate outcomes. We follow the same setting in Section 5 and implement method in [13] via similar cross-fitting techniques as in Algorithm 2 except that there is no need to fit \(\rho(a,\mathbf{x})= \mathbb{P}(R=1 \mid A=a, \mathbf{X}= \mathbf{x})\) since their method only uses labeled data. All the nuisance functions are estimated by correctly specified parametric models. We generate samples with sample size \(n \in \{ 10^{2.6}, 10^{2.8}, \dots, 10^{4.6}\}\), apply supervised and semi-supervised methods to estimate \({\theta}(1)\) and repeat the process \(M=500\) times. The results are summarized in Figure 7.

Figure 7 shows the semi-supervised method incorporating unlabeled data with surrogate outcomes has a smaller estimation error compared with the supervised method only using labeled data. The improvement comes from many aspects: As discussed in Section 4.3.2, the asymptotic variance of our semi-supervised method is smaller than that of the supervised method only using labeled data; Also in the semi-supervised setting, both labeled and unlabeled data are used to estimate the nuisance functions, making nuisance estimation more accurate than supervised method solely based on labeled data since the effective sample size of semi-supervised method is larger.

10 Real Data Analysis↩︎

In this section we apply the proposed method to the Job Corps study, conducted in the 1990s to evaluate the effects of the publicly funded U.S. Job Corps program. The Job Corps program targets a population between 16 and 24 years old living in the U.S. and coming from low-income households, where participants received an average of 1,200 hours of vocational training over approximately eight months. [49] and [50] discuss the study design in detail and analyze the effects of the Jobs Corps program on various outcomes. They found the Job Corps program effectively increased educational attainment, prevented arrests, and increased employment and earnings. The effects of the Job Corps program have been extensively studied under different causal inference frameworks [51]–[53].

However, these previous studies on the Job Corps program mainly considered binary treatment definitions. In this work, we are interested in estimating the effects of different doses of participation in the program on future involvement in the criminal justice system, namely the number of arrests in the fourth year after the program (outcome \(Y\)). Specifically, our treatment variable \(A\) is defined as the total hours spent either in academic or vocational classes of the program. The short-term surrogate outcome \(S\) is the proportion of weeks employed in the second year after the program. [54] used generalized propensity score weighting to estimate the continuous treatment effects of time spent in the Job Corps program under a mediation analysis framework. We re-analyze their dataset publically available on Harvard dataverse [55], using the doubly robust estimator proposed as an illustration of our method.

We follow [54] and focus on \(n=4000\) samples with a positive treatment (i.e., \(A_i>0, 1 \leq i \leq n\)). To identify the causal estimand, we invoke the conditional exchangeability in Section 3, where a set of covariates is conditioned on to adjust for confounding bias. The covariate set \(\mathbf{V}\) we adjust for confounding consists of age, gender, ethnicity, education, marital status, previous employment status and income, welfare receipt during childhood, and family background (e.g., parents’ education). Missing dummies are created for covariates in \(\mathbf{V}\) containing missing values. The readers are referred to Table 4 in [54] for descriptive statistics of the pretreatment covariates as well as the treatment, surrogate, and outcome variables in the data. Conditioning on a rich set of covariates is important since it enables us to identify the causal estimand by making the conditional exchangeability assumption plausible. In our analysis, the outcome \(Y\) (number of arrests in year 4) is omitted for 25% of the samples randomly to mimic the setting where the primary outcome is missing and a surrogate outcome is used as auxiliary information. We estimate the treatment effects \(\theta(a)\) for each of \(a \in \{100,150,200,\dots,2000\}\) using both plug-in-style estimator and doubly robust estimator in Algorithm 2, where the nuisance functions \(\rho, \mu, \tau\) are estimated by superlearner [56] combining generalized linear model and random forests, conditional density \(\pi\) is estimated by kernel density estimator. The estimated dose response curve is plotted in Figure 8.

Figure 8 shows that, as participants spend more hours in the program, the expected number of arrests in year 4 has a decreasing trend, which confirms the conclusion that such training programs effectively reduce involvement in the criminal justice system. Importantly, the dose-response fitted by the doubly robust estimator is very similar to the results in [54]: the shape of the function is similar to their Figure 2. (Note that the specific values on the y-axis are different since they plot a contrast effect and we instead plot the expectation of potential outcome \(Y^a\).) As pointed out in [54], the treatment effect is highly nonlinear, which is further verified by the curve estimated from the doubly robust estimator. By contrast, the curve estimated by a simple plug-in estimator in 2 fails to capture this non-linearity, possibly because it suffers from a large first-order bias.

References↩︎