# A Meta-Learning Method for Estimation of Causal Excursion Effects to Assess
Time-Varying Moderation

Jieru Shi, Walter Dempsey
Department of Biostatistics, University of Michigan

###### Abstract

Twin revolutions in wearable technologies and smartphone-delivered digital
health interventions have significantly expanded the accessibility and uptake
of mobile health (mHealth) interventions across various health science
domains. Sequentially randomized experiments called micro-randomized trials
(MRTs) have grown in popularity to empirically evaluate the effectiveness of
these mHealth intervention components. MRTs have given rise to a new class of
causal estimands known as “causal excursion effects”, which enable health
scientists to assess how intervention effectiveness changes over time or is
moderated by individual characteristics, context, or responses in the past.
However, current data analysis methods for estimating causal excursion effects
require pre-specified features of the observed high-dimensional history to
construct a working model of an important nuisance parameter. While machine
learning algorithms are ideal for automatic feature construction, their naive
application to causal excursion estimation can lead to bias under model
misspecification, potentially yielding incorrect conclusions about
intervention effectiveness. To address this issue, this paper revisits the
estimation of causal excursion effects from a meta-learner perspective, where
the analyst remains agnostic to the choices of supervised learning algorithms
used to estimate nuisance parameters. The paper presents asymptotic properties
of the novel estimators and compares them theoretically and through extensive
simulation experiments, demonstrating relative efficiency gains and supporting
the recommendation for a doubly robust alternative to existing methods.
Finally, the practical utility of the proposed methods is demonstrated by
analyzing data from a multi-institution cohort of first-year medical residents
in the United States (NeCamp et al., 2020).

Keywords: Debiased/Orthogonal Estimation, Machine Learning, Double Robustness,
Causal Excursion Effect, Mobile Health, Time-Varying Treatment.

##  1 Introduction

The use of smart devices (e.g., smartphones, smartwatches) and other wearables
to deliver digital interventions to improve health outcomes has grown
significantly in the past few years. Low-cost, accessible digital
interventions can be delivered everywhere, anytime, and in any amount, even to
reticent or hard-to-reach populations. Interventions of this type are
hypothesized to result in meaningful short- and long-term behavior changes.
The assessment of such time-varying effects prompted the development of micro-
randomized trials (MRTs), in which individuals are randomized to receive
notifications at hundreds or thousands of decision points.

The MRT enables estimation of proximal or lagged effects of push notifications
on pre-specified outcomes of interest, referred to as “causal excursion
effects” (Boruvka et al., 2018; Qian et al., 2020; Dempsey et al., 2020; Shi
et al., 2022). Semiparametric inference of the causal excursion effects can be
conducted via a weighted, centered least squares (WCLS) criterion (Boruvka et
al., 2018). A key feature of implementing the WCLS criterion is that health
scientists must pre-specify features from the high-dimensional observed
history to formulate a linear working model for a critical nuisance parameter,
which is a challenging, non-trivial task.

Machine learning (ML) algorithms offer powerful tools to automatically
construct features for the nuisance components, but their naive application to
semiparametric inference can lead to bias in estimation. Chernozhukov et al.
(2018) shows that Neyman-orthogonal moments and cross-fitting can remove the
impact of regularization bias and overfitting caused by naive application of
the ML methods. Later, several meta-learner algorithms emerged that can take
advantage of supervised learning or regression methods in ML and statistics
—such as random forests, Bayesian Additive Regression Trees (BART) or neural
networks— to estimate the nuisance components (Hill, 2011; Semenova and
Chernozhukov, 2021; Künzel et al., 2019; Nie and Wager, 2021; Kennedy, 2020).
These papers provide flexible, well-performing methods for estimating
conditional average treatment effects (CATE) in randomized controlled trials
and observational studies, and illustrate how ML methods can be applied for
semiparametric inference.

While meta-learning approaches have been developed extensively for CATE
estimation, their application to longitudinal studies has been relatively
limited. Viviano and Bradic (2021) proposed a dynamic covariate balancing
method when high-dimensional covariates are present. Bodory et al. (2022) used
this approach to examine effects under dynamic treatment regimes using Double
Machine Learning (DML) (Chernozhukov et al., 2018) and semiparametrically
efficient estimation. In this setting, DML is used to control for observed and
time-varying covariates in a data-driven way, across treatment sequences in
different time periods. Lewis and Syrgkanis (2020) proposed a new DML approach
for estimating causal effect under dynamic treatment regimes by using
g-estimation – a sequential residualization approach that uses supervised
learning of debiased outcomes on debiased treatments over a specific time
period based on linear parameterization of blip functions.

While prior studies on DML methods for longitudinal analysis have mainly
focused on estimating the average treatment effect (ATE) on a distal outcome,
often under pre-specified or dynamic non-random treatment sequences, our
current work takes a different perspective. Specifically, we propose a meta-
learning framework that assesses time-varying causal effect moderation in
MRTs, providing a novel solution for estimating moderated treatment effects on
proximal outcomes using DML methods. The proposed method can help health
scientists improve their ability to answer critical scientific questions
regarding time-varying effect moderation, and find out when, in what context,
and what intervention content to deliver to each person to make the
intervention most effective (Qian et al., 2022).

###  1.1 Outline

The rest of the paper proceeds as follows. Section 2 reviews existing analytic
techniques for MRTs and ML methods in causal inference. We then summarize our
main contributions in Section 3 and explain why allowing for the use of ML
approaches in WCLS is challenging. We propose two new inferential methods in
Section 4, which leverage supervised learning methods to improve the
efficiency of time-varying treatment effect estimation. In Section 5.1, we
outline the inferential algorithms. Section 5.2 and 5.3 present the asymptotic
theory as well as discuss the relative efficiency gain of the proposed
methods. Section 6 discusses the extension of the proposed methods to settings
such as missing data, lagged effects, and binary outcomes. Section 7 uses
simulation studies to compare various estimators and standard errors. Section
8 illustrates the efficiency improvement using our proposed methods with a
recent MRT: the Intern Health Study (NeCamp et al., 2020). The paper concludes
with a brief discussion in Section 9. All technical proofs are collected in
the Supplementary Materials.

##  2 Preliminaries

###  2.1 Micro-Randomized Trials (MRT)

An MRT consists of a sequence of within-subject decision times $t=1,\ldots,T$
at which treatment options are randomly assigned (Liao et al., 2016).
Individual-level data can be summarized as
$\\{O_{0},O_{1},A_{1},O_{2},A_{2},\ldots,O_{T},A_{T},O_{T+1}\\}$ where $t$
indexes a sequence of decision points, $O_{0}$ is the baseline information,
$O_{t}$ is the information collected between time $t-1$ and $t$, and $A_{t}$
is the treatment option provided at time $t$; here we consider binary
treatment options, i.e., $A_{t}\in\\{0,1\\}$. In an MRT, $A_{t}$ is randomized
with randomization probabilities that may depend on the complete observed
history $H_{t}:=\\{O_{0},O_{1},A_{1},\ldots,A_{t-1},O_{t}\\}$, denoted
$\mathbf{p}=\\{p_{t}(A_{t}\,|\,H_{t})\\}_{t=1}^{T}$. Treatment options are
designed to impact a proximal response, denoted by $Y_{t+1}$, which is a
function of the observed history and the latest treatment, i.e.,
$Y_{t+1}=y(H_{t},A_{t})$ (Dempsey et al., 2020).

###  2.2 Estimands and Inferential Methods: A Review

The class of estimands, referred to as “causal excursion effects”, was
developed to assess whether mobile health interventions influence the proximal
health outcomes they were designed to impact (Heron and Smyth, 2010). These
time-varying effects are a function of the decision point $t$ and a set of
moderators $S_{t}$ and marginalize over all other observed and unobserved
variables (Dempsey et al., 2020; Qian et al., 2020). We provide formal
definitions using potential outcomes (Rubin, 1978; Robins, 1986).

Let $Y_{t+1}(\bar{a}_{t-1})$ denote the potential outcome for the proximal
response under treatment sequence $\bar{a}_{t-1}$. Let $O_{t}(\bar{a}_{t-1})$
denote the potential information collected between time $t-1$ and $t$. Let
$S_{t}(\bar{a}_{t-1})$ denote the potential outcome for a time-varying effect
moderator which is a deterministic function of the potential history up to
time $t$, $H_{t}(\bar{a}_{t-1})$. We consider the setting in which the
potential outcomes are i.i.d over users according to a distribution
$\mathcal{P}$, i.e.,
$\left\\{O_{t}(\bar{a}_{t-1}),Y_{t+1}(\bar{a}_{t-1})\right\\}_{t=1}^{T}\overset{\text{i.i.d}}{\sim}\mathcal{P}$.
The causal excursion effect estimand is defined as:

$\beta_{\mathbf{p}}(t;s)=\mathbb{E}_{{\mathbf{p}}}\left[Y_{t+1}\left(\bar{A}_{t-1},A_{t}=1\right)-Y_{t+1}\left(\bar{A}_{t-1},A_{t}=0\right)\,|\,S_{t}(\bar{A}_{t-1})=s\right].$
(1)

Equation (1) is defined with respect to a reference distribution $\mathbf{p}$,
i.e., the joint distribution of treatments
$\bar{A}_{t-1}:=\\{A_{1},A_{2},\dots,A_{t-1}\\}$. We follow common practice in
observational mobile health studies where analyses such as GEEs (Liang and
Zeger, 1986) are conducted marginally over $\mathbf{p}$. To express the
proximal response in terms of the observed data, we assume positivity,
consistency, and sequential ignorability (Robins, 1994, 1997):

######  Assumption 2.1.

We assume consistency, positivity, and sequential ignorability:

  * •

Consistency: For each $t\leq T$,
$\\{Y_{t+1}(\bar{A}_{t}),O_{t}(\bar{A}_{t-1}),A_{t}(\bar{A}_{t-1})\\}=\\{Y_{t+1},O_{t},A_{t}\\}$,
i.e., observed values equal the corresponding potential outcomes;

  * •

Positivity: if the joint density $\\{A_{t}=a_{t},H_{t}=h_{t}\\}$ is greater
than zero, then $P(A_{t}=a_{t}\,|\,H_{t}=h_{t})>0$;

  * •

Sequential ignorability: For each $t\leq T$, the potential outcomes
$\\{Y_{t+1}(\bar{a}_{t}),O_{t+1}(\bar{a}_{t}),A_{t+1}(\bar{a}_{t}),\dots,Y_{T+1}(\bar{a}_{T})\\}$
are independent of $A_{t}$ conditional on the observed history $H_{t}$.

Under Assumption 2.1, (1) can be re-expressed in terms of observable data:

$\beta_{\mathbf{p}}(t;s)=\mathbb{E}\left[\mathbb{E}_{{\mathbf{p}}}\left[Y_{t+1}\mid
A_{t}=1,H_{t}\right]-\mathbb{E}_{{\mathbf{p}}}\left[Y_{t+1}\mid
A_{t}=0,H_{t}\right]\mid S_{t}=s\right].$ (2)

The causal excursion effect is typically assumed to take a known linear form,
i.e., $\beta_{{\mathbf{p}}}(t;s)=f_{t}(s)^{\top}\beta^{\star}$, where
$f_{t}(s)\in\mathbb{R}^{q}$ is a feature vector comprised of a $q$-dimensional
summary of observed information depending only on state $s$ and decision point
$t$. MRTs are experimental studies with pre-specified randomization
probabilities. It is therefore common to impose the following condition:

######  Assumption 2.2.

The randomization probability $p_{t}(A_{t}|H_{t})$ is known or correctly
specified via a parametric model $p_{t}(A_{t}|H_{t};\theta)$ for
$\theta\in\mathbb{R}^{d}$.

A consistent estimator $\hat{\beta}_{n}$ can be obtained by minimizing a
weighted and centered least squares (WCLS) criterion (Boruvka et al., 2018):

$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}(Y_{t+1}-g_{t}(H_{t})^{\top}\alpha-(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\beta)^{2}\right],$
(3)

where $\mathbb{P}_{n}$ is an operator denoting the sample average,
$W_{t}=\tilde{p}_{t}(A_{t}|S_{t})/p_{t}(A_{t}|H_{t})$ is a weight where the
numerator is an arbitrary function with range $(0,1)$ that only depends on
$S_{t}$, and $g_{t}(H_{t})\in\mathbb{R}^{p}$ are $p$ control variables.
Important to this paper, the linear term $g_{t}(H_{t})^{\top}\alpha$ is a
working model for $\mathbb{E}[W_{t}Y_{t+1}|H_{t}]$, which can be viewed as a
nuisance function. A high-quality model of the nuisance function can help
reduce variance and construct more powerful test statistics. See Boruvka et
al. (2018) for more details on the estimand formulation and consistency,
asymptotic normality, and robustness properties of this method.

######  Remark 2.3.

Correct causal effect specification, i.e., $\beta_{{\bf
p}}(t;s)=f_{t}(s)^{\top}\beta^{\star}$ is not required. Instead, we can follow
prior literature (Dempsey et al., 2020; Shi et al., 2022) and interpret the
proposed linear form as a working model. Specifically, $\hat{\beta}$ is a
consistent and asymptotically normal estimator for

$\beta^{\star}=\arg\min_{\beta}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))(\beta(t;S_{t})-f_{t}(S_{t})^{\top}\beta)^{2}\right].$

Therefore, the working model can be interpreted as an $L_{2}$ projection of
the true causal excursion effect onto the space spanned by a $q$-dimensional
feature vector that only includes $t$ and $s$, denoted by
$f_{t}(s)^{\top}\beta^{\star}$ (Dempsey et al., 2020). Interpretation as a
projection or as a correctly specified causal effect can be viewed as a bias-
variance trade-off. The projection interpretation guarantees well-defined
parameter interpretation in practice.

###  2.3 Machine Learning and Causal Effect Estimation

Chernozhukov et al. (2018) provided a generic DML approach for obtaining valid
inferential statements about focal parameters, using Neyman-orthogonal scores
and cross-fitting, in settings where nuisance parameters are estimated using
ML methods. As a motivating example, here we consider the following partially
linear regression (PLR) model as in Robinson (1988):

$\displaystyle Y$
$\displaystyle=A\beta_{0}+m_{0}(X)+U,~{}~{}~{}~{}~{}~{}~{}\mathbb{E}[U|X,A]=0,$
$\displaystyle A$
$\displaystyle=p_{0}(X)+V,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathbb{E}[V|X]=0.$

Here, $Y$ is the outcome variable, $A$ represents the treatment indicator,
$X=(X_{1},...,X_{p})$ consists of baseline covariates, and $U$ and $V$ are
noise variables. The parameter of interest is $\beta_{0}$, i.e., the treatment
effect.

In many applications, the dimension of baseline covariates, $dim(X)=p$, is
large relative to the sample size $N$. To capture this, modern theoretical
analyses consider $p$ increasing with sample size. Traditional analyses that
limit the complexity of the nuisance functions $g_{0}=(m_{0},p_{0})$ will fail
in these settings. Chernozhukov et al. (2018) apply Neyman orthogonality and
sample splitting to overcome the failures of traditional methods. Suppose, for
the sake of clarity, that we randomly split the sample into two parts: a main
part of size n, with observation numbers indexed by $i\in\mathcal{I}$, and an
auxiliary part of size $N-n$, with observations indexed by
$i\in\mathcal{I}^{\complement}$. Let $\hat{g}_{0}$ be the estimator of the
nuisance parameter $g_{0}$, which is obtained using the auxiliary sample, and
the estimator of $\beta_{0}$ is obtained using the main sample and satisfying:

$\frac{1}{n}\sum_{i\in\mathcal{I}}\psi(W;\hat{\beta}_{0},\hat{g}_{0})=0,$

where $\psi$ is an orthogonalized or debiased score function, i.e., it
satisfies the property that the Gateaux derivative operator with respect to
$g$ vanishes when evaluated at the true parameter values:

$\partial_{g}\mathbb{E}\left[\psi(W;\beta_{0},g_{0})\right](g-g_{0})=0.$ (4)

Using moment conditions that satisfy Equation (4) to construct estimators and
inference procedures that are robust to mistakes in nuisance parameters has a
long history in statistics. We refer to property (4) as _Neyman orthogonality_
and to $\psi$ as the _Neyman orthogonal score function_ due to the fundamental
contributions in Neyman (1979), where this notion was introduced. The score
functions $\psi$ are not sensitive to biased estimation of $g_{0}$ in the
sense that (4) holds. Neyman orthogonality and sample splitting are the main
tools for establishing good behavior of an estimator for $\beta_{0}$.

##  3 Contributions

The objective of this paper is to advance the understanding of optimal methods
for estimating causal excursion effects. The main challenges in achieving this
are two-fold: First, in MRTs where treatments, responses, and moderators all
vary with time, there is a need to effectively incorporate supervised learning
methods and sample splitting into the causal excursion effect estimation
without introducing bias. Second, it is increasingly common for Assumption 2.2
to be violated due to unknown randomization probabilities or implemented
probabilities not matching the recorded probabilities. In these settings, the
current WCLS approach can only provide consistent estimates if the working
_linear_ outcome regression model for $\mathbb{E}[Y_{t+1}|H_{t},A_{t}]$ is
correctly specified. It is therefore important to develop an inferential
procedure that is robust to model misspecification.

This paper makes three original contributions: First, we introduce two
inferential procedures for estimating causal excursion effects that
incorporate DML techniques, called “R-WCLS” and “DR-WCLS”. To mitigate
overfitting bias, cross-fitting is employed to estimate the first stage of
data-driven plug-ins. Second, we demonstrate the $\sqrt{n}$-consistency and
asymptotic normality of the proposed estimators under regularity conditions,
while remaining agnostic to the specific supervised learning algorithm used to
learn the plug-in models. Third, we provide theoretical guarantees of double-
robustness and gains in estimation efficiency relative to the WCLS approach.

##  4 A Supervised Learning Algorithm Agnostic Approach to Moderation Analysis

