# Doubly Robust Targeted Estimation of Conditional Average Treatment Effects
for Time-to-event Outcomes with Competing Risks

Runjia Li1 Victor B. Talisa2 and Chung-Chou H. Chang3,∗
1Department of Biostatistics<EMAIL_ADDRESS>University of Pittsburgh
Pittsburgh  Pennsylvania  U.S.A.
2Department of Critical Care  University of Pittsburgh  Pittsburgh
Pennsylvania  U.S.A.
3Department of Medicine  University of Pittsburgh  Pittsburgh  Pennsylvania
U.S.A

###### Abstract

In recent years, precision treatment strategy have gained significant
attention in medical research, particularly for patient care. We propose a
novel framework for estimating conditional average treatment effects (CATE) in
time-to-event data with competing risks, using ICU patients with sepsis as an
illustrative example. Our approach, based on cumulative incidence functions
and targeted maximum likelihood estimation (TMLE), achieves both asymptotic
efficiency and double robustness. The primary contribution of this work lies
in our derivation of the efficient influence function for the targeted causal
parameter, CATE. We established the theoretical proofs for these properties,
and subsequently confirmed them through simulations. Our TMLE framework is
flexible, accommodating various regression and machine learning models, making
it applicable in diverse scenarios. In order to identify variables
contributing to treatment effect heterogeneity and to facilitate accurate
estimation of CATE, we developed two distinct variable importance measures
(VIMs). This work provides a powerful tool for optimizing personalized
treatment strategies, furthering the pursuit of precision medicine.

###### keywords:

Competing Risks; Conditional average treatment effect; Efficient influence
function; Survival Analysis; Targeted maximum likelihood estimation.

##  1 Introduction

In the realm of medical research and patient care, the demand for precision
medicine and personalized treatment approaches is on the rise. This need is
particularly critical when dealing with life-threatening conditions such as
sepsis in intensive care unit (ICU) patients, characterized by time-to-event
outcomes in the presence of competing risks. To estimate the heterogeneous
treatment effect across diverse patient populations and enable tailored
interventions, we need a robust approach for estimating the conditional
average treatment effect (CATE) for time-to-event data in the presence of
competing risks.

The subdistribution hazards model by Fine and Gray (1999) is the most widely
used method for analyzing time-to-event data in the presence of competing
risks. To address issues of collapsibility (Gail et al., 1984), we focus on
the estimand of differences in the cumulative incidence function (CIF) to
infer causality. Unlike an alternative approach based on cause-specific
hazards of all competing events, the subdistribution hazards model directly
connects covariate effects to the CIF, simplifying the interpretation of
covariates impacts on outcomes. Therefore, our targeted estimand becomes a
preferred choice for estimating patient prognosis and predicting CIF (Lau et
al., 2009; Koller et al., 2012; Austin and Fine, 2017).

Existing methods for estimating causal treatment effects in the presence of
competing risks for time-to-event outcomes often fall short due to the
challenges of adjusting for both censoring and confounding bias. Parametric
methods, heavily reliant on distributional assumptions, are at great risk of
bias due to misspecification. Nonparametric models for the outcome or
treatment assignment can result in the slow convergence rates, and/or a lack
of statistical inference. For these reasons, doubly robust estimators, based
on semiparametric efficiency theory, have gained popularity. The estimators
have been developed for binary and continuous outcome variables (Robins et
al., 1994; van der Laan and Rubin, 2006; Kennedy et al., 2017). Doubly robust
estimators provide consistent effect estimates when either the outcome model
or the treatment model is correctly specified, and yield efficient estimates
when both are correctly specified. Among these, targeted maximum likelihood
estimation (TMLE), proposed by van der Laan and Rubin (2006), stands out for
its advantageous properties in finite samples by the nature of being a
substitution estimator, especially when the positivity assumption is nearly
violated. For example, the CATE, defined as the difference between CIFs,
should lie between $-1$ and $1$, which could be ensured by TMLE but not by
non-substitution estimators. TMLE has also been extended to survival outcomes,
with additional work for adjusting for censoring (Moore and van der Laan,
2009; Cai and van der Laan, 2020; Benkeser et al., 2018; Rytgaard et al.,
2023). However, while some studies within the TMLE framework have aimed at CIF
estimation, they rely on estimating cause-specific hazards functions, where
covariate effects do not provide direct interpretations on the effects of the
CIF. Furthermore, these efforts generally focus on average treatment effects,
providing limited guidance for treatment decisions across diverse patient
population.

In this study, we propose a comprehensive framework that leverages TMLE to
consistently estimate CATE in time-to-event data with competing risks. The
primary contribution of this work lies in our derivation of the efficient
influence function (EIF) for the target causal parameter, CATE, which
quantifies the difference in potential outcomes defined by CIFs under distinct
treatment conditions and baseline characteristics. Our proposed estimator
posses doubly robustness and achieves asymptotic efficiency under mild
conditions. We provide theoretical proofs and empirical validation through
simulation studies. Leveraging the empirical variance of EIF, we construct
confidence intervals and conduct hypothesis tests to identify overall
treatment effects and treatment effect heterogeneity. Moreover, our flexible
framework can incorporate various machine learning algorithms for modeling
outcome and treatment mechanisms.

This paper is organized as follows. In Section 2, we define the problem
setting and notation for competing risks data, introducing the causal
parameters of interest. Section 3 reviews the underlying theory and TMLE
framework, illustrates the construction process and the properties of the
proposed estimator. In Section 4, we demonstrate the practical performance and
properties through simulations across scenarios and settings. In Section 5,
the proposed method is applied to data from ICU patients with sepsis to
estimate the CATE of steroids. Finally, in Section 6, we provide a discussion
of the findings.

##  2 Setting and Notation

###  2.1 Time-to-event data with competing risks

We consider a data setting that involves time-to-event outcomes in the
presence of competing risks and encompasses two treatment arms. Let the
observed data be denoted as
$O_{i}=\left(X_{i},A_{i},\widetilde{T}_{i},\widetilde{\Delta}_{i}\right)_{i=1}^{N}\sim^{i.i.d}P$
for subjects $i=1,\ldots,N$, where $X\in\mathcal{X}\subseteq\mathbb{R}^{p}$
represents a vector of baseline covariates, $A\in\\{0,1\\}$ is a binary
treatment assignment, $T$ is the failure time, and $C$ is the right-censoring
time. $P$ denotes the unknown true probability distribution of $O$, where our
observed sample drawn from. Consider $\widetilde{T}=T\wedge C$ as the observed
event time, $\Delta\in\\{1,2\\}$ as the event type, and
$\widetilde{\Delta}=\mathbbm{1}(T\leq C)\Delta$ as the observed event type.
Without loss of generality, let event 1 ($\Delta=1$) be the primary event of
interest, and event 2 ($\Delta=2$) represent the competing events. The
cumulative incidence function (CIF) for the primary event at time $t_{0}$ is
defined as $F_{1}(t_{0})=\operatorname{Pr}(T\leq t_{0},\Delta=1)$. Let $T^{a}$
and $\Delta^{a}$ be the potential event time and event type, respectively, had
the subject been assigned treatment $A=a$ $(a=0,1)$. Then, the potential CIF
for the primary event at time $t_{0}$ had the subject received treatment $a$
is $F_{1}^{a}(t_{0})=\operatorname{Pr}(T^{a}\leq t_{0},\Delta^{a}=1)$.

###  2.2 Causal estimand

Our parameter of interest is the conditional average treatment effect (CATE)
in prespecified subgroups by a set of predictive variables
$V\in\mathcal{V}\subset\mathcal{X}$, where treatment effects differ among
patients with various values of $V$. $V$ are usually selected based on
clinical knowledge, interests, or by a data-driven approach, e.g., by
assessing the significance of interaction terms between treatment and
candidate covariates through a preliminary full model of subdistribution
hazards model, which include all covariates and their interactions with
treatment. The CATE is the average treatment effect conditional on $V=v_{m}$
within each subgroup $m$ $(m=1,2,...,M)$, defined as the difference in
potential CIF of main event:

$\displaystyle\Psi_{t_{0},m}(P)$
$\displaystyle=F^{1}_{1,m}(t_{0}|V=v_{m})-F^{0}_{1,m}(t_{0}|V=v_{m})$
$\displaystyle=E\left\\{F^{1}_{1}(t_{0}|V=v_{m},L)-F^{0}_{1}(t_{0}|V=v_{m},L)|V=v_{m}\right\\}$

where $F^{1}_{1,m}(t_{0}|V=v_{m})$ is the potential CIF had the subject
received treatment $a$, possibly contrary to fact, conditional on the
prognostic variables
$L\in\mathcal{L}\subseteq\mathcal{X}\backslash\mathcal{V}$ and predictive
variables $V=v_{m}$. When $\mathcal{V}=\varnothing$, $M=1$ and
$\Psi_{t_{0}}(P)$ is the average treatment effect.

Within the potential outcome framework by Rosenbaum and Rubin (1983), this
parameter can be causally identified from observed data under the following
assumptions: 1) Consistency: $T=T^{a}$, $C=C^{a}$ and $\Delta=\Delta^{a}$; 2)
No unmeasured confounder: $T^{a},\Delta^{a}\perp\\!\\!\\!\perp{A}|V,L$; 3)
Positivity: $0<\operatorname{Pr}(A=1|V,L)<1$ and
$\operatorname{Pr}(C>t_{0}|V,L,A)>0$ with probability 1; and 4) Independent
censoring: $T,\Delta\perp\\!\\!\\!\perp{C}|V,L$. In addition, we assume that
the potential outcomes are not affected by the treatment assignments of other
subjects, and only one version of the treatment (known as the stable unit
treatment value assumption, SUTVA). The target causal parameter within each
subgroup can be then identified as the following statistical parameter:

$\Psi_{t_{0},m}(P)=E\left\\{F_{1}(t_{0}|A=1,L,V=v_{m})-F_{1}(t_{0}|A=0,L,V=v_{m})|V=v_{m}\right\\}.$
(1)

For simplicity, the $\Psi_{t_{0},m}(P)$ is denoted as $\Psi(P)$, and
conditioning on $V=v_{m}$ and subscript $m$ are dropped:

$\Psi(P)=E\left\\{F_{1}(t_{0}|A=1,L)-F_{1}(t_{0}|A=0,L)\right\\}.$ (2)

##  3 Targeted MLE Framework

We developed an estimation framework to construct the targeted estimator
$\Psi(\hat{P}_{n}^{*})$ with doubly robust feature. This framework is based on
the targeted maximum likelihood estimation (TMLE) introduced by van der Laan
and Rubin (2006), along with methods extended to estimating average treatment
effect for time-to-event outcomes (Rytgaard et al., 2023; Moore and van der
Laan, 2009; Cai and van der Laan, 2020). When subjects are stratified
according to a predetermined set of predictive covariates, our estimation
procedure follows these steps: 1) Initial estimation: We begin by generating
initial estimates for the nuisance parameters; 2) Targeting step: The
estimators are updated iteratively to solve the corresponding efficient
influence function (EIF) equations, thereby ensuring asymptotic properties;
and 3) CATE estimation: Finally, we compute the CATE by averaging the updated
estimates on the baseline covariates that are not used for stratification. We
will delve into the details in the following subsections. The step-by-step
procedure for estimating CATE is summarized in Supplementary Materials, Web
Appendix A Algorithm 1.

###  3.1 Initial estimation

To solve the EIF equations in the targeting step, we need initial estimates
for the nuisance functions: (i) CIF of the main event
$F_{1}(t|A,L)=1-\exp\\{-\int_{0}^{t}\lambda_{1}(s|A,L)ds\\}$, where
$\lambda_{1}(s|A,L)$ is the subdistribution hazard function of the main event;
(ii) survival function for censoring $G(t|A,L)=\operatorname{Pr}(C>t|A,L)$;
and (iii) the treatment mechanism $\pi(a|L)=\operatorname{Pr}(A=a|L)$.

To estimate CIF, let’s start with the subdistribution hazard function (Gray,
1988)

$\displaystyle\lambda_{1}\left(t|A,L\right)=\lim_{dt\rightarrow
0}\frac{1}{dt}\operatorname{Pr}$ $\displaystyle\\{t\leq T<t+dt,{\Delta}=1\mid$
(3) $\displaystyle T\geq t\cup(T\leq t\cap{\Delta}\neq 1),A,L\\},$

where the risk set includes subjects who are currently free of any event and
those who have previously experienced a competing event. Because there is a
one-to-one relationship between the subdistribution hazard function and the
CIF, we chose to use the subdistribution hazards model proposed by Fine and
Gray (1999) to predict CIF. To simultaneously estimate unknown parameters and
perform variable selection, researchers developed penalized subdistribution
hazards model and related algorithms (Fu et al., 2017). Kawaguchi et al.
(2020) enhanced computational efficient estimators in R package fastcmprsk. We
adopted the algorithm in our parameter estimation and variable selection
procedure.

