High-dimensional Inference for Dynamic Treatment Effects

High-dimensional inference for dynamic treatment effects

Introduction

The complexity of a given disease or economic policy often manifests in the diversity and size of the personal characteristics pertaining to each individual or economy under consideration, causing a considerable degree of heterogeneity in observed outcomes. However, the utility of randomized control trials (RCTs), especially over time, is frequently curtailed by prohibitive costs or ethical concerns. In contrast, the accessibility of time-varying observational studies has burgeoned of late. The ubiquity of data-driven decision-making is evident in various aspects of daily life, such as the continuous monitoring of individuals' health using mobile devices and consequential medical interventions, tracking of online presence and real-time measurement of economic and social policies implemented to enhance public health. The present study contributes novel insights to the literature by proposing a novel framework to construct confidence intervals pertaining to dynamic treatment effects amid high-dimensional observations. In a Job Corps real-data analysis, our novel framework provides more accurate estimates of the long-term impact of additional schooling over time on wages, which has important practical implications for designing effective policies aimed at increasing educational attainment and improving economic outcomes.

In light of intricate notational complexities, we exemplify our ideas and findings for two-stage trials while affirming that the same theoretical framework and methodology developed are extensible to multiple-stage trials; see, e.g., Section (ref). Consider a two-stage series of binary treatment assignments, denoted by $A_1$ and $A_2$, and an outcome of interest, $Y \in \mathbb{R}$. Alongside this, a set of possibly high-dimensional sequential pre-treatment covariates $\mathbf{S}_1 \in \mathbb{R}^{d_1}$ and $\mathbf{S}_2 \in \mathbb{R}^{d_2}$, possibly of different dimensions, are also observed. The potential or counterfactual outcomes, $Y(a)$, refer to the outcome that a participant would have experienced had they followed a particular treatment sequence, $a = (a_1,a_2) \in \{0,1\}^2$, which may differ from the treatment they were observed with. Our parameter of interest is the dynamic treatment effect (DTE) between two treatment paths, $a$ and $a' $, which is defined as follows:

Estimating the DTE is a challenging task when there are multiple exposures involved. The influence of past treatments on future confounders and treatment choices complicates the identifiability of $\theta$ rosenbaum1983central. Adjusting for confounders may not have a causal interpretation, even when all confounders are measured and the regression is correctly specified daniel2013methods. In this context, alternative methods such as Sequential Multiple Randomized Control Trials (SMART) hernan2016specifying, Structural Nested Mean (SNM) Robins1997causal, and Marginal Structural Mean (MSM) models murphy2001marginal have become the gold standard for addressing these challenges. This paper contributes to the field by establishing robust MSM model estimations with new effective rates.

The doubly robust representations

Throughout this work, we assume that any treatment-specific variable can only be affected by past treatments or past covariates; and not the future. This is sometimes called temporal ordering. We also assume a “no interference” setting and Assumption (ref) below robins1987addendum, robins2000marginal,murphy2003optimal.

The following lemma provides a doubly robust (DR) representation of $\theta_a$. This result is consistent with previous studies in the literature, including works by van2012targeted,orellana2010dynamic,murphy2001marginal, bang2005doubly. We consider the MSM models where we adjust for confounding variables that may affect both the treatment assignment and the outcome of interest. In an MSM, the treatment assignment and the outcome of interest are modeled separately using propensity scores $\pi_{a} (\mathbf{s}_{1})$ and $\rho_{a}(\mathbf{s} )$ together with the first-time and second-time conditional means, $\mu_{a}(\mathbf{s}_{1}):= E[Y(a)|\mathbf{S}_{1}=\mathbf{s}_{1}]$ and $ \nu_{a}(\mathbf{s} ):= E[Y(a)|\mathbf{S} =\mathbf{s},A_1=a_1]$. Throughout this work, we use $\pi_{a}^* (\cdot)$ and $\rho_{a}^*(\cdot)$ as well as $\mu_a^*(\cdot)$ and $\nu_a^*(\cdot)$ to refer to the working models, i.e., the population-level approximations of the propensity scores and conditional means, respectively.

Based on Lemma (ref), consistent estimates of $\theta_{a}$ are expected as long as at least one nuisance model is correctly parametrized at each exposure time. However, this goal has not been achieved yet; see babino2019multiple for an overview. The main obstacle is the estimation of interlocking nuisance functions, especially the first-time conditional mean, as it cannot be identified directly through the observable variables as $\mu_{a}(\mathbf{s}_{1})=E[Y(a)|\mathbf{S}_{1}=\mathbf{s}_{1}]\neq E[Y|\mathbf{S}_{1}=\mathbf{s}_{1},A_1=a_1]$. Under Assumption (ref), existing DTE literature typically considers the following nested representation of $\mu_{a}(\cdot)$,

and suggests a nested regression (NR) of the conditional means -- as long as an estimate $\widehat{\nu}_a(\cdot)$ of $\nu_a(\cdot)$ is obtained, one can use $\widehat{\nu}_a(\mathbf{S}_i)$ as the imputed outcomes and perform regression to construct $\widehat{\mu}_{a,{\mbox{\tiny NR}}}(\cdot)$; see, e.g., murphy2001marginal. We formalize these properties under high-dimensional linear working models, naming the resulting DTE estimator the “dynamic treatment Lasso” (DTL) estimator. We show that the nested-regression approach faces certain limitations and fails to attain the DR property equivalent to Lemma (ref). Among the multiple factors contributing to this, the biggest one is arising from a peculiar model misspecification that we identified arising from the nested representation in Equation (ref). In the event of a misspecified linear working model $\nu_a^*(\cdot)$, the corresponding $\mu_a^*(\cdot)$ will inevitably be misspecified as well, leading to $\mu_a^*(\cdot) \neq \mu_a(\cdot)$, even when $\mu_a(\cdot)$ is itself linear. Besides the linearity of $\mu_a(\cdot)$, additional conditions on $\nu_a(\cdot)$ are necessary for the correctness of the nested-regression-based linear working model, as discussed in Section (ref).

This issue necessitates the use of specialized methods for which we propose a new DR representation of the first-time conditional mean function $\mu_a(\cdot)$; see (ref) below. It provides tools to quantify the DR property of the resulting DTE estimate and to develop correction techniques that can mitigate the DR gap by achieving the estimation under model conditions equivalent to Lemma (ref).

Utilizing the two DR representations (ref) and (ref) simultaneously, we propose a sequential doubly robust Lasso (S-DRL) estimator. The proposed estimator is consistent as long as either the conditional mean function is truly linear or the propensity score function is truly logistic (or both) for each exposure time. To the best of our knowledge, this is the first estimator that matches Lemma (ref) conditions empirically. The inverse probability weighting (IPW) methods robins1986new,robins2000marginal,hernan2001marginal,robins2004optimal require all the propensity score models to be correctly parametrized. The covariate balancing methods kallus2018optimal,yiu2018covariate,viviano2021dynamic require all the conditional mean models to be correctly parametrized. Perhaps unexpectedly, the standard low-dimensional DR methods robins2000robust,murphy2001marginal,bang2005doubly,yu2006double and the targeted maximum likelihood estimation (TMLE) van2012targeted require either all the propensity score functions or all the conditional mean (or density) functions to be correctly parametrized. The “multiply robust” (MR) estimator of babino2019multiple reaches better robustness than all of the aforementioned methods. In general $T$-stage trials, they allow for the first $t$ conditional mean models and the last $T-t$ propensity score models to be correctly parametrized for any $t$. The DTL estimator allows the first $t$ propensity score models and the last $T-t$ conditional mean models to be correctly parametrized. Our S-DRL estimator is strictly more robust in terms of consistency; see Table (ref) and Remark (ref) for further details.