We begin with a moment-based estimand, as defined in Equation (2), and assume
a parametric model for the causal excursion effect, denoted as
$\beta(t;s)=f_{t}(S_{t})^{\top}\beta^{\star}$, where
$\beta^{\star}\in\mathbb{R}^{q}$. The WCLS criterion provides a set of
estimating equations used to perform inference about the causal parameter
$\beta^{\star}$. This approach suggests that the nuisance parameter can be
expressed as a sequence of expectations
$\mathbf{g}=\\{g_{t}(H_{t})=\mathbb{E}[W_{t}Y_{t+1}|H_{t}]\\}_{t=1}^{T}$, with
a population value of $\mathbf{g}^{\star}$. To estimate these quantities, the
WCLS criterion only considers linear working models
$\\{g_{t}(H_{t})^{\top}\alpha\\}_{t=1}^{T}$.

Constructing linear working models, however, can pose a significant challenge
as researchers must pre-specify features from the high-dimensional observed
history $H_{t}$, which is a non-trivial task. To address this challenge and
increase the modeling flexibility of nuisance functions, we reframe the
estimating equation (3) in a more general form that puts no modeling
assumptions in $g_{t}(H_{t})$ and allows its dimensions to grow with sample
size. Here, we assume that $\beta^{\star}$ satisfies the moment conditions:
$\mathbb{E}[\psi(\beta^{\star};\mathbf{g}^{\star})]=0$, where
$\psi(\beta;\mathbf{g})$ is the estimating equation for $\beta$:

$\psi(\beta;\mathbf{g})=\sum_{t=1}^{T}W_{t}\Big{[}Y_{t+1}-g_{t}(H_{t})-(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\beta\Big{]}(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t}).$
(5)

We can recover WCLS by replacing $g_{t}(H_{t})$ with a linear working model
with fixed dimension, i.e., $g(H_{t})^{\top}\alpha$ for
$\alpha\in\mathbb{R}^{p}$. To ensure robustness and valid inference for
$\beta$, we require Neyman orthogonality for the estimating equation
$\psi(\beta;\mathbf{g})$ (Chernozhukov et al., 2015). The Gateaux derivative
operator with respect to $\mathbf{g}$ is:

$G(\mathbf{g})=\mathbb{E}\left[\sum_{t=1}^{T}W_{t}g_{t}(H_{t})(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right](\mathbf{g}-\mathbf{g}^{\star})=0,$
(6)

thus (5) satisfies Neyman orthogonality. Intuitively, Neyman orthogonality
implies that the moment conditions used to identify $\beta^{\star}$ are
sufficiently insensitive to the nuisance parameter estimates, allowing us to
directly plug in estimates of $\mathbf{g}^{\star}$ while still obtaining high-
quality inference for $\beta$.

We now consider estimating $g_{t}(H_{t})$. By definition,
$g_{t}(H_{t})=\mathbb{E}[W_{t}Y_{t+1}|H_{t}]$ is the conditional expectation
of the weighted proximal outcome. Let $g_{t}(H_{t},A_{t})$ denote a working
model for $\mathbb{E}[Y_{t+1}|H_{t},A_{t}]$, then we have the decomposition:
$g_{t}(H_{t})=\tilde{p}_{t}(1|S_{t})g_{t}(H_{t},1)+(1-\tilde{p}_{t}(1|S_{t}))g_{t}(H_{t},0)$.
Based on this, we propose the following _R-WCLS criterion_ , which minimizes:

$\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-\Big{(}\tilde{p}_{t}(1|S_{t})g_{t}(H_{t},1)+(1-\tilde{p}_{t}(1|S_{t}))g_{t}(H_{t},0)\Big{)}-(A_{t}-\tilde{p}_{t}(1|S_{t})))f_{t}(S_{t})^{\top}\beta\right)^{2}\right].$
(7)

The proposed approach (7) allows us to learn these conditional expectations
without pre-specifying features to construct a parametric working model. This
reflects the fact that in practice, only a small subset of the parameters are
of key scientific interest, and an analyst may prefer to be agnostic about the
nuisance parameters. As an example, we can use supervised learning methods to
train the working model $g_{t}(H_{t},A_{t})$, along with cross-fitting to
obtain informative error analysis, avoid overfitting, and optimize out-of-
sample prediction accuracy, which is particularly useful when the dimension of
the complete observed history is high. Theorem 5.2 in Section 5.2 shows that
the estimator $\hat{\beta}^{(R)}_{n}$ obtained by minimizing Equation (7) is a
consistent estimator of the parameter of interest $\beta^{\star}$ under the
assumption that the randomization probability $p_{t}(A_{t}|H_{t})$ is either
known or correctly specified.

######  Remark 4.1 (Connection to the WCLS criterion).

The R-WCLS criterion replaces $g_{t}(H_{t})^{\top}\alpha$ in the WCLS
criterion, which was a linear working model for
$\mathbb{E}[W_{t}Y_{t+1}|H_{t}]$, with a general choice of working models.
Setting $g_{t}(H_{t},A_{t})$ to be the linear working model
$g_{t}(H_{t})^{\top}\alpha+(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\beta$,
the R-WCLS criterion recovers the original WCLS criterion. Thus, (7) is a
strict generalization of (3).

######  Remark 4.2 (Connection to the R-learner).

In traditional causal inference with a single treatment $A$, fully-observed
set of confounders $X$, and outcome $Y$, a two-stage estimator, referred to as
the R-Learner, was previously proposed by Nie and Wager (2021). Beyond our
extension to the time-varying setting, there are two key distinguishing
features of R-WCLS in (7) when compared with the R-Learner. First, we focus on
estimating a low-dimensional target parameter, whereas R-learner seeks to
estimate the conditional average treatment effect and allows it to be a
complex function of baseline covariates. Second, the weight $W_{t}$ in R-WCLS
criterion implicitly depends on the propensity $\tilde{p}_{t}(1|S_{t})$, we
thereby replace the R-learner data-adaptive model for
$\mathbb{E}[W_{t}Y_{t+1}|H_{t}]$ with one for each
$\mathbb{E}[Y_{t+1}|H_{t},a]$, $a\in\\{0,1\\}$, which is invariant to
different choices of moderators $S_{t}$.

###  4.1 A Doubly-Robust Alternative

The above discussion relies on Assumption 2.2 being held. In many MRTs, one
may fail to correctly implement or collect the desired randomization
probabilities, leading to unknown randomization probabilities or uncertainty
in their recorded values. In such cases, the R-WCLS criterion in (7) can only
provide consistent estimates of $\beta^{\star}$ if the fully conditional model
for $\mathbb{E}[Y_{t+1}|H_{t},A_{t}]$ has been correctly specified:

$\mathbb{E}[Y_{t+1}|H_{t},A_{t}]=\left(\tilde{p}_{t}(1|S_{t})g_{t}^{\star}(H_{t},1)+(1-\tilde{p}_{t}(1|S_{t}))g_{t}^{\star}(H_{t},0)\right)+(A_{t}-\tilde{p}_{t}(1|S_{t})))f_{t}(S_{t})^{\top}\beta^{\star}.$

This implies that
$g^{\star}_{t}(H_{t},1)-g^{\star}_{t}(H_{t},0)=f_{t}(S_{t})^{\top}\beta^{\star}$,
where the fully conditional treatment effect depends only on the specified
moderators $S_{t}$ and the linear model is correctly specified. However, in
practice, $S_{t}$ is often a subset of the potential moderators, so this
assumption is not expected to hold. Therefore, an estimating procedure that
does not rely on a correct model specification will be preferred.

In this section, we present a novel derivation of an alternative, doubly
robust estimator called the _DR-WCLS criterion_ , which is given by:

$\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\left(\frac{W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g_{t}(H_{t},A_{t}))}{\tilde{\sigma}^{2}_{t}(S_{t})}+\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta\right)f_{t}(S_{t})\right]=0,$
(8)

where $\beta(t;H_{t})\coloneqq g_{t}(H_{t},1)-g_{t}(H_{t},0)$ is the causal
excursion effect under the fully observed history $H_{t}$, and
$\tilde{\sigma}^{2}_{t}(S_{t})\coloneqq\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))$.
Theorem 5.3 below shows that the estimator $\hat{\beta}^{(DR)}_{n}$ obtained
from solving (8) is doubly-robust, i.e., (8) will yield a consistent estimator
of $\beta^{\star}$ if either the randomization probability
$p_{t}(A_{t}|H_{t})$ _or_ the conditional expectation $g_{t}(H_{t},A_{t})$ is
correctly specified.

######  Remark 4.3 (Connection to the DR-learner).

The general DR-learner approach was first proposed by Van Der Laan and Rubin
(2006). In later research, the DR-learner, a two-stage doubly robust estimator
with a single treatment $A$, a fully observed set of confounders $X$, and an
outcome $Y$, was proposed (Kennedy, 2020). Beyond our extension to the time-
varying setting, there are two key distinguishing features of (8) when
compared with existing variants of the DR-learner. First, the causal excursion
effect is a marginal effect, so the weights $W_{t}$ and the treatment-
centering probability are dependent on the moderators, whereas the DR-learner
estimates a fully conditional causal effect. Second, time-varying treatments
and the projection interpretation (see Dempsey et al. (2020) for a detailed
discussion) in the feature space $f_{t}(S_{t})$ require the additional weights
$\tilde{\sigma}^{2}_{t}(S_{t})$ in the DR-WCLS estimating equations.

###  4.2 Connection Between R-WCLS and DR-WCLS

In recent work from Morzywolek et al. (2023), a unified framework was
presented for estimating heterogeneous treatment effects, resulting in a class
of weighted loss functions with nuisance parameters. They showed that the
R-Learner (Nie and Wager, 2021) and the DR-Learner (Kennedy, 2020) can be seen
as special cases resulting from particular weighting choices. Here, we present
a complementary viewpoint by showing a simple relationship between the two
proposed methods. We begin by adding and subtracting
$g_{t}(H_{t},A_{t})=A_{t}g_{t}(H_{t},1)+(1-A_{t})g_{t}(H_{t},0)$ from Equation
(7):

$\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-g_{t}(H_{t},A_{t})+\left(A_{t}-\tilde{p}_{t}(1|S_{t})\right)\left(\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta\right)\right)^{2}\right].$

One can then obtain an estimate of $\beta^{\star}$ by solving the following
estimating equation:

$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-g_{t}(H_{t},A_{t})\right)\left(A_{t}-\tilde{p}_{t}(1|S_{t})\right)f_{t}(S_{t})\right]+$
(9)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}\left(A_{t}-\tilde{p}_{t}(1|S_{t})\right)^{2}\left(\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta\right)f_{t}(S_{t})\right].$
(10)

Under the correct specification of the randomization probabilities, the
Gateaux derivative with respect to $\mathbf{g}$ of both terms (9) and (10)
will be 0. However, if the randomization probabilities are not specified
correctly, term (10) may not have a Gateaux derivative of 0. To address this,
we replace the stochastic term $W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))^{2}$ in
(10) with its expectation under the correct randomization probability:

$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})(\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta^{\star})f_{t}(S_{t})\right].$

After this substitution, we recover (8). And by doing so, the Gateaux
derivative with respect to $\mathbf{g}$ of both terms will no longer be
affected by the randomization probability specification. The above derivation
links the R-WCLS and DR-WCLS, showing that the doubly-robust estimators can be
constructed from R-learner methods. Finally, (7) and (8) yield estimation
procedures that are presented in Section 5.1.

##  5 Algorithm and Asymptotic Theory

This section presents the inferential algorithms for R-WCLS and DR-WCLS. In
addition, we provide the corresponding asymptotic theory that can be used for
hypothesis testing and constructing confidence intervals.

###  5.1 Algorithm

Both R-WCLS and DR-WCLS algorithms exploit the structure of (7) and (8) to
characterize the problem as a two-stage weighted regression estimation that
regresses estimated _pseudo-outcomes_ on a feature vector. Cross-fitting is
employed to obtain asymptotic theories and convergence results that are
agnostic of the supervised learning algorithm used for the estimation of the
nuisance parameters, avoiding the Donsker conditions that were prevalent in
the classic semi-parametric inference literature.

##### Step I

Let $K$ be a fixed integer. Form a $K$-fold random partition of
$\\{1,2,\dots,n\\}$ by dividing it to equal parts, each of size $n/k$,
assuming $n$ is a multiple of $k$. From each set $I_{k}$, let
$I_{k}^{\complement}$ denote the observation indices that are not in $I_{k}$.

##### Step II

Learn the appropriate working models for each fold $I_{k}$ using the
individuals in $I_{k}^{\complement}$. Let $\hat{g}_{t}^{(k)}(H_{t},A_{t})$,
$\hat{p}_{t}^{(k)}(1|H_{t})$, and $\hat{\tilde{p}}_{t}^{(k)}(1|S_{t})$ denote
the estimates for $\mathbb{E}[Y_{t+1}|H_{t},A_{t}]$,
$\mathbb{E}[A_{t}|H_{t}]$, and $\mathbb{E}[A_{t}|S_{t}]$ respectively, i.e.,
estimates of the nuisance parameters in the $k$th fold. Note that when
randomization probabilities are known, $\hat{p}_{t}(A_{t}|H_{t})$ is set equal
to $p_{t}(A_{t}|H_{t})$.

##### Step III

Construct the pseudo-outcomes and perform weighted regression estimation:

  * •

R-WCLS: For individual $j$ at time $t$, define the pseudo-outcome:

$\tilde{Y}^{(R)}_{t+1,j}:=Y_{t+1,j}-\hat{g}^{(k)}_{t}(H_{t,j},A_{t,j})+\left(A_{t,j}-\hat{\tilde{p}}_{t}^{(k)}(1|S_{t,j})\right)\left(\hat{g}_{t}^{(k)}(H_{t,j},1)-\hat{g}_{t}^{(k)}(H_{t,j},0)\right),$

where $j\in I_{k}$. Then regress $\tilde{Y}^{(R)}_{t+1}$ on
$(A_{t}-\hat{\tilde{p}}_{t}^{(k)}(1|S_{t}))f_{t}(S_{t})^{\top}\beta$ with
weights
$\hat{W}_{t}^{(k)}=\hat{\tilde{p}}_{t}^{(k)}(A_{t}|S_{t})/\hat{p}_{t}^{(k)}(A_{t}|H_{t})$
to obtain estimate $\hat{\beta}_{n}^{(R)}$.

  * •

DR-WCLS: For individual $j$ at time $t$, define the pseudo-outcome:

$\tilde{Y}^{(DR)}_{t+1,j}:=\frac{\hat{W}^{(k)}_{t,j}(A_{t,j}-\hat{\tilde{p}}^{(k)}_{t}(1|S_{t,j}))(Y_{t+1,j}-\hat{g}_{t}^{(k)}(H_{t,j},A_{t,j}))}{\hat{\tilde{p}}^{(k)}_{t}(1|S_{t,j})(1-\hat{\tilde{p}}^{(k)}_{t}(1|S_{t,j}))}+\left(\hat{g}_{t}^{(k)}(H_{t,j},1)-\hat{g}_{t}^{(k)}(H_{t,j},0)\right),$

where $j\in I_{k}$. Then regress $\tilde{Y}^{(DR)}_{t+1}$ on
$f_{t}(S_{t})^{\top}\beta$ with weights
$\hat{\tilde{p}}^{(k)}_{t}(1|S_{t})(1-\hat{\tilde{p}}^{(k)}_{t}(1|S_{t}))$ to
obtain $\hat{\beta}_{n}^{(DR)}$.

######  Remark 5.1.

Without sample splitting, the estimated nuisance functions are now correlated,
hence introducing spurious correlation in Step III. Typical approaches would
constrain the function $g(\cdot)$ in Step II to belong to a function class
with relatively simple statistical complexity, typically referred to as a
Donsker function class. Then the Step III estimate is $\sqrt{n}$-consistent
and asymptotically normal. Chen et al. (2022) shows that neither the sample
splitting nor Donsker property is required if the estimate $\hat{g}(\cdot)$
satisfies leave-one-out stability properties and the moment function satisfies
the weak mean-squared-continuity property of Chernozhukov et al. (2021). This
allows for sample reuse, which can benefit moderately sized sample regimes.
Here we aim to stay agnostic about the choice of $g(\cdot)$, but we consider
extensions that do not require sample splitting as important future work.

###  5.2 Asymptotic Properties

Here, we demonstrate the asymptotic theory for both R-WCLS and DR-WCLS
estimators obtained using the algorithm described above. All asymptotic
statements assume $\hat{p}_{t}(A_{t}|H_{t})$ is bounded away from 0 and 1, $T$
and $K$ both finite and fixed, and $n$ increasing to infinity.

######  Theorem 5.2 (Asymptotic property of R-WCLS estimator).

Under Assumption 2.2, and given invertibility and moment conditions, the
estimator $\hat{\beta}^{(R)}_{n}$ that minimizes (7) is consistent and
asymptotically normal such that
$\sqrt{n}(\hat{\beta}^{(R)}_{n}-\beta^{\star})\rightarrow\mathcal{N}(0,\Sigma_{R})$,
where $\Sigma_{R}$ is defined in Appendix A. In particular, with the algorithm
outlined in Section 5.1, $\Sigma_{R}$ can be consistently estimated by:

$\left[\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left\\{\dot{m}(\hat{\beta},\hat{\eta}_{k})\right\\}\right]^{-1}\times\left[\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left\\{m(\hat{\beta},\hat{\eta}_{k})m(\hat{\beta},\hat{\eta}_{k})^{\top}\right\\}\right]\times\left[\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left\\{\dot{m}(\hat{\beta},\hat{\eta}_{k})\right\\}\right]^{-1},$
(11)

