# Automated, efficient and model-free inference for randomized clinical trials
via data-driven covariate adjustment

Kelly Van Lancker1, Iván Díaz2, and Stijn Vansteelandt1

(1Department of Applied Mathematics, Computer Science and Statistics, Ghent
University, Ghent, Belgium
2Division of Biostatistics, Department of Population Health, New York
University School of Medicine, NY, United States)

###### Abstract

In May 2023, the U.S. Food and Drug Administration (FDA) released guidance for
industry on “Adjustment for Covariates in Randomized Clinical Trials for Drugs
and Biological Products”. Covariate adjustment is a statistical analysis
method for improving precision and power in clinical trials by adjusting for
pre-specified, prognostic baseline variables. Though recommended by the FDA
and the European Medicines Agency (EMA), many trials do not exploit the
available information in baseline variables or make use only of the baseline
measurement of the outcome. This is likely (partly) due to the regulatory
mandate to pre-specify baseline covariates for adjustment, leading to
challenges in determining appropriate covariates and their functional forms.

We will explore the potential of automated data-adaptive methods, such as
machine learning algorithms, for covariate adjustment, addressing the
challenge of pre-specification. Specifically, our approach allows the use of
complex models or machine learning algorithms without compromising the
interpretation or validity of the treatment effect estimate and its
corresponding standard error, even in the presence of misspecified outcome
working models. This contrasts the majority of competing works which assume
correct model specification for the validity of standard errors. Our proposed
estimators either necessitate ultra-sparsity in the outcome model (which can
be relaxed by limiting the number of predictors in the model) or necessitate
integration with sample splitting to enhance their performance. As such, we
will arrive at simple estimators and standard errors for the marginal
treatment effect in randomized clinical trials, which exploit data-adaptive
outcome predictions based on prognostic baseline covariates, and have low (or
no) bias in finite samples even when those predictions are themselves biased.

We provide a detailed methodology overview as well as empirical study results.
The findings offer a promising avenue for improving the statistical power of
trial analyses through automated covariate adjustment.

##  1 Introduction

For many years, the debate on whether or not to adjust for baseline covariates
in the primary intention-to-treat analysis of randomized clinical trials has
been obscured by the lack of focus on a single pre-specified estimand.
Unadjusted comparisons deliver marginal measures of treatment effect, but tend
to be inefficient as a result of ignoring potential imbalances in prognostic
baseline covariates. Correction for such covariates happens in adjusted
comparisons delivered by standard regression analyses; however, at the risk of
bias resulting from misspecification of the regression model that is needed
for adjustment (with the exception of targeting marginal effects with linear
regression). Moreover, the potential efficiency gain delivered by those
analyses is difficult to quantify (Robinson and Jewell, 1991) as they,
conversely, deliver conditional measures of treatment effect, whose meaning
and magnitude may depend on precisely what covariates have been adjusted for
(Daniel et al., 2021). Tsiatis et al. (2008) brought much clarity to the
debate. He focused on estimation of the marginal treatment effect and found
(semi-parametric) efficient estimators to exploit covariate adjustment for
baseline covariates in order to boost precision, whilst being insulated - in
sufficiently large samples - against the bias that might otherwise result from
model misspecification.

Despite these advances, practical implementation of (semi-parametric)
efficient estimators of the marginal treatment effect in randomized clinical
trials has been hindered by the need to pre-specify the baseline covariates
for adjustment, as mandated by regulatory agencies such as the U.S. Food and
Drug Administration (FDA) and the European Medicine Agency (EMA) (FDA and EMA,
1998; EMA, 2015). This poses challenges in determining the appropriate
covariates and their functional forms (e.g., which interactions between
covariates will be included in the model for the outcome) prior to seeing the
data, often leading to difficult or even infeasible tasks in many trials. To
address this, there seems promise in the use of automated data-adaptive
methods, which ‘adapt’ to the data by making use of automated variable
selection strategies (e.g., stepwise AIC-based selection, the lasso, …) or
even machine learning methods (e.g., random forest regression, gradient
boosting, neural networks, …). Such algorithms allow for objective, guided
selection of variables and construction of statistical models, and can be pre-
specified (Tsiatis et al., 2008; Van der Laan et al., 2011; Wager et al.,
2016; Wu and Gagnon-Bartsch, 2018; Vansteelandt, 2021).

Nonetheless, the application of data-adaptive methods also introduces
potential concerns. These methods are generally optimized for predictive
accuracy, and may be sub-optimal for the estimation of treatment effects. In
particular, they may deliver finite-sample biases that are non-negligible to
the extent that they exceed the magnitude of the standard error at all sample
sizes, thereby distorting subsequent inferences (Chernozhukov et al., 2018).
In the literature on variable and model selection, this is often referred to
as the problem of post-selection inference (Leeb and Pötscher, 2008; Belloni
et al., 2014a, 2017). It stems from the estimator’s distribution being
different depending on the chosen model, with standard confidence intervals
disregarding the resulting uncertainty in model selection.

Debiasing procedures have been developed to remedy this, thereby ensuring more
robust treatment effect estimates. Targeted learning procedures (Moore and van
der Laan, 2009; Van der Laan et al., 2011) work by cleverly tuning the
predictions delivered by data-adaptive algorithms in a way that removes bias
in the resulting treatment effect estimator. Debiased machine learning methods
(Chernozhukov et al., 2018) instead make use of default predictions, but
subsequently debias the obtained treatment effect estimate. Interestingly, the
(semi-parametric) efficient estimators of Tsiatis et al. (2008) incorporate
such debiasing by default. Unlike many other covariate-adjusted estimators of
treatment effect, they can thus readily be used in conjunction with data-
adaptive procedures, such as automated variable selection or machine learning.
They thereby deliver treatment effect estimates that invoke covariate
adjustment in a manner that can be pre-specified - as it can be automated -
and have low bias even when the employed data-adaptive predictions are
themselves biased (e.g., as a result of model misspecification or the use of
poorly fitting/converging machine learning predictions). Wager et al. (2016)
shows the resulting estimators to be exactly unbiased, even in finite samples,
provided linear models are used and a specific sample splitting procedure
(known as cross-fitting) is employed, whereby the model building process runs
on a different part of the data than the part on which the estimator is
calculated. This article expands the unbiasedness result for linear models to
data-driven covariate adjustment in a broader context. We thereby demonstrate
that achieving exact unbiasedness is even feasible when using machine learning
predictions, provided cross-fitting is employed, across various types of
endpoints.

Unfortunately, however, this robustness to model misspecification or poor
convergence of machine learning predictions does not generally extend to
standard errors, whether obtained using sandwich estimators or the bootstrap.
In particular, the standard errors delivered by the proposal of Tsiatis et al.
(2008) or by more general targeted learning procedures (Van der Laan et al.,
2011) may be biased when based on poor data-adaptive predictions. Vermeulen
and Vansteelandt (2015) remedied this problem by employing additional
debiasing procedures for parametric, pre-specified models, which were later
extended to include automated variable selection in parametric models (Avagyan
and Vansteelandt, 2022) and even to arbitrary regression learners (Benkeser et
al., 2017); however, the estimators of Benkeser et al. (2017) are non-regular
and may therefore show erratic behaviour under certain data-generating
distributions (Dukes et al., 2021). Williams et al. (2021) showed that this
problem is nonetheless easy to remedy for targeted learning procedures in
randomized clinical trials w.r.t. ordinal and time-to-event outcomes. Their
results focus on targeted learning (1) with a propensity score estimated
independently of baseline covariates and (2) with outcome models obtained via
cross-fitting. Bannick et al. (2023) extended these results to covariate-
adaptive randomization.

Wu and Gagnon-Bartsch (2018) explore similar concepts under the Neyman–Rubin
model, assuming fixed potential outcomes and baseline covariates, with
randomness solely in treatment assignment. This contrasts with the approach
taken in this article, which assumes participants are drawn from a super-
population. Specifically, Wu and Gagnon-Bartsch (2018) proposed a ‘leave-one-
out-potential outcomes’ estimator with a conservative variance estimator where
automatic prediction can also be performed using various regression techniques
such as linear regression or random forests. The drawback of this estimator
lies in its complexity of implementation and the ambiguity surrounding the
appropriateness of Wald-type confidence intervals derived from the normal
distribution.

In this paper, we extend the results of Williams et al. (2021) for data-driven
covariate adjustment, without necessitating use of the more general, but more
technically demanding strategies listed above. Our results support parametric
estimation of the propensity score (e.g., using a logistic regression model)
to improve precision (see Section 6). Specifically, we show that valid
inference can be obtained in the non-trivial case where the propensity score
may become a complicated expression of the estimated parameters, even when the
outcome model is misspecified. We moreover illustrate that exact unbiasedness
can be attained, even when leveraging machine learning predictions, as long as
cross-fitting is utilized. In addition, we show that, with some variable
selection approaches, valid inference can be attained without the need for
sample splitting. As such, we will arrive at simple estimators and standard
errors for the marginal treatment effect in randomized clinical trials, which
exploit data-adaptive outcome predictions based on prognostic baseline
covariates, and have low (or no) bias even when those predictions are
themselves biased. These features are appealing as they allow for the
utilization of exceedingly complex models, such as those involving splines or
machine learning algorithms, without distorting the interpretation or validity
of the treatment effect estimate. This flexibility guarantees that the
analysis remains solely evidence-based and leverages covariate information to
maximize statistical power (Tsiatis et al., 2008).

The remainder of this paper is organized as follows. In Section 2 we introduce
some notation and the estimand of interest, Section 3 presents a simple
covariate-adjusted estimator, Section 4 provides a detailed overview of the
proposed data-adaptive methodology, Section 5 discusses the results of our
empirical study, Section 6 touches on the results when estimating the
propoensity score parametrically, and finally, Section 7 concludes with a
summary of our findings and potential future research directions.

##  2 Notation and Estimand

For participant $i$ ($i=1,\dots n$), denote by $Z_{i}$ the study arm
assignment (1 for experimental treatment and 0 for control), by $X_{i}$ a
vector of baseline variables measured before treatment assignment, and by
$Y_{i}$ the primary outcome of interest which can be continuous, binary or
ordinal. When followed up completely, participant $i$ contributes full data
$\mathbf{O}_{i}=(X_{i},Z_{i},Y_{i})$.

We are interested in making inferences about marginal / unconditional
intention-to-treat effects, which we will denote by $\theta$ and refer to as
the policy estimand of interest. Examples include the difference in means of
the primary outcome between study arms for continuous and binary outcomes,
relative risk and odds ratio for binary outcomes, and Mann-Whitney estimand
for ordinal outcomes. In what follows, without loss of generality, we focus on
estimation and inference for

$\theta=E\left(Y|Z=1\right)-E\left(Y|Z=0\right),$

based on different (covariate-adjusted) estimators. According to the law of
iterated expectation, this expression is equivalent to
$E\\{E\left(Y|Z=1,X\right)-E\left(Y|Z=0,X\right)\\}$. This provides insight
into the covariate adjustment technique discussed in Section 3. It indicates
the procedure of first estimating subgroup effects (e.g., via
logistic/linear/…regression models) and then determining the standardized
value by averaging these subgroup effects across all participants.

##  3 Regular Covariate-Adjusted Estimation

Unadjusted estimators of the estimands are readily obtained by replacing
population by sample means in the estimand definition. For example, for the
difference in (population) mean outcomes estimand (i.e.,
$\theta=E\left(Y|Z=1\right)-E\left(Y|Z=0\right)$), the unadjusted estimator is
the difference in sample means.
To take full advantage of prognostic baseline covariates, we focus on a
generalization of the covariate adjusted estimator suggested in the recent
U.S. Food and Drug Administration guidance draft on covariate adjustment (FDA,
2021). To focus ideas, we first consider a model-based covariate-adjusted
estimator of $\mu_{1}=E(Y|Z=1)$:

  * Step 1:

Model fitting
Regress the outcome $Y$ on pre-specified baseline covariates $X$ among
participants with $Z=1$, by fitting a generalized linear regression model
(e.g., logistic, linear, …) with a canonical link (e.g., logit, identity, …)
and intercept via maximum likelihood. For example, one could model
$E(Y|Z=1,X)$ by
$h_{1}(X;\boldsymbol{\gamma})=logit^{-1}(\gamma_{0}+\gamma_{1}X)$ for a binary
endpoint $Y$ among participants with $Z=1$.

  * Step 2:

Predicting
Use the fitted regression model in Step 1 to compute a prediction of the
outcome under $Z=1$ for each participant, that is,
$\hat{Y}_{1,i}=h_{1}(X_{i};\hat{\boldsymbol{\gamma}})=logit^{-1}(\hat{\gamma}_{0}+\hat{\gamma}_{1}X_{i})$
where
$\hat{\boldsymbol{\gamma}}=(\hat{\gamma}_{0},\hat{\gamma}_{1})^{\prime}$.

  * Step 3:

Averaging
Take the average of these predicted outcomes over all study participants to
obtain an estimate for the average outcome under $Z=1$:
$\hat{\mu}_{1}=n^{-1}\sum_{i=1}^{n}\hat{Y}_{1,i}$.

Similarly, a covariate-adjusted estimator $\hat{\mu}_{0}$ for
$\mu_{0}=E(Y|Z=0)$ can be obtained by

$\hat{\mu}_{0}=n^{-1}\sum_{i=1}^{n}\hat{Y}_{0,i},$

where the predicted values $h_{0}(X_{i},\hat{\boldsymbol{\eta}})$ for each
participant are obtained via regression of $Y$ on $X$ among the participants
with $Z=0$ using a canonical generalized linear working model
$h_{0}(X,\boldsymbol{\eta})$ with intercept for $E(Y|Z=0,X)$, fitted via
maximum likelihood. For example, one could specify
$h_{0}(X,\boldsymbol{\eta})=logit^{-1}(\eta_{0}+\eta_{1}X)$ for a binary
endpoint $Y$. Then,
$\hat{\boldsymbol{\eta}}=(\hat{\eta}_{0},\hat{\eta}_{1})^{\prime}$ and
$\hat{Y}_{0,i}=h_{0}(X_{i};\hat{\boldsymbol{\eta}})$. A covariate-adjusted
estimator for the estimand of interest $\theta$ can then be obtained as
$\hat{\theta}=\hat{\mu}_{1}-\hat{\mu}_{0}$. It is known as a standardization,
G-computation or imputation estimator (Tsiatis et al., 2008).

Importantly, this covariate-adjusted estimator is robust to model
misspecification (i.e., the estimator will be unbiased in large samples even
if the outcome working models are misspecified) and has a high potential to
improve precision (i.e., to deliver tighter confidence intervals and higher
power), making it a low-risk, high-reward method (Tsiatis et al., 2008;
Benkeser et al., 2020). This robustness to model misspecification (in large
samples) only holds under a prediction unbiasedness condition:

$\sum_{i:Z_{i}=z}\hat{Y}_{z,i}=\sum_{i:Z_{i}=z}Y_{i},$ (1)

where $\hat{Y}_{z,i}$ is a prediction under treatment $z\in\\{0,1\\}$ for
participant $i$ ($i=1,\dots,n$). This holds for the considered estimator as it
is guaranteed for canonical generalized linear models fitted via maximum
likelihood. For other models, additional debiasing is needed to restore model
robustness (Robins et al., 1994). With concern for finite-sample bias, the
estimator can even be made finite-sample-unbiased by using sample splitting
(see later in Section 4.3) and employing the known randomization probabilities
(see Appendix A.2.2 for a proof).

##  4 Data-Adaptive Covariate-Adjusted Estimation

Pre-specification of the optimal baseline covariates that will be adjusted for
along with the functional form of the outcome model, as mandated by regulatory
agencies (FDA and EMA, 1998; EMA, 2015), is a difficult, if not impossible,
task in many trials. There thus seems promise in pre-specifying how to use the
collected data to guide the selection of variables and a statistical model,
along with respect to what covariate set. In particular, data-adaptive methods
provide a better opportunity to explore covariate-outcome relationships. The
flexibility afforded by the opportunity to freely exploit modeling methods and
expertise, allows covariate information to be optimally used to gain
efficiency (Tsiatis et al., 2008).

Specifically, in this article, we show that similar methods to the one
proposed in the FDA guidance remain valid when the selection of baseline
covariates is based on the trial data, so long as one pre-specifies how the
selection will be done, that is, based on what covariate set and algorithm.
For instance, a valid data-adaptive method can be obtained by replacing Step 1
by Step 1a and 1b, and similarly for the model fitting step for estimation of
$\mu_{0}$ (Tsiatis et al., 2008):

  * Step 1a:

Model selection Fit a generalized linear regression model with a canonical
link (e.g., logistic, linear, …) that regresses the outcome $Y$ on pre-
specified baseline covariates $X$ among the treated participants (i.e., with
$Z=1$), using a data-adaptive method (e.g., stepwise regression or Lasso).
Select the covariates — which can for instance include interaction terms or
higher order terms — with non-zero estimated coefficients.

  * Step 1b:

Model fitting Fit a generalized linear regression model with a canonical link
(e.g., logistic, linear, …) and intercept via maximum likelihood that
regresses the outcome $Y$ on all covariates selected in Step 1a among the
treated participants. This regression may also include additional variables
(such as stratification factors) that were not selected in the previous steps,
but that were identified a priori as being important.