The S-DRL estimator demonstrates superior estimation rates in high-dimensional contexts when compared to the DTL estimator; see Table (ref) as well as Remark (ref). Root-sample-size inference based on the S-DRL estimator is guaranteed when two product-sparsity conditions are satisfied, whereas the DTL method requires three product-sparsity conditions, as demonstrated in Theorems (ref) and (ref). The errors in nuisance estimation at different stages have a parallel effect on the estimation; see the consistency rate in Theorem (ref).

The estimation of the in-between outcome models is intrinsically linked to regression with imputed outcomes. We have developed a novel cone-set analysis of imputed Lasso estimates that is of independent interest to other imputed, high-dimensional regressions. Existing Lasso proof techniques provide conservative bounds only; see Section (ref). Our results are adaptive to the imputation error and can be used to guide the selection of tuning parameters in high-dimensional regression models with imputed outcomes.

In the multi-stage exposure setting, we extend our method and develop DR representations to identify both the expected potential outcomes and conditional means, as shown in Section (ref). While the consistency rate and asymptotic normality require intricate proofs, we anticipate they hold analogously to those in the two-stage case. It is worth noting that Theorem (ref) provides new DR representations that are independent of any specific parametric models, allowing the sequential doubly robust (S-DR) method to be utilized with non-parametric nuisance estimates, which enhances its versatility.

Organization of the paper

In Section (ref), we introduce the DR estimators of the DTE, including the proposed S-DRL estimator, the DTL estimator, and a general DR estimator. The theoretical properties of the considered DTE estimators are established in Section (ref). In Section (ref), we formalize the supporting theoretical discoveries, including a general theory for imputed Lasso estimation and the consistency results of the nuisance estimates. We further extend our setting to the case of multi-stage treatments and provide general DR representations for the intermediate conditional means in Section (ref). Section (ref) presents numerical results, including simulation studies and an application to the National Job Corps Study. Further discussion is provided in Section (ref).

Notation

For any $\alpha>0$, let $\psi_{\alpha}(\cdot)$ denote the function given by $\psi_{\alpha}(x):=\exp(x^\alpha)-1$, $\forall x>0$. Then the $\psi_{\alpha}$-Orlicz norm $\|\cdot\|_{\psi_{\alpha}}$ of a random variable $X$ is defined as $ \|X\|_{\psi_{\alpha}} :=\inf\{c>0:E[\psi_{\alpha}(|X|/c)]\leq 1\}. $ Two special cases of finite $\psi_{\alpha}-$Orlicz norm are given by $\psi_{2}(x)=\exp(x^2)-1$ and $\psi_{1}(x)=\exp(x)-1$, which correspond to sub-Gaussian and sub-exponential random variables, respectively. The notation $a_N\ll b_N$ denotes $a_N=o(b_N)$, and $a_N\gg b_N$ denotes $b_N\ll a_N$ as $N \to \infty$. The notation $a_N\asymp b_N$ denotes $cb_N\leq a_N \leq Cb_N$ for all $N \geq 1$ and with constants $c,C>0$. Define $g(u)= \exp(u)/\{1+\exp(u)\}$ as the logistic function and $\phi(u)=\log(1+\exp(u))$ as the corresponding link function throughout.

The doubly robust estimators

We observe a collection of independent and identically distributed (i.i.d.) samples $\mathcal D:=\{W_{i}\}_{i=1}^{N}=(Y_{i},\mathbf{S}_{1i},A_{1i},\mathbf{S}_{2i},A_{2i})_{i=1}^{N}$, drawn from the same distribution as $(Y,\mathbf{S}_1,A_1,\mathbf{S}_2,A_2)$. In the following subsections, we present three DTE estimators: the new sequential doubly robust Lasso (S-DRL) estimator, the dynamic treatment Lasso (DTL) estimator, and the general DR estimator.

The sequential doubly robust Lasso (S-DRL) estimator

We focus on the high-dimensional scenario, and consider linear (working) models for the conditional means $\mu_a(\cdot)$ and $\nu_a(\cdot)$, along with logistic (working) models for the propensities $\pi_a(\cdot)$ and $\rho_a(\cdot)$. The population minimizer approximating $\pi_{a} (\mathbf{s}_1)$ is defined as $\pi_{a}^{*}(\mathbf{s}_{1}) = g(\mathbf{v}^{\top}\bm{\gamma}_{a}^{*})$ with $\mathbf{v}=(1, \mathbf{s}_{1}^{\top})^{\top}$, whereas that of approximating $\rho_{a}(\mathbf{s} )$ is $ \rho_{a}^{*}(\mathbf{s} )=g(\mathbf{u}^{\top}\bm{\delta}_{a}^{*})$ with $\mathbf{u}=(1, \mathbf{s}^{\top})^{\top}$. Here

One can also consider a feature map $\varphi(\mathbf{s}_1)$ (e.g., a polynomial basis) and a working model $\pi_a^*(\mathbf{s}_1)=g(\varphi(\mathbf{s}_1)^\top\boldsymbol{\gamma}_a^*)$ with some $\boldsymbol{\gamma}_a^*$ defined correspondingly. We focus on $\varphi(\mathbf{s}_1)=\mathbf{v}$, although the results apply more broadly. The above working models can be estimated with many regularizations. Throughout this work, we focus on the $\ell_1$-regularization, albeit the theoretical developments apply more broadly. With a subset of training data $\mathcal D_\mathcal J=\{W_i\}_{i\in\mathcal J}\subset\mathcal D$, where $\mathcal J\subset\{1,\dots,N\}$, we define

with tuning parameters $ \lambda_{\bm{\gamma}},{\lambda}_{\bm{\delta}} >0$. Observe that for $ \widehat{\bm{\gamma}}_{a}$, we utilize all of the observations regardless of its treatment path, whereas for $ \widehat{\bm{\delta}}_{a}$, only those whose treatment path matches $a_1$ regardless of what $a_2$ is. The best linear working model for the second-time conditional mean $\nu_a(\cdot)=E[Y|\mathbf{S},A_1=a_1,A_2=a_2]$ is denoted as

An estimator of (ref) can be obtained similarly with $ {\lambda}_{\bm{\alpha}}> 0$:

\paragraph*{Estimation of the first-time conditional mean $\mu_a(\cdot)$}

Recall that $\boldsymbol{\delta}_a^*$ and $\boldsymbol{\alpha}_a^*$ are defined in Equations (ref) and (ref), respectively. We propose the following DR imputed outcome

With this in mind, we consider a linear working model for the first-time conditional mean

To estimate the best linear slope $\boldsymbol{\beta}_a^*$ based on a subset of training data $\mathcal D_\mathcal J\subset\mathcal D$, we consider an additional sample splitting with $\mathcal D_\mathcal J=\mathcal D_{\mathcal J_1}\cup\mathcal D_{\mathcal J_2}$, where $\mathcal J_1$ and $\mathcal J_2$ are disjoint subsets of $\mathcal J$. Using the first half of the subsamples $\mathcal D_{\mathcal J_1}$, we first obtain the second-time nuisance estimates $\widetilde{\boldsymbol{\delta}}_a:=\widehat{\boldsymbol{\delta}}_{a}(\mathcal D_{\mathcal J_1})$ and $\widetilde{\boldsymbol{\alpha}}_a:=\widehat{\boldsymbol{\alpha}}_a(\mathcal D_{\mathcal J_1})$ as (ref) and (ref), respectively. Then, for each $i\in\mathcal J_2$, we construct a DR imputed outcome

