Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
94,751 characters · 14 sections · 52 citation commands
Dynamic treatment effects: high-dimensional doubly robust inference under model misspecification
Statistical inference and estimation of causal relationships have a long-standing tradition. In various applications, data is collected dynamically over time, and individuals undergo treatments at multiple stages. Examples include mobile health datasets, electronic health records, and a broad range of studies from biomedicine and public health to political science. Conducting randomized controlled trials, especially those with multiple treatment stages, is often time-consuming, and results are not immediately available. Additionally, financial and ethical constraints frequently result in small sample sizes with unrepresentative populations. In contrast, observational studies provide a more accessible alternative, generating large-scale dynamic datasets with rich information. While observational studies have advantages in terms of economic efficiency, sample representativeness, and timely results, statistical analysis based on them is more challenging. Over time, confounding variables at each time point simultaneously affect future treatments and final outcomes. As a result, methods that are effective in randomized controlled trials, such as the two-sample t-test, generally suffer from non-negligible bias in observational studies.
Bias from dynamic confounders is a major challenge in causal inference. In large-scale dynamic studies, confounders often outnumber treatment-specific samples due to exponentially decreasing sizes across stages, leading to high-dimensional settings. Additionally, model misspecification complicates inference, particularly when earlier counterfactual models depend on later ones babino2019multiple. We demonstrate that this challenge can be overcome by introducing new estimates of the underlying causal effects.
Consider a dynamic setting with binary treatments at two exposure times, $A_1$ and $A_2$, although our results extend to any finite number of exposures. We observe independent and identically distributed samples $\mathcal{S} = \{\mathbf{W}_i\}_{i=1}^N =\{ (Y_i, A_{1i}, A_{2i}, \mathbf{S}_{1i}, \mathbf{S}_{2i})\}_{i=1}^N$. Here, $Y \in \mathbb{R}$ denotes the observed outcome at the final stage. We assume the existence of potential outcome variables $Y(a_1, a_2)$ for each $(a_1, a_2) \in \{0, 1\}^2$, representing the outcome an individual would have experienced if exposed to a treatment path $(a_1, a_2)$. Before each exposure, we also collect confounders (or covariates), denoted by $\mathbf{S}_1 \in \mathbb{R}^{d_1}$ and $\mathbf{S}_2 \in \mathbb{R}^{d_2}$, respectively. Covariate history up to the second exposure is denoted by $\bar{\mathbf{S}}_2 := (\mathbf{S}_1^\top, \mathbf{S}_2^\top)^\top \in \mathbb{R}^d$, where the dimensions $d_1$ and $d := d_1 + d_2$ are potentially much larger than $N$. We consider observational studies that allow all variables evaluated at previous stages to potentially influence later ones, without relying on any Markov assumptions, as illustrated in Figure (ref).
In this work, we concentrate on estimating a causal effect known as the dynamic treatment effect (DTE), defined as $\mathrm{DTE}:=\mathbb{E}\{Y(a_1,a_2)-Y(a_1’,a_2’)\}$. We focus on counterfactual mean \(\theta_{1,1} := \mathbb{E}\{Y(1,1)\}\) because the same method extends to any \(\mathbb{E}\{Y(a_1,a_2)\}\) and thus to the DTE. Because \(Y_i(1,1)\) is observed only when \(A_{1i} = A_{2i} = 1\), we cannot simply average across all potential outcomes \(Y_i(1,1)\) since many remain unobserved. Additionally, due to the presence of confounding, we generally have $$\theta_{1,1} \neq \mathbb{E}\{Y(1,1) \mid A_1 = A_2 = 1\}.$$ Marginal Structural Mean (MSM) models are widely used in causal inference to assess the impact of time-dependent treatments on outcomes, allowing for time-dependent covariates affected by previous treatments robins2000marginal. MSMs are determined by a score function that identifies $\theta_{1,1}$ and a number of nuisance parameters. Our goal is to ensure correct optimal inference — even if nuisance models are misspecified and cannot be estimated at the usual \(\sqrt{N}\)-rate. We begin with the minimal set of assumptions outlined in Assumption 1.
In the standard conditions above (see, e.g., robins1987addendum,robins2000marginal,murphy2003optimal), we include two propensity score (PS) models — \(\pi(\mathbf{s}_1) := \mathbb{P}(A_1 = 1 \mid \mathbf{S}_1 = \mathbf{s}_1)\) and \(\rho(\bar{\mathbf{s}}_2) := \mathbb{P}(A_2 = 1 \mid \bar{\mathbf{S}}_2 = \bar{\mathbf{s}}_2, A_1 = 1)\) — and two conditional, counterfactual, outcome regression (OR) models — \(\mu(\mathbf{s}_1) := \mathbb{E}\{Y(1,1) \mid \mathbf{S}_1 = \mathbf{s}_1\}\) and \(\nu(\bar{\mathbf{s}}_2) := \mathbb{E}\{Y(1,1) \mid \bar{\mathbf{S}}_2 = \bar{\mathbf{s}}_2, A_1 = 1\}\). These four models are the true, yet unknown population processes that will factor into the identification of $\theta_{1,1}$. We introduce “working” models \(\pi^*\), \(\rho^*\), \(\mu^*\), and \(\nu^*\), which need not coincide with the true processes but are used to guide estimation the necessary nuisances. Given this framework, \(\theta_{1,1}\) can be identified in multiple ways using only observable variables. We focus on the doubly robust identification approach, as outlined in murphy2001marginal, bang2005doubly, yu2006double, where
as long as Assumption (ref) holds.
However, modeling the counterfactual mean given covariates and treatment history up to a certain time point inherently imposes constraints on the counterfactual mean given earlier histories. For instance, in the representation $ \mu(\mathbf{s}_1) \;=\; \mathbb{E}\{\nu(\bar{\mathbf{S}}_2) \mid \mathbf{S}_1 = \mathbf{s}_1, A_1 = 1\} $ from murphy2001marginal, the correctness of \(\mu\) hinges on the correctness of \(\nu\). We show that the above double-robust representation is not sufficient to guarantee model double robustness. A method is deemed doubly robust if it yields consistent estimates under Assumption (ref). In contrast, a method is considered model doubly robust (or said to provide doubly robust inference) if its inference remains valid under the same Assumption (ref), as long as the (possibly misspecified) working models can be estimated with $o(N^{-1/4})$ rates smucler2019unifying. New estimates of the nuisances are required to accommodate non-\(\sqrt{N}\) convergence rate while still guaranteeing \(\sqrt{N}\)-rate inference for $\theta_{1,1}$ under the same Assumption (ref).
We illustrate the existing results under the following four settings all encompassed within Assumption (ref): \vskip -10pt \begingroup {0pt} {0pt} {0pt} {0pt}
\endgroup
In low dimensions, inverse probability weighting (IPW) robins1986new,hernan2001marginal,robins2004optimal provide valid inference allowing only (ref). On the other hand, covariate balancing methods kallus2018optimal,yiu2018covariate only allow for (ref). Doubly robust methods bang2005doubly, yu2006double,orellana2010dynamic allow (ref) or (ref), but do not accommodate (ref) and (ref). The multiple robust estimator proposed by babino2019multiple allows for (ref), (ref), or (ref), but does not accommodate (ref). Using the following double-robust imputation step,
luedtke2017sequential propose a consistent estimator and rotnitzky2017multiply develop an asymptotically normal estimator of the DTE under conditions (ref)--(ref), requiring that the nuisance estimates belong to a Donsker class van2000asymptotic. However, Donsker conditions place bounded complexity on the functional class, ensuring \(\sqrt{N}\)-rate convergence of nuisance estimates. These constraints are unsuitable for many nonparametric or high-dimensional models. rotnitzky2017multiply,bodory2022evaluating,diaz2023nonparametric adopt cross-fitting techniques and double machine learning chernozhukov2017double to relax the need for Donsker conditions, thereby allowing more flexible, nonparametric methods. However, valid inference still requires all nuisance estimates to converge sufficiently quickly to their true underlying models, effectively ruling out model misspecification.
Similar picture persists with high-dimensional working models -- in order to guarantee \(\sqrt{N}\)-inference, all working nuisance models must be correctly specified. This is achieved in a sequence of papers across different models. For structural nested mean models Robins1997causal, lewis2021double achieve inferential guarantees only when the “blip functions”—differences in outcome regression across treatment paths—are low-dimensional and correctly specified. viviano2021dynamic require both OR models to be correct, covering only (ref). Dynamic Treatment Lasso (DTL) bodory2022evaluating achieves consistency under (ref), (ref), and (ref), and Sequential Doubly Robust Lasso (S-DRL) bradic2024high extends this to (ref) via (ref). However, both DTL and S-DRL require all nuisance models to be correctly specified for valid inference, thus excluding misspecification in (ref)--(ref).
All the doubly robust methods discussed above share a common limitation: while estimators based on the doubly robust score (ref) reduce bias to product (quadratic) terms of the nuisance estimation errors, this reduction critically depends on correctly specifying all nuisance models. Once any model is misspecified, the bias reduction fails. To address this, we propose new moment conditions that directly control bias under Assumption (ref). To be concrete, we employ linear OR and logistic PS working models: \[ \nu^*(\bar\mathbf{s}_2) = \bar\mathbf{s}_2^\top\boldsymbol{\alpha}^*, \quad \mu^*(\mathbf{s}_1) = \mathbf{s}_1^\top\boldsymbol{\beta}^*, \quad \pi^*(\mathbf{s}_1)=g(\mathbf{s}_1^\top\boldsymbol{\gamma}^*), \quad \rho^*(\bar\mathbf{s}_2)=g(\bar\mathbf{s}_2^\top\boldsymbol{\delta}^*), \] where \(g(u)=\exp(u)/\{\exp(u)+1\}\). These represent the “best” linear or logistic models approximating the true underlying processes buja2019models. With \(\boldsymbol{\eta}^* := (\boldsymbol{\alpha}^{*\top}, \boldsymbol{\beta}^{*\top}, \boldsymbol{\gamma}^{*\top}, \boldsymbol{\delta}^{*\top})^\top\), \(\widehat{\theta}_{1,1}=N^{-1}\sum_{i=1}^N\psi(\mathbf{W}_i;\widehat{\boldsymbol{\eta}})\) for $\widehat{\boldsymbol{\eta}}=(\widehat{\boldsymbol{\alpha}}^\top,\widehat{\boldsymbol{\beta}}^\top,\widehat{\boldsymbol{\gamma}}^\top,\widehat{\boldsymbol{\delta}}^\top)^\top$. Then,
In the above, \(\Delta_1 := N^{-1} \sum_{i=1}^N \psi(\mathbf{W}_i;\boldsymbol{\eta}^*) - \theta_{1,1}\) is \(O_p(N^{-1/2})\) and asymptotically normal under Assumption (ref) and standard moment conditions. \(\Delta_3\) depends quadratically on \(\widehat{\boldsymbol{\eta}} - \boldsymbol{\eta}^*\) and is typically negligible. There are two main strategies to control the main bias term $ \Delta_2 :=N^{-1}\sum_{i=1}^N \nabla_{\boldsymbol{\eta}} \psi(\mathbf{W}_i;\boldsymbol{\eta}^*)^\top (\widehat{\boldsymbol{\eta}} - \boldsymbol{\eta}^*)$. One is to constrain the nuisance models to be either Donsker or low-dimensional, ensuring \(\sqrt{N}\)-rate convergence of $\hat \boldsymbol{\eta} -\boldsymbol{\eta}^*$ -- even if those models are misspecified -- and yielding \(\Delta_2 = o_p(N^{-1/2})\). The other is to apply cross-fitting, which induces independence between the summands in \(N^{-1}\sum_{i=1}^N \nabla_{\boldsymbol{\eta}} \psi(\mathbf{W}_i;\boldsymbol{\eta}^*)\) and thereby allows a weak law of large numbers argument to deliver \(\Delta_2 = o_p(N^{-1/2})\). However, cross-fitting requires \(\mathbb{E}\{\nabla_{\boldsymbol{\eta}}\psi(\mathbf{W};\boldsymbol{\eta}^*)\}=\mathbf{0}\), a “Neyman orthogonality” property chernozhukov2017double that holds only when the nuisance models are correctly specified—thus excluding the misspecification scenario. Whenever misspecification occurs, \(\mathbb{E}\{\nabla_{\boldsymbol{\eta}}\psi(\mathbf{W};\boldsymbol{\eta}^*)\}\neq\mathbf{0}\) and \[ \Delta_2 = \mathbb{E}\{\nabla_{\boldsymbol{\eta}}\psi(\mathbf{W};\boldsymbol{\eta}^*)\}^\top(\widehat{\boldsymbol{\eta}} - \boldsymbol{\eta}^*) \;+\; o_p(N^{-1/2}), \] meaning \(\Delta_2\) depends linearly on the estimation error \(\widehat{\boldsymbol{\eta}} - \boldsymbol{\eta}^*\) even if cross-fitting is used and \(\widehat{\boldsymbol{\eta}}\) is estimated on a separate dataset. Instead, we design new nuisance estimates, named moment-targeted estimates that directly guarantee the following moment condition
This reduction remains effective even when $\widehat{\boldsymbol{\eta}} - \boldsymbol{\eta}^*$ does not converge at the $\sqrt{N}$ rate. We leverage two key components: (i) the standard doubly robust score (ref), and (ii) new moment-targeted nuisance estimates and new loss functions. Both components are needed. The former addresses bias when all models are correctly specified. The latter specifically targets the bias of model misspecification. As shown in Table (ref), this approach transforms linear bias into quadratic, achieving faster convergence rates and robust inference even under model misspecification; see Table (ref). Our approach applies to all cases (ref)-(ref).
In this section, we introduce he sequential model doubly robust (SMDR) estimator. For any $\boldsymbol{\omega} \in \{\boldsymbol{\gamma}, \boldsymbol{\delta}, \boldsymbol{\alpha}, \boldsymbol{\beta}\}$, let $\Delta_{2,\boldsymbol{\omega}}:=\mathbb{E}\{\boldsymbol{\nabla}_{\boldsymbol{\omega}}\psi(\mathbf{W}; \boldsymbol{\eta}^*)\}^\top (\widehat{\boldsymbol{\omega}} - \boldsymbol{\omega}^*)$ be the bias resulting from the nuisance estimation of $\boldsymbol{\omega}^*$. Then, we also have \(\Delta_2 \approx \sum_{\boldsymbol{\omega} \in \{\boldsymbol{\gamma}, \boldsymbol{\delta}, \boldsymbol{\alpha}, \boldsymbol{\beta}\}} \Delta_{2,\boldsymbol{\omega}}\) with
The above equations are listed in an order such that the right-hand side of each equation involves progressively more nuisance parameters. For instance, (ref) only involves $\boldsymbol{\gamma}^*$, while (ref) involves both $\boldsymbol{\gamma}^*$ and $\boldsymbol{\delta}^*$.
We first propose nuisance estimates \(\widehat{\boldsymbol{\gamma}}\), \(\widehat{\boldsymbol{\delta}}\), \(\widehat{\boldsymbol{\alpha}}\), and \(\widehat{\boldsymbol{\beta}}\) in a sequential manner. Define the index sets as \(\mathcal{I}_{\boldsymbol{\gamma}}, \mathcal{I}_{\boldsymbol{\delta}}, \mathcal{I}_{\boldsymbol{\alpha}}, \mathcal{I}_{\boldsymbol{\beta}} \subseteq \{1, \dots, N\}\), see Algorithm (ref). We define the estimate for the first PS, \(\pi^*(\mathbf{s}_1)\), with \(\lambda_{\boldsymbol{\gamma}} > 0\), as
The loss function (ref) is designed to achieve covariate balancing as \[ \mathbb{E}\{w_1 \mathbf{S}_1\} = \mathbb{E}(\mathbf{S}_1) , \qquad w_1:= A_1g^{-1}(\mathbf{S}_1^\top\boldsymbol{\gamma}^*). \] Strong covariate balancing has been used for estimation of single-stage average treatment effect ning2020robust. Moreover, even when the PS model is misspecified, it still strictly controls the first term in the above decomposition by ensuring \(\Delta_{2,\beta}=\mathbb{E}\{\boldsymbol{\nabla}_{\boldsymbol{\beta}}\psi(\mathbf{W}; \boldsymbol{\eta}^*)\} = \mathbf{0}\) for \(\boldsymbol{\gamma}^* := \arg\min_{\boldsymbol{\gamma} \in \mathbb{R}^{d_1}} \mathbb{E}\{\ell_1(\mathbf{W}; \boldsymbol{\gamma})\}\). This in turn, effectively reduces the bias induced by the estimation error of \(\widehat{\boldsymbol{\beta}}\), the nuisance estimate for the first OR model. On the other hand, if the PS model is logistic, i.e., \(\pi(\mathbf{S}_1) = g(\mathbf{S}_1^\top\boldsymbol{\gamma}^0)\) for some \(\boldsymbol{\gamma}^0 \in \mathbb{R}^{d_1}\), then, \(\boldsymbol{\gamma}^* = \boldsymbol{\gamma}^0\); see Section (ref) of the Supplementary Material.
Next, we construct the estimate for the second PS model, $\rho^*(\bar\mathbf{s}_2)$, with \(\lambda_{\boldsymbol{\delta}}>0\) and
This loss function achieves a new kind of covariate balancing where \[ \mathbb{E}\{w_2 w_1\bar\mathbf{S}_2\} = \mathbb{E}(w_1\bar\mathbf{S}_2), \qquad w_2= A_2g^{-1}(\bar\mathbf{S}_2^\top\boldsymbol{\delta}^*). \] One can interpret the above as conditional covariate balancing: given the information from the first time period, we achieve classical balancing at the second time exposure. The weights \(w_1\) align with those of \(\pi^*\) but are essential for correctly targeting the bias term \(\Delta_2\). Specifically, this conditional covariate balancing ensures that $ \Delta_{2,\alpha} =\mathbb{E}\left\{\nabla_{\boldsymbol{\alpha}} \psi(\mathbf{W}; \boldsymbol{\eta}^*)\right\} = \nabla_{\boldsymbol{\delta}} \mathbb{E}\left\{\ell_2(\mathbf{W}; \boldsymbol{\gamma}^*, \boldsymbol{\delta}^*)\right\} = \mathbf{0}, $ for $ \boldsymbol{\delta}^* := \arg\min_{\boldsymbol{\delta} \in \mathbb{R}^{d}} \mathbb{E}\left\{\ell_2(\mathbf{W}; \boldsymbol{\gamma}^*, \boldsymbol{\delta})\right\}, $ even when the models are misspecified. This in turn, effectively reduces the bias induced by the estimation error of \(\widehat{\boldsymbol{\alpha}}\), the nuisance estimate for the second OR model.
For the remaining two OR models, we define moment-targeted nuisance estimators with properly chosen \(\lambda_{\boldsymbol{\alpha}},\lambda_{\boldsymbol{\beta}}>0\) as
The corresponding loss functions are defined as:
The above (ref)-(ref) mitigate the estimation bias by ensuring \[ \mathbb{E}\left\{\boldsymbol{\nabla}_{\boldsymbol{\gamma}}\psi(\mathbf{W}; \boldsymbol{\eta}^*)\right\} = \boldsymbol{\nabla}_{\boldsymbol{\beta}}\mathbb{E}\{\ell_4(\mathbf{W}; \boldsymbol{\gamma}^*, \boldsymbol{\delta}^*, \boldsymbol{\alpha}^*, \boldsymbol{\beta}^*)\}/2 = \mathbf{0}\] \[ \mathbb{E}\left\{\boldsymbol{\nabla}_{\boldsymbol{\delta}}\psi(\mathbf{W}; \boldsymbol{\eta}^*)\right\} = \boldsymbol{\nabla}_{\boldsymbol{\alpha}}\mathbb{E}\{\ell_3(\mathbf{W}; \boldsymbol{\gamma}^*, \boldsymbol{\delta}^*, \boldsymbol{\alpha}^*)\}/2 = \mathbf{0}, \] leading to \(\Delta_{2,\boldsymbol{\gamma}} = \Delta_{2,\boldsymbol{\delta}} = 0\) for population slopes $\boldsymbol{\alpha}^*:=\arg\min_{\boldsymbol{\alpha}\in\mathbb{R}^d}\mathbb{E}\{\ell_3(\mathbf{W};\boldsymbol{\gamma}^*,\boldsymbol{\delta}^*,\boldsymbol{\alpha})\}$ and $\boldsymbol{\beta}^*:= \arg\min_{\boldsymbol{\beta}\in\mathbb{R}^{d_1}}\mathbb{E}\{\ell_4(\mathbf{W};\boldsymbol{\gamma}^*,\boldsymbol{\delta}^*,\boldsymbol{\alpha}^*,\boldsymbol{\beta})\}$, respectively. Uniqueness of \(\boldsymbol{\gamma}^*\), \(\boldsymbol{\delta}^*\), \(\boldsymbol{\alpha}^*\), and \(\boldsymbol{\beta}^*\) are discussed in Section (ref) of the Supplementary Material. The introduced loss functions are performing (imputed) residual covariate balancing of the outcome regressions as \[ \mathbb{E}\left\{w_1 w_2' ( Y-\bar\mathbf{S}_2^\top\boldsymbol{\alpha}) \right\}=0, \qquad \mathbb{E}\left\{w_1'( Y^{\mbox{\tiny DR}}-\mathbf{S}_1^\top\boldsymbol{\beta} ) \right\} =0 \] where $w_1'=\boldsymbol{\nabla}_{\boldsymbol{\gamma}} w_1$ and $w_2'=\boldsymbol{\nabla}_{\boldsymbol{\delta}} w_2$ and $Y^{\mbox{\tiny DR}}=\bar\mathbf{S}_2^\top\boldsymbol{\alpha}+ {A_2(Y-\bar\mathbf{S}_2^\top\boldsymbol{\alpha})}/{g(\bar\mathbf{S}_2^\top\boldsymbol{\delta})}$ is the double robust imputation based of (ref). Here, we are ensuring that the (imputed) residuals are uncorrelated with the adjustments made to the second and first PS estimations.
The loss functions (ref), (ref), (ref), and (ref) are termed moment-targeting loss functions. The sequential model doubly robust (SMDR) estimator of \(\theta_{1,1}\) utilizing a cross-fitting technique can be found in Algorithm (ref).
Off-the-shelf methods cannot achieve model double robustness, even with doubly robust or Neyman orthogonal scores, for estimating average treatment effects (ATE) with a single time exposure smucler2019unifying, tan2020model, dukes2020doubly, avagyan2021high, dukes2021inference, bradic2019sparsity. Our introduced loss functions reduce to \(\ell_2 \) and \(\ell_4 \) in the single time exposure setting, but our dynamic problem introduces more complex challenges compared to the static case. We believe our work is the first to achieve model double robustness fully under Assumption (ref).
In dynamic settings, luedtke2017sequential, rotnitzky2017multiply, bradic2024high, diaz2023nonparametric employ doubly robust imputation, \( Y^{\mbox{\tiny DR}} \), to improve the rate conditions of the final estimator. However, beyond this step, they rely solely on off-the-shelf methods, such as (regularized) maximum likelihood estimation. We show that this approach alone is insufficient for robustness against model misspecification and that our newly introduced loss functions are essential. While the classical logistic loss \[ A_1\{(1-A_2)\bar\mathbf{S}_2^\top\boldsymbol{\delta} + A_2 h(\bar\mathbf{S}_2^\top\boldsymbol{\delta})\} \] with \( h(u) = -\log g(u) \) is commonly used, we propose a covariate (conditional) balancing-inspired modification, replacing \( h(\bar\mathbf{S}_2^\top\boldsymbol{\delta}) \) with \( \exp(-\bar\mathbf{S}_2^\top\boldsymbol{\delta}) \) and introducing an additional weight \( w_2 \). Furthermore, while classical least squares loss is typically applied in OR estimation, we identify the weights \( w_1 w_2' \) and \( w_1' \) as necessary to ensure score orthogonality under PS model misspecification. Our numerical experiments confirm that, in finite samples, our approach consistently achieves better bias control than existing off-the-shelf methods (see Tables (ref)-(ref)).
In the following, we choose tuning parameters $\lambda_{\boldsymbol{\gamma}}\asymp\sqrt{\log d_1/N}$, $\lambda_{\boldsymbol{\delta}}\asymp\sqrt{\log d/N}$, $\lambda_{\boldsymbol{\alpha}}\asymp\sqrt{\log d/N}$, $\lambda_{\boldsymbol{\beta}}\asymp\sqrt{\log d_1/N}$. Define $s_{\boldsymbol{\gamma}}:=\|\boldsymbol{\gamma}^*\|_0$, $s_{\boldsymbol{\delta}}:=\|\boldsymbol{\delta}^*\|_0$, $s_{\boldsymbol{\alpha}}:=\|\boldsymbol{\alpha}^*\|_0$, and $s_{\boldsymbol{\beta}}:=\|\boldsymbol{\beta}^*\|_0$ as the sparsity levels of the population nuisance parameters.
The sparsity conditions of the form $s=o\left(N/\log d\right)$ are very common in the high-dimensional statistics literature and guarantee estimation consistency. The additional condition $(s_{\boldsymbol{\gamma}}+s_{\boldsymbol{\alpha}})\log d_1\log d+s_{\boldsymbol{\delta}}\log^2d=O(N)$ is necessary since, in general, the imputed outcomes considered in the Lasso problem (ref) do not have a bounded $\psi_\alpha$-Orlicz norm. However, this condition is no longer required if we further assume that $\|\bar\mathbf{S}_2\|_\infty<C$, as in, e.g., bradic2019sparsity,tan2020model,smucler2019unifying.
The following assumption imposes some standard moment conditions where $\|X\|_{\psi_2}:=\inf\{c>0:\mathbb{E}[\psi_{2}(\lvert X\rvert/c)]\leq 1\}$, with $\psi_2(x)=\exp(x^2)-1$.
Theorem (ref) characterizes the convergence rate of the SMDR estimator. When all nuisance models are correctly specified, we have \[ \widehat{\theta}_{1,1} - \theta_{1,1} = O_p\left(\sigma N^{-1/2} + r_{\boldsymbol{\gamma}}r_{\boldsymbol{\beta}} + r_{\boldsymbol{\delta}}r_{\boldsymbol{\alpha}}\right), \] which matches the S-DRL estimator bradic2024high and outperforms the DTL estimator bradic2024high, bodory2022evaluating. Under model misspecification, DTL and S-DRL exhibit convergence rates with additional linear terms (see Table (ref)), whereas our result in (ref) involves only quadratic terms—products of nuisance estimation errors. When a nuisance model is misspecified, the convergence rate of S-DRL (and DTL) includes an additional term that is linearly dependent on the estimation error of the other nuisance model at the same exposure. In contrast, the proposed SMDR method mitigates such model misspecification errors by introducing a multiplicity factor that incorporates estimation errors from nuisance estimates constructed prior to the misspecified model. For instance, when \(\mu^*\) is misspecified, the convergence rates of DTL and S-DRL both involve \(r_{\boldsymbol{\gamma}}\), the nuisance estimation rate for \(\nu^*\). In comparison, the SMDR method reduces this term to a product of \(r_{\boldsymbol{\gamma}}\) and \(r_{\boldsymbol{\gamma}} + r_{\boldsymbol{\delta}} + r_{\boldsymbol{\alpha}}\). The sequentially designed loss functions in SMDR effectively leverage the structures of previously estimated models to downstream the impact of model misspecification when estimating subsequent models, thereby enhancing overall robustness and accuracy in the final DTE estimation.
Whenever all nuisance models are correctly specified, we have the following result.
When all nuisance functions are correctly specified, the result coincides with that of bradic2024high while also achieving the semi-parametric efficiency of bang2005doubly. Hence, we do not lose accuracy when the nuisance models are correctly specified. As shown in Theorem (ref), root-$N$ inference requires product sparsity conditions between the nuisance parameters' sparsity levels at each exposure, i.e., (ref); we name such a property as “sequential rate double robustness”. This condition is weaker than the DTL estimator where an additional product sparsity condition $s_{\boldsymbol{\gamma}}s_{\boldsymbol{\alpha}}=o(N/(\log d_1\log d))$ is imposed.
In the following, we develop theoretical properties of the proposed moment-targeted nuisance estimators $\widehat{\boldsymbol{\gamma}}$, $\widehat{\boldsymbol{\delta}}$, $\widehat{\boldsymbol{\alpha}}$, and $\widehat{\boldsymbol{\beta}}$, defined in (ref)-(ref). The analysis of nuisance estimation is non-trivial since the nuisance estimates are constructed in a sequential manner. Section (ref) shows the nuisance estimators' consistency despite potential model misspecification, while Section (ref) presents their faster consistency rates when certain models are correctly specified. Our findings show that the accuracy of nuisance models affects the estimation errors.
Our first focus is on the asymptotic behavior of moment-targeted nuisance estimators with possibly inaccurate models. In determining convergence rates, we confront the complexities of RSC conditions caused by dependent loss functions in Lemma (ref), and manage the increased gradient variability in Lemma (ref), as expanded upon in the Supplementary Material.
Among the results in Theorem (ref), part (b) is the most challenging to show. Notice that $\widehat{\boldsymbol{\delta}}$ is constructed based on a first-stage estimate $\widehat{\boldsymbol{\gamma}}$. Due to the occurrence of the imputation error $\widehat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^*$, the estimation error $\widehat{\boldsymbol{\delta}}-\boldsymbol{\delta}^*$ no longer belongs to the usual cone set $\mathbb{C}(S,k):=\{\boldsymbol{\Delta}\in\mathbb{R}^d:\|\boldsymbol{\Delta}_{S^c}\|_1\leq k\|\boldsymbol{\Delta}_S\|_1\}$. A similar problem has been recently studied by bradic2024high, where their Theorem 8 provides consistency rates of imputed Lasso estimates. The problem we consider here is even more technically challenging in that the loss function (ref) is non-quadratic with respect to $\boldsymbol{\delta}$. We consider a cone set $\widetilde{\mathbb{C}}(s,k):=\{\boldsymbol{\Delta}\in\mathbb{R}^d:\|\boldsymbol{\Delta}\|_1\leq k\sqrt{s}\|\boldsymbol{\Delta}\|_2\}$ that is “larger” than the usual $\mathbb{C}(S,k)$ and also different from the cone set studied by bradic2024high. We show that $\widehat{\boldsymbol{\delta}}-\boldsymbol{\delta}^*\in\widetilde{\mathbb{C}}(s,k)$ with high probability and some $k,s>0$; see details in Lemma (ref). Together with some empirical process results as in Lemma (ref), we control the imputation error's effect and finally reach the consistency rates introduced above; see Lemma (ref) and the proof of Theorem (ref). Although we focus on a specific loss function (ref), the results of part (b) in fact apply more broadly to other smooth and convex loss functions. Since the nuisance estimators $\widehat{\boldsymbol{\gamma}},\widehat{\boldsymbol{\delta}},\widehat{\boldsymbol{\alpha}},\widehat{\boldsymbol{\beta}}$ are constructed sequentially, and the later estimators depend on all the previous ones, the estimation errors of the nuisance parameters are cumulative, i.e., the consistency rate depends on the sparsity levels of all the nuisance parameters up to the current one.
If we have additional information that some of the nuisance models are correctly specified, we are able to achieve better consistency rates.
The new convergence rates in Theorem (ref) are established through Lemmas (ref) and (ref) of the Supplementary Material. Assuming certain nuisance models being correct, unlike Theorem (ref) and Lemma (ref), we can control the gradients involving the estimated nuisance parameters and control the imputation errors from the previous steps' nuisances in a more efficient way; see more details in Lemma (ref). As a result, we obtain faster convergence rates than Theorem (ref) given additional model correctness information.
First, we see that the subsequent estimates, with correct model specifications, lose a factor of $r_{\boldsymbol{\gamma}}$ in their rates of estimation. Secondly, specific estimates demonstrate asymptotic decoupling: (i) the convergence rate of $\widehat{\boldsymbol{\delta}}$ depends only on $s_{\boldsymbol{\delta}}$ when $\rho^*$ is correctly specified; (ii) the convergence rate of $\widehat{\boldsymbol{\alpha}}$ depends only on $s_{\boldsymbol{\alpha}}$ when $\nu^*$ is correctly specified. Lastly, the convergence of \(\widehat{\boldsymbol{\beta}}\) relies on the model correctness of \(\mu^*\), as well as the preceding models \(\rho^*\) and \(\nu^*\). Specifically, if either $\rho^*$ or $\nu^*$ is correctly specified, as explored in cases (c) and (d), the consistency rate of $\widehat{\boldsymbol{\beta}}$ depends on $s_{\boldsymbol{\beta}}$ and the sparsity level of the correctly specified model, be it $\rho^*$ or $\nu^*$. When both $\rho^*$ and $\nu^*$ are accurate, as in case (e), the consistency rate of $\widehat{\boldsymbol{\beta}}$ depends on $s_{\boldsymbol{\beta}}$ and a product sparsity $s_{\boldsymbol{\delta}}s_{\boldsymbol{\alpha}}$. When a product sparsity condition, $s_{\boldsymbol{\delta}}s_{\boldsymbol{\alpha}}=o(N/\log^2d)$, is assumed as in (ref) of Theorem (ref), $\widehat{\boldsymbol{\beta}}$ also becomes asymptotically decoupled from the other three estimates.
We illustrate the finite sample properties of the introduced estimator on a number of simulated experiments. We focus on the estimation of $\theta=\theta_a-\theta_{a'}$ where $a=(a_1,a_2)=(1,1)$ and $a'=(a_1',a_2')=(0,0)$. We describe the considered data generating processes below. The outcome variables are generated as $Y_i=A_{1i}A_{2i}Y_i(1,1)+(1-A_{1i})(1-A_{2i})Y_i(0,0)$.
Setting (a): Non-linear $\mu$ and non-logistic $\rho$. Generate covariates at the first exposure: for each $i\leq N$, $\mathbf{S}_{1i}\sim^\mathrm{iid} N_{d_1}(\mathbf{0},\mathbf{I}_{d_1}).$ The treatment indicators of the first exposure are generated as $A_{1i}\mid\mathbf{S}_{1i}\sim\mathrm{Bernoulli}(g(\mathbf{S}_{1i}^\top\boldsymbol{\gamma}))$. A $d$ dimensional vector of all ones and zeros are denoted with $\mathbf{1}_{(d)}$ and $\mathbf{0}_{(d)}$, respectively. Covariates at the second exposure satisfy $\mathbf{S}_{2i}=0.5Q(A_{1i})(\mathbf{S}_{1i}^2-1)+Q(A_{1i})\mathbf{S}_{1i}+A_{i}(1+\delta_{1i})\mathbf{1}_{(d_2)}+\boldsymbol{\delta}_{1i},$ where $\mathbf{S}_{1i}^2\in\mathbb{R}^{d_1}$ is the coordinate-wise square of $\mathbf{S}_{1i}$, $\boldsymbol{\delta}_{1i}\sim^\mathrm{iid} N_{d_2}(0,\mathbf{I}_{d_2})$, and a matrix $Q$ is defined with $\{Q(1)\}_{i,j}=0.8^{|i-j|}\mathbbm1\{|i-j|\leq1\}$ and $\{Q(0)\}_{i,j}=0.7^{|i-j|}\mathbbm1\{|i-j|\leq2\}$ for $i\leq d_2$ and $j\leq d_1$. The treatment indicators at the second exposure are generated as $A_{2i}\mid(\bar\mathbf{S}_{2i},A_{1i})\sim\mathrm{Bernoulli}(A_{1i}\tilde g(\bar\mathbf{S}_{2i}^\top\boldsymbol{\delta})+(1-A_{1i})\tilde g(-\bar\mathbf{S}_{2i}^\top\boldsymbol{\delta}))$, where $\tilde g(u):=(|u+1|+0.1)/(|u+1|+1)$. Lastly, $Y_i(1,1)=\bar\mathbf{S}_{2i}^\top\boldsymbol{\alpha}+1+\epsilon_i$, $Y_i(0,0)=-\bar\mathbf{S}_{2i}^\top\boldsymbol{\alpha}-1+\epsilon_i$ and $\epsilon_i\sim^\mathrm{iid}N(0,1)$. We consider $\boldsymbol{\alpha}=(1,\mathbf0_{(d_1-1)},0.5,0.5,0.5,0.5,\mathbf0_{(d_2-4)})^\top$, $\boldsymbol{\gamma}=(1,1,\mathbf0_{(d_1-2)})^\top$ and $\boldsymbol{\delta}=(1,\mathbf0_{(d_1-1)},0.5,0.5,0.5,0.5,\mathbf0_{(d_2-4)})^\top$.
Setting (b): Non-linear $\mu$ and non-linear $\nu$. At the first exposure, generate covariates from a centered Beta distribution, i.e., $\mathbf{S}_{1ij}\sim^\mathrm{iid} \mathrm{Beta}(1,2)-1/3$ for each $i\leq N$ and $j\leq d_1$; generate $A_{1i}\mid\mathbf{S}_{1i}\sim\mathrm{Bernoulli}(g(\mathbf{S}_{1i}^\top\boldsymbol{\gamma}))$. At the second exposure, generate $\mathbf{S}_{2i}=W(A_{1i})\mathbf{S}_{1i}+A_{2i}\mathbf{1}_{(d_2)}+\boldsymbol{\delta}_i$ and $A_{2i}\mid(\bar\mathbf{S}_{2i},A_{1i})\sim\mathrm{Bernoulli}(A_{1i}g(\bar\mathbf{S}_{2i}^\top\boldsymbol{\delta})+(1-A_{1i})g(-\bar\mathbf{S}_{2i}^\top\boldsymbol{\delta}))$, where $\boldsymbol{\delta}_{ij}\sim^\mathrm{iid} \mathrm{Beta}(1,4)-1/5$, $\{W(1)\}_{i,j}=0.2^{|i-j|}\mathbbm1\{|i-j|\leq1\}$ and $\{W(0)\}_{i,j}=0.2^{|i-j|}\mathbbm1\{|i-j|\leq2\}+0.1\mathbbm1\{|i-j|=2\}$ for each $i\leq d_2$ and $j\leq d_1$. Here, $Y_i(1,1)=\bar\mathbf{S}_{2i}^\top\boldsymbol{\alpha}-1+2r_i+\epsilon_i$, $Y_i(0,0)=-\bar\mathbf{S}_{2i}^\top\boldsymbol{\alpha}+1-2r_i+\epsilon_i$ and $\epsilon_i\sim^\mathrm{iid}N(0,1)$. Here, we consider non-linear signals with $r_i$ as the standardized version of $\mathbf{S}_{1i1}\mathbf{S}_{1i2}\mathbbm1\{\mathbf{S}_{1i2}>0.3\}+\mathbf{S}_{1i1}\mathbf{S}_{1i3}\mathbbm1\{\mathbf{S}_{1i1}>0.3\}+\mathbf{S}_{1i2}\mathbf{S}_{1i3}\mathbbm1\{\mathbf{S}_{1i1}>0.3\}$. The parameters are $\boldsymbol{\alpha}=(-1,0,0,1/18,\mathbf0_{(d_1-4)},-1,-1,-1,\mathbf0_{(d_2-3)})^\top$, $\boldsymbol{\gamma}=(1,1,\mathbf0_{(d_1-2)})^\top$ and $\boldsymbol{\delta}=(-2,-2,\mathbf0_{(d_1+d_2-2)})^\top$.
For each setting, we consider dimensions $d_1=100$ and $d_2=50$ (resulting in $d=d_1+d_2=150$), with total sample sizes $N$ ranging from $400$ to $16,000$. The experiments are repeated 200 times. Our proposed SMDR estimator is denoted as SMDR1 (see Algorithm (ref) with $\mathbb{K}=5$). Additionally, we present a slightly modified version, SMDR2, which constructs all the nuisances on the entire sub-sample of $\mathcal{I}_{-k}$ in Steps 4-7 of Algorithm (ref).
For comparison, we include several existing estimators: the inverse probability weighting (IPW) estimator, where propensity score models are estimated using $\ell_1$-regularized logistic regression without cross-fitting; the sequential doubly robust (SDR) estimator by Luedtke et al. (2017), where nuisance functions are estimated through linear and logistic regression without $\ell_1$-regularization; a cross-fitted version of the SDR estimator rotnitzky2017multiply,diaz2023nonparametric using random forest nuisance estimates, denoted as SDR-RF; the sequential doubly robust Lasso (S-DRL) estimator proposed by Bradic et al. (2021); and two versions of the dynamic treatment Lasso (DTL) estimator proposed by Bradic et al. (2021) and Bodory et al. (2022), denoted as DTL2 and DTL1, respectively. DTL2's nuisances are estimated using samples in $\mathcal{I}_{-k}$, while DTL1's nuisances use different sub-samples in $\mathcal{I}_{\boldsymbol{\gamma}}$, $\mathcal{I}_{\boldsymbol{\delta}}$, $\mathcal{I}_{\boldsymbol{\alpha}}$, and $\mathcal{I}_{\boldsymbol{\beta}}$. Here, DTL1 and SMDR1 share the same type of sample splitting, while DTL2 and SMDR2 share the same type of sample splitting. The tuning parameters are chosen through 5-fold cross-validations. Additionally, we report the performance of a naive empirical difference estimator (empdiff), $\widehat{\theta}_\mathrm{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})$, as well as an oracle doubly robust estimator, $\widehat{\theta}_\mathrm{oracle}$, which uses the doubly robust score with correct nuisance functions. The results are reported in Tables (ref)-(ref).
Due to the confounding factors, the naive empirical difference estimator $\widehat{\theta}_\mathrm{empdiff}$ is not consistent with large biases and poor coverage; see Tables (ref)-(ref). The IPW estimator also has very large biases (especially when $N$ is large) and provides bad coverage results under Setting (a), where the PS model at the second exposure is misspecified. In Setting (b) where both PS models are correctly specified, surprisingly, the IPW estimator provides acceptable coverages although there is no theoretical guarantees from existing work in high dimensions. However, the RMSEs of IPW are comparable with SMDR2 only when $N=12000$, and worse than SMDR2 in other cases; see Table (ref).
In scenarios where the sample size is relatively small compared to the dimensionality, the SDR estimator exhibits substantial estimation errors due to the absence of regularization in the nuisance estimation process. As the sample size increases, the SDR estimator tends to yield more acceptable estimation and inference results; however, its efficiency remains notably inferior to that of the proposed SMDR1 and SMDR2 estimators. On the other hand, the SDR-RF estimator typically yields satisfactory results in scenarios with relatively small sample sizes. However, as the sample size increases, the inferential outcomes, as well as the estimation results in Setting (a), tend to deteriorate. This phenomenon occurs due to the relatively slow convergence rates of random forests for nuisance estimation, leading to notable biases that cannot be ignored and consequently compromising the accuracy of inference.
The DTL1 estimator exhibits relatively poor performance overall, with biases often close to RMSE, and coverages far below the desired $95\%$. This suboptimal performance arises from two main factors: (i) The DTL estimators are only proven to be consistent when model misspecification occurs bradic2024high and are not necessarily $\sqrt N$-consistent nor asymptotically normal; (ii) The sample splitting method of DTL1 is inefficient in finite samples, as only $1/5$ of the samples are used to obtain each nuisance estimator when $\mathbb{K}=5$. DTL2 is constructed using a more efficient sample splitting, leading to smaller biases than DTL1. However, it fails to achieve satisfactory coverage guarantees even with a large sample size.
The S-DRL estimator is constructed similarly to DTL2, except with a different doubly robust estimation strategy for the first OR model. When the sample size is large enough, the S-DRL method provides RMSEs similar to (see Table (ref)) or smaller than (see Table (ref)) the DTL2 estimator. However, the coverages based on S-DRL remain below the desired $95\%$, even with a large sample size. Interestingly, the DTL and S-DRL methods yield coverages closer to $95\%$ when the sample size $N$ is small; however, an increase in the sample size does not lead to better coverage. This is due to the biases of these methods decaying slower than the parametric rate under model misspecification, making the normal approximation inaccurate.
For sufficiently large sample sizes, the proposed SMDR1 estimator consistently outperforms IPW, SDR, SDR-RF, DTL1, DTL2, and S-DRL in terms of estimation, exhibiting smaller biases and RMSEs across all considered settings (refer to Tables (ref)-(ref)). Moreover, SMDR1 provides satisfactory coverages that closely approach the desired $95\%$. However, its performance can be occasionally inferior to that of SDR-RF, DTL2, and S-DRL when the total sample size $N$ is small, as observed in Table (ref) for $N\in\{400,1000\}$. Notably, SMDR1 consistently outperforms the DTL1 method, which utilizes the same type of sample splitting. This observation suggests the inadequacy of the sample splitting technique introduced in SMDR1. Indeed, when $N=1000$, only approximately $N_0\approx341$ samples are observed with the treatment path $(1,1)$ under Setting (a). Consequently, Step 4 of Algorithm (ref) results in only $N_0(\mathbb{K}-1)/(4\mathbb{K})=N_0/5\approx68$ training samples for nuisance estimation, while the nuisance parameters have dimensions $d_1=100$ and $d=150$ for the first and second exposures, respectively. In contrast, the SMDR2 method employs $N_0(\mathbb{K}-1)/\mathbb{K}=0.8N_0\approx273$ training samples under the same setup. Notably, the SMDR2 method yields more stable results when $N$ is small, especially under Setting (a); see Table (ref). Therefore, while we acknowledge that the SMDR2 estimator may demand more stringent sparsity conditions compared to SMDR1 from a theoretical perspective, we recommend also considering the more efficient sample splitting technique of SMDR2, particularly when the sample size is small.
In this section, we compare the estimation and inference performance of the DTE estimators through semi-synthetic experiments. We consider a dataset from the National Job Corps Study, which is the largest and most comprehensive job training program in the US established in 1964, and serves approximately 50,000 disadvantaged youths aged 16-24 each year by providing vocational training and academic education. A detailed description of the original design and main effects can be found in schochet2008does,schochet2001national.
We consider a dataset of 11,313 individuals, with 6,828 assigned to the Job Corps and 4,485 not. Treatments, denoted as $Z_{ti}\in\{0,1,2,3\}$ ($t\in{1,2}$), are assigned to the $i$th individual in the first and second years after the initial randomization, where $Z_{ti}=0$ represents non-enrollment, $Z_{ti}=1$ enrollment without program participation, $Z_{ti}=2$ high-school-level education, and $Z_{ti}=3$ vocational training. The baseline covariate vector, $\mathbf{S}_{1i}$, has 909 characteristics, while $\mathbf{S}_{2i}$ includes 1,427 characteristics that are evaluated before the second-year treatment assignment. We exclude 2,610 individuals whose treatment stages are missing completely at random Schochet2003national, resulting in a final sample of 8,703 individuals. We also exclude the binary characteristics, if the 0/1 groups are extremely unbalanced in that the minority group's size is less than 10 within the 8703 individuals, resulting in the final $\mathbf{S}_{1i}$ with 891 characteristics and $\mathbf{S}_{2i}$ with 1350 characteristics. After standardizing the covariates, we generate the potential outcomes $\widetilde{Y}_i(z)$ corresponding to treatment paths $z=(z_1,z_2)\in\{0,1,2,3\}^2$ based on Settings (a) and (b) below. The observed outcome is generated as $Y_i=Y_i(Z_{i1},Z_{i2})=\sum_{z\in\{0,1,2,3\}^2}\mathbbm1_{\{(Z_{1i},Z_{2i})=z\}}Y_i(z)$.
We consider estimation of the DTE, $\theta=E\{\widetilde{Y}_i(z)\}-E\{\widetilde{Y}_i(z')\}$, focusing on the treatment path $z=(z_1,z_2)=(3,3)$ and the control path $z'=(z_1',z_2')=(1,1)$. To estimate the expected potential outcome $E\{\widetilde{Y}_i(z)\}$, we set $A_{1i}=\mathbbm1_{\{Z_{1i}=z_1\}}$, $A_{2i}=\mathbbm1_{\{Z_{2i}=z_2\}}$, and $Y_i(1,1)=\widetilde{Y}_i(z)$. Then $E\{\widetilde{Y}_i(z)\}=E\{Y_i(1,1)\}$ can be estimated using Algorithm (ref). The control arm $E\{\widetilde{Y}_i(z')\}=E\{Y_i(0,0)\}$ can be estimated analogously, and the final DTE estimator is constructed as the difference of the obtained estimates. Let $\epsilon_i\sim^\mathrm{iid}N(0,1)$. The potential outcomes $Y_i(1,1)=\widetilde{Y}_i(z)$ and $Y_i(0,0)=\widetilde{Y}_i(z')$ are generated as below.
Setting (a): Linear $\nu$. Let $Y_i(1,1)=\bar\mathbf{S}_{2i}^\top\boldsymbol{\alpha}+\epsilon_i$ and $Y_i(0,0)=-\bar\mathbf{S}_{2i}^\top\boldsymbol{\alpha}+\epsilon_i$, where $\boldsymbol{\alpha}=0.5\cdot(\alpha_0,\mathbf{1}_{(8)},\mathbf{0}_{(d_1-8)}, \mathbf{1}_{(4)},\mathbf0_{(d_2-4)})^\top$ with $\alpha_0$ varying from 0.15 to 0.3.
Setting (b): Non-linear $\nu$. Let $Y_i(1,1)=\mathbf{S}_{1i}^\top\boldsymbol{\alpha}_1+ (\mathbf{S}_{2i}^2-1)^\top\boldsymbol{\alpha}_2+\epsilon_i$ and $Y_i(0,0)=-\mathbf{S}_{1i}^\top\boldsymbol{\alpha}_1-(\mathbf{S}_{2i}^2-1)^\top\boldsymbol{\alpha}_2+\epsilon_i$, where $\mathbf{S}_{2i}^2 $ is the coordinate-wise square of $\mathbf{S}_{2i}$, $\boldsymbol{\alpha}_1=0.5\cdot(\alpha_0,\mathbf{1}_{(8)},\mathbf{0}_{(d_1-8)})^\top$, $\boldsymbol{\alpha}_2=0.05\cdot (\mathbf{1}_{(4)},\mathbf0_{(d_2-4)})^\top$, and $\alpha_0$ varies from 0.25 to 0.5.
For each setting, we implement the DTL2, S-DRL, and the proposed SMDR1 estimators (see Section (ref)). The results are reported in Table (ref), where the biases are calculated based on the oracle difference-in-mean estimate $\widehat{\theta}_O:=N^{-1}\sum_{i=1}^N\{Y_i(1,1)-Y_i(0,0)\}=\alpha_0$. Recall that under our simulated outcome setting, we get to see all potential outcomes: two per individual. Under both Settings (a) and (b), the proposed SMDR1 method provides smaller absolute biases than the DTL2 and S-DRL estimators. In addition, under Setting (a), where the OR model at the second exposure is truly linear, all the constructed confidence intervals contain the oracle estimate $\widehat{\theta}_O$. However, when the potential outcome is generated through a quadratic function (under Setting (b)), the oracle estimate $\widehat{\theta}_O$ does not lie in the confidence intervals based on the DTL2 and S-DRL methods; on the other hand, the proposed SMDR1 method leads to confidence intervals containing the oracle estimate. Moreover, considering the hypothesis testing problem with the null $H_0:\theta=0$ and the alternative $H_1:\theta\neq0$, the reported p-values decay as $\widehat{\theta}_O=\alpha_0$ grows; see Figure (ref). When $\alpha_0$ is large enough, all the methods return p-values smaller than $0.05$; however, different methods require different signal levels to detect the causal effect and reject the null successfully. Under Setting (a), the proposed SMDR1 method is able to detect the causal effect with a significance level of $95\%$ when $\alpha_0=0.2$; however, under the same signal level, both the DTL2 and S-DRL methods fail to reject the null as the corresponding p-values are larger than $0.05$. Similarly, under Setting (b) with $\widehat{\theta}_O=\alpha_0=0.3$, the proposed SMDR1 method is able to detect the causal effect, whereas the p-value based on the DTL2 and S-DRL methods are both very large. Therefore, we observed a significantly better power in the SMDR1 method than both DTL2 and S-DRL.
This paper introduces new techniques to enable statistical inference for treatment effects in dynamic, high-dimensional, and potentially misspecified settings. By proposing a set of novel loss functions for nuisance models, we develop a sequential model doubly robust (SMDR) method that achieves root-\(N\) inference under minimal requirements. Our findings highlight the critical role of nuisance model estimation—naive, off-the-shelf estimators fail to achieve the desired robustness, even within doubly robust frameworks. While some nuisance models can be estimated independently, our results demonstrate that adopting a sequential estimation approach with nested designs significantly reduces the final estimation error for causal parameters. This observation raises an intriguing question: does this phenomenon persist in other statistical estimation problems, particularly in complex longitudinal settings requiring multi-stage estimation?
In the context of dynamic treatment regimes, a related but distinct doubly robust (DR) property has been explored. Existing methods for consistently estimating the optimal regime often require correctly specified contrast models at all later stages of estimation schulte2014q, shi2018high. This condition is highly restrictive, especially in settings with multiple exposure occasions, and is not required in our framework. A natural question arises: How should decisions be made if contrast models cannot be accurately specified at later stages? Our results, outlined in Theorem (ref), suggest that outcome regression (OR) models, such as \(\mathbb{E}\{Y(a_1,a_2) \mid \mathbf{S}_1 = \mathbf{s}_1\}\), can still be consistently estimated under these conditions. Consequently, optimizing the estimated OR functions at the first exposure over possible treatment paths offers a conservative yet viable strategy for newly arriving individuals. Notably, this approach eliminates the need for additional covariate evaluations at later stages, making it particularly useful in applications where accessing longitudinal covariates is costly or impractical.
The proposed algorithms can also be implemented using linear or logistic forms with basis functions, such as B-splines. However, further theoretical analysis is required for such approaches, as well as for the application of other non-parametric methods, including random forests and boosting. Additionally, while our methods leverage sparse structures in the models, future research should explore strategies that accommodate dense models with robust guarantees, broadening the applicability of our framework.
Sections (ref)-(ref) contain additional discussions, justifications, and proofs of the main results. Additional notations used in the supplementary material are introduced in Section (ref). Section (ref) discusses the uniqueness of the moment-targeted parameters; the justification of their identification is provided in Section (ref). We introduce some useful auxiliary lemmas in Section (ref). The proofs of main results and auxiliary lemmas are in Sections (ref) and (ref), respectively.