The estimated treatment effect is denoted by
$\hat{\theta}_{DA}=\hat{\mu}_{1,DA}-\hat{\mu}_{0,DA}$, where ‘DA’ stands for
‘data-adaptive’.

###  4.1 Characteristics

The estimator discussed in Section 4 inherits the appealing properties
discussed in Section 3, provided that canonical generalized linear models are
used and fitted via maximum likelihood (which ensures the prediction
unbiasedness condition (1) to hold); for other models, additional debiasing is
needed (see Section 4.3). In particular, misspecification of the outcome model
or poor covariate selection then does not introduce bias in the estimated
treatment effect $\hat{\theta}_{DA}$ in large samples (and even in finite
samples when making use of the known (instead of estimated) randomization
probabilities and when, moreover, sample splitting is used; see Section 4.3).
This is the case even when the set of covariates to choose from is large
relative to the sample size, as long as there is ultra-sparsity (see Appendix
A.1), in the sense that the number of predictors with non-zero coefficients in
the model is relatively small. This condition is relatively weak since the
user has full freedom what (and how many) predictors to include in the model.

Unfortunately, this robustness to model misspecification or (imperfect)
covariate selection does not generally extend to standard error estimators. In
particular, the available theory on inference for targeted or debiased
learning estimators, like the estimator of Section 4 (achieved by fitting a
canonical generalized linear model in Step 1b, guaranteeing condition (1)),
assumes correctly specified models (Van der Laan et al., 2011; Wager et al.,
2016; Chernozhukov et al., 2018), with some exceptions that invoke additional
debiasing procedures (Vermeulen and Vansteelandt, 2015; Benkeser et al., 2017;
Avagyan and Vansteelandt, 2022). In Appendix A.3, we show that, surprisingly,
correct model specification is not a requirement for the validity of standard
errors when randomization is independent of baseline covariates (Williams et
al., 2021; Bannick et al., 2023). In particular, we show that the use of data-
adaptive model selection techniques to determine which covariates to include
in the outcome working model does not affect the validity of large sample
inferences, even when the outcome model is misspecified. This is even so when
the set of covariates to choose from is large relative to the sample size,
provided that there is ultra-sparsity. In that case, inference can proceed as
if the covariates were pre-selected before seeing the data. Specifically, the
asymptotic variance of $\hat{\theta}_{DA}=\hat{\mu}_{1,DA}-\hat{\mu}_{0,DA}$
can then be calculated as $1/n$ times the sample variance of the values

$\displaystyle\frac{Z_{i}}{\hat{\pi}}(Y_{i}-\hat{Y}_{1,i})+\hat{Y}_{1,i}-\left[\frac{1-Z_{i}}{1-\hat{\pi}}(Y_{i}-\hat{Y}_{0,i})+\hat{Y}_{0,i}\right],$
(2)

with $\hat{\pi}=\hat{P}(Z=1)=n^{-1}\sum_{i=1}^{n}Z_{i}$ the empirical
randomization probability and $\hat{Y}_{z,i}$ ($z\in\\{0,1\\}$, $i=1,\dots,n$)
the predictions obtained in the data-adaptive procedure in Section 4.

This is a strong result as it implies that, besides the aforementioned
sparsity conditions, only weak assumptions (see Appendix A.1) are needed to
ensure that the data-adaptive approach in Section 4 is insensitive to
selection mistakes in large samples (i.e., preserves Type I error). A proof is
given in Appendix A.3.

When the outcome model is correctly specified, a consistent variable selection
procedure is used (such as the Lasso or stepwise regression) and the sample
size is sufficiently large, the above data-adaptive covariate-adjusted
estimator has the further advantage of being at least as efficient as all
other estimators that are unbiased under simple randomization (such as the
unadjusted estimator) (Yang and Tsiatis, 2001; Tsiatis, 2006).

Even though the misspecification of outcome working models might result in a
loss of precision (in comparison to using correctly specified models), in our
experience the proposal is rarely less efficient than an unadjusted analysis,
though this cannnot be theoretically guaranteed. However, by relating our
estimators to the literature on empirical efficiency maximization (Rubin and
van der Laan, 2008), some form of guaranteed efficiency over a broad class of
estimators (which includes the unadjusted estimator; see Appendix A.6) is
attainable, even when the outcome model is misspecified. Concretely, it
follows from Rubin and van der Laan (2008), Cao et al. (2009), Gruber and van
der Laan (2012) and others (see also Appendix A.6) that the estimator with
minimal variance can be obtained by selecting and fitting the outcome model to
have minimal mean squared error (MSE). This estimator is then moreover at
least as efficient as the unadjusted estimator. However, to attain this
efficiency guarantee, logistic regression models must be fitted via least
squares instead of maximum likelihood estimation. As other appealing
properties – such as the robustness against model misspecification – will no
longer hold when using least squares (as the prediction unbiasedness condition
(1) is no longer met), we discuss a slightly different approach (using sample
splitting) in Section 4.3 which has all the aforementioned properties
including this efficiency property when the outcome model is misspecified.

All these appealing characteristics – which are summarized in Table 1 – also
hold when the association between $Z$ and $X$ is parametrically estimated
(e.g., using a logistic regression model), as demonstrated in Appendix A.5.2.
This extension can be accomplished by fitting a weighted regression model in
Step 1b with weights equal to the inverse of the probability of receiving the
assigned treatment conditional on baseline covariates. Further elaboration on
this is provided in Section 6.

Table 1:  Overview of important characteristics of the different covariate-adjusted estimators. Sample splitting: whether sample splitting is considered; Algorithm: the prediction algorithms that are considered [Can. GLMs: only allowing for automated variable selection approaches (e.g., lasso, stepwise AIC-based selection, …) for canonical generalized linear models fitted via maximum likelihood; ML: also allowing for more geberal machine learning algorithms; GLMs∗: automated variable selection approach for parametric and non-parametric models, for which the prediction unbiasedness condition (1) does not necessarily hold]; SP efficient: whether the estimators are semi-parametric efficient; Sample bounded: whether estimates of $\mu_{1}$ and $\mu_{0}$ are bounded to the outcomes in the sample; Robust estimation: whether estimators are unbiased in large samples when the outcome working models are misspecified [Asympt.: estimator is (only) asymptotically unbiased, when models are misspecified; Exact∗∗: estimator is exactly unbiased - also in finite samples - as long as true $\pi$ has been used, even when models are misspecified]; SN w/o SS: can one avoid sample splitting when testing strong null (without any additional assumptions)? [${\dagger}$: We refer to Section 4.4 as it requires an adaptation of the regular and data-adaptive covariate-adjusted estimators with canonical GLMs.] Method | Sample | Algorithm | SP | Sample | Robust | SN
---|---|---|---|---|---|---
| splitting |  | efficient | bounded | estimation | w/o SS
Unadjusted | No | / | / | $\surd$ | Exact∗∗ | $\surd$
Regular cov. adj. | No | Can. GLMs | $\surd$ | $\surd$ | Asympt. | $\surd^{\dagger}$
Data-adaptive cov. adj. | No | Can. GLMs | $\surd$ | $\surd$ | Asympt. | $\surd^{\dagger}$
| Yes | ML | $\surd$ | / | Exact∗∗ | $\surd$
TMLE | No | GLMs∗ | $\surd$ | $\surd$ | Asympt. | /
| Yes | ML | $\surd$ | $\surd$ | Asympt. | /

###  4.2 Connection with Targeted Maximum Likelihood Estimation

The proposed approach is closely related to targeted maximum likelihood
estimation (TMLE) (see e.g., Rosenblum and Van Der Laan, 2010), which conducts
Step 1b in a slightly different way in order to enforce the prediction
unbiasedness condition (1). In particular, the initial predictions from Step
1a (see Section 4), denoted by $\hat{Y}_{1,i}^{init}$, are updated by fitting
the following regression model with single parameter $\epsilon$ and offset
based on $\hat{Y}_{1,i}^{init}$:

$\displaystyle
E(Y|Z=1,X)=g^{-1}\left\\{g\left(\hat{Y}_{1,i}^{init}\right)+\epsilon\right\\},$
(3)

with $g$ the considered canonical link function (e.g., identity for continuous
outcome and logit for binary outcome). This can be done using standard
statistical software for fitting a generalized linear regression (with
canonical link) of $Y$ with an intercept (denoted as $\epsilon$) using offset
$g\left(\hat{Y}_{1,i}^{init}\right)$ among participants with $Z=1$. For
example, for a binary endpoint $Y$, one can fit a logistic regression model
(via maximum likelihood) with only an intercept and an offset
$logit\left(\hat{Y}_{1,i}^{init}\right)$. The final estimate for
$\hat{\mu}_{1}$ is then obtained by

$n^{-1}\sum_{i=1}^{n}g^{-1}\left\\{g\left(\hat{Y}_{1,i}^{init}\right)+\hat{\epsilon}\right\\},$

where $\hat{\epsilon}$ is obtained via maximum likelihood. While the primary
attention has been on parametric models in Step 1a, TMLE also allows for the
utilization of sophisticated data-adaptive methods, including machine learning
algorithms. However, this necessitates sample splitting, as outlined in
Section 4.3.

In Appendix A.4, we show that the standard errors of Section 4.1 are valid for
the TMLE, under the same weak conditions as discussed in that section. This
has been conjectured before (Moore and van der Laan, 2009) and recently been
shown for sample splitting (w.r.t. ordinal and time-to-event outcomes) by
Williams et al. (2021), but to the best of our knowledge has not been formally
shown in full generality (e.g., for parametrically estimated propensity scores
as is done in Section 6).

###  4.3 Role of sample splitting

The sparsity assumptions required for the method(s) proposed in the previous
section(s) justify the use of certain data-adaptive approaches such as the
(post-) Lasso or stepwise selection. They can be relaxed by ensuring that the
predictions $\hat{Y}_{1,i}$ and $\hat{Y}_{0,i}$ (and similarly
$\hat{Y}_{1,i}^{init}$ and $\hat{Y}_{0,i}^{init}$ for the targeted maximum
likelihood approach) are estimated from a separate sample. This can be
achieved via sample splitting (see next paragraph; Zheng and van der Laan
(2011); Chernozhukov et al. (2018)). Its use is appealing as it enables the
use of more complex data-adaptive methods, such as machine learning
algorithms.

Sample splitting is performed as follows. First, we split the data into $K$
non-overlapping folds. Let $\mathcal{V}_{1},\dots,\mathcal{V}_{K}$ denote the
partition of the index set $\\{1,\dots,n\\}$ into $K$ non-overlapping subsets.
Next, for the participants belonging to a certain fold $k=1,\dots,K$, the
prediction models (under both treatments) are estimated based on the $K-1$
remaining folds, that is, $\\{1,\dots,n\\}\backslash\mathcal{V}_{k}$ (Step 1).
This might involve more complex data-adaptive methods such as machine learning
algorithms as well as variable selection approaches for parametric models.
Then, for each participant $i$ in fold $k$, predictions $\hat{Y}_{1,k,i}$
(under $Z=1$) and $\hat{Y}_{0,k,i}$ (under $Z=0$) are obtained based on these
models (Step 2). Define
$\hat{\pi}_{k}=\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}$ as
an estimate of $\pi=P(Z=1)$ based on fold $k$, with $|\mathcal{V}_{k}|$ the
number of participants in fold $k$ (i.e., participants with index in
$\mathcal{V}_{k}$). As the prediction unbiasedness condition (1) no longer
holds under sample splitting, an estimate $\hat{\theta}_{k}$ of $\theta$, is
then obtained by taking the mean of

$\displaystyle\frac{Z_{i}}{\hat{\pi}_{k}}\left\\{Y_{i}-\hat{Y}_{1,k,i}\right\\}+\hat{Y}_{1,k,i}-\left[\frac{1-Z_{i}}{1-\hat{\pi}_{k}}\left\\{Y_{i}-\hat{Y}_{0,k,i}\right\\}+\hat{Y}_{0,k,i}\right],$
(4)

over all participants $i$ in fold $k$ (Robins et al., 1994). The final
estimate $\hat{\theta}_{SS}$ (where ‘SS’ stands for ‘sample splitting’) is
then obtained by taking the mean over these $K$ estimates
$\hat{\theta}_{1},\dots,\hat{\theta}_{K}$. This estimator as well as the
variant where the association between $Z$ and $X$ is parametrically estimated
(see Section 6) inherits the appealing properties discussed in the previous
sections.

In addition, $\hat{\theta}_{SS}$ has some other attractive properties under
model misspecification. First, as mentioned in Section 4.1, some form of
efficiency is attainable when the outcome model is misspecified. Concretely,
the estimator based on Equation (4) with minimal variance (over a broad class
of estimators, including the unadjusted estimator; see Appendix A.6) can be
obtained by selecting and fitting the outcome model to have minimal MSE over
that class (see Appendix A.6). This is especially appealing as the approach in
this section allows the use of machine learning algorithms, many of which
allow to choose the MSE as the loss function of interest. Moreover,
surprisingly, when the known randomization probability $\pi$ (e.g., 0.5 for a
1:1 randomized trial) is plugged-in for the propensity score (instead of the
estimated $\hat{\pi}_{k}$’s), $\hat{\theta}_{SS}$ is exactly unbiased, even in
finite samples (see Appendix A.2.2). This does not require correctly specified
models. This also means that, when using estimated propensity scores
$\hat{\pi}_{k}$, any potential bias is solely attributed to the estimation of
these randomization probabilities, which we can reasonably expect to be small
in magnitude. Indeed, note that unadjusted estimators (e.g., mean differences)
likewise exhibit finite sample bias due to the denominator (representing
number of people in each arm) being subject to randomness. The weak
assumptions needed to obtain valid inference with this approach are described
in Appendix A.1. A possible issue of sample splitting is that it may introduce
additional variation, which is a potential concern (see e.g., Guo and Shah,
2023). An explanation and a proof for the validity of this approach under
misspecified working models is given in Appendix A.2.

Sample splitting can also be combined with TMLE (see Appendix A.4.1 for
technical details); commonly referred to as cross-validated targeted minimum-
loss-based estimation (CV-TMLE). Remarkably, CV-TMLE no longer has the
attractive finite-sample-unbiasedness property due to the estimation of the
intercept $\epsilon$ (over all folds). Nevertheless, we still expect at most
small finite sample bias as any potential bias is solely attributed to the
estimation of the randomization probabilities and $\epsilon$.

###  4.4 Avoiding sample splitting when testing the strong null – without
additional assumptions

Interestingly, sample splitting can be avoided — even when using machine
learning predictions — if interest lies solely in testing the (strong) null
hypothesis that $Z$ has no effect on $Y$ (i.e., $Z$ is independent of $Y$). In
particular, one may then estimate $E(Y|Z=1,X)$ and $E(Y|Z=0,X)$ under the null
hypothesis that $Z$ is jointly independent of $Y$ and $X$. Specifically, we
assume $E(Y|Z=1,X)=E(Y|Z=0,X)=E(Y|X)$ by fitting one model in all participants
only including the covariates $X$ (and not the study arm assignment indicator
$Z$). The prediction for participant $i$ ($i=1,\dots,n$) under this model is
denoted by $\hat{Y}_{i}$. Under the assumption that $Z$ is jointly independent
of $Y$ and $X$, $\mu_{1}=E(Y|Z=1)$ can be estimated as
$\hat{\mu}_{1,SN}=\frac{1}{n}\sum_{i=1}^{n}\frac{Z_{i}}{\hat{\pi}}\left\\{Y_{i}-\hat{Y}_{i}\right\\}+\hat{Y}_{i}$.
Here, ‘SN’ stands for ‘strong null’. Note again that there are no restrictions
on how the predictions are obtained, meaning that they are potentially
obtained by machine learning algorithms. Similarly,
$\hat{\mu}_{0,SN}=\frac{1}{n}\sum_{i=1}^{n}\frac{1-Z_{i}}{1-\hat{\pi}}\left\\{Y_{i}-\hat{Y}_{i}\right\\}+\hat{Y}_{i}$,
which uses the same predictions $\hat{Y}_{i}$ as $\hat{\mu}_{1,SN}$.

In Appendix A.2.1, we give a variance estimator and prove the validity of this
approach. Under the strong null, the resulting estimator,
$\hat{\theta}_{SN}=\hat{\mu}_{1,SN}-\hat{\mu}_{0,SN}$, also has the advantage
of being exactly unbiased, even in finite samples as long as the estimated
probability $\hat{\pi}$ is replaced with the known randomization probability
$\pi$ (e.g., 0.5 for a 1:1 randomized trial).

##  5 Simulation Study

Will be added in a later version.

##  6 Enhancing precision: Estimating propensity scores parametrically

Up to this point, we have estimated the propensity score $\pi=P(Z=1|X)$ as the
proportion of treated individuals. Even so, under misspecification of the
outcome model, more efficiency can be attained by estimating the propensity
score as a function of covariates $X$ (despite the fact that it is constant by
design). The results from the previous sections remain valid in that case, so
long as a parametric propensity score model is used (e.g., logistic
regression).