The conditional survival function for censoring events can be estimated using
various methods, including the Kaplan-Meier estimator or regression models
like the Cox proportional hazards model. Similarly, the treatment mechanism
can be estimated using a range of statistical techniques for binary outcome,
such as logistic regression and generalized additive model (GAM) (van der Laan
et al., 2011).

###  3.2 Targeting step: iterative updates

To achieve asymptotic efficiency and doubly robust properties for the
estimator, the key step in TMLE involves refining the initial estimates from
Step 1 by deriving the EIF expression and solving the estimating equation. We
will first explain the underlying theory and then illustrate the process of
obtaining the asymptotically efficient and doubly robust estimator.

####  3.2.1 Efficient influence function (EIF)

The EIF has been widely used in constructing efficient and doubly robust
estimators (Robins et al., 1994; van der Laan and Rubin, 2006; Kennedy et al.,
2017) based on the semiparametric and empirical process theories (Bickel et
al., 1993; van der Vaart, 1998). Under regularity conditions, the targeted
estimator $\Psi(\hat{P}_{n}^{*})$ is asymptotically linear if

$\sqrt{n}\\{\Psi(\hat{P}_{n}^{*})-\Psi(P)\\}=\sqrt{n}\mathbb{P}_{n}\\{D(P)\\}+o_{P}(1),$
(4)

where $D(P)$ is an influence function of the estimator, $\mathbb{P}_{n}$
denote the empirical average that
$\mathbb{P}_{n}(\cdot)=\frac{1}{n}\sum_{i=1}^{n}(\cdot)$, $E\\{D(P)\\}=0$, and
$E\\{D(P)^{2}\\}<\infty$. A regular estimator in a semiparametric model is
asymptotically efficient with the lowest asymptotic variance if it is
asymptotically linear with $D(P)$ being the EIF (Bickel et al., 1993). To
obtain the asymptotically efficient estimator and calculate the lower bound of
the variance, we need to: 1) derive the expression of $D(P)$, and 2) find the
targeted distribution
$\hat{P}_{n}^{*}=(\hat{\lambda}_{1}^{*},\hat{\pi},\hat{G}$) that solving the
efficient influence function equation:

$\mathbb{P}_{n}\\{D(\hat{P}_{n}^{*})\\}=0,$ (5)

where $\hat{\lambda}^{*}$ is the updated estimated subdistribution hazard
function.

A general approach for deriving the EIF is to calculate the pathwise
derivative $\frac{d}{d\epsilon}\Psi(P_{\epsilon})\Bigr{|}_{\epsilon=0}$ for a
submodel $P_{\epsilon}$. Different from conventional concept of parametric
model, a submodel is a useful conceptual tool in semiparametric theory,
instead of being used for analysis (Tsiatis, 2006). It is a smooth parametric
model that is contained in the larger model of interest and is equal to the
true distribution when evaluated at $\epsilon=0$. A commonly chosen submodel
has the form: $dP_{\epsilon}=(1+\epsilon b)dP$ (Kennedy, 2022), where
$\epsilon\in\mathbb{R}$ is the index for perturbation; $b=b(T,\Delta,A,L)$ is
a zero-mean score function and for this submodel,
$b=\frac{d}{d\epsilon}\log(dP_{\epsilon})\Bigr{|}_{\epsilon=0}$. The EIF can
be obtained by expressing the pathwise derivative of the target estimand as an
inner product of two terms:
$\frac{d}{d\epsilon}\Psi(P_{\epsilon})\Bigr{|}_{\epsilon=0}=E\\{D(P)\cdot
b\\}$. A detailed description is provided by Levy (2019). Per the definition
of CIF, our target parameter can be expressed as:

$\displaystyle\Psi(P)$
$\displaystyle=E\left\\{F_{1}(t_{0}|A=1,L)-F_{1}(t_{0}|A=0,L)\right\\}$ (6)
$\displaystyle=\int_{\mathcal{L}}\sum_{a=0,1}\left[1-\exp\left\\{-\int_{0}^{t_{0}}\Lambda_{1}(dt|a,l)\right\\}\right]\left\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\right\\}\mu(l)dl,$

where $\mu(l)$ is the density of prognostic variables $L$.

Let $P^{\delta}(dt|a,l)=\operatorname{Pr}(t\leq{T}<t+dt,\Delta=\delta|a,l)$,
then from the definition of subdistribution hazard function (3) we have
$\Lambda_{1,\epsilon}(dt|a,l)=P_{\epsilon}^{1}(dt|a,l)\\{P_{\epsilon}({T}\geq
t\cup({T}\leq t\cap{\Delta}\neq 1)|a,l)\\}^{-1}$. The derivation relies on the
differentiation of $P^{\delta}_{\epsilon}$ with respect to $\epsilon$
evaluated at $\epsilon=0$, i.e.,
$\frac{d}{d\epsilon}P^{\delta}_{\epsilon}(dt|a,l)\Bigr{|}_{\epsilon=0}=\frac{d}{d\epsilon}P^{\delta}(dt|a,l)(1+\epsilon
b)\Bigr{|}_{\epsilon=0}=P^{\delta}(dt|a,l)b(t,\delta,a,l)$.
Theorem 1. (Efficient influence function of $\Psi(P)$). Under Assumptions 1 -
4 in Section 2.2, the efficient influence function is in the form of:

$\displaystyle D(P)=$
$\displaystyle\int_{0}^{t_{0}}\left\\{h(t,A,L)w(t)M(dt|A,L)\right\\}+\left\\{F_{1}(t_{0}|1,L)-F_{1}(t_{0}|0,L)\right\\}-\Psi(P),$
(7)

where

$\displaystyle h(t,A,L)=$
$\displaystyle\frac{\\{\mathbbm{1}(A=1)-\mathbbm{1}(A=0)\\}}{\pi(A|L)G(t^{-}|A,L)}\cdot\frac{1-F_{1}\left(t_{0}|A,L\right)}{1-F_{1}(t|A,L)}$
(8)

is defined as the clever covariate, which plays a pivotal role in the
fluctuation model for estimating CATE under the TMLE framework (will be
described in detail in next subsection); $w(t)=\mathbbm{1}(C\geq T\wedge
t)G(t^{-}|A,L)/G\\{(\tilde{T}\wedge t)^{-}|A,L\\}$ is the weight determined by
adapting the inverse probability of censoring weighting (IPCW) technique (Fine
and Gray, 1999); and $M(dt|A,L)=N(dt)-Y(t)\Lambda(dt|A,L)$, where
$N(t)=\mathbbm{1}(T\leq t,\Delta=1)$ and $Y(t)=\\{1-N(t^{-})\\}$. Derivation
of EIF is in Supplementary Materials, Web Appendix B.

####  3.2.2 Solving the efficient influence function equation

The targeted estimator $\Psi(\hat{P}_{n}^{*})$ is obtained by solving the
efficient influence function equation (5) via a fluctuation model. The
fluctuation model is selected to have its corresponding mean score function
equal to the empirical mean of first part of the EIF expressed in (7) when the
fluctuation parameter $\epsilon=0$. In our case, we consider the following
fluctuation model:

$\lambda_{1,\epsilon}\left(t|A,L\right)=\lambda_{1}\left(t|A,L\right)\exp\left\\{\epsilon
h\left(t,A,L\right)\right\\},$ (9)

where $h\left(t,A,L\right)$ is the clever covariate defined in (8), and
$\epsilon\in\mathbb{R}$ is the fluctuation parameter. The score function with
respect to this model is:

$\displaystyle\frac{d\log\mathscr{L}\left(\lambda_{1,\epsilon}\right)}{d\epsilon}$
$\displaystyle=\sum_{i=1}^{n}\int_{0}^{t_{0}}w_{i}(t)h(t,A_{i},L_{i})\left[N_{i}(dt)-Y_{i}(t)\lambda_{1}(t|A_{i},L_{i})\exp\\{\epsilon
h(t,A_{i},L_{i})\\}dt\right],$ (10)

where $\mathscr{L}$ is the likelihood function of the model. We obtain the
maximum likelihood estimator $\hat{\epsilon}^{(k)}$ by solving
$\frac{d\log\mathscr{L}\left(\lambda_{1,\epsilon}\right)}{d\epsilon}=0$. At
each iteration (index $k=0,1,2,...$), we use $\hat{\epsilon}^{(k)}$ to update
$\hat{P}_{n}$, specifically,
$\hat{\lambda}_{1}^{(k+1)}(t|A,L)=\hat{\lambda}_{1}^{(k)}(t|A,L)\exp\left\\{\hat{\epsilon}^{(k)}h(t,A,L)\right\\}$,
and
$\hat{F}_{1}^{(k+1)}\left(t|A,L\right)=1-\exp\left\\{-\int_{0}^{t}\hat{\lambda}_{1}^{(k+1)}\left(t|A,L\right)dt\right\\}$.
It will iterate until $\hat{\epsilon}$ becomes very small and approaches 0.
Consequently, the mean score function (10) converges to the empirical mean of
the first part of the EIF (7):

$\frac{1}{n}\frac{d\log\mathscr{L}\left(\lambda_{1,\epsilon}\right)}{d\epsilon}\Bigr{|}_{\epsilon=0}=\mathbb{P}_{n}\left[\int_{0}^{t_{0}}\left\\{h(t,A,L)w(t)M(dt|A,L)\right\\}\right].$
(11)

Note that the empirical mean of the remaining part of the EIF is zero by
definition,
$\mathbb{P}_{n}[\\{\hat{F}_{1}^{*}(t_{0}|1,L)-\hat{F}_{1}^{*}(t_{0}|0,L)\\}-\Psi(\hat{P}_{n}^{*})]=0$.
The property of (11) enables us to solve the efficient influence function
equation (5) around $\epsilon=0$ and obtain the solution $\hat{P}_{n}^{*}$.
For implementation, we can set a small number $s_{n}=o_{p}(n^{-\frac{1}{2}})$
and stop the iteration when $\left|\epsilon\right|\leq s_{n}$. Derivation of
the score function can be found in Supplementary Materials, Web Appendix C.
Subsequently, we can substitute this solution to obtain the targeted estimator
$\Psi(\hat{P}_{n}^{*})$, which is introduced in the next subsection.

###  3.3 Asymptotic properties

With the estimator $\hat{P}_{n}^{*}$ being the solution of the efficient
influence function equation, $\Psi(\hat{P}_{n}^{*})$ is asymptotically linear
with respect to $D(P)$ if following conditions are satisfied: 1) the second-
order remainder
$R(\hat{P}_{n}^{*},P)=\Psi(\hat{P}_{n}^{*})-\Psi(P)+E_{P}\\{D(\hat{P}_{n}^{*})\\}$
is of the order $o_{p}(n^{-1/2}$); and 2) $D(\hat{P}_{n}^{*})$ belongs to the
Donsker class (van der Vaart, 1998; van der Laan and Rubin, 2006; van der Laan
et al., 2011).

For the first condition, the second-order remainder can be decomposed as

$\displaystyle R(\hat{P}_{n}^{*},P)=\int_{\mathcal{L}}\sum_{a=0,1}$
$\displaystyle\int_{0}^{t_{0}}\left[\frac{\pi(a|l)}{\hat{\pi}(a|l)}\frac{G\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}{\hat{G}\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}-1\right]\cdot\left\\{\Lambda_{1}(dt|a,l)-\hat{\Lambda}_{1}^{*}(dt|a,l)\right\\}$
(12)
$\displaystyle\cdot\left\\{\frac{1-\hat{F}_{1}^{*}(t_{0}|a,l)}{1-\hat{F}_{1}^{*}(t|a,l)}\right\\}\left\\{1-F_{1}(t^{-}|a,l)\right\\}\cdot\left\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\right\\}\mu(l)dl,$

where the double-robust structure is demonstrated. The remainder is equal to
$o_{p}(n^{-1/2})$ when either both the treatment and censoring mechanisms are
consistently estimated or the CIF is consistently estimated, i.e., either
$||\\{\pi(a|l)G(t|a,l)\\}-\\{\hat{\pi}(a|l){\hat{G}(t|a,l)}\\}||=o_{p}(n^{-1/2})$,
or $||\Lambda_{1}(dt|a,l)-\hat{\Lambda}_{1}^{*}(dt|a,l)||=o_{p}(n^{-1/2})$.
$R(\hat{P}_{n}^{*},P)=o_{p}(n^{-1/2})$ can also be achieved if both
$\hat{\pi}_{0}(a|l)\hat{G}(t|a,l)$ and $\hat{\lambda}_{1}^{*}(t|a,l)$ converge
at a rate faster than $o_{P}(n^{-1/4})$. See Supplementary Materials, Web
Appendix D for details.