where
$m(\hat{\beta},\hat{\eta}_{k})=\sum_{t=1}^{T}\hat{W}^{(k)}_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\hat{\tilde{p}}_{t}^{(k)}(1|S_{t}))f_{t}(S_{t})^{\top}\hat{\beta}^{(R)}_{n}\right)(A_{t}-\hat{\tilde{p}}_{t}^{(k)}(1|S_{t}))f_{t}(S_{t})$,
$\dot{m}(\hat{\beta},\hat{\eta}_{k})=\sum_{t=1}^{T}\hat{\tilde{p}}_{t}^{(k)}(1|S_{t})(1-\hat{\tilde{p}}_{t}^{(k)}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}$,
and $\mathbb{P}_{n,k}\\{\bullet\\}$ refers to the empirical average within
fold $k$.

If Assumption 2.2 is violated, using the R-WCLS criterion in equation (7) will
only produce consistent estimates if the fully conditional model for
$\mathbb{E}[Y_{t+1}|H_{t},A_{t}]$ is correctly specified, which is difficult
to achieve in practice due to the complexity of the true data-generating
process. Therefore, the DR-WCLS estimator is especially valuable for
safeguarding against model misspecification. The asymptotic properties of the
DR-WCLS estimator are as follows:

######  Theorem 5.3 (Asymptotic property of DR-WCLS estimator).

Given invertibility and moment conditions, the estimator
$\hat{\beta}^{(DR)}_{n}$ that solves (8) is subject to an error term, which
(up to a multiplicative constant) is bounded above by:

$\mathbf{\hat{B}}=\mathbb{E}\left[\sum_{t=1}^{T}\sum_{a\in\\{0,1\\}}\left\|p_{t}(A_{t}=1|H_{t})-\hat{p}_{t}(A_{t}=1|H_{t})\right\|\left\|g_{t}(H_{t},a)-\hat{g}_{t}(H_{t},a)\right\|\right],$
(12)

where $\left\|X\right\|:=(\mathbb{E}[X^{2}])^{1/2}$. If
$\mathbf{\hat{B}}=o_{p}(n^{-1/2})$, then $\hat{\beta}^{(DR)}_{n}$ is
consistent and asymptotically normal such that
$\sqrt{n}(\hat{\beta}^{(DR)}_{n}-\beta^{\star})\rightarrow\mathcal{N}(0,\Sigma_{DR})$,
where $\Sigma_{DR}$ is defined in Appendix B. In particular, with the
algorithm outlined in Section 5.1, $\Sigma_{DR}$ can be consistently estimated
by formula (11) with
$m(\hat{\beta},\hat{\eta}_{k})=\sum_{t=1}^{T}\hat{\tilde{p}}_{t}^{(k)}(1|S_{t})(1-\hat{\tilde{p}}_{t}^{(k)}(1|S_{t}))(\tilde{Y}^{(DR)}_{t+1}-f_{t}(S_{t})\hat{\beta}^{(DR)}_{n})f_{t}(S_{t})$,
and
$\dot{m}(\hat{\beta},\hat{\eta}_{k})=\sum_{t=1}^{T}\hat{\tilde{p}}_{t}^{(k)}(1|S_{t})(1-\hat{\tilde{p}}_{t}^{(k)}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}$.

It follows that $\hat{\beta}^{(DR)}_{n}$ is doubly robust since it is
consistent when either (1) the treatment model is correctly specified or (2)
the conditional model is correctly specified. Importantly, the model-agnostic
error bound applies to arbitrary first-stage estimators. The bound
$\mathbf{\hat{B}}$ on the DR-WCLS estimator error shows that it can only
deviate from $\beta^{\star}$ by at most a (smoothed) product of errors in the
estimation of treatment propensities and conditional expectation of outcomes,
thus allowing faster rates for estimating the causal effect even when the
nuisance estimates converge at slower rates. For detailed proofs of Theorems
5.2 and 5.3, please refer to Appendices A and B respectively.

######  Remark 5.4.

A variety of flexible options for nuisance estimates are available to attain
the convergence rate of $o_{p}(n^{-1/2})$. For example, if
$\left\|p_{t}(1|H_{t})-\hat{p}_{t}(1|H_{t})\right\|=o_{p}(n^{-1/4})$ and
$\left\|g_{t}(H_{t},a)-\hat{g}_{t}(H_{t},a)\right\|=o_{p}(n^{-1/4})$, then the
product term is $o_{p}(n^{-1/2})$ and is thus asymptotically negligible. This
occurs when both $\hat{g}(H_{t},a)$ and $\hat{p}_{t}(a|H_{t})$ are based on
correctly specified parametric models, but also achievable for many ML methods
under structured assumptions on the nuisance parameters, for example,
regularized estimators such as the Lasso, and random forest (Chernozhukov et
al., 2018; Athey et al., 2018). Worth noticing, in this setting, completely
nonparametric estimators are usually not an option as they tend to converge at
rates slower than $n^{-1/4}$ unless there are strong smoothness or sparsity
assumptions in place (Kennedy, 2016).

###  5.3 Estimation Efficiency Comparison

Previous research (Qian et al., 2020; Shi et al., 2022) has demonstrated that
a locally efficient, semiparametric estimator for the fully conditional causal
effect (i.e., $S_{t}=H_{t}$) can be derived based on semiparametric efficiency
theory (Robins, 1994; Newey, 1990; Tsiatis, 2007). These findings provide the
motivation for the development of the proposed methods described above. In
this section, we investigate the relative efficiency between the proposed
estimators and WCLS, with a focus on the situation where $S_{t}\subset H_{t}$.

The term “more efficient” here means one method achieves lower asymptotic
variance than another method for any linear combination of the causal
parameter estimates, i.e., the asymptotic variance of $c^{\top}\hat{\beta}$ is
smaller for any $c\in\mathbb{R}^{q}$. This is equivalent to the difference
between the asymptotic variance matrices being negative semidefinite.

####  5.3.1 An Augmented R-WCLS Estimator

To make full use of the estimates $g_{t}(H_{t},A_{t})$, we propose an
augmented version of the R-WCLS criterion. According to Equation (7),
$g_{t}(H_{t},A_{t})$ is only included to construct the plug-in estimator
$g_{t}(H_{t})$, but the difference between $g_{t}(H_{t},1)$ and
$g_{t}(H_{t},0)$, i.e., the causal excursion effect under fully observed
history, is not incorporated into the estimating equation. Therefore, an
augmented R-WCLS criterion can efficiently use this information in the
following manner:

$0=\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-g_{t}(H_{t})-\left(A_{t}-\tilde{p}_{t}(1|S_{t})\right)\left(f_{t}(S_{t})^{\top}\beta+\Delta_{t}^{\perp}\right)\right)\left(A_{t}-\tilde{p}_{t}(1|S_{t})\right)f_{t}(S_{t})\right],$
(13)

where $\Delta_{t}^{\perp}$ denotes the projection of the fully conditional
causal effect $\beta(t;H_{t})$ onto the orthogonal complement of
$f_{t}(S_{t})$. The definition of the orthogonal complement can be found in
(32), while the details on constructing $\Delta_{t}^{\perp}$ are provided in
Appendix D.1 for further reference. The estimator obtained from solving
Equation (13) has similar asymptotic properties as R-WCLS and it leads to the
following lemma. Proof can be found in Appendix D.2.

######  Lemma 5.5.

Let $\hat{\beta}^{(AR)}_{n}$ denote the augmented R-WCLS estimator obtained
from solving Equation (13). Under Assumption 2.2, given invertibility and
moment conditions, $\hat{\beta}^{(AR)}_{n}$ is consistent and asymptotically
normal such that
$\sqrt{n}(\hat{\beta}^{(AR)}_{n}-\beta^{\star})\rightarrow\mathcal{N}(0,\Sigma_{AR})$,
where $\Sigma_{AR}$ is defined in Appendix D.2.

####  5.3.2 Efficiency Improvement

We compare the estimation efficiency of the meta-learning estimators proposed
with that of the existing WCLS methods. Before presenting the theorem, we
first list the additional required assumptions :

######  Assumption 5.6.

In addition to Assumption 2.2, we make the following assumptions:

  1. (a)

The residual $e_{t}\coloneqq Y_{t+1}-g(H_{t},A_{t})$ is uncorrelated with
future states given history $H_{t}$ and treatment $A_{t}$, i.e.,
$\mathbb{E}[e_{t}f_{t^{\prime}}(S_{t^{\prime}})\Delta^{\perp}_{t^{\prime}}|H_{t},A_{t}]=0,~{}t<t^{\prime}$;

  2. (b)

The estimator $\hat{g}(H_{t},A_{t})$ is consistent for the true conditional
expectation $g(H_{t},A_{t})$.

Assumption 5.6 a implies that the residuals do not convey any additional
information beyond that which is already contained in the observed history.
Such an assumption would hold when considering Markov Decision Processes
(MDPs), where the current state and treatment determine future states.
Assumption 5.6b is derived from the discussion in Theorem 5.3, where we
require the error term to converge at a rate of $o_{p}(n^{-1/2})$. This
assumption, in conjunction with Assumption 2.2, ensures the consistency of the
DR-WCLS estimator.

######  Theorem 5.7 (Efficiency gain over WCLS estimator).

Under Assumption 2.2 and given invertibility and moment conditions:

  1. (a)

if Assumption 5.6a also holds, the augmented R-WCLS estimator is guaranteed to
be at least as efficient as WCLS;

  2. (b)

if Assumption 5.6b holds, the DR-WCLS estimator is guaranteed to be at least
as efficient as the R-WCLS estimator.

Theorem 5.7a requires Assumption 5.6a as a sufficient condition when
considering a smoothed model for the causal excursion effect
$f_{t}(S_{t})^{\top}\beta^{\star}$ over time. However, Assumption 5.6a is not
necessary, since efficiency gains can always be guaranteed if the causal
excursion effect is modeled nonparametrically over time, i.e.,
$f_{t}(S_{t})^{\top}\beta_{t}^{\star}$. Theorem 5.7b indicates further
asymptotic efficiency gains by employing our proposed doubly-robust
alternative. Bang and Robins (2005) do comment that doubly-robust estimators
(i.e., DR-WCLS) may be less efficient in finite samples than the inverse
probability of treatment weighting (IPTW) estimators (i.e., R-WCLS) under
extremely strong model misspecification. In the current context, if
$\hat{p}(1|H_{t})$ is based on a correctly specified parametric model so that
$\left\|p_{t}(1|H_{t})-\hat{p}_{t}(1|H_{t})\right\|=O_{p}(n^{-1/2})$, then
$\hat{g}_{t}(H_{t},a_{t})$ need only to be consistent, i.e,
$\left\|g_{t}(H_{t},a_{t})-\hat{g}_{t}(H_{t},a_{t})\right\|=o_{p}(1)$ for the
DR-WCLS estimator to be least asymptotically as efficient as the R-WCLS
estimator. Our model-agnostic approach reduces the risk of severe model
misspecification. The detailed proof is provided in Appendix E.1 and E.2.

##  6 Extensions

###  6.1 Missing Data

In mHealth studies, it is common for both the proximal outcome $Y_{t+1}$ and
elements of the history $H_{t}$ to be missing. In the case study of a 6-month
MRT on medical interns presented in Section 8, for example, the proximal
outcomes are self-reported mood score and step count. Self-reports are often
missing due to non-response, while step count can be missing due to
individuals not wearing the wrist sensors. Current methods are lacking for
addressing missing data in the context of MRT (Boruvka et al., 2018; Dempsey
et al., 2020; Qian et al., 2020).

Here we extend the DR-WCLS criterion to be robust to missing data.
Specifically, we consider two types of missingness: (1) in the outcome
$Y_{t+1}$ and (2) in the observed history $H_{t}$ being used in the supervised
learning algorithm. We do _not_ consider missingness in the moderator set
$S_{t}$ which we assume is completely observed, but consider this important
future work. Let $R_{t}$ be the binary indicator of whether the proximal
outcome $Y_{t+1}$ is observed ($R_{t}=1$) or not ($R_{t}=0$) at decision time
$t$, and $R_{t}(\bar{a}_{t})$ denotes the potential observation status.
Clearly, missingness is a post-treatment variable and therefore we require
additional assumptions:

######  Assumption 6.1.

We assume consistency, missing at random, and positivity:

  * •

Consistency: For each $t\leq T$, $R_{t}(\bar{A}_{t})=R_{t}$, i.e., the
observed missing data indicator is equal to the corresponding potential
outcome observation status;

  * •

Missing at random: For each $t\leq T$, $R_{t}(\bar{a}_{t})$ is independent of
$A_{t}$ conditional on the observed history $H_{t}$;

  * •

Positivity: if the joint density $\\{R_{t}=r_{,}H_{t}=h_{t},A_{t}=a_{t}\\}$ is
greater than zero, then $p(R_{t}=1|H_{t},A_{t})=p(R_{t}|H_{t})>0$.

Under Assumption 6.1, we can derive a doubly robust extension for missing data
by augmenting the DR-WCLS criterion:

$\begin{split}\mathbb{P}_{n}\bigg{[}\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\bigg{(}&\frac{{\bf
1}(R_{t}=1)}{p(R_{t}|H_{t})}\frac{W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g_{t}(H_{t},A_{t}))}{\tilde{\sigma}^{2}_{t}(S_{t})}+\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta\bigg{)}f_{t}\bigg{]}=0.\end{split}$
(14)

Equation (14) is equal to (8) except that we multiply the first term by the
inverse probability of missing data. As the data-missing mechanism is a
complex nuisance function, it can also be considered part of the meta-learning
algorithm. Theorem 5.3 can be extended to the current setting, leading to
Corollary 6.1.1. See Appendix F for the proofs.

######  Corollary 6.1.1.

(Asymptotic property for DR-WCLS estimator with missing data) Under Assumption
6.1, given invertibility and moment conditions, the estimator
$\hat{\beta}_{n}$ that solves (14) is subject to an error term, which (up to a
multiplicative constant) is bounded above by:

$\mathbf{\hat{B}}^{R}=\mathbb{E}\left[\sum_{t=1}^{T}\sum_{a\in\\{0,1\\}}\left\|p(R_{t}=1|H_{t})p_{t}(a|H_{t})-\hat{p}(R_{t}=1|H_{t})\hat{p}_{t}(a|H_{t})\right\|\left\|g_{t}(H_{t},a)-\hat{g}_{t}(H_{t},a)\right\|\right].$
(15)

If we further assume $\mathbf{\hat{B}}^{R}=o_{p}(n^{-1/2})$, $\hat{\beta}_{n}$
is consistent and asymptotically normal such that
$\sqrt{n}(\hat{\beta}_{n}-\beta^{\star})\rightarrow\mathcal{N}(0,\Sigma^{R}_{DR})$,
where $\Sigma^{R}_{DR}$ is defined in Appendix F.

###  6.2 Lagged Effects

Beyond the interest on proximal outcomes, additional attention has been paid
to lagged outcomes defined over future decision points with a fixed window
length $\Delta>1$, denoted as $Y_{t,\Delta}$, which is a known function of the
observed history and latest treatment:
$Y_{t,\Delta}=y(H_{t+\Delta-1},A_{t+\Delta-1})$. In practice, $\Delta$ is
explicitly chosen to avoid the curse of the horizon problem (Dempsey et al.,
2020). While this has been true to date, we acknowledge that larger $\Delta$
will be more common as MRT data sets grow in size and these longer-term
outcomes become of primary interest. Under Assumption 2.1, the causal estimand
for lagged effect can be expressed in terms of observable data (Shi et al.,
2022):

$\beta_{\mathbf{p},\pi}(t+\Delta;s)=\mathbb{E}\left[\mathbb{E}_{{\mathbf{p}}}\left[W_{t,\Delta-1}Y_{t,\Delta}\mid
A_{t}=1,H_{t}\right]-\mathbb{E}_{{\mathbf{p}}}\left[W_{t,\Delta-1}Y_{t,\Delta}\mid
A_{t}=0,H_{t}\right]\mid S_{t}=s\right],$ (16)

where
$W_{t,u}=\prod_{s=1}^{u}\pi_{t}(A_{t+s}|H_{t+s})/p_{t}(A_{t+s}|H_{t+s})$, with
$W_{t,0}=1$. Here, we assume the reference distribution for treatment
assignments from $t+1$ to $t+\Delta-1$ ($\Delta>1$) is given by a
randomization probability generically represented by
$\\{\pi_{u}(a_{u}|H_{u})\\}_{u=t+1}^{t+\Delta-1}$. This generalization
contains previous definitions such as lagged effects (Boruvka et al., 2018)
where $\pi_{u}=p_{u}$ and deterministic choices such as
$a_{t+1:(t+\Delta-1)}={\bf 0}$ (Dempsey et al., 2020; Qian et al., 2020),
where $\pi_{u}={\bf 1}\\{a_{u}=0\\}$ and ${\bf 1}\\{\cdot\\}$ is the indicator
function.

A brief discussion in Shi et al. (2022) presented an approach to improve the
estimation efficiency of the lagged effect and alleviate the curse of the
horizon (Liu et al., 2018). Specifically, it was shown that an optimal
estimating function will be orthogonal to the score functions for the
treatment selection probabilities (Bickel et al., 1993). This implies the
estimator may be improved by replacing the estimating equation by itself minus
its projection on the score functions for the treatment selection
probabilities (Murphy et al., 2001). This can be done in the case of the DR-
WCLS estimating equation as follows:

$\begin{split}\mathbb{P}_{n}\Bigg{[}\sum_{t=1}^{T-\Delta+1}\Bigg{[}&W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))\Bigg{(}W_{t,\Delta-1}\left(Y_{t,\Delta}-g_{t+\Delta-1}(H_{t+\Delta-1},A_{t+\Delta-1})\right)\\\
&-\sum_{u=0}^{\Delta-2}W_{t,u}\left(g_{t+u}(H_{t+u},A_{t+u})-\sum_{a_{t+u+1}}\pi(a_{t+u+1}|H_{t+u+1})g_{t+u+1}(H_{t+u+1},a_{t+u+1})\right)\Bigg{)}\\\
&+\tilde{\sigma}_{t}^{2}(S_{t})\left(\beta(t+\Delta,H_{t})-f_{t}(S_{t})^{\top}\beta\right)\Bigg{]}f_{t}(S_{t})^{\top}\Bigg{]}=0,\end{split}$
(17)