Based on the DR imputed outcomes $\widehat Y^{\mbox{\tiny DR}}_{\mathcal J_2}:=\{\widehat Y^{\mbox{\tiny DR}}_i\}_{i\in\mathcal J_2}$, we propose a DR estimate:

where $ {\lambda}_{\bm{\beta}}>0$. To regain full sample size efficiency, we can always swap the samples $\mathcal D_{\mathcal J_1}$ and $\mathcal D_{\mathcal J_2}$, repeat the procedure, and average the results.

\paragraph*{The S-DRL estimator of the DTE} For each $c\in\{a,a'\}$ and for any $\eta=(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\delta})$, define the DR score function based on the DR representation (ref):

We propose the sequential doubly robust Lasso (S-DRL) estimator of $\theta$:

where $\hat\eta_c:=(\widehat{\boldsymbol{\alpha}}_c,\widehat{\boldsymbol{\beta}}_c,\widehat{\boldsymbol{\gamma}}_c, \widehat{\boldsymbol{\delta}}_c)$ are the nuisance estimates of (ref), (ref), (ref), and (ref), respectively. A cross-fitting technique is used. The details are provided in Algorithm (ref); see chernozhukov2018double,smucler2019unifying where the cross-fitting leads to weaker sparsity restrictions than those without it, such as farrell2015robust,tan2020model.

The dynamic treatment Lasso (DTL) estimator

In this section, we formally define a dynamic treatment Lasso (DTL) estimator based on the DR score (ref) of $\theta_a$ and the nested representation (ref) of $\mu_a(\cdot)$. Here, $\ell_1$-regularized nuisance estimates $\widehat{\boldsymbol{\gamma}}_a$, $\widehat{\boldsymbol{\delta}}_a$, and $\widehat{\boldsymbol{\alpha}}_a$ will be the same as before; see (ref), (ref), and (ref) above. Estimation of the first-time conditional mean model is different. Based on the best linear approximation $\nu_{a}^{*}(\mathbf{s})=\mathbf{u}^{\top} \bm{\alpha}_{a}^{*}$, (ref), of $\nu_a(\cdot)$, we introduce the following nested “best linear working model”:

Note that the two linear working models $\mu_{a,{\mbox{\tiny NR}}}^*(\cdot)$ and $\mu_a^*(\cdot)$ are not necessarily the same; see Section (ref) for detailed comparisons. We consider the following imputed Lasso estimate of $\bm{\beta}_{a,{\mbox{\tiny NR}}}^{*}$, defined as $\widehat{\bm{\beta}}_{a,{\mbox{\tiny NR}}}:=\widehat{\bm{\beta}}_{a,{\mbox{\tiny NR}}}(\mathcal D_\mathcal J,\widehat{\boldsymbol{\alpha}}_a)$ with

Now we introduce the dynamic treatment Lasso (DTL) estimator of $\theta$:

where $\psi_c(\cdot;\cdot)$ is defined in (ref) and $\hat\eta_{c,{\mbox{\tiny NR}}}:=(\widehat{\boldsymbol{\alpha}}_c,\widehat{\boldsymbol{\beta}}_{c,{\mbox{\tiny NR}}},\widehat{\boldsymbol{\gamma}}_c, \widehat{\boldsymbol{\delta}}_c)$ are the nuisance estimates as in (ref), (ref), (ref), and (ref), respectively; see Algorithm (ref) for details.

Comparisons between the first-time working models $\mu_a^*(\cdot)$ and $\mu_{a,{\mbox{\tiny NR}}}^*(\cdot)$

In the dynamic treatment setting, the relationship between the linear conditional mean function, $\mu_a(\cdot)$, and its corresponding approximations, $\mu_a^*(\cdot)$ and $\mu_{a,{\mbox{\tiny NR}}}^*(\cdot)$, obtained via different identification strategies, is not straightforward. Specifically, a linear $\mu_a(\cdot)$ is only a necessary condition for correctly specified linear working models; it does not guarantee equality. Additional conditions are required to ensure that $\mu_a^*(\cdot)=\mu_a(\cdot)$ and $\mu_{a,{\mbox{\tiny NR}}}^*(\cdot)=\mu_a(\cdot)$. In the following, we will discuss these necessary conditions in detail.

\paragraph*{The working model $\mu_a^*(\cdot)$} The proposed S-DRL approach utilizes the following identification $\mu_a(\mathbf{s}_1)=E[Y^{\mbox{\tiny DR}}\mid\mathbf{S}_1=\mathbf{s}_1,A_1=a_1]$. This representation remains valid for a linear $\mu_a(\cdot)$ as long as either $\rho_a(\cdot)$ is truly logistic or $\nu_a(\cdot)$ is truly linear and not necessarily both. In this case, we also have the following equivalent expressions of $\boldsymbol{\beta}_a^*$:

That is, $\mu_a^*(\mathbf{s}_1)=\mathbf{v}^{\top}\boldsymbol{\beta}_a^*$ satisfies $\mu_a^*(\cdot)=\mu_a(\cdot)$. Hence, $\mu_a^*(\cdot)=\mu_a(\cdot)$ whenever (a) $\mu_a(\cdot)$ is a linear function and (b) either $\rho_a(\cdot)$ is a logistic function or $\nu_a(\cdot)$ is a linear function. It is worth noting that Condition (b) is already a prerequisite for the identification of DTE, as stated in Lemma (ref). Consequently, there is no need to introduce any other conditions beyond those outlined in Condition (a).

\paragraph*{The working model $\mu_{a,{\mbox{\tiny NR}}}^*(\cdot)$} The DTL estimator relies on the nested-regression identification for which Conditions (a) and (b) above are insufficient for $\mu_{a,{\mbox{\tiny NR}}}^*(\cdot)=\mu_a(\cdot)$; see Example (ref) below. Additional conditions are needed. For instance, $\mu_{a,{\mbox{\tiny NR}}}^*(\cdot)=\mu_a(\cdot)$ if we further assume the following:

One sufficient but not necessary condition for (c) is that the second-time conditional mean function $\nu_a(\cdot)$ is also linear; see further justifications in Section A of Supplementary Materials bradic2023supplement.

\paragraph*{The misspecification errors of $\mu_a^*(\cdot)$ and $\mu_{a,{\mbox{\tiny NR}}}^*(\cdot)$} Now we consider the case where $\mu_a(\cdot)$ is possibly non-linear and compare the misspecification (approximation) errors of the linear working models $\mu_a^*(\cdot)$ and $\mu_{a,{\mbox{\tiny NR}}}^*(\cdot)$. As long as Condition (b) above holds, we have

Hence, we have the following conclusions. (1) When $\nu_a(\cdot)$ is linear, $\mu_{a,{\mbox{\tiny NR}}}^*(\cdot)=\mu_a^*(\cdot)$, and both of them are the best linear approximations of the true conditional mean $\mu_a(\cdot)$ among the group $A_1=a_1$, i.e., $\mathrm{Err}_{\mbox{\tiny NR}}=\mathrm{Err}_{\mbox{\tiny DR}}$. (2) When $\nu_a(\cdot)$ is non-linear and $\rho_a(\cdot)$ is logistic, $\mu_{a,{\mbox{\tiny NR}}}^*(\cdot)\neq\mu_a^*(\cdot)$ and $\boldsymbol{\beta}_{a,{\mbox{\tiny NR}}}^*\neq\boldsymbol{\beta}_{a}^*$ in general. If Assumption (ref) holds, the strict inequality above holds in that $\mathrm{Err}_{\mbox{\tiny NR}}>\mathrm{Err}_{\mbox{\tiny DR}}$ as long as $\boldsymbol{\beta}^*_{a,{\mbox{\tiny NR}}}\neq\boldsymbol{\beta}_{a}^*$. That is, the nested “best linear approximation”, $\mu_{a,{\mbox{\tiny NR}}}^*(\cdot)$, is in general sub-optimal. If we further consider the case where $\mu_a(\cdot)$ is linear, then we have $\mathrm{Err}_{\mbox{\tiny DR}}=0$ and possibly $\mathrm{Err}_{\mbox{\tiny NR}}>0$; see also an illustration in Example (ref) below.

The general DR DTE estimator

In this section, we present a general doubly robust (DR) estimator of the DTE. We assume that we have access to estimators $\widehat{\nu}_a(\cdot)$, $\widehat{\mu}_a(\cdot)$, $\widehat{\pi}_a(\cdot)$, and $\widehat{\rho}_a(\cdot)$ of $\nu_a(\cdot)$, $\mu_a(\cdot)$, $\pi_a(\cdot)$, and $\rho_a(\cdot)$, respectively. The functions $\nu_a(\mathbf{s}) $, $\pi_a(\mathbf{s}_1) $, and $\rho_a(\mathbf{s}) $ can be directly estimated using observable variables, while the remaining nuisance function $\mu_a(\cdot)$ can be identified using either the proposed DR representation (ref) or the usual nested representation (ref). We consider flexible estimation strategies for all nuisance functions, including both parametric and non-parametric methods. Using the DR representation of $\theta_a$ given by (ref), we propose a general DR estimator of the DTE through a cross-fitting procedure. For any $K\geq2$, randomly split $\mathcal I=\{1,\dots,N\}$ into $K$ equal-sized parts with $|\mathcal I_{k}|=n=N/K$. For the sake of simplicity, we consider $n$ as an integer. Based on the training samples $\mathcal W_{-k}$, construct $\widehat{\nu}_{c,-k}(\cdot)$, $\widehat{\mu}_{c,-k}(\cdot)$, $\widehat{\pi}_{c,-k}(\cdot)$, and $\widehat{\rho}_{c,-k}(\cdot)$ as estimates of the nuisance functions $\nu_c(\cdot)$, $\mu_c(\cdot)$, $\pi_c(\cdot)$, and $\rho_c(\cdot)$, respectively. For each $c\in\{a,a'\}$, let

The general DR DTE estimator and the corresponding variance estimate are then defined with $\widehat \Delta_{-k}(\cdot) = \widehat{\psi}_{a,-k} (\cdot) -\widehat{\psi}_{a',-k} (\cdot) $ and $\mathcal{K}=\{1,\cdots, K\}$ as

Asymptotic properties

Here we establish consistency and asymptotic normality of the S-DRL, DTL, and the general DR estimator.

Properties of the S-DRL estimator

We use $s_{\bm{\alpha}_{a}}:=\|\boldsymbol{\alpha}_a^*\|_0$, $s_{\bm{\beta}_{a}}:=\|\boldsymbol{\beta}_a^*\|_0$, $s_{\bm{\gamma}_{a}}:=\|\boldsymbol{\gamma}_a^*\|_0$, and $s_{\bm{\delta}_{a}}:=\|\boldsymbol{\delta}_a^*\|_0$ to denote sparsity levels of the nuisance parameters as defined in (ref) , (ref), (ref) and (ref), respectively. The number of covariates, $d_1$ and $d$, are possibly much larger than $N$; for simplicity, we consider $d_1\asymp d_2\asymp d:=d_1+d_2$.

Assumptions (ref) and (ref) are fairly general even among the high-dimensional literature. As $N \to \infty$, we allow $\psi_2$-norm bounds of $\zeta$ and $\varepsilon$ to diverge or to shrink to zero. When all the nuisance models are correctly specified, under the overlap condition in Assumption (ref), $\sigma^2\asymp E[\zeta^2]+E[\varepsilon^2]+E[\xi^2]\geq\max\{E[\zeta^2],E[\varepsilon^2]\},$ where $\xi :=\mu_{a}(\mathbf{S}_{1})-\mu_{a'}(\mathbf{S}_{1})-\theta$ denotes the centered conditional effect at the first exposure. A sufficient condition for Assumption (ref) is $\|\zeta/\sqrt{E[\zeta^2]}\|_{\psi_2} \leq\sigma_\zeta \mbox{ and } \|\varepsilon/\sqrt{E[\varepsilon^2]}\|_{\psi_2}\leq\sigma_\varepsilon,$ i.e., the “normalized” residuals have constant ${\psi_2}$-norms. Note that, we allow $\sigma=\sigma_N$ to be dependent on $N$ while assuming $\sigma_\zeta$ and $\sigma_\varepsilon$ to be constants independent of $N$; $\sigma\to0$ and $\sigma\to\infty$ are both allowed as $N\to\infty$. The following Assumption (ref) is an overlap condition for the working propensity score models, which is additionally required only when model misspecification occurs.

The following theorem characterizes the consistency rate of the S-DRL estimator of $\theta$.

We categorize bounded and unbounded covariate support and add a $\log(d)$ restriction to $s_{\boldsymbol{\delta}_a}$ for the latter. The DR imputed outcome (ref) has an unbounded $\psi_\alpha$-Orlicz norm for any $\alpha>0$. Yet, if $\mathbf{S}_1$ has bounded support, no extra sparsity condition is required as the inverse probability weighting is stable and the DR imputation has a well-behaved tail distribution.

When all the nuisance models are correctly specified, we further establish asymptotic normality results and the corresponding rate DR property of the S-DRL estimator.

Theorem (ref), as per bang2005doubly, indicates that the S-DRL estimator achieves the semiparametric efficiency bound when all nuisance models are correctly specified.

Properties of the DTL estimator

With slight abuse of notation, let $s_{\boldsymbol{\beta}_a}=\|\boldsymbol{\beta}_{a,{\mbox{\tiny NR}}}^*\|_0$. \vskip -10pt

Properties of the general DR estimator

In this section, we provide a new consistency result of the general DR DTE estimator. Here we consider arbitrary working models $\pi_a^*(\cdot)$, $\rho_a^*(\cdot)$, $\mu_a^*(\cdot)$, and $\nu_a^*(\cdot)$, which may not follow the logistic or linear forms as before. For each $c\in\{a,a'\}$, define the corresponding DR score function as

with

The probability measures and corresponding expectations above are with respect to a fresh draw $\mathbf{S}$ (or $\mathbf{S}_1$). Note that Assumption (ref) allows for $\rho_a^*(\cdot)$ to differ from $\rho_a(\cdot)$ while requiring a overlap condition consistent with the existing literature; see, e.g., chernozhukov2018double. The $\max$ condition, satisfied by sub-Gaussian random variables, controls the tails of $\zeta$, $\varepsilon$, and $\xi$. The last two conditions of Assumption (ref) aim to ensure the interpretability of the results by bounding the “normalized” conditional second moments.

The aforementioned theorem yields two distinct conclusions that warrant discussion. The first pertains to the conditions that are necessary for achieving root-$N$ consistency, while the second relates to the issue of consistency under model misspecification. If all the models are correctly specified, $\widehat{\theta}_\mathrm{gen}-\theta=O_p(b_Nc_N+a_Nd_N+\sigma N^{-1/2})$ and root-$N$ consistency happens as long as $b_Nc_N+a_Nd_N=O(N^{-1/2})$ and $\sigma=O(1)$.

If, on the other hand, at least one of the nuisance models is correctly specified at each exposure time, $\widehat{\theta}_\mathrm{gen}$ is a consistent estimator as long as $\sigma=O(1)$. Model misspecification can take an asymmetric form in terms of estimation rates. Specifically, while $q_N$ is symmetric in the rates themselves, the dependence of $b_N$ on $a_N$ and/or $d_N$ can introduce potential asymmetries. For example, when performing $\ell_1$-regularized nested regression, Theorem (ref) indicates that $b_N$ depends additively on $a_N$ as $b_N=b_N^*+a_N$, where $\{b_N^*\}^2= {s_{\boldsymbol{\beta}_a}\log(d)/N}$ is the estimation error of $\mu_a(\cdot)$ when $\nu_a(\cdot)$ is known. As a result, the consistency rate of $\widehat{\theta}_\mathrm{gen}$ includes an additional term $a_Nc_N+a_N\mathbbm1_{\{\pi_a^*\neq\pi_a\}}$, as illustrated by the DTL estimator in (ref).

On the other hand, if we consider a new DR approach based on the DR representation (ref) to estimate $\mu_a(\cdot)$, the corresponding $b_N$ will depend on both $a_N$ and $d_N$. For instance, when all the nuisance models are correctly specified, Theorem (ref) indicates that $\ell_1$-regularized DR estimation leads to a symmetric rate with $b_N=b_N^*+a_Nd_N$, resulting in $\widehat{\theta}_\mathrm{gen}-\theta=O_p(b_N^*c_N+a_Nd_N+1/\sqrt N)$ if $\sigma\asymp1$ and $a_N,b_N^*,c_N,d_N=o(1)$. The approach used to estimate the first-time conditional mean $\mu_a(\cdot)$ determines the persistence of the symmetry.

Supporting theoretical discoveries

This section presents supplementary findings that, while not the primary focus of the research, may nonetheless be informative or valuable.

An adaptive theory for imputed Lasso with high-dimensional covariates

Let $\mathbb{S}:=(Y_i^*,\mathbf{X}_i)_{i=1}^M$ be i.i.d. observations and let $(Y^*,\mathbf{X})$ be an independent copy with $Y^*\in\mathbb{R}$ and $\mathbf{X}\in\mathbb{R}^d$. Suppose that there exists, possibly random, $\widehat{Y}_i\in\mathbb{R}$. Note that for some, and possibly all observations, outcomes $Y^*$ are imputed, i.e., estimated using $\hat Y_i$. The true population slope is defined as $\boldsymbol{\beta}^*=\mathrm{argmin}_{\boldsymbol{\beta}\in\mathbb{R}^d}E [Y^*-\mathbf{X}^\top\boldsymbol{\beta} ]^2.$ Then its estimator is

for $\lambda_M>0$. The following result delineates properties of such imputed-Lasso estimator, $\widehat{\boldsymbol{\beta}}$.

The above result is of independent interests as it provides a general theory for any Lasso estimators based on imputed outcomes. It contributes to the literature in three specific aspects: (a) The “imputation error”, $\widehat{Y}_i-Y_i^*$, can be dependent on and even possibly correlated with covariates $\mathbf{X}_i$; (b) We allow every $\widehat{Y}_i$ to be fitted using the same set of observations $(X_i,Y_i)_{i=1}^M$, i.e., $\widehat{Y}_i$s are also possibly dependent on each other; and (c) The tuning parameter $\lambda_M$ is of the same order as the one chosen for the fully observed data and is independent of the imputation error or any sparsity parameter.

Compared with the existing literature, Theorem (ref) requires weaker sparsity assumptions and provides better rates of estimation. Imputed Lasso of zhu2019proper requires $s=o(M)$, $\log(d)=o(\sqrt M)$, and $\sqrt s\delta_M=o(1)$. That of lewis2021double requires an ultra-sparse setup $s^2\log(d)=o(M)$ and $s\delta_M=o(1)$. In contrast, we only require $s\log(d)=o(M)$ and $\delta_M=o(1)$. Additionally, zhu2019proper choose a tuning parameter $\lambda_M\gg\sqrt s\delta_M$ and provide $\| \widehat{\boldsymbol{\beta}} -\boldsymbol{\beta}^* \|_2 =O_p(\sqrt s\delta_M+\sqrt{s/M})$ upon requiring strict conditions for assuring model selection consistency; see Theorem 2 therein. lewis2021double take $\lambda_M\asymp\sqrt{\log(d)/M}+\delta_M$ and establish $\| \widehat{\boldsymbol{\beta}} -\boldsymbol{\beta}^* \|_2 =O_p(s\sqrt{\log (d)/M}+ s\delta_M)$; see Theorem 13 therein. In contrast, we allow $\lambda_M\asymp\sqrt{\log(d)/M}$. The imputation error $\delta_M$ only appears in our final estimation rate (ref) additively, and its effect does not explode as the sparsity level grows.

Theorem (ref) requires development of new proof techniques: the standard Lasso inequality followed by the cone-set reduction are not valid in this instance. In fact, the error, $\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^* $, no longer belongs to the accustomed cone set, $\mathcal C(S,k):=\{\boldsymbol{\Delta}\in\mathbb{R}^d:\|\boldsymbol{\Delta}_{S^c}\|_1\leq k\|\boldsymbol{\Delta}_S\|_1\}$. Instead, we identify a new set, $\widetilde{\mathcal C}(S):=\{\boldsymbol{\Delta}\in\mathbb{R}^d:\|\boldsymbol{\Delta}_{S^c}\|_1\leq 16 \lambda_M^{-1} \delta_M^2, \|\boldsymbol{\Delta}_{S}\|_1\leq 4\lambda_M^{-1} \delta_M^2\}$, and show that the error vector belongs to the union of the above two sets. This enables us to avoid choosing a tuning parameter dependent on the imputation error, as is done in the above literature. Moreover, our results are adaptive to the imputation error in that when there is no imputation, i.e., $\delta_M=0$, our result reaches the standard consistency rate in the high-dimensional statistics literature, e.g., bickel2009simultaneous,negahban2012unified,wainwright2019high.

Theoretical characteristics of nuisance estimators with imputed outcomes

As a result of constraints on the length of the main file, we have included the theoretical properties of the nuisance estimates $\widehat{\boldsymbol{\alpha}}_a$, $\widehat{\boldsymbol{\gamma}}_a$, and $\widehat{\boldsymbol{\delta}}_a$ as defined by equations (ref), (ref), and (ref) respectively, in the Supplementary Materials {bradic2023supplement where we show $\|\widehat{\bm{\alpha}}_a-\bm{\alpha}_a^*\|_2=O_p(\sigma\sqrt{s_{\bm{\alpha}_{a}}\log(d)/N})$, $\|\widehat{\bm{\gamma}}_a-\bm{\gamma}_a^*\|_2=O_p(\sqrt{s_{\bm{\gamma}_{a}}\log(d_1)/N})$, and $\|\widehat{\bm{\delta}}_a-\bm{\delta}_a^*\|_2=O_p(\sqrt{s_{\bm{\delta}_{a}}\log(d)/N})$. Now we establish the properties of the first-time conditional mean model estimates, where imputation is required. We first consider the DR-imputation-based estimator $\widehat{\boldsymbol{\beta}}_a$ defined as (ref) and the corresponding conditional mean estimate $\widehat{\mu}_a(\mathbf{s}_1)=\mathbf{v}^{\top}\widehat{\boldsymbol{\beta}}_a$.

Theorem (ref) elucidates that the consistency rate of $\widehat{\boldsymbol{\beta}}_a$ is subject to the fidelity of the second-time nuisance models $\rho_a^*(\cdot)$ and $\nu_a^*(\cdot).$ More specifically, when both models are accurately specified, the DR imputation error contributes multiplicatively to the consistency rate in Theorem (ref)(a). In contrast, when only one of $\rho_a^*(\cdot)$ and $\nu_a^*(\cdot)$ is correctly specified, the estimation error of the correctly specified model contributes additively to the rates presented in Theorem (ref)(b) and Theorem (ref)(c). It is noteworthy that these results do not rely on the correctness of the first-time conditional mean model per se.

In the following, we also provide the consistency results of the nested-regression-based estimator, $\widehat{\boldsymbol{\beta}}_{a,{\mbox{\tiny NR}}}$, defined in Equation (ref), and the corresponding conditional mean estimate $\widehat{\mu}_{a,{\mbox{\tiny NR}}}(\mathbf{s}_1)=\mathbf{v}^{\top}\widehat{\boldsymbol{\beta}}_{a,{\mbox{\tiny NR}}}$.

In general, $\mu_a(\cdot)-\mu_{a'}(\cdot)$ can be seen as a conditional average treatment effect (CATE) parameter through the well established nested representation (ref). Outside of dynamic settings, DR approaches for CATE estimation typically rely on DR influence function representation of the conditional means. When those conditional means independently are not smooth enough, kennedy2020towards proposes to instead use DR imputations for the joint estimation of the difference of the conditional means. Here, the nested structure of $\mu_a(\cdot)$, where the true outcome is never observed, prevents direct influence function approaches. Instead, our approach leverages cases when $\mu_a(\cdot)$ has {\it sparser} structure than $\nu_a(\cdot)$.

Advancing multi-stage treatment estimation with DR methods

The objective of this section is to expand upon the methodology of sequential doubly robust estimation by considering its application in multi-stage settings. Consider $T\geq2$ exposure times and suppose that we observe i.i.d. samples $\{W_{T,i}\}_{i=1}^N= (\mathbf{S}_{1i},A_{1i},\dots,\mathbf{S}_{Ti},A_{Ti},Y)_{i=1}^N$. Let $W _T:= (\mathbf{S}_1,A_{1},\dots,\mathbf{S}_T,A_T,Y)$ be an independent copy of $W_{T,i}$. For each $t\leq T$, let $\mathbf{S}_t\in\mathbb{R}^{d_t}$ and $A_t\in\{0,1\}$ denote the covariate vector and the treatment assignment at the $t$-th exposure time, respectively. Let $Y \in \mathbb{R}$ denote the observed outcome variable at the final stage. Denote $\overline{\mathbf{S}}_t:=(\mathbf{S}_1,\dots,\mathbf{S}_t)$ and $\overline{A}_t:=(A_1,\dots,A_t)$ for any $1\leq t\leq T$. Let $Y(\overline{a}_T)$ be the counterfactual outcome corresponding to the treatment path $\overline{a}_T=(a_1,\dots,a_T)\in\{0,1\}^T$. The DTE between any treatment paths $\overline{a}_T,\overline{a}_T'\in\{0,1\}^T$ is now defined as

We define the conditional mean and propensity score functions as

where for the sake of simplicity, we denote with $\bar A_0=\overline{a}_0=\emptyset$ and $\overline{\mathbf{S}}_{T+1}:=(\mathbf{S}_1,\dots,\mathbf{S}_T,Y)$. For each $1\leq t\leq T$, we denote $\mu_{t}^*(\overline{\mathbf{s}}_{t},\overline{a}_{T})$ and $\pi_{t}^*(\overline{\mathbf{s}}_{t},\overline{a}_{t})$ as the working models for the conditional mean and propensity score, respectively. Additionally, with $\overline{\mathbf{S}}_0=\overline{\mathbf{s}}_0=\emptyset$, we set $\mu_{0}(\overline{\mathbf{s}}_{0},\overline{a}_{T}):=E[Y(\overline{a}_T)|\overline{\mathbf{S}}_0=\overline{\mathbf{s}}_0]=\theta_{\overline{a}_T}$ and $\mu_{T+1}^{*}(\overline{\mathbf{s}}_{T+1},\overline{a}_{T}):=s_{T+1}$. Note that, under the Assumption (ref)(b) below, we have $\mu_{T+1}^*(\overline{\mathbf{S}}_{T+1},\overline{a}_{T})=\mu_{T+1}(\overline{\mathbf{S}}_{T+1},\overline{a}_{T})=Y$. To identify $\theta_{\overline{a}_T}=E[Y(\overline{a}_T)]$ for any $\overline{a}_T\in\{0,1\}^T$, we assume a multi-stage version of Assumption (ref) in the following; see also, e.g., murphy2003optimal,robins2000marginal,robins1987addendum.

The following proposition presents a well-known DR representation of $E[Y(\overline{a}_t)]$ under the multi-stage dynamic setting; see, e.g., bang2005doubly, murphy2001marginal.

According to Proposition (ref), a consistent estimate of $\theta_{\overline{a}_T}$ should be achievable as long as we can consistently estimate at least one of the nuisance functions $\mu_{t}(\cdot,\overline{a}_{T})$ or $\pi_{t}(\cdot,\overline{a}_{t})$ for each exposure time $t$. In the present context, the propensity score functions of (ref) are identifiable via observable variables. Additionally, by Assumption (ref), $\mu_T(\overline{\mathbf{s}}_T,\overline{a}_{T})= E[Y|\overline{\mathbf{S}}_T=\overline{\mathbf{s}}_T,\overline{A}_{T}=\overline{a}_{T}]$, thereby facilitating its estimation using the corresponding samples. However, the remaining conditional mean functions for stages $ t\leq T-1$ cannot be identified directly. To address this challenge, we propose DR representations of these intermediate conditional means, as an alternative to the conventional nested representation of (ref).

Theorem (ref) can be regarded as an overarching, comprehensive, umbrella result that subsumes a range of components, particularly encompassing Proposition (ref) as a specific case when $t=0$. Indeed, $\theta_{\overline{a}_T}=E[Y(\overline{a}_T)|\overline{\mathbf{S}}_0=\overline{\mathbf{s}}_0]=\mu_0(\overline{\mathbf{s}}_0,\overline{a}_T)$ is a conditional mean function at “stage zero”. Theorem (ref) indicates that $\mu_{t}(\cdot,\overline{a}_{T})$ can be identified through a DR representation using all the conditional means and propensity scores at later stages. Therefore, $\mu_{t}(\cdot,\overline{a}_{T})$ can be estimated sequentially backward in time based on the DR imputations.

For example, if we use linear working models for the conditional means and logistic models for the propensity scores, and either the true conditional mean $\mu_{r}(\cdot,\overline{a}_{T})$ is linear or the true propensity score $\pi_{r}(\cdot,\overline{a}_{r})$ is logistic at each later stage $r\geq t+1$, we can get a consistent estimate of the $\mu_{t}(\cdot,\overline{a}_{T})$ using DR imputed linear regression. By repeating this process backwards, we conclude that if either the conditional means or propensity scores are correctly parametrized at every stage $t$, we can estimate all nuisance functions consistently, leading to a consistent estimate of $\theta_{\overline{a}_T}$. An alternative approach to our proposed sequential doubly robust method is the nested estimator murphy2001marginal. This approach represents all conditional means using the following equation:

However, in order to ensure the consistency of the nested estimator for $\mu_{t}(\cdot,\overline{a}_{T})$, it is essential that all subsequent conditional mean functions exhibit true linearity. Interestingly, even the multiply robust approach presented by babino2019multiple falls short in achieving the same level of robustness as the S-DRL method. Further insights can be found in the comments following Theorem (ref) and Table (ref). Our method demonstrates a growing advantage as the number of exposure times increases. For instance, in the case of $T$ exposure times, consider all the cases including correctly or incorrectly parametrized $\pi_t^*(\cdot,\overline{a}_t)$ and $\mu_t^*(\cdot,\overline{a}_T)$ for each $t$, our method enables $3^T$ out of $4^T$ possible cases. In contrast, the nested-regression-based and multiply robust approaches only allow for $(T+2)2^{T-1}$ cases. This conclusion is independent of the particular parametrization used -- nonparametric, smooth models are permissible -- and extends to the difference of means $\theta=\theta_{\overline{a}_T}-\theta_{\overline{a}'_T}$ as well.

Numerical Experiments

Simulation studies

We illustrate the finite sample properties of the introduced estimators in several simulated experiments; auxiliary settings are relegated to Section C of the Supplementary Material bradic2023supplement. We consider $a=(1,1)$ and $a'=(0,0)$, and use $\mathbf{1}_{(q)}:=(1,\dots,1)^\top\in\mathbb{R}^{q}$ as well as $\mathbf0_{(q)}:=(0,\ldots,0)\in\mathbb{R}^q$. Below we use $\zeta_i\sim^\mathrm{iid}\mathrm{Uniform}(-1,1)$ and $\{\boldsymbol{\delta}_{i}\}_j\sim^\mathrm{iid}\mathrm{Uniform}(-1,1)$. We decompose $\boldsymbol{\alpha}_a$ into two components as $\boldsymbol{\alpha}_a =(\boldsymbol{\alpha}_{a,1}^\top,\boldsymbol{\alpha}_{a,2}^\top)^\top$ and consider $\boldsymbol{\alpha}_{a'}=-\boldsymbol{\alpha}_a$ and $\boldsymbol{\eta}_{a'}=-\boldsymbol{\eta}_a$ unless specified differently. For each $c\in\{a,a'\}$, $\rho_c(\mathbf{S}_{i})=g(\mathbf{U}_i^\top\boldsymbol{\eta}_c)$ and $A_{2i}|(\mathbf{S}_{i},A_{1i}=c_1)\sim\mathrm{Bernoulli}(\rho_{c}(\mathbf{S}_{i}))$.

\paragraph*{M1: Correctly parametrized models} Consider $\mathbf{S}_{1i}\sim^\mathrm{iid} N_{d_1}(\mathbf{0},\mathbf{I}_{d_1})$ and $A_{1i}|\mathbf{S}_{1i}\sim\mathrm{Bernoulli}$ \ $(\pi_{a} (\mathbf{S}_{1i}))$, with $\pi_{a} (\mathbf{S}_{1i})=g(\mathbf{V}_i^\top\boldsymbol{\gamma}_{a})$. Let $\delta_{1i}\sim^\mathrm{iid} N(0,1)$, $\boldsymbol{\delta}_{1i}\sim^\mathrm{iid} N_{d_1}(0,\mathbf{I}_{d_1})$, and $$\mathbf{S}_{2i}=\mathbf{S}_{1i}+A_{1i}(1+\delta_{1i})\mathbf{1}_{(d_1)}+\boldsymbol{\delta}_{1i}.$$ The outcomes are $Y_i(c)=\mathbf{U}_i^\top\boldsymbol{\alpha}_c+ N(0,1)$ with parameters $\boldsymbol{\alpha}_{a} =(-1,-1,1,-1,\mathbf0_{(d_1-3)},$ $-1,-1,1,\mathbf0_{(d_2-3)})^\top,$ $\boldsymbol{\alpha}_{a'} =(1,1,1,-1,\mathbf0_{(d_1-3)},1,1,1,\mathbf0_{(d_2-3)})^\top,$ $\boldsymbol{\gamma}_{a} =(0,1,1,1,\mathbf0_{(d_1-3)})^\top,$ $\boldsymbol{\eta}_a =(0,1,1,\mathbf0_{(d_1-2)},1,-1,\mathbf0_{(d_2-2)})^\top,$ and $\boldsymbol{\eta}_{a'} =(0,0.5,0,-0.5,\mathbf0_{(d_1-3)},0.5,0,0.5,\mathbf0_{(d_2-3)})^\top$.

\paragraph*{M2: Weakly sparse $\nu_c(\cdot)$ and dense $\pi_c(\cdot)$} Let $D_i\sim^\mathrm{iid} \mathrm{Bernoulli}(0.5)$ and $$\{\mathbf{S}_{1i}\}_{j}\sim D_i\cdot\mathrm{Uniform}(-1,-0.5)+(1-D_i)\cdot\mathrm{Uniform}(0.5,1).$$ Define $\{\overline{W}\}_j=0.5\cdot 0.9^j$ and let $\mathbf{S}_{2i}=\mathbf{S}_{1i}^\top\overline{W}\mathbf{1}_{(d_2)}+\boldsymbol{\delta}_{i}$. Let $Y_i(c)=\mathbf{U}_i^\top\boldsymbol{\alpha}_c+\zeta_i$. The parameter $\boldsymbol{\alpha}_{a,1} =(-1,\mathbf0_{(d_1)})^\top$ and $\{\boldsymbol{\alpha}_{a,2}\}_j=-0.3\cdot 0.99^{j-1}$, $\boldsymbol{\eta}_a =(3,0.1,\mathbf0_{(d_1+d_2-1)})^\top$.

\paragraph*{M3: Non-linear $\nu_c(\cdot)$ and non-logistic $\pi_{c}(\cdot)$} Consider $\{\mathbf{S}_{1i}\}_{j}\sim^\mathrm{iid}\mathrm{Uniform}(-1,1)$. Let $\pi_{a}(\mathbf{S}_{1i})=\bar g(\mathbf{V}_i^\top\boldsymbol{\gamma}_{a})$, where $$\bar g (u)=(|u|/(|u|+1))\mathbbm{1}_{\{u>0\}}+(1/(|u|+1))\mathbbm{1}_{\{u<0\}}.$$ Define $\{\widetilde{W}(a)\}_j=0.7 \cdot 0.8^j$ and $\{\widetilde{W}(a')\}_j=0.5\cdot 0.9^j$. Let $\mathbf{W}_{2i}=\{\widetilde W(A_{1i})\}^\top\mathbf{S}_{1i}\mathbf{1}_{(d_2)}+\boldsymbol{\delta}_{i}$ and $\{\mathbf{S}_{2i}\}_j=\sqrt{|\{\mathbf{W}_{2i}\}_j|}$. The parameter $\boldsymbol{\alpha}_a$ has the same $\boldsymbol{\alpha}_{a,1}$ as M2 and $\{\boldsymbol{\alpha}_{a,2}\}_j=-0.3\cdot 0.9^{j-1}$. The parameter $\boldsymbol{\gamma}_{a} =5\cdot\mathbf1_{(10)}^\top,$ and $\boldsymbol{\eta}_a =(2,0.1,\mathbf0_{(d_1-1)},0.1,\mathbf0_{(d_2-1)})^\top$. Let $$Y_i(c)=\mathbf{V}_i^\top\boldsymbol{\alpha}_{c,1}+\sum_{j=1}^{d_2}\{\alpha_{c,2}\}_j\mathrm{sgn}(\{\mathbf{W}_{2i}\}_j)\{\mathbf{S}_{2i}\}_j^2+\zeta_i.$$

For M1, $d_1=d_2=100$ and $N\in\{1000,4000\}$; for M2, $d_2=500$, $(N,d_1)$ are chosen from $(2000,20)$, $(4000,20)$, and $(4000,50)$; for M3, $d_1=20$, $(N,d_2)$ are chosen from $(500,500)$, $(1000,500)$, $(2000,500)$, $(4000,500)$, $(1000,1000)$, and $(2000,1000)$. We replicate settings 500 times. We report S-DRL as well as a version S-DRL', which has $\widehat{\bm{\beta}}_{c}$ constructed with $\widehat{\boldsymbol{\delta}}_c$ and $\widehat{\boldsymbol{\alpha}}_c$ build on the whole sub-sample $\mathcal W_{-k}=\mathcal W_{-k,1}\cup\mathcal W_{-k,2}$.} We also present (a) DTL, Algorithm (ref), (b) IPW with $\ell_1$-regularized logistic PS, (c) an empirical difference estimator (empdiff), $\widehat{\theta}_{\mbox{\tiny empdiff}}:=\sum_{i=1}^NA_{1i}A_{2i}Y_i/\sum_{i=1}^NA_{1i}A_{2i}-\sum_{i=1}^N(1-A_{1i})(1-A_{2i})Y_i/\sum_{i=1}^N(1-A_{1i})(1-A_{2i})$, and (d) an oracle DR estimator constructed with the true nuisances. All methods use $10$-fold cross validation for selection of tuning parameters.

Tables (ref) and (ref) show the estimation and inference results for the DTE estimators, while Table (ref) focuses on estimation performances, as valid inference is unlikely with misspecified models. Our summarized findings, shown in Tables (ref)-(ref), reveal that the naive empirical difference estimator has large biases due to confounding between outcome and treatment assignments. The IPW method also performs poorly, with large biases and RMSEs, and confidence interval coverages far from the desired $95\%$. The DTL, S-DRL, and S-DRL' estimators behave similarly in Table (ref) (under M1), with correctly specified nuisance models and relatively low sparsity levels. The S-DRL method's additional sample splitting in Algorithm (ref) (Steps 5-7) leads to larger estimation errors in the first-time conditional mean estimates than those in the DTL and S-DRL' methods. Consequently, when $N=1000$, the S-DRL estimator's RMSE is slightly larger than that of the DTL and S-DRL' estimators, but they have similar RMSEs when $N=4000$. In terms of inference behaviors, the corresponding confidence interval coverages are below the desired $95\%$ when $N=1000$. However, increasing the total sample size to $N=4000$ brings the coverages closer to $95\%$. Estimating $\nu_c(\cdot)$ under M2 is more challenging than estimating $\rho_c(\cdot)$. As a result, the DR estimates of $\mu_c(\cdot)$ in the S-DRL and S-DRL' methods have significantly smaller estimation errors compared to the nested regression used in the DTL method, as shown in Table (ref), leading to smaller RMSEs and coverages closer to $95\%$. Moving on to M3, both $\nu_c(\cdot)$ and $\pi_c(\cdot)$ are misspecified. Table (ref) shows that the estimation errors of $\mu_c(\cdot)$ with the S-DRL and S-DRL' are substantially smaller than those of the DTL. Consequently, we see an improvement of the RMSEs in the S-DRL and S-DRL' estimators.

Application to National Job Corps Study (NJCS)

Job Corps (JC) is the largest and most comprehensive federal job training program in the US for disadvantaged youth between 16 and 24 years old. Each year, about 50,000 participants receive vocational training and academic education at JC centers to improve their job prospects. On average, a JC student spends 8 months at a local center, completing around 1,100 hours of instruction, which is roughly equivalent to one year of high school. For a more detailed description, refer to schochet2008does and schochet2001national.

Numerous studies have investigated the effects of Job Corps on wages. lee2009training highlighted sample selection issues in their analysis. zhang2008evaluating separated the causal effects of JC enrollment on wages from those on employment. flores2012estimating found that longer exposure to JC training is associated with higher future earnings. chen2015bounds separated the effects of sample selection from noncompliance, while huber2020direct distinguished between the causal direct and indirect effects in the presence of mediators. In addition to studying the effects of Job Corps in single-time treatment settings, researchers have also explored the dynamic treatment setting offered by Job Corps. bodory2022evaluating investigated the effects of JC's educational and training programs and found positive impacts on fourth-year employment compared to no program participation. Meanwhile, singh2021kernel analyzed the total, direct, and indirect dynamic dose response of job training on employment. Their study concluded that a few class hours in the first and second years significantly increase employment in the fourth year. In this section, we will evaluate the effects of sequential job training programs on wages using the S-DRL and DTL methods, as defined in Algorithms (ref) and (ref).

We analyze a dataset of 11,313 individuals, with 6,828 assigned to the Job Corps and 4,485 not. They are interviewed 1, 2, and 4 years post-randomization. For each year $t\in{1,2}$, $Z_t\in\{0,1,2,3\}$ represents the treatment assignment in the $t$-th year. We assign $Z_t=0$ for non-enrollment, $Z_t=1$ for enrollment without program participation, $Z_t=2$ for high-school-level education, and $Z_t=3$ for vocational training. The baseline covariate vector, $\mathbf{S}_1$, has 909 characteristics, while $\mathbf{S}_2$ includes 1,427 characteristics. In total, there are 2,336 covariates. The outcome is the log-transformed wage $\widetilde{Y}=\log(\mathrm{wage}+1)\in\mathbb{R}$. We exclude 2,610 individuals with missing treatment stages that are missing completely at random Schochet2003national and an additional 133 with missing covariates or outcomes, resulting in a final sample of 8,570 individuals.

Table (ref) shows estimated DTEs between treatment paths $(3,3)$ vs. $(1,1)$, $(3,3)$ vs. $(2,2)$, and $(2,2)$ vs. $(1,1)$. Both the S-DRL and DTL methods suggest that vocational training has a positive impact on achieving higher wages by showing non-zero effects between the first two paths. On the other hand, the estimates between paths $(2,2)$ and $(1,1)$ are negative, and their corresponding confidence intervals contain zero, making it impossible to determine if academic education is beneficial or detrimental. However, our analysis does suggest that individuals seeking higher-paying jobs would benefit more from vocational training compared to academic education, which only provides high school-level education without any significant vocational training. The S-DRL estimates have a slightly greater distance from zero compared to DTL's, with similar standard errors, leading to slightly smaller p-values. Figure (ref) examines the overlap of estimated propensity scores, displaying mirror histograms of estimated propensity scores within the treatment and control groups. The substantial overlap seen in the mirror histograms indicates that the inverse propensity score weights are relatively stable. Figure (ref)(d) displays bimodal patterns in the histograms due to a binary confounder variable that significantly influences the propensity score estimate of $P[Z_{2}=2\mid\mathbf{S},Z_{1}=2]$, with a close association between a participant's decision to enroll in second-year education and their attendance in the class during the final weeks of the first year.

Discussion

This paper aims to enhance the understanding of estimating causal parameters in multi-stage settings. While prior DR literature has recognized the importance of the stage-zero DR representation for the expected potential outcome, it has overlooked the fact that all the intermediate conditional mean functions can also be identified in a DR manner. This approach leads to better theoretical guarantees and greater flexibility in modeling dynamic dependencies, which can be complex and involve multiple time exposures. Furthermore, our findings have significant practical implications beyond parametric models, especially in situations where doctors or policymakers cannot rely on randomized treatments or simplistic treatment rules. With the ability to model dynamic treatment effects using robust principles, new avenues of discovery are emerging, including optimal treatment rules and determining the best treatment times. Our approach also enables the exploration of further important issues, such as mitigating network spillover effects through robustness perspectives and enriching balancing methods with better robustness properties. The significance of our work lies in the fact that it allows researchers to estimate treatment effects in complex settings more accurately and provides a valuable tool for policymakers seeking to make informed decisions based on robust causal inference methods.

Acknowledgement

This work was supported in part by NSF awards CNS-1730158, ACI-1540112, ACI-1541349, OAC-1826967, the University of California Office of the President, and the University of California San Diego's California Institute for Telecommunications and Information Technology/Qualcomm Institute. Jelena Bradic's work has been supported by the NSF grand DMS-1712481. The majority of this work was done while Yuqian Zhang was with the Department of Mathematics, University of California San Diego.