For the second condition, the nuisance functions $\hat{F}_{1}$,
$\hat{\lambda}_{1}$, $\hat{G}$, and $\hat{\pi}$ of $D(\hat{P}_{n})$ , which
are all components of the EIF function, belong to the Donsker class. These
functions can be parameterized by the Càdlàg functions, which are right-
continuous with left limits, ensuring that the corresponding estimators are
not too complex. Because all components of the EIF function belong to the
Donsker class, the EIF itself also falls within the Donsker class (van der
Vaart, 1998; van der Laan and Rubin, 2006; Rytgaard et al., 2023).

When these two conditions are met, the targeted substitution estimator
$\Psi(\hat{P}^{*}_{n})$ exhibits asymptotic linearity with respect to the EIF
$D(P)$. Moreover, by the Central Limit Theorem, $\Psi(\hat{P}^{*}_{n})$
possesses asymptotic normality:

$\sqrt{n}\left\\{\Psi(\hat{P}_{n}^{*})-\Psi(P)\right\\}\xrightarrow{\mathcal{D}}N\left[0,E_{P}\\{D(P)^{2}\\}\right].$

The asymptotic variance of this estimator can be estimated by
$\hat{\sigma}^{2}=\frac{1}{n}\sum_{i=1}^{n}\\{D(\hat{P}_{n}^{*})\\}^{2}$. The
Wald confidence interval can be constructed by $\Psi(\hat{P}_{n}^{*})\pm
z_{1-\alpha/2}\sqrt{\mathbb{P}_{n}\\{D(\hat{P}_{n}^{*})\\}^{2}}$, where
$z_{1-\alpha/2}$ is the $(1-\alpha/2)^{th}$ percentile of the standard normal
distribution.

###  3.4 Variable importance

To investigate the key predictive variables contributing to treatment effect
heterogeneity, as well as the importance of prognostic variables for treatment
effect estimation within each subgroup, we calculate two variable importance
measures (VIM) after constructing the estimator.

To evaluate the variable importance of predictive variables in terms of the
contribution to the treatment effect heterogeneity, our VIM is calculated
based on the variance of the estimated treatment effect, inspired by Hines et
al. (2022). For $j^{th}$ predictive variable ($j=1,2,...,d,d<p$, recall that
$V\in\mathcal{V}\subset\mathcal{X}\subseteq\mathbb{R}^{p}$):

$\displaystyle\text{VIM}_{1,j}$
$\displaystyle=\frac{\Bigr{|}\text{var}\left\\{\hat{\tau}(V_{i},L_{i})\right\\}-\text{var}\left\\{\hat{\tau}(V_{i,-j},L_{i})\right\\}\Bigr{|}}{\text{var}\\{\hat{\tau}(V_{i},L_{i})\\}},$
(13)

where $\hat{\tau}(V_{i},L_{i})$ is the estimated treatment effect of
individual $i$ using all predictive and prognostic variables, while
$\hat{\tau}(V_{i,-j},L_{i})$ is the one with $j^{th}$ predictive variable
excluded. We calculate ${\text{VIM}}_{1,j}$ as the difference in variance of
the estimated treatment effect when the variable is removed.

Additionally, to guide physicians in identifying critical variables for
treatment effect estimation, we estimate VIM for each prognostic variable.
This measure, akin to the one proposed by Pirracchio et al. (2018), evaluates
the change in estimated CATE within each subgroup when one prognostic variable
is omitted at a time:

${\text{VIM}}_{2,k}=\Psi(\hat{P}_{n}^{*})-\Psi(\hat{P}^{*}_{n,-k}),$ (14)

where $\Psi(\hat{P}_{n}^{*})$ and $\Psi(\hat{P}^{*}_{n,-k})$ represent the
estimated treatment effects using all variables and using all variables except
the $k^{th}$ prognostic variable ($k=1,2,...q,q\leq(p-d)$, recall that
$L\in\mathcal{L}\subseteq\mathcal{X}\backslash\mathcal{V}$), respectively.

##  4 Simulation studies

###  4.1 Simulation setup

To assess the performance of our proposed TMLE framework and illustrate the
characteristics of the targeted estimator, we conducted simulation studies
across various scenarios in R (code are available on Github as the link
provided in Supplementary Materials). We generated a two-cause competing risks
data following a scheme similar to Fine and Gray (1999), simulating observed
data akin to sepsis patients.

The baseline predictive covariates $V_{1}$, $V_{2}$, and the prognostic
covariate $L_{1}$ were independently generated from $Ber(0.5)$, and the other
three prognostic covariates $L_{2}$, $L_{3}$, and $L_{4}$, were independently
generated from $N(0,1)$. The binary treatment variable $A$ was generated from
$Ber([1+\exp\\{-(-0.2V_{1}-0.1V_{2}+1.5L_{1}+0.1L_{3})\\}]^{-1})$. To generate
the main event time $T_{1}$, we used a unit exponential mixture with the
failure function taking the form:
$F_{1}(t|A,V,L)=1-\\{1-p(1-\exp(-t))\\}^{\exp\\{Y_{1}(A,V,L)\\}}$, where
$Y_{1}(A,V,L)=(0.2V_{1}-0.9L_{1}-0.1L_{2})+(-0.8V_{1}+0.6V_{2})A+0.5(A-0.5)$,
$p=0.7$. The competing event time $T_{2}$ was generated from exponential
distribution with mean $[\exp\\{Y_{2}(A,V,L)\\}]^{-1}$, where
$Y_{2}(A,V,L)=0.5Y_{1}(A,V,L)$. The event indicator $\Delta\sim
1+Ber\\{(1-p)^{Y_{1}(A,V,L)}\\}$. Finally, the noninformative censoring times
$C$ were generated from exponential distribution with mean
$[\lambda_{0}\exp\\{0.1V_{1}-0.2V_{2}-0.1L_{1}+0.05L_{4}\\}]^{-1}$, resulting
in approximately 42% for the main event, 33% for the competing event, and 25%
for noninformative censoring. The interaction terms between treatment $A$ and
covariates $V_{1}$ and $V_{2}$ in $Y_{1}$ played a critical role in
introducing treatment effect heterogeneity across the simulation scenarios.

To investigate the double robustness and asymptotic properties of our proposed
estimator, we included scenarios with correctly specified and misspecified
models for the outcome, treatment and censoring. In the misspecified outcome
models, $L_{1}$ and $L2$ were deliberately excluded, with additional
covariates $L_{3}$ and $L_{4}$ introduced in the subdistribution hazard model.
The misspecified treatment models excluded $L_{1}$ and $L3$, and introduced
additional covariates, $L_{2}$ and $L_{4}$. Similarly, for the misspcification
in censoring models, we deliberately omitted $L_{1}$ and $L4$, and extra
covariates, $L_{2}$ and $L3$, were introduced. To provide a benchmark, we
implemented existing meta-algorithms for CATE, specifically the S-learner and
T-learner as described in (Künzel et al., 2019), and compared their
performance to our proposed TMLE framework. The S-learner utilizes a single
model to predict the outcome, with treatment as a covariate and including its
interaction terms with covariates, and calculates the difference in
predictions between different treatment assignments. Conversely, the T-learner
uses two separate models, one for treated and one for untreated individuals,
and subsequently computes the difference in predictions. The S-learner and
T-learner were also incorporated in the first step in our TMLE framework to
generate the initial estimates of the outcome. Furthermore, our simulations
included adjustments to the sample size ($n$) with values of 800, 1,500, and
3,000, as well as the time points of interests at the 25th, 50th and 75th
percentile event time.

In addition to evaluate the properties of our proposed estimator and compare
with existing methods, we estimated CATE by incorporating various combinations
of models in our TMLE framework, to investigate the performance in high-
dimensional data. For the outcome model, we employed subdistribution hazards
model or LASSO-regularized subdistribution hazards model, implemented using
the R package fastcmprsk (Kawaguchi et al., 2020). The treatment model was
built using logistic regression or a super-learner that encompasses
generalized linear models, generalized additive models, and LASSO-regularized
logistic regression using the R package SuperLearner (Van der Laan et al.,
2007). The censoring model was built using Cox model, or a super-learner that
encompasses Cox model, exponential model, log-logistic model, and LASSO-
regularized Cox model using the R package survSuperLearner (Golmakani and
Polley, 2020). We applied these methods in another set of simulated data to
mimic the situation with high dimensional data with large amount of noises.
This set of data are generated by adding more prognostic covariates and noise
variables, which could be selected in the model but are not involved in the
true models. The number of prognositic variables was increased to 28, where
$L_{1},L_{2},...,L_{8}\sim Ber(0.5)$ and $L_{9},L_{10},...,L_{28}\sim N(0,1)$.
The data generating models are updated with two additional variables, $L_{9}$
and $L_{10}$: $A\sim
Ber([1+\exp\\{-(-0.2V_{1}-0.1V_{2}+1.5L_{1}-0.1L_{2}+0.4L_{3}-0.1L_{9}-0.1L_{10})\\}]^{-1})$,
$Y_{1}(A,V,L)=(0.2V_{1}-0.9L_{1}-0.15L_{2}+0.2L_{9}-0.2L_{10})+(-0.8V_{1}+0.6V_{2})A+0.5(A-0.5)$,
and the censoring times $C$ were generated from exponential distribution with
mean
$[\lambda_{0}\exp\\{0.1V_{1}-0.2V_{2}-0.1L_{1}+0.1L_{9}-0.1L_{10}\\}]^{-1}$.

We conducted $B=500$ simulations under the aforementioned settings. Estimator
is constructed for each one of the four subgroup defined by $V_{1}$ and
$V_{2}$: group $m=1$ with $V_{1}=0$ and $V_{2}=0$, group $m=2$ with $V_{1}=0$
and $V_{2}=1$, group $m=3$ with $V_{1}=1$ and $V_{2}=0$, and group $m=4$ with
$V_{1}=1$ and $V_{2}=1$. Bias and root-mean-square error (RMSE) were computed
for each subgroup ($m=1,2,3,4$) defined by $V_{1}$ and $V_{2}$ as
$\text{Bias}_{m}(t_{0})=\Psi_{m}(\hat{P}^{*}_{n})-\Psi_{m}(P)\text{ and RMSE
}_{m}(t_{0})=\sqrt{\frac{1}{B}\sum^{B}_{b=1}\\{\Psi_{m}(\hat{P}^{*}_{n})-\Psi_{m}(P)\\}^{2}}$,
where $\Psi_{m}(\hat{P}^{*}_{n})$ and $\Psi_{m}(P)$ represent the estimated
and true CATE, respectively for subgroup $m$ at time $t_{0}$. Additionally, we
calculate the 95% confidence interval (CI) for the estimated CATE and the
proportion of true CATE covered by the 95% CI. Hypothesis tests are conducted
to make inference about whether there are significant treatment effects:
$H_{0}:\Psi_{m}(P)=0$.

###  4.2 Simulation results

Our simulation results confirm asymptotic normality and doubly robustness of
our CATE estimator. Figure 1 shows that our TMLE CATE estimator outperforms S-
and T-learners, especially when outcome or treatment models are misspecified,
with lower bias and RMSE.

Figure 1: Bias (upper panels) and RMSE (lower panels) of CATE estimations at
the median follow-up time point under different model specifications. Each
panel corresponds to one of the four subgroups determined by two binary
predictive variables, $V_{1}$ and $V_{2}$. The red line represents zero bias.
$\checkmark$ and $\times$ indicate that the corresponding models were
correctly specified and misspecified, respectively. Light grey and dark grey
represent S-learner and T-learner, respectively, while light blue and dark
blue represent TMLE with S-learner and T-learner employed for the initial
subdistribution hazard model in the first step, respectively. Treatment and
censoring models are not involved in S-learner and T-learner. Sample size is
$n=3000$ and each scenario was replicated 500 times.

In Scenario 1, where all models are correctly specified, all methods perform
well with minimal bias and RMSE. However, in Scenario 2, with a misspecified
outcome model, S- and T-learners exhibit a 200% increase in bias and RMSE
compared to Scenario 1. Remarkably, our TMLE framework maintains bias and RMSE
similar to Scenario 1. This resilience stems from the correction applied in
the targeting step of the TMLE framework, effectively mitigating the impact of
outcome model misspecification. Similarly, in Scenarios 3, 4 and 5, where the
treatment model and/or the censoring model are misspecified, respectively, our
TMLE methods show no inflation in bias and RMSE compared to Scenarios 1 and 2
(S- and T- learner do not involve treatment model so the results are not
reported for Scenarios 3 and 4). We also evaluated the bias and RMSE of the
same scenarios in data with smaller sample size $n=1500$ (Web Figure 1) and
$n=800$ (Web Figure 2), as well as the results at different time points, see
details in Supplementary Materials, Web Appendix E. As anticipated, as sample
sizes decrease, bias and RMSE increase across all methods and scenarios, and
our methods still perform better and remain stable despite misspecification.
Maintaining low bias and RMSE reflect the double robustness of our method that
can protect against misspecification for outcome model, or for treatment and
censoring models, indicating the valuable advantages of our method in
estimating CATE and guiding treatment decisions.

