Dynamic treatment effects: high-dimensional inference under model misspecification

The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.

94,718 characters

Dynamic treatment effects: high-dimensional doubly robust inference under model misspecification



\date{}

\title{\bf Dynamic treatment effects: high-dimensional doubly robust inference under model misspecification}

\author{Yuqian Zhang\thanks{Institute of Statistics and Big Data, Renmin University of China} \and Weijie Ji\thanks{School of Statistics and Management, Shanghai University of Finance and Economics} \and Jelena Bradic\thanks{Department of Mathematics and Halicioglu Data Science Institute, University of California, San Diego, E-mail: [email removed]}}

\maketitle

\begin{abstract}
Estimating dynamic treatment effects is crucial across various disciplines, providing insights into the time-dependent causal impact of interventions. However, this estimation poses challenges due to time-varying confounding, leading to potentially biased estimates. Furthermore, accurately specifying the growing number of treatment assignments and outcome models with multiple exposures appears increasingly challenging to accomplish. Double robustness, which permits model misspecification, holds great value in addressing these challenges. This paper introduces a novel “sequential model doubly robust” estimator. We develop novel moment-targeting estimates to account for confounding effects and establish that root-$N$ inference can be achieved as long as at least one nuisance model is correctly specified at each exposure time, despite the presence of high-dimensional covariates. Although the nuisance estimates themselves do not achieve root-\(N\) rates, the carefully designed loss functions in our framework ensure final root-\(N\) inference for the causal parameter of interest. Unlike off-the-shelf high-dimensional methods, which fail to deliver robust inference under model misspecification even within the doubly robust framework, our newly developed loss functions address this limitation effectively.
\end{abstract}

\begin{keyword}
Causal Inference, High-dimensional Statistics, Robust Inference, Longitudinal Data
\end{keyword}

\section{Introduction}\label{sec:intro}
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 \citep{babino2019multiple}. We demonstrate that this challenge can be overcome by introducing new estimates of the underlying causal effects.

\subsection{Estimation of causal effects under dynamic setups}

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{fig:diag}.



\begin{figure}[h]
	\scriptsize
	\centering
	\begin{tikzpicture}
		[
		roundnode/.style={circle, draw=red!20, fill=red!5, very thick, minimum size=7mm},
		circlenode/.style={circle, draw=blue!20, fill=blue!5, very thick, minimum size=7mm},
		];
		\node[circlenode] (s1) at (0,0) {$\mathbf{S}_1$};
		\node[roundnode] (a1) at (1.5,2) {$A_1$};
		\node[circlenode] (s2) at (3,0.5) {$\mathbf{S}_2$};
		\node[roundnode] (a2) at (4.5,2) {$A_2$};
		\node[circlenode] (y) at (6,0) {$Y$};
		\draw[-{Stealth[length=2mm, width=1mm]}] (s1) -- (a1);
		\draw[-{Stealth[length=2mm, width=1mm]}] (a1) -- (s2);
		\draw[-{Stealth[length=2mm, width=1mm]}] (s1) -- (s2);
		\draw[-{Stealth[length=2mm, width=1mm]}] (s2) -- (a2);
		\draw[-{Stealth[length=2mm, width=1mm]}] (s1) -- (a2);
		\draw[-{Stealth[length=2mm, width=1mm]}] (a1) -- (a2);
		\draw[-{Stealth[length=2mm, width=1mm]}] (s1) -- (y);
		\draw[-{Stealth[length=2mm, width=1mm]}] (a1) -- (y);
		\draw[-{Stealth[length=2mm, width=1mm]}] (s2) -- (y);
		\draw[-{Stealth[length=2mm, width=1mm]}] (a2) -- (y);
	\end{tikzpicture}
	\caption{\centering Causal diagrams for dynamic settings with two exposure occasions.}
	\label{fig:diag}