Let us define $P(Z=1|X)=p(X;\beta_{0})$, where $p(X;\beta)$ represents a known
function of an unknown parameter $\beta$, the population value of which is
denoted as $\beta_{0}$. Let $\hat{\beta}$ be an consistent estimator of
$\beta_{0}$ under this model. The data-adaptive covariate-adjustment approach
using canonical generalized linear models (see Section 4) can then be extended
by retaining all steps (i.e., Step 1a, Step 2 and Step 3), but replacing Step
1b by:

  * Step 1b:

Model fitting Fit a weighted generalized linear regression model with a
canonical link (e.g., logistic, linear, …) and intercept via maximum
likelihood that regresses the outcome $Y$ on all covariates selected in Step
1a among the treated participants (i.e., with $Z=1$), with weights
$1/p(X;\hat{\beta})$.

A similar approach can be followed to estimate $\mu_{0}$, but by using weights
$1/\\{1-p(X;\hat{\beta})\\}$ in Step 1b.

Accordingly, the targeted maximum likelihood estimator in Section 4.2 can be
adapted by replacing the updating step in Equation (3) by

$\displaystyle
E(Y|Z=1,X)=g^{-1}\left\\{g\left(\hat{Y}_{1}^{init}\right)+\epsilon\cdot\frac{1}{p(X;\hat{\beta}_{k})}\right\\},$

where $\hat{\beta}_{k}$ is an estimate of $\beta_{0}$ obtained via estimation
within a parametric model (e.g., using logistic regression) in fold $k$.
Specifically, we can fit a generalized linear regression (with canonical link)
of $Y$ without an intercept but with a clever covariate
$\frac{1}{p(X;\hat{\beta}_{k})}$ using offset
$g\left(\hat{Y}_{1}^{init}\right)$ among participants with $Z=1$. For
$\mu_{0}$, the clever covariate can be replaced by
$\frac{1}{1-p(X;\hat{\beta}_{k})}$. The estimation approach using sample
splitting (see Section 4.3), can be extended by replacing $\hat{\pi}_{k}$ by
$p(X;\hat{\beta}_{k})$ in Equation (4).

A proof of the validity of these approaches under possibly misspecified
working models is given in Appendix A.5. As we do not allow variable or model
selection for the propensity score, we need to pre-specify the covariates that
will be adjusted for along with the functional form of the propensity score
model. In order to reduce variance, we suggest the propensity score model
should include variables stratified during randomization (if not a priori
included in the outcome model) and a pre-specified selection of a small number
of key prognostic variables measured at baseline.

##  7 Discussion

It is imperative that the analysis of randomized experiments is based on
statistical procedures that are completely pre-specified (FDA, 2021; EMA,
2015). This requirement is nonetheless difficult to reconcile with the need of
most statistical procedures to be based on correctly specified models. This
has kept trialists away from considering extensive covariate adjustment beyond
simple adjustment for a handful of discrete covariates, in turn leading trial
analyses to have suboptimal power. Given the rapidly growing interest in the
use of data-adaptive methods in statistical inference in general (Van der Laan
et al., 2011; Belloni et al., 2014a; Chernozhukov et al., 2018; Van Lancker et
al., 2021a, b), we wanted to address the conditions under which these methods
can be safely used to estimate (marginal) treatment effects in clinical
trials. In particular, we have shown that (a generalization of) the covariate
adjusted estimator suggested in the recent U.S. Food and Drug Administration
guidance draft on covariate adjustment (FDA, 2021) can be combined with data-
adaptive methods.

The independence between treatment and baseline covariates renders these
covariate-adjusted estimators immune to model misspecification as well as
variable selection uncertainty in randomized clinical trials. It moreover
ensures that valid inference can be obtained, based on simple analytic
expressions, despite these complications. This is attractive because it
enables the use of arbitrarily complex models (e.g., including splines) or
(machine learning) algorithms to be used without distorting the interpretation
or validity of the treatment effect estimate. This flexibility ensures that
the resulting analysis is purely evidence-based and allows covariate
information to be optimally used to gain power (Tsiatis et al., 2008). Our
derivations relied on ultra-sparsity, which we could relax either via sample
splitting or by limiting the number of predictors in the model.

By enabling the use of data-adaptive analysis strategies, while at the same
time guaranteeing a valid inference that is insulated against model
misspecification, the approaches outlined in this article offer the
possibility to incorporate covariate adjustment in an automated manner and on
a large scale. They simplify the difficulty of pre-specifying the statistical
analysis to pre-specification of the full covariate set, along with the
algorithm for variable selection / model building that will be employed. This
form of pre-specification is what we should strive for in the future. Note
that without such pre-specification of the algorithm, there is a general
danger for data fishing, e.g., conducting multiple analyses and reporting only
the one with largest estimated effect, which can lead to bias and / or
inflated $p$-values.

When using covariate-adjusted estimators – even without data-adaptive methods,
different methods exist for determining the required sample size. The first
approach takes a conservative standpoint by assuming that covariate adjustment
will not yield any increase in precision. Any actual precision gain in this
case would then enhance the statistical power. An alternative is to assess the
potential precision improvement achievable through external trial data while
calculating the sample size. This tends to deliver smaller sample sizes but
comes with the risk of underpowering the trial when erroneously projecting
covariate’s prognostic value (Li et al., 2023).

Due to the concern of underestimated standard errors for large sample methods
such as standardization, we suggested multiplying the variance estimators by
the small sample correction factor proposed by Tsiatis et al. (2008). For
data-adaptive approaches, this is, however, conservative. Further research
into refining approximations of the asymptotic distribution might be
advantageous.

Finally, a drawback of covariate-adjusted methods is that they may be prone to
additional missing (baseline covariate) data in the statistical analysis.
However, if covariates have missing values, this does not necessarily pose
problems for the analysis as it is justified to adopt missing indicator and
mean imputation procedures (Zhao and Ding, 2024). These imputed covariates,
along with the missing indicators, can be considered as a separate set of
baseline covariates. As our procedures demand treatment assignment to be
independent of all baseline covariates, and thus also of the imputed
covariates, it is then crucial not to utilize any treatment or outcome
information during the imputation process (except if correct imputation models
are used).

## References

  * Robinson and Jewell [1991] Laurence D Robinson and Nicholas P Jewell.  Some surprising results about covariate adjustment in logistic regression models.  _International Statistical Review/Revue Internationale de Statistique_ , pages 227–240, 1991.
  * Daniel et al. [2021] Rhian Daniel, Jingjing Zhang, and Daniel Farewell.  Making apples from oranges: Comparing noncollapsible effect estimators and their standard errors after adjustment for different covariate sets.  _Biometrical Journal_ , 63(3):528–557, 2021.
  * Tsiatis et al. [2008] Anastasios A Tsiatis, Marie Davidian, Min Zhang, and Xiaomin Lu.  Covariate adjustment for two-sample treatment comparisons in randomized clinical trials: a principled yet flexible approach.  _Statistics in medicine_ , 27(23):4658–4677, 2008\.
  * FDA and EMA [1998] FDA and EMA.  E9 statistical principles for clinical trials. U.S. Food and Drug Administration: CDER/CBER. European Medicines Agency: CPMP/ICH/363/96.  https://www.ema.europa.eu/en/documents/scientific-guideline/ich-e-9-statistical-principles-clinical-trials-step-5$\\_$en.pdf, 1998\.  Last checked: 2021-02-03.
  * EMA [2015] EMA.  Guideline on adjustment for baseline covariates in clinical trials., 2015.  URL www.ema.europa.eu.  Last checked: 2022-05-30.
  * Van der Laan et al. [2011] Mark J Van der Laan, Sherri Rose, et al.  _Targeted learning: causal inference for observational and experimental data_ , volume 10.  Springer, 2011.
  * Wager et al. [2016] Stefan Wager, Wenfei Du, Jonathan Taylor, and Robert J Tibshirani.  High-dimensional regression adjustments in randomized experiments.  _Proceedings of the National Academy of Sciences_ , 113(45):12673–12678, 2016.
  * Wu and Gagnon-Bartsch [2018] Edward Wu and Johann A Gagnon-Bartsch.  The loop estimator: Adjusting for covariates in randomized experiments.  _Evaluation review_ , 42(4):458–488, 2018.
  * Vansteelandt [2021] Stijn Vansteelandt.  Statistical modelling in the age of data science.  _Observational Studies_ , 7(1):217–228, 2021\.
  * Chernozhukov et al. [2018] Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins.  Double/debiased machine learning for treatment and structural parameters.  _The Econometrics Journal_ , 21(1):C1–C68, 01 2018.  ISSN 1368-4221.  doi: 10.1111/ectj.12097.
  * Leeb and Pötscher [2008] Hannes Leeb and Benedikt M Pötscher.  Can one estimate the unconditional distribution of post-model-selection estimators?  _Econometric Theory_ , 24(2):338–376, 2008.
  * Belloni et al. [2014a] Alexandre Belloni, Victor Chernozhukov, and Christian Hansen.  Inference on treatment effects after selection among high-dimensional controls.  _The Review of Economic Studies_ , 81(2):608–650, 2014a.
  * Belloni et al. [2017] Alexandre Belloni, Victor Chernozhukov, Iván Fernández-Val, and Christian Hansen.  Program evaluation and causal inference with high-dimensional data.  _Econometrica_ , 85(1):233–298, 2017.
  * Moore and van der Laan [2009] Kelly L Moore and Mark J van der Laan.  Covariate adjustment in randomized trials with binary outcomes: targeted maximum likelihood estimation.  _Statistics in medicine_ , 28(1):39–64, 2009\.
  * Vermeulen and Vansteelandt [2015] Karel Vermeulen and Stijn Vansteelandt.  Bias-reduced doubly robust estimation.  _Journal of the American Statistical Association_ , 110(511):1024–1036, 2015.
  * Avagyan and Vansteelandt [2022] Vahe Avagyan and Stijn Vansteelandt.  High-dimensional inference for the average treatment effect under model misspecification using penalized bias-reduced double-robust estimation.  _Biostatistics & Epidemiology_, 6(2):221–238, 2022.
  * Benkeser et al. [2017] David Benkeser, Marco Carone, MJ Van Der Laan, and Peter B Gilbert.  Doubly robust nonparametric inference on the average treatment effect.  _Biometrika_ , 104(4):863–880, 2017.
  * Dukes et al. [2021] Oliver Dukes, Stijn Vansteelandt, and David Whitney.  On doubly robust inference for double machine learning.  _arXiv preprint arXiv:2107.06124_ , 2021.
  * Williams et al. [2021] Nicholas Williams, Michael Rosenblum, and Iván Díaz.  Optimizing precision and power by machine learning in randomized trials, with an application to covid-19, 2021.  URL https://arxiv.org/abs/2109.04294.
  * Bannick et al. [2023] Marlena S Bannick, Jun Shao, Jingyi Liu, Yu Du, Yanyao Yi, and Ting Ye.  A general form of covariate adjustment in randomized clinical trials.  _arXiv preprint arXiv:2306.10213_ , 2023.
  * FDA [2021] FDA.  Adjusting for Covariates in Randomized Clinical Trials for Drugs and Biological Products. Guidance for Industry.  https://www.fda.gov/media/148910/download, 2021.  Last checked: 2021-10-20.
  * Benkeser et al. [2020] David Benkeser, Iván Díaz, Alex Luedtke, Jodi Segal, Daniel Scharfstein, and Michael Rosenblum.  Improving precision and power in randomized trials for COVID-19 treatments using covariate adjustment, for binary, ordinal, and time-to-event outcomes.  _Biometrics_ , 77:1467–1481, 2020.
  * Robins et al. [1994] James M Robins, Andrea Rotnitzky, and Lue Ping Zhao.  Estimation of regression coefficients when some regressors are not always observed.  _Journal of the American statistical Association_ , 89(427):846–866, 1994.
  * Yang and Tsiatis [2001] Li Yang and Anastasios A Tsiatis.  Efficiency study of estimators for a treatment effect in a pretest–posttest trial.  _The American Statistician_ , 55(4):314–321, 2001\.
  * Tsiatis [2006] Anastasios A Tsiatis.  _Semiparametric theory and missing data_.  Springer, 2006.
  * Rubin and van der Laan [2008] Daniel B Rubin and Mark J van der Laan.  Empirical efficiency maximization: Improved locally efficient covariate adjustment in randomized experiments and survival analysis.  _The International Journal of Biostatistics_ , 4(1), 2008\.
  * Cao et al. [2009] Weihua Cao, Anastasios A Tsiatis, and Marie Davidian.  Improving efficiency and robustness of the doubly robust estimator for a population mean with incomplete data.  _Biometrika_ , 96(3):723–734, 2009.
  * Gruber and van der Laan [2012] Susan Gruber and Mark J van der Laan.  Targeted minimum loss based estimator that outperforms a given estimator.  _The International Journal of Biostatistics_ , 8(1), 2012\.
  * Rosenblum and Van Der Laan [2010] Michael Rosenblum and Mark J Van Der Laan.  Simple, efficient estimators of treatment effects in randomized trials using generalized linear models to leverage baseline variables.  _The international journal of biostatistics_ , 6(1), 2010\.
  * Zheng and van der Laan [2011] Wenjing Zheng and Mark J van der Laan.  Cross-validated targeted minimum-loss-based estimation.  In _Targeted Learning_ , pages 459–474. Springer, 2011.
  * Guo and Shah [2023] F Richard Guo and Rajen D Shah.  Rank-transformed subsampling: inference for multiple data splitting and exchangeable p-values.  _arXiv preprint arXiv:2301.02739_ , 2023.
  * Van Lancker et al. [2021a] Kelly Van Lancker, Oliver Dukes, and Stijn Vansteelandt.  Ensuring valid inference for hazard ratios after variable selection, 2021a.  URL https://arxiv.org/abs/2112.00172.
  * Van Lancker et al. [2021b] Kelly Van Lancker, Oliver Dukes, and Stijn Vansteelandt.  Principled selection of baseline covariates to account for censoring in randomized trials with a survival endpoint.  _Statistics in Medicine_ , 40(18):4108–4121, 2021b.
  * Li et al. [2023] Xiudi Li, Sijia Li, and Alex Luedtke.  Estimating the efficiency gain of covariate-adjusted analyses in future clinical trials using external data.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 85(2):356–377, 2023.
  * Zhao and Ding [2024] Anqi Zhao and Peng Ding.  To adjust or not to adjust? estimating the average treatment effect in randomized experiments with missing covariates.  _Journal of the American Statistical Association_ , 119(545):450–460, 2024.
  * Farrell [2015] Max H Farrell.  Robust inference on average treatment effects with possibly more covariates than observations.  _Journal of Econometrics_ , 189(1):1–23, 2015\.
  * Belloni et al. [2014b] Alexandre Belloni, Victor Chernozhukov, and Christian Hansen.  High-dimensional methods and inference on structural and treatment effects.  _Journal of Economic Perspectives_ , 28(2):29–50, 2014b.
  * Bradic et al. [2019] Jelena Bradic, Stefan Wager, and Yinchu Zhu.  Sparsity double robust inference of average treatment effects.  _arXiv preprint arXiv:1905.00744_ , 2019.
  * Belloni and Chernozhukov [2013] Alexandre Belloni and Victor Chernozhukov.  Least squares after model selection in high-dimensional sparse models.  _Bernoulli_ , 19(2):521–547, 05 2013.  doi: 10.3150/11-BEJ410.
  * Billingsley [1999] Patrick Billingsley.  _Convergence of probability measures_.  Wiley Series in Probability and Statistics: Probability and Statistics. John Wiley & Sons Inc., New York, second edition, 1999.  ISBN 0-471-19745-9.  A Wiley-Interscience Publication.
  * de la Peña et al. [2009] Victor de la Peña, Tze Lai, and Qi-Man Shao.  _Self-normalized processes. Limit theory and statistical applications_.  Probability and Its Applications. Springer Berlin Heidelberg, 2009.  ISBN 9783540856368.  doi: 10.1007/978-3-540-85636-8.
  * Belloni et al. [2012] Alexandre Belloni, Daniel Chen, Victor Chernozhukov, and Christian Hansen.  Sparse models and methods for optimal instruments with an application to eminent domain.  _Econometrica_ , 80(6):2369–2429, 2012.
  * Krantz and Parks [2002] Steven G Krantz and Harold R Parks.  _A primer of real analytic functions_.  Springer Science & Business Media, 2002.
  * Newey and McFadden [1994] Whitney K Newey and Daniel McFadden.  Large sample estimation and hypothesis testing.  _Handbook of econometrics_ , 4:2111–2245, 1994.

##  Appendix A Theoretical derivations

###  A.1 Assumptions