The percentage of true CATE covered by the 95% confidence interval (CI)
constructed using our proposed method is reported in Table 1. In all
scenarios, whether with or without model misspecification, our method
consistently achieves approximately 95% coverage of the true CATE. These CIs,
established by the asymptotic normality of our proposed estimator, demonstrate
their reliability by providing valuable information and facilitating accurate
inference.

Table 1: Coverage probability of CATE estimations at the median follow-up time
point across various sample sizes and simulation scenariosa.

Subgroupb | $\begin{array}[]{c}\textbf{Scenario 1${}^{c}$}\\\ \text{Outcome \checkmark}\\\ \text{Treatment \checkmark}\\\ \text{Censoring \checkmark}\end{array}$ | $\begin{array}[]{c}\textbf{Scenario 2${}^{c}$}\\\ \text{Outcome }\times\\\ \text{Treatment }\checkmark\\\ \text{Censoring }\checkmark\end{array}$ | $\begin{array}[]{c}\textbf{Scenario 3${}^{c}$}\\\ \text{Outcome }\checkmark\\\ \text{Treatment }\times\\\ \text{Censoring \checkmark}\end{array}$ | $\begin{array}[]{c}\textbf{Scenario 4${}^{c}$}\\\ \text{Outcome }\checkmark\\\ \text{Treatment }\checkmark\\\ \text{Censoring }\times\end{array}$ | $\begin{array}[]{c}\textbf{Scenario 5${}^{c}$}\\\ \text{Outcome }\checkmark\\\ \text{Treatment }\times\\\ \text{Censoring }\times\end{array}$
---|---|---|---|---|---
$n=3,000$
$V_{1}=0,V_{2}=0$ | 95.2% | 95.8% | 96.4% | 95.0% | 96.4%
$V_{1}=0,V_{2}=1$ | 96.0% | 95.8% | 96.6% | 95.6% | 96.6%
$V_{1}=1,V_{2}=0$ | 96.6% | 97.0% | 96.4% | 96.8% | 96.4%
$V_{1}=1,V_{2}=1$ | 96.4% | 96.8% | 95.2% | 96.4% | 95.2%
$n=1,500$
$V_{1}=0,V_{2}=0$ | 95.2% | 95.6% | 94.6% | 95.4% | 94.6%
$V_{1}=0,V_{2}=1$ | 94.4% | 94.8% | 96.8% | 94.4% | 96.8%
$V_{1}=1,V_{2}=0$ | 95.8% | 96.0% | 95.4% | 95.8% | 95.4%
$V_{1}=1,V_{2}=1$ | 95.0% | 96.2% | 95.4% | 95.0% | 95.4%
$n=800$
$V_{1}=0,V_{2}=0$ | 93.4% | 94.0% | 95.6% | 93.8% | 95.6%
$V_{1}=0,V_{2}=1$ | 94.2% | 95.2% | 97.4% | 94.8% | 97.4%
$V_{1}=1,V_{2}=0$ | 94.4% | 94.8% | 94.6% | 94.4% | 94.6%
$V_{1}=1,V_{2}=1$ | 94.4% | 94.0% | 95.0% | 94.4% | 95.0%

a$B=500$ Monte-Carlo repetitions were implemented.
b Subgroups were defined by the two binary predictive variables $V_{1}$ and
$V_{2}$.
c$\checkmark$ indicates that the corresponding model was correctly specified,
while $\times$ indicates that the model was misspecified.

Figure 2 displays the bias and RMSE of our proposed estimator in higher-
dimensional data, and compares the results of employing various models for the
outcome, treatment and censoring. As shown in the upper panel, our proposed
estimator provides unbiased estimates no matter which model was employed. The
number of outliers (upper panel) and the RMSE (lower panel) are decreased by
utilizing super-learner methods for treatment and censoring. Using LASSO
regularized subdistribution hazards model for CIF do not significantly lower
bias or RMSE.

Figure 2: Bias (upper panels) and RMSE (lower panels) of CATE estimations at
the median follow-up time point in high-dimensional data with 30 covariates.
Each panel corresponds to one of the four subgroups determined by two binary
predictive variables, $V_{1}$ and $V_{2}$. The red line represents zero bias.
”+S” (light colors) and ”+T” (dark colors) represent TMLE with S-learner and
T-learner employed for the initial subdistribution hazard model in the first
step, respectively; ”+LASSO” represents using LASSO penalized subdistribution
hazards model for the outcome; ”+SL” represents utilizing super learner for
treatment and censoring. Sample size is $n=3000$ and each scenario was
replicated 500 times.

For medical decision-making, relying on S- or T-learners might lead to
incorrect treatment recommendations. Our TMLE framework ensures more reliable
results, guarding against misspecification and inappropriate treatment
decisions based on inconsistent findings.

##  5 Analysis of Sepsis Dataset

###  5.1 Data source and analysis

We analyzed electronic medical records data from 2012 to 2018, encompassing
all ICU patients with sepsis in a healthcare system (N=60,886). We included
3,245 patients who met the eligibility criteria of a recent landmark clinical
trial of steroids in septic shock, ADRENAL (Venkatesh et al., 2018), focusing
on short-term steroid effects. Our primary aim was to assess the impact of
steroids on reducing ICU mortality among various sepsis patients. To focus on
short-term treatment effects and mitigate potential violations of the
positivity assumption, we considered follow-up was capped at 10 days post-ICU
admission, with discharge by this time considered a competing event. We also
identified important biomarkers and indices for treatment strategy
development.

Within the cohort, 885 patients (27.3%) died within 10 ICU days, with 274
(31.0%) receiving steroids. Meanwhile, 1,557 (48.0%) were discharged within
the first 10 days, including 188 (12.1%) who received steroids. Treated
patients had a median time to death of 1.83 days, compared to 2.79 days for
untreated patients.

As presented in Web Table 1 in Supplementary Materials, we analyzed baseline
characteristics of eligible patients, including demographics, comorbidities,
clinical indices, and medication status. To identify potential predictive
variables, a preliminary subdistribution hazards model was fitted with
treatment and all covariates, revealing age, surgical history, and Glasgow
Coma Scale (GCS) as predictive variables with treatment effect heterogeneity.
As shown in Web Figure 3 in Supplementary Materials, cumulative incidence
plots illustrated this heterogeneity, prompting us to categorize age ($<65$
vs. $\geq 65$) and GCS ($<7$ vs. $7-15$). Conditional on these three
predictive variables, we considered gender, BMI, Sequential Organ Failure
Assessment (SOFA) score, Charlson Comorbidity Index (CCI), and serum lactate
for estimating treatment effects in each subgroup. R code for data analysis is
available on Github as the link provided in Supplementary Materials.

###  5.2 Data analysis results

We estimated CATE using our proposed TMLE framework and compared them with
S-learner. For the TMLE, we reported 95% CI based on the empirical EIF, while
S-learner relied on bootstrapping.

In Figure 3, we present the results of our study regarding treatment effect
heterogeneity, which provides valuable insights into distinct patient
subgroups. Positive estimated CATEs indicate that administering steroids led
to a higher cumulative incidence of death at day 7 compared to the placebo,
and negative ones suggest potential reducing in incidence. As shown in the
sixth and eighth row, patients who did not undergo surgery and had a GCS of 7
or higher (relative higher consciousness level) showed estimated CATEs
significantly greater than 0, indicating that the use of steroids are not
recommended for these patient subgroups. Among these people, those aged 65 or
older showed an estimated CATE of 0.118 (95% CI: 0.018-0.217), while those
younger than 65 years old had an estimated CATE of 0.173 (95% CI:
0.072-0.273). Similarly, for patients younger than 65 who met the criteria for
undergoing surgery and had a GCS below 7, the estimated CATE was 0.144 (95%
CI: 0.010-0.278), also indicating a significant increase in the cumulative
incidence of death at day 7 with steroid treatment compared to the placebo.
The use of steroids is not advised for this patient subgroup either.

Figure 3: Conditional average treatment effect (CATE) defined using the
cumulative incidence function (CIF) at day 7 for patient subgroups categorized
by age ($<$65 or $\geq$65 years old), surgical history (with or without
surgery), and Glasgow Coma Scale (GCS). CATE values were estimated by the
proposed targeted maximum likelihood estimation (TMLE) framework (black
rectangle) and the S-learner (grey circle). The 95% confidence intervals are
provided using the efficient influence function for TMLE, and through
bootstrapping for the S-learner. The dashed line indicates no treatment
effect.

However, for patients aged 65 or older with a GCS below 7 (relatively low
consciousness level) who may have poorer health condition, regardless of
surgery status, the estimated CATEs were negative, with the 95% CI including
zero, suggesting no significant treatment effect. Although there is some
uncertainty, it suggests the possibility of a beneficial effect of steroid
treatment in these patient subgroups. Further investigations are needed to
confirm potential benefits of steroids in this context.

The estimates obtained from our TMLE framework differed from S-learner’s
estimates in various groups. This discrepancy can be attributed to the absence
of the doubly robust feature in the S-learner, and the uncertainty in whether
the outcome, treatment, and censoring models were all correctly specified.

To understand the factors contributing to treatment effect heterogeneity and
identify variables for accurate CATE estimation, we presented Variable
Importance Measures (VIMs) for both predictive and prognostic variables within
each subgroup. In Figure Web Appendix F: Additional Analysis Results of Sepsis
Dataset(A), we illustrate the variable importance of predictive variables,
calculated using the difference in empirical treatment effect variance
($\text{VIM}_{1}$ from equation 13). Notably, GCS emerges as the most
influential predictive variable, explaining the majority of variance in
treatment effects. Omitting GCS from the analysis would reduce treatment
effect variance by 63.7%. Following GCS, age was a critical variable for
detecting heterogeneous treatment effects. We also checked the contribution of
prognostic variables to treatment effect heterogeneity and compared with the
predictive variables, which are negligible as expected, see details in Web
Figure 4 in Supplementary Materials.

Figure 4: Variable importance for predictive (A) and prognostic (B) variables.
In (A), variable importance as estimated by $\text{VIM}_{1}$, represents the
difference in empirical variance of the treatment effect when each variable is
ignored during analysis. In (B), variable importance as estimated by
$\text{VIM}_{2}$, shows the difference in the estimated treatment effect when
each variable is removed from the analysis.

In Figure Web Appendix F: Additional Analysis Results of Sepsis Dataset(B), we
demonstrate the variable importance of prognostic variables, calculated using
the differences in estimated treatment effects ($\text{VIM}_{2}$ in equation
14). For instance, in patients under 65 who underwent surgery and had a GCS
$<7$, omitting serum lactate would result in a 0.0513 difference in treatment
effect estimates, a significant figure given our narrow treatment effect scale
($-1$ to 1). Serum lactate proves vital for estimating treatment effects in
patients with lower consciousness levels (GCS $<7$), while SOFA is more
important for patients with relatively higher consciousness levels (GCS $\geq
7$). This underscores the importance of collecting these variables for
different groups to obtain accurate CATE estimates.

##  6 Discussion

We developed a targeted maximum likelihood estimation (TMLE) framework for
estimating conditional average treatment effects (CATE) in time-to-event data
with competing risks. Instead of using cause-specific hazards, our framework
relies on the cumulative incidence function (CIF), which is based on
subdistribution hazards. This approach facilitates a more straightforward
interpretation of how covariates impact the CIF. Specifically, it allows
clinical researchers to more easily understand and apply the findings to make
informed treatment decisions.

Our proposed estimator exhibits a valuable doubly robust feature, ensuring the
maintenance of consistent results even when the outcome model or both
treatment and censoring mechanisms is misspecified. This feature has been
rigorously demonstrated through both theoretical proofs and extensive
simulations across diverse settings and scenarios. Using our TMLE framework
for treatment decisions enhances strategy reliability, even when modeling
specification are imperfect due to limited knowledge or data, which is the
most cases in real-world problems.

We offer statistical inference with confidence intervals and hypothesis tests
based on asymptotic linearity and normality properties, using estimated
variances from the efficient influence function (EIF) we developed. Our TMLE
framework is flexible, accommodating various regression and machine learning
models for outcomes, treatment, and censoring mechanisms, making it applicable
in diverse scenarios, including the high-dimensional data.