where $g_{t+u}(H_{t+u},A_{t+u})$ is a working model for
$\mathbb{E}[W_{t+u+1:t+\Delta-1}Y_{t,\Delta}|H_{t+u},A_{t+u}]$. Specifically,
$g_{t+\Delta-1}(H_{t+\Delta-1},A_{t+\Delta-1})=\mathbb{E}[Y_{t,\Delta}|H_{t+\Delta-1},A_{t+\Delta-1}]$,
and
$\mathbb{E}[g_{t+u-1}(H_{t+u-1},A_{t+u-1})]=\mathbb{E}\left[\sum_{a_{t+u}}\pi_{t+u}(a_{t+u}|H_{t+u})g_{t+u}(H_{t+u},a_{t+u})\right]$.
The parameterized linear working model of the conditional expectation
$g_{t+u}(H_{t+u},A_{t+u})$ in Murphy et al. (2001) can be improved by
leveraging supervised learning algorithms to construct data-adaptive
estimates. Based on the estimating Equation (17), the estimator
$\hat{\beta}^{\Delta}_{n}$ has the following property:

######  Corollary 6.1.2 (Asymptotic property for DR-WCLS estimator for lagged
outcomes).

Given invertibility and moment conditions, the estimation
$\hat{\beta}^{\Delta}_{n}$ obtained by solving Equation (17) is subject to an
error term, which is (up to a multiplicative constant) bounded above by
$\sum_{u=0}^{\Delta-1}\mathbf{\hat{B}}_{u}$, where

$\begin{split}\mathbf{\hat{B}}_{u}=\mathbb{E}\Bigg{[}\sum_{t=1}^{T-u}\sum_{a_{t+u}}&\left\|\hat{p}_{t+u}(a_{t+u}|H_{t+u})-p_{t+u}(a_{t+u}|H_{t+u})\right\|\left\|\hat{g}_{t+u}(H_{t+u},a_{t+u})-g_{t+u}(H_{t+u},a_{t+u})\right\|\Bigg{]}.\end{split}$
(18)

If we assume that $\mathbf{\hat{B}}_{u}=o_{p}(n^{-1/2})$, the estimator
$\hat{\beta}^{\Delta}_{n}$ is consistent and asymptotically normal such that
$\sqrt{n}(\hat{\beta}^{\Delta}_{n}-\beta^{\star})\rightarrow\mathcal{N}(0,\Sigma^{\Delta}_{DR})$,
where $\Sigma^{\Delta}_{DR}$ is defined in Appendix G.

When $\Delta$ is large, correctly specifying the conditional expectation model
$g_{t}(H_{t},A_{t})$ is particularly useful to avoid the variance estimation
growing exponentially due to the weight $W_{t,\Delta}$, thus offering a remedy
for the curse of the horizon.

###  6.3 Binary Outcomes

Qian et al. (2020) proposed an estimator of the marginal excursion effect
(EMEE) by adopting a log relative risk model to examine whether a particular
time-varying intervention has an effect on a longitudinal binary outcome. The
causal excursion effect is defined by:

$\displaystyle\beta_{\mathbf{p}}(t;s)$
$\displaystyle=\log\frac{\mathbb{E}\left[Y_{t+1}(\bar{A}_{t-1},1)\,|\,S_{t}(\bar{A}_{t-1})=s\right]}{\mathbb{E}\left[Y_{t+1}(\bar{A}_{t-1},0)\,|\,S_{t}(\bar{A}_{t-1})=s\right]}$
(19)
$\displaystyle=\log\frac{\mathbb{E}\left[\mathbb{E}\left[Y_{t+1}\,|\,{A_{t}=1,H_{t}}\right]\,|\,S_{t}=s\right]}{\mathbb{E}\left[\mathbb{E}\left[Y_{t+1}\,|\,{A_{t}=0,H_{t}}\right]\,|\,S_{t}=s\right]}.$
(20)

Assuming $\beta_{{\bf p}}(t;s)=f_{t}(s)^{\top}\beta^{\star}$, where
$f_{t}(s)\in\mathbb{R}^{q}$ is a feature vector of a $q$-dimension and only
depends on state $s$ and decision point $t$, a consistent estimator for
$\beta^{*}$ can be obtained by solving a set of weighted estimating equations:

$\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}e^{-A_{t}f_{t}(S_{t})^{\top}\beta}\left(Y_{t+1}-e^{g_{t}(H_{t})^{\top}\alpha+A_{t}f_{t}(S_{t})^{\top}\beta}\right)\left(\begin{array}[]{c}g_{t}(H_{t})\\\
(A_{t}-\tilde{p}_{t}(1\mid S_{t}))f_{t}(S_{t})\end{array}\right)\right]=0.$
(21)

See Qian et al. (2020) for more details on the estimand formulation and
consistency, asymptotic normality, and robustness properties of the EMEE
estimation method.

Based on Equation (21), we propose a doubly-robust alternative of EMEE, termed
“DR-EMEE”. A doubly robust estimator for the log-relative risk is constructed
by solving the following set of estimating equations:

$\begin{split}\mathbb{P}_{n}\bigg{[}\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\bigg{(}&\frac{W_{t}e^{-A_{t}f_{t}(S_{t})^{\top}\beta}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g_{t}(H_{t},A_{t}))}{\tilde{\sigma}^{2}_{t}(S_{t})}+\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}e^{-f_{t}(S_{t})^{\top}\beta}g_{t}(H_{t},1)-g_{t}(H_{t},0)\bigg{)}f_{t}\bigg{]}=0.\end{split}$
(22)

######  Corollary 6.1.3 (Asymptotic property for DR-EMEE estimator).

Upon correctly specifying either conditional expectation model
$g_{t}(H_{t},A_{t})$ or treatment randomization probability
$p_{t}(A_{t}|H_{t})$, given invertibility and moment conditions, the estimator
$\hat{\beta}_{n}$ obtained from solving Equation (22) is consistent and
asymptotically normal such that
$\sqrt{n}(\hat{\beta}_{n}-\beta^{\star})\rightarrow\mathcal{N}(0,\Sigma^{b}_{DR})$,
where $\Sigma^{b}_{DR}$ is defined in Appendix H.

##  7 Simulation

###  7.1 Simulation Setup

To evaluate the proposed estimator, we extend the simulation setup in Boruvka
et al. (2018). We first present a base data generation model. Consider an MRT
with a known randomization probability, and $g(H_{t})$ in the generative model
is a complex function of high-dimensional history information $H_{t}$. Let
$S_{t}\in\\{-1,1\\}$ denote a single state variable, which is an effect
moderator, and $S_{t}\subset H_{t}$. We have the generative model as:

$Y_{t,j}=g_{t}(H_{t})+\big{(}A_{t,j}-p_{t}(1|H_{t})\big{)}(\beta_{10}+\beta_{11}S_{t,j})+e_{t,j}.$
(23)

The randomization probability is
$p_{t}(1|H_{t})=\text{expit}(\eta_{1}A_{t-1,j}+\eta_{2}S_{t,j})$ where
$\text{expit}(x)=(1+\exp(-x))^{-1}$; the state dynamics are given by
$\mathbb{P}(S_{t,j}=1|A_{t-1},H_{t-1})=1/2$ with $A_{0}=0$, and the
independent error term satisfies $e_{t,j}\sim\mathcal{N}(0,1)$ with
$\text{Corr}(e_{u,j},e_{t,j^{\prime}})=1[j=j^{\prime}]0.5^{|u-t|/2}$. As in
Boruvka et al. (2018), we set $\eta_{1}=-0.8,\eta_{2}=0.8,\beta_{10}=-0.2$,
and $\beta_{11}=0.2$. The marginal proximal effect is equal to
$\beta_{10}+\beta_{11}\mathbb{E}\left[S_{t,j}\right]=\beta_{10}=-0.2$. The
marginal treatment effect is thus constant in time and is given by
$\beta_{0}^{\star}=\beta_{10}=-0.2$.

In the following, we set the complex function $g(H_{t})$ as a decision tree,
and the flow chart Figure 2 in Appendix I visualizes the decision-making
process as well as the outcomes. We consider the estimation of the fully
marginal proximal treatment effect, thus $f_{t}(S_{t})=1$ in Equation (23)
(i.e., $S_{t}=\emptyset$). The results below report average point estimate
(Est), standard error (SE), and 95% confidence interval coverage probabilities
(CP) across 1000 replicates. Here, we report results with $N=100$ showing the
relative advantage of R-WCLS and DR-WCLS over WCLS.

###  7.2 Estimation Methods

##### WCLS

The estimation model in Boruvka et al. (2018) assumes a linear function for
the control variables, i.e., $g(H_{t};\alpha)=g(H_{t})^{\top}\alpha$. This
method is guaranteed to produce an unbiased estimate with a valid confidence
interval. Thus, it will be used as a reference for comparison of the
estimation results from the following method.

##### R-WCLS

This estimation model incorporates modern supervised learning techniques to
construct a $\hat{g}(H_{t})$ for the control variables. The plug-in estimators
$\hat{g}(H_{t},1)$ and $\hat{g}(H_{t},0)$ are learned separately with
corresponding treatment assignments in training data. As the number of control
variables is relatively small in this simulation study, we will use random
forests. The plug-in estimator is then
$\hat{g}(H_{t})=\tilde{p}_{t}(1|S_{t})\hat{g}(H_{t},1)+(1-\tilde{p}_{t}(1|S_{t}))\hat{g}(H_{t},0)$.

##### DR-WCLS

The plug-in estimators $\hat{g}(H_{t},1)$ and $\hat{g}(H_{t},0)$ can be
obtained the same way as above, and using these two estimators, we can get
$\hat{\beta}(t;H_{t})=\hat{g}(H_{t},1)-\hat{g}(H_{t},0)$.

###  7.3 Simulation Results

Table 1 reports the simulation results. “%RE gain” indicates the percentage of
times we achieve an efficiency gain out of 1000 Monte Carlo replicates. “mRE”
stands for the average relative efficiency, and “RSD” represents the relative
standard deviation between two estimates. The proposed R-WCLS and DR-WCLS
methods significantly improve the efficiency of the WCLS when estimating the
fully marginal causal effect. In addition, we find that mRE varies with
$\beta_{11}$. R-WCLS has higher mRE than DR-WCLS when $\beta_{11}$ is small,
and this reverses when $\beta_{11}$ increases. In our simulation, $\beta_{11}$
being large indicates that an important moderator $S_{t,j}$ was not included
in the causal effect model (i.e., $f_{t}(S_{t})^{\top}\beta=\beta_{0}$).
Therefore, when model misspecification occurs, DR-WCLS shows better
performance than R-WCLS.

Table 1: Fully marginal causal effect estimation efficiency comparison.  Method | $\beta_{11}$ | Est | SE | CP | %RE gain | mRE | RSD
---|---|---|---|---|---|---|---
WCLS | 0.2 | -0.198 | 0.049 | 0.946 | - | - | -
0.5 | -0.195 | 0.050 | 0.945 | - | - | -
0.8 | -0.193 | 0.053 | 0.951 | - | - | -
R-WCLS | 0.2 | -0.200 | 0.044 | 0.950 | 100% | 1.231 | 1.260
0.5 | -0.199 | 0.045 | 0.944 | 100% | 1.218 | 1.255
0.8 | -0.200 | 0.048 | 0.956 | 99.9% | 1.203 | 1.236
DR-WCLS | 0.2 | -0.200 | 0.045 | 0.954 | 99.7% | 1.216 | 1.249
0.5 | -0.199 | 0.045 | 0.947 | 99.9% | 1.228 | 1.261
0.8 | -0.200 | 0.047 | 0.954 | 99.7% | 1.254 | 1.282

##  8 Intern Health Study: A Worked Example

The Intern Health Study (IHS) is a 6-month micro-randomized trial on medical
interns (NeCamp et al., 2020), which aimed to investigate when to provide
mHealth interventions to individuals in stressful work environments to improve
their behavior and mental health. In this section, we evaluate the
effectiveness of targeted notifications in improving individuals’ moods and
step counts. The exploratory and MRT analyses conducted in this paper focus on
weekly randomization, thus, an individual was randomized to receive mood,
activity, sleep, or no notifications with equal probability ($1/4$ each) every
week. We choose the outcome $Y_{t+1,j}$ as the self-reported mood score (a
Likert scale taking values from 1 to 10) and step count (cubic root) for
individual $j$ in study week $t$.

Missing data occurred throughout the trial when interns did not complete the
self-reported mood survey or were not wearing their assigned Fitbit wrist-worn
device; thus, multiple imputation was originally used to impute missing daily
data. See NeCamp et al. (2020) for further details. The following analysis is
based on one of the imputed data sets. The data set used in the analyses
contains 1562 participants. The average weekly mood score when a notification
is delivered is 7.14, and 7.16 when there is no notification; The average
weekly step count (cubic root) when a notification is delivered is 19.1, and
also 19.1 when there is no notification. In Section 8.1 and 8.2, we evaluate
the targeted notification treatment effect for medical interns using our
proposed methods and WCLS.

###  8.1 Comparison of the Marginal Effect Estimation

First, we are interested in assessing the fully marginal excursion effect
(i.e., $\beta(t)=\beta_{0}^{\star}$). For an individual $j$, the study week is
coded as a subscript $t$. $Y_{t+1,j}$ is the self-reported mood score or step
count (cubic root) for individual $j$ in study week $t+1$. $A_{t}$ is defined
as the specific type of notification that targets improving the outcome. For
example, if the outcome is the self-reported mood score, sending mood
notifications would be the action, thus, $\mathbb{P}(A_{t}=1)=0.25$. We
analyze the marginal causal effect $\beta_{0}$ of the targeted notifications
on self-reported mood score and step count using the following model for WCLS:

$Y_{t+1,j}\sim g_{t}(H_{t,j})^{\top}\alpha+(A_{t,j}-\tilde{p}_{t})\beta_{0}.$

The term $g_{t}(H_{t})^{\top}\alpha$ represents a linear working model of
prognostic control variables including two baseline characteristics, study
week $t$, and the prior week’s outcome $Y_{t,j}$. For R-WCLS and DR-WCLS
methods, we include a total of 12 control variables and use random forests to
construct the plug-in estimators $\hat{g}(H_{t},A_{t})$ as described in
Section 7.2. For a detailed description of the control variables, see Appendix
J.

Table 2: IHS Study: Fully marginal treatment effect estimation.  Outcome | Method | Estimation | Std.err | P-value | RE
---|---|---|---|---|---
Mood | WCLS | -0.016 | $9.03\times 10^{-3}$ | 0.078 | -
R-WCLS | -0.017 | $8.14\times 10^{-3}$ | 0.038 | 1.23
DR-WCLS | -0.017 | $8.18\times 10^{-3}$ | 0.042 | 1.22
Steps | WCLS | 0.070 | $2.41\times 10^{-2}$ | 0.004 | -
R-WCLS | 0.065 | $2.34\times 10^{-2}$ | 0.005 | 1.06
DR-WCLS | 0.070 | $2.37\times 10^{-2}$ | 0.003 | 1.03

We report various estimators in Table 2 and present more details in Appendix
J. In comparison with WCLS, the estimations using R-WCLS and DR-WCLS have a
tangible improvement in the standard error estimates. We conclude that sending
activity notifications can increase (the cubic root of) step counts by 0.07,
and sending mood notifications can negatively affect users’ moods by -0.017,
with statistical significance at level 95%. By comparison, R-WCLS and DR-WCLS
have enough power to detect a causal relationship between sending mobile
prompts and lower mood scores, while WCLS does not.

###  8.2 Time-varying Treatment Effect Estimation

For further analysis, we include study week in the moderated treatment effect
model: $\beta(t)=\beta_{0}^{\star}+\beta_{1}^{\star}t$, and examine how
treatment effect varies over time. Estimated time-varying treatment moderation
effects and their relative efficiency are shown in Figure 1. The shaded area
in Figure 1 represents the 95% confidence band of the moderation effects as a
function of study week. Narrower confidence bands were observed for estimators
constructed using both R-WCLS and DR-WCLS methods. Relative efficiencies
between 1.2 and 1.3 were observed over the study week.

Figure 1: Causal effects estimates with confidence intervals of R-WCLS (left)
and DR-WCLS (middle), and their relative efficiency in comparisons with WCLS
(right).

Based on the results above, we can conclude that sending notifications doesn’t
have a significant impact on mood scores for the first 12 weeks. Nevertheless,
sending notifications later in the study is less likely to improve the mood of
participants. In light of this, it might not be ideal to overburden
participants for an extended period if the notifications don’t serve any
therapeutic purpose. Additionally, we assessed the time-varying treatment
effect on step count, which is detailed in Appendix J.

###  8.3 Treatment Effect Estimation with Missing Data

Next, we apply our proposed methods to evaluate the treatment effect based on
the raw observed data rather than the imputed dataset. To maintain consistency
with previous analyses, we still use the weekly average mood score and step
count (cubic root) as outcomes. Self-report mood scores and step counts are
collected every day, so if no records were observed for the entire week, we
indicate the weekly outcome as missing. Otherwise, the average mood score and
step count (cubic root) are calculated as outcomes. For mood outcome, there is
in total 31.3% person/week missing, and for step count outcome, there is in
total 48.1% person/week missing.

We carried out the same analysis as above for marginal treatment effects.
Inverse probability weighting is used when implementing estimation using WCLS
and R-WCLS criteria. Estimated treatment effects and their relative efficiency
are shown in Table 3. It is no longer evident that mood notifications have a
significant overall impact on participants’ moods, but the step count analysis
still indicates a positive effect of sending activity notifications on
participants’ physical activity levels.