\end{figure}



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 \citep{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.
\begin{assumption} \label{cond:basic}
(a) Sequential ignorability:
$
Y(a_1,a_2)\protect\mathpalette{\protect\independenT}{\perp} A_1\mid\mathbf{S}_1$, $Y(a_1,a_2)\protect\mathpalette{\protect\independenT}{\perp} A_2\mid(\mathbf{S}_1,\mathbf{S}_2,A_1=a_1).
$
(b) Consistency:
$Y=Y(A_1,A_2).$
(c) Overlap:
let $\mathbb{P}(c_0<\pi(\mathbf{S}_1)<1-c_0)=1,\ \mathbb{P}(c_0<\rho(\bar\mathbf{S}_2)<1-c_0)=1$ with some constant $c_0\in(0,1)$. Additionally, let $\pi^*$ and $\rho^*$ be some functions satisfying $\mathbb{P}(c_0<\pi^*(\mathbf{S}_1)<1-c_0)=1,\ \mathbb{P}(c_0<\rho^*(\bar\mathbf{S}_2)<1-c_0)=1.$
\end{assumption}

In the standard conditions above (see, e.g., \cite{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 \cite{murphy2001marginal, bang2005doubly, yu2006double}, where
\begin{equation}\label{def:score}
\theta_{1,1}=\mathbb{E}\{\psi(\mathbf{W};\boldsymbol{\eta}^*)\}, \ \psi(\mathbf{W};\boldsymbol{\eta}^*):=\mu^*(\mathbf{S}_1)+\frac{A_1\{\nu^*(\bar\mathbf{S}_2)-\mu^*(\mathbf{S}_1)\}}{\pi^*(\mathbf{S}_1)}+\frac{A_1A_2\{Y-\nu^*(\bar\mathbf{S}_2)\}}{\pi^*(\mathbf{S}_1)\rho^*(\bar\mathbf{S}_2)},
\end{equation}
as long as Assumption \ref{cond:mis} holds.
\begin{assumption}[Sequential model double robustness]\label{cond:mis}
Let (a) either $\pi=\pi^*$ or $\mu=\mu^*$ hold, but not necessarily both; and (b) either $\rho=\rho^*$ or $\nu=\nu^*$, but not necessarily both.
\end{assumption}

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 \cite{murphy2001marginal}, the correctness of \(\mu\) hinges on the correctness of \(\nu\). We show that the above double-robust representation is not sufficient to guarantee \emph{model double robustness}. A method is deemed \emph{doubly robust} if it yields consistent estimates under Assumption \ref{cond:mis}.
In contrast, a method is considered \emph{model doubly robust} (or said to provide \emph{doubly robust inference}) if its inference remains valid under the same Assumption \ref{cond:mis}, as long as the (possibly misspecified) working models can be estimated with $o(N^{-1/4})$ rates \citep{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{cond:mis}.

\subsection{Existing results}

We illustrate the existing results under the following four settings all encompassed within Assumption \ref{cond:mis}:
\vskip -10pt
\begingroup
\setlength{\abovedisplayskip}{0pt}
\setlength{\belowdisplayskip}{0pt}
\setlength{\abovedisplayshortskip}{0pt}
\setlength{\belowdisplayshortskip}{0pt}
\begin{align}
&\text{Only the OR models are correctly parametrized;}\label{CAN_a}
\end{align}
\begin{align}
&\text{Only the PS models are correctly parametrized;}\label{CAN_b}
\end{align}
\begin{align}
&\text{Only the first OR and the second PS models are correctly parametrized;}\label{CAN_c}
\end{align}
\begin{align}
&\text{Only the second OR and first the PS models are correctly parametrized.}\label{CAN_d}
\end{align}
\endgroup


\noindent In low dimensions, inverse probability weighting (IPW) \citep{robins1986new,hernan2001marginal,robins2004optimal} provide valid inference allowing only \eqref{CAN_b}. On the other hand, covariate balancing methods \citep{kallus2018optimal,yiu2018covariate} only allow for \eqref{CAN_a}. Doubly robust methods \citep{bang2005doubly, yu2006double,orellana2010dynamic} allow \eqref{CAN_a} or \eqref{CAN_b}, but do not accommodate \eqref{CAN_c} and \eqref{CAN_d}. The multiple robust estimator proposed by \cite{babino2019multiple} allows for \eqref{CAN_a}, \eqref{CAN_b}, or \eqref{CAN_c}, but does not accommodate \eqref{CAN_d}. Using the following double-robust imputation step,
\begin{equation}\label{rep:DR-mu}
\mu(\mathbf{s}_1) \;=\; \mathbb{E}\!\Biggl[\nu^*(\bar{\mathbf{S}}_2) \;+\; \frac{A_2 \bigl\{Y - \nu^*(\bar{\mathbf{S}}_2)\bigr\}}{\rho^*(\bar{\mathbf{S}}_2)} \;\mid\; \mathbf{S}_1 = \mathbf{s}_1, A_1 = 1\Biggr],
\end{equation}
\citet{luedtke2017sequential} propose a consistent estimator and \citet{rotnitzky2017multiply} develop an asymptotically normal estimator of the DTE under conditions \eqref{CAN_a}--\eqref{CAN_d}, requiring that the nuisance estimates belong to a Donsker class \citep{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. \cite{rotnitzky2017multiply,bodory2022evaluating,diaz2023nonparametric} adopt cross-fitting techniques and double machine learning \citep{chernozhukov2017double} to relax the need for Donsker conditions, thereby allowing more flexible, nonparametric methods. However, valid inference still requires \emph{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 \citep{Robins1997causal}, \citet{lewis2021double} achieve inferential guarantees only when the “blip functions”—differences in outcome regression across treatment paths—are low-dimensional and correctly specified. \citet{viviano2021dynamic} require both OR models to be correct, covering only \eqref{CAN_a}. Dynamic Treatment Lasso (DTL) \citep{bodory2022evaluating} achieves \emph{consistency} under \eqref{CAN_a}, \eqref{CAN_b}, and \eqref{CAN_d}, and Sequential Doubly Robust Lasso (S-DRL) \citep{bradic2024high} extends this to \eqref{CAN_c} via \eqref{rep:DR-mu}. However, both DTL and S-DRL require all nuisance models to be correctly specified for valid \emph{inference}, thus excluding misspecification in \eqref{CAN_a}--\eqref{CAN_d}.

\subsection{Addressing misspecification: direct bias control}

All the doubly robust methods discussed above share a common limitation: while estimators based on the doubly robust score \eqref{def:score} reduce bias to product (quadratic) terms of the nuisance estimation errors, this reduction critically depends on correctly specifying \emph{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{cond:mis}.
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 \citep{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,
\begin{equation}\label{eq:taylor}
\widehat{\theta}_{1,1}-\theta_{1,1}=\Delta_1+\Delta_2+\Delta_3.
\end{equation}
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{cond:mis} 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 \citep{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 \emph{moment-targeted estimates} that directly guarantee the following moment condition
\begin{equation}
\mathbb{E}\{\nabla_{\boldsymbol{\eta}}\psi(\mathbf{W};\boldsymbol{\eta}^*)\}=\mathbf{0}, \;\; \text{even under model misspecification}. \label{eq:SMDR}
\end{equation}
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 \eqref{def:score}, 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{table:rate}, this approach transforms linear bias into quadratic, achieving faster convergence rates and robust inference even under model misspecification; see Table \ref{table:sparsity}.
Our approach applies to all cases \eqref{CAN_a}-\eqref{CAN_d}.




\section{The sequential model doubly robust estimator}\label{sec:DR-DTE}

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
\begin{align}
 \mathbb{E}\{\boldsymbol{\nabla}_{\boldsymbol{\beta}}\psi(\mathbf{W}; \boldsymbol{\eta}^*)\} &= \mathbb{E}\left[\left\{1 - \frac{A_1}{g(\mathbf{S}_1^\top \boldsymbol{\gamma}^*)}\right\} \mathbf{S}_1\right],\label{eq:moment-beta}
 \\
 \mathbb{E}\{\boldsymbol{\nabla}_{\boldsymbol{\alpha}}\psi(\mathbf{W}; \boldsymbol{\eta}^*)\} &= \mathbb{E}\left[\frac{A_1}{g(\mathbf{S}_1^\top \boldsymbol{\gamma}^*)}\left\{1 - \frac{A_2}{g(\bar{\mathbf{S}}_2^\top \boldsymbol{\delta}^*)}\right\} \bar{\mathbf{S}}_2\right],\label{eq:moment-alpha}\\
 \mathbb{E}\{\boldsymbol{\nabla}_{\boldsymbol{\delta}}\psi(\mathbf{W}; \boldsymbol{\eta}^*)\} &= \mathbb{E}\left\{-\frac{A_1 A_2 \exp(-\bar{\mathbf{S}}_2^\top \boldsymbol{\delta}^*) (Y - \bar{\mathbf{S}}_2^\top \boldsymbol{\alpha}^*)}{g(\mathbf{S}_1^\top \boldsymbol{\gamma}^*)} \bar{\mathbf{S}}_2\right\},\nonumber
 \\
 \mathbb{E}\{\boldsymbol{\nabla}_{\boldsymbol{\gamma}}\psi(\mathbf{W}; \boldsymbol{\eta}^*)\} &= \mathbb{E}\left[-A_1 \exp(-\mathbf{S}_1^\top \boldsymbol{\gamma}^*)\left\{\frac{A_2 (Y - \bar{\mathbf{S}}_2^\top \boldsymbol{\alpha}^*)}{g(\bar{\mathbf{S}}_2^\top \boldsymbol{\delta}^*)} + \bar{\mathbf{S}}_2^\top \boldsymbol{\alpha}^* - \mathbf{S}_1^\top \boldsymbol{\beta}^*\right\} \mathbf{S}_1\right].\nonumber
\end{align}
The above equations are listed in an order such that the right-hand side of each equation involves progressively more nuisance parameters. For instance, \eqref{eq:moment-beta} only involves $\boldsymbol{\gamma}^*$, while \eqref{eq:moment-alpha} 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{alg:BRDR}. We define the estimate for the first PS, \(\pi^*(\mathbf{s}_1)\), with \(\lambda_{\boldsymbol{\gamma}} > 0\), as
\begin{align}
\widehat{\boldsymbol{\gamma}} & := \arg\min_{\boldsymbol{\gamma} \in \mathbb{R}^{d_1}} \left\{ \lvert\mathcal{I}_{\boldsymbol{\gamma}}\rvert^{-1} \sum_{i \in \mathcal{I}_{\boldsymbol{\gamma}}} \ell_1(\mathbf{W}_i; \boldsymbol{\gamma}) + \lambda_{\boldsymbol{\gamma}} \|\boldsymbol{\gamma}\|_1 \right\},\;\;\mbox{where}\label{def:alphahat}\\
\ell_1(\mathbf{W}; \boldsymbol{\gamma}) & := (1 - A_1) \mathbf{S}_1^\top \boldsymbol{\gamma} + A_1 \exp(-\mathbf{S}_1^\top \boldsymbol{\gamma}). \label{def:alpha}
\end{align}

The loss function \eqref{def:alpha} 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 \citep{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{sec:exist_unique} 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
\begin{align}
\widehat{\boldsymbol{\delta}}=\widehat{\boldsymbol{\delta}}(\widehat{\boldsymbol{\gamma}})~:=~&\arg\min_{\boldsymbol{\delta}\in\mathbb{R}^{d}}\biggr\{M^{-1}\sum_{i\in\mathcal I_{\boldsymbol{\delta}}}\ell_2(\mathbf{W}_i;\widehat{\boldsymbol{\gamma}},\boldsymbol{\delta})+\lambda_{\boldsymbol{\delta}}\|\boldsymbol{\delta}\|_1\biggr\},\;\;\mbox{where}\label{def:betahat}\\
\ell_2(\mathbf{W};\boldsymbol{\gamma},\boldsymbol{\delta})~:=~&\frac{A_1}{g(\mathbf{S}_1^\top\boldsymbol{\gamma})}\left\{(1-A_2)\bar\mathbf{S}_2^\top\boldsymbol{\delta}+A_2\exp(-\bar\mathbf{S}_2^\top\boldsymbol{\delta})\right\}.\label{def:l2}
\end{align}
 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
\begin{align}
\widehat{\boldsymbol{\alpha}}=\widehat{\boldsymbol{\alpha}}(\widehat{\boldsymbol{\gamma}},\widehat{\boldsymbol{\delta}})~:=~&\arg\min_{\boldsymbol{\alpha}\in\mathbb{R}^d}\left\{M^{-1}\sum_{i\in\mathcal I_{\boldsymbol{\alpha}}}\ell_3(\mathbf{W}_i;\widehat{\boldsymbol{\gamma}},\widehat{\boldsymbol{\delta}},\boldsymbol{\alpha})+\lambda_{\boldsymbol{\alpha}}\|\boldsymbol{\alpha}\|_1\right\},\label{def:gammahat}
\\
\widehat{\boldsymbol{\beta}}=\widehat{\boldsymbol{\beta}}(\widehat{\boldsymbol{\gamma}},\widehat{\boldsymbol{\delta}},\widehat{\boldsymbol{\alpha}})~:=~&\arg\min_{\boldsymbol{\beta}\in\mathbb{R}^{d_1}}\left\{M^{-1}\sum_{i\in\mathcal I_{\boldsymbol{\beta}}}\ell_4(\mathbf{W}_i;\widehat{\boldsymbol{\gamma}},\widehat{\boldsymbol{\delta}},\widehat{\boldsymbol{\alpha}},\boldsymbol{\beta})+\lambda_{\boldsymbol{\beta}}\|\boldsymbol{\beta}\|_1\right\}.\label{def:deltahat}
\end{align}
The corresponding loss functions are defined as:
\begin{align}
\ell_3(\mathbf{W};\boldsymbol{\gamma},\boldsymbol{\delta},\boldsymbol{\alpha})&~:=~\frac{A_1A_2\exp(-\bar\mathbf{S}_2^\top\boldsymbol{\delta})}{g(\mathbf{S}_1^\top\boldsymbol{\gamma})} \left(Y-\bar\mathbf{S}_2^\top\boldsymbol{\alpha} \right)^2,\label{def:l3}
\end{align}
\vspace{-2em}
\begin{align}
\ell_4(\mathbf{W};\boldsymbol{\gamma},\boldsymbol{\delta},\boldsymbol{\alpha},\boldsymbol{\beta})&~:=~A_1\exp(-\mathbf{S}_1^\top\boldsymbol{\gamma})\left\{\bar\mathbf{S}_2^\top\boldsymbol{\alpha}+\frac{A_2(Y-\bar\mathbf{S}_2^\top\boldsymbol{\alpha})}{g(\bar\mathbf{S}_2^\top\boldsymbol{\delta})}-\mathbf{S}_1^\top\boldsymbol{\beta}\right\}^2.\label{def:l4}
\end{align}


The above \eqref{def:l3}-\eqref{def:l4} 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{sec:exist_unique} 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 \eqref{rep:DR-mu}. Here, we are ensuring that the (imputed) residuals are uncorrelated with the adjustments made to the second and first PS estimations.



 The loss functions \eqref{def:alpha}, \eqref{def:l2}, \eqref{def:l3}, and \eqref{def:l4} are termed \emph{moment-targeting loss functions}. The \emph{sequential model doubly robust} (SMDR) estimator of \(\theta_{1,1}\) utilizing a cross-fitting technique can be found in Algorithm \ref{alg:BRDR}.

\begin{algorithm}[ht]
\caption{The sequential model doubly robust (SMDR) estimator of $\theta_{1,1}$}\label{alg:BRDR}
\begin{algorithmic}[1]
\Require Observations $\mathbb{S}=(\mathbf{W}_i)_{i=1}^N $ and the treatment path $(a_1,a_2)=(1,1)$.
\State Let $\mathcal I = \{1,2,\dots,N\}=\cup_{k=1}^{\mathbb{K}}\mathcal I_k$ with equal sized splits $n=N/\mathbb{K}$ and $\mathbb{K}\geq2$.
\For{$k=1,2,...,\mathbb{K}$}
\State $\mathcal I_{-k}\leftarrow\mathcal I\setminus\mathcal I_k$
\State $\mathcal I_{\boldsymbol{\gamma}},\mathcal I_{\boldsymbol{\delta}}, \mathcal I_{\boldsymbol{\alpha}}, \mathcal I_{\boldsymbol{\beta}}$ $\leftarrow$ size $M$ disjoint partition of $\mathcal I_{-k}$ with $M=N(\mathbb{K}-1)/(4\mathbb{K})$.
\State Propensity estimate at the first exposure $\widehat{\boldsymbol{\gamma}}_{-k}\leftarrow\widehat{\boldsymbol{\gamma}}$ as in \eqref{def:alphahat}.
\State Propensity estimate at the last/second exposure $\widehat{\boldsymbol{\delta}}_{-k} \leftarrow \widehat{\boldsymbol{\delta}}$ as in \eqref{def:betahat}.
\State Outcome estimate at the last exposure $\widehat{\boldsymbol{\alpha}}_{-k} \leftarrow\widehat{\boldsymbol{\alpha}}$ as in \eqref{def:gammahat}.
\State Outcome estimate at the first exposure $\widehat{\boldsymbol{\beta}}_{-k} \leftarrow \widehat{\boldsymbol{\beta}}$ as in \eqref{def:deltahat}.
\EndFor\\
\Return The SMDR estimator is
\begin{equation}
\widehat{\theta}_{1,1}=N^{-1}\sum_{k=1}^{\mathbb{K}}\sum_{i\in\mathcal I_k}\psi(\mathbf{W}_i;\widehat{\boldsymbol{\eta}}_{-k}),\;\;\mbox{where}\;\;\widehat{\boldsymbol{\eta}}_{-k} = (\widehat{\boldsymbol{\gamma}}_{-k}^\top,\widehat{\boldsymbol{\delta}}_{-k}^\top,\widehat{\boldsymbol{\alpha}}_{-k}^\top,\widehat{\boldsymbol{\beta}}_{-k}^\top)^\top\label{def:thetahat}
\end{equation}
and $\psi(\mathbf{W}_i;\widehat{\boldsymbol{\eta}}_{-k})$ is defined through \eqref{def:score} replacing $\boldsymbol{\eta}^*$ with $\widehat{\boldsymbol{\eta}}_{-k}$.
\end{algorithmic}
\end{algorithm}


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 \citep{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{cond:mis}.

In dynamic settings, \citep{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{table:settinga1}-\ref{table:lin}).


\begin{remark}[Correctness of nuisance models]\label{remark:correctness}
We now discuss the correct specification of nuisance models. The first outcome regression \(\mu\), is inherently challenging to interpret in dynamic settings \citep{babino2019multiple}. This complexity arises from the dynamic nature of the problem, not from our representation.
 Below, we introduce the specific meaning of a ``correctly specified model'':
\begin{enumerate}
\item[(a)] We say $\pi^*$ is correctly specified when $\pi^*=\pi$, which occurs if and only if (iff) there exists some $\boldsymbol{\gamma}^0\in\mathbb{R}^{d_1}$, such that $\pi(\mathbf{s}_1)=g(\mathbf{s}_1^\top\boldsymbol{\gamma}^0)$ holds. Additionally, $\boldsymbol{\gamma}^*=\boldsymbol{\gamma}^0$.
\item[(b)]
We say $\rho^*$ is correctly specified when $\rho^*=\rho$, which occurs iff there exists some $\boldsymbol{\delta}^0\in\mathbb{R}^d$, such that $\rho(\bar\mathbf{s}_2)=g(\bar\mathbf{s}_2^\top\boldsymbol{\delta}^0)$ holds. Additionally, $\boldsymbol{\delta}^*=\boldsymbol{\delta}^0$.
\item[(c)]
We say $\nu^*$ is correctly specified when $\nu^*=\nu$, which occurs iff there exists some $\boldsymbol{\alpha}^0\in\mathbb{R}^d$, such that $\nu(\bar\mathbf{s}_2)=\bar\mathbf{s}_2^\top\boldsymbol{\alpha}^0$ holds. Additionally, $\boldsymbol{\alpha}^*=\boldsymbol{\alpha}^0$.
\item[(d)]
We say $\mu^*$ is correctly specified when $\mu^*=\mu$, which occurs if there exists some $\boldsymbol{\beta}^0\in\mathbb{R}^d$, such that $\mu(\mathbf{s}_1)=\mathbf{s}_1^\top\boldsymbol{\beta}^0$ and, furthermore, either case (b) or (c) holds. Additionally, $\boldsymbol{\beta}^*=\boldsymbol{\beta}^0$.
\end{enumerate}


Note that $\boldsymbol{\delta}^*$ depends on $\boldsymbol{\gamma}^*$. However, condition (b) ensures that the correctness of $\rho^*$ is independent of $\boldsymbol{\gamma}^*$. Similarly, conditions (a)-(c) establish that the correctness of $\pi^*$, $\rho^*$, and $\nu^*$ are mutually independent.
 However, this is not the case for the OR model at the first exposure, $\mu^*$. Specifically, if $\mu$ is linear with $\mu(\mathbf{s}_1) = \mathbf{s}_1^\top\boldsymbol{\beta}^0$ for some $\boldsymbol{\beta}^0$, this does not imply that $\mu^*$ is correctly specified, as $\boldsymbol{\beta}^*$ may differ from $\boldsymbol{\beta}^0$. In particular, $\mu^*$ is correctly specified if, additionally, either $\rho^*$ or $\nu^*$ is (or both are) correctly specified. That is, under Assumption \ref{cond:mis}, $\mu^*$ is correctly specified if and only if $\mu$ is indeed a linear function -- no further constraints are needed. In contrast, the correctness of a linear $\mu^*$ defined through nested approaches \citep{murphy2001marginal, bodory2022evaluating} typically relies additionally on the linearity of $\nu$; see Section \ref{sec:just} of the Supplementary Materials.\end{remark}




\section{Sequential model doubly robust estimation and inference}\label{sec:DTE}

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.
\begin{assumption}[Sparsity]\label{cond:sparse}
Let $s_{\boldsymbol{\gamma}}+s_{\boldsymbol{\beta}}=o(N/\log d_1)$, $s_{\boldsymbol{\delta}}+s_{\boldsymbol{\alpha}}=o(N/\log d)$, and $(s_{\boldsymbol{\gamma}}+s_{\boldsymbol{\alpha}})\log d_1\log d+s_{\boldsymbol{\delta}}\log^2d=O(N)$.
\end{assumption}

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 \eqref{def:deltahat} 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., \cite{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$.

\begin{assumption}[Sub-Gaussianity]\label{cond:subG}
Let $\bar\mathbf{S}_2$ be a sub-Gaussian random vector with $\|\mathbf{v}^\top\bar\mathbf{S}_2\|_{\psi_2}\leq\sigma_{\mathbf{S}}\|\mathbf{v}\|_2$ for all $\mathbf{v}\in\mathbb{R}^d$. Let $\varepsilon:=Y(1,1)-\bar\mathbf{S}_2^\top\boldsymbol{\alpha}^*$ and $\zeta:=\bar\mathbf{S}_2^\top\boldsymbol{\alpha}^*-\mathbf{S}_1^\top\boldsymbol{\beta}^*$ be sub-Gaussian with $\|\varepsilon\|_{\psi_2}\leq\sigma_\varepsilon$ and $\|\zeta\|_{\psi_2}\leq\sigma_\zeta$.
In addition, let $\mathbb{E}[A_1A_2\{Y(1,1)-\nu(\bar\mathbf{S}_2)\}^2]>c_Y$ and the smallest eigenvalue of $\mbox{Cov}(A_1\bar\mathbf{S}_2)$ is bounded below by $c_{\min}$. Here, $\sigma_{\mathbf{S}},\sigma_\varepsilon,\sigma_\zeta,c_Y,c_{\min}$ are some positive constants.
\end{assumption}


\begin{theorem}[Convergence rates]\label{thm:rate}
Let Assumptions \ref{cond:basic}-\ref{cond:subG} hold. Define
\begin{equation}
r_{\boldsymbol{\gamma}}:=\sqrt\frac{s_{\boldsymbol{\gamma}}\log d_1}{N},\;\;r_{\boldsymbol{\delta}}:=\sqrt\frac{s_{\boldsymbol{\delta}}\log d}{N},\;\;r_{\boldsymbol{\alpha}}:=\sqrt\frac{s_{\boldsymbol{\alpha}}\log d}{N},\;\;r_{\boldsymbol{\beta}}:=\sqrt\frac{s_{\boldsymbol{\beta}}\log d_1}{N}.\label{def:rs}
\end{equation}
Then, as $N,d_1,d_2\to\infty$, $\sigma^2:=\mathbb{E}\{\psi(\mathbf{W};\boldsymbol{\eta}^*)-\theta_{1,1}\}^2\asymp\|\boldsymbol{\beta}^*\|_2+1$ and
\begin{align}
\widehat{\theta}_{1,1}-\theta_{1,1}&=O_p(\sigma N^{-1/2}+r_{\boldsymbol{\gamma}}r_{\boldsymbol{\beta}}+r_{\boldsymbol{\delta}}r_{\boldsymbol{\alpha}})+\mathbbm1_{\rho\neq\rho^*}O_p(r_{\boldsymbol{\gamma}}r_{\boldsymbol{\alpha}})\nonumber\\
&\qquad+\mathbbm1_{\nu\neq\nu^*}O_p(r_{\boldsymbol{\gamma}}r_{\boldsymbol{\delta}}+r_{\boldsymbol{\delta}}^2)+\mathbbm1_{\mu\neq\mu^*}O_p(r_{\boldsymbol{\gamma}}^2+r_{\boldsymbol{\gamma}}r_{\boldsymbol{\delta}}+r_{\boldsymbol{\gamma}}r_{\boldsymbol{\alpha}}).\label{rate:thetahat}
\end{align}
\end{theorem}

Theorem~\ref{thm:rate} 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 \citep{bradic2024high} and outperforms the DTL estimator \citep{bradic2024high, bodory2022evaluating}. Under model misspecification, DTL and S-DRL exhibit convergence rates with additional linear terms (see Table~\ref{table:rate}), whereas our result in \eqref{rate:thetahat} 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.


\begin{table}[h]
\caption{Convergence rates of DTL, S-DRL, and SMDR estimators under model misspecification situations in high dimensions. The sequences $r_{\boldsymbol{\gamma}}, r_{\boldsymbol{\delta}}, r_{\boldsymbol{\alpha}}, r_{\boldsymbol{\beta}}$ are defined in \eqref{def:rs}. For simplicity, we consider $\sigma \asymp 1$, and denote $r_0^2 := N^{-1/2} + r_{\boldsymbol{\gamma}}r_{\boldsymbol{\beta}} + r_{\boldsymbol{\delta}}r_{\boldsymbol{\alpha}}$. The quantities in \textcolor{red}{red} denote the additional linear terms in the convergence rates of DTL and S-DRL. The quantities in {\color{ao(english)} green} denote quadratic terms that decay faster than the corresponding {\color{red} red} terms on the same row. The quantities in {\color{orange} orange} denote the additional terms that appear in DTL's convergence rate only.} \label{table:rate}
\begin{center}
\begin{tabular}{| c | c | c | c | c | c | c |}
\hline
\multicolumn{4}{| c |}{Model correctness}&\multicolumn{3}{ c |}{Convergence raets}\\
\hline
$\pi^*$&$\rho^*$&$\nu^*$&$\mu^*$&DTL&S-DRL&SMDR\\
\hline
\ding{51}&\ding{51}&\ding{51}&\ding{51}&$r_0^2+{\color{orange} r_{\boldsymbol{\gamma}}r_{\boldsymbol{\alpha}}}$&$r_0^2$&$r_0^2$\\
\hline
\ding{51}&\ding{51}&\ding{51}&\ding{55}&$r_0^2+{\color{red} r_{\boldsymbol{\gamma}}}$&$r_0^2+{\color{red} r_{\boldsymbol{\gamma}}}$&$r_0^2+{\color{ao(english)} r_{\boldsymbol{\gamma}}(r_{\boldsymbol{\gamma}}+r_{\boldsymbol{\delta}}+r_{\boldsymbol{\alpha}})}$\\
\hline
\ding{51}&\ding{51}&\ding{55}&\ding{51}&$r_0^2+{\color{orange} r_{\boldsymbol{\gamma}}r_{\boldsymbol{\alpha}}}+{\color{red} r_{\boldsymbol{\delta}}}$&$r_0^2+{\color{red} r_{\boldsymbol{\delta}}}$&$r_0^2+{\color{ao(english)} r_{\boldsymbol{\delta}}(r_{\boldsymbol{\gamma}}+r_{\boldsymbol{\delta}})}$\\
\hline
\ding{51}&\ding{55}&\ding{51}&\ding{51}&$r_0^2+{\color{red} r_{\boldsymbol{\alpha}}}$&$r_0^2+{\color{red} r_{\boldsymbol{\alpha}}}$&$r_0^2+{\color{ao(english)} r_{\boldsymbol{\alpha}}r_{\boldsymbol{\gamma}}}$\\
\hline
\ding{55}&\ding{51}&\ding{51}&\ding{51}&$r_0^2+{\color{orange} r_{\boldsymbol{\alpha}}}+{\color{red} r_{\boldsymbol{\beta}}}$&$r_0^2+{\color{red} r_{\boldsymbol{\beta}}}$&$r_0^2$\\
\hline
\ding{51}&\ding{51}&\ding{55}&\ding{55}&$r_0^2+{\color{red} r_{\boldsymbol{\gamma}}+r_{\boldsymbol{\delta}}}$&$r_0^2+{\color{red} r_{\boldsymbol{\gamma}}+r_{\boldsymbol{\delta}}}$&$r_0^2+{\color{ao(english)} r_{\boldsymbol{\gamma}}(r_{\boldsymbol{\gamma}}+r_{\boldsymbol{\delta}}+r_{\boldsymbol{\alpha}})+r_{\boldsymbol{\delta}}^2}$\\
\hline
\ding{51}&\ding{55}&\ding{51}&\ding{55}&$r_0^2+{\color{red} r_{\boldsymbol{\alpha}}+r_{\boldsymbol{\gamma}}}$&$r_0^2+{\color{red} r_{\boldsymbol{\alpha}}+r_{\boldsymbol{\gamma}}}$&$r_0^2+{\color{ao(english)} r_{\boldsymbol{\gamma}}(r_{\boldsymbol{\gamma}}+r_{\boldsymbol{\delta}}+r_{\boldsymbol{\alpha}})}$\\
\hline
\ding{55}&\ding{51}&\ding{55}&\ding{51}&$r_0^2+{\color{orange} r_{\boldsymbol{\alpha}}}+{\color{red} r_{\boldsymbol{\beta}}+r_{\boldsymbol{\delta}}}$&$r_0^2+{\color{red} r_{\boldsymbol{\beta}}+r_{\boldsymbol{\delta}}}$&$r_0^2+{\color{ao(english)} r_{\boldsymbol{\delta}}(r_{\boldsymbol{\gamma}}+r_{\boldsymbol{\delta}})}$\\
\hline
\ding{55}&\ding{55}&\ding{51}&\ding{51}&$r_0^2+{\color{red} r_{\boldsymbol{\alpha}}+r_{\boldsymbol{\beta}}}$&$r_0^2+{\color{red} r_{\boldsymbol{\alpha}}+r_{\boldsymbol{\beta}}}$&$r_0^2+{\color{ao(english)} r_{\boldsymbol{\alpha}}r_{\boldsymbol{\gamma}}}$\\
\hline
\end{tabular}
\end{center}
\end{table}






\begin{theorem}[Inference under model misspecification]\label{thm:main}
Let Assumptions \ref{cond:basic}-\ref{cond:subG} hold. Let the following product sparsity conditions hold
\begin{equation}\label{cond:s2}
r_{\boldsymbol{\gamma}}r_{\boldsymbol{\beta}}=o(N^{-1/2})\quad\mbox{and}\quad r_{\boldsymbol{\delta}}r_{\boldsymbol{\alpha}}=o(N^{-1/2}),
\end{equation}
where sequences $r_{\boldsymbol{\gamma}}$, $r_{\boldsymbol{\delta}}$, $r_{\boldsymbol{\alpha}}$, and $r_{\boldsymbol{\beta}}$ are defined in \eqref{def:rs}.
We assume the following additional conditions if model misspecification occurs:
\begin{align}
&\text{if}\;\;\rho\neq\rho^*,\;\;\text{further let}\;\;r_{\boldsymbol{\gamma}}r_{\boldsymbol{\alpha}}=o(N^{-1/2});\label{cond:s3}\\
&\text{if}\;\;\nu\neq\nu^*,\;\;\text{further let}\;\;r_{\boldsymbol{\gamma}}r_{\boldsymbol{\delta}}=o(N^{-1/2}),\;\;r_{\boldsymbol{\delta}}=o(N^{-1/4});\label{cond:s4}\\
&\text{if}\;\;\mu\neq\mu^*,\;\;\text{further let}\;\;r_{\boldsymbol{\gamma}}=o(N^{-1/4}),\;\;r_{\boldsymbol{\gamma}}r_{\boldsymbol{\delta}}+r_{\boldsymbol{\gamma}}r_{\boldsymbol{\alpha}}=o(N^{-1/2}).\label{cond:s5}
\end{align}
Then, as $N,d_1,d_2\to\infty$,
$$\sigma^{-1}N^{1/2}(\widehat{\theta}_{1,1}-\theta_{1,1})\to\mathcal N(0,1)$$ in distribution and $\widehat{\sigma}^2=\sigma^2\{1+o_p(1)\}$, where $\widehat{\sigma}^2:=N^{-1}\sum_{k=1}^{\mathbb{K}}\sum_{i\in\mathcal I_k}\{\psi(\mathbf{W}_i;\widehat{\boldsymbol{\eta}}_{-k})-\widehat{\theta}_{1,1}\}^2$.
\end{theorem}


\begin{remark}[Sequential model double robustness]\label{remark:MDR}
In Theorem \ref{thm:main}, we establish the ``sequential model double robustness'' (SMDR) property of our proposed estimator, ensuring root-\(N\) inference as long as at least one nuisance model is correctly specified at each exposure (see Assumption \ref{cond:mis}), and some of the working models are estimated with \(o(N^{-1/4})\) rates. The specific nuisance estimation conditions (or equivalently, sparsity conditions) will be detailed in Remark \ref{remark:sparsity} and Table \ref{table:sparsity}.

The SMDR property differs from the sequential double robustness (SDR) established by \cite{luedtke2017sequential}, which guarantees only \emph{consistency} of the DTE estimates under model misspecification. Root-\(N\) inference under model misspecification has been achieved in low-dimensional settings \citep{rotnitzky2017multiply}, where nuisance estimates are assumed to satisfy Donsker conditions and exhibit parametric convergence rates, i.e., \(r_{\boldsymbol{\alpha}} \asymp r_{\boldsymbol{\beta}} \asymp r_{\boldsymbol{\gamma}} \asymp r_{\boldsymbol{\delta}} \asymp N^{-1/2}\). Due to the fast convergence rates they require, their method also fails to achieve the SMDR property that we aim to establish. In fact, provided these parametric convergence rates, the required conditions \eqref{cond:s2}-\eqref{cond:s5} above are automatically satisfied.

Our sequentially designed loss functions expand the complexity of the functional class to which the nuisance functions belong, requiring only some models to be estimated with a rate of \(o(N^{-1/4})\), while others can be estimated with even slower rates. For instance, when $\rho^*$ is misspecified, we allow for scenarios where \(r_{\boldsymbol{\alpha}} \asymp r_{\boldsymbol{\beta}} \asymp N^{-1/3}\) and \(r_{\boldsymbol{\gamma}} \asymp r_{\boldsymbol{\delta}} \asymp N^{-1.01/6}\). This is equivalent to \(s_{\boldsymbol{\alpha}} \log d \asymp s_{\boldsymbol{\beta}} \log d_1 \asymp N^{1/3}\) and \(s_{\boldsymbol{\gamma}} \log d \asymp s_{\boldsymbol{\delta}} \log d_1 \asymp N^{1.99/3}\), accommodating growing sparsity levels and dimensions that violate Donsker conditions. The proposed method extends root-\(N\) inference from low-dimensional to more challenging high-dimensional settings and achieves the SMDR property.

In high dimensions, the existing DTL and S-DRL methods \citep{bodory2022evaluating, bradic2024high} only ensure consistent DTE estimates without inferential guarantees, as long as at least one nuisance model is misspecified. Only \cite{viviano2021dynamic} provided valid inference without relying on specific parametric forms for the PS functions; however, they always require all OR models to be correctly specified, i.e., \eqref{CAN_a} holds. Our proposed strategy accommodates all cases \eqref{CAN_a}-\eqref{CAN_d}, achieving the best model robustness in high dimensions.
\end{remark}




\begin{remark}[Required sparsity conditions under model misspecification]\label{remark:sparsity}
We now discuss the sparsity conditions required for root-$N$ inference in Theorem \ref{thm:main}.

First, we examine a simpler static scenario and compare the sparsity conditions with existing literature on settings with a single exposure. These settings can be viewed as a special case of our framework, where the exposure $A_1$ is independent of $\bar{\mathbf{S}}_2$, $A_2$, and $Y$. In such cases, the DTE reduces to the average treatment effect (ATE), and robust inference under model misspecification has been established by \cite{smucler2019unifying}, \cite{tan2020model}, \cite{avagyan2021high}, \cite{ning2020robust}, among others. For these setups, the required sparsity conditions are $r_{\boldsymbol{\delta}} r_{\boldsymbol{\alpha}} = o(N^{-1/2})$, and if the outcome regression (OR) model $\nu^*$ is misspecified, we additionally require $r_{\boldsymbol{\delta}} = o(N^{-1/4})$, as given in \eqref{cond:s2} and \eqref{cond:s4}. These conditions match those in \cite{smucler2019unifying}, but are weaker than the conditions in \cite{tan2020model}, \cite{avagyan2021high}, and \cite{ning2020robust}, where a stronger condition $r_{\boldsymbol{\delta}} + r_{\boldsymbol{\alpha}} = o(N^{-1/4})$ is required.

Next, we consider the more complex dynamic scenarios. Similar to the static case, the correctness of one PS model, $\pi^*$, does not affect the required sparsity conditions. However, the correctness of the other PS model, $\rho^*$, along with both OR models, $\nu^*$ and $\mu^*$, influences the sparsity requirements. The more models that are misspecified, the more stringent the sparsity conditions become. When $\rho^*$, $\nu^*$, and $\mu^*$ are all correctly specified, we require Assumption \ref{cond:sparse} and \eqref{cond:s2}. If any model at exposure $t \in \{1, 2\}$ is misspecified, we impose a product condition between (i) the sparsity level of the other (correctly specified) model at the same exposure $t$ and (ii) the summation of sparsity levels corresponds to all the nuisance estimators that such a misspecified estimator is constructed based on. Recall that we estimate the nuisance models sequentially in the order: $\widehat{\boldsymbol{\gamma}}$, then $\widehat{\boldsymbol{\delta}}$, followed by $\widehat{\boldsymbol{\alpha}}$ and $\widehat{\boldsymbol{\beta}}$. For instance, when the OR model $\mu^*$ is misspecified, as shown in \eqref{cond:s5}, we require a product condition between (i) $s_{\boldsymbol{\gamma}}$ and (ii) $s_{\boldsymbol{\gamma}} + s_{\boldsymbol{\delta}} + s_{\boldsymbol{\alpha}}$. Moreover, if the OR model at the $t$-th exposure is misspecified, an ultra-sparse PS parameter is required at that exposure, as the OR models are estimated based on the PS estimates.

Our results provide a clear framework for understanding the additional challenges posed by model misspecification. Since achieving the SMDR property requires sequential estimation of the PS models followed by backward estimation of the OR models, misspecification of early-stage OR models imposes the most stringent conditions on the model structures. In contrast, misspecification of the first-stage PS model does not impact the sparsity conditions. The sparsity conditions required in \cite{smucler2019unifying} can be viewed as a special case of the more general phenomenon we identify for single-exposure settings. Further details are provided in Table \ref{table:sparsity}.
\end{remark}




Whenever all nuisance models are correctly specified, we have the following result.
\begin{theorem}[Inference under correctly specified models]\label{cor:correct}
Suppose all the nuisance models are correctly specified. Let Assumptions \ref{cond:basic}, \ref{cond:sparse}, and \ref{cond:subG} hold, as well as the product sparsity \eqref{cond:s2}. Then, as $N,d_1,d_2\to\infty$,
$$\sigma^{-1}N^{1/2}(\widehat{\theta}_{1,1}-\theta_{1,1})~\to~\mathcal N(0,1)$$ in distribution and $\widehat{\sigma}^2=\sigma^2\{1+o_p(1)\}$.
\end{theorem}

When all nuisance functions are correctly specified, the result coincides with that of \cite{bradic2024high} while also achieving the semi-parametric efficiency of \cite{bang2005doubly}. Hence, we do not lose accuracy when the nuisance models are correctly specified.
As shown in Theorem \ref{cor:correct}, root-$N$ inference requires product sparsity conditions between the nuisance parameters' sparsity levels at each exposure, i.e., \eqref{cond:s2}; 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.





\begin{table}[h]
\caption{Sparsity conditions required for the SMDR estimator asymptotically normal. For simplicity, $\|\bar\mathbf{S}_2\|_\infty<C$, $d_1\asymp d$, and $s_{\boldsymbol{\gamma}}+s_{\boldsymbol{\delta}}+s_{\boldsymbol{\alpha}}+s_{\boldsymbol{\beta}}=o\left(N/\log d\right)$. } \label{table:sparsity}
\begin{center}
\begin{tabular}{| c | c | c | c | c |}
\hline
\multicolumn{4}{| c |}{Model correctness}&\multirow{2}{*}{Required sparsity conditions}\\
\cline{1-4}
$\pi^*$&$\rho^*$&$\nu^*$&$\mu^*$&\\
\hline
\ding{51}&\ding{51}&\ding{51}&\ding{51}&$s_{\boldsymbol{\gamma}}s_{\boldsymbol{\beta}}+s_{\boldsymbol{\delta}}s_{\boldsymbol{\alpha}}=o\left(\frac{N}{\log^2d}\right)$\\
\hline
\ding{51}&\ding{51}&\ding{51}&\ding{55}&$s_{\boldsymbol{\gamma}}=o\left(\frac{\sqrt N}{\log d}\right),\;s_{\boldsymbol{\gamma}}s_{\boldsymbol{\delta}}+s_{\boldsymbol{\gamma}}s_{\boldsymbol{\alpha}}+s_{\boldsymbol{\gamma}}s_{\boldsymbol{\beta}}+s_{\boldsymbol{\delta}}s_{\boldsymbol{\alpha}}=o\left(\frac{N}{\log^2d}\right)$\\
\hline
\ding{51}&\ding{51}&\ding{55}&\ding{51}&$s_{\boldsymbol{\delta}}=o\left(\frac{\sqrt N}{\log d}\right),\;s_{\boldsymbol{\gamma}}s_{\boldsymbol{\delta}}+s_{\boldsymbol{\gamma}}s_{\boldsymbol{\beta}}+s_{\boldsymbol{\delta}}s_{\boldsymbol{\alpha}}=o\left(\frac{N}{\log^2d}\right)$\\
\hline
\ding{51}&\ding{55}&\ding{51}&\ding{51}&$s_{\boldsymbol{\gamma}}s_{\boldsymbol{\alpha}}+s_{\boldsymbol{\gamma}}s_{\boldsymbol{\beta}}+s_{\boldsymbol{\delta}}s_{\boldsymbol{\alpha}}=o\left(\frac{N}{\log^2d}\right)$\\
\hline
\ding{55}&\ding{51}&\ding{51}&\ding{51}&$s_{\boldsymbol{\gamma}}s_{\boldsymbol{\beta}}+s_{\boldsymbol{\delta}}s_{\boldsymbol{\alpha}}=o\left(\frac{N}{\log^2d}\right)$\\
\hline
\ding{51}&\ding{51}&\ding{55}&\ding{55}&$s_{\boldsymbol{\gamma}}+s_{\boldsymbol{\delta}}=o\left(\frac{\sqrt N}{\log d}\right),\;s_{\boldsymbol{\gamma}}s_{\boldsymbol{\alpha}}+s_{\boldsymbol{\gamma}}s_{\boldsymbol{\beta}}+s_{\boldsymbol{\delta}}s_{\boldsymbol{\alpha}}=o\left(\frac{N}{\log^2d}\right)$\\
\hline
\ding{51}&\ding{55}&\ding{51}&\ding{55}&$s_{\boldsymbol{\gamma}}=o\left(\frac{\sqrt N}{\log d}\right),\;s_{\boldsymbol{\gamma}}s_{\boldsymbol{\delta}}+s_{\boldsymbol{\gamma}}s_{\boldsymbol{\alpha}}+s_{\boldsymbol{\gamma}}s_{\boldsymbol{\beta}}+s_{\boldsymbol{\delta}}s_{\boldsymbol{\alpha}}=o\left(\frac{N}{\log^2d}\right)$\\
\hline
\ding{55}&\ding{51}&\ding{55}&\ding{51}&$s_{\boldsymbol{\delta}}=o\left(\frac{\sqrt N}{\log d}\right),\;s_{\boldsymbol{\gamma}}s_{\boldsymbol{\delta}}+s_{\boldsymbol{\gamma}}s_{\boldsymbol{\beta}}+s_{\boldsymbol{\delta}}s_{\boldsymbol{\alpha}}=o\left(\frac{N}{\log^2d}\right)$\\
\hline
\ding{55}&\ding{55}&\ding{51}&\ding{51}&$s_{\boldsymbol{\gamma}}s_{\boldsymbol{\alpha}}+s_{\boldsymbol{\gamma}}s_{\boldsymbol{\beta}}+s_{\boldsymbol{\delta}}s_{\boldsymbol{\alpha}}=o\left(\frac{N}{\log^2d}\right)$\\
\hline
\end{tabular}
\end{center}
\end{table}




\section{Theoretical results for the nuisance estimators}\label{sec:nuis}

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 \eqref{def:alphahat}-\eqref{def:deltahat}. The analysis of nuisance estimation is non-trivial since the nuisance estimates are constructed in a sequential manner. Section \ref{sec:asymp_nuisance} shows the nuisance estimators' consistency despite potential model misspecification, while Section \ref{sec:asymp_nuisance'} presents their faster consistency rates when certain models are correctly specified. Our findings show that the accuracy of nuisance models affects the estimation errors.


\subsection{Results for misspecified models}\label{sec:asymp_nuisance}

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{lemma:RSC}, and manage the increased gradient variability in Lemma \ref{lemma:gradient}, as expanded upon in the Supplementary Material.

\begin{theorem}\label{thm:nuisance}
Let Assumptions \ref{cond:basic} and \ref{cond:subG} hold. Define sequences $r_{\boldsymbol{\gamma}}$, $r_{\boldsymbol{\delta}}$, $r_{\boldsymbol{\alpha}}$, and $r_{\boldsymbol{\beta}}$ as in \eqref{def:rs}. Then, as $N,d_1,d_2\to\infty$, the following holds:
\begin{enumerate}
\item[(a)] If $r_{\boldsymbol{\gamma}}=o(1)$, then
$\|\widehat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^*\|_2=O_p(r_{\boldsymbol{\gamma}}).$
\item[(b)] In addition to (a), if $r_{\boldsymbol{\delta}}=o(1)$, then
$\|\widehat{\boldsymbol{\delta}}-\boldsymbol{\delta}^*\|_2= O_p(r_{\boldsymbol{\gamma}}+r_{\boldsymbol{\delta}}).$
\item[(c)] In addition to (a) and (b), if $r_{\boldsymbol{\alpha}}=o(1)$, then
$\|\widehat{\boldsymbol{\alpha}}-\boldsymbol{\alpha}^*\|_2= O_p(r_{\boldsymbol{\gamma}}+r_{\boldsymbol{\delta}}+r_{\boldsymbol{\alpha}}).$
\item[(d)] In addition to (a), (b), and (c), if $r_{\boldsymbol{\beta}}=o(1)$, then $
\|\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^*\|_2= O_p(r_{\boldsymbol{\gamma}}+r_{\boldsymbol{\delta}}+r_{\boldsymbol{\alpha}}+r_{\boldsymbol{\beta}}).$
\end{enumerate}
\end{theorem}

 Among the results in Theorem \ref{thm:nuisance}, 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 \cite{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 \eqref{def:l2} 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 \cite{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{lemma:beta2}. Together with some empirical process results as in Lemma \ref{lemma:beta1}, we control the imputation error's effect and finally reach the consistency rates introduced above; see Lemma \ref{lemma:beta3} and the proof of Theorem \ref{thm:nuisance}. Although we focus on a specific loss function \eqref{def:l2}, 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.


\subsection{Results for correctly specified models}\label{sec:asymp_nuisance'}

If we have additional information that some of the nuisance models are correctly specified, we are able to achieve better consistency rates.

\begin{theorem}\label{thm:nuisance'}
Let Assumptions \ref{cond:basic} and \ref{cond:subG} hold. Suppose that the sequences defined in \eqref{def:rs} satisfy $r_{\boldsymbol{\gamma}}+r_{\boldsymbol{\delta}}+r_{\boldsymbol{\alpha}}+r_{\boldsymbol{\beta}}=o(1)$. Then, as $N,d_1,d_2\to\infty$, the following holds:
\begin{enumerate}
\item[(a)] Let $\rho=\rho^*$ and $s_{\boldsymbol{\gamma}}=O(N/(\log d_1\log d))$, then
 $\|\widehat{\boldsymbol{\delta}}-\boldsymbol{\delta}^*\|_2=O_p(r_{\boldsymbol{\delta}}).$
\item[(b)] Let $\nu=\nu^*$, $s_{\boldsymbol{\gamma}}$ as in (a) and $s_{\boldsymbol{\delta}}=O(N/\log^2d)$, then
 $\|\widehat{\boldsymbol{\alpha}}-\boldsymbol{\alpha}^*\|_2=O_p(r_{\boldsymbol{\alpha}}).$
\item[(c)] Let $\nu=\nu^*$, $\mu=\mu^*$, $s_{\boldsymbol{\gamma}}$ and $s_{\boldsymbol{\delta}}$ are as in (b), then
 $\|\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^*\|_2 =O_p(r_{\boldsymbol{\alpha}}+r_{\boldsymbol{\beta}}).$
\item[(d)] Let $\rho=\rho^*$, $\mu=\mu^*$, and $s_{\boldsymbol{\gamma}}+s_{\boldsymbol{\delta}}+s_{\boldsymbol{\alpha}}=O(N/(\log d_1\log d))$, then
 $\|\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^*\|_2 =O_p(r_{\boldsymbol{\delta}}+r_{\boldsymbol{\beta}}).$
 \item[(e)] Let $\rho=\rho^*$, $\nu=\nu^*$, $\mu=\mu^*$, $s_{\boldsymbol{\gamma}}$ and $s_{\boldsymbol{\delta}}$ are as in (b) and $s_{\boldsymbol{\alpha}}=O(N/(\log d_1\log d))$, then
 $\|\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^*\|_2 =O_p(r_{\boldsymbol{\delta}}r_{\boldsymbol{\alpha}}+r_{\boldsymbol{\beta}}).$
\end{enumerate}
\end{theorem}

The new convergence rates in Theorem \ref{thm:nuisance'} are established through Lemmas \ref{lemma:RSC} and \ref{lemma:gradient'} of the Supplementary Material. Assuming certain nuisance models being correct, unlike Theorem \ref{thm:nuisance} and Lemma \ref{lemma:gradient}, we can control the gradients involving the \emph{estimated} nuisance parameters and control the imputation errors from the previous steps' nuisances in a more efficient way; see more details in Lemma \ref{lemma:gradient'}. As a result, we obtain faster convergence rates than Theorem \ref{thm:nuisance} 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 \eqref{cond:s2} of Theorem \ref{thm:main}, $\widehat{\boldsymbol{\beta}}$ also becomes asymptotically decoupled from the other three estimates.


\section{Numerical Experiments}\label{sec:num}
\subsection{Simulation studies}\label{sec:sim}
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{alg:BRDR} 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{alg:BRDR}.

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 \citep{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{table:settinga1}-\ref{table:settingb1}.

\begin{table}[h!]
\centering
\caption{Simulation under Setting (a) with $d_1=100$, $d=150$. Bias: empirical bias; RMSE: root mean square error; Length: average length of the $95\%$ confidence intervals; Coverage: average coverage of the $95\%$ confidence intervals; ESD: empirical standard deviation; ASD: average of estimated standard deviations. All the reported values (except Coverage) are based on robust (median-type) estimates. $N_1$ and $N_0$ denote the expected numbers of observations in the treatment groups $(1,1)$ and $(0,0)$, respectively.} \label{table:settinga1}

\begin{tabular}{lcccccccccc}
\toprule
Method&Bias&RMSE&Length&Coverage&&Bias&RMSE&Length&Coverage\\
\hline
\multicolumn{1}{c}{ } & \multicolumn{4}{c}{\cellcolor{gray!50} $N=400,N_1\approx136,N_0\approx68$}&&\multicolumn{4}{c}{ \cellcolor{gray!50} $N=1000,N_1\approx341,N_0\approx169$}\\
\cline{2-5}\cline{7-11}
oracle&0.037&0.298&1.627&0.955&&0.004&0.190&1.056&0.975\\
\cdashline{2-5}\cdashline{7-11}
empdiff&-0.409&0.460&0.491&0.325&&-0.392&0.392&0.310&0.120\\
\cdashline{2-5}\cdashline{7-11}
IPW&0.553&0.582&1.941&0.770&&0.730&0.730&1.254&0.390\\
\cdashline{2-5}\cdashline{7-11}
SDR&0.597&4.067&14.304&0.800&&0.126&0.479&2.396&0.930\\
\cdashline{2-5}\cdashline{7-11}
SDR-RF&0.294&0.353&1.564&0.885&&0.423&0.423&0.945&0.595\\
\cdashline{2-5}\cdashline{7-11}
DTL1&0.643&0.775&2.086&0.775&&0.488&0.500&1.048&0.580\\
\cdashline{2-5}\cdashline{7-11}
DTL2&0.466&0.499&1.597&0.775&&0.278&0.286&1.044&0.775\\
\cdashline{2-5}\cdashline{7-11}
S-DRL&0.507&0.520&1.632&0.710&&0.311&0.313&1.055&0.785\\
\cdashline{2-5}\cdashline{7-11}
SMDR1&0.664&0.666&1.633&0.600&&0.462&0.462&0.975&0.530\\
\cdashline{2-5}\cdashline{7-11}
SMDR2&0.307&0.359&1.546&0.840&&0.091&0.193&0.998&0.890\\
\hline
\multicolumn{1}{c}{ } & \multicolumn{4}{c}{\cellcolor{gray!50} $N=12000,N_1\approx4103,N_0\approx2033$}&&\multicolumn{4}{c}{ \cellcolor{gray!50} $N=16000,N_1\approx5471,N_0\approx2710$}\\
\cline{2-5}\cline{7-11}
oracle&-0.002&0.053&0.317&0.945&&0.007&0.056&0.276&0.955\\
\cdashline{2-5}\cdashline{7-11}
empdiff&-0.401&0.401&0.090&0.000&&-0.397&0.397&0.078&0.000\\
\cdashline{2-5}\cdashline{7-11}
IPW&0.946&0.946&0.418&0.000&&0.957&0.957&0.369&0.000 \\
\cdashline{2-5}\cdashline{7-11}
SDR&0.012&0.090&0.571&0.955&&0.011&0.073&0.508&0.945\\
\cdashline{2-5}\cdashline{7-11}
SDR-RF&0.560&0.560&0.276&0.000&&0.567&0.567&0.241&0.000\\
\cdashline{2-5}\cdashline{7-11}
DTL1&0.243&0.243&0.329&0.260&&0.212&0.212&0.293&0.235\\
\cdashline{2-5}\cdashline{7-11}
DTL2&0.137&0.141&0.355&0.670&&0.122&0.122&0.314&0.655\\
\cdashline{2-5}\cdashline{7-11}
S-DRL&0.143&0.143&0.356&0.650&&0.123&0.123&0.313&0.670 \\
\cdashline{2-5}\cdashline{7-11}
SMDR1&0.053&0.075&0.311&0.890&&0.048&0.069&0.269&0.920\\
\cdashline{2-5}\cdashline{7-11}
SMDR2&0.020&0.058&0.319&0.935&&0.013&0.053&0.277&0.925\\
\bottomrule
\end{tabular}
\end{table}

\begin{table}[h]
\centering
\caption{Simulation under Setting (b) with $d_1=100$, $d=150$. The rest of the caption details remain the same as those in Table \ref{table:settinga1}.} \label{table:settingb1}
\begin{tabular}{lcccccccccc}
\toprule
Method&Bias&RMSE&Length&Coverage&&Bias&RMSE&Length&Coverage\\
\hline
\multicolumn{1}{c}{ } & \multicolumn{4}{c}{\cellcolor{gray!50} $N=400,N_1\approx109,N_0\approx92$}&&\multicolumn{4}{c}{ \cellcolor{gray!50} $N=1000,N_1\approx271,N_0\approx227$}\\
\cline{2-5}\cline{7-11}
oracle&-0.016&0.140&1.038&0.980&&-0.007&0.110&0.663&0.940\\
\cdashline{2-5}\cdashline{7-11}
empdiff&-0.066&0.220&0.369&0.450&&-0.097&0.173&0.239&0.330\\
\cdashline{2-5}\cdashline{7-11}
IPW&-0.170&0.266&1.581&0.950&&-0.044&0.153&0.941&0.980\\
\cdashline{2-5}\cdashline{7-11}
SDR&-0.259&3.874&14.351&0.740&&-0.156&0.386&1.303&0.875\\
\cdashline{2-5}\cdashline{7-11}
SDR-RF&-0.109&0.219&1.382&0.940&&-0.121&0.169&0.844&0.925\\
\cdashline{2-5}\cdashline{7-11}
DTL1&-0.078&0.366&1.908&0.940&&-0.161&0.220&0.972&0.895\\
\cdashline{2-5}\cdashline{7-11}
DTL2&-0.161&0.252&1.410&0.905&&-0.160&0.194&0.825&0.865\\
\cdashline{2-5}\cdashline{7-11}
S-DRL&-0.160&0.265&1.418&0.915&&-0.153&0.187&0.815&0.895\\
\cdashline{2-5}\cdashline{7-11}
SMDR1&-0.083&0.223&1.642&0.925&&-0.154&0.193&0.954&0.880\\
\cdashline{2-5}\cdashline{7-11}
SMDR2&-0.156&0.232&1.393&0.910&&-0.114&0.146&0.816&0.925\\
\hline
\multicolumn{1}{c}{ } & \multicolumn{4}{c}{\cellcolor{gray!50} $N=12000,N_1\approx3264,N_0\approx2726$}&&\multicolumn{4}{c}{ \cellcolor{gray!50} $N=16000,N_1\approx4352,N_0\approx3634$}\\
\cline{2-5}\cline{7-11}
oracle&-0.005&0.034&0.192&0.950&&0.002&0.030&0.166&0.960\\
\cdashline{2-5}\cdashline{7-11}
empdiff&-0.103&0.103&0.069&0.135&&-0.094&0.094&0.060&0.125\\
\cdashline{2-5}\cdashline{7-11}
IPW&0.034&0.048&0.274&0.945&&0.033&0.045&0.237&0.965\\
\cdashline{2-5}\cdashline{7-11}
SDR&-0.010&0.050&0.266&0.940&&-0.009&0.042&0.231&0.965\\
\cdashline{2-5}\cdashline{7-11}
SDR-RF&-0.090&0.090&0.219&0.600&&-0.087&0.087&0.189&0.565\\
\cdashline{2-5}\cdashline{7-11}
DTL1&-0.123&0.123&0.237&0.475&&-0.100&0.100&0.210&0.540\\
\cdashline{2-5}\cdashline{7-11}
DTL2&-0.072&0.073&0.246&0.775&&-0.060&0.062&0.217&0.880\\
\cdashline{2-5}\cdashline{7-11}
S-DRL&-0.060&0.069&0.246&0.790&&-0.048&0.051&0.215&0.815\\
\cdashline{2-5}\cdashline{7-11}
SMDR1&-0.035&0.051&0.227&0.890&&-0.018&0.037&0.198&0.950\\
\cdashline{2-5}\cdashline{7-11}
SMDR2&-0.013&0.047&0.228&0.940&&0.000&0.034&0.199&0.935\\
\bottomrule
\end{tabular}
\end{table}

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{table:settinga1}-\ref{table:settingb1}. 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{table:settingb1}.

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 \citep{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{table:settinga1}) or smaller than (see Table \ref{table:settingb1}) 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{table:settinga1}-\ref{table:settingb1}). 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{table:settinga1} 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{alg:BRDR} 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{table:settinga1}. 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.


\subsection{A semi-synthetic analysis based on the National Job Corps Study (NJCS)}

\begin{table}[h!]
\centering
\caption{Semi-synthetic analysis. Bias: empirical bias; CI: the $95\%$ confidence interval; p-value: the p-value of $H_0: \theta=0$ v.s. $H_1: \theta\neq0$.} \label{table:lin}
\resizebox{15cm}{!}{
\begin{tabular}{lcccccccccccc}
\toprule
Method&$\widehat{\theta}_O$&$\widehat{\theta}$&Bias&CI&p-value&&$\widehat{\theta}_O$&$\widehat{\theta}$&Bias&CI&p-value&\\
\hline
\multicolumn{1}{c}{ } & \multicolumn{12}{c}{\cellcolor{gray!50} Setting (a)}\\
\cline{2-13}
DTL2&\multirow{3}{*}{0.15}&0.103&-0.047&[-0.090, 0.297]&0.295&&\multirow{3}{*}{0.20}&0.153&-0.047&[-0.040, 0.347]&0.120\\
\cdashline{3-6}\cdashline{9-12}
S-DRL&&0.092&-0.058&[-0.101, 0.308]&0.378&&&0.142&-0.058&[-0.051, 0.358]&0.173\\
\cdashline{3-6}\cdashline{9-12}
SMDR1&&0.180&0.030&[-0.019, 0.380]&0.076&&&0.230&0.030&[0.031, 0.430]&0.024\\
\hline
DTL2&\multirow{3}{*}{0.25}&0.203&-0.047&[0.010, 0.397]&0.039&&\multirow{3}{*}{0.30}&0.253&-0.047&[0.060, 0.447]&0.010\\
\cdashline{3-6}\cdashline{9-12}
S-DRL&&0.192&-0.058&[-0.001, 0.408]&0.066&&&0.241&-0.059&[0.049, 0.458]&0.021\\
\cdashline{3-6}\cdashline{9-12}
SMDR1&&0.280&0.030&[0.080, 0.480]&0.005&&&0.330&0.030&[0.131, 0.530]&0.001\\
\hline
\multicolumn{1}{c}{ } & \multicolumn{12}{c}{\cellcolor{gray!50} Setting (b)}\\
\cline{2-13}
DTL2&\multirow{3}{*}{0.25}&-0.027&-0.277&[-0.223,0.168]&0.783&&\multirow{3}{*}{0.30}&0.023&-0.277&[-0.172,0.218]&0.817\\
\cdashline{3-6}\cdashline{9-12}
S-DRL&&0.013&-0.237&[-0.232, 0.178]&0.902&&&0.063&-0.237&[-0.182,0.228]&0.548\\
\cdashline{3-6}\cdashline{9-12}
SMDR1&&0.141&-0.109&[-0.038, 0.320]&0.123&&&0.192&-0.108&[0.013, 0.371]&0.035\\
\hline
DTL2&\multirow{3}{*}{0.35}&0.072&-0.278&[-0.123,0.268]&0.468&&\multirow{3}{*}{0.40}&0.122&-0.278&[-0.074,0.317]&0.222\\
\cdashline{3-6}\cdashline{9-12}
S-DRL&&0.113&-0.237&[-0.133,0.277]&0.281&&&0.162&-0.238&[-0.083,0.327]&0.121\\
\cdashline{3-6}\cdashline{9-12}
SMDR1&&0.242&-0.108&[0.063, 0.421]&0.008&&&0.291&-0.109&[0.111, 0.470]&0.001\\
\hline
DTL2&\multirow{3}{*}{0.45}&0.172&-0.278&[-0.023,0.367]&0.084&&\multirow{3}{*}{0.50}&0.223&-0.277&[0.028,0.418]&0.025\\
\cdashline{3-6}\cdashline{9-12}
S-DRL&&0.212&-0.238&[-0.033,0.377]&0.043&&&0.263&-0.237&[0.018,0.428]&0.012\\
\cdashline{3-6}\cdashline{9-12}
SMDR1&&0.340&-0.110&[0.161, 0.519]&0.000&&&0.406&-0.094&[0.225, 0.586]&0.000\\
\bottomrule
\end{tabular}}
\end{table}

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 \cite{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 \citep{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{alg:BRDR}. 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.

\begin{figure*}
	\captionsetup[subfloat]{labelformat=empty}
	\subfloat[Setting (a)]{
	\includegraphics[height=0.35\linewidth,width=0.45\linewidth]{Rplot_1.pdf}
	}
 \hfill
 \subfloat[Setting (b)]{
 \includegraphics[height=0.35\linewidth,width=0.45\linewidth]{Rplot_2.pdf}
 }
\caption{\centering The p-values for the null $H_0: \theta = 0$ as $\hat{\theta}_O$ varies. Since the true $\theta$ is unknown, the x-axis denotes the oracle difference-in-mean estimate of $\theta$.}\label{fig_JC}
\end{figure*}


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{sec:sim}). The results are reported in Table \ref{table:lin}, 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{fig_JC}. 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.


















\section{Discussion}\label{sec:dis}

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 \citep{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{thm:nuisance}, 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.



\section*{Supplementary Material}\label{supp_mat}

Sections \ref{sec:notation}-\ref{sec:proof_lemmas} contain additional discussions, justifications, and proofs of the main results. Additional notations used in the supplementary material are introduced in Section \ref{sec:notation}. Section \ref{sec:exist_unique} discusses the uniqueness of the moment-targeted parameters; the justification of their identification is provided in Section \ref{sec:just}. We introduce some useful auxiliary lemmas in Section \ref{sec:lemmas}. The proofs of main results and auxiliary lemmas are in Sections \ref{sec:proof_DTE} and \ref{sec:proof_lemmas}, respectively.