We denote the regression coefficients obtained by the data-adaptive procedure
in Section 4 as
$\tilde{\boldsymbol{\gamma}}=(\tilde{\gamma}_{0},\tilde{\gamma}_{1})^{\prime}$
for estimation of $\mu_{1}$, and
$\tilde{\boldsymbol{\eta}}=(\tilde{\eta}_{0},\tilde{\eta}_{1})^{\prime}$ for
estimation of $\mu_{0}$. To allow for model misspecification, we introduce the
limiting value of the estimated value $\tilde{\boldsymbol{\gamma}}$ as
$\boldsymbol{\gamma}^{*}$ and the limiting value of the estimated value
$\tilde{\boldsymbol{\eta}}$ as $\boldsymbol{\eta}^{*}$. Let $p$ denote the
number of considered terms in the model (including possible interactions,
transformations, higher order terms, …). Define the sets of terms that have
coefficients that are ‘truly’ non-zero (i.e., active set) as
$S_{\gamma}=\text{support}(\boldsymbol{\gamma}^{*})$ and
$S_{\eta}=\text{support}(\boldsymbol{\eta}^{*})$, and define
$s_{\gamma}=|S_{\gamma}|$ and $s_{\eta}=|S_{\eta}|$.

######  Assumption 1 (Ultra-sparsity)

$n^{-1/2}(s_{\gamma})\log(p\lor n)=o(1)$ $n^{-1/2}(s_{\eta})\log(p\lor
n)=o(1)$

This assumption frequently finds application in the literature [see e.g.,
Farrell, 2015, Belloni et al., 2014b]. Although in general this might be a
strong assumption that may be difficult to validate in practice [Bradic et
al., 2019], in our setting this seems rather a weak assumption. The rationale
behind this is our allowance for a misspecified outcome working model,
enabling us to keep the number of covariates ($p$) limited.

We will also make use of the following rate conditions:

######  Assumption 2 (Rates on error of estimated coefficients)

  1. 1.

$||\boldsymbol{\gamma}^{*}-\tilde{\boldsymbol{\gamma}}||_{1}=O_{P_{n}}\left(s_{\gamma}\sqrt{\log(p\lor
n)/n}\right)$

  2. 2.

$||\boldsymbol{\gamma}^{*}-\tilde{\boldsymbol{\gamma}}||_{2}=O_{P_{n}}\left(\sqrt{s_{\gamma}\log(p\lor
n)/n}\right)$

  3. 3.

$||\boldsymbol{\eta}^{*}-\tilde{\boldsymbol{\eta}}||_{1}=O_{P_{n}}\left(s_{\eta}\sqrt{\log(p\lor
n)/n}\right)$

  4. 4.

$||\boldsymbol{\eta}^{*}-\tilde{\boldsymbol{\eta}}||_{2}=O_{P_{n}}\left(\sqrt{s_{\eta}\log(p\lor
n)/n}\right)$

where $\\{P_{n}\\}$ refers to a sequence of data-generating processes that
obey, for each $n$, the assumptions above. We refer the reader to Belloni and
Chernozhukov [2013] and Belloni et al. [2017] who show that these rates hold
for the Lasso estimator in the context of linear and logistic regression.

Besides the aforementioned sparsity conditions, we only need the following
weak assumptions to ensure that the data-adaptive approach in Section 4 is
asymptotically insensitive to selection mistakes (i.e., preserves Type I
error).

######  Assumption 3 (Positivity assumption)

$1>1-\rho\geq P_{n}(Z=1)\geq\rho>0,$

with $\rho$ a constant.

######  Assumption 4 (Convergence in probability of predictions)

The imputations $h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})$ and
$h_{0}(X_{i},\tilde{\boldsymbol{\eta}})$ converge in probability to (possibly
misspecified) limits $h_{1}(X_{i};\boldsymbol{\gamma}^{*})$ and
$h_{0}^{*}(X_{i};\boldsymbol{\eta}^{*})$, respectively, in the sense that the
$L_{2}$-norms

$\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left\\{h_{1}(X_{i};\boldsymbol{\gamma}^{*})-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}^{2}}=o_{P_{n}}(1),$

and

$\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left\\{h_{0}(X_{i};\boldsymbol{\eta}^{*})-h_{0}(X_{i};\tilde{\boldsymbol{\eta}})\right\\}^{2}}=o_{P_{n}}(1).$

Here, $o_{P_{n}}(1)$ denotes a stochastic variable that converges in
probability to zero as $n\rightarrow\infty$.

We make similar assumptions for
$\hat{Y}_{1,i}^{init}=h_{1}(X_{i};\tilde{\tilde{\boldsymbol{\gamma}}})$ and
$\hat{Y}_{0,i}^{init}=h_{0}(X_{i},\tilde{\tilde{\boldsymbol{\eta}}})$.

We denote the imputations in Section 4.3, $\hat{Y}_{1,k,i}$ (under $Z=1$) and
$\hat{Y}_{0,k,i}$ (under $Z=0$), which can be based on machine learning
algorithms, as $\hat{h}_{1,k}(X_{i})$ and $\hat{h}_{0,k}(X_{i})$ (for
$k=1,\dots,K$), respectively. The following assumption focuses on these
imputations $\hat{h}_{1,k}(X_{i})$ and $\hat{h}_{0,k}(X_{i})$.

###### Assumption 4′ (Convergence in probability of predictions for sample
splitting)

The imputations $\hat{h}_{1,k}(X_{i})$ and $\hat{h}_{0,k}(X_{i})$ (for
$k=1,\dots,K$) converge in probability to (possibly misspecified) limits
$h_{1}^{*}(X_{i})$ and $h_{0}^{*}(X_{i})$, respectively, in the sense that the
$L_{2}$-norms

$\sqrt{\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}^{2}}$

and

$\sqrt{\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{h_{0}^{*}(X_{i})-\hat{h}_{0,k}(X_{i})\right\\}^{2}}$

converge in probability to zero.

######  Assumption 5 (Uniform integrability without sample splitting)

$\frac{1}{n}\sum_{i=1}^{n}\left\\{h^{*}_{1}(X_{i})-\hat{h}_{1}(X_{i})\right\\}^{2}$

is uniformly integrable.
A collection of random variables $\\{G_{n}:n\geq n_{0}>0\\}$ with
$G_{n}=\frac{1}{|\mathcal{V}_{k}|}\sum_{i=\in\mathcal{V}_{k}}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1}(X_{i})\right\\}^{2}$
is uniformly integrable if $sup_{n\geq
n_{0}}E\left\\{|G_{n}|I(|G_{n}|>g)\right\\}=0$, as $g\rightarrow\infty$. Here,
$n_{0}$ can be seen as a minimal (total) sample size.

We make a similar assumption for the imputations $\hat{h}_{1,k}(X_{i})$ and
$\hat{h}_{0,k}(X_{i})$:

###### Assumption 5′ (Uniform integrability with sample splitting)

$\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}^{2}$

is uniformly integrable.

In the next section, we first explain why the data-adaptive approach with
sample splitting in Section 4.3 gives valid inference. Only in Appendix A.3,
we will give a proof to show that we can avoid cross-fitting when a parametric
outcome regression is fitted, possibly data-adaptive (using consistent
variable selection, regardless of loss function).

###  A.2 Data-adaptive covariate adjustment with sample splitting

We denote the imputations in Section 4.3 (possibly obtained by machine
learning), $\hat{Y}_{1,k,i}$ (under $Z=1$) and $\hat{Y}_{0,k,i}$ (under
$Z=0$), as $\hat{h}_{1,k}(X_{i})$ and $\hat{h}_{0,k}(X_{i})$ (for
$k=1,\dots,K$), respectively. In this appendix we first show that the
estimator based on Equation (4),

$\hat{\theta}_{k}=\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\frac{Z_{i}}{\hat{\pi}_{k}}\left\\{Y_{i}-\hat{h}_{1,k}(X_{i})\right\\}+\hat{h}_{1,k}(X_{i})-\left[\frac{1-Z_{i}}{1-\hat{\pi}_{k}}\left\\{Y_{i}-\hat{h}_{0,k}(X_{i})\right\\}+\hat{h}_{0,k}(X_{i})\right],$

is asymptotically linear, assuming the predictions $\hat{h}_{1,k}(X_{i})$ and
$\hat{h}_{0,k}(X_{i})$ are obtained from a separate sample. This can be
accommodated by sample splitting [Zheng and van der Laan, 2011, Chernozhukov
et al., 2018]. In particular, for the participants belonging to a certain fold
$k=1,\dots,K$, the predictions $\hat{h}_{1,k}(X)$ and $\hat{h}_{0,k}(X)$ are
estimated based on the $K-1$ remaining folds, that is,
$\\{1,\dots,n\\}\backslash\mathcal{V}_{k}$.

For simplicity, we focus on estimation of $\mu_{1}=E(Y|Z=1)$, which is
estimated by
$\hat{\mu}_{1,k}=\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\frac{Z_{i}}{\hat{\pi}_{k}}\left\\{Y_{i}-\hat{h}_{1,k}(X_{i})\right\\}+\hat{h}_{1,k}(X_{i})$
in fold $k$. We then have

$\displaystyle\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{\hat{\pi}_{k}}\left\\{Y_{i}-\hat{h}_{1,k}(X_{i})\right\\}+\hat{h}_{1,k}(X_{i})\right]$
$\displaystyle=\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-\hat{h}_{1,k}(X_{i})\right\\}+\hat{h}_{1,k}(X_{i})\right]$
(5)
$\displaystyle+\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{Y_{i}-\hat{h}_{1,k}(X_{i})\right\\}\left(\frac{1}{\hat{\pi}_{k}}-\frac{1}{\pi}\right)\sqrt{|\mathcal{V}_{k}|}.$
(6)

For term (6), we have

$\displaystyle\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}$
$\displaystyle
Z_{i}\left\\{Y_{i}-\hat{h}_{1,k}(X_{i})\right\\}\left(\frac{1}{\hat{\pi}_{k}}-\frac{1}{\pi}\right)\sqrt{|\mathcal{V}_{k}|}$
$\displaystyle=\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{Y_{i}-\hat{h}_{1,k}(X_{i})+h_{1}^{*}(X_{i})-h_{1}^{*}(X_{i})\right\\}\left(\frac{1}{\hat{\pi}_{k}}-\frac{1}{\pi}\right)\sqrt{|\mathcal{V}_{k}|}$
$\displaystyle=\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left[Y_{i}-h_{1}^{*}(X_{i})+\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}\right]\left(\frac{1}{\hat{\pi}_{k}}-\frac{1}{\pi}\right)\sqrt{|\mathcal{V}_{k}|}$
$\displaystyle=\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{Y_{i}-h_{1}^{*}(X_{i})\right\\}\left(\frac{1}{\hat{\pi}_{k}}-\frac{1}{\pi}\right)\sqrt{|\mathcal{V}_{k}|}$
(7)
$\displaystyle+\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}\left(\frac{1}{\hat{\pi}_{k}}-\frac{1}{\pi}\right)\sqrt{|\mathcal{V}_{k}|}.$
(8)

Then, following a Taylor expansion for term (7) with respect to $\pi$, we have

$\displaystyle\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{Y_{i}-h_{1}^{*}(X_{i})\right\\}\left(\frac{1}{\hat{\pi}_{k}}-\frac{1}{\pi}\right)\sqrt{|\mathcal{V}_{k}|}=-\frac{1}{\pi^{2}}(\hat{\pi}_{k}-\pi)\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{Y_{i}-h_{1}^{*}(X_{i})\right\\}\sqrt{|\mathcal{V}_{k}|}+o_{P_{n}}(1),$

where the $o_{P_{n}}(1)$ term follows from the fact that
$|\hat{\pi}_{k}-\pi|^{2}=O_{P_{n}}(n^{-1})$. Note that we use “$n$” instead of
$|\mathcal{V}_{k}|$, which is justified as both quantities are proportional.
By the weak law of large numbers, this equals

$\displaystyle\sqrt{|\mathcal{V}_{k}|}\left\\{-\frac{1}{\pi^{2}}(\hat{\pi}_{k}-\pi)\right\\}E[Z\left\\{Y-h_{1}^{*}(X)\right\\}]+o_{P_{n}}(1).$
(9)

Here, the first term is $O_{P_{n}}(1)$ when
$E[Z\left\\{Y-h_{1}^{*}(X)\right\\}]$ is bounded as
$|\hat{\pi}_{k}-\pi|=O_{P_{n}}(n^{-1/2})$. The boundedness of
$E[Z\left\\{Y-h_{1}^{*}(X)\right\\}]$ is often generally satisfied; e.g., for
learners that provide predictions in range of the data.

Next, applying a Taylor expansion as for term (7), term (8) yields

$\displaystyle\sqrt{|\mathcal{V}_{k}|}$
$\displaystyle\left|\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}\left(\frac{1}{\hat{\pi}_{k}}-\frac{1}{\pi}\right)\right|$
$\displaystyle\leq\sqrt{|\mathcal{V}_{k}|}\rho^{-2}\left|\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[Z_{i}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}\right](\hat{\pi}_{k}-\pi)\right|+o_{P_{n}}(1)$
$\displaystyle\leq\sqrt{|\mathcal{V}_{k}|}\rho^{-2}\sqrt{\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[Z_{i}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}\right]^{2}}|\hat{\pi}_{k}-\pi|+o_{P_{n}}(1)$
$\displaystyle\leq\sqrt{|\mathcal{V}_{k}|}\rho^{-2}\sqrt{\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}^{2}}|\hat{\pi}_{k}-\pi|+o_{P_{n}}(1)$
$\displaystyle=\sqrt{|\mathcal{V}_{k}|}\rho^{-2}o_{P_{n}}(1)O_{P_{n}}(n^{-1/2})+o_{P_{n}}(1)$
$\displaystyle=o_{P_{n}}(1),$

where the second inequality follows from the Cauchy-Schwarz inequality, the
third inequality from the fact that $Z$ is bounded (by 1) and the first
equality from Assumption 4′ (Convergence in probability of predictions for
sample splitting). We conjecture that this term can probably be handled under
weaker assumptions.

In order to handle term (5), we define
$\hat{\mu}_{1,k}^{(\pi)}=\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-\hat{h}_{1,k}(X_{i})\right\\}+\hat{h}_{1,k}(X_{i})\right]$.
Then, it holds that

$\displaystyle\sqrt{|\mathcal{V}_{k}|}(\hat{\mu}_{1,k}^{(\pi)}-\mu_{1})$
$\displaystyle=\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}^{*}(X_{i})\right\\}+h_{1}^{*}(X_{i})-\mu_{1}\right]$
$\displaystyle+\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}\right].$
(10)

Assuming the predictions $\hat{h}_{1,k}(X_{i})$ for fold $k$ are obtained from
a “separate” sample
$\mathbf{O}^{\bar{k}}=\\{\mathbf{O}_{j}:j\in\\{1,\dots,n\\}\backslash\mathcal{V}_{k}\\}$
(i.e., independent from the sample
$\mathbf{O}^{k}=\\{\mathbf{O}_{j}:j\in\mathcal{V}_{k}\\}$ which is used for
calculation of $\hat{\mu}_{1,k}$), we can reason conditional on the sample
$\mathbf{O}^{\bar{k}}$ to handle term (10). First, note that term (10) has
mean zero as

$\displaystyle
E\left(\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}\right]\right)$
$\displaystyle=E\left\\{E\left(\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}\right]\bigg{|}\mathbf{O}^{\bar{k}}\right)\right\\}$
$\displaystyle=E\left(E\left[\left(\frac{Z}{\pi}-1\right)\left\\{h_{1}^{*}(X)-\hat{h}_{1,k}(X)\right\\}\bigg{|}\mathbf{O}^{\bar{k}}\right]\right)$
$\displaystyle=E\left(\frac{Z}{\pi}-1\right)E\left[E\left\\{h_{1}^{*}(X)-\hat{h}_{1,k}(X)\bigg{|}\mathbf{O}^{\bar{k}}\right\\}\right]$
$\displaystyle=0,$

where the third equation follows from the randomization (i.e., $Z$ being
independent of $X$) and the last equation from
$E\left(\frac{Z}{\pi}-1\right)=0$. By the law of total variance, it follows
from the above that

$\displaystyle Var$
$\displaystyle\left(\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}\right]\right)$
$\displaystyle=E\left\\{Var\left(\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}\right]\bigg{|}\mathbf{O}^{\bar{k}}\right)\right\\}$
$\displaystyle+Var\left\\{E\left(\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}\right]\bigg{|}\mathbf{O}^{\bar{k}}\right)\right\\}$
$\displaystyle=E\left\\{Var\left(\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}\right]\bigg{|}\mathbf{O}^{\bar{k}}\right)\right\\}.$

By randomization, this can be rewritten as

$\displaystyle\frac{1}{|\mathcal{V}_{k}|^{2}}E\left(\sum_{i\in\mathcal{V}_{k}}Var\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}\bigg{|}\mathbf{O}^{\bar{k}}\right]\right)$
$\displaystyle=\frac{1}{|\mathcal{V}_{k}|^{2}}E\left(\sum_{i\in\mathcal{V}_{k}}E\left[\left(\frac{Z_{i}}{\pi}-1\right)^{2}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}^{2}\bigg{|}\mathbf{O}^{\bar{k}}\right]\right)$
$\displaystyle=\frac{1}{|\mathcal{V}_{k}|^{2}}E\left\\{\left(\frac{Z}{\pi}-1\right)^{2}\right\\}E\left(\sum_{i\in\mathcal{V}_{k}}E\left[\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}^{2}\bigg{|}\mathbf{O}^{\bar{k}}\right]\right)$
$\displaystyle=\frac{1}{|\mathcal{V}_{k}|}E\left\\{\left(\frac{Z}{\pi}-1\right)^{2}\right\\}E\left[\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}^{2}\right]$