Table 3: IHS Study: Fully marginal treatment effect estimation with missing outcomes.  Outcome | Method | Estimation | Std.err | P-value
---|---|---|---|---
Mood | WCLS | $7.71\times 10^{-3}$ | $1.73\times 10^{-2}$ | 0.655
R-WCLS | $1.81\times 10^{-3}$ | $1.62\times 10^{-2}$ | 0.911
DR-WCLS | $3.00\times 10^{-3}$ | $1.68\times 10^{-2}$ | 0.858
Steps | WCLS | $6.71\times 10^{-2}$ | $3.94\times 10^{-2}$ | 0.088
R-WCLS | $7.43\times 10^{-2}$ | $4.05\times 10^{-2}$ | 0.067
DR-WCLS | 0.104 | $4.09\times 10^{-2}$ | 0.011

##  9 Discussion

Scientists wish to leverage the large volume of data generated by mobile
health systems to better answer scientific questions regarding the time-
varying intervention effects. The application of ML algorithms can make good
use of high-dimensional mobile health data to provide high-quality predictions
of proximal outcomes. In this paper, we proposed two methods, termed R-WCLS
and DR-WCLS respectively, to estimate time-varying treatment effects and
illustrated both the theoretical and practical benefits of incorporating
machine learning algorithms into the inference procedure. Particularly, both
R-WCLS and DR-WCLS criteria provide sufficient flexibility in their nuisance
model specifications and improve estimation efficiency over existing
approaches. A crucial feature of the DR-WCLS criterion that is not shared by
the WCLS and R-WCLS criteria is that it is doubly robust. This feature is
critical even in MRTs where the treatment model may be misspecified and
missing data are common. The DR-WCLS criterion is especially powerful when
both the treatment randomization probability and conditional expectation model
are accurately specified, which results in the highest relative asymptotic
efficiency.

Although this work represents a major step forward in analyzing MRT data,
there are still some interesting questions to explore in the future. For
example, all the discussion above relies on sample splitting having little to
no impact on the estimation asymptotically, but the effect of the particular
random split on the estimate can be important in finite samples. Chernozhukov
et al. (2018) introduced a finite-sample adjustment to incorporate uncertainty
induced by sample splitting, which can be a straightforward and useful
extension to our proposed methods. In addition, the accuracy of different
supervised learning algorithms in estimating nuisance components might differ,
and it is worth researching whether there are methods that work best in
practice. In relation to the same topic, we can learn the nuisance component
using non-parametric models, such as kernel regression. These methods enable
us to obtain an explicit error bound for the DR-WCLS estimator, which serves
as sufficient conditions for the nuisance model’s smoothness and dimension for
the estimator to be $\sqrt{n}$-consistent (Kennedy, 2020). Last but not least,
Han and Wang (2013) proposed a “multiply robust” estimator that is more robust
than doubly robust estimators, which allows multiple models for both the
propensity score and the outcome regression. This estimator is consistent if
any of the multiple models are correctly specified. It is worth considering
extending our proposed method to a multiply robust version.

## References

  * Athey et al. (2018) Athey, S., G. W. Imbens, and S. Wager (2018, 02).  Approximate Residual Balancing: Debiased Inference of Average Treatment Effects in High Dimensions.  Journal of the Royal Statistical Society Series B: Statistical Methodology 80(4), 597–623.
  * Bang and Robins (2005) Bang, H. and J. M. Robins (2005).  Doubly robust estimation in missing data and causal inference models.  Biometrics 61(4), 962–973.
  * Bickel et al. (1993) Bickel, P. J., C. A. Klaassen, P. J. Bickel, Y. Ritov, J. Klaassen, J. A. Wellner, and Y. Ritov (1993).  Efficient and adaptive estimation for semiparametric models, Volume 4.  Springer.
  * Bodory et al. (2022) Bodory, H., M. Huber, and L. Lafférs (2022).  Evaluating (weighted) dynamic treatment effects by double machine learning.  The Econometrics Journal 25(3), 628–648.
  * Boruvka et al. (2018) Boruvka, A., D. Almirall, K. Witkiewitz, and S. A. Murphy (2018).  Assessing time-varying causal effect moderation in mobile health.  Journal of the American Statistical Association 113(523), 1112–1121.
  * Chen et al. (2022) Chen, Q., V. Syrgkanis, and M. Austern (2022).  Debiased machine learning without sample-splitting for stable estimators.  arXiv preprint arXiv:2206.01825.
  * Chernozhukov et al. (2018) Chernozhukov, V., D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, and J. Robins (2018).  Double/debiased machine learning for treatment and structural parameters.
  * Chernozhukov et al. (2015) Chernozhukov, V., C. Hansen, and M. Spindler (2015).  Post-selection and post-regularization inference in linear models with many controls and instruments.  American Economic Review 105(5), 486–90.
  * Chernozhukov et al. (2021) Chernozhukov, V., W. K. Newey, and R. Singh (2021).  A simple and general debiased machine learning theorem with finite sample guarantees.  arXiv preprint arXiv:2105.15197.
  * Dempsey et al. (2020) Dempsey, W., P. Liao, S. Kumar, and S. A. Murphy (2020).  The stratified micro-randomized trial design: sample size considerations for testing nested causal effects of time-varying treatments.  The annals of applied statistics 14(2), 661.
  * Han and Wang (2013) Han, P. and L. Wang (2013).  Estimation with missing data: beyond double robustness.  Biometrika 100(2), 417–430.
  * Heron and Smyth (2010) Heron, K. E. and J. M. Smyth (2010).  Ecological momentary interventions: incorporating mobile technology into psychosocial and health behaviour treatments.  British journal of health psychology 15(1), 1–39.
  * Hill (2011) Hill, J. L. (2011).  Bayesian nonparametric modeling for causal inference.  Journal of Computational and Graphical Statistics 20(1), 217–240.
  * Kennedy (2016) Kennedy, E. H. (2016).  Semiparametric theory and empirical processes in causal inference.  In Statistical causal inferences and their applications in public health research, pp. 141–167. Springer.
  * Kennedy (2020) Kennedy, E. H. (2020).  Optimal doubly robust estimation of heterogeneous causal effects.  arXiv preprint arXiv:2004.14497.
  * Künzel et al. (2019) Künzel, S. R., J. S. Sekhon, P. J. Bickel, and B. Yu (2019).  Metalearners for estimating heterogeneous treatment effects using machine learning.  Proceedings of the national academy of sciences 116(10), 4156–4165.
  * Lewis and Syrgkanis (2020) Lewis, G. and V. Syrgkanis (2020).  Double/debiased machine learning for dynamic treatment effects via g-estimation.  arXiv preprint arXiv:2002.07285.
  * Liang and Zeger (1986) Liang, K.-Y. and S. L. Zeger (1986).  Longitudinal data analysis using generalized linear models.  Biometrika 73(1), 13–22.
  * Liao et al. (2016) Liao, P., P. Klasjna, A. Tewari, and S. Murphy (2016).  Micro-randomized trials in mhealth.  Statistics in Medicine 35(12), 1944–71.
  * Liu et al. (2018) Liu, Q., L. Li, Z. Tang, and D. Zhou (2018).  Breaking the curse of horizon: Infinite-horizon off-policy estimation.  Advances in Neural Information Processing Systems 31.
  * Morzywolek et al. (2023) Morzywolek, P., J. Decruyenaere, and S. Vansteelandt (2023).  On a general class of orthogonal learners for the estimation of heterogeneous treatment effects.  arXiv preprint arXiv:2303.12687.
  * Murphy et al. (2001) Murphy, S. A., M. J. van der Laan, J. M. Robins, and C. P. P. R. Group (2001).  Marginal mean models for dynamic regimes.  Journal of the American Statistical Association 96(456), 1410–1423.
  * NeCamp et al. (2020) NeCamp, T., S. Sen, E. Frank, M. A. Walton, E. L. Ionides, Y. Fang, A. Tewari, and Z. Wu (2020).  Assessing real-time moderation for developing adaptive mobile health interventions for medical interns: micro-randomized trial.  Journal of medical Internet research 22(3), e15033.
  * Newey (1990) Newey, W. K. (1990).  Semiparametric efficiency bounds.  Journal of applied econometrics 5(2), 99–135.
  * Neyman (1979) Neyman, J. (1979).  C ($\alpha$) tests and their use.  Sankhyā: The Indian Journal of Statistics, Series A, 1–21.
  * Nie and Wager (2021) Nie, X. and S. Wager (2021).  Quasi-oracle estimation of heterogeneous treatment effects.  Biometrika 108(2), 299–319.
  * Qian et al. (2022) Qian, T., A. E. Walton, L. M. Collins, P. Klasnja, S. T. Lanza, I. Nahum-Shani, M. Rabbi, M. A. Russell, M. A. Walton, H. Yoo, et al. (2022).  The microrandomized trial for developing digital interventions: Experimental design and data analysis considerations.  Psychological methods.
  * Qian et al. (2020) Qian, T., H. Yoo, P. Klasnja, D. Almirall, and S. A. Murphy (2020).  Estimating time-varying causal excursion effect in mobile health with binary outcomes.
  * Robins (1986) Robins, J. (1986).  A new approach to causal inference in mortality studies with a sustained exposure period-application to control of the healthy worker survivor effect.  Mathematical Modelling 7(9), 1393–1512.
  * Robins (1994) Robins, J. M. (1994).  Correcting for non-compliance in randomized trials using structural nested mean models.  Communications in Statistics-Theory and methods 23(8), 2379–2412.
  * Robins (1997) Robins, J. M. (1997).  Causal inference from complex longitudinal data.  In Latent variable modeling and applications to causality, pp. 69–117. Springer.
  * Robinson (1988) Robinson, P. M. (1988).  Root-n-consistent semiparametric regression.  Econometrica: Journal of the Econometric Society, 931–954.
  * Rubin (1978) Rubin, D. (1978).  Bayesian inference for causal effects: The role of randomization.  The Annals of Statistics 6(1), 34–58.
  * Semenova and Chernozhukov (2021) Semenova, V. and V. Chernozhukov (2021).  Debiased machine learning of conditional average treatment effects and other causal functions.  The Econometrics Journal 24(2), 264–289.
  * Shi et al. (2022) Shi, J., Z. Wu, and W. Dempsey (2022).  Assessing time-varying causal effect moderation in the presence of cluster-level treatment effect heterogeneity and interference.  Biometrika.
  * Shi et al. (2023) Shi, J., Z. Wu, and W. Dempsey (2023).  Incorporating auxiliary variables to improve the efficiency of time-varying treatment effect estimation.  Manuscript in preparation.
  * Tsiatis (2007) Tsiatis, A. (2007).  Semiparametric Theory and Missing Data.  Springer Science & Business Media.
  * Van Der Laan and Rubin (2006) Van Der Laan, M. J. and D. Rubin (2006).  Targeted maximum likelihood learning.  The international journal of biostatistics 2(1).
  * Van der Vaart (2000) Van der Vaart, A. W. (2000).  Asymptotic statistics, Volume 3.  Cambridge university press.
  * Viviano and Bradic (2021) Viviano, D. and J. Bradic (2021).  Dynamic covariate balancing: estimating treatment effects over time.  arXiv preprint arXiv:2103.01280.

##  Appendix A Proof of Theorem 5.2

Assume $\hat{\beta}^{(R)}_{n}$ minimizes the R-WCLS criterion:

$\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-\left(\tilde{p}_{t}(1|S_{t})g_{t}(1,H_{t})+(1-\tilde{p}_{t}(1|S_{t}))g_{t}(0,H_{t})\right)-(A_{t}-\tilde{p}_{t}(1|S_{t})))f_{t}(S_{t})^{\top}\beta\right)^{2}\right].$
(24)

Denote
$g^{\star}_{t}(H_{t})=\tilde{p}_{t}(1|S_{t})g_{t}(H_{t},1)+(1-\tilde{p}_{t}(1|S_{t}))g_{t}(H_{t},0)=\mathbb{E}[W_{t}Y_{t+1}|H_{t}]$,
where we applied a supervised-learning algorithm and get an estimator
$\hat{g}_{t}(H_{t})$. The asymptotic properties of the R-WCLS estimator follow
from the expansion:

$\displaystyle 0$
$\displaystyle=\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-\hat{g}(H_{t})-(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\hat{\beta}^{(R)}_{n}\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right]$
$\displaystyle=\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-g^{\star}_{t}(H_{t})-(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\beta^{\star}\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right]$
$\displaystyle~{}~{}~{}~{}~{}-\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))^{2}f_{t}(S_{t})f_{t}(S_{t})^{\top}\right](\hat{\beta}^{(R)}_{n}-\beta^{\star})$
$\displaystyle~{}~{}~{}~{}~{}+\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(g^{\star}_{t}(H_{t})-\hat{g}(H_{t}))f_{t}(S_{t})\right]$
(25)

Because

$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))^{2}f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]\overset{P}{\to}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right],$
$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})(g^{\star}(H_{t})-\hat{g}(H_{t}))\right]\overset{P}{\to}0.~{}~{}~{}~{}\text{(by
design)}$

The second property holds true for any $\hat{g}(H_{t})$ because:

$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})(g^{\star}(H_{t})-\hat{g}(H_{t}))\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\Bigg{[}\sum_{t=1}^{T}\mathbb{E}[\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})(g^{\star}(H_{t})-\hat{g}(H_{t}))|A_{t}=1]$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\mathbb{E}[(1-\tilde{p}_{t}(1|S_{t}))(0-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})(g^{\star}(H_{t})-\hat{g}(H_{t}))|A_{t}=0]\Bigg{]}$
$\displaystyle=$ $\displaystyle 0.$

By solving (A), we obtain:

$\displaystyle n^{1/2}(\hat{\beta}^{(R)}_{n}-\beta^{\star})$
$\displaystyle=\left\\{\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))^{2}f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]\right\\}^{-1}\times$
$\displaystyle~{}~{}~{}~{}~{}~{}n^{1/2}\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-g^{\star}_{t}(H_{t})-(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\beta^{\star}\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right].$

Denote $\tilde{Y}^{(R)}_{t+1}=Y_{t+1}-g^{\star}_{t}(H_{t})$, and we obtain

$\displaystyle n^{1/2}(\hat{\beta}^{(R)}_{n}-\beta^{\star})$
$\displaystyle=n^{1/2}~{}\mathbb{P}_{n}\Bigg{[}\sum_{t=1}^{T}~{}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]^{-1}\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\beta^{\star}\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\Bigg{]}+o_{p}(1).$

By definition of $\beta^{\star}$:

$\mathbb{E}\left[\sum_{t=1}^{T}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\beta^{\star}\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right]=0$

Consequently, under regularity conditions, the estimator
$\hat{\beta}^{(R)}_{n}\overset{P}{\to}\beta^{\star}$; that is,
$\hat{\beta}^{(R)}_{n}$ is a consistent estimator of $\beta^{\star}$. The
influence function for $\hat{\beta}^{(R)}_{n}$ is:

$\displaystyle\sum_{t=1}^{T}~{}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]^{-1}\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\beta^{\star}\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t}).$
(26)

Then under moment conditions, we have asymptotic normality with variance given
by $\Sigma_{R}=Q^{-1}WQ^{-1}$, where

$\displaystyle Q$
$\displaystyle=\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right],$
$\displaystyle W$
$\displaystyle=\mathbb{E}\left[\left(\sum_{t=1}^{T}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\beta^{\star}\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right)^{2}\right],$

due to space constraints, we use $\mathbb{E}[X^{2}]$ to denote
$\mathbb{E}[XX^{\top}]$. In conclusion, we establish that the estimator
minimizing the R-WCLS criterion $\hat{\beta}^{(R)}_{n}$ is consistent and
asymptotically normal:

$n^{1/2}(\hat{\beta}^{(R)}_{n}-\beta^{\star})\sim\mathcal{N}(0,\Sigma_{R}).$

We further prove that the variance is consistently estimated by Equation (11).
Using sample splitting, the estimating equation can be written as:

$\displaystyle 0$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-\hat{g}^{(k)}(H_{t})-(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\hat{\beta}^{(R)}_{n}\right)(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))f_{t}(S_{t})\right]$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-g^{\star}_{t}(H_{t})-(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\beta^{\star}\right)(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))f_{t}(S_{t})\right]$
$\displaystyle~{}~{}~{}~{}~{}-\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))^{2}f_{t}(S_{t})f_{t}(S_{t})^{\top}\right](\hat{\beta}^{(R)}_{n}-\beta^{\star})$
$\displaystyle~{}~{}~{}~{}~{}+\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(g^{\star}_{t}(H_{t})-\hat{g}^{(k)}(H_{t}))f_{t}(S_{t})\right]$
(27)

Assume $K$ is finite and fixed, and we have the same reasoning as above:

$\displaystyle\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))^{2}f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]\overset{P}{\to}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}^{(k)}_{t}(1|S_{t})(1-\tilde{p}^{(k)}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right],$
$\displaystyle\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))f_{t}(S_{t})(g^{\star}(H_{t})-\hat{g}^{(k)}(H_{t}))\right]\overset{P}{\to}0.~{}~{}~{}~{}\text{(by
design)}$

Then we obtain the following:

$\displaystyle n^{1/2}(\hat{\beta}^{(R)}_{n}-\beta^{\star})$
$\displaystyle=\left\\{\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))^{2}f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]\right\\}^{-1}\times$
$\displaystyle
n^{1/2}\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-g^{\star}_{t}(H_{t})-(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\beta^{\star}\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right]$

Further,

$\displaystyle n^{1/2}(\hat{\beta}^{(R)}_{n}-\beta^{\star})$
$\displaystyle=n^{1/2}~{}\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\Bigg{[}\sum_{t=1}^{T}~{}\left\\{\frac{1}{K}\sum_{k=1}^{K}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}^{(k)}_{t}(1|S_{t})(1-\tilde{p}^{(k)}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]\right\\}^{-1}\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\beta^{\star}\right)(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))f_{t}(S_{t})\Bigg{]}+o_{p}(1).$

By definition of $\beta^{\star}$:

