key: cord-0627691-ofuhi0p5 authors: D'iaz, Iv'an; Hoffman, Katherine L; Hejazi, Nima S. title: Causal survival analysis under competing risks using longitudinal modified treatment policies date: 2022-02-07 journal: nan DOI: nan sha: ae6d676179f1a83445e57cd22cd522c042b2c8d9 doc_id: 627691 cord_uid: ofuhi0p5 Longitudinal modified treatment policies (LMTP) have been recently developed as a novel method to define and estimate causal parameters that depend on the natural value of treatment. LMTPs represent an important advancement in causal inference for longitudinal studies as they allow the non-parametric definition and estimation of the joint effect of multiple categorical, numerical, or continuous exposures measured at several time points. We extend the LMTP methodology to problems in which the outcome is a time-to-event variable subject to right-censoring and competing risks. We present identification results and non-parametric locally efficient estimators that use flexible data-adaptive regression techniques to alleviate model misspecification bias, while retaining important asymptotic properties such as $sqrt{n}$-consistency. We present an application to the estimation of the effect of the time-to-intubation on acute kidney injury amongst COVID-19 hospitalized patients, where death by other causes is taken to be the competing event. In survival analysis, a competing event is one whose realization precludes the occurrence of the event of interest. As an example, consider a study on the effect of mechanical ventilation of hospitalized COVID-19 patients on the onset of acute kidney injury, where the death of a patient by extraneous causes may be considered a competing event. Commonly used methods for analyzing time-to-event outcomes under competing events include estimating the causespecific hazard function (Prentice et al., 1978) , the sub-distribution function (also known as cause-specific cumulative incidence; Fine & Gray, 1999) , or the marginal cumulative incidence that treats the competing event as a censoring event. Early techniques for the estimation of these parameters were developed in the context of the Cox proportional hazards models (Cox, 1972; Andersen & Gill, 1982) . While these methods have allowed much progress in applied research, they have two important limitations: (i) they lack a formal framework that clarifies the conditions required for a causal interpretation of the estimates, and (ii) they typically rely upon parametric modeling assumptions that can induce non-negligible bias in the effect estimates, even when the assumptions required for causal identification hold. Recent developments in causal inference and semiparametric estimation have yielded important progress in addressing these significant limitations. For example, when the exposure is binary, Young et al. (2020) present a study on the definition and identification of several causal effect parameters in the presence of competing risks, and Benkeser et al. (2018) and Rytgaard & van der Laan (2021) describe arXiv:2202.03513v1 [stat.ME] 7 Feb 2022 semiparametric-efficient estimators that leverage the flexibility of data-adaptive regression procedures in order to mitigate model misspecification bias while retaining asymptotic properties critical to reliable statistical inference (e.g., √ n-consistency). Among other insights, Young et al. (2020) clarify that identification of the marginal cumulative incidence entails considering interventions that eliminate the competing event, and therefore require no unmeasured confounding of the competing event-outcome relation. This assumption is not required to identify the cause-specific cumulative incidence, which is the focus of this work. Defining causal effects requires considering counterfactual outcomes under hypothetical interventions on the exposure of interest. In the case of binary exposures, causal effects are often defined as differences between the expected outcomes under hypothetical interventions that assign the exposure to all units or to none -that is, the average treatment effect. While useful, causal effects defined through such deterministic interventions are limited in many applications, including problems where the exposure of interest is numerical rather than categorical, where the exposure is multivariate, where the exposure is a time-to-event-variable, or when exposure assignment is deterministic conditional on potential confounders (cf. positivity/overlap assumption). As a solution to these problems, recent work has focused on stochastic interventions and modified treatment policies (Stock, 1989; Robins et al., 2004; Díaz & van der Laan, 2012; Haneuse & Rotnitzky, 2013; Young et al., 2014) . Stochastic interventions inquire what would have happened in a hypothetical world where the post-intervention exposure is a random draw from a user-specified distribution that may itself depend on the observed data distribution. Examples of stochastic interventions include interventional mediation effects (VanderWeele et al., 2014; Hejazi et al., 2022) , as well as incremental propensity score interventions that address violations of the positivity assumption (Kennedy, 2019; Wen et al., 2021) . Modified treatment policies (MTPs), on the other hand, are defined as interventions that depend on the natural value of treatment (Young et al., 2014) . MTPs allow the definition and estimation of parameters defined in terms of natural research questions, such as inquiring what the incidence of wheezing among children with asthma would have been had the concentration of PM 2.5 particles been 5% lower, or what the mortality rate due to opioid overdose would have been had naloxone access laws been implemented earlier (Rudolph et al., 2021) . Identifying assumptions for causal effects of longitudinal modified treatment policies were first articulated by Richardson & Robins (2013) and Young et al. (2014) . Recently, Díaz et al. (2021) developed sequential regression and sequentially doubly robust estimators for longitudinal modified treatment policies (LMTP). These estimators are generalizations of estimators developed by Díaz & van der Laan (2018) for MTPs in the single time-point case, and of estimators developed by Luedtke et al. (2017) and Rotnitzky et al. (2017) for static interventions in the multiple timepoint case. Here, we extend the LMTP methodology for the estimation of the cause-specific cumulative incidence under competing events, illustrating how LMTPs can be used to solve problems whose solution has remained elusive in causal inference, such as the definition and estimation of effects for continuous, multi-valued, time-valued, or multivariate exposures. In our motivating example we illustrate estimation of the effect of time-to-intubation on time to acute kidney injury, with death by other causes acting as a competing risk. Our estimation strategy is rooted in recent developments allowing the estimation of causal effects using flexible regression techniques from the machine learning literature (van der Laan & Rose, 2011 Chernozhukov et al., 2017) . Central to this theory is the notion of von-Mises expansion (e.g., von Mises, 1947; van der Vaart, 1998; Robins et al., 2009) , which together with sample-splitting techniques (Bickel, 1982; Klaassen, 1987; Zheng & van der Laan, 2011) allows the use of flexible regression methods for the estimation of nuisance parameters while still allowing the construction of efficient and √ n-consistent estimators. Assume that at each time point t = 1, . . . τ , we measure (D t , Y t , L t , A t , C t ) on each of i = 1, . . . , n study participants, where L t denotes time-varying covariates, A t denotes a vector of exposure variables at time t, C t denotes an indicator of loss-to-follow-up (equal to one if the unit remains on study at time t + 1 and equal to zero otherwise), D t is an indicator of the competing event occurring on or before time t, and Y t denotes an indicator of the event of interest occurring on or before time t. At the end of the study, we measure the outcome of interest Y τ +1 . We let D 1 = Y 1 = 0, meaning that all participants have not experienced the event of interest or the competing event at the beginning of the study. We let R t = 1{D t = 0, Y t = 0} denote an indicator that neither the event of interest nor the competing event have occurred by t. Assume monotone loss-to-follow-up so that C t = 0 implies C k = 0 for all k > t, in which case all data become degenerate for k > t. In this work we are interested in a scenario where occurrence of the competing event precludes occurrence of the event of interest. For example, in a study looking to assess the effect of hypertension medications on the likelihood of a stroke, death for reasons other than a stroke is a competing event. In this notation, we have Y t = 1 if Y t−1 = 1, D t = 1 if D t−1 = 1, and Y t = 0 if D t = 1 and Y t−1 = 0, since the competing event precludes occurrence of the event of interest. Let X i denote the observed data for subject i. We let Pf = f (x) dP(x) for a given function f (x). We use P n to denote the empirical distribution of X 1 , . . . X n , and assume P is an element of the nonparametric statistical model defined as all continuous densities on X with respect to a dominating measure ν. We let E denote the expectation with respect to P, i.e., E{f (X)} = f (x) dP(x). We also let ||f || 2 denote the L 2 (P) norm f 2 (x) dP(x). We useW t = (W 1 , . . . , W t ) to denote the past history of a variable W , useW t = (W t , . . . , W τ ) to denote the future of a variable, and use H t = (L t ,Ā t−1 ) to denote the history of the time-varying covariates and the treatment process up until just before time t.W denotes the entire history (W 1 , . . . , W τ ) of a variable. By convention, variables with an index t ≤ 0 are defined as the null set, expectations conditioning on a null set are marginal, products of the type We formalize the definition of the causal effects using a non-parametric structural equation model (Pearl, 2000) . Specifically, for each time point t, we assume the existence of deterministic functions , . . . , τ + 1}) is a vector of exogenous variables with unrestricted joint distribution. Without loss of generality, we assume that variables which are undefined are equal to zero (e.g., L t = 0 if C k = 1 for any k < t). We are concerned with the definition and estimation of the causal effect of an intervention on the treatment processĀ on the event indicator Y τ +1 in a hypothetical world where there is no loss to follow-up (i.e.,C = 0). The interventions on the treatment process are defined in terms of longitudinal modified treatment policies (Díaz & van der Laan, 2012; Haneuse & Rotnitzky, 2013) , which are hypothetical interventions where the exposure is assigned as a new random variable A d t instead of being assigned according to the structural equation model. An intervention that sets the exposures up to time t − 1 toĀ d t−1 generates a counterfactual variable A t (Ā d t−1 ), which is referred to as the natural value of treatment (Richardson & Robins, 2013; Young et al., 2014) , and represents the value of treatment that would have been observed at time t under an intervention carried out up until time t − 1 but discontinued thereafter. An intervention on all the treatment and censoring variables up to t = τ generates a counterfactual outcome Y τ +1 (Ā d ). Causal effects are defined as contrasts between the counterfactual probability P[Y τ +1 (Ā d ) = 1] of a failure event by the end of the study under interventions implied by different interventions d, as defined below. The causal effects we define are characterized by a user-given function d(a t , h t , ε t ) that maps a given treatment value a t (i.e., the "natural value") and a history h t into a new exposure value. The function d is also allowed to depend on a randomizer ε t drawn independently across units and independently of U , and whose distribution does not depend on P. For fixed valuesā t , andl t , we recursively define a d The intervention is undefined if h d t is such that d d s = 1 or y d s = 1 for any s ≤ t. The LMTPs that we study are thus defined as To ground the ideas and illustrate the usefulness of this framework, we now Example 1 (Additive shift LMTP). Let A t denote a vector of numerical exposures, such as a drug dose or energy expenditure through physical activity. Assume that A t is supported as P Then, for a user-given vector δ, we let (1) ). This intervention has been discussed by Díaz & van der Laan (2012), Díaz & van der Laan (2018), and Haneuse & Rotnitzky (2013) in the context of single time-point exposures; by and Hejazi et al. (2022) in the context of causal mediation analysis; and by Hejazi et al. (2020) in the context of studies with two-phase sampling designs. This intervention considers hypothetical worlds in which the natural exposure at time t is increased by a user-given value δ, whenever such an increase is feasible for a unit with history H t (A d t−1 ). Example 2 (Multiplicative shift LMTP). Let A t denote a vector of numerical exposures such as pollutant concentrations for various pollutants. To define this intervention, assume that A t is supported as P The intervention implied by this function considers hypothetical worlds where the natural exposure at time t is increased by a user-given factor δ, whenever such increase is feasible for a unit with history H t (A d t−1 ). This intervention seems useful, for example, to solve current problems in environmental epidemiology related to the effect of environmental mixtures (Gibson et al., 2019) . Example 3 (Incremental propensity score interventions based on the risk ratio). Kennedy (2019) proposed an intervention tailored to binary exposures in which the post-intervention exposure is a draw from a user-given distribution g δ (a t | h t ), where this distribution is chosen such that the odds ratio comparing the likelihood of exposure under intervention versus the actual exposure is a user-given value δ. As identifying the effect of this intervention does not require the positivity assumption, these incremental propensity score interventions are very useful in settings where violations of the positivity assumption preclude identification of static effects such as the average treatment effect. Recently, Wen et al. (2021) proposed a similar intervention that instead uses the risk ratio to quantify the likelihood of exposure under the intervention. This intervention may be expressed as an LMTP as follows where ε t is a random draw from the uniform distribution on (0, 1). As before, we make use of the definition Note that the probability of exposure under the intervention is equal to g δ (1 | h t ) = δ×g(1 | h t ); thus, the risk ratio of exposure under the intervention versus the natural exposure value is δ. Importantly, a unit with zero probability of exposure conditional on their history also has zero probability of exposure under the intervention. This weakens the positivity assumptions required for identification and estimation, as we will discuss below. Example 4 (Dynamic treatment initiation strategies with grace period). Let A t denote a binary treatment variable. Consider an intervention of the form "for a user-given grace period m, if a condition for treatment initiation is met at or before t − m then start treatment at t; otherwise, do not intervene." Letting L t denote an indicator that the condition has been met, this intervention can be operationalized as an LMTP as follows: This type of intervention has been considered in the context of treatment initiation for HIV (van der Laan et al., 2005; Cain et al., 2010) , where the condition L t is that CD4+ cell count drops below a pre-specified value. In this type of application, realistic treatment policies may involve a delay of treatment initiation for logistical reasons. Accordingly, allowing for a grace period quantifies an effect that is arguably of more practical relevance. Remark 1. In contrast to LMTPs, stochastic interventions proceed by specifying an exposure density g d t (a t | h t ) a priori and considering a sequence of draws Counterfactual outcomes under stochastic interventions can then be defined as Y (Q d ). Wen et al. (2021) study stochastic interventions where the post-intervention density can be expressed as a mixture between a known density f t (a t | h t ) and the density of the observed data g t (a t | h t ). Any such intervention has an LMTP counterpart defined by where ε t is a random draw from a uniform distribution in [0, 1], 0 < c t (h t ) < 1 is a mixture constant, and q t represents a draw from f t (a t | h t ). Efficient estimators for this type of LMTP have been previously proposed by Díaz et al. (2021) . Wen et al. (2021) discuss the stochastic interventions implied by Examples 3 and 4, but incorrectly claim that the estimators they develop are not particular cases of the methods of Díaz et al. (2021) . The following assumptions will be sufficient to prove identification: Assumptions A1 and A2 are required so that the interventions considered are supported in the data. In words, Assumption A1 requires that if there is an individual with treatment level a t and covariate values h t who is event-free and uncensored at time t, there must also be an individual with treatment level d(a t , h t , ε t ) and covariate values h t who is event-free and uncensored at time t. Assumptions A2 states that for every possible value of the covariates h t that has positive probability of occurring, there are uncensored individuals. Assumption A3 is required to identify the effect of an LMTP. It states that, among event-free uncensored individuals, H t contains all the common causes of exposure at time t and future events, exposures, and covariates. Assumption A4 is required to identify the effect of a static intervention that prevents loss to follow-up. It states that, among event-free uncensored individuals, H t contains all the common causes of censoring at time t and future events and covariates. Proposition 1 (Identification of the effect of LMTPs with time-to-event outcomes subject to competing events). Set q τ +1 = Y τ +1 and R τ +1 = 1. For t = τ, . . . , 1, recursively define The identification result above, proved in the supplementary materials, extends the identification result of Young et al. (2020) to LMTPs, and the identification result of Díaz et al. (2021) to the case of competing events. In the next section we present efficiency theory and efficient estimators for θ, based on an extension of our prior work (Díaz et al., 2021) . These estimators allow the use of slow-converging data-adaptive regression techniques to obtain estimators of θ that converge at parametric rates. The general approach involves finding a von-Mises-type approximation (von Mises, 1947; van der Vaart, 1998; Robins et al., 2009; Bickel et al., 1997) for the parameter q L,t , which can be intuitively understood to be first-order expansions with second-order error (remainder) terms. As the errors in the expansion are second-order, the resulting estimator of θ will be √ n-consistent and asymptotically normal as long as the second-order error terms converge to zero at rate √ n. This would be satisfied, for example, if all the regression functions used for estimation converge at rate n 1/4 (see §3). The reader interested in this general theory is encouraged to consult the works of Buckley & James (1979) As noted by Díaz et al. (2021) , √ n-consistent estimation is not possible for arbitrary functions d when the exposure is continuous. Thus, we impose the following assumption, which is a generalization of an assumption first used by (Haneuse & Rotnitzky, 2013) . In what follows, we will assume that the function d does not depend on the distribution of the observed data P. 1 (Piecewise smooth invertibility). Assume that the cardinality of A t is p. For each h t , assume the support of A t conditional on H t = h t may be partitioned into p-dimensional subintervals I t,j (h t ) : In what follows, we assume that either (i) (i) A t is a discrete random variable for all t, or (ii) A t is a continuous random variable and the modified treatment policy d satisfies Assumption 1. Next, for continuous A t we define 1t,j{bj(at, ht, εt), ht} × gA,t{bj(at, ht, εt) | ht} × det{b j (at, ht, εt)} dP(εt), where 1 t,j {u, h t } = 1 if u ∈ I t,j (h t ) and 1 t,j {u, h t } = 0 otherwise. Under Assumption 1, it is easy to show that the p.d.f. of A d t conditional on the history h t is g d A,t (a t | h t ). In the case of Example 2, the post-intervention , which shows that piecewise smoothness is sufficient to handle interventions such as (2) which are not smooth in the full range of the exposure. Assumption 1 and expression (6) ensure that we can use the change of variable formula when computing integrals over A t for continuous exposures. This is useful for studying properties of the parameter and estimators we propose. For discrete exposure variables, we let Evaluated using the definition of the incremental propensity score intervention based on the risk ratio given in Example 3, this expression yields , which clarifies that the risk ratio for receiving the exposure under the intervention vs the observed data is equal to δ. In what follows, it will be useful to define the parameter and the function for t = τ, . . . , 1. When necessary, we use the notation ϕ t (z; η) or ϕ t (z; η t ) to highlight the dependence of ϕ t on η t = (w t , m t , . . . , w τ , m τ ). We also use η to denote (w 1 , m 1 , . . . , w τ , m τ ), and define ϕ τ +1 (Z; η) = Y . Estimation of θ will proceed by regression of the transformation ϕ t+1 sequentially from t = τ to 0 = 1 to obtain estimates m t . We say that ϕ t+1 is a multiply robust unbiased transformation for m t due to the following proposition: Proposition 2 (Doubly robust transformation). Let η be such that either m s = m s or w s = w s for all s > t. Then, in the event R t = 1, we have This proposition is a straightforward application of the von-Mises expansion given in Lemma S1 in the supplementary materials. Another consequence of Lemma S1 is that the efficient influence function for estimating θ = E[m 1 (A d , L 1 )] in the non-parametric model is given by ϕ 1 (Z) − θ (see Díaz et al., 2021) . Proposition 2 motivates the construction of the sequential regression estimator by iteratively regressing an estimate of the data transformation ϕ t+1 (Z; η) on (A t , H t ) among individuals with C t = R t = 1, starting at ϕ τ +1 (Z; η) = Y . In order to avoid imposing entropy conditions on the initial estimators, we use sample splitting and crossfitting (Klaassen, 1987; Zheng & van der Laan, 2011; Chernozhukov et al., 2018) . Let V 1 , . . . , V J denote a random partition of the index set {1, . . . , n} into J prediction sets of approximately the same size. That is, V j ⊂ {1, . . . , n}; J j=1 V j = {1, . . . , n}; and V j ∩ V j = ∅. In addition, for each j, the associated training sample is given by T j = {1, . . . , n} \ V j . We letη j denote the estimator of η obtained by training the corresponding prediction algorithm using only data in the sample T j . Further, we let j(i) denote the index of the validation set containing unit i. For preliminary cross-fitted estimatesŵ 1,j(i) , . . . ,ŵ τ,j(i) , the estimator is defined as follows: Step 1: Initialize ϕ τ +1 (Z i ;η τ,j(i) ) = Y i for i = 1, . . . , n. Step 2: For t = τ, . . . , 1: • Compute the pseudo-outcomeY t+1,i = R t+1 × ϕ t+1 (Z i ;η t,j(i) ) for all i = 1, . . . , n. • For j = 1, . . . , J: -RegressY t+1,i on (A t,i, H t,i ) using any regression technique and using only data points i ∈ T j . -Letq t,j denote the output, updateη t,j = (ŵ t,j ,m t,j , . . . ,ŵ τ,j ,m τ,j ), and iterate. Step 3: Define the sequential regression estimator aŝ θ = 1 n n i=1 ϕ 1 (Z i ;η j(i) ). Implementation of the above estimator requires an algorithm for estimation of the parameters w t . This parameter may be estimated component-wise as follows. First, the censoring mechanism g C,t may be estimated through any regression or classification procedure by regressing the indicator C t on data (A t , H t ) among individuals with R t = C t = 1. Thus, it only remains to estimate the density ratio g d t (a t | h t )/g t (a t | h t ). An option to estimate this density ratio is to first obtain an estimate of the density g t (a t | h t ), and then to obtain an estimator of g d t (a t | h t ) by applying formula (6) or (7). This approach is highly problematic in the case of continuous or multivariate exposures due both to the curse of dimensionality and to the fact that flexible and computationally efficient estimation of conditional densities is an under-developed area of machine learning (relative to regression). This approach can also be problematic for discrete exposures of high cardinality. Since the conditional densities themselves are unnecessary for estimation, and direct estimation of the ratio g d t (a t | h t )/g t (a t | h t ) suffices, we propose to use Díaz et al. (2021) 's method for direct estimation of density ratios. Briefly, this approach starts by constructing an augmented artificial dataset of size 2n in which each observation is duplicated. For each observation, one duplicate gets assigned the post-intervention exposure value d(A t , H t ) and the other gets the actual observed value A t . Then, the task of density ratio estimation reduces to the task of classifying which record belongs to the intervened-upon unit, based on data on the exposure and the history H t . Specifically, Díaz et al. (2021) show that the odds of this classification problem is equal to the density ratio. The above estimation algorithm requires the implementation of several regression procedures. Many algorithms from the machine learning literature are available for such regression problems, including those based on ensembles of regression trees, splines, etc. As stated in the asymptotic normality result below, it is important that the regression procedures used approximate as well as possible the relevant components of the true data-generating mechanism. Increasing the predictive performance of these regression fits can entail estimator selection from a large list of candidate estimators. If the sample size is large enough, a principled solution to this problem is the use of Super Learning , a form of model stacking (Wolpert, 1992; Breiman, 1996) that constructs an ensemble model as a convex combination of candidate algorithms from a user-supplied regression library, where the weights in the combination are computed based on cross-validation and are such that they minimize the prediction error of the resulting estimator. The asymptotic distribution of θ is Gaussian with variance equal to the efficiency bound, under conditions. To state the conditions, we first define the data-dependent parameteř where the outer expectation is only with respect to the distribution P of Z (i.e.,η is fixed). The following proposition is a consequence of the results of Díaz et al. (2021) : Proposition 3 (Consistency and weak convergence of SDR estimator). We have where σ 2 = Var{ϕ 1 (Z; η)} is the non-parametric efficiency bound. Note that this means thatθ is sequentially doubly robust in the sense that, at each time point, it is robust to misspecification of at most one of the two regression procedures used. This property is also known in the literature under the name of 2 τ -multiple robustness. We prefer the term sequentially doubly robust for two reasons. First, contrasted to the double robustness property which requires one out of two nuisance estimators to be consistent, the term 2 τ -multiple robustness erroneously conveys the message that consistency would follow from one out of 2 τ nuisance estimators being consistent. Second, retaining the term "double" in this description conveys a fundamental property of the estimation procedure: it is based on an expansion (given in Lemma S1 in the supplementary material) with second-order remainder terms. In this nomenclature, a triply robust estimating procedure would be one based on an expansion with third-order remainder terms, and so forth. 1 Under the conditions of the above theorem, a Wald-type procedure, in which the standard error may be estimated based on the empirical variance of ϕ 1 (Z;η j(i) ), yields asymptotically valid confidence intervals and hypothesis tests. A drawback of the above estimation approach is that there is no guarantee that the parameter will remain within the bounds of the parameter space. That is, it is possible, in principle, that the sequential regression estimator will provide time-to-event probabilities outside of the closed unit interval [0, 1]. Furthermore, when the survival function must be estimated at several points in time, there is no guarantee that the resulting point estimates will be antitonic (i.e., monotonic decreasing). In the supplementary materials, we present a targeted minimum loss estimator (TMLE) capable of producing estimates that always lie within the unit interval. However, boundedness of the TMLE comes with a robustness tradeoff. In particular, the robustness properties of the TMLE are inferior to those of the sequentially doubly robust estimator. Specifically, consistency of the TMLE requires that either ŵ t − w t = o P (1) or q t −q t = o P (1), for an estimatorq t constructed based on sequential regression using (5). As a result, consistent estimation of q t in the TMLE procedure requires consistent estimation of q s : s > t. By contrast, the sequentially doubly robust estimator provides consistent estimation of q t when either q s or w s is consistently estimated for all s > t. Both Luedtke et al. (2017) and Díaz et al. (2021) provide further in-depth discussion on this topic. To resolve the above issues, we leverage a set of techniques proposed by . Under their approach, out-of-bounds estimates are truncated to remain within bounds, and the possibly non-monotonic estimate of the survival curve, alongside any confidence region limits, is projected onto the space of monotonic functions by way of isotonic regression. Importantly, show that the resulting projected estimator retains the same asymptotic properties as the original estimator. In our case, this implies that the projected estimator retains the properties of asymptotic normality and consistency under the assumptions of Proposition 3. To illustrate the proposed methods, we aim to answer the following question: what is the effect of invasive mechanical ventilation (IMV) on acute kidney injury (AKI) among COVID-19 patients? This question has been recently identified by an expert panel on lung-kidney interactions in critically ill patients (Joannidis et al., 2020) as an area in need of further research. AKI is a common condition in the ICU, complicating about 30% of ICU admissions, and causing increased risk of in-hospital mortality and long-term morbidity and mortality (Kes & Jukić, 2010) . AKI is often viewed as a tolerable consequence of interventions to support other failing organs, such as IMV (Husain-Syed et al., 2016 ). To answer this question, we will use data on approximately 3,300 patients hospitalized with COVID-19 at the New York Presbyterian Cornell, Queens, or Lower Manhattan hospitals. Our analysis focuses on patients hospitalized between March 3 and May 15, 2020, who did not have a previous history of Chronic Kidney Disease (CKD). This choice of time frame is rooted in there being a lack of treatment guidelines in the initial months of the COVID-19 pandemic, which resulted in a large variability in provider practice regarding mechanical ventilation. This clinical heterogeneity in the deployment of IMV makes the problem particularly well-suited to assessment via the causal effects of LMTPs. The analytical dataset was built from two distinct databases. Demographics, comorbidity, intubation, death, and discharge data were abstracted from electronic health records by trained medical professionals into a secure RedCAP database (Goyal et al., 2020) . These data were supplemented with the Weill Cornell Critical carE Database for Advanced Research (WC-CEDAR), a comprehensive data repository housing laboratory, procedure, diagnosis, medication, microbiology, and flowsheet data documented as part of standard care . This dataset is unique in that it contains a large sample of patients treated in the same hospital system within a short period of time, while still exhibiting extensive treatment variability. The most frequent reason for morbidity and mortality among patients with COVID-19 is pneumonia leading to acute respiratory distress syndrome (ARDS). Respiratory support via devices such as nasal cannulae, face masks, or endotracheal tubes (i.e., intubation and subsequent mechanical ventilation) is a primary treatment for ARDS . For severe ARDS patients, an important decision must be made regarding the optimal timing of intubation based on clinical factors such as oxygen saturation, dyspnea, respiratory rate, chest radiograph, and, occasionally, external resource considerations such as ventilator or provider availability (Tobin, 2020; Thomson & Calligaro, 2021) . At the outset of the COVID-19 pandemic, multiple international guidelines recommended early intubation in an effort to protect healthcare workers from infection and to avoid complications stemming from "crash" intubations (Papoutsi et al., 2021) . Yet, as the pandemic progressed and multiple reports documented high mortality for mechanically ventilated patients, guidance changed to delaying intubation with the rationale that IMV increases risk of secondary infections, ventilator-induced lung injury, damage to other organs, and death (Bavishi et al., 2021; Tobin, 2006) . The lungs and kidneys are physiologically closely tied, and it has been hypothesized that mechanical ventilation may cause AKI in COVID-19 patients via oxygen toxicity and capillary endothelial damage leading to inflammation, hypotension, and sepsis (Durdevic et al., 2020) . Despite these hypotheses, a causal link between intubation and AKI in patients with COVID-19 has not been definitively established. Furthermore, such a link is difficult to study in a randomized fashion due to the clinically subjective nature of intubation assignment -that is, it is impossible to pick an oxygen saturation breakpoint to intubate across all patients and hospital settings Perkins et al., 2020) . For each patient, study time begins on the day of hospitalization. The exposure of interest is daily maximum respiratory support, which can take three levels: no supplemental oxygen, supplemental oxygen not including IMV, and IMV. Our goal is to estimate the overall causal effect on daily AKI rates of an intervention that would delay intubation for IMV by a single day among 449 patients who received IMV in the first 14 days of hospitalization, where death by other causes is a competing risk for AKI development. In notation, this LMTP may be expressed as follows. Let A t ∈ {0, 1, 2}, where 0 indicates no supplemental oxygen at end of day t, 1 indicates supplemental oxygen not including IMV, and 2 indicates IMV. Consider the following intervention: Under this intervention, a patient who is first intubated at day t would receive non-invasive oxygen support at day t. Otherwise, the patient would receive the oxygen support actually (i.e., "naturally") observed at day t. This intervention assesses what would have happened in a hypothetical world where intubation procedures were more conservative than they were at the beginning of the pandemic in New York City's surge conditions, when limited information and treatments were available for COVID-19. Figure 1 shows the observed supplemental oxygen levels and observed AKI and death rates across days 0-14. Baseline confounders include age, sex, race, ethnicity, body mass index (BMI), hospital location, home oxygen status, and comorbidities (e.g., hypertension, history of stroke, Diabetes Mellitus, Coronary Artery Disease, active cancer, cirrhosis, asthma, Chronic Obstructive Pulmonary Disease, Interstitial Lung Disease, HIV infection, and immunosuppression). Time-dependent confounders include vital signs (e.g., highest and lowest respiratory rate, oxygen saturation, temperature, heart rate, and blood pressure), laboratory results (e.g., C-Reactive Protein, BUN Creatinine Ratio, Activated Partial Thromboplastin time, Creatinine, lymphocytes, neutrophils, bilirubin, platelets, D-dimer, glucose, arterial partial pressure of oxygen, and arterial partial pressure of carbon dioxide), and an indicator for daily corticosteroid administration greater than 0.5 mg/kg body weight. In cases of missing baseline confounders, multiple imputation using chained equations (MICE) (van Buuren & Groothuis-Oudshoorn, 2011) is used. For missing data at later time points, the last observation is carried forward. Patients are censored at their day of hospital discharge, as AKI and vital status were unknown after this point. We estimate the effect of the LMTP with the SDR estimator (on account of its enhanced robustness relative to TMLE), using a version of the lmtp R package (Williams et al., 2022) modified to support the competing risks setting. The nuisance functions r t and m t are estimated using the Super Learner ensemble modeling algorithm via the sl3 R package (Coyle et al., 2022) . The candidate library considered by the Super Learner included a total of 32 learning algorithms, inclusive of hyperparameter variations. The algorithms included multivariate adaptive regression splines (MARS; Friedman et al., 1991) ; random forests (Breiman, 2001; Wright & Ziegler, 2017) ; extreme gradient boosted trees (Friedman, 2001; Chen & Guestrin, 2016) ; lasso, ridge, and elastic net (with equal weights, i.e., α = 0.5) regularized regression (Friedman et al., 2022) ; main terms generalized linear models (McCullagh & Nelder, 1989; Enea, 2009) ; generalized linear models with Bayesian priors on parameter estimates (Gelman & Hill, 2006) ; and an unadjusted outcome mean. The estimators of r t and m t were computed under a Markov assumption with lag r = 2 to reflect the collection of laboratory results at minimum 48-hour intervals. Figure 2 displays the SDR-based cumulative event curve estimates for each treatment condition (panel A) after the correction for monotonicity proposed by , and the incidence difference between the two (panel B), together with 95% simultaneous confidence bands. Examination of Figure 2 panel A reveals a slight but increasingly widening gap (across days) between the SDR-based incidence estimates for each arm. In agreement with this, inspection of panel B of Figure 2 makes apparent the protective effect of the delay-in-intubation policy against AKI. The estimates of the incidence difference are consistently negative and decrease steadily in time; moreover, the simultaneous confidence bands do not overlap with the null value Days since hospital admission Estimated incidence difference Estimated AKI Incidence Difference via SDR B Figure 2 : Panel A: Estimated AKI incidence across the arms defined by the two treatment strategies, 1-day delay-ofintubation and no intervention. Panel B: The time-specific difference between the two treatment conditions. In both panels, 95% simultaneous confidence bands cover the sets of point estimates. of zero except for on the first day. Hypothesis testing of the incidence difference yields adjusted p-values < 0.001 for all days after the first, indicating a statistically significant protective effect of the LMTP; the Bonferroni correction was used to adjust for testing multiplicity. Based on these estimates, the delay-in-intubation policy would be expected to result in a decrease in 14-day AKI rates by approximately 3.5% (simultaneous 95% CI: [2.7%, 4.3%]) relative to no intervention. Modified treatment policies represent a recent but powerful advance in modern causal inference. These intervention regimes have a variety of applications, including (i) practical solutions to violations of the positivity assumption; (ii) definition and estimation of the effects of continuous, time-valued, and multivariate exposures; and (iii) definition of effects that depend on the natural value of treatment. Here, we have advanced the LMTP methodology, extending it to the case of survival outcomes with competing events. We have presented identification assumptions for the causal parameters, as well as a sequentially doubly robust and efficient estimation algorithm capable of leveraging the many data-adaptive regression procedures available for nuisance parameter estimation. We ensure that the estimates of the survival function remain within the parameter space (antitonic and bounded in the closed unit interval) by leveraging the monotonic projection framework of and . We have illustrated the scientific utility of our approach through application to the study of negative side effects of invasive mechanical ventilation on acute kidney injury in patients hospitalized with COVID-19, for which death by other causes serves as a competing event. Our Proposition 2 shows that doubly robust estimation is possible for any intervention that is defined in terms of an LMTP satisfying the stated conditions. Namely, any randomized intervention defined by a known function d(a t , h t , ε t ) through either an LMTP or a general stochastic intervention accommodates doubly robust estimation, provided that: (1) either (i) A t is a discrete random variable for all t, or (ii) A t is a continuous random variable and the modified treatment policy d satisfies Assumption 1, and (2) the randomizer ε t is drawn independently across units and independently of U , and its distribution does not depend on P. Notably, this rich class contains all interventions considered by Theorem 2 of Wen et al. (2021) . Cox's regression model for counting processes: a large sample study Doubly robust estimation in missing data and causal inference models Timing of intubation in coronavirus disease 2019: A study of ventilator mechanics, imaging, findings, and outcomes Improved estimation of the cumulative incidence of rare outcomes On adaptive estimation Efficient and Adaptive Estimation for Semiparametric Models Stacked regressions Random forests Linear regression with censored data When to start treatment? a systematic approach to the comparison of dynamic regimes using observational data Xgboost: A scalable tree boosting system Double/debiased/neyman machine learning of treatment effects Double/debiased machine learning for treatment and structural parameters Regression models and life-tables sl3: Modern pipelines for machine learning and Super Learning Causal mediation analysis for stochastic interventions Non-parametric efficient causal mediation with intermediate confounders Population intervention causal effects based on stochastic interventions Stochastic treatment regimes Nonparametric causal effects based on longitudinal modified treatment policies Progressive renal failure in patients with covid-19 after initiating mechanical ventilation: A case series Fitting linear models and generalized linear models with large data sets in R. In Statistical Methods for the Analysis of Large Datasets A proportional hazards model for the subdistribution of a competing risk glmnet: Lasso and elastic-net regularized generalized linear models Greedy function approximation: a gradient boosting machine Multivariate adaptive regression splines Data Analysis Using Regression and Multilevel/Hierarchical Models An overview of methods to address distinct research questions on environmental mixtures: an application to persistent organic pollutants and leukocyte telomere length Clinical characteristics of covid-19 in new york city Estimation of the effect of interventions that modify the received treatment Mortality in covid-19 patients with acute respiratory distress syndrome and corticosteroids use: a systematic review and meta-analysis Nonparametric causal mediation analysis for stochastic interventional (in) direct effects Efficient nonparametric inference on the effects of stochastic interventions under two-phase sampling, with applications to vaccine efficacy trials Lung-kidney cross-talk in the critically ill patient Lung-kidney interactions in critically ill patients: consensus report of the acute disease quality initiative (adqi) 21 workgroup Nonparametric causal effects based on incremental propensity score interventions Acute kidney injury in the intensive care unit Consistent estimation of the influence function of locally asymptotically linear estimators Sequential double robustness in right-censored longitudinal models Generalized Linear Models Effect of timing of intubation on clinical outcomes of critically ill patients with covid-19: a systematic review and meta-analysis of non-randomized cohort studies Causality: Models, Reasoning, and Inference Recovery-respiratory support: Respiratory strategies for patients with suspected or proven covid-19 respiratory failure; continuous positive airway pressure, high-flow nasal oxygen, and standard care: A structured summary of a study protocol for a randomised controlled trial The analysis of failure times in the presence of competing risks Single world intervention graphs (SWIGs): A unification of the counterfactual and graphical approaches to causality Quadratic semiparametric von Mises calculus Robust estimation in sequentially ignorable missing data and causal inference models Effects of multiple interventions. Comparative quantification of health risks: global and regional burden of disease attributable to selected major risk factors 1 Estimation of regression coefficients when some regressors are not always observed On the multiply robust estimation of the mean of the g-functional A doubly robust censoring unbiased transformation When effects cannot be estimated: redefining estimands to understand the effects of naloxone access laws One-step tmle for targeting cause-specific absolute risks and survival curves Critical care database for advanced research (cedar): An automated method to support intensive care units with electronic health record data Nonparametric policy analysis Timing of intubation in COVID-19: Not just location, location, location? Principles and practice of mechanical ventilation Basing respiratory management of COVID-19 on physiological principles Caution about early intubation and mechanical ventilation in covid-19 mice: Multivariate imputation by chained equations in R The cross-validated adaptive epsilon-net estimator History-adjusted marginal structural models and statically-optimal dynamic treatment regimens Super learner Unified Methods for Censored Longitudinal Data and Causality Targeted Learning: Causal Inference for Observational and Experimental Data Targeted Learning in Data Science: Causal Inference for Complex longitudinal Studies Targeted maximum likelihood learning Asymptotic Statistics Effect decomposition in the presence of an exposure-induced mediator-outcome confounder On the asymptotic distribution of differentiable statistical functions Intervention treatment distributions that depend on the observed treatment process and model double robustness in causal survival analysis A unified study of nonparametric inference for monotone functions Correcting an estimator of a multivariate monotone function with isotonic regression lmtp: Non-parametric causal effects of feasible interventions based on modified treatment policies Stacked generalization ranger: A fast implementation of random forests for high dimensional data in C++ and R Identification, estimation and approximation of risk under interventions that depend on the natural value of treatment using observational data A causal framework for classical statistical estimands in failure-time settings with competing events Cross-validated targeted minimum-loss-based estimation The authors declare that they have no conflict of interest.