We explored the identification of predictive covariates and proposed two
methods for ranking the importance of predictive and prognostic variables in
accurately estimating CATE.

In our current research, we employed a two-stage approach: first, identifying
predictive variables and subgroups based on these variables, and then
estimating CATE for each subgroup. In our future research, we will focus on
integrating these two stages into one. When the data set is high-dimensional,
LASSO does not significantly help with variable selection and modeling CIF, so
we may need to incorporate more advanced machine learning techniques (e.g.,
deep neural network) to enhance CIF estimation, especially when it is
challenging to correctly specify the outcome and treatment models concurrently
while still maintaining the double robustness feature of CATE estimation.

## Data Availability

Data sharing is not applicable to this article as no new data were created or
analyzed in this paper.

## References

  * Austin and Fine (2017) Austin, P. C. and Fine, J. P. (2017).  Practical recommendations for reporting Fine-Gray model analyses for competing risk data.  Statistics in Medicine, 36, 4391–4400.
  * Bellach et al. (2019) Bellach, A., Kosorok, M. R., Rüschendorf, L., and Fine, J. P. (2019).  Weighted npmle for the subdistribution of a competing risk.  Journal of the American Statistical Association, 114, 259–270.
  * Benkeser et al. (2018) Benkeser, D., Carone, M., and Gilbert, P. B. (2018).  Improved estimation of the cumulative incidence of rare outcomes.  Statistics in Medicine, 37, 280–293.
  * Bickel et al. (1993) Bickel, P. J., Klaassen, C. A., Bickel, P. J., Ritov, Y., Klaassen, J., Wellner, J. A., and Ritov, Y. (1993).  Efficient and Adaptive Estimation for Semiparametric Models, volume 4.  Baltimore, MD: Johns Hopkins University Press.
  * Cai and van der Laan (2020) Cai, W. and van der Laan, M. J. (2020).  One-step targeted maximum likelihood estimation for time-to-event outcomes.  Biometrics, 76, 722–733.
  * Fine and Gray (1999) Fine, J. P. and Gray, R. J. (1999).  A proportional hazards model for the subdistribution of a competing risk.  Journal of the American Statistical Association, 94, 496–509.
  * Fu et al. (2017) Fu, Z., Parikh, C. R., and Zhou, B. (2017).  Penalized variable selection in competing risks regression.  Lifetime Data Analysis, 23, 353–376.
  * Gail et al. (1984) Gail, M. H., Wieand, S., and Piantadosi, S. (1984).  Biased estimates of treatment effect in randomized experiments with nonlinear regressions and omitted covariates.  Biometrika, 71, 431–444.
  * Gill (1994) Gill, R. D. (1994).  Lectures on survival analysis.  In Lectures on Probability Theory, pages 115–241. Springer Berlin Heidelberg, Berlin, Heidelberg.
  * Golmakani and Polley (2020) Golmakani, M. K. and Polley, E. C. (2020).  Super learner for survival data prediction.  The International Journal of Biostatistics, 16, 20190065.
  * Gray (1988) Gray, R. J. (1988).  A class of $K$-sample tests for comparing the cumulative incidence of a competing risk.  The Annals of Statistics, 16, 1141–1154.
  * Hines et al. (2022) Hines, O., Diaz-Ordaz, K., and Vansteelandt, S. (2022).  Variable importance measures for heterogeneous causal effects.  arXiv preprint arXiv:2204.06030.
  * Kawaguchi et al. (2020) Kawaguchi, S., E., Shen, I., J., Li, G., and Suchard, A., M. (2020).  A fast and scalable implementation method for competing risks data with the R package fastcmprsk.  The R Journal, 12, 163–172.
  * Kennedy (2022) Kennedy, E. H. (2022).  Semiparametric doubly robust targeted double machine learning: a review.  arXiv preprint arXiv:2203.06469.
  * Kennedy et al. (2017) Kennedy, E. H., Ma, Z., McHugh, M. D., and Small, D. S. (2017).  Non-parametric methods for doubly robust estimation of continuous treatment effects.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 79, 1229–1245.
  * Koller et al. (2012) Koller, M. T., Raatz, H., Steyerberg, E. W., and Wolbers, M. (2012).  Competing risks and the clinical community: irrelevance or ignorance?  Statistics in Medicine, 31, 1089–1097.
  * Künzel et al. (2019) Künzel, S. R., Sekhon, J. S., Bickel, P. J., and Yu, B. (2019).  Meta-learners for Estimating Heterogeneous Treatment Effects using Machine Learning.  Proceedings of the National Academy of Sciences, 116, 4156–4165.
  * Lau et al. (2009) Lau, B., Cole, S. R., and Gange, S. J. (2009).  Competing risk regression models for epidemiologic data.  American Journal of Epidemiology, 170, 244–256.
  * Levy (2019) Levy, J. (2019).  Tutorial: Deriving the efficient influence curve for large models.  arXiv preprint arXiv:1903.01706.
  * Moore and van der Laan (2009) Moore, K. L. and van der Laan, M. J. (2009).  Application of time-to-event methods in the assessment of safety in clinical trials.  In Peace, K. P., editor, Design and Analysis of Clinical Trials with Time-to-Event Endpoints, pages 473–500. London, UK: Taylor & Francis.
  * Pirracchio et al. (2018) Pirracchio, R., Yue, J. K., Manley, G. T., van Der Laan, M. J., and Hubbard, A. E. (2018).  Collaborative targeted maximum likelihood estimation for variable importance measure: Illustration for functional outcome prediction in mild traumatic brain injuries.  Statistical Methods in Medical Research, 27, 286–297.
  * Robins et al. (1994) Robins, J. M., Rotnitzky, A., and Zhao, L. P. (1994).  Estimation of regression coefficients when some regressors are not always observed.  Journal of the American Statistical Association, 89, 846–866.
  * Rosenbaum and Rubin (1983) Rosenbaum, P. R. and Rubin, D. B. (1983).  The central role of the propensity score in observational studies for causal effects.  Biometrika, 70, 41–55.
  * Rytgaard et al. (2023) Rytgaard, H. C., Eriksson, F., and van der Laan, M. J. (2023).  Estimation of time-specific intervention effects on continuously distributed time-to-event outcomes by targeted maximum likelihood estimation.  Biometrics, 79, 3038–3049.
  * Tsiatis (2006) Tsiatis, A. A. (2006).  Semiparametric Theory and Missing Data.  New York, NY: Springer.
  * Van der Laan et al. (2007) Van der Laan, M. J., Polley, E. C., and Hubbard, A. E. (2007).  Super learner.  Statistical Applications in Genetics and Molecular Biology, 6,.
  * van der Laan et al. (2011) van der Laan, M. J., Rose, S., et al. (2011).  Targeted Learning: Causal Inference for Observational and Experimental Data.  New York, NY: Springer.
  * van der Laan and Rubin (2006) van der Laan, M. J. and Rubin, D. (2006).  Targeted maximum likelihood learning.  The International Journal of Biostatistics, 2, 1–40.
  * van der Vaart (1998) van der Vaart, A. W. (1998).  Asymptotic Statistics.  Cambridge, UK: Cambridge University Press.
  * Venkatesh et al. (2018) Venkatesh, B., Finfer, S., Cohen, J., Rajbhandari, D., Arabi, Y., Bellomo, R., Billot, L., Correa, M., Glass, P., Harward, M., et al. (2018).  Adjunctive glucocorticoid therapy in patients with septic shock.  New England Journal of Medicine, 378, 797–808.

## Supporting Information

Web Appendices, Tables, and Figures, and R code referenced in Sections 3,
3.2.1, 3.2.2, 3.3, 4, and 5 are available with this paper at the Biometrics
website on Wiley Online Library.

Supplementary Materials for Doubly Robust Targeted Estimation of Conditional
Average Treatment Effects for Time-to-event Outcomes with Competing Risks

## Web Appendix A: Algorithm of TMLE for Time-to-Event Outcome with Competing
Risks

Algorithm 1

Preliminary step:

  1. (1)

Fit a subdistribution hazards model that includes all baseline covariates,
treatment, and interaction terms between treatment and covariates.

  2. (2)

Identify predictive covariates based on the fitted regression results and
clinical expertise.

  3. (3)

Split data into $M$ subgroups based the predictive variables.

Within each subgroup $m=1,2,...,M$:

Initial estimation step

  1. (1)

Subdistribution hazard $\hat{\lambda}_{1,m}^{(0)}(t|A,L)$ using unpenalized or
penalized subdistribution hazards model, and CIF as
$\hat{F}_{1,m}^{(0)}\left(t|A,L\right)=1-\exp\\{-\int^{t}_{0}\hat{\lambda}^{(0)}_{1,m}(t|A,L)dt\\}$;

  2. (2)

Treatment assignment $\hat{\pi}_{m}(t|A,L)$ by logistic regression, or super
learner;

  3. (3)

Censoring probability $\hat{G}_{m}(t|A,L)$ by Cox regression, or super
learner.

Targeting step

Update estimated subdistribution hazard and CIF. Let $k=1,2,\ldots$ denote the
iteration index.

Initiate $\hat{\epsilon}=1$. Let $s_{n}$ be a small number, e.g., 0.001.

while $|\hat{\epsilon}|>s_{n}$ do

(1) Obtain $\hat{\epsilon}$ by maximizing the partial likelihood with respect
to the model:
$\lambda_{1,m,\epsilon}^{(k)}\left(t|A,L\right)=\lambda_{1,m}^{(k)}\left(t|A,L\right)\exp\left\\{\epsilon
h_{t_{0},m}^{(k)}\left(t,A,L\right)\right\\},$ where the clever covariate
$h_{m}^{(k)}(t,A,L)=\frac{\\{\mathbb{I}(A=1)-\mathbb{I}(A=0)\\}}{\hat{\pi}_{m}(A|L)}\cdot\frac{\\{1-\hat{F}_{1,m}^{(k)}\left(t_{0}|A,L\right)\\}}{\\{1-\hat{F}_{1,m}^{(k)}(t|A,L)\\}\hat{G}_{m}(t|A,L)}.$

(2) Update estimators:
$\hat{\lambda}_{1,m}^{(k+1)}\left(t|A,L\right)=\hat{\lambda}_{1,m}^{(k)}\left(t|A,L\right)\exp\left\\{\hat{\epsilon}h_{m}^{(k)}\left(t,A,L\right)\right\\}\text{,
and}$
$\hat{F}_{1,m}^{(k+1)}\left(t|A,L\right)=1-\exp\left\\{-\int_{0}^{t}\hat{\lambda}_{1,m}^{(k+1)}\left(t|A,L\right)dt\right\\}.$
$k=k+1$

end while

CATE estimation:

Denote the final estimators obtained from Step 2 as
$\hat{F}_{1,m}^{*}(t|A,L)$. Plug in the estimates and calculate the average
treatment effect within each subgroup $m$ to obtain the CATE:
$\Psi_{m}(\hat{P}_{n}^{*})=\frac{1}{n_{m}}\sum_{i=1}^{n_{m}}\left\\{(\hat{F}_{1,m}^{*}(t_{0}|L_{i},A_{i}=1)-\hat{F}_{1,m}^{*}(t_{0}|L_{i},A_{i}=0)\right\\}.$

## Web Appendix B: Derivation of the Efficient Influence Function

The EIF is derived by calculating the pathwise derivative
$\frac{d}{d\epsilon}\Psi(P_{\epsilon})\Bigr{|}_{\epsilon=0}$ for a submodel
$P_{\epsilon}$: $dP_{\epsilon}=\\{1+\epsilon b\\}dP$. Our goal is to express
the pathwise derivative of the target estimand as an inner product of two
terms, i.e.,
$\frac{d}{d\epsilon}\Psi(P_{\epsilon})\Bigr{|}_{\epsilon=0}=E\\{D(P)\cdot
b\\}$.

Under Assumptions 1-4, our target estimand, CATE, can be identified and
expressed as

$\displaystyle\Psi(P)$
$\displaystyle=\int_{\mathcal{L}}\sum_{a=0,1}F_{1}\left(t_{0}|A=a,L=l\right)\left\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\right\\}\mu(l)dl$
$\displaystyle=\int_{\mathcal{L}}\sum_{a=0,1}\left\\{1-\exp\left(-\int_{0}^{t_{0}}\Lambda_{1}(dt|a,l)\right)\right\\}\left\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\right\\}\mu(l)dl,$

where $\mu(l)$ is the density of prognostic variables $L$, and