$\mathbb{E}\left[\sum_{t=1}^{T}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\beta^{\star}\right)(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))f_{t}(S_{t})\right]=0$

Consequently, under regularity conditions, the estimator
$\hat{\beta}^{(R)}_{n}\overset{P}{\to}\beta^{\star}$; that is,
$\hat{\beta}^{(R)}_{n}$ is consistent. The influence function for
$\hat{\beta}^{(R)}_{n}$ is:

$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\sum_{t=1}^{T}~{}\left\\{\frac{1}{K}\sum_{k=1}^{K}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}^{(k)}_{t}(1|S_{t})(1-\tilde{p}^{(k)}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]\right\\}^{-1}\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))f_{t}(S_{t})^{\top}\beta^{\star}\right)(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))f_{t}(S_{t}).$
(28)

Recall

$\displaystyle m$
$\displaystyle=\sum_{t=1}^{T}\hat{W}^{(k)}_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\hat{\tilde{p}}_{t}^{(k)}(1|S_{t}))f_{t}(S_{t})^{\top}\hat{\beta}^{(R)}_{n}\right)(A_{t}-\hat{\tilde{p}}_{t}^{(k)}(1|S_{t}))f_{t}(S_{t}),$
$\displaystyle\dot{m}$ $\displaystyle=\frac{\partial
m(\beta,\eta)}{\partial\beta}=\sum_{t=1}^{T}\hat{\tilde{p}}_{t}^{(k)}(1|S_{t})(1-\hat{\tilde{p}}_{t}^{(k)}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}.$

Then the variance can be consistently estimated by:

$\left[\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\dot{m}(\hat{\beta},\hat{\eta}_{k})\right]^{-1}\times\left[\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}m(\hat{\beta},\hat{\eta}_{k})m(\hat{\beta},\hat{\eta}_{k})^{\top}\right]\times\left[\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\dot{m}(\hat{\beta},\hat{\eta}_{k})\right]^{-1}.$

##  Appendix B Proof of Theorem 5.3

###  B.1 Double robustness property

The following is proof of the double robustness of the DR-WCLS estimator.
Assume $\hat{\beta}^{(DR)}_{n}$ minimizes the _DR-WCLS_ criterion:

$\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\left(\frac{W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g_{t}(H_{t},A_{t}))}{\tilde{\sigma}^{2}_{t}(S_{t})}+\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta\right)^{2}\right].$

Here the true randomization probability is $p_{t}(A_{t}|H_{t})$, and the
conditional expectation (also known as the outcome regression):

$\mathbb{E}[Y_{t+1}|H_{t},A_{t}]=g_{t}^{\star}(H_{t},A_{t}).$

Denote $\beta(t;H_{t})=g_{t}^{\star}(H_{t},1)-g_{t}^{\star}(H_{t},0)$. The
corresponding ML estimators are denoted as $\hat{g}_{t}(H_{t},A_{t})$ and
$\hat{\beta}(t;H_{t})$. We consider the estimating equation of the objective
function above:

$\displaystyle 0$
$\displaystyle=\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\left(\frac{W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g_{t}(H_{t},A_{t}))}{\tilde{\sigma}^{2}_{t}(S_{t})}+\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta^{(DR)}\right)f_{t}(S_{t})\right]$
$\displaystyle=\mathbb{E}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g_{t}(H_{t},A_{t}))f_{t}(S_{t})\right]+\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\left(\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta^{(DR)}\right)f_{t}(S_{t})\right].$

If the conditional expectation $g_{t}^{\star}(H_{t},A_{t})$ is correctly
specified, we have:

$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g_{t}^{\star}(H_{t},A_{t}))f_{t}(S_{t})\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))\mathbb{E}[Y_{t+1}-g_{t}^{\star}(H_{t},A_{t})|H_{t},A_{t}]f_{t}(S_{t})\right]$
$\displaystyle=$ $\displaystyle 0,$

and we are left to solve:

$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\left(\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta^{(DR)}\right)f_{t}(S_{t})\right]=0.$

Under regularity conditions, the estimator
$\hat{\beta}^{(DR)}_{n}\overset{P}{\to}\beta^{\star}$; that is,
$\hat{\beta}^{(DR)}_{n}$ is consistent. Another case is when the treatment
randomization probability is correctly specified. Then we have:

$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-\hat{g}_{t}(H_{t},A_{t}))f_{t}(S_{t})\right]+\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\left(\hat{\beta}(t;H_{t})-f_{t}(S_{t})^{\top}\beta^{(DR)}\right)f_{t}(S_{t})\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))\Big{(}\mathbb{E}[Y_{t+1}|H_{t},A_{t}=1]-\mathbb{E}[Y_{t+1}|H_{t},A_{t}=0]-\hat{\beta}(H_{t},A_{t})\Big{)}f_{t}(S_{t})\right]$
$\displaystyle+\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\left(\hat{\beta}(t;H_{t})-f_{t}(S_{t})^{\top}\beta^{(DR)}\right)f_{t}(S_{t})\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\left(\mathbb{E}[Y_{t+1}|H_{t},A_{t}=1]-\mathbb{E}[Y_{t+1}|H_{t},A_{t}=0]-f_{t}(S_{t})^{\top}\beta^{(DR)}\right)f_{t}(S_{t})\right]$

Under regularity conditions, the estimator
$\hat{\beta}^{(DR)}_{n}\overset{P}{\to}\beta^{\star}$; that is,
$\hat{\beta}^{(DR)}_{n}$ is consistent.

###  B.2 Asymptotic properties for DR-WCLS estimators

Assume $\hat{\beta}^{(DR)}_{n}$ minimizes the _DR-WCLS_ criterion:

$\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\left(\frac{\hat{W}_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-\hat{g}_{t}(H_{t},A_{t}))}{\tilde{\sigma}^{2}_{t}(S_{t})}+\hat{\beta}(t;H_{t})-f_{t}(S_{t})^{\top}\beta^{(DR)}\right)^{2}\right].$

The estimated treatment randomization probability is denoted as
$\hat{p}_{t}=\hat{p}(A_{t}|H_{t})$, thus we have the weight $W_{t}$ estimated
by $\hat{W}_{t}=\tilde{p}_{t}(A_{t}|S_{t})/\hat{p}_{t}(A_{t}|H_{t})$. And the
estimating equation is:

$\displaystyle 0$
$\displaystyle=\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\left(\frac{\tilde{p}_{t}(A_{t}|S_{t})(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-\hat{g}_{t}(H_{t},A_{t}))}{\hat{p}_{t}(A_{t}|H_{t})\tilde{\sigma}^{2}_{t}(S_{t})}+\hat{\beta}(t;H_{t})-f_{t}(S_{t})^{\top}\beta^{(DR)}\right)f_{t}(S_{t})\right].$
(29)

Expand the right-hand side, we have:

$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\left(\frac{\tilde{p}_{t}(A_{t}|S_{t})(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-\hat{g}_{t}(H_{t},A_{t}))}{\hat{p}_{t}(A_{t}|H_{t})\tilde{\sigma}^{2}_{t}(S_{t})}+\hat{\beta}(t;H_{t})-f_{t}(S_{t})^{\top}\beta^{(DR)}\right)f_{t}(S_{t})\right]$
$\displaystyle=$
$\displaystyle\mathbb{P}_{n}\Bigg{[}\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\Big{(}\frac{\tilde{p}_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g^{\star}_{t}(H_{t},A_{t})+g^{\star}_{t}(H_{t},A_{t})-\hat{g}_{t}(H_{t},A_{t}))}{\tilde{\sigma}^{2}_{t}(S_{t})}\left(\frac{1}{\hat{p}_{t}}-\frac{1}{p_{t}}+\frac{1}{p_{t}}\right)+$
$\displaystyle\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta^{\star}+(\hat{\beta}(t;H_{t})-\beta(t;H_{t}))-f_{t}(S_{t})^{\top}(\beta^{(DR)}-\beta^{\star})\Big{)}f_{t}(S_{t})\Bigg{]}$
$\displaystyle=$
$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\left(\frac{W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g^{\star}_{t}(H_{t},A_{t}))}{\tilde{\sigma}^{2}_{t}(S_{t})}+\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta^{\star}\right)f_{t}(S_{t})\right]+$
$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{p}_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g^{\star}_{t}(H_{t},A_{t}))\left(\frac{1}{\hat{p}_{t}}-\frac{1}{p_{t}}\right)f_{t}(S_{t})\right]+$
$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{p}_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(g^{\star}_{t}(H_{t},A_{t})-\hat{g}_{t}(H_{t},A_{t}))\left(\frac{1}{\hat{p}_{t}}-\frac{1}{p_{t}}\right)f_{t}(S_{t})\right]+$
$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(g^{\star}(H_{t},A_{t})-\hat{g}_{t}(H_{t},A_{t}))f_{t}(S_{t})\right]+\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})(\hat{\beta}(t;H_{t})-\beta(t;H_{t}))f_{t}(S_{t})\right]$
$\displaystyle-\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})f_{t}(S_{t})f^{\top}_{t}(S_{t})\right](\hat{\beta}^{(DR)}_{n}-\beta^{\star})$

Because

$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{p}_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g^{\star}_{t}(H_{t},A_{t}))\left(\frac{1}{\hat{p}_{t}}-\frac{1}{p_{t}}\right)f_{t}(S_{t})\right]$
$\displaystyle\overset{P}{\to}0,~{}~{}~{}~{}\text{(correct model
specification)}$
$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})f_{t}(S_{t})f^{\top}_{t}(S_{t})\right]$
$\displaystyle\overset{P}{\to}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})f_{t}(S_{t})f^{\top}_{t}(S_{t})\right],$

and

$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(g^{\star}(H_{t},A_{t})-\hat{g}_{t}(H_{t},A_{t}))f_{t}(S_{t})\right]+$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})(\hat{\beta}(t;H_{t})-\beta(t;H_{t}))f_{t}(S_{t})\right]\overset{P}{\to}0.~{}~{}~{}~{}\text{(terms
cancellation)}$

Apart from the nicely-behaved term above, the only term that might be
problematic is:

$\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{p}_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(g^{\star}_{t}(H_{t},A_{t})-\hat{g}_{t}(H_{t},A_{t}))\left(\frac{1}{\hat{p}_{t}}-\frac{1}{p_{t}}\right)f_{t}(S_{t})\right],$

which will converge to:

$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}\left(\underbrace{\sum_{a\in\\{0,1\\}}\frac{\tilde{\sigma}^{2}_{t}(S_{t})}{a\hat{p}_{t}(1|H_{t})+(1-a)(1-\hat{p}_{t}(1|H_{t}))}(p_{t}(1|H_{t})-\hat{p}_{t}(1|H_{t}))(g(H_{t},a)-\hat{g}(H_{t},a))}_{\text{(I)}}\right)f_{t}(S_{t})\right].$

In our context, $T$ is finite and fixed. Therefore, by the fact that
$\hat{p}_{t}(1|H_{t}))$ is bounded away from zero and one, along with the
Cauchy–Schwarz inequality, we have that (up to a multiplicative constant) term
(I) is bounded above by:

$\mathbf{\hat{B}}=\mathbb{E}\left[\sum_{t=1}^{T}\sum_{a\in\\{0,1\\}}\left\|p_{t}(1|H_{t})-\hat{p}_{t}(1|H_{t})\right\|\left\|g(H_{t},a)-\hat{g}(H_{t},a)\right\|\right].$
(30)

Assuming we have nuisance estimates that can make $\mathbf{\hat{B}}$
asymptotically negligible, then the DR-WCLS estimator satisfies:

$\displaystyle n^{1/2}(\hat{\beta}^{(DR)}_{n}-\beta^{\star})$
$\displaystyle=n^{1/2}~{}\mathbb{P}_{n}\Bigg{[}\sum_{t=1}^{T}~{}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})f_{t}(S_{t})f^{\top}_{t}(S_{t})\right]^{-1}\tilde{\sigma}^{2}_{t}(S_{t})(\beta(t;H_{t})-f_{t}(S_{t})\beta^{\star})f_{t}(S_{t})\Bigg{]}+o_{p}(1),$

and it is efficient with influence function:

$\displaystyle\sum_{t=1}^{T}~{}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]^{-1}\tilde{\sigma}^{2}_{t}(S_{t})(\beta(t;H_{t})-f_{t}(S_{t})\beta^{\star})f_{t}(S_{t}).$

Under moment conditions, we have asymptotic normality with variance given by
$\Sigma_{DR}=Q^{-1}WQ^{-1}$, where

$\displaystyle Q$
$\displaystyle=\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})f_{t}(S_{t})f_{t}(S_{t})^{\top}\right],$
$\displaystyle W$
$\displaystyle=\mathbb{E}\left[\left(\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})(\beta(t;H_{t})-f_{t}(S_{t})\beta^{\star})f_{t}(S_{t})\right)^{2}\right].$

###  B.3 Asymptotic variance using sample splitting

Built on the previous doubly robust property, we know that if either the
conditional expectation model $g_{t}(H_{t},A_{t})$ or the treatment
randomization probability $p_{t}(A_{t}|H_{t})$ is correctly specified, we can
obtain a consistent estimator of $\beta^{\star}$. In this section, we provide
the asymptotic variance estimation under sample splitting. Without loss of
generality, we assume that the treatment randomization probability
$p_{t}(A_{t}|H_{t})$ is correctly specified. For simplicity, we use
$\tilde{\sigma}_{t}^{2(k)}$ to denote $\tilde{\sigma}^{2}_{t}(S_{t})^{(k)}$
The asymptotic properties of the DR-WCLS estimator follow from the expansion:

$\displaystyle 0$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}\tilde{\sigma}_{t}^{2(k)}\left(\frac{W_{t}(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))(Y_{t+1}-\hat{g}^{(k)}_{t}(H_{t},A_{t}))}{\tilde{\sigma}_{t}^{2(k)}}+\hat{\beta}^{(k)}(t;H_{t})-f_{t}(S_{t})^{\top}\hat{\beta}^{(DR)}_{n}\right)f_{t}(S_{t})\right]$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\Bigg{[}\sum_{t=1}^{T}\tilde{\sigma}_{t}^{2(k)}\Big{(}\frac{W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g^{\star}(H_{t},A_{t})+g^{\star}(H_{t},A_{t})-\hat{g}^{(k)}_{t}(H_{t},A_{t}))}{\tilde{\sigma}_{t}^{2(k)}}+$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}\beta(t;H_{t})+(\hat{\beta}^{(k)}(t;H_{t})-\beta(t;H_{t}))-f_{t}(S_{t})^{\top}\beta^{\star}-f_{t}(S_{t})^{\top}(\hat{\beta}^{(DR)}_{n}-\beta^{\star})\Big{)}f_{t}(S_{t})\Bigg{]}$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}\tilde{\sigma}_{t}^{2(k)}\left(\frac{W_{t}(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))(Y_{t+1}-g^{\star}(H_{t},A_{t}))}{\tilde{\sigma}_{t}^{2(k)}}+\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta^{\star}\right)f_{t}(S_{t})\right]$
$\displaystyle+\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))(g^{\star}(H_{t},A_{t})-\hat{g}^{(k)}_{t}(H_{t},A_{t}))f_{t}(S_{t})\right]$
$\displaystyle+\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}\tilde{\sigma}_{t}^{2(k)}(\hat{\beta}^{(k)}(t;H_{t})-\beta(t;H_{t}))f_{t}(S_{t})\right]$
$\displaystyle-\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}\tilde{\sigma}_{t}^{2(k)}f_{t}(S_{t})f^{\top}_{t}(S_{t})\right](\hat{\beta}^{(DR)}_{n}-\beta^{\star})$

Because

$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))(g^{\star}(H_{t},A_{t})-\hat{g}^{(k)}_{t}(H_{t},A_{t}))f_{t}(S_{t})\right]+$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}\tilde{\sigma}_{t}^{2(k)}(\hat{\beta}^{(k)}(t;H_{t})-\beta(t;H_{t}))f_{t}(S_{t})\right]\overset{P}{\to}0,$

and

$\displaystyle\mathbb{P}_{n,k}\left[\sum_{t=1}^{T}\tilde{\sigma}_{t}^{2(k)}f_{t}(S_{t})f^{\top}_{t}(S_{t})\right]$
$\displaystyle\overset{P}{\to}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}_{t}^{2(k)}f_{t}(S_{t})f^{\top}_{t}(S_{t})\right].$

Denote

$\tilde{Y}^{(DR)}_{t+1}=\frac{W_{t}(A_{t}-\tilde{p}^{(k)}_{t}(1|S_{t}))(Y_{t+1}-g^{\star}(H_{t},A_{t}))}{\tilde{\sigma}_{t}^{2(k)}}+\beta(t;H_{t}),$

we obtain:

$\displaystyle n^{1/2}(\hat{\beta}^{(DR)}_{n}-\beta^{\star})$
$\displaystyle=n^{1/2}~{}\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}_{n,k}\Bigg{[}\sum_{t=1}^{T}~{}\left\\{\frac{1}{K}\sum_{k=1}^{K}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})f_{t}(S_{t})f^{\top}_{t}(S_{t})\right]\right\\}^{-1}\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\tilde{\sigma}_{t}^{2(k)}(\tilde{Y}^{(DR)}_{t+1}-f_{t}(S_{t})\beta^{\star})f_{t}(S_{t})\Bigg{]}+o_{p}(1).$

By definition of $\beta^{\star}$:

$\frac{1}{K}\sum_{k=1}^{K}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}_{t}^{2(k)}\left(\tilde{Y}^{(DR)}_{t+1}-f_{t}(S_{t})^{\top}\beta^{\star}\right)f_{t}(S_{t})\right]=0.$

Consequently, the influence function for $\hat{\beta}^{(DR)}_{n}$ is:

