Estimation of Heterogeneous Restricted Mean Survival Time Using Random Forest

Estimation and prediction of heterogeneous restricted mean survival time (hRMST) is of great clinical importance, which can provide an easily interpretable and clinically meaningful summary of the survival function in the presence of censoring and individual covariates. The existing methods for the modeling of hRMST rely on proportional hazards or other parametric assumptions on the survival distribution. In this paper, we propose a random forest based estimation of hRMST for right-censored survival data with covariates and prove a central limit theorem for the resulting estimator. In addition, we present a computationally efficient construction for the confidence interval of hRMST. Our simulations show that the resulting confidence intervals have the correct coverage probability of the hRMST, and the random forest based estimate of hRMST has smaller prediction errors than the parametric models when the models are mis-specified. We apply the method to the ovarian cancer data set from The Cancer Genome Atlas (TCGA) project to predict hRMST and show an improved prediction performance over the existing methods. A software implementation, srf using R and C++, is available at https://github.com/lmy1019/SRF.


INTRODUCTION
In epidemiological and biomedical studies, time to an event or survival time T is often the primary outcome of interest. Important quantities related to survival time include hazard rate (HR), tyear survival probability, and the mean survival time. Among these, HR is one of the most commonly used quantity due to its strong connection to the proportional hazards regression model or Cox model. Cox model is a very popular regression model for censored survival data due to its computational feasibility and theoretical properties (Cox, 1972(Cox, , 1975Andersen and Gill, 1982;Gill and Gill, 1984;Huang et al., 2013;Fang et al., 2017). However, when there is a departure from the proportional hazards assumption, the connection between HR and survival function is lost and it is difficult to interpret HR (Wang and Schaubel, 2018). The t-year survival probability is the probability of survival time greater than a pre-specified time t. It is not suitable for summarizing the global profile of T over the duration of a study (Tian et al., 2014). In contrast, mean survival time is an alternative quantity since it takes the whole distribution of T into account. However, the mean of T may not always be estimable in the presence of censoring. For example, let C denotes the FIGURE 1 | Training data are simulated from Equation (2), with n = 600 training points, dimension p = 20 and errors ǫ ∼ N(0, 10 2 ). Random forests are trained based using R package grf. Truth is shown as red curve, with green curve corresponding to the random forest predictions, and upper and lower bounds of the point-wise confidence intervals connected in the black lines. Brown curve and blue curve are based on the approaches of Wang and Schaubel (2018) with Identity and Exp link functions.
censoring time, and C max = inf c {P(C ≤ c) = 1} be the upper limit of the censoring distribution, If the survival time T satisfies P(T > C max ) > 0, then we cannot estimate E T [T], since we never observe any event after C max .
The restricted mean survival time (RMST) (Royston and Parmar, 2013) summarizes the survival process and provides an attractive alternative to the proportional hazards regression model (Tian et al., 2014). The restricted survival time of T up to a fixed point L is defined as T ∧ L, and the restricted mean survival time is defined as the expectation of the restricted survival time. (1) If L is chosen to be less than C max , hRMST is estimable since P(T ∧ L > C max ) = 0. RMST also plays a role in the context of inverse probability censoring weighting (IPCW). A key assumption for applying IPCW is P(T < C max ) = 1, making 1/(1 − G(T)) well-defined, where G(T) = P(C ≤ T|T). If we set L properly such that P(T ∧ L < C max ) = 1, then G(T ∧ C ∧ L|X) < 1 and the IPCW is well-defined under the restricted survival time context. There are two main approaches for hRMST regression. One approach is to estimate hRMST indirectly through hazard regression (Zucker, 1998;Chen and Tsiatis, 2001;Zhang and Schaubel, 2011). This approach starts by estimating the regression parameters and the baseline hazard from a Cox model, calculating the cumulative baseline hazard, transforming it to obtain the survival function and, finally, obtaining the hRMST through Equation (1). Such an indirect hRMST estimation is inconvenient and computationally cumbersome for obtaining a point estimate and its corresponding asymptotic standard error. An alternative approach is to model hRMST with the baseline covariates X directly via some parametric assumptions, eg. g[µ L (X i )] = β ′ 0 X i , where g is a strictly monotone link function with a continuous derivative within an open neighborhood (Tian et al., 2014;Wang and Schaubel, 2018). A major weakness of this approach, however, is their inability to choose a proper link function, which may lead to the model misspecification. As an example, we simulate x 1 , . . . , x n independently from the uniform distribution on [0, 1] 20 with a survival time model where we assume that the censoring time C and the restricted time L satisfy P(C ≤ T ∧L) = 33% and P(L ≤ T ∧C) = 11%. Our goal is to estimate µ L (x). Figure 1 shows a set of predictions on an artificially generated data set from Equation (2). Compared with other methods, the random forest is able to estimate the target function closely, especially when µ L (x) approaches L.
For the continuous outcomes without censoring, random forest (Breiman, 2001(Breiman, , 2004) is a popular method of nonparametric regression that has shown effectiveness in many applications (Svetnik et al., 2003;Díaz-Uriarte and Alvarez de Andrés, 2006;Cutler et al., 2007). It is invariant under scaling and various other transformations of feature values, robust to inclusion of irrelevant features (Hastie et al., 2001), and versatile enough to be applied to large-scale problems (Biau and Scornet, 2016). Besides strong empirical results, theoretical results such as consistency (Meinshausen, 2006;Biau et al., 2008;Biau, 2012;Denil et al., 2014) and asymptotic normality (Wager and Athey, 2015;Mentch and Hooker, 2016;Athey et al., 2018;Friedberg et al., 2018) have also been obtained for regression models without censoring. Extending random forest to censored survival data has been proposed in several recent papers (Ishwaran et al., 2008;Steingrimsson et al., 2019), focusing on implementations and algorithms. However, there has been little theoretical work in statistical inference of such random survival forest. Ishwaran and Kogalur (2011) proved the consistency of the random survival forest by showing that the forest ensemble survival function converges uniformly to the true population survival function.
Instead of focusing on predicting the survival function or the survival probability as the algorithms implemented by Ishwaran et al. (2008) and Steingrimsson et al. (2019), we develop in this paper a random forest framework to model the hRMST directly given the baseline covariates in the presence of possibly covariate-dependent censoring. This approach provides a nonparametric estimation of hRMST adjusting for covariates. Due to the complex relationship between the survival time and the covariates, it is desirable to have more flexible methods to estimate the hRMST than the approaches that a certain link function has to be assumed. Our construction of random forest is based on the estimated IPCW. We show that the resulting survival random forest estimates of hRMST has the asymptotic normality property that can be used to obtain the point-wise confidence interval with theoretical guarantees. To the best of our knowledge, it is the first asymptotic normality result for the predictions in the context of censored survival data using random forest.
The remainder of the paper is organized as follows. In section 2, we describe the proposed random forest estimator. Asymptotic properties are given in section 3. In section 4, we conduct simulation studies to evaluate the accuracy of the proposed method in the finite sample settings. In section 5, we apply our method to an ovarian cancer data set of The Cancer Genome Atlas (TCGA) project (http://cancergenome.nih.gov/abouttcga) to evaluate the predictions of the hRMST for ovarian cancer patients using their acylcarnitine measurements and clinical variables. We conclude this chapter with a brief discussion in section 6.

RANDOM FOREST FOR ESTIMATING THE hRMST
We begin with some notation. Let X i be the baseline covariates for subject i from a cohort of sample size n and T i be the survival time for subject i. Let C i be the censoring time, which is independent of T i conditional on the baseline covariates X i . The observation time for subject i is Z i = T i ∧ C i , where a ∧ b = min{a, b}. The indicator for censoring is denoted by Let L be a pre-specified time point of interest, before the maximum follow-up time τ = max{Z i : i = 1, . . . , n}. As in Wang and Schaubel (2018), L is normally chosen as a time point of clinical relevance or, at least, of particular interest to the investigators, respecting the bound at the maximum follow-up time. Denote the restricted observation time as Z L i = Z i ∧L and its corresponding indicator δ L i = 1 {T i ∧L≤C i } . Our goal is to estimate covariate-adjusted RMST or hRMST µ L (x) = E(Z L |X = x) and to construct its confidence interval.

Forest-Based Local Estimating Equation for hRMST
, and a restriction threshold L, we first present a random forest method to estimate µ L (x). The idea of the approach is to solve a weighted estimating equation for µ L (x), where the estimating equation functions of the observations whose covariates closer to x will have larger weights. Specifically, let w i = δ L i /(1 − G(Z L i |X i )) be the IPCW of the ith data point under the true censoring distribution G(·|X i ). The (infeasible) estimating equation function are also known, the solution to the empirical estimating equation for µ L (x) is given as which provides a good candidate of estimator for µ L (x). However we do not know the censoring distribution G and the local weights {α i (x)} n i=1 , which need to be estimated from the data. We assume censoring distribution G follows a Cox model, a natural choice for modeling censoring times in the context of IPCW. Let be the estimated IPCW for ith observation withĜ(·|X i ) derived from the data through Cox model. We define the estimating equation function for ith observation with its corresponding estimated IPCW as Our approach to derive the local weights {α i (x)} n i=1 is through the random forest, which is an ensemble of survival trees constructed by Algorithm 1.

Algorithm 1: Survival tree
SurvivalTree (set of observations J, domain X); IPCW←CoxModel(J); Root P 0 ←CreateNode(J, X); Queue Q →InitializeQueue(P 0 ); while Q is NotNull do node P ← Pop(Q); if split succeeds then AddQueue(C 1 ); AddQueue(C 2 ); end end It can be shown that ρ i is the influence function of the ith observation forμ L P . Let F n be the empirical distribution of the observations in node P, and let F n,i = (1 − ǫ)F n + ǫν i , with ν i be the Dirac delta function at ith observation. Set µ L P,i =μ L P + i , whereμ L P,i = argmin µ | ψ µ (X, Z L , δ L )dF n,i |.
By Taylor expansion, where µ * is a value betweenμ L P andμ L P,i . The above equation implies , and therefore the influence function of ith observation forμ L P is Athey et al. (2018) shows that maximizing the splitting criterion˜ (C 1 , C 2 ) is approximately equivalent to minimizing the weighted mean squared error err(C 1 , In order to achieve consistency and asymptotic normality, we split the tree and make predictions in an honest way as introduced in Wager and Athey (2015). Specifically, each tree in an honest forest is grown using two non-overlapping subsamples of the training data. For the bth tree, given I b and J b , we first choose the tree structure T b using only the data in J b , and write x ↔ b x ′ as the boolean indicator for whether the points x and x ′ fall into the same leaf of T b . In a second step, we define the set of neighbors of The weights of point x from a survival forest with B trees can be written as The empirical locally weighted estimating equation forμ L (x) is then defined as FIGURE 2 | Simulation results of the coverage probability for Model 1 with three different link functions, sample size of n = 1, 000, 2, 000, 5, 000, and p = 2, 4, 6, 8. For each case, prediction coverage probability is calculated over the samples in the testing data set.
Frontiers in Genetics | www.frontiersin.org and the random forest estimator for the hRMST is the solution of Equation (4), which iŝ We emphasize the difference between the IPCW used in building the survival trees and IPCW used to deriveμ L (x). The IPCW used in building survival trees is estimated only by the data points from J b so that the resulting survival forest is honest. The IPCW used to deriveμ L (x) is estimated from all data points.

Asymptotic Normality
We derive a central limit theorem for survival forest estimate of hRMST. We first give three common assumptions that required for the most of the theoretical analysis of random forests.
Assumption 2. There exists a restricted time threshold L, such that P(C > t ∧ L|X = x) ≥ ǫ L > 0 for any x, t.
As mentioned in the previous section, we model the conditional survival function of censoring distribution G given baseline covariates. Because of its flexibility and popularity in practice, we adopt the proportional hazards model for hazard function of censoring distribution.
We make additional regularity assumptions that are widely used in analysis of estimates from the proportional hazards models. These assumptions are needed in order to quantify the difference between the estimated IPCW and true IPCW.
FIGURE 3 | Simulation results of coverage probability for Model 2 with three different link functions, sample size of n = 1, 000, 2, 000, 10, 000, and p = 2, 4, 6, 8. For each case, prediction coverage probability is calculated over the samples in the testing data set.
for some positive constant and for any t, x. This assumption implies that Following Wager and Athey (2015) and Athey et al. (2018), we assume that all trees are symmetric, in that their output is invariant to permuting the indices of Estimation-Part in training examples (see Corollary 6 of Wager and Athey (2015) for more details about this symmetry). They also require balanced splits in the sense that every split puts at least a fraction ω of the observations in the parent node into each child, for some ω > 0. Finally, the trees are randomized in such a way that, at every split, the probability that the tree splits on the jth feature is bounded from below by some π > 0. The forest is honest and built via subsampling with subsample size s satisfying s/n → 0 and s → ∞.
Under the assumptions listed above, we have the following asymptotic distribution result for the random forest-based estimate of the hRMST. Theorem 1. Under Assumptions 1, 2, 3, 4, 5, 6, 7, 8, for each fixed test point x, there is a sequence σ 2 where ω > 0 is the low-bound fraction for observations in the parent node into each child, and π > 0 is the lower-bound of the probability that the tree splits on any features.
We give a consistent estimate of σ 2 n (x) based on half-sampling (Efron, 1980) and the method of Sexton and Laake (2009).

Estimation of the Variance
Following Athey et al. (2018), we use the random forest delta method to develop a variance estimate of the survival forest predictionμ L (x). Athey et al. (2018) provides a consistent estimate of σ 2 n (x) using s 2 In our context, V(x) = −1, then simply we have s 2 n (x) = H n (x). A consistent estimator for H n (x) can be obtained using halfsampling estimator (Efron, 1980;Athey et al., 2018). Let H be the average of the empirical estimating equation functions averaged over the trees that only use the data from the halfsample H, denoted by S H , where L b (x) contains neighbors of x in the bth tree. An ideal half-sampling estimator is then defined aŝ  Wang and Schaubel (2018) with identity link; Wang:exp, method of Wang and Schaubel (2018) with exponential link.
where is the randomness in building honest tree, including splitting data into random halves and randomness in selecting variables to split.Ĥ HS n (x) is similar to classic bootstrap estimator for the standard error, except that the sampling distribution forĤ HS n (x) is the half sampling distribution instead of the bootstrap sampling. Denote E ss and Var ss as the expectation and variance under the half sampling distribution, thenĤ HS Since carrying out the full half-sampling computation and expectation with respect to are impractical, Sexton and Laake (2009) pointed out thatĤ HS n (x) can be efficiently approximated by the following law of total variance: which leads to a Monte Carlo approximation ofĤ HS n (x) bŷ In order to approximate random forest randomness quantity Var and sampling randomness quantities Var ss ,Ê ss , we split B trees in G groups and each group has l trees, and the trees in the same group have the same half sample. The final consistent estimatorσ 2 n (x) can be written aŝ The following diagram summarizes the procedure of estimating the variance σ 2 n (x).
where from left to right, the first arrow is based on Theorem 5 of Athey et al. (2018), the second arrow is based on half-sampling of Efron (1980), and the third arrow is supported by Equations (5) and (6)

SIMULATION STUDIES
We present simulations to evaluate the performance of the proposed method in finite sample setting. Two different models for the survival time are considered . , X ip are independently generated from Unif (−1, 1), α 0 = 5, α 1 = α 2 = 0.25 and α i = 0 for i > 2, and ǫ ∼ N(0, σ 2 ). The variance σ 2 is chosen to have proper signal-noise ratio (SNR), We generate the independent censoring time C i from a Cox model with the following hazard λ = λ C exp(X 1 log 2) and λ C is chosen to have a proper un-censoring rate. The link function g can have the following form • Identity link: g −1 (x) = x; • Exp link: g −1 (x) = exp(x); • Log-exp link: g −1 (x) = log(exp(x) + 1).

Evaluation of Coverage Probability of Predictions
To evaluate the asymptotic results in Theorem 1, we generate five training data sets and one testing data set with the same sample size. The coverage probability performance is evaluated on the testing data set with predictions and confidence intervals derived from 5 independent training data sets. More specifically, for each observation in the testing sample, we obtain the 95% confidence intervals and record how many times a hRMST observation in test sample is within five estimated 95% confidence intervals. The coverage probability of an observation is defined by the its proportion of being covered, and the overall coverage probability of the testing sample is defined by the average of coverage probability of each of its observation. We present the coverage probability results with sample size n = 1, 000, 2, 000, 5, 000 for Model 1, and n = 1, 000, 2, 000, 10, 000 for Model 2. By choosing the proper λ C , we control the un-censoring rate around 60-70% for different link functions: λ C ∼ 0.08 for Identity link and Logexp link, and λ C ∼ 0.003 for Exp link. The truncation time L is  Tian et al. (2014) with identity link; Lu.exp, method of Tian et al. (2014) with exponential link; Wang.id, method of Wang and Schaubel (2018) with identity link; Wang:exp, method of Wang and Schaubel (2018) with exponential link.
chosen to make the truncation rate fall into 2% − 5%. Specifically, L ∼ 5.4 for Identity link and Log-exp link, and L ∼ 220 for Exp link. Figures 2, 3 present the results for Model 1 and Model 2 under three different link functions. We see that the coverage probability approaches to nominal level 95% when the sample size gets larger. If p is smaller, the coverage probability is closer to 95%. This corresponds to the result of Theorem 3 in Wager and Athey (2015), which states that the rate of convergence of the bias of random forest estimator is O(n K p ) for some constant K. When the sample size n is fixed, bigger p leads to larger bias in the estimates of hRMST, and under-coverage of the confidence interval. On the other hand, when p is fixed, bigger n results in a smaller bias and leads to a better coverage of the confidence interval.

Comparison of Prediction Performance With Existing Methods
We compare our proposed method with several existing methods for hRMST estimation, including • Naive.km: using Kaplan-Meier estimator for survival function and computing hRMST by Equation (1). Covariates are not adjusted. • Naive.Cox: using proportational hazards estimator for the survival function and computing hRMST by Equation (1). The censoring distribution is assumed to follow the proportional hazards assumption. • Lu.method: using some parametric forms of hRMST and computing hRMST by solving a weighted estimating equation. The censoring distribution is assumed to be independent of the covariates (Tian et al., 2014). We consider Identity link and Exp link in the simulations. • Wang.method: using some parametric forms of hRMST and computing hRMST by solving a weighted estimating equation. The censoring distribution is assumed to follow the proportional hazards assumption. We consider Identity link and Exp link in the simulations (Wang and Schaubel, 2018).
We compare all these methods under Model 1 and Model 2, and use the Mean-Absolute-Error (MAE) and Rooted-Mean-Squared-Error (RMSE), introduced in Davison and Hinkley (1997), Tian et al. (2007), and Wang and Schaubel (2018), to measure the performance of these methods.
We set n = 3, 000, SNR = 0.3. For Identity link and Log-exp link, λ C = 0.08, L = 5.3. For Exp link λ C = 0.0026, L = 190. We calculate the MAE and RMSE for our method and four existing methods(both Lu.method and Wang.method have two link functions) under Model 1 and Model 2 and p = 5, 10, 20. Among all the considered models, our method in general has a better performance. As an example, Figure 4 visualizes the observed hRMST generated from Log-exp link and predicted hRMST from our method and Wang.method, showing that the random forest can give better predictions. Tables 1, 2 show the MAE and RMSE for Model 1 and Model 2, respectively. For Model 1, the parametric models are correctly specified using the methods of Tian et al. (2014), Wang and Schaubel (2018), we expect that both methods perform well, and our method can have a comparable performance. For Model 2, our proposed method dominates all other methods. Increasing the number of non-predictive covariates does not have a big impact on the performance of our method.
When the censoring distribution does not follow PH assumption, we may expect a difference in the prediction performance because of the bias of IPCW from mis-specification. To check whether our method can still outperform the existing methods, we conduct additional numerical studies. In particular, we simulate the censoring time from the following gamma distributions C ∼ Ŵ(α, β), β = 1 λ C exp(X 1 log 2) , and α ∈ {0.5, 1.5} When α = 1, the gamma distribution degenerates to the exponential distribution we used for Tables 1, 2. Tables 3, 4 show the MAE and RMSE for Model 1 and Model 2 when α = 0.5, and Tables 5, 6 show the MAE and RMSE for Model 1 and Model 2 when α = 1.5. Results of α ∈ {0.5, 1.5} are not very different from the results of α = 1. Under Model 1, our method performs comparably well as methods of Tian et al. (2014), Wang and Schaubel (2018), and it dominates the others under Model 2. When feature dimension is low(p = 5), the error metrics of our method when α = 1 are in general lower than the error metrics when α = 0.5, 1.5 for both Model 1 and Model 2. The additional errors can be regarded as the bias induced from the violation of PH assumption of the censoring distribution. When feature dimension is high(p = 10, 20), bias from large p may dominate the bias from the violation of PH assumption of the censoring distribution.

APPLICATION TO THE TCGA OVARIAN CANCER DATA SET
We apply the proposed method to The Cancer Genome Atlas (TCGA) ovarian cancer functional proteomics data set (Akbani et al., 2015) that is publicly available (http://gdac.broadinstitute. org). The data sets include proteomic characterization of tumors using reverse-phase protein arrays (RPPA). Specifically, Akbani et al. (2015) reported an RPPA-based proteomic analysis using 195 high-quality antibodies that target total, cleaved, acetylated and phosphorylated forms of proteins in 412 high-grade serous ovarian cystadenocarcinoma (OVCA) samples. The function space covered by the antibodies used in the RPPA analysis emcompasses major functional and signaling pathways of relevance to human cancer, including proliferation, DNA damage, polarity, vesicle function, EMT, invasiveness, hormone signaling, apoptosis, metabolism, immunological, and stromal function as well as transmembrane receptors, integrin, TGFβ, LKB1/AMPK, TSC/mTOR, PI3K/Akt, Ras/MAPK, Hippo, Notch, and Wnt/beta-catenin signaling (Akbani et al., 2015).
After removing a few samples with missing data, the final data set includes 407 OVCA samples with a mean/median follow-up of 3.20/2.79 years and a total of 242 deaths and 40% censoring. To assess how different methods predict the hRMST, we performed the following cross-validation analysis. For a given L, we did 10-fold cross-validation on the data set. For each training data set in the cross-validation, we perform a univariate analysis to select top 5 most significant features based on univariate Cox regression analysis. We then estimate the hRMST on the test set using the training data sets with these 5 features as the predictors. We apply 7 different methods, including estimate based on the KM estimator, estimate based on the Cox model, the method of Tian et al. (2014) and the method of Wang and Schaubel (2018). We report the average  Wang and Schaubel (2018) with identity link; Wang:exp, method of Wang and Schaubel (2018) with exponential link.  Tian et al. (2014) with exponential link; Wang.id, method of Wang and Schaubel (2018) with identity link; Wang:exp, method of Wang and Schaubel (2018) with exponential link. of MAE and RMSE on the samples in the testing sets over the 10-fold cross-validation. The results are shown in Table 7 and Figure 5 for L = 3, 4, 5 (see Supplementary Material for L = 6, 7, 8). There are 45.9, 31.2, 19.4, 11.8, 8.1, 4.4% of the observations larger than L for L = 3, 4, 5, 6, 7, 8 correspondingly. For different choices of L, our proposed random forest based method dominates the other methods in MAE and RMSE. The methods of Tian et al. (2014) and Wang and Schaubel (2018) are based on parametric form of hRMST. Cox model is heavily dependent on the proportional hazard assumption, and the Kaplan-Meier approach does not take the covariates into account. We also notice that the method of Wang and Schaubel (2018) always performs better than the method of Tian et al. (2014), possibly due to the fact that the censoring mechanism in the data depends on the covariates.

DISCUSSION
In this paper, we have developed a non-parametric random forest-based method for estimation of hRMST. Compared with traditional Cox model, which gets hRMST estimates by transforming the estimated hazard functions, directly modeling hRMST would be more preferable for computation and feature importance analysis. The proposed estimator can relax the parametric assumptions imposed on the survival time used in Tian et al. (2014) and Wang and Schaubel (2018), and can achieve better prediction performance. We have derived the asymptotic distribution of the random forest estimator using IPCW approach, and presented a procedure based on bags of little bootstraps to obtain the variance of the estimator. Our simulation results and analysis of TCGA data sets have shown promising performance in predicting hRMST as compared to the other available methods, even when the dimension is high and the covariates include irrelevant variables. The method is implemented by R and C++, and is available at https://github. com/lmy1019/SRF.
The proposed method can be used to estimate the heterogeneous treatment effects in randomized clinical trials when the outcome is censored. One can simply apply the method separately to the treated group and the placebo group and take the difference. However, for the observational studies, one needs to account for the fact that the treatment assignments might not be completely at random. Wager and Athey (2015) developed a non-parametric causal forest for estimating heterogeneous treatment effects that extends Breiman's random forest algorithm. In the potential outcomes framework with non-confounding, they showed that causal forest are pointwise consistent for the true treatment effect and have an asymptotically Gaussian and centered sampling distribution. For the observational studies with censored survival outcomes, it is also possible to combine the methods proposed here and the method of Wager and Athey (2015) in order to estimate the treatment effect on the restricted mean survival time.
The proposed methods can also be extended to take into account possible competing risk. This can be done by introducing an additional inverse probability weight (IPCW) to differentiate the non-informative censoring and competing risk censoring. In this case, the estimation equation ψ function with covariates historyX =x under true G C and G R becomes where under competing risk scenario, δ L = 1 {T∧L≤C∧R} . The method proposed in this paper can be automatically adapted to the competing risk case and the asymptotic normality result can be derived similarly.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author/s.