$\displaystyle\quad\quad\quad\Lambda_{1}(dt|a,l)$
$\displaystyle=\frac{\operatorname{Pr}\\{t\leq{T}<t+dt,{\Delta}=1,C\geq
t|a,l\\}}{\operatorname{Pr}[\\{{T}\geq t\cup(T\leq t\cap{\Delta}\neq 1)\\}\cap
C\geq t|a,l]}$
$\displaystyle=\frac{\operatorname{Pr}\\{t\leq{T}<t+dt,{\Delta}=1|a,l\\}\operatorname{Pr}\\{C\geq
t|a,l\\}}{\operatorname{Pr}\\{{T}\geq t\cup({T}\leq t\cap\Delta\neq
1)|a,l\\}\operatorname{Pr}\\{C\geq t|a,l\\}}$
$\displaystyle=\frac{P^{1}(dt|a,l)}{\sum_{\delta}\\{\int_{0}^{\infty}P^{\delta}(ds|a,l)-\int_{0}^{t}\mathbbm{1}(\delta=1)P^{\delta}(ds|a,l)\\}},$

where
$P^{\delta}(dt|a,l)=\operatorname{Pr}(t\leq{T}<t+dt,{\Delta}=\delta|a,l),\delta\in\\{1,2\\}$,
assuming two causes for simplicity.

The pathwise derivative of the target estimand

$\displaystyle\frac{d}{d\epsilon}\Psi(P_{\epsilon})\Bigr{|}_{\epsilon=0}$
$\displaystyle=\frac{d}{d\epsilon}\int_{\mathcal{L}}\sum_{a=0,1}\left\\{1-\exp\left(-\int_{0}^{t_{0}}\Lambda_{1,\epsilon}(dt|a,l)\right)\right\\}\left\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\right\\}\mu_{\epsilon}(l)dl\Bigr{|}_{\epsilon=0}$
$\displaystyle=\frac{d}{d\epsilon}\int_{\mathcal{L}}\sum_{a=0,1}\left[1-\exp\left\\{-\int_{0}^{t_{0}}\Lambda_{1,\epsilon}(dt|a,l)\right\\}\right]\left\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\right\\}\mu_{\epsilon}(l)dl\Bigr{|}_{\epsilon=0}$
$\displaystyle=\int_{\mathcal{L}}\sum_{a=0,1}\exp\left\\{-\int_{0}^{t_{0}}\Lambda_{1}(dt|a,l)\right\\}\left\\{\int_{0}^{t_{0}}\frac{d}{d\epsilon}\Lambda_{1,\epsilon}(dt|a,l)\Bigr{|}_{\epsilon=0}\right\\}\left\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\right\\}\mu(l)dl$
$\displaystyle\quad\quad+\int_{\mathcal{L}}\sum_{a=0,1}\left[1-\exp\left\\{-\int_{0}^{t_{0}}\Lambda_{1}(dt|a,l)\right\\}\right]\left\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\right\\}\frac{d}{d\epsilon}\mu_{\epsilon=0}\Bigr{|}_{\epsilon=0}dl$
$\displaystyle:=(*),$

where

$\displaystyle\frac{d}{d\epsilon}\Lambda_{1,\epsilon}(dt|a,l)\Bigr{|}_{\epsilon=0}=\frac{d}{d\epsilon}\frac{P_{\epsilon}^{1}(dt|a,l)}{P_{\epsilon}\\{{T}\geq
t\cup({T}\leq t\cap\Delta\neq 1)|a,l\\}}\Bigr{|}_{\epsilon=0}$
$\displaystyle=P\\{{T}\geq t\cup({T}\leq t\cap\Delta\neq 1)|a,l\\}^{-2}\cdot$
$\displaystyle\quad\quad\left[P^{1}(dt|a,l)b(t,1,a,l)P\\{{T}\geq t\cup({T}\leq
t\cap\Delta\neq 1)|a,l\\}\right.$
$\displaystyle\quad\quad-P^{1}(dt|a,l)\sum_{\delta=1,2}\int_{0}^{\infty}\\{P^{\delta}(ds|a,l)-\mathbbm{1}(T<t,\delta=1)\left.P^{\delta}(ds|a,l)\\}b(s,\delta,a,l)\right]$
$\displaystyle=\sum_{\delta=1,2}\left[\frac{\mathbbm{1}(\delta=1)P^{\delta}(dt|a,l)b(t,\delta,a,l)}{P\\{{T}\geq
t\cup({T}\leq t\cap\Delta\neq 1)|a,l\\}}\right.$
$\displaystyle\quad\quad-\left.\frac{\Lambda_{1}(dt|a,l)}{P\left\\{{T}\geq
t\cup({T}\leq t\cap\Delta\neq
1)|a,l\right\\}}\int_{0}^{\infty}Y(t)P^{\delta}(ds|a,l)b(s,\delta,a,l)\right]\text{,
and}$

by the key identity (1) in Levy (2019),

$\frac{d}{d\epsilon}\mu_{\epsilon}(l)\Bigr{|}_{\epsilon=0}=\left[E_{P}\left\\{b(T,\Delta,A,L)|l\right\\}-E_{P}\left\\{b(T,\Delta,A,L)\right\\}\right]\mu(l).$

Let (1)
$\displaystyle=\int_{0}^{t_{0}}\frac{d}{d\epsilon}\Lambda_{1,\epsilon}(dt|a,l)\Bigr{|}_{\epsilon=0}$
$\displaystyle=\int_{0}^{\infty}\sum_{\delta=1,2}\mathbbm{1}(t\leq
t_{0})\left[\frac{\mathbbm{1}(\delta=1)P^{\delta}(dt|a,l)b(t,\delta,a,l)}{P\left\\{{T}\geq
t\cup({T}\leq t\cap\Delta\neq 1)|a,l\right\\}}\right.$
$\displaystyle\quad\quad-\left.\frac{\Lambda_{1}(dt|a,l)}{P\\{{T}\geq
t\cup({T}\leq t\cap\Delta\neq
1)|a,l\\}}\int_{0}^{\infty}Y(t)P^{\delta}(ds|a,l)b(s,\delta,a,l)\right]$
$\displaystyle=\int_{0}^{\infty}\sum_{\delta=1,2}\left[\frac{\mathbbm{1}(t\leq
t_{0},\delta=1)}{P\\{{T}\geq t\cup({T}\leq t\cap\Delta\neq
1)|a,l\\}}P^{\delta}(dt|a,l)b(t,\delta,a,l)\right.$
$\displaystyle\quad\quad-\int_{0}^{\infty}\left.\frac{\Lambda_{1}(dt|a,l)\mathbbm{1}(t\leq
t_{0})}{P\\{{T}\geq t\cup({T}\leq t\cap\Delta\neq
1)|a,l\\}}Y(t)P^{\delta}(ds|a,l)b(s,\delta,a,l)\right]$
$\displaystyle=\int_{0}^{\infty}\sum_{\delta=1,2}\left[\int^{t_{0}}_{0}\frac{N(dt)}{P\\{{T}\geq
t\cup({T}\leq t\cap\Delta\neq
1)|a,l\\}}P^{\delta}(ds|a,l)b(s,\delta,a,l)\right.$
$\displaystyle\quad\quad-\left.\int_{0}^{t_{0}}\frac{\Lambda_{1}(dt|a,l)Y(t)}{P\\{{T}\geq
t\cup({T}\leq t\cap\Delta\neq
1)|a,l\\}}P^{\delta}(ds|a,l)b(s,\delta,a,l)\right]$
$\displaystyle=\int_{0}^{\infty}\sum_{\delta=1,2}\left\\{\int^{t_{0}}_{0}\frac{N(dt)-\Lambda_{1}(dt|a,l)Y(t)}{1-F_{1}(t|a,l)}\right\\}P^{\delta}(ds|a,l)b(s,\delta,a,l);\text{
and}$

(2)
$\displaystyle=\int_{\mathcal{L}}\sum_{a=0,1}\left[1-\exp\left\\{-\int_{0}^{t_{0}}\Lambda_{1}(dt|a,l)\right\\}\right]\left\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\right\\}\frac{d}{d\epsilon}\mu_{\epsilon=0}\Bigr{|}_{\epsilon=0}\mu(l)dl$
$\displaystyle=\int_{\mathcal{L}}\sum_{a=0,1}\left[1-\exp\left\\{-\int_{0}^{t_{0}}\Lambda_{1}(dt|a,l)\right\\}\right]\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\\}\cdot$
$\displaystyle\quad\quad\left[E_{P}\left\\{b(T,\Delta,A,L)|l\right\\}-E_{P}\left\\{b(T,\Delta,A,L)\right\\}\right]\mu(l)dl$
$\displaystyle=\int_{\mathcal{L}}\left[\left\\{1-\exp\left(\int_{0}^{t_{0}}\Lambda_{1}(dt|1,l)\right)\right\\}\right.$
$\displaystyle\quad\quad-\left.\left\\{1-\exp\left(\int_{0}^{t_{0}}\Lambda_{1}(dt|0,l)\right)\right\\}\right]E_{P}\\{b(T,\Delta,A,L)|l\\}\mu(l)dl$
$\displaystyle\quad\quad-\int_{\mathcal{L}}\left[\left\\{1-\exp\left(\int_{0}^{t_{0}}\Lambda_{1}(dt|1,l)\right)\right\\}\right.$
$\displaystyle\quad\quad-\left.\left\\{1-\exp\left(\int_{0}^{t_{0}}\Lambda_{1}(dt|0,l)\right)\right\\}\right]E_{P}\\{b(T,\Delta,A,L)\\}\mu(l)dl$
$\displaystyle=E_{P}\left[\\{F_{1}(t_{0}|1,L)-F_{1}(t_{0}|0,L)\\}b({T},\Delta,A,L)\right]-E_{P}\\{\Psi(P)b({T},\Delta,A,L)\\}.$

Plug (1) and (2) back into (*), we have

$\displaystyle\frac{d}{d\epsilon}\Psi(P_{\epsilon})\Bigr{|}_{\epsilon=0}=$
$\displaystyle\int_{\mathcal{L}}\sum_{a=0,1}\exp\left\\{-\int_{0}^{t_{0}}\Lambda_{1}(dt|a,l)\right\\}\left\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\right\\}$
$\displaystyle\quad\quad\left[\int_{0}^{\infty}\sum_{\delta=1,2}\left\\{\int^{t_{0}}_{0}\frac{N(dt)-\Lambda_{1}(dt|a,l)Y(t)}{1-F_{1}(t|a,l)}\right\\}P^{\delta}(ds|a,l)b(s,\delta,a,l)\right]\mu(l)dl$
$\displaystyle\quad\quad+E_{P}\left[\\{F_{1}(t_{0}|1,L)-F_{1}(t_{0}|0,L)\\}b({T},\Delta,A,L)\right]-E_{P}\\{\Psi(P)b({T},\Delta,A,L)\\}$
$\displaystyle=$
$\displaystyle\int_{\mathcal{L}}\sum_{a=0,1}\left\\{1-F_{1}(t_{0}|a,l)\right\\}\left\\{\frac{\mathbbm{1}(a=1)}{\pi(a|l)}-\frac{\mathbbm{1}(a=0)}{\pi(a|l)}\right\\}\pi(a|l)$
$\displaystyle\quad\quad\left[\int_{0}^{\infty}\sum_{\delta=1,2}\left\\{\int^{t_{0}}_{0}\frac{N(dt)-\Lambda_{1}(dt|a,l)Y(t)}{1-F_{1}(t|a,l)}\right\\}P^{\delta}(ds|a,l)b(s,\delta,a,l)\right]\mu(l)dl$
$\displaystyle\quad\quad+E_{P}\left[\\{F_{1}(t_{0}|1,L)-F_{1}(t_{0}|0,L)\\}b({T},\Delta,A,L)\right]-E_{P}\\{\Psi(P)b({T},\Delta,A,L)\\}$
$\displaystyle=$
$\displaystyle\int_{\mathcal{L}}\sum_{a=0,1}\pi(a|l)\int_{0}^{\infty}\sum_{\delta=1,2}\left\\{\int^{t_{0}}_{0}\left[\frac{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)}{\pi(a|l)}\cdot\frac{1-F_{1}(t_{0}|a,l)}{1-F_{1}(t|a,l)}\cdot\right.\right.$
$\displaystyle\quad\quad\left.\left.\left\\{N(dt)-\Lambda_{1}(dt|a,l)Y(t)\right\\}\right]P^{\delta}(ds|a,l)b(s,\delta,a,l)\right\\}\mu(l)dl$
$\displaystyle\quad\quad+E_{P}\left[\\{F_{1}(t_{0}|1,L)-F_{1}(t_{0}|0,L)\\}b({T},\Delta,A,L)\right]-E_{P}\\{\Psi(P)b({T},\Delta,A,L)\\}$