$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\sum_{t=1}^{T}~{}\left\\{\frac{1}{K}\sum_{k=1}^{K}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}_{t}^{2(k)}f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]\right\\}^{-1}\tilde{\sigma}_{t}^{2(k)}(\tilde{Y}^{(DR)}_{t+1}-f_{t}(S_{t})\beta^{\star})f_{t}(S_{t}).$
(31)

Then under moment conditions, we have asymptotic normality with variance given
by Equation (11) where:

$\displaystyle m$
$\displaystyle=\sum_{t=1}^{T}\hat{\tilde{p}}_{t}^{(k)}(1|S_{t})(1-\hat{\tilde{p}}_{t}^{(k)}(1|S_{t}))(\tilde{Y}^{(DR)}_{t+1}-f_{t}(S_{t})\hat{\beta}^{(DR)}_{n})f_{t}(S_{t}),$
$\displaystyle\dot{m}$ $\displaystyle=\frac{\partial
m(\beta,\eta)}{\partial\beta}=\sum_{t=1}^{T}\hat{\tilde{p}}_{t}^{(k)}(1|S_{t})(1-\hat{\tilde{p}}_{t}^{(k)}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}.$

In conclusion, we establish that the estimator minimizing the DR-WCLS
criterion $\hat{\beta}^{(DR)}_{n}$ is consistent and asymptotically normal.

##  Appendix C Sample Splitting

The estimation technique developed in this paper relies on K-fold cross-
validation obtained by randomly partitioning the sample, i.e., estimate the
nuisance models $\hat{g}_{t}(H_{t}),\hat{g}_{t}(H_{t},A_{t})$ and
$\hat{p}_{t}(t;H_{t})$ on one part of the data (training data) and estimate
the parameter of interest $\hat{\beta}$ on the other part of the data (test
data).

Cross-fitting plays an important role here: the regression procedure as
defined estimates the pseudo-outcome on a separate sample, independent from
the one used in the second-stage regression (Kennedy, 2020), which allows for
informative error analysis while being agnostic about the first- and second-
stage methods.

##  Appendix D Proof of Lemma 5.5

###  D.1 Implementation of the augmented R-WCLS criterion

Let $f_{t}(S_{t})^{\perp}$ denote the orthogonal complement of $f_{t}(S_{t})$
in $H_{t}$, which refers to the set of random variables that are uncorrelated
with $f_{t}(S_{t})$ smoothing over time. Here, a rigorous definition of the
orthogonal complement of $f_{t}(S_{t})$ is given below (Shi et al., 2023):

$\displaystyle f_{t}(S_{t})^{\perp}\coloneqq\left\\{X_{t}\in
H_{t}:\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1-\tilde{p}_{t})X_{t}f_{t}(S_{t})\right]=0\right\\}.$
(32)

To construct $\Delta_{t}^{\perp}$, i.e., the projection of $\beta(t;H_{t})$
onto $f_{t}(S_{t})^{\perp}$, we can apply a linear working model as follows:

$\beta(t;H_{t})\sim(f_{t}(S_{t})^{\perp})^{\top}\eta+f_{t}(S_{t})^{\top}\beta.$

Therefore, $\Delta_{t}^{\perp}=(f_{t}(S_{t})^{\perp})^{\top}\eta$. This
approach allows us to effectively leverage the information from the nuisance
functions $\beta(t;H_{t})=g_{t}(H_{t},1)-g_{t}(H_{t},0)$, which can be
decomposed into $\Delta_{t}^{\perp}$ and $f_{t}(S_{t})^{\top}\beta$. Most
importantly, the inclusion of $\Delta_{t}^{\perp}$ in the estimating equation
does not compromise the consistency of the estimator $\hat{\beta}_{n}^{(AR)}$.

###  D.2 Asymptotic properties

Assume $\hat{\beta}^{(AR)}_{n}$ minimizes the augmented R-WCLS criterion and
therefore it solves:

$0=\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-g_{t}(H_{t})-\left(A_{t}-\tilde{p}_{t}(1|S_{t})\right)\left(f_{t}(S_{t})^{\top}\beta+\Delta_{t}^{\perp}\right)\right)\left(A_{t}-\tilde{p}_{t}(1|S_{t})\right)f_{t}(S_{t})\right].$

The asymptotic properties of the augmented R-WCLS estimator follow from the
expansion:

$\displaystyle 0$
$\displaystyle=\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-\hat{g}(H_{t})-(A_{t}-\tilde{p}_{t}(1|S_{t}))\left(f_{t}(S_{t})^{\top}\hat{\beta}^{(AR)}_{n}+\hat{\Delta}_{t}^{\perp}\right)\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right]$
$\displaystyle=\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-g^{\star}_{t}(H_{t})-(A_{t}-\tilde{p}_{t}(1|S_{t}))\left(f_{t}(S_{t})^{\top}\beta^{\star}+\Delta_{t}^{\perp}\right)\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right]$
$\displaystyle~{}~{}~{}~{}~{}+\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(g^{\star}_{t}(H_{t})-\hat{g}(H_{t}))f_{t}(S_{t})\right]$
$\displaystyle~{}~{}~{}~{}~{}-\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))^{2}(\hat{\Delta}_{t}^{\perp}-\Delta_{t}^{\perp})f_{t}(S_{t})\right]$
$\displaystyle~{}~{}~{}~{}~{}-\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))^{2}f_{t}(S_{t})f_{t}(S_{t})^{\top}\right](\hat{\beta}^{(AR)}_{n}-\beta^{\star})$

Because

$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))^{2}f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]\overset{P}{\to}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right],$
$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})(g^{\star}(H_{t})-\hat{g}(H_{t}))\right]\overset{P}{\to}0,~{}~{}~{}~{}\text{(by
design)}$
$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))^{2}(\hat{\Delta}_{t}^{\perp}-\Delta_{t}^{\perp})f_{t}(S_{t})\right]\overset{P}{\to}0.~{}~{}~{}~{}\text{(orthogonal
projection)}$

We obtain:

$\displaystyle n^{1/2}(\hat{\beta}^{(AR)}_{n}-\beta^{\star})$
$\displaystyle=\left\\{\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))^{2}f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]\right\\}^{-1}\times$
$\displaystyle
n^{1/2}\mathbb{P}_{n}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-g^{\star}_{t}(H_{t})-(A_{t}-\tilde{p}_{t}(1|S_{t}))\left(f_{t}(S_{t})^{\top}\beta^{\star}+\Delta_{t}^{\perp}\right)\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right]$

Recall $\tilde{Y}^{(R)}_{t+1}=Y_{t+1}-g^{\star}_{t}(H_{t})$, and we obtain

$\displaystyle n^{1/2}(\hat{\beta}^{(AR)}_{n}-\beta^{\star})$
$\displaystyle=n^{1/2}~{}\mathbb{P}_{n}\Bigg{[}\sum_{t=1}^{T}~{}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]^{-1}\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}_{t}(1|S_{t}))\left(f_{t}(S_{t})^{\top}\beta^{\star}+\Delta_{t}^{\perp}\right)\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\Bigg{]}+o_{p}(1).$

By definition of $\beta^{\star}$:

$\mathbb{E}\left[\sum_{t=1}^{T}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}_{t}(1|S_{t}))\left(f_{t}(S_{t})^{\top}\beta^{\star}+\Delta_{t}^{\perp}\right)\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right]=0.$

Consequently, under regularity conditions, the estimator
$\hat{\beta}^{(AR)}_{n}\overset{P}{\to}\beta^{\star}$; that is,
$\hat{\beta}^{(AR)}_{n}$ is consistent. The influence function for
$\hat{\beta}^{(AR)}_{n}$ is:

$\displaystyle\sum_{t=1}^{T}~{}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]^{-1}\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}_{t}(1|S_{t}))\left(f_{t}(S_{t})^{\top}\beta^{\star}+\Delta_{t}^{\perp}\right)\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t}).$
(33)

Then under moment conditions, we have asymptotic normality with variance given
by $\Sigma_{R}=Q^{-1}WQ^{-1}$, where

$\displaystyle Q$
$\displaystyle=\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right],$
$\displaystyle W$
$\displaystyle=\mathbb{E}\left[\left(\sum_{t=1}^{T}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}_{t}(1|S_{t}))\left(f_{t}(S_{t})^{\top}\beta^{\star}+\Delta_{t}^{\perp}\right)\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right)^{2}\right].$

In conclusion, we establish that the estimator minimizing the augmented R-WCLS
criterion $\hat{\beta}^{(AR)}_{n}$ is consistent and asymptotically normal.
Under sample splitting, the asymptotic variance can be estimated by Equation
(11) with:

$\displaystyle m$
$\displaystyle=\sum_{t=1}^{T}\hat{W}^{(k)}_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\hat{\tilde{p}}_{t}^{(k)}(1|S_{t}))\left(f_{t}(S_{t})^{\top}\hat{\beta}^{(AR)}_{n}+\hat{\Delta}_{t}^{\perp(k)}\right)\right)(A_{t}-\hat{\tilde{p}}_{t}^{(k)}(1|S_{t}))f_{t}(S_{t}),$
$\displaystyle\dot{m}$ $\displaystyle=\frac{\partial
m(\beta,\eta)}{\partial\beta}=\sum_{t=1}^{T}\hat{\tilde{p}}_{t}^{(k)}(1|S_{t})(1-\hat{\tilde{p}}_{t}^{(k)}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}.$

##  Appendix E Proof of Theorem 5.7

###  E.1 Efficiency gain over WCLS estimator

To reconcile the notations, we write the estimating equation in a general
form, from which can obtain a consistent estimate of $\beta^{\star}$ by
solving:

$\mathbb{E}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-g_{t}(H_{t})-\left(A_{t}-\tilde{p}_{t}(1|S_{t})\right)\left(f_{t}(S_{t})^{\top}\beta+\Delta_{t}^{\perp}\right)\right)\left(A_{t}-\tilde{p}_{t}(1|S_{t})\right)f_{t}(S_{t})\right]=0.$

For WCLS, denote the linear working model for
$\mathbb{E}[Y_{t+1}|H_{t},A_{t}]$ as $\tilde{g}_{t}(H_{t},A_{t})$. We can then
write the estimating equation as:

$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-\tilde{g}_{t}(H_{t})-\left(A_{t}-\tilde{p}_{t}(1|S_{t})\right)(f_{t}(S_{t})^{\top}\beta+\tilde{\Delta}_{t}^{\perp})\right)\left(A_{t}-\tilde{p}_{t}(1|S_{t})\right)f_{t}(S_{t})\right]=0.$

For an augmented R-WCLS estimator, recall
$\tilde{Y}^{(R)}_{t+1}=Y_{t+1}-g(H_{t})$, the estimating equation can be
written as:

$\displaystyle 0$
$\displaystyle=\mathbb{E}\left[\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-g_{t}(H_{t})-\left(A_{t}-\tilde{p}_{t}(1|S_{t})\right)\left(f_{t}(S_{t})^{\top}\beta+\Delta_{t}^{\perp}\right)\right)\left(A_{t}-\tilde{p}_{t}(1|S_{t})\right)f_{t}(S_{t})\right]$
$\displaystyle=\mathbb{E}\left[\sum_{t=1}^{T}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-\left(A_{t}-\tilde{p}_{t}(1|S_{t})\right)\left(f_{t}(S_{t})^{\top}\beta+\Delta_{t}^{\perp}\right)\right)\left(A_{t}-\tilde{p}_{t}(1|S_{t})\right)f_{t}(S_{t})\right]$

Since both methods can produce consistent estimations, now we consider the
comparison between their asymptotic variances. For WCLS, the asymptotic
variance can be calculated as:

$\displaystyle\Sigma=$
$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]^{-1}\times$
$\displaystyle~{}~{}~{}~{}~{}\mathbb{E}\left[\left(\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-\tilde{g}_{t}(H_{t})-(A_{t}-\tilde{p}_{t}(1|S_{t}))(f_{t}(S_{t})^{\top}\beta^{\star}+\tilde{\Delta}_{t}^{\perp})\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right)^{2}\right]\times$
$\displaystyle~{}~{}~{}~{}~{}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]^{-1}$

and for the augmented R-WCLS, the asymptotic variance can be calculated as:

$\displaystyle\Sigma^{(AR)}=$
$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]^{-1}\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathbb{E}\left[\left(\sum_{t=1}^{T}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}_{t}(1|S_{t}))(f_{t}(S_{t})^{\top}\beta^{\star}+\Delta_{t}^{\perp})\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right)^{2}\right]\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]^{-1}.$

Denote $\epsilon(H_{t})=g(H_{t})-\tilde{g}(H_{t})$, we have the following
derivation:

$\displaystyle\mathbb{E}\left[\left(\sum_{t=1}^{T}W_{t}\left(Y_{t+1}-\tilde{g}_{t}(H_{t})-(A_{t}-\tilde{p}_{t}(1|S_{t}))(f_{t}(S_{t})^{\top}\beta^{\star}+\tilde{\Delta}_{t}^{\perp})\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right)^{2}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\left(\sum_{t=1}^{T}W_{t}\left(\tilde{Y}^{(R)}_{t+1}+\epsilon(H_{t})-(A_{t}-\tilde{p}_{t}(1|S_{t}))(f_{t}(S_{t})^{\top}\beta^{\star}+\tilde{\Delta}_{t}^{\perp}+\Delta_{t}^{\perp}-\Delta_{t}^{\perp})\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right)^{2}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\left(\sum_{t=1}^{T}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}_{t}(1|S_{t}))(f_{t}(S_{t})^{\top}\beta^{\star}+\Delta_{t}^{\perp})\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right)^{2}\right]$
$\displaystyle~{}~{}~{}~{}~{}~{}+\mathbb{E}\Bigg{[}\left(\sum_{t=1}^{T}W_{t}\epsilon(H_{t},A_{t})(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right)^{2}\Bigg{]}$
$\displaystyle\geq$
$\displaystyle\mathbb{E}\left[\left(\sum_{t=1}^{T}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}_{t}(1|S_{t}))(f_{t}(S_{t})^{\top}\beta^{\star}+\Delta_{t}^{\perp})\right)(A_{t}-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})\right)^{2}\right]$

where
$\epsilon(H_{t},A_{t})=\epsilon(H_{t})+(A_{t}-\tilde{p}_{t}(1|S_{t}))(\Delta_{t}^{\perp}-\tilde{\Delta}_{t}^{\perp})$.
The interaction term is 0 because:

$\displaystyle\mathbb{E}\left[\sum_{t,t^{\prime}}^{T}W_{t}\left(\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}_{t})(f_{t}(S_{t})^{\top}\beta^{\star}+\Delta_{t}^{\perp})\right)(A_{t}-\tilde{p}_{t})f_{t}(S_{t})W_{t^{\prime}}\epsilon(H_{t^{\prime}},A_{t^{\prime}})(A_{t^{\prime}}-\tilde{p}_{t^{\prime}})f_{t^{\prime}}(S_{t^{\prime}})^{\top}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{t,t^{\prime}}^{T}W_{t}\left(Y_{t+1}-g^{\star}_{t}(H_{t},A_{t})\right)(A_{t}-\tilde{p}_{t})f_{t}(S_{t})W_{t^{\prime}}\epsilon(H_{t^{\prime}},A_{t^{\prime}})(A_{t^{\prime}}-\tilde{p}_{t^{\prime}})f_{t^{\prime}}(S_{t^{\prime}})^{\top}\right]$

Here the first to second line uses the fact that
$f_{t}(S_{t})^{\top}\beta^{\star}+\Delta_{t}^{\perp}=g_{t}(H_{t},1)-g_{t}(H_{t},0)$
so we can then get
$\tilde{Y}^{(R)}_{t+1}-(A_{t}-\tilde{p}_{t})(f_{t}(S_{t})^{\top}\beta^{\star}+\Delta_{t}^{\perp})=Y_{t+1}-g^{\star}_{t}(H_{t},A_{t})$.
For $t\geq t^{\prime}$, by iterated expectation, we have:

$\displaystyle\mathbb{E}\left[\sum_{t,t^{\prime}}^{T}W_{t}\underbrace{\mathbb{E}\left[\left(Y_{t+1}-g^{\star}_{t}(H_{t},A_{t})\right)|H_{t},A_{t}\right]}_{=0}W_{t^{\prime}}\epsilon(H_{t^{\prime}},A_{t^{\prime}})(A_{t^{\prime}}-\tilde{p}_{t^{\prime}})f_{t^{\prime}}(S_{t^{\prime}})^{\top}(A_{t}-\tilde{p}_{t})f_{t}(S_{t})\right]$
$\displaystyle=0.$

For $t<t^{\prime}$, by iterated expectation, we have:

$\displaystyle\mathbb{E}\left[\sum_{t,t^{\prime}}^{T}W_{t}\left(Y_{t+1}-g^{\star}_{t}(H_{t},A_{t})\right)\mathbb{E}\left[W_{t^{\prime}}\epsilon(H_{t^{\prime}},A_{t^{\prime}})(A_{t^{\prime}}-\tilde{p}_{t^{\prime}})f_{t^{\prime}}(S_{t^{\prime}})^{\top}|H_{t^{\prime}}\right](A_{t}-\tilde{p}_{t})f_{t}(S_{t})\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{t,t^{\prime}}^{T}W_{t}\left(Y_{t+1}-g^{\star}_{t}(H_{t},A_{t})\right)\mathbb{E}\left[\tilde{p}_{t^{\prime}}(1-\tilde{p}_{t^{\prime}})\Delta_{t^{\prime}}^{\perp}f_{t^{\prime}}(S_{t^{\prime}})^{\top}|H_{t^{\prime}}\right](A_{t}-\tilde{p}_{t})f_{t}(S_{t})\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{t,t^{\prime}}^{T}W_{t}\mathbb{E}\left[\underbrace{\left(Y_{t+1}-g^{\star}_{t}(H_{t},A_{t})\right)\tilde{p}_{t^{\prime}}(1-\tilde{p}_{t^{\prime}})\Delta_{t^{\prime}}^{\perp}f_{t^{\prime}}(S_{t^{\prime}})^{\top}|A_{t},H_{t}}_{\text{conditionally
independent}}\right](A_{t}-\tilde{p}_{t})f_{t}(S_{t})\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{t,t^{\prime}}^{T}W_{t}\underbrace{\mathbb{E}\left[Y_{t+1}-g^{\star}_{t}(H_{t},A_{t})|A_{t},H_{t}\right]}_{=0}\mathbb{E}\left[\tilde{p}_{t^{\prime}}(1-\tilde{p}_{t^{\prime}})\Delta_{t^{\prime}}^{\perp}f_{t^{\prime}}(S_{t^{\prime}})^{\top}|A_{t},H_{t}\right](A_{t}-\tilde{p}_{t})f_{t}(S_{t})\right]$
$\displaystyle=$ $\displaystyle 0$