where $E\left\\{\left(\frac{Z}{\pi}-1\right)^{2}\right\\}$ is bounded and
$E\left[\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}^{2}\right]$
converges to zero assuming that
$\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}^{2}$
is uniformly integrable (Assumption 5′ (Uniform integrability with sample
splitting)). The latter follows from Theorem 3.5 in Billingsley [1999] which
states that convergence in distribution (which is implied by convergence in
probability) implies convergence of the means under extra regularity
conditions, that is, under uniform integrability. Specifically, assuming that
$\frac{1}{|\mathcal{V}_{k}|}\sum_{i=\in\mathcal{V}_{k}}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}^{2}$
is uniformly integrable, it follows from Theorem 3.5 in Billingsley [1999]
that
$\sqrt{\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}^{2}}=o_{P_{n}}(1)$
(see Assumption 4′ (Convergence in probability of predictions for sample
splitting)) implies that
$E\left[\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}^{2}\right]$
converges to zero. Application of Chebyshev’s inequality then shows that term
(10) converges in probability to zero: for any constant $l>0$,

$P\left(\left|\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}\right]\right|\geq
l\right)\leq\frac{E\left\\{\left(\frac{Z}{\pi}-1\right)^{2}\right\\}E\left[\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}^{2}\right]}{l^{2}},$

where the right hand side converges to zero.

Consequently,

$\displaystyle\sqrt{|\mathcal{V}_{k}|}(\hat{\mu}_{1,k}^{(\pi)}-\mu_{1})$
$\displaystyle=\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}^{*}(X_{i})\right\\}+h_{1}^{*}(X_{i})-\mu_{1}\right]+o_{P_{n}}(1),$

and thus

$\displaystyle\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-\hat{h}_{1,k}(X_{i})\right\\}+\hat{h}_{1,k}(X_{i})\right]=\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}^{*}(X_{i})\right\\}+h_{1}^{*}(X_{i})\right]+o_{P_{n}}(1).$

Finally, we have

$\displaystyle\frac{1}{\sqrt{|\mathcal{V}_{k}|}}$
$\displaystyle\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{\hat{\pi}_{k}}\left\\{Y_{i}-\hat{h}_{1,k}(X_{i})\right\\}+\hat{h}_{1,k}(X_{i})\right]$
$\displaystyle=\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left(\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}^{*}(X_{i})\right\\}+h_{1}^{*}(X_{i})-E\left[\frac{Z}{\pi^{2}}\left\\{Y-h_{1}^{*}(X)\right\\}\right](Z_{i}-\pi)\right)+o_{P_{n}}(1),$

which leads to

$\displaystyle\begin{split}&\sqrt{|\mathcal{V}_{k}|}(\hat{\mu}_{1,k}-\mu_{1})\\\
&=\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left(\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}^{*}(X_{i})\right\\}+h_{1}^{*}(X_{i})-E\left[\frac{Z}{\pi^{2}}\left\\{Y-h_{1}^{*}(X)\right\\}\right](Z_{i}-\pi)-\mu_{1}\right)+o_{P_{n}}(1).\end{split}$
(11)

A similar reasoning can be followed to prove that

$\displaystyle\begin{split}&\sqrt{|\mathcal{V}_{k}|}(\hat{\mu}_{0,k}-\mu_{0})\\\
&=\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left(\frac{1-Z_{i}}{1-\pi}\left\\{Y_{i}-h_{0}^{*}(X_{i})\right\\}+h_{0}^{*}(X_{i})+E\left[\frac{1-Z}{(1-\pi)^{2}}\left\\{Y-h_{0}^{*}(X)\right\\}\right](Z_{i}-\pi)-\mu_{0}\right)+o_{P_{n}}(1),\end{split}$
(12)

where
$\hat{\mu}_{0,k}=\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\frac{1-Z_{i}}{1-\hat{\pi}_{k}}\left\\{Y_{i}-\hat{h}_{0,k}(X_{i})\right\\}+\hat{h}_{0,k}(X_{i})$.

The final estimate $\hat{\mu}_{1}$ for $\mu_{1}$ is then obtained by taking
the mean over these $K$ estimates $\hat{\mu}_{1,1},\dots,\hat{\mu}_{1,K}$.
Then,

$\displaystyle\sqrt{n}(\hat{\mu}_{1,SS}-\mu_{1})$
$\displaystyle=\sqrt{n}\left(\frac{1}{K}\sum_{k=1}^{K}\hat{\mu}_{1,k}-\mu_{1}\right)$
$\displaystyle=\sqrt{n}\frac{1}{K}\sum_{k=1}^{K}(\hat{\mu}_{1,k}-\mu_{1})$
$\displaystyle=\sqrt{n}\frac{1}{K}\sum_{k=1}^{K}\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sqrt{|\mathcal{V}_{k}|}(\hat{\mu}_{1,k}-\mu_{1})$
$\displaystyle=\sqrt{n}\frac{1}{K}\sum_{k=1}^{K}\frac{1}{\sqrt{|\mathcal{V}_{k}|}}$
$\displaystyle\hskip
7.11317pt\left\\{\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left(\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}^{*}(X_{i})\right\\}+h_{1}^{*}(X_{i})-E\left[\frac{Z}{\pi^{2}}\left\\{Y-h_{1}^{*}(X)\right\\}\right](Z_{i}-\pi)-\mu_{1}\right)+o_{P_{n}}(1)\right\\}$
$\displaystyle=\sqrt{n}\frac{1}{K}\sum_{k=1}^{K}\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left(\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}^{*}(X_{i})\right\\}+h_{1}^{*}(X_{i})-E\left[\frac{Z}{\pi^{2}}\left\\{Y-h_{1}^{*}(X)\right\\}\right](Z_{i}-\pi)-\mu_{1}\right)$
$\displaystyle+o_{P_{n}}(1).$

Assuming that $n$ is a multiple of $K$, this can be simplified as

$\displaystyle\sqrt{n}(\hat{\mu}_{1,SS}-\mu_{1})$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{k=1}^{K}\sum_{i\in\mathcal{V}_{k}}\left(\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}^{*}(X_{i})\right\\}+h_{1}^{*}(X_{i})-E\left[\frac{Z}{\pi^{2}}\left\\{Y-h_{1}^{*}(X)\right\\}\right](Z_{i}-\pi)-\mu_{1}\right)+o_{P_{n}}(1).$

Following a similar reasoning for $\hat{\mu}_{0,SS}$, we have (for the general
case)

$\displaystyle\sqrt{n}(\hat{\theta}_{SS}-\theta)$
$\displaystyle=\sqrt{n}\frac{1}{K}\sum_{k=1}^{K}\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{\left(\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}^{*}(X_{i})\right\\}+h_{1}^{*}(X_{i})-E\left[\frac{Z}{\pi^{2}}\left\\{Y-h_{1}^{*}(X)\right\\}\right](Z_{i}-\pi)-\mu_{1}\right)\right.$
$\displaystyle-\left.\left(\frac{1-Z_{i}}{1-\pi}\left\\{Y_{i}-h_{0}^{*}(X_{i})\right\\}+h_{0}^{*}(X_{i})+E\left[\frac{1-Z}{(1-\pi)^{2}}\left\\{Y-h_{0}^{*}(X)\right\\}\right](Z_{i}-\pi)-\mu_{0}\right)\right\\}+o_{P_{n}}(1).$

The asymptotic variance of
$\hat{\theta}_{SS}=\hat{\mu}_{1,SS}-\hat{\mu}_{0,SS}$ can be calculated as
$1/n$ times the sample variance of the values

$\displaystyle\left\\{\left(\frac{Z_{i}}{\hat{\pi}_{k(i)}}\left\\{Y_{i}-\hat{h}_{1,k(i)}(X_{i})\right\\}+\hat{h}_{1,k(i)}(X_{i})-\frac{1}{|\mathcal{V}_{k(i)}|}\sum_{j\in\mathcal{V}_{k(i)}}\left[\frac{Z_{j}}{\hat{\pi}_{k(i)}^{2}}\left\\{Y_{j}-\hat{h}_{1,k(i)}(X_{j})\right\\}\right](Z_{i}-\hat{\pi}_{k(i)})\right)\right.$
$\displaystyle-\left.\left(\frac{1-Z_{i}}{1-\hat{\pi}_{k(i)}}\left\\{Y_{i}-\hat{h}_{0,k(i)}(X_{i})\right\\}+\hat{h}_{0,k(i)}(X_{i})+\frac{1}{|\mathcal{V}_{k(i)}|}\sum_{j\in\mathcal{V}_{k(i)}}\left[\frac{1-Z_{j}}{(1-\hat{\pi}_{k(i)})^{2}}\left\\{Y_{j}-\hat{h}_{0,k(i)}(X_{j})\right\\}\right](Z_{i}-\hat{\pi}_{k(i)})\right)\right\\},$

where $k(i)$ indicates to which fold participant $i$ belongs. It may seem to
the reader that one need to weight for the unequal fold sizes, but this has no
asymptotic impact.

The uncertainty in the randomization probabilities is taken into account via
the terms

$\frac{1}{|\mathcal{V}_{k(i)}|}\sum_{j\in\mathcal{V}_{k(i)}}\left[\frac{Z_{j}}{\hat{\pi}_{k(i)}^{2}}\left\\{Y_{j}-\hat{h}_{1,k(i)}(X_{j})\right\\}\right](Z_{i}-\hat{\pi}_{k(i)}),$

and

$\frac{1}{|\mathcal{V}_{k(i)}|}\sum_{j\in\mathcal{V}_{k(i)}}\left[\frac{1-Z_{j}}{(1-\hat{\pi}_{k(i)})^{2}}\left\\{Y_{j}-\hat{h}_{0,k(i)}(X_{j})\right\\}\right](Z_{i}-\hat{\pi}_{k(i)}).$

Note that we do not have such a term for the outcome working models as term
(10) converges to zero (due to the randomization probability being a
constant). In Appendix A.4.1, we prove that TMLE would not need to account for
these terms. For time-to-event and ordinal outcomes, Williams et al. [2021]
have already shown this result, which can be extended straightforwardly to
binary and continuous outcomes. In addition, when using the true randomization
probability $\pi$ (see Appendix A.2.2), even the above error terms for the
estimation of the randomization probability can be removed.

####  A.2.1 Avoiding sample splitting when testing the strong null (without
additional assumptions)

In this section, we show how sample splitting can be avoided if interest lies
solely in testing the (strong) null hypothesis, that is, $Z$ has no effect on
$Y$ (i.e., $Z$ is independent of $Y$). As discussed in Section 4.4, we should
then estimate $E(Y|Z=1,X)$ and $E(Y|Z=0,X)$ under the null hypothesis that $Z$
is jointly independent of $Y$ and $X$. Specifically, we assume
$E(Y|Z=1,X)=E(Y|Z=0,X)=E(Y|X)$ by fitting one model $h(X)$ in all participants
only including the covariates $X$ (and not the treatment indicator $Z$).

Under the assumption that $Z$ is independent of $Y$ and $X$,
$\mu_{1}=E(Y|Z=1)$ can be estimated as
$\hat{\mu}_{1,SN}=\frac{1}{n}\sum_{i=1}^{n}\frac{Z_{i}}{\hat{\pi}}\left\\{Y_{i}-\hat{h}(X_{i})\right\\}+\hat{h}(X_{i})$.
Here, ‘SN’ stands for ‘strong null’. Note that there are no restrictions –
besides blinding – on how the predictions are obtained, meaning that they are
potentially obtained by machine learning algorithms. Following the same
reasoning as in Appendix A.2,

$\displaystyle\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\hat{\pi}}\left\\{Y_{i}-\hat{h}(X_{i})\right\\}+\hat{h}(X_{i})\right]$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-\hat{h}(X_{i})\right\\}+\hat{h}(X_{i})\right]$
(13)
$\displaystyle+\sqrt{n}\left\\{-\frac{1}{\pi^{2}}(\hat{\pi}-\pi)\right\\}E[Z\left\\{Y-h^{*}(X)\right\\}]$
$\displaystyle+o_{P_{n}}(1),$

where $h^{*}(X)$ is the probability limit of $\hat{h}(X)$. Then, in order to
handle term (13), we define
$\hat{\mu}_{1,SN}^{(\pi)}=\frac{1}{n}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-\hat{h}(X_{i})\right\\}+\hat{h}(X_{i})\right]$.
Then, it holds that

$\displaystyle\sqrt{n}(\hat{\mu}_{1,SN}^{(\pi)}-\mu_{1})$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h^{*}(X_{i})\right\\}+h^{*}(X_{i})-\mu_{1}\right]$
$\displaystyle+\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}\right].$
(14)

Then, it holds for term (14) that

$\displaystyle
E\left(\frac{1}{n}\sum_{i=1}^{n}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}\right]\right)$
$\displaystyle=E\left\\{E\left(\frac{1}{n}\sum_{i=1}^{n}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}\right]\bigg{|}X_{1},\dots,X_{n},Y_{1},\dots,Y_{n}\right)\right\\}$
$\displaystyle=E\left(\frac{1}{n}\sum_{i=1}^{n}\left[\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}E\left(\frac{Z_{i}}{\pi}-1\bigg{|}X_{1},\dots,X_{n},Y_{1},\dots,Y_{n}\right)\right]\right)$
$\displaystyle=E\left(\frac{1}{n}\sum_{i=1}^{n}\left[\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}E\left(\frac{Z_{i}}{\pi}-1\right)\right]\right)$
$\displaystyle=0,$

where the third equality is a consequence of $Z$ being jointly independence of
($X$, $Y$) and the last equality follows from
$E\left(\frac{Z}{\pi}-1\right)=0$.

Its variance equals

$\displaystyle
Var\left(\frac{1}{n}\sum_{i=1}^{n}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}\right]\right)$
$\displaystyle=E\left\\{Var\left(\frac{1}{n}\sum_{i=1}^{n}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}\right]\bigg{|}X_{1},\dots,X_{n},Y_{1},\dots,Y_{n}\right)\right\\}$
$\displaystyle+Var\left\\{E\left(\frac{1}{n}\sum_{i=1}^{n}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}\right]\bigg{|}X_{1},\dots,X_{n},Y_{1},\dots,Y_{n}\right)\right\\}$
$\displaystyle=E\left\\{Var\left(\frac{1}{n}\sum_{i=1}^{n}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}\right]\bigg{|}X_{1},\dots,X_{n},Y_{1},\dots,Y_{n}\right)\right\\}$
$\displaystyle=\frac{1}{n^{2}}E\left(\sum_{i=1}^{n}Var\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}\bigg{|}X_{1},\dots,X_{n},Y_{1},\dots,Y_{n}\right]\right)$
$\displaystyle=\frac{1}{n^{2}}E\left(\sum_{i=1}^{n}E\left[\left(\frac{Z_{i}}{\pi}-1\right)^{2}\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}^{2}\bigg{|}X_{1},\dots,X_{n},Y_{1},\dots,Y_{n}\right]\right)$
$\displaystyle=\frac{1}{n^{2}}E\left[\sum_{i=1}^{n}\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}^{2}E\left\\{\left(\frac{Z_{i}}{\pi}-1\right)^{2}\bigg{|}X_{1},\dots,X_{n},Y_{1},\dots,Y_{n}\right\\}\right]$
$\displaystyle=\frac{1}{n}E\left\\{\left(\frac{Z}{\pi}-1\right)^{2}\right\\}E\left[\frac{1}{n}\sum_{i=1}^{n}\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}^{2}\right],$

where $E\left\\{\left(\frac{Z}{\pi}-1\right)^{2}\right\\}$ is bounded and
$E\left[\frac{1}{n}\sum_{i=1}^{n}\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}^{2}\right]$
converges to zero, assuming that
$\frac{1}{n}\sum_{i=1}^{n}\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}^{2}$ is
uniformly integrable (see Assumption 5), as
$\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left\\{h^{*}(X_{i})-\hat{h}(X_{i})\right\\}^{2}}=o_{P_{n}}(1)$.
Therefore, application of Chebyshev’s inequality shows that term (14)
converges in probability to zero. Consequently,

$\displaystyle\sqrt{n}(\hat{\mu}_{1,SN}^{(\pi)}-\mu_{1})$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h^{*}(X_{i})\right\\}+h^{*}(X_{i})-\mu_{1}\right]+o_{P_{n}}(1).$

Following a similar reasoning as in Appendix A.2, we have

$\displaystyle\sqrt{n}(\hat{\theta}_{SN}-\theta)$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left\\{\left(\frac{Z_{i}}{\pi}\left\\{Y_{i}-h^{*}(X_{i})\right\\}+h^{*}(X_{i})-E\left[\frac{Z}{\pi^{2}}\left\\{Y-h^{*}(X)\right\\}\right](Z_{i}-\pi)-\mu_{1}\right)\right.$
$\displaystyle-\left.\left(\frac{1-Z_{i}}{1-\pi}\left\\{Y_{i}-h^{*}(X_{i})\right\\}+h^{*}(X_{i})+E\left[\frac{1-Z}{(1-\pi)^{2}}\left\\{Y-h^{*}(X)\right\\}\right](Z_{i}-\pi)-\mu_{0}\right)\right\\}+o_{P_{n}}(1),$