$\displaystyle=$ $\displaystyle
E_{P}\left(\int^{t_{0}}_{0}\left[\frac{\mathbbm{1}(A=1)-\mathbbm{1}(A=0)}{\pi(A|L)}\cdot\frac{1-F_{1}(t_{0}|A,L)}{1-F_{1}(t|A,L)}\left\\{N(dt)-\Lambda_{1}(dt|A,L)Y(t)\right\\}\right]b(T,\Delta,A,L)\right)$
$\displaystyle\quad\quad+E_{P}\left[\left\\{F_{1}(t_{0}|1,L)-F_{1}(t_{0}|0,L)\right\\}b({T},\Delta,A,L)\right]-E_{P}\left\\{\Psi(P)b({T},\Delta,A,L)\right\\}$
$\displaystyle=$ $\displaystyle
E_{P}\left\\{\left[\int^{t_{0}}_{0}\frac{\mathbbm{1}(A=1)-\mathbbm{1}(A=0)}{\pi(A|L)}\cdot\frac{1-F_{1}(t_{0}|A,L)}{1-F_{1}(t|A,L)}\left\\{N(dt)-\Lambda_{1}(dt|A,L)Y(t)\right\\}\right.\right.$
$\displaystyle\quad\quad+\left.\left.\\{F_{1}(t_{0}|1,L)-F_{1}(t_{0}|0,L)\\}-\Psi(P)\right]b({T},\Delta,A,L)\right\\}.$

Now we have shown that the pathwise derivative of the targeted estimand can be
expressed as an inner product of $D(P)$ and $b({T},\Delta,A,L)$. Therefore,
the EIF of $\Psi(P)$ for complete data (no censoring) is

$\displaystyle D(P)$
$\displaystyle=\int^{t_{0}}_{0}\frac{\mathbbm{1}(A=1)-\mathbbm{1}(A=0)}{\pi(A|L)}\cdot\frac{1-F_{1}(t_{0}|A,L)}{1-F_{1}(t|A,L)}\left\\{N(dt)-\Lambda_{1}(dt|A,L)Y(t)\right\\}$
$\displaystyle\quad\quad+\\{F_{1}(t_{0}|1,L)-F_{1}(t_{0}|0,L)\\}-\Psi(P),$

with the clever covariate is

$h(t,A,L)=\frac{\mathbbm{1}(A=1)-\mathbbm{1}(A=0)}{\pi(A|L)}\cdot\frac{1-F_{1}(t_{0}|A,L)}{1-F_{1}(t|A,L)}.$

For right-censored data, $N(t)$ and $\left\\{1-Y(t)\right\\}$ are not
observable when $r(t)=\mathbbm{1}(C\geq T\wedge t)=0$, but $r(t)N(t)$ and
$r(t)\left\\{1-Y(t)\right\\}$ are computable for $r(t)=0,1$ (Fine and Gray,
1999). The term $\frac{r(t)}{G\\{(\tilde{T}\wedge t)^{-}|A,L\\}}$ has an
expectation 1 conditional of $T,\Delta,$ and $L$. We multiply the first part
of EIF by $\frac{w(t)}{G(t^{-}|A,L)}$, then the EIF of $\Psi(P)$ for the
right-censored data has the form:

$\displaystyle D(P)$
$\displaystyle=\int^{t_{0}}_{0}\frac{\mathbbm{1}(A=1)-\mathbbm{1}(A=0)}{\pi(A|L)}\cdot\frac{1-F_{1}(t_{0}|A,L)}{1-F_{1}(t|A,L)}\cdot\frac{w(t)}{G(t^{-}|A,L)}\left\\{N(dt)-\Lambda_{1}(dt|A,L)Y(t)\right\\}$
$\displaystyle\quad\quad+\\{F_{1}(t_{0}|1,L)-F_{1}(t_{0}|0,L)\\}-\Psi(P),$

where $w(t)=\frac{\mathbbm{1}(C\geq T\wedge
t)G(t^{-}|A,L)}{G\\{(\tilde{T}\wedge t)^{-}|A,L\\}}$ and the clever covariate
being

$h(t,A,L)=\frac{\mathbbm{1}(A=1)-\mathbbm{1}(A=0)}{\pi(A|L)G(t^{-}|A,L)}\cdot\frac{1-F_{1}(t_{0}|A,L)}{1-F_{1}(t|A,L)}.$

## Web Appendix C: Derivation of the Score Function of the Fluctuation Model

The fluctuation model that we use for updating the estimated target parameter
is

$\lambda_{1,\epsilon}(t|A,L)=\lambda_{1}(t|A,L)\exp\\{\epsilon h(t,A,L)\\}.$

The weighted likelihood (the pseudo likelihood (page 5, section 2.3) in
Bellach et al. (2019)) has the form:

$\mathscr{L}_{n}(\lambda_{1,\epsilon})=\prod_{i=1}^{n}\left[\left\\{w_{i}(t)Y_{i}(t)\lambda_{1,\epsilon}(t|A_{i},L_{i})^{\mathbbm{1}(C_{i}\geq
t)N_{i}(dt)}\right\\}\exp\left\\{-\int_{0}^{t_{0}}w_{i}(t)Y_{i}(t)\lambda_{1,\epsilon}(t|A_{i},L_{i})dt\right\\}\right];$

and the weighted log-likelihood function (page 8, section 3.2) in Bellach et
al. (2019))

$\displaystyle\log\mathscr{L}_{n}(\lambda_{1,\epsilon})$
$\displaystyle=\sum_{i=1}^{n}\left[\int_{0}^{t_{0}}\mathbbm{1}(C_{i}\geq
t)Y_{i}(t)w_{i}(t)\log\left\\{\lambda_{1,\epsilon}(t|A_{i},L_{i})\right\\}N_{i}(dt)-\int_{0}^{t_{0}}w_{i}(t)Y_{i}(t)\lambda_{1,\epsilon}\left(t|A_{i},L_{i}\right)dt\right].$

Since $N_{i}(dt)=N_{i}(t)-N_{i}(t-)=0$ as long as $\mathbbm{1}(C_{i}\geq
t)Y_{i}(t)=0$, we can ignore $\mathbbm{1}(C_{i}\geq t)Y_{i}(t)$ in the first
part of the likelihood. Therefore, the score function has the form:

$\displaystyle\frac{d}{d\epsilon}\log\mathscr{L}_{n}\left(\lambda_{1,\epsilon}\right)$
$\displaystyle=\sum_{i=1}^{n}\left[\int_{0}^{t_{0}}w_{i}(t)\frac{d}{d\epsilon}\log\left\\{\lambda_{1,\epsilon}(t|A_{i},L_{i})\right\\}N_{i}(dt)-\int_{0}^{t_{0}}w_{i}(t)Y_{i}(t)\frac{d}{d\epsilon}\lambda_{1,\epsilon}(t|A_{i},L_{i})dt\right]$
$\displaystyle=\sum_{i=1}^{n}\left[\int_{0}^{t_{0}}w_{i}(t)\frac{\lambda_{1}(t|A_{i},L_{i})\exp\\{\epsilon
h(t,A_{i},L_{i})\\}}{\lambda_{1}(t|A_{i},L_{i})\exp\\{\epsilon
h(t,A_{i},L_{i})\\}}h(t,A_{i},L_{i})N_{i}(dt)\right.$
$\displaystyle\quad\quad-\left.\int_{0}^{t_{0}}w_{i}(t)Y_{i}(t)\lambda_{1}(t|A_{i},L_{i})\exp\\{\epsilon
h(t,A_{i},L_{i})\\}h(t,A_{i},L_{i})dt\right]$
$\displaystyle=\sum_{i=1}^{n}\int_{0}^{t_{0}}w_{i}(t)h(t,A_{i},L_{i})\left[N_{i}(dt)-Y_{i}(t)\lambda_{1}(t|A_{i},L_{i})\exp\\{\epsilon
h(t,A_{i},L_{i})\\}dt\right].$

Then,

$\displaystyle\frac{1}{n}\cdot\frac{d}{d\epsilon}\log\mathscr{L}_{n}$
$\displaystyle\left(\lambda_{1,\epsilon}\right)\Bigr{|}_{\epsilon=0}=\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{t_{0}}w_{i}(t)h(t,A_{i},L_{i})\left[N_{i}(dt)-Y_{i}(t)\lambda_{1}(t|A_{i},L_{i})dt\right]$

is exactly the same as the empirical average of the first part of EIF.
Empirically, the maximum likelihood estimator $\hat{\epsilon}$ is obtained by
using R package nleqslv.

## Web Appendix D: Asymptotic Properties: Derivation of the Second-order
Remainder

By the von Mises expansion (distributional Taylor expansion),

$\displaystyle\Psi(\hat{P}^{*}_{n})-\Psi(P)$
$\displaystyle=-E(D(\hat{P}^{*}_{n}))+R(\hat{P}^{*}_{n},P)$ (15)
$\displaystyle=\left[\mathbb{P}_{n}\left\\{D(\hat{P}^{*}_{n})-D(P)\right\\}-E\left\\{D(\hat{P}^{*}_{n})-D(P)\right\\}\right]$
$\displaystyle\quad-\mathbb{P}_{n}(D(\hat{P}^{*}_{n}))+\mathbb{P}_{n}(D(P))$
$\displaystyle\quad+R(\hat{P}^{*}_{n},P),$

The first term is $o_{p}(n^{-1/2})$ if $D(\hat{P}^{*}_{n})$ belongs to the
Donsker class by the empirical process theory (van der Vaart, 1998; van der
Laan et al., 2011). The second term, the mean of EIF, is $o_{p}(n^{-1/2})$ by
solving the efficient influence function equation in the targeting step (see
(5)). The third term is $o_{p}(n^{-1/2})$ by the Central Limit Theorem. To
show the asymptotic properties of the targeted estimator
$\Psi(\hat{P}_{n}^{*})$, we need to show that the last term, the second order
remainder $R(\hat{P}_{n}^{*},P_{0})=o_{P}(n^{-1/2})$.

The second-order remainder
$R(\hat{P}^{*}_{n},P)=\Psi(\hat{P}^{*}_{n})-\Psi(P)+E\left\\{D(\hat{P}^{*}_{n})\right\\},$
where

$\displaystyle E\left\\{D(\hat{P}^{*}_{n})\right\\}$
$\displaystyle=E\left[\int_{0}^{t_{0}}\hat{h}(t,A,L)\hat{w}(t)\left\\{N(dt)-\hat{\Lambda}_{1}^{*}(dt|A,L)Y(t)\right\\}+\left\\{\hat{F}_{1}^{*}(t_{0}|1,L)-\hat{F}^{*}_{1}(t_{0}|0,L)\right\\}-\Psi(\hat{P})\right]$
$\displaystyle=E\left[\int_{0}^{t_{0}}\frac{\\{\mathbbm{1}(A=1)-\mathbbm{1}(A=0)\\}}{\hat{\pi}(A|L)}\cdot\frac{\\{1-\hat{F}^{*}_{1}\left(t_{0}|A,L\right)\\}}{\\{1-\hat{F}^{*}_{1}(t|A,L)\\}\hat{G}(t^{-}|A,L)}\right.\cdot$
$\displaystyle\quad\quad\left.\frac{\mathbbm{1}(C\geq T\wedge
t)\hat{G}(t^{-}|A,L)}{\hat{G}\\{(\tilde{T}\wedge
t)^{-}|A,L\\}}\left\\{N(dt)-\hat{\Lambda}_{1}^{*}(dt|A,L)Y(t)\right\\}\right]$
$\displaystyle\quad\quad+E\left\\{\hat{F}_{1}^{*}(t_{0}|1,L)-\hat{F}_{1}^{*}(t_{0}|0,L)\right\\}-E\left\\{\Psi(\hat{P}^{*}_{n})\right\\}$
$\displaystyle=\int_{\mathcal{L}}\Sigma_{a=0,1}\int_{0}^{t_{0}}\frac{\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\\}\pi(a|l)}{\hat{\pi}(a|l)}\cdot$
$\displaystyle\quad\quad\frac{\\{1-\hat{F}_{1}^{*}\left(t_{0}|a,l\right)\\}\\{1-F_{1}(t^{-}|a,l)\\}G\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}{\\{1-\hat{F}_{1}^{*}(t|a,l)\\}\hat{G}\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}\left\\{\Lambda_{1}(dt|a,l)-\hat{\Lambda}_{1}^{*}(dt|a,l)\right\\}\mu(l)dl$
$\displaystyle\quad\quad+E\left\\{\hat{F}_{1}^{*}(t_{0}|1,L)-\hat{F}_{1}^{*}(t_{0}|0,L)\right\\}-E\left\\{\Psi(\hat{P}^{*}_{n})\right\\}$
$\displaystyle\quad\quad({\color[rgb]{.5,.5,.5}\text{Because
}E\left\\{f(Y)\mathbbm{1}(A=1)\right\\}=E(f(Y)|A=1)P(A=1),}$ where $f(Y)$ is a
function of $Y$, and $A$ is a random variable
$\displaystyle=\int_{\mathcal{L}}\Sigma_{a=0,1}\int_{0}^{t_{0}}\left[\frac{\pi(a|l)}{\hat{\pi}(a|l)}\frac{G\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}{\hat{G}\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}\frac{\left\\{1-\hat{F}_{1}^{*}\left(t_{0}|a,l\right)\right\\}}{\left\\{1-\hat{F}_{1}^{*}(t|a,l)\right\\}}\left\\{\Lambda_{1}(dt|a,l)-\hat{\Lambda}_{1}^{*}(dt|a,l)\right\\}\right.$
$\displaystyle\quad\quad\quad\left.\left\\{1-F_{1}(t^{-}|a,l)\right\\}\right]\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\\}\mu(l)dl$
$\displaystyle\quad\quad+E\left\\{\hat{F}_{1}^{*}(t_{0}|1,L)-\hat{F}_{1}^{*}(t_{0}|0,L)\right\\}-E\left\\{\Psi(\hat{P}^{*}_{n})\right\\}.$