Therefore, the derivation above shows that $\Sigma^{(AR)}-\Sigma$ is negative
semidefinite. This indicates using the augmented R-WCLS to estimate treatment
effect $\beta^{\star}$ is more efficient than WCLS. In the case when we
estimate $\beta_{t}^{\star}$ nonparametrically rather than smoothing over
time, the interaction terms for $t\neq t^{\prime}$ don’t exist, therefore the
conclusion holds without the conditional independence assumption.

###  E.2 Efficiency comparison between R-WCLS and DR-WCLS

As stated in Section 4.2, the difference between R-WCLS and DR-WCLS criterion
is DR-WCLS replaces one term in R-WCLS with its expectation. Denote the
following:

$\displaystyle
M_{DR,t}=\tilde{\sigma}^{2}_{t}(S_{t})\left(\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta\right)f_{t}(S_{t}),$
$\displaystyle
M_{R,t}=W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))\left(\tilde{Y}_{t+1}^{(R)}-(A_{t}-\tilde{p}_{t})f_{t}(S_{t})^{\top}\beta\right)f_{t}(S_{t}).$

Here we apply the Conditional Jensen’s Inequality: since $g(x)=x^{2}$ is
convex, $\mathbb{E}[X^{2}|Z]\geq(\mathbb{E}[X|Z])^{2}$ holds. Then we have the
following inequality hold:

$\displaystyle\mathbb{E}\left[\left(\sum_{t=1}^{T}M_{R,t}\right)^{2}\,|\,H_{t}\right]\geq\left(\mathbb{E}\left[\sum_{t=1}^{T}M_{R,t}\,|\,H_{t}\right]\right)^{2}$

By iterative expectation, we can easily show:

$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}M_{R,t}\,|\,H_{t}\right]=\sum_{t=1}^{T}M_{DR,t},$

thus, it’s equivalently saying:

$\displaystyle\mathbb{E}\left[\left(\sum_{t=1}^{T}M_{R,t}\right)^{2}\,|\,H_{t}\right]\geq\left(\sum_{t=1}^{T}M_{DR,t}\right)^{2}$

Both sides of the inequality are functions of random variables collected in
$H_{t}$, thus the inequality still holds after taking expectation w.r.t.
$H_{t}$, which yields:

$\displaystyle\mathbb{E}\left[\left(\sum_{t=1}^{T}M_{R,t}\right)^{2}\right]\geq\mathbb{E}\left[\left(\sum_{t=1}^{T}M_{DR,t}\right)^{2}\right]$

Following this, the asymptotic variance of the DR-WCLS estimator can be
calculated as:

$\displaystyle\Sigma^{(DR)}=$
$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]^{-1}\times\mathbb{E}\left[\left(\sum_{t=1}^{T}M_{DR,t}\right)^{2}\right]\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]^{-1}$

Recall $\Sigma^{(R)}$ from previous proof:

$\displaystyle\Sigma^{(R)}=$
$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]^{-1}\times\mathbb{E}\left[\left(\sum_{t=1}^{T}M_{R,t}\right)^{2}\right]\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{p}_{t}(1|S_{t})(1-\tilde{p}_{t}(1|S_{t}))f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]^{-1}$

Therefore $\Sigma^{(DR)}-\Sigma^{(R)}$ is negative semidefinite. In
conclusion, we proved that the DR-WCLS estimator is more efficient than R-WCLS
estimator when we estimate the time-varying treatment effect.

##  Appendix F Proof of Corollary 6.1.1

###  F.1 Double robustness property

To derive the DR-WCLS criterion (14) with missing indicator $R_{t}$. Under
Assumption 2.1 and 6.1, the pseudo outcome $\tilde{Y}_{t+1}^{(DR)}$ can be
written as:

$\displaystyle\tilde{Y}_{t+1}^{(DR)}=$
$\displaystyle\beta(t;H_{t})+\frac{A_{t}R_{t}(Y_{t+1}-g(H_{t},A_{t}))}{p_{t}(A_{t},R_{t}|H_{t})}-\frac{(1-A_{t})R_{t}(Y_{t+1}-g(H_{t},A_{t}))}{p_{t}(A_{t},R_{t}|H_{t})}$
$\displaystyle=$
$\displaystyle\beta(t;H_{t})+\frac{A_{t}R_{t}(Y_{t+1}-g(H_{t},A_{t}))}{p(R_{t}|H_{t})p(A_{t}|H_{t})}-\frac{(1-A_{t})R_{t}(Y_{t+1}-g(H_{t},A_{t}))}{p(R_{t}|H_{t})p(A_{t}|H_{t})}$
$\displaystyle=$
$\displaystyle\beta(t;H_{t})+\frac{R_{t}}{p(R_{t}|H_{t})}\left[\frac{A_{t}(Y_{t+1}-g(H_{t},A_{t}))}{p(A_{t}|H_{t})}-\frac{(1-A_{t})(Y_{t+1}-g(H_{t},A_{t}))}{p(A_{t}|H_{t})}\right]$
$\displaystyle=$ $\displaystyle\beta(t;H_{t})+\frac{{\bf
1}(R_{t}=1)}{p(R_{t}|H_{t})}\frac{W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g_{t}(H_{t},A_{t}))}{\tilde{\sigma}^{2}_{t}(S_{t})}$

and the corresponding estimating equation is:

$\displaystyle\mathbb{P}_{n}\Bigg{[}\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\Big{(}\frac{{\bf
1}(R_{t}=1)}{p(R_{t}|H_{t})}\frac{W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g_{t}(H_{t},A_{t}))}{\tilde{\sigma}^{2}_{t}(S_{t})}+\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta\Big{)}f_{t}\Bigg{]}$

Furthermore, based on previous proofs, we can conclude that the
$\hat{\beta}_{n}$ obtained by solving the above estimating equation is doubly
robust.

###  F.2 Asymptotic normality

We use the same notation as the previous section, assume
$\hat{\beta}^{(DR)}_{n}$ solves:

$\displaystyle 0$
$\displaystyle=\mathbb{P}_{n}\Bigg{[}\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\Big{(}\frac{{\bf
1}(R_{t}=1)}{p(R_{t}|H_{t})}\frac{W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g_{t}(H_{t},A_{t}))}{\tilde{\sigma}^{2}_{t}(S_{t})}+\beta(t;H_{t})-f_{t}(S_{t})^{\top}\hat{\beta}^{(DR)}_{n}\Big{)}f_{t}\Bigg{]}$

When the true randomization probability $p_{t}=p(A_{t}|H_{t})$ and missing
mechanism $p_{t}^{R}=p(R_{t}|H_{t})$ are unknown, we have the weight $W_{t}$
estimated by
$\hat{W}_{t}=\tilde{p}_{t}(A_{t}|S_{t})/\hat{p}_{t}(A_{t}|H_{t})$, and missing
mechanism estimated by $\hat{p}(R_{t}|H_{t})$. Then the estimating equation
can be decomposed as:

$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\left(\frac{R_{t}}{\hat{p}(R_{t}|H_{t})}\frac{\tilde{p}_{t}(A_{t}|S_{t})(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-\hat{g}_{t}(H_{t},A_{t}))}{\hat{p}_{t}(A_{t}|H_{t})\tilde{\sigma}^{2}_{t}(S_{t})}+\hat{\beta}(t;H_{t})-f_{t}(S_{t})^{\top}\beta^{(DR)}\right)f_{t}(S_{t})\right]$
$\displaystyle=$
$\displaystyle\mathbb{P}_{n}\Bigg{[}\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\Big{(}\frac{R_{t}\tilde{p}_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g^{\star}_{t}(H_{t},A_{t})+g^{\star}_{t}(H_{t},A_{t})-\hat{g}_{t}(H_{t},A_{t}))}{\tilde{\sigma}^{2}_{t}(S_{t})}\left(\frac{1}{\hat{p}_{t}\hat{p}_{t}^{R}}-\frac{1}{p_{t}p_{t}^{R}}+\frac{1}{p_{t}p_{t}^{R}}\right)$
$\displaystyle+\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta^{\star}+(\hat{\beta}(t;H_{t})-\beta(t;H_{t}))-f_{t}(S_{t})^{\top}(\beta^{(DR)}-\beta^{\star})\Big{)}f_{t}(S_{t})\Bigg{]}$
$\displaystyle=$
$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})\left(\frac{R_{t}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g^{\star}_{t}(H_{t},A_{t}))}{p_{t}^{R}\tilde{\sigma}^{2}_{t}(S_{t})}+\beta(t;H_{t})-f_{t}(S_{t})^{\top}\beta^{\star}\right)f_{t}(S_{t})\right]+$
$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}R_{t}\tilde{p}_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g^{\star}_{t}(H_{t},A_{t}))\left(\frac{1}{\hat{p}_{t}^{R}\hat{p}_{t}}-\frac{1}{p_{t}^{R}p_{t}}\right)f_{t}(S_{t})\right]+$
$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}R_{t}\tilde{p}_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(g^{\star}_{t}(H_{t},A_{t})-\hat{g}_{t}(H_{t},A_{t}))\left(\frac{1}{\hat{p}_{t}^{R}\hat{p}_{t}}-\frac{1}{p_{t}^{R}p_{t}}\right)f_{t}(S_{t})\right]+$
$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}\frac{R_{t}}{p_{t}^{R}}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(g^{\star}(H_{t},A_{t})-\hat{g}_{t}(H_{t},A_{t}))f_{t}(S_{t})\right]+$
$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})(\hat{\beta}(t;H_{t})-\beta(t;H_{t}))f_{t}(S_{t})\right]-$
$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})f_{t}(S_{t})f^{\top}_{t}(S_{t})\right](\hat{\beta}^{(DR)}_{n}-\beta^{\star})$

Because

$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}R_{t}\tilde{p}_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(Y_{t+1}-g^{\star}_{t}(H_{t},A_{t}))\left(\frac{1}{\hat{p}_{t}^{R}\hat{p}_{t}}-\frac{1}{p_{t}^{R}p_{t}}\right)f_{t}(S_{t})\right]\overset{P}{\to}0,$
$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})f_{t}(S_{t})f^{\top}_{t}(S_{t})\right]\overset{P}{\to}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})f_{t}(S_{t})f^{\top}_{t}(S_{t})\right],$

and

$\displaystyle\mathbb{P}_{n}\left[\sum_{t=1}^{T}\frac{R_{t}}{p_{t}^{R}}W_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(g^{\star}(H_{t},A_{t})-\hat{g}_{t}(H_{t},A_{t}))f_{t}(S_{t})\right]+$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathbb{P}_{n}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})(\hat{\beta}(t;H_{t})-\beta(t;H_{t}))f_{t}(S_{t})\right]\overset{P}{\to}0.$

Apart from the nicely-behaved terms above, the only term that might be
problematic is:

$\mathbb{P}_{n}\left[\sum_{t=1}^{T}R_{t}\tilde{p}_{t}(A_{t}-\tilde{p}_{t}(1|S_{t}))(g^{\star}_{t}(H_{t},A_{t})-\hat{g}_{t}(H_{t},A_{t}))\left(\frac{1}{\hat{p}_{t}^{R}\hat{p}_{t}}-\frac{1}{p_{t}^{R}p_{t}}\right)f_{t}(S_{t})\right],$

which will converge to:

$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}\left(\underbrace{\sum_{a\in\\{0,1\\}}{\bf
c}(a)(p_{t}^{R}(1|H_{t})p_{t}(a|H_{t})-\hat{p}_{t}^{R}(1|H_{t})\hat{p}_{t}(a|H_{t}))(g^{\star}(H_{t},a)-\hat{g}(H_{t},a))}_{\text{(II)}}\right)f_{t}(S_{t})\right],$

where ${\bf
c}(a)=\frac{\tilde{\sigma}^{2}_{t}(S_{t})/\hat{p}_{t}^{R}(1|H_{t})}{a\hat{p}_{t}(1|H_{t})+(a-1)(1-\hat{p}_{t}(1|H_{t}))}$.
In our context, $T$ is finite and fixed. Therefore, by the fact that
$\hat{p}_{t}(1|H_{t}))$ is bounded away from zero and one, along with the
Cauchy–Schwarz inequality, we have that (up to a multiplicative constant) term
(II) is bounded above by:

$\mathbf{\hat{B}}^{R}=\mathbb{E}\left[\sum_{t=1}^{T}\sum_{a\in\\{0,1\\}}\left\|p_{t}^{R}(1|H_{t})p_{t}(a|H_{t})-\hat{p}_{t}^{R}(1|H_{t})\hat{p}_{t}(a|H_{t})\right\|\left\|g^{\star}(H_{t},a)-\hat{g}(H_{t},a)\right\|\right].$
(34)

Same argument as in the previous section, if $\hat{p}(a|H_{t})$ and
$\hat{p}_{t}^{R}(1|H_{t})$ are based on a correctly specified parametric
model, so that
$\left\|\hat{p}_{t}^{R}(1|H_{t})\hat{p}_{t}(a|H_{t})-p_{t}^{R}(1|H_{t})p_{t}(a|H_{t})\right\|=O_{p}(n^{-1/2})$,
then we only need $\hat{g}(H_{t},a)$ to be
consistent,$\left\|g^{\star}(H_{t},a)-\hat{g}(H_{t},a)\right\|=o_{p}(1)$, to
make $\mathbf{\hat{B}}^{R}$ asymptotically negligible. Thus if we know the
treatment and data missingness mechanism, the outcome model can be very
flexible. Another way to achieve efficiency is if we have both
$\left\|\hat{p}_{t}^{R}(1|H_{t})\hat{p}_{t}(a|H_{t})-p_{t}^{R}(1|H_{t})p_{t}(a|H_{t})\right\|=o_{p}(n^{-1/4})$
and $\left\|g^{\star}(H_{t},a)-\hat{g}(H_{t},a)\right\|=o_{p}(n^{-1/4})$, so
that their product term is $o_{p}(n^{-1/2})$ and asymptotically negligible
(Kennedy, 2016). This of course occurs if both $\hat{g}(H_{t},a)$ and
$\hat{p}_{t}^{R}(1|H_{t})\hat{p}_{t}(a|H_{t})$ are based on correctly
specified models, but it can also hold even for estimators that are very
flexible and not based on parametric models.

Assuming we have nuisance estimates that can make $\mathbf{\hat{B}}^{R}$
asymptotically negligible, then the DR-WCLS estimator satisfies:

$\displaystyle n^{1/2}(\hat{\beta}^{(DR)}_{n}-\beta^{\star})$
$\displaystyle=n^{1/2}~{}\mathbb{P}_{n}\Bigg{[}\sum_{t=1}^{T}~{}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})f_{t}(S_{t})f^{\top}_{t}(S_{t})\right]^{-1}\tilde{\sigma}^{2}_{t}(S_{t})(\beta(t;H_{t})-f_{t}(S_{t})\beta^{\star})f_{t}(S_{t})\Bigg{]}+o_{p}(1),$

and it is efficient with influence function:

$\displaystyle\sum_{t=1}^{T}~{}\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})f_{t}(S_{t})f_{t}(S_{t})^{\top}\right]^{-1}\tilde{\sigma}^{2}_{t}(S_{t})(\beta(t;H_{t})-f_{t}(S_{t})\beta^{\star})f_{t}(S_{t}).$

Under moment conditions, we have asymptotic normality with variance given by
$\Sigma^{R}_{DR}=Q^{-1}WQ^{-1}$, where

$\displaystyle Q$
$\displaystyle=\mathbb{E}\left[\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})f_{t}(S_{t})f_{t}(S_{t})^{\top}\right],$
$\displaystyle W$
$\displaystyle=\mathbb{E}\left[\left(\sum_{t=1}^{T}\tilde{\sigma}^{2}_{t}(S_{t})(\beta(t;H_{t})-f_{t}(S_{t})\beta^{\star})f_{t}(S_{t})\right)^{2}\right].$

###  F.3 Algorithm

##### Step I

Let $K$ be a fixed integer. Form a K-fold random partition of
$\\{1,2,\dots,N\\}$ by dividing it to equal parts, each of size $n:=N/K$,
assuming $N$ is a multiple of $K$. Form each set $I_{k}$, let
$I_{k}^{\complement}$ denote the observation indices that are not in $I_{k}$.

##### Step II

For each fold, use any supervised learning algorithm to estimate the
appropriate working models. Let $\hat{g}_{t}^{(k)}(H_{t},A_{t})$,
$\hat{p}_{t}^{(k)}(1|H_{t})$, $\hat{p}_{t}^{(k)}(R_{t}|H_{t})$ and
$\hat{\tilde{p}}_{t}^{(k)}(1|S_{t})$ denote the estimates for
$\mathbb{E}[Y_{t+1}|H_{t},A_{t}]$,
$\mathbb{E}[A_{t}|H_{t}]$,$\mathbb{E}[R_{t}|H_{t}]$, and
$\mathbb{E}[A_{t}|S_{t}]$ respectively using individuals in
$I_{k}^{\complement}$, i.e., estimates of the nuisance parameters the $k$th
fold.

##### Step III

Construct the pseudo-outcomes and perform weighted regression estimation:

$\tilde{Y}^{(DR)}_{t+1}:=\frac{{\bf