which can be simplified to

$\displaystyle\sqrt{n}(\hat{\theta}_{SN}-\theta)$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(\frac{Z_{i}-\pi}{\pi(1-\pi)}\left\\{Y_{i}-h^{*}(X_{i})\right\\}-E\left[\left(\frac{Z}{\pi^{2}}+\frac{1-Z}{(1-\pi)^{2}}\right)\left\\{Y-h^{*}(X)\right\\}\right](Z_{i}-\pi)-\left(\mu_{1}-\mu_{0}\right)\right)$
$\displaystyle+o_{P_{n}}(1).$

Inference can be conducted by calculating the standard error of
$\hat{\theta}_{SN}$ as $1/n$ times the sample variance of

$\frac{Z_{i}-\hat{\pi}}{\hat{\pi}(1-\hat{\pi})}\left\\{Y_{i}-\hat{h}(X_{i})\right\\}-\frac{1}{n}\sum_{j=1}^{n}\left[\left(\frac{Z_{j}}{\hat{\pi}^{2}}+\frac{1-Z_{j}}{(1-\hat{\pi})^{2}}\right)\left\\{Y_{j}-\hat{h}(X_{j})\right\\}\right](Z_{i}-\hat{\pi}).$

####  A.2.2 Finite sample unbiasedness

We prove that $\hat{\theta}_{SS}$ is exactly unbiased, even in finite samples,
when the the known randomization probability $\pi$ (e.g., 0.5 for a 1:1
randomized trial) is plugged-in for the propensity score. Thus we show that

$\hat{\theta}_{k}=\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\frac{Z_{i}}{\pi}\left\\{Y_{i}-\hat{h}_{1,k}(X_{i})\right\\}+\hat{h}_{1,k}(X_{i})-\left[\frac{1-Z_{i}}{1-\pi}\left\\{Y_{i}-\hat{h}_{0,k}(X_{i})\right\\}+\hat{h}_{0,k}(X_{i})\right],$

has mean $E\left\\{Y(1)-Y(0)\right\\}$ for each $k=1,\dots,K$ with $Y(j)$
denoting the counterfactual outcome under $Z=j$ for $j=0,1$. Following a
similar proof as for term (10), we have

$\displaystyle
E\left(\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\frac{Z_{i}}{\pi}\left\\{Y_{i}-\hat{h}_{1,k}(X_{i})\right\\}+\hat{h}_{1,k}(X_{i})-\left[\frac{1-Z_{i}}{1-\pi}\left\\{Y_{i}-\hat{h}_{0,k}(X_{i})\right\\}+\hat{h}_{0,k}(X_{i})\right]\right)$
$\displaystyle=E\left[\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\frac{Z_{i}}{\pi}Y_{i}-\frac{1-Z_{i}}{1-\pi}Y_{i}+\left\\{1-\frac{Z_{i}}{\pi}\right\\}\hat{h}_{1,k}(X_{i})-\left\\{1-\frac{1-Z_{i}}{1-\pi}\right\\}\hat{h}_{0,k}(X_{i})\right]$
$\displaystyle=E\left[\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\frac{Z_{i}}{\pi}Y_{i}-\frac{1-Z_{i}}{1-\pi}Y_{i}\right]+E\left[\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{1-\frac{Z_{i}}{\pi}\right\\}\hat{h}_{1,k}(X_{i})-\left\\{1-\frac{1-Z_{i}}{1-\pi}\right\\}\hat{h}_{0,k}(X_{i})\right]$
$\displaystyle=E\left[Y(1)-Y(0)\right]+E\left[\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{1-\frac{Z_{i}}{\pi}\right\\}\hat{h}_{1,k}(X_{i})-\left\\{1-\frac{1-Z_{i}}{1-\pi}\right\\}\hat{h}_{0,k}(X_{i})\right],$

because of consistency and randomization. Relying on sample splitting, the
second expectation can be rewritten as

$\displaystyle
E\left[\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{1-\frac{Z_{i}}{\pi}\right\\}\hat{h}_{1,k}(X_{i})-\left\\{1-\frac{1-Z_{i}}{1-\pi}\right\\}\hat{h}_{0,k}(X_{i})\right]$
$\displaystyle=E\left(E\left[\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{1-\frac{Z_{i}}{\pi}\right\\}\hat{h}_{1,k}(X_{i})-\left\\{1-\frac{1-Z_{i}}{1-\pi}\right\\}\hat{h}_{0,k}(X_{i})\bigg{|}\mathbf{O}^{\bar{k}}\right]\right)$
$\displaystyle=E\left(E\left[\left\\{1-\frac{Z}{\pi}\right\\}\hat{h}_{1,k}(X)-\left\\{1-\frac{1-Z}{1-\pi}\right\\}\hat{h}_{0,k}(X)\bigg{|}\mathbf{O}^{\bar{k}}\right]\right)$
$\displaystyle=E\left(1-\frac{Z}{\pi}\right)E\left[E\left\\{\hat{h}_{1,k}(X)\bigg{|}\mathbf{O}^{\bar{k}}\right\\}\right]-E\left(1-\frac{1-Z}{1-\pi}\right)E\left[E\left\\{\hat{h}_{0,k}(X)\bigg{|}\mathbf{O}^{\bar{k}}\right\\}\right]$
$\displaystyle=0,$

as $E\left(\frac{Z}{\pi}-1\right)=0$ and $E\left(1-\frac{1-Z}{1-\pi}\right)=0$
because of randomization.

Under the strong null, sample splitting can even be avoided. Specifically,
under the strong null, $\hat{\theta}_{SN}=\hat{\mu}_{1,SN}-\hat{\mu}_{0,SN}$
also has the advantage of being exactly unbiased, even in finite samples as
long as the estimated probability $\hat{\pi}$ is replaced with the known
randomization probability $\pi$ (e.g., 0.5 for a 1:1 randomized trial). A
proof for this result follows a similar reasoning as in Appendix A.2.1 by
using the independence of $Z$ and $\hat{h}(X)$.

###  A.3 Data-adaptive covariate adjustment without sample splitting: Outcome
regression fitted in a parametric model

In this section we prove that we can avoid cross-fitting when a generalized
linear model with canonical link for the outcome is fitted via maximum
likelihood, possibly in a data-adaptive way (using consistent variable
selection, regardless of loss function). We denote the regression coefficients
obtained by the data-adaptive procedure in Section 4 as
$\tilde{\boldsymbol{\gamma}}=(\tilde{\gamma}_{0},\tilde{\gamma}_{1})^{\prime}$
for estimation of $\mu_{1}$, and
$\tilde{\boldsymbol{\eta}}=(\tilde{\eta}_{0},\tilde{\eta}_{1})^{\prime}$ for
estimation of $\mu_{0}$.

The “Model fitting” step (Step 1b) in Section 4 ensures that

$\displaystyle\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\hat{\pi}}\left\\{Y_{i}-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}+h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right]$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}+h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right]$
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}Z_{i}\left\\{Y_{i}-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}\left(\frac{1}{\hat{\pi}}-\frac{1}{\pi}\right)\sqrt{n}$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}+h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right],$
(15)

as
$\frac{1}{n}\sum_{i=1}^{n}Z_{i}\left\\{Y_{i}-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}=0$
for generalized linear models with a canonical link function and intercept,
fitted via maximum likelihood.

In order to handle term (15), we define
$\hat{\mu}_{1,DA}^{(\pi)}=\frac{1}{n}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}+h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right]$.
Then, it holds that

$\displaystyle\sqrt{n}(\hat{\mu}_{1,DA}^{(\pi)}-\mu_{1})$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}(X_{i};\boldsymbol{\gamma}^{*})\right\\}+h_{1}(X_{i};\boldsymbol{\gamma}^{*})-\mu_{1}\right]$
$\displaystyle+\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h_{1}(X_{i};\boldsymbol{\gamma}^{*})-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}\right].$
(16)

Next, following a Taylor expansion, term (16) can be rewritten as

$\displaystyle\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\left(\frac{Z_{i}}{\pi}-1\right)\frac{\partial}{\partial\boldsymbol{\gamma}}h_{1}(X_{i};\boldsymbol{\gamma})\biggr{\rvert}_{\boldsymbol{\gamma}=\boldsymbol{\gamma}^{*}}\right](\tilde{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^{*})+O_{P_{n}}(\sqrt{n}||\boldsymbol{\gamma}^{*}-\tilde{\boldsymbol{\gamma}}||^{2}_{2}),$

assuming that $\max_{j\leq p,j^{\prime}\leq p,i\leq
n}\left|\frac{\partial^{2}}{\partial\gamma_{j}\partial\gamma_{j^{\prime}}}h_{1}(X_{i};\boldsymbol{\gamma})\biggr{\rvert}_{(\gamma_{j}=\gamma^{*}_{j},\gamma_{j^{\prime}}=\gamma^{*}_{j^{\prime}})}\right|$
is bounded, with $\gamma_{j}$ and $\gamma^{*}_{j}$ denoting component $j$
($j\in 1,\dots,p$) of $\boldsymbol{\gamma}$ and $\boldsymbol{\gamma}^{*}$,
respectively.

Define
$\tilde{\Sigma}_{j}=\frac{1}{n}\sum_{i=1}^{n}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left(\frac{\partial}{\partial\gamma_{j}}h_{1}(X_{i};\boldsymbol{\gamma})\biggr{\rvert}_{\gamma_{j}=\gamma^{*}_{j}}\right)\right]^{2}$.
Then, by Hölder’s inequality

$\displaystyle\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\left(\frac{Z_{i}}{\pi}-1\right)\frac{\partial}{\partial\boldsymbol{\gamma}}h_{1}(X_{i};\boldsymbol{\gamma})\biggr{\rvert}_{\boldsymbol{\gamma}=\boldsymbol{\gamma}^{*}}\right](\tilde{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^{*})\right|$
$\displaystyle\leq\left(\max_{j\leq
p}\tilde{\Sigma}_{j}^{1/2}\right)\max_{j\leq
p}\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\frac{\left(\frac{Z_{i}}{\pi}-1\right)\left(\frac{\partial}{\partial\gamma_{j}}h_{1}(X_{i};\boldsymbol{\gamma})\biggr{\rvert}_{\gamma_{j}=\gamma^{*}_{j}}\right)}{\tilde{\Sigma}_{j}^{1/2}}\right|\left|\left|\tilde{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^{*}\right|\right|_{1}$
$\displaystyle\leq
O(1)O_{P_{n}}(\sqrt{\log(p)})O_{P_{n}}\left(s_{\gamma}\sqrt{\log(p\lor
n)/n}\right)$ $\displaystyle=o_{P_{n}}(1),$

where the second inequality follows from the bound on the center factor by
applying the moderate deviation theory for self-normalized sums of de la Peña
et al. [2009], and in particular Lemma 5 in Belloni et al. [2012] and the
verification of Assumption 3 in Farrell [2015].

Consequently,

$\displaystyle\sqrt{n}(\hat{\mu}_{1,DA}^{(\pi)}-\mu_{1})$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}(X_{i};\boldsymbol{\gamma}^{*})\right\\}+h_{1}(X_{i};\boldsymbol{\gamma}^{*})-\mu_{1}\right]+o_{P_{n}}(1),$

and thus

$\displaystyle\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}+h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right]$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}(X_{i};\boldsymbol{\gamma}^{*})\right\\}+h_{1}(X_{i};\boldsymbol{\gamma}^{*})\right]+o_{P_{n}}(1).$

Finally, we have

$\displaystyle\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\hat{\pi}}\left\\{Y_{i}-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}+h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right]$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}(X_{i};\boldsymbol{\gamma}^{*})\right\\}+h_{1}(X_{i};\boldsymbol{\gamma}^{*})\right]+o_{P_{n}}(1),$
(17)

which leads to

$\displaystyle\sqrt{n}(\hat{\mu}_{1,DA}-\mu_{1})$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}(X_{i};\boldsymbol{\gamma}^{*})\right\\}+h_{1}(X_{i};\boldsymbol{\gamma}^{*})-\mu_{1}\right]+o_{P_{n}}(1).$
(18)

Following a similar reasoning for $\hat{\mu}_{0,DA}$, we have

$\displaystyle\sqrt{n}(\hat{\theta}_{DA}-\theta)$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}(X_{i};\boldsymbol{\gamma}^{*})\right\\}+h_{1}(X_{i};\boldsymbol{\gamma}^{*})-\mu_{1}\right]-\left[\frac{1-Z_{i}}{1-\pi}\left\\{Y_{i}-h_{0}(X_{i};\boldsymbol{\eta}^{*})\right\\}+h_{0}(X_{i};\boldsymbol{\eta}^{*})-\mu_{0}\right]\right)$
$\displaystyle+o_{P_{n}}(1).$

The asymptotic variance of
$\hat{\theta}_{DA}=\hat{\mu}_{1,DA}-\hat{\mu}_{0,DA}$ can be calculated as
$1/n$ times the sample variance of the values

$\displaystyle\frac{Z_{i}}{\hat{\pi}}(Y_{i}-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}}))+h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})-\left[\frac{1-Z_{i}}{1-\hat{\pi}}(Y_{i}-h_{0}(X_{i};\tilde{\boldsymbol{\eta}}))+h_{0}(X_{i};\tilde{\boldsymbol{\eta}})\right].$

As opposed to the asymptotic variance for
$\hat{\theta}_{SS}=\hat{\mu}_{1,SS}-\hat{\mu}_{0,SS}$, there is no need to
take into account the uncertainty in the randomization probabilities as there
are no terms like
$\frac{1}{n}\sum_{j=1}^{n}\left[\frac{Z_{j}}{\hat{\pi}^{2}}\left\\{Y_{j}-\hat{h}_{1}(X_{j};\tilde{\boldsymbol{\gamma}})\right\\}\right](Z_{i}-\hat{\pi})$
and
$\frac{1}{n}\sum_{j=1}^{n}\left[\frac{1-Z_{j}}{(1-\hat{\pi})^{2}}\left\\{Y_{j}-\hat{h}_{0}(X_{j};\tilde{\boldsymbol{\eta}})\right\\}\right](Z_{i}-\hat{\pi})$
for participant $i$.

###  A.4 Targeted maximum likelihood estimation

As described in Section 4.2, the data-adaptive approach discussed in Section 4
is closely related to targeted maximum likelihood estimation (TMLE) [see e.g.,
Rosenblum and Van Der Laan, 2010]. In this appendix, we show that the validity
of the inference based on this TMLE approach is asymptotically insensitive to
selection mistakes. In particular, standard errors obtained using the
aforementioned calculation are unbiased in large samples, even when the chosen
model is misspecified. This has been conjectured before by Moore and van der
Laan [2009], and has more recently been shown by Williams et al. [2021] for
time-to-event and ordinal outcomes. Their result can be extended directly to
continuous and binary outcomes. Such a general proof is given in Appendix
A.4.1. This is also in line with the findings in Benkeser et al. [2017]; their
approach to obtain doubly robust non-parametric inference for the average
treatment effect reduces to standard TMLE – as described here – in randomized
clinical trials.

####  A.4.1 Cross-validated targeted minimum-loss-based estimation

In this appendix, we discuss cross-validated targeted minimum-loss-based
estimation (CV-TMLE), another imputation estimator that ensures valid
inference when data-adaptive prediction models are used, irrespective of model
misspecification. Its simplicity makes it an appealing estimator for valid
inference and estimation in randomized trials. This approach first constructs
an initial prediction $\hat{h}_{1,k}(X_{i})$ for participant $i$ belonging to
a certain fold $k=1,\dots,K$ based on the $K-1$ remaining folds, which is then
updated in a clever way to $\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})$ so that
[Zheng and van der Laan, 2011]

$\frac{1}{n}\sum_{i=1}^{n}\left[Z_{i}\left\\{Y_{i}-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right\\}\right]=0.$

This ensures that

$\frac{1}{K}\sum_{k=1}^{K}\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[Z_{i}\left\\{Y_{i}-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right\\}\right]=0,$

if $n/K=|\mathcal{V}_{k}|$ for all $k\leq K$, and

$\frac{1}{K}\sum_{k=1}^{K}\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[Z_{i}\left\\{Y_{i}-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right\\}\right]=o_{P_{n}}(1),$

otherwise (assuming $K$ does not grow with $n$). Specifically, the outcome
model fit $\hat{h}_{1,k}(X)$ is updated by fitting the following regression
model with single parameter $\epsilon$ and offset based on the initial
predictions $\hat{h}_{1,k}(X)$:

$E(Y|Z=1,X)=g^{-1}\left\\{g(\hat{h}_{1,k}(X))+\epsilon\right\\},$

with $g$ the considered canonical link function (e.g., identity for continuous
outcome and logit for binary outcome). This can be done using standard
statistical software for fitting a generalized linear regression (with
canonical link) of $Y$ with an intercept (denoted as $\epsilon$) using offset
$g(\hat{h}_{1,k}(X))$ among participants with $Z=1$ pooled over all folds. In
fact, the updating step can also be performed on the validation sample, which
leads to $K$ different estimates $\hat{\epsilon}_{1},\dots,\hat{\epsilon}_{k}$
for $\epsilon$ (i.e., one for each fold).