Therefore, the second-order remainder

$\displaystyle R(\hat{P}^{*}_{n},P)$
$\displaystyle=\Psi(\hat{P}^{*}_{n})-\Psi(P)+E\left\\{D(P^{*}_{n})\right\\}$
$\displaystyle=\Psi(\hat{P}^{*}_{n})-\Psi(P)$
$\displaystyle\quad+\int_{\mathcal{L}}\Sigma_{a=0,1}\int_{0}^{t_{0}}\left[\frac{\pi(a|l)}{\hat{\pi}(a|l)}\frac{G\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}{\hat{G}\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}\frac{\left\\{1-\hat{F}_{1}^{*}\left(t_{0}|a,l\right)\right\\}}{\left\\{1-\hat{F}_{1}^{*}(t|a,l)\right\\}}\left\\{\Lambda_{1}(dt|a,l)-\hat{\Lambda}_{1}^{*}(dt|a,l)\right\\}\right.$
$\displaystyle\quad\quad\quad\left.\left\\{1-F_{1}(t^{-}|a,l)\right\\}\right]\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\\}\mu(l)dl$
$\displaystyle\quad\quad+E\left\\{\hat{F}_{1}^{*}(t_{0}|1,L)-\hat{F}_{1}^{*}(t_{0}|0,L)\right\\}-E\left\\{\Psi(\hat{P}^{*}_{n})\right\\}$
$\displaystyle=\int_{\mathcal{L}}\Sigma_{a=0,1}\int_{0}^{t_{0}}\left[\frac{\pi(a|l)}{\hat{\pi}(a|l)}\frac{G\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}{\hat{G}\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}\frac{\left\\{1-\hat{F}_{1}^{*}\left(t_{0}|a,l\right)\right\\}}{\left\\{1-\hat{F}_{1}^{*}(t|a,l)\right\\}}\left\\{\Lambda_{1}(dt|a,l)-\hat{\Lambda}_{1}^{*}(dt|a,l)\right\\}\right.$
$\displaystyle\quad\quad\quad\left.\left\\{1-F_{1}(t^{-}|a,l)\right\\}\right]\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\\}\mu(l)dl$
$\displaystyle\quad\quad+E\left\\{\hat{F}_{1}^{*}(t_{0}|1,L)-\hat{F}_{1}^{*}(t_{0}|0,L)\right\\}-\Psi(P)$
$\displaystyle=\int_{\mathcal{L}}\Sigma_{a=0,1}\int_{0}^{t_{0}}\left[\frac{\pi(a|l)}{\hat{\pi}(a|l)}\frac{G\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}{\hat{G}\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}\frac{\left\\{1-\hat{F}_{1}^{*}\left(t_{0}|a,l\right)\right\\}}{\left\\{1-\hat{F}_{1}^{*}(t|a,l)\right\\}}\left\\{\Lambda_{1}(dt|a,l)-\hat{\Lambda}_{1}^{*}(dt|a,l)\right\\}\right.$
$\displaystyle\quad\quad\quad\left.\left\\{1-F_{1}(t^{-}|a,l)\right\\}\right]\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\\}\mu(l)dl$
$\displaystyle\quad\quad\quad+\int_{\mathcal{L}}\Sigma_{a=0,1}\left\\{\hat{F}_{1}^{*}(t_{0}|a,l)-F_{1}(t_{0}|a,l)\right\\}\left\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\right\\}\mu(l)dl$
$\displaystyle=\int_{\mathcal{L}}\Sigma_{a=0,1}\int_{0}^{t_{0}}\left[\frac{\pi(a|l)}{\hat{\pi}(a|l)}\frac{G\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}{\hat{G}\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}\frac{\left\\{1-\hat{F}_{1}^{*}\left(t_{0}|a,l\right)\right\\}}{\left\\{1-\hat{F}_{1}^{*}(t|a,l)\right\\}}\left\\{\Lambda_{1}(dt|a,l)-\hat{\Lambda}_{1}^{*}(dt|a,l)\right\\}\right.$
$\displaystyle\quad\quad\quad\left.\left\\{1-F_{1}(t^{-}|a,l)\right\\}\right]\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\\}\mu(l)dl$
$\displaystyle\quad\quad\quad+\int_{\mathcal{L}}\Sigma_{a=0,1}\left\\{1-F_{1}(t^{-}|a,l)\right\\}\left\\{\Lambda_{1}(dt|a,l)-\hat{\Lambda}_{1}^{*}(dt|a,l)\right\\}\frac{\left\\{1-\hat{F}_{1}^{*}(t_{0}|a,l)\right\\}}{\left\\{1-\hat{F}_{1}^{*}(t|a,l)\right\\}}$
$\displaystyle\quad\quad\quad\left\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\right\\}\mu(l)dl{\color[rgb]{.5,.5,.5}\text{
\quad\quad by the Duhamel equation
\cite[citep]{(\@@bibref{AuthorsPhrase1Year}{gill_lectures_1994}{\@@citephrase{,
}}{})},}}$
$\displaystyle{\color[rgb]{.5,.5,.5}\quad\hat{F}_{1}^{*}(t_{0}|a,L)-F_{1}(t_{0}|a,l)=-\int_{0}^{t_{0}}\left\\{1-F_{1}(t^{-}|a,l)\right\\}\left\\{\Lambda_{1}(dt|a,l)-\hat{\Lambda}_{1}^{*}(dt|a,l)\right\\}\frac{\left\\{1-\hat{F}_{1}^{*}\left(t_{0}|a,l\right)\right\\}}{\left\\{1-\hat{F}_{1}^{*}(t|a,l)\right\\}}\text{}}$
$\displaystyle=\int_{\mathcal{L}}\Sigma_{a=0,1}\int_{0}^{t_{0}}\left(\left[\frac{\pi(a|l)}{\hat{\pi}(a|l)}\frac{G\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}{\hat{G}\\{(\tilde{T}\wedge
t)^{-}|a,l\\}}-1\right]\left\\{\Lambda_{1}(dt|a,l)-\hat{\Lambda}_{1}^{*}(dt|a,l)\right\\}\right.$
$\displaystyle\quad\quad\quad\left.\frac{\left\\{1-\hat{F}_{1}^{*}\left(t_{0}|a,l\right)\right\\}}{\left\\{1-\hat{F}_{1}^{*}(t|a,l)\right\\}}\left\\{1-F_{1}(t^{-}|a,l)\right\\}\right)\\{\mathbbm{1}(a=1)-\mathbbm{1}(a=0)\\}\mu(l)dl$
$\displaystyle=O_{p}\left(\Sigma_{a=0,1}\left\lVert\int_{0}^{t_{0}}\left[\pi(a|l)G\\{(\tilde{T}\wedge
t)^{-}|a,l\\}-\hat{\pi}(a|l)\hat{G}\\{(\tilde{T}\wedge
t)^{-}|a,l\\}\right]dt\right\rVert\right.\cdot$
$\displaystyle\quad\quad\quad\left.\left\lVert\int_{0}^{t_{0}}\left\\{\Lambda_{1}(dt|a,l)-\hat{\Lambda}_{1}^{*}(dt|a,l)\right\\}\right\rVert\right),{\color[rgb]{.5,.5,.5}\text{by
the Cauchy-Schwarz inequality.}}$

The remainder is equal to $o_{p}(n^{-1/2})$ when either both the treatment and
censoring mechanisms are consistently estimated or the CIF is consistently
estimated, i.e., either
$\left\lVert\\{\pi(a|l)G(t|a,l)\\}-\\{\hat{\pi}(a|l){\hat{G}(t|a,l)}\\}\right\rVert=o_{p}(n^{-1/2})$,
or
$\left\lVert\Lambda_{1}(dt|a,l)-\hat{\Lambda}_{1}^{*}(dt|a,l)\right\rVert=o_{p}(n^{-1/2})$.
$R(\hat{P}_{n}^{*},P)=o_{p}(n^{-1/2})$ can also be achieved if both
$\hat{\pi}_{0}(a|l)\hat{G}(t|a,l)$ and $\hat{\lambda}_{1}^{*}(t|a,l)$ converge
at a rate faster than $o_{P}(n^{-1/4})$.

## Web Appendix E: Additional Simulation Results

Web Figure 1.  Bias (upper panels) and RMSE (lower panels) of CATE estimations
at the median follow-up time point under different model specifications. Each
panel corresponds to one of the four subgroups determined by two binary
predictive variables, $V_{1}$ and $V_{2}$. The red line represents zero bias.
$\checkmark$ and $\times$ indicate that the corresponding models were
correctly specified and misspecified, respectively. Light grey and dark grey
represent S-learner and T-learner, respectively, while light blue and dark
blue represent TMLE with S-learner and T-learner employed for the initial
subdistribution hazard model in the first step, respectively. Treatment and
censoring models are not involved in S-learner and T-learner. Sample size is
$n=1500$ and each scenario was replicated 500 times.

Web Figure 2.  Bias (upper panels) and RMSE (lower panels) of CATE estimations
at the median follow-up time point under different model specifications. Each
panel corresponds to one of the four subgroups determined by two binary
predictive variables, $V_{1}$ and $V_{2}$. The red line represents zero bias.
$\checkmark$ and $\times$ indicate that the corresponding models were
correctly specified and misspecified, respectively. Light grey and dark grey
represent S-learner and T-learner, respectively, while light blue and dark
blue represent TMLE with S-learner and T-learner employed for the initial
subdistribution hazard model in the first step, respectively. Treatment and
censoring models are not involved in S-learner and T-learner. Sample size is
$n=800$ and each scenario was replicated 500 times.

## Web Appendix F: Additional Analysis Results of Sepsis Dataset

Web Table 1. Baseline characteristics of patients, grouped by untreated and
treated by steroidsa.
$\begin{array}[]{l}\textbf{Baseline}\\\ \textbf{Variable}\end{array}$
$\begin{array}[]{c}\textbf{No Steroids}\\\ \textbf{(N=2,610)}\end{array}$
$\begin{array}[]{c}\textbf{Steroids}\\\ \textbf{(N=635)}\end{array}$ P-valuea
Age (Year) 63.8 (16.6) 62.9 (14.8) 0.193 Gender - Female 1063 (40.7%) 313
(49.3%) 0.04 Race White 2081 (79.7%) 505 (79.5%) 0.538 Black 289 (11.1%) 78
(12.3%) Other 240 (9.2%) 52 (8.2%) BMI 29.9 (8.92) 29.4 (9.14) 0.158 Surgery -
Yes 1003 (38.4%) 190 (29.9%) $<$0.001 Sequential Organ Failure Assessment
(SOFA) 10.1 (3.26) 11.7 (3.85) $<$0.001 Charlson Comorbidity Index (CCI) 2.89
(2.33) 3.22 (2.37) 0.001 Glasgow Coma Scale (GCS) 7.07 (3.88) 7.08 (4.11) 0.97
Serum Lactate 4.82 (4.04) 6.83 (5.82) $<$0.001
a For continuous variables, we report mean (SD) and compare using t-test. For
categorical variables, we report frequency (%) and compare using Chi-squared
test.

Web Figure 3. Cumulative incidence plot of death and discharge from ICU within
10 days post-ICU admission, grouped by treatment assignment and (A) age group,
(B) surgery status, and (C) Glasgow Coma Scale (GCS).

Web Figure 4. Variable importance for predictive and prognostic variables,
estimated by $\text{VIM}_{1}$, represents the difference in empirical variance
of the treatment effect when each variable is ignored during analysis.