This results in the imputation estimator

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon}).$
(19)

Alternatively, one can consider
$\hat{\mu}_{1,TMLE}=\frac{1}{K}\sum_{k=1}^{K}\hat{\mu}_{1,k,TMLE}$ with
$\hat{\mu}_{1,k,TMLE}=\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})$
for fold $k$, which equals (19) if $n/K=|\mathcal{V}_{k}|$ for all $k\leq K$
and is equal up to a term that’s asymptotically negligible otherwise. In what
follows, we therefore focus on
$\hat{\mu}_{1,TMLE}=\frac{1}{K}\sum_{k=1}^{K}\hat{\mu}_{1,k,TMLE}$ to be
consistent with previous sections.

The update of the predictions is the reason why the proof in Appendix A.2 no
longer holds for CV-TMLE as the predictions
$\hat{h}_{1,k}^{(u)}(X;\hat{\epsilon})$ are no longer obtained from a
‘separate sample’. In particular, an update of the proof for term (10) is
required.

To this end, we assume that $\hat{\epsilon}$ converges in probability to a
limit $\epsilon^{*}$ in the sense that
$|\hat{\epsilon}-\epsilon^{*}|=o_{P_{n}}(1)$ since $\hat{\epsilon}$ is
obtained with maximum likelihood estimation, and that
$\hat{h}_{1,k}^{(u)}(X_{i};\epsilon^{*})$ converges in probability to
$h_{1}^{*}(X_{i};\epsilon^{*})$ in the sense that
$\sqrt{\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{h_{1}^{*}(X_{i};\epsilon^{*})-\hat{h}_{1,k}^{(u)}(X_{i};\epsilon^{*})\right\\}^{2}}=o_{P_{n}}(1)$,
$\forall k\in\\{1,\dots,K\\}$ (see also Assumption 4′ (Convergence in
probability of predictions for sample splitting)).

First, the ‘clever’ update ensures that

$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{\hat{\pi}}\left\\{Y_{i}-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right\\}+\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right]$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right\\}+\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right]$
$\displaystyle+\frac{1}{K}\sum_{k=1}^{K}\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{Y_{i}-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right\\}\left(\frac{1}{\hat{\pi}}-\frac{1}{\pi}\right)$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right\\}+\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right],$

if $n/K=|\mathcal{V}_{k}|$ for all $k\leq K$; otherwise, the left hand side
equals

$\frac{1}{K}\sum_{k=1}^{K}\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right\\}+\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right]+o_{P_{n}}(1).$

In line with the previous proofs, we focus on

$\displaystyle\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right\\}+\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right].$
(20)

In order to handle term (20), we define
$\hat{\mu}_{1,k,TMLE}^{(\pi)}=\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right\\}+\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right]$.
Then, it holds that

$\displaystyle\sqrt{|\mathcal{V}_{k}|}(\hat{\mu}_{1,k,TMLE}^{(\pi)}-\mu_{1})$
$\displaystyle=\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}^{*}(X_{i};\epsilon^{*})\right\\}+h^{*}_{1}(X_{i};\epsilon^{*})-\mu_{1}\right]$
$\displaystyle+\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h^{*}_{1}(X_{i};\epsilon^{*})-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right\\}\right].$
(21)

Next, for term (21) it holds that

$\displaystyle\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h^{*}_{1}(X_{i};\epsilon^{*})-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right\\}\right]$
$\displaystyle=\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{h^{*}_{1}(X_{i};\epsilon^{*})-\hat{h}_{1,k}^{(u)}(X_{i};\epsilon^{*})\right\\}\right]$
(22)
$\displaystyle+\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{\pi}-1\right)\left\\{\hat{h}_{1,k}^{(u)}(X_{i};\epsilon^{*})-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon})\right\\}\right].$
(23)

Under the assumption that
$\sqrt{\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{h_{1}^{*}(X_{i};\epsilon^{*})-\hat{h}_{1,k}^{(u)}(X_{i};\epsilon^{*})\right\\}^{2}}=o_{P_{n}}(1)$,
one can prove that term (22) converges in probability to zero by following a
similar approach as for term (10). Moreover, following a Taylor expansion for
term (23) gives

$\displaystyle\sqrt{|\mathcal{V}_{k}|}(\epsilon^{*}-\hat{\epsilon})\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left(\frac{Z_{i}}{\pi}-1\right)\frac{\partial}{\partial\epsilon}\hat{h}_{1,k}^{(u)}(X_{i};\epsilon)\biggr{\rvert}_{\epsilon=\epsilon^{*}}+O_{P_{n}}\left(\sqrt{|\mathcal{V}_{k}|}|\epsilon^{*}-\hat{\epsilon}|^{2}\right).$
(24)

For the main term in Equation (24), we have

$\displaystyle
E\left[E\left\\{\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left(\frac{Z_{i}}{\pi}-1\right)\frac{\partial}{\partial\epsilon}\hat{h}_{1,k}^{(u)}(X_{i};\epsilon)\biggr{\rvert}_{\epsilon=\epsilon^{*}}\bigg{|}\mathbf{O}^{\bar{k}}\right\\}\right]$
$\displaystyle=E\left[E\left\\{\left(\frac{Z}{\pi}-1\right)\frac{\partial}{\partial\epsilon}\hat{h}_{1,k}^{(u)}(X;\epsilon)\biggr{\rvert}_{\epsilon=\epsilon^{*}}\bigg{|}\mathbf{O}^{\bar{k}}\right\\}\right]$
$\displaystyle=E\left(\frac{Z}{\pi}-1\right)E\left[E\left\\{\frac{\partial}{\partial\epsilon}\hat{h}_{1,k}^{(u)}(X;\epsilon)\biggr{\rvert}_{\epsilon=\epsilon^{*}}\bigg{|}\mathbf{O}^{\bar{k}}\right\\}\right]$
$\displaystyle=0.$

By the law of total variance

$\displaystyle
Var\left[\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left(\frac{Z_{i}}{\pi}-1\right)\frac{\partial}{\partial\epsilon}\hat{h}_{1,k}^{(u)}(X_{i};\epsilon)\biggr{\rvert}_{\epsilon=\epsilon^{*}}\right]$
$\displaystyle=E\left[Var\left\\{\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left(\frac{Z_{i}}{\pi}-1\right)\frac{\partial}{\partial\epsilon}\hat{h}_{1,k}^{(u)}(X_{i};\epsilon)\biggr{\rvert}_{\epsilon=\epsilon^{*}}\bigg{|}\mathbf{O}^{\bar{k}}\right\\}\right]$
$\displaystyle=E\left[\frac{1}{|\mathcal{V}_{k}|}Var\left\\{\left(\frac{Z}{\pi}-1\right)\frac{\partial}{\partial\epsilon}\hat{h}_{1,k}^{(u)}(X;\epsilon)\biggr{\rvert}_{\epsilon=\epsilon^{*}}\bigg{|}\mathbf{O}^{\bar{k}}\right\\}\right]$
$\displaystyle=\frac{1}{|\mathcal{V}_{k}|}E\left(E\left[\left\\{\left(\frac{Z}{\pi}-1\right)\frac{\partial}{\partial\epsilon}\hat{h}_{1,k}^{(u)}(X;\epsilon)\biggr{\rvert}_{\epsilon=\epsilon^{*}}\right\\}^{2}\bigg{|}\mathbf{O}^{\bar{k}}\right]\right)$
$\displaystyle=\frac{1}{|\mathcal{V}_{k}|}E\left\\{\left(\frac{Z}{\pi}-1\right)^{2}\right\\}E\left(E\left[\left\\{\frac{\partial}{\partial\epsilon}\hat{h}_{1,k}^{(u)}(X;\epsilon)\biggr{\rvert}_{\epsilon=\epsilon^{*}}\right\\}^{2}\bigg{|}\mathbf{O}^{\bar{k}}\right]\right)$
$\displaystyle=\frac{1}{|\mathcal{V}_{k}|}E\left\\{\left(\frac{Z}{\pi}-1\right)^{2}\right\\}E\left[\left\\{\frac{\partial}{\partial\epsilon}\hat{h}_{1,k}^{(u)}(X;\epsilon)\biggr{\rvert}_{\epsilon=\epsilon^{*}}\right\\}^{2}\right],$

Under the assumption that
$E\left[\left\\{\frac{\partial}{\partial\epsilon}\hat{h}_{1,k}^{(u)}(X;\epsilon)\biggr{\rvert}_{\epsilon=\epsilon^{*}}\right\\}^{2}\right]$
is bounded (which is satisfied for generalized linear models with common
canonical link when $X$ is bounded), application of Chebyshev’s inequality
shows that

$\displaystyle\left|\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left(\frac{Z_{i}}{\pi}-1\right)\frac{\partial}{\partial\epsilon}\hat{h}_{1,k}^{(u)}(X_{i};\epsilon)\biggr{\rvert}_{\epsilon=\epsilon^{*}}\right|=O_{P_{n}}\left(\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\right).$

Then,

$\displaystyle\left|\sqrt{|\mathcal{V}_{k}|}(\epsilon^{*}-\hat{\epsilon})\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left(\frac{Z_{i}}{\pi}-1\right)\frac{\partial}{\partial\epsilon}\hat{h}_{1,k}^{(u)}(X_{i};\epsilon)\biggr{\rvert}_{\epsilon=\epsilon^{*}}\right|$
$\displaystyle\leq\sqrt{|\mathcal{V}_{k}|}\left|\epsilon^{*}-\hat{\epsilon}\right|\left|\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left(\frac{Z_{i}}{\pi}-1\right)\frac{\partial}{\partial\epsilon}\hat{h}_{1,k}^{(u)}(X_{i};\epsilon)\biggr{\rvert}_{\epsilon=\epsilon^{*}}\right|$
$\displaystyle=\sqrt{|\mathcal{V}_{k}|}o_{P_{n}}(1)O_{P_{n}}\left(\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\right)$
$\displaystyle=o_{P_{n}}(1),$

as $\left|\epsilon^{*}-\hat{\epsilon}\right|=o_{P_{n}}(1)$ (see e.g., Zheng
and van der Laan [2011]). Moreover, following a similar reasoning, it holds
that
$O_{P_{n}}\left(\sqrt{|\mathcal{V}_{k}|}|\epsilon^{*}-\hat{\epsilon}|^{2}\right)=o_{P_{n}}(1)$
in Equation (24), as $|\hat{\epsilon}-\epsilon^{*}|=o_{P_{n}}(1)$ and assuming
$\hat{h}_{1,k}^{(u)}(X_{i};\epsilon)$ is an analytic function which ensures
that
$\frac{1}{k!}\left|\frac{\partial^{k}}{\partial\epsilon^{k}}\hat{h}_{1,k}^{(u)}(X_{i};\epsilon)\biggr{\rvert}_{\epsilon=\epsilon^{*}}\right|$
is bounded for every compact set [Krantz and Parks, 2002]. Note that this
assumption is satisfied for generalized linear models with common canonical
link such as linear regression, logistic regression, and Poisson regression
over a compact parameter space.

Then,

$\displaystyle\sqrt{n}(\hat{\mu}_{1,TMLE}-\mu_{1})$
$\displaystyle=\sqrt{n}\frac{1}{K}\sum_{k=1}^{K}\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{\left(\frac{Z_{i}}{\pi}\left\\{Y_{i}-h_{1}^{*}(X_{i};\epsilon^{*})\right\\}+h_{1}^{*}(X_{i};\epsilon^{*})-\mu_{1}\right)\right\\}+o_{P_{n}}(1).$

We can follow a similar reasoning for $\hat{\mu}_{0,TMLE}$.

####  A.4.2 Targeted maximum likelihood estimation without sample splitting:
Outcome regression fitted in a parametric model

A proof for the validity of TMLE without sample splitting can be obtained by
following a similar reasoning as in Appendix A.3 and A.4.1, under similar
assumptions as in A.3 (e.g., Assumption 1 and Assumption 2).

###  A.5 Parametric estimation of propensity score

Till now, we have estimated the propensity score $\pi=P(Z=1)$ independently of
baseline variables $X$. In this section, we prove that the results above hold
in more generality when the treatment probability is estimated within a
parametric model (e.g., using a logistic regression model). This might be of
interest in practice as estimating the propensity score in a parametric model
may yield larger efficiency gains. Specifically, we define
$P(Z=1|X)=p(X;\beta_{0})$ where $p(X;\beta)$ is a known function evaluated at
a parameter $\beta$ with population value $\beta_{0}$. For example, as $Z$ is
binary, its association with $X$ can be modeled via a logistic regression.

####  A.5.1 Data-adaptive covariate adjustment with sample splitting

Since $\pi=p(X;\beta_{0})$ we can then follow a similar reasoning as before,
except for Term (7) and Term (8) as these involve $p(X;\hat{\beta})$, with
$\hat{\beta}$ an estimate of $\beta_{0}$ obtained via estimation within a
parametric model (e.g., using logistic regression).

First, for Term (7), we have

$\displaystyle\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}$
$\displaystyle
Z_{i}\left\\{Y_{i}-h_{1}^{*}(X_{i})\right\\}\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\sqrt{|\mathcal{V}_{k}|}$
$\displaystyle=-\sqrt{|\mathcal{V}_{k}|}(\hat{\beta}-\beta_{0})\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{Y_{i}-h_{1}^{*}(X_{i})\right\\}\frac{1}{p(X_{i};\beta_{0})^{2}}\frac{\partial}{\partial\beta}p(X_{i};\beta)\biggr{\rvert}_{\beta=\beta_{0}}+o_{P_{n}}(1),$

where the $o_{P_{n}}(1)$ term follows from the fact that
$|\hat{\beta}-\beta_{0}|^{2}=O_{P_{n}}(n^{-1})$ when estimating the propensity
score via maximum likelihood estimation. Under suitable regularity conditions
[see e.g., Newey and McFadden, 1994], it follows from the uniform weak law of
large numbers that this equals

$\displaystyle-\sqrt{|\mathcal{V}_{k}|}(\hat{\beta}-\beta_{0})E\left[Z\left\\{Y-h_{1}^{*}(X)\right\\}\frac{1}{p(X;\beta_{0})^{2}}\frac{\partial}{\partial\beta}p(X;\beta)\biggr{\rvert}_{\beta=\beta_{0}}\right]+o_{P_{n}}(1).$

Here, the first term is $O_{P_{n}}(1)$ as
$|\hat{\beta}-\beta_{0}|=O_{P_{n}}(n^{-1/2})$.

Next, we define $s(X;\beta)$ as the score function of $\beta$. It then follows
from Tsiatis [2006] that

$\displaystyle\sqrt{|\mathcal{V}_{k}|}(\hat{\beta}-\beta_{0})=-\left[E\left\\{\frac{\partial
s(X;\beta_{0})}{\partial\beta}\right\\}\right]^{-1}\left\\{\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}s(X_{i};\beta_{0})\right\\}+o_{P_{n}}(1).$

Thus, Term (7) can be rewritten as

$\displaystyle
E\left[Z\left\\{Y-h_{1}^{*}(X)\right\\}\frac{1}{p(X;\beta_{0})^{2}}\frac{\partial}{\partial\beta}p(X;\beta)\biggr{\rvert}_{\beta=\beta_{0}}\right]\left[E\left\\{\frac{\partial
s(X;\beta_{0})}{\partial\beta}\right\\}\right]^{-1}\left\\{\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}s(X_{i};\beta_{0})\right\\}+o_{P_{n}}(1).$

Next, from the the Cauchy Schwarz inequality, it follows that Term (8) yields

$\displaystyle\sqrt{|\mathcal{V}_{k}|}$
$\displaystyle\left|\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}\left\\{\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right\\}\right|$
$\displaystyle\leq\sqrt{|\mathcal{V}_{k}|}\left(\max_{i\in\mathcal{V}_{k}}\frac{1}{p(X_{i};\beta_{0})p(X_{i};\hat{\beta})}\right)\sqrt{\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{h_{1}^{*}(X_{i})-\hat{h}_{1,k}(X_{i})\right\\}^{2}}\sqrt{\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{p(X_{i};\hat{\beta})-p(X_{i};\beta_{0})\right\\}^{2}}$
$\displaystyle=o_{P_{n}}(1),$

as $1>1-\rho\geq p(X_{i};\beta_{0})\geq\rho>0$ for $i=1,\dots n$, Assumption
4′ (Convergence in probability of predictions for sample splitting) and the
fact that $\hat{\beta}$ is a root-$n$ consistent maximum likelihood estimator
of a parameter in a parametric model.

Consequently,

$\displaystyle\sqrt{|\mathcal{V}_{k}|}(\hat{\mu}_{1,k}-\mu_{1})$
$\displaystyle=\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left(\vphantom{E\left[A\left\\{Y-m^{*}(L)\right\\}\frac{1}{p(L;\beta_{0})^{2}}\frac{\partial}{\partial\beta}p(L;\beta)\biggr{\rvert}_{\beta=\beta_{0}}\right]}\frac{Z_{i}}{p(X_{i};\beta_{0})}\left\\{Y_{i}-h_{1}^{*}(X_{i})\right\\}+h_{1}^{*}(X_{i})+\right.$
$\displaystyle\left.E\left[Z\left\\{Y-h_{1}^{*}(X)\right\\}\frac{1}{p(X;\beta_{0})^{2}}\frac{\partial}{\partial\beta}p(X;\beta)\biggr{\rvert}_{\beta=\beta_{0}}\right]\left[E\left\\{\frac{\partial
s(X;\beta_{0})}{\partial\beta}\right\\}\right]^{-1}s(X_{i};\beta_{0})-\mu_{1}\right)+o_{P_{n}}(1).$

We can follow a similar reasoning for $\hat{\mu}_{0,k}$.

####  A.5.2 Outcome regression fitted in a parametric model with canonical
link

Recall that, when parametrically modeling the association of $Z$ with $X$,
“Step 1b: Model fitting” is performed as follows:

  * Step 1b:

Model fitting Fit a generalized linear regression model with a canonical link
(e.g., logistic, linear, …) and intercept via maximum likelihood that
regresses the outcome $Y$ on all covariates selected in Step 1a among the
treated participants with $Z=1$, using weights $1/p(X;\hat{\beta})$.

This step then ensures that

$\displaystyle\frac{1}{\sqrt{n}}\sum_{i=1}^{n}$
$\displaystyle\left[\frac{Z_{i}}{p(X_{i};\hat{\beta})}\left\\{Y_{i}-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}+h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right]$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{p(X_{i};\beta_{0})}\left\\{Y_{i}-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}+h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right]$
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}Z_{i}\left\\{Y_{i}-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\sqrt{n}$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{p(X_{i};\beta_{0})}\left\\{Y_{i}-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}+h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right]$
(25)
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}Z_{i}\left\\{Y_{i}-h_{1}(X_{i};\boldsymbol{\gamma}^{*})\right\\}\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\sqrt{n}$
(26)
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}Z_{i}\left\\{h_{1}(X_{i};\boldsymbol{\gamma}^{*})-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\sqrt{n}.$
(27)

Term (25) is equal to Term (15) as $\pi=p(X;\beta_{0})$, and can therefore be
handled in the same way. It then follows that Term (25) equals

$\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\frac{Z_{i}}{p(X_{i};\beta_{0})}\left\\{Y_{i}-h_{1}(X_{i};\boldsymbol{\gamma}^{*})\right\\}+h_{1}(X_{i};\boldsymbol{\gamma}^{*})\right]+o_{P_{n}}(1).$

Following a similar reasoning as in Appendix A.5.1, Term (26) can be rewritten
as

$\displaystyle
E\left[Z\left\\{Y-h_{1}(X_{i};\boldsymbol{\gamma}^{*})\right\\}\frac{1}{p(X;\beta_{0})^{2}}\frac{\partial}{\partial\beta}p(X;\beta)\biggr{\rvert}_{\beta=\beta_{0}}\right]\left[E\left\\{\frac{\partial
s(X;\beta_{0})}{\partial\beta}\right\\}\right]^{-1}\left\\{\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}s(X_{i};\beta_{0})\right\\}+o_{P_{n}}(1),$

which converges to zero as
$E\left[Z\left\\{Y-h_{1}(X_{i};\boldsymbol{\gamma}^{*})\right\\}\frac{1}{p(X;\beta_{0})^{2}}\frac{\partial}{\partial\beta}p(X;\beta)\biggr{\rvert}_{\beta=\beta_{0}}\right]=0$.

Next, from the the Cauchy Schwarz inequality, it follows that Term (27) yields

$\displaystyle\sqrt{n}$
$\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}Z_{i}\left\\{h_{1}(X_{i};\boldsymbol{\gamma}^{*})-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\right|$
$\displaystyle\leq\sqrt{n}\left(\max_{i\leq
n}\frac{1}{p(X_{i};\beta_{0})p(X_{i};\hat{\beta})}\right)\sqrt{\frac{1}{n}\sum_{i=1}^{n}Z_{i}\left\\{h_{1}(X_{i};\boldsymbol{\gamma}^{*})-h_{1}(X_{i};\tilde{\boldsymbol{\gamma}})\right\\}^{2}}\sqrt{\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{p(X_{i};\hat{\beta})-p(X_{i};\beta_{0})\right\\}^{2}}$
$\displaystyle=o_{P_{n}}(1),$

as $1>p(X_{i};\beta_{0})\geq\rho>0$ for $i=1,\dots n$, Assumption 4 and the
fact that $\hat{\beta}$ is a root-$n$ consistent maximum likelihood estimator
of a parameter in a parametric model.

####  A.5.3 Targeted maximum likelihood estimation

To allow for the propensity score, and thus the ‘clever covariate’, to depend
on baseline covariates $X$, we extend the notation of our outcome models. In
particular, we add an argument for the parameter $\beta$ which indexes the
propensity score model. For each $k\in{1,\dots,K}$, we then consider

$\displaystyle\frac{1}{\sqrt{|\mathcal{V}_{k}|}}$
$\displaystyle\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{p(X_{i};\hat{\beta})}\left\\{Y_{i}-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\hat{\beta})\right\\}+\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\hat{\beta})\right]$
$\displaystyle=\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{p(X_{i};\beta_{0})}\left\\{Y_{i}-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\hat{\beta})\right\\}+\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\hat{\beta})\right]$
(28)
$\displaystyle+\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{Y_{i}-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\hat{\beta})\right\\}\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\sqrt{|\mathcal{V}_{k}|}.$
(29)

For Term (28), we follow a similar reasoning as for Term (20). In particular,
defining
$\hat{\mu}_{1,k,TMLE}^{(\pi)}=\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{p(X_{i};\beta_{0})}\left\\{Y_{i}-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\hat{\beta})\right\\}+\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\hat{\beta})\right]$,
it holds that

$\displaystyle\sqrt{|\mathcal{V}_{k}|}(\hat{\mu}_{1,k,TMLE}^{(\pi)}-\mu_{1})$
$\displaystyle=\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{p(X_{i};\beta_{0})}\left\\{Y_{i}-h_{1}^{*}(X_{i};\epsilon^{*},\beta_{0})\right\\}+h^{*}_{1}(X_{i};\epsilon^{*},\beta_{0})-\mu_{1}\right]$
$\displaystyle+\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{p(X_{i};\beta_{0})}-1\right)\left\\{h^{*}_{1}(X_{i};\epsilon^{*},\beta_{0})-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\hat{\beta})\right\\}\right],$

where we have for the last term that

$\displaystyle\frac{1}{\sqrt{|\mathcal{V}_{k}|}}$
$\displaystyle\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{p(X_{i};\beta_{0})}-1\right)\left\\{h^{*}_{1}(X_{i};\epsilon^{*},\beta_{0})-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\hat{\beta})\right\\}\right]$
$\displaystyle=\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{p(X_{i};\beta_{0})}-1\right)\left\\{h^{*}_{1}(X_{i};\epsilon^{*},\beta_{0})-\hat{h}_{1,k}^{(u)}(X_{i};\epsilon^{*},\beta_{0})\right\\}\right]$
$\displaystyle+\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{p(X_{i};\beta_{0})}-1\right)\left\\{\hat{h}_{1,k}^{(u)}(X_{i};\epsilon^{*},\beta_{0})-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta_{0})\right\\}\right]$
$\displaystyle+\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{p(X_{i};\beta_{0})}-1\right)\left\\{\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta_{0})-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\hat{\beta})\right\\}\right].$

We can show that the first two terms converge to zero in a similar way as for
Term (22) and Term (23), respectively. For the last term, we have

$\displaystyle\frac{1}{\sqrt{|\mathcal{V}_{k}|}}$
$\displaystyle\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{p(X_{i};\beta_{0})}-1\right)\left\\{\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta_{0})-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\hat{\beta})\right\\}\right]$
$\displaystyle=\sqrt{|\mathcal{V}_{k}|}(\beta_{0}-\hat{\beta})\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{p(X_{i};\beta_{0})}-1\right)\frac{\partial}{\partial\beta}\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta)\biggr{\rvert}_{\boldsymbol{\beta}=\beta_{0}}\right]+\sqrt{|\mathcal{V}_{k}|}||\beta_{0}-\hat{\beta}||^{2}_{2}O_{P_{n}}(1),$

where $||\cdot||^{2}_{2}$ denotes the Euclidian norm. The remainder term
$\sqrt{|\mathcal{V}_{k}|}||\beta_{0}-\hat{\beta}||^{2}_{2}O_{P_{n}}(1)$ tends
to be of lower order than the other term. The first term converges to zero as
$\sqrt{|\mathcal{V}_{k}|}(\beta_{0}-\hat{\beta})=O_{p}(1)$ and
$\left|\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{Z_{i}}{p(X_{i};\beta_{0})}-1\right)\frac{\partial}{\partial\beta}h_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta)\biggr{\rvert}_{\beta=\beta_{0}}\right]\right|=o_{p}(1)$
by following a similar reasoning as for Term (23) (applying Chebyshev’s
inequality for a finite-dimensional vector) assuming
$E\left(\left|\left|\frac{\partial}{\partial\beta}h_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta)\biggr{\rvert}_{\beta=\beta_{0}}\right|\right|_{2}\right)$
is bounded.

Consequently,

$\displaystyle\sqrt{|\mathcal{V}_{k}|}(\hat{\mu}_{1,k,TMLE}^{(\pi)}-\mu_{1})$
$\displaystyle=\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}\left[\frac{Z_{i}}{p(X_{i};\beta_{0})}\left\\{Y_{i}-h_{1}^{*}(X_{i};\epsilon^{*},\beta_{0})\right\\}+h^{*}_{1}(X_{i};\epsilon^{*},\beta_{0})-\mu_{1}\right]+o_{P_{n}}(1).$

Term (29) can be rewritten as follows

$\displaystyle\frac{1}{|\mathcal{V}_{k}|}$
$\displaystyle\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{Y_{i}-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\hat{\beta})\right\\}\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\sqrt{|\mathcal{V}_{k}|}$
$\displaystyle=\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{Y_{i}-h^{*}_{1}(X_{i};\epsilon^{*},\beta_{0})\right\\}\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\sqrt{|\mathcal{V}_{k}|}$
(30)
$\displaystyle+\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{h^{*}_{1}(X_{i};\epsilon^{*},\beta_{0})-\hat{h}_{1,k}^{(u)}(X_{i};\epsilon^{*},\beta_{0})\right\\}\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\sqrt{|\mathcal{V}_{k}|}$
(31)
$\displaystyle+\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{\hat{h}_{1,k}^{(u)}(X_{i};\epsilon^{*},\beta_{0})-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta_{0})\right\\}\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\sqrt{|\mathcal{V}_{k}|}$
(32)
$\displaystyle+\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta_{0})-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\hat{\beta})\right\\}\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\sqrt{|\mathcal{V}_{k}|}.$
(33)

Following a similar reasoning as in Appendix A.5.1, Term (30) can be rewritten
as

$\displaystyle
E\left[Z\left\\{Y-h^{*}_{1}(X_{i};\epsilon^{*},\beta_{0})\right\\}\frac{1}{p(X;\beta_{0})^{2}}\frac{\partial}{\partial\beta}p(X;\beta)\biggr{\rvert}_{\beta=\beta_{0}}\right]\left[E\left\\{\frac{\partial
s(X;\beta_{0})}{\partial\beta}\right\\}\right]^{-1}\left\\{\frac{1}{\sqrt{|\mathcal{V}_{k}|}}\sum_{i\in\mathcal{V}_{k}}s(X_{i};\beta_{0})\right\\}+o_{P_{n}}(1),$

which converges to zero as
$E\left[Z\left\\{Y-h^{*}_{1}(X_{i};\epsilon^{*},\beta_{0})\right\\}\frac{1}{p(X;\beta_{0})^{2}}\frac{\partial}{\partial\beta}p(X;\beta)\biggr{\rvert}_{\beta=\beta_{0}}\right]=0$.
Term (31) also converges to zero by following a similar reasoning as for Term
(8) with parametrically estimated propensity scores.

Next, from the the Cauchy Schwarz inequality, it follows that Term (32) yields

$\displaystyle\sqrt{|\mathcal{V}_{k}|}\left|\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{\hat{h}_{1,k}^{(u)}(X_{i};\epsilon^{*},\beta_{0})-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta_{0})\right\\}\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\right|$
$\displaystyle\leq\sqrt{|\mathcal{V}_{k}|}\left(\max_{i\in\mathcal{V}_{k}}\frac{1}{p(X_{i};\beta_{0})p(X_{i};\hat{\beta})}\right)\sqrt{\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{\hat{h}_{1,k}^{(u)}(X_{i};\epsilon^{*},\beta_{0})-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta_{0})\right\\}^{2}}\sqrt{\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left\\{p(X_{i};\hat{\beta})-p(X_{i};\beta_{0})\right\\}^{2}}$
$\displaystyle=o_{P_{n}}(1),$

as $1>p(X_{i};\beta_{0})\geq\rho>0$ for $i=1,\dots n$,
$|\hat{\epsilon}-\epsilon^{*}|=o_{P_{n}}(1)$, and the fact that $\hat{\beta}$
are root-$n$ consistent maximum likelihood estimators of parameters in a
parametric model.

For Term (33), we have

$\displaystyle\frac{1}{|\mathcal{V}_{k}|}$
$\displaystyle\sum_{i\in\mathcal{V}_{k}}Z_{i}\left\\{\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta_{0})-\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\hat{\beta})\right\\}\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\sqrt{|\mathcal{V}_{k}|}$
$\displaystyle=\sqrt{|\mathcal{V}_{k}|}(\beta_{0}-\hat{\beta})\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\frac{\partial}{\partial\beta}\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta)\biggr{\rvert}_{\beta=\beta_{0}}\right]$
$\displaystyle+\sqrt{|\mathcal{V}_{k}|}||\beta_{0}-\hat{\beta}||^{2}_{2}O_{P_{n}}(1),$

where $||\cdot||^{2}_{2}$ denotes the Euclidian norm. The remainder term
$\sqrt{|\mathcal{V}_{k}|}||\beta_{0}-\hat{\beta}||^{2}_{2}O_{P_{n}}(1)$ tends
to be of lower order than the other term. Denote $\beta_{j}$ and $\beta_{0,j}$
for component $j$ ($j\in 1,\dots,p$) of $\beta$ and $\beta_{0}$, respectively.
Next, by Hölder’s inequality

$\displaystyle\left|\sqrt{|\mathcal{V}_{k}|}(\beta_{0}-\hat{\beta})\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\frac{\partial}{\partial\beta}\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta)\biggr{\rvert}_{\beta=\beta_{0}}\right]\right|$
$\displaystyle\leq\sqrt{|\mathcal{V}_{k}|}\left|\left|\beta_{0}-\hat{\beta}\right|\right|_{1}\left|\left|\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left[\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\frac{\partial}{\partial\beta}\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta)\biggr{\rvert}_{\beta=\beta_{0}}\right]\right|\right|_{\infty}$
$\displaystyle\leq\sqrt{|\mathcal{V}_{k}|}\left|\left|\beta_{0}-\hat{\beta}\right|\right|_{1}\max_{j\leq
p,i\in\mathcal{V}_{k}}\left|\frac{\partial}{\partial\beta_{j}}\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta)\biggr{\rvert}_{\beta_{j}=\beta_{0,j}}\right|\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left|\left(\frac{1}{p(X_{i};\hat{\beta})}-\frac{1}{p(X_{i};\beta_{0})}\right)\right|$
$\displaystyle\leq\sqrt{|\mathcal{V}_{k}|}\left|\left|\beta_{0}-\hat{\beta}\right|\right|_{1}\max_{j\leq
p,i\in\mathcal{V}_{k}}\left|\frac{\partial}{\partial\beta_{j}}\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta)\biggr{\rvert}_{\beta_{j}=\beta_{0,j}}\right|\left(\max_{i\in\mathcal{V}_{k}}\frac{1}{p(X_{i};\beta_{0})p(X_{i};\hat{\beta})}\right)\frac{1}{|\mathcal{V}_{k}|}\sum_{i\in\mathcal{V}_{k}}\left|p(X_{i};\hat{\beta})-p(X_{i};\beta_{0})\right|$
$\displaystyle\leq\sqrt{|\mathcal{V}_{k}|}\left|\left|\beta_{0}-\hat{\beta}\right|\right|_{1}\max_{j\leq
p,i\in\mathcal{V}_{k}}\left|\frac{\partial}{\partial\beta_{j}}\hat{h}_{1,k}^{(u)}(X_{i};\hat{\epsilon},\beta)\biggr{\rvert}_{\beta_{j}=\beta_{0,j}}\right|\left(\max_{i\in\mathcal{V}_{k}}\frac{1}{p(X_{i};\beta_{0})p(X_{i};\hat{\beta})}\right)\frac{\sqrt{|\mathcal{V}_{k}|}}{|\mathcal{V}_{k}|}\sqrt{\sum_{i\in\mathcal{V}_{k}}\left\\{p(X_{i};\hat{\beta})-p(X_{i};\beta_{0})\right\\}^{2}}$