{'fraction_non_alphanumeric': 0.03840977907367949, 'fraction_numerical': 0.01307895031961422, 'mean_word_length': 4.532144241954246, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 17, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ensemble methods for survival function estimation with time-varying covariates

Survival data with time-varying covariates are common in practice. If relevant, they can improve on the estimation of survival function. However, the traditional survival forests - conditional inference forest, relative risk forest and random survival forest - have accommodated only time-invariant covariates. We generalize the conditional inference and relative risk forests to allow time-varying covariates. We also propose a general framework for estimation of a survival function in the presence of time-varying covariates. We compare their performance with that of the Cox model and transformation forest, adapted here to accommodate time-varying covariates, through a comprehensive simulation study in which the Kaplan-Meier estimate serves as a benchmark, and performance is compared using the integrated L2 difference between the true and estimated survival functions. In general, the performance of the two proposed forests substantially improves over the Kaplan-Meier estimate. Taking into account all other factors, under the proportional hazard (PH) setting, the best method is always one of the two proposed forests, while under the non-PH setting, it is the adapted transformation forest. K-fold cross-validation is used as an effective tool to choose between the methods in practice.


Introduction
Methodology for survival data often assumes that covariate information is time-invariant; that is, only values measured at time 0 are used. In this situation, survival analysis models can provide an estimate of the survival function (the probability of surviving past time t) for a subgroup of the population (i.e., a subpopulation) with a specific set of values for the covariates. Time-varying covariates, however, are common in practice and play an important role in the analysis of censored time-to-event data. For example, in a study of the effect of heart transplant on survival for heart patients, the occurrence of a transplant can be modeled as a time-varying binary covariate (Crowley and Hu, 1977), and in a study of the effect of CD4+ T-cell counts on the occurrence of AIDS or death for HIV-infected patients, the cell count was used as a time-varying numerical covariate, measured longitudinally (Tsiatis et al., 1995).
The Cox proportional hazards model (Cox, 1972) has a long history of being used to model and analyze censored survival data. As a semi-parametric model, it assumes that the time-invariant covariates have a proportional effect on the hazard function. The Cox model was extended to fit time-varying covariates using a counting process formulation as follows (Andersen and Gill, 1982). Consider continuous-time survival data with time-varying covariates, where each subject may have multiple records of measurements of risk factors at multiple time points. In practice, as the subjects are observed intermittently, the time-varying covariates are assumed constant between observation times. One can then reformat the data structure using the counting process approach by which a data record of a subject becomes a list of pseudo-subjects, that are treated as being independent, left-truncated and right-censored observations.
It is important to recognize that survival regression models, like all regression models, can be used for the two distinct purposes of estimation and prediction. These two purposes, while related (prediction almost always involves estimation of some kind as a first step), are distinctly different from each other. Prediction is meaningful at the level of an individual. In the context of survival data, this would correspond to an estimated survival function for a particular individual. For such an individual, if a change in the value of a covariate at time t * can potentially impact the future probabilities of survival of that individual, then the estimated survival function for that individual must be 1 for any time t < t * (since the covariate couldn't impact the future probabilities of survival unless the individual was alive at that time). A related point is the existence of so-called internal covariates, in which a variable (for example, blood pressure) can only be measured when an individual is alive, meaning that the act of measuring the covariate implies that the survival function value for that individual must equal 1. These facts, and their implications for prediction, are well-known, and are the reason for the prominence of joint modeling methods (Rizopoulos, 2012) in the prediction of survival data with time-varying covariates.
The focus of this paper, however, is exclusively on estimation, rather than prediction. The goal here is not prediction of an individual's survival function, but rather estimation of a population-level survival function. This is an important problem in any situation that involves strategic decision making, such as public health policy or business operations and planning. Civic policymakers are interested in developing strategies that change over time, in response to changing conditions; this does not involve prediction at the individual level, but rather estimation of population probabilities.
Such estimation can answer "what if" questions, in order to either be prepared for what might happen, or to try to control what will happen. Recognizing that various aspects of society can change, it can also examine what the effects of such changes might be. Examples of this could include modeling the incidence of COVID-19 infection in the population as vaccination and testing rates change, or estimating what proportion of deliveries will be made by Christmas as various characteristics of the supply chain change. These problems are fundamentally different from the problem of making a prediction for a particular individual, and issues related to joint modeling and the distinction between internal and external covariates are no longer relevant.
When time-varying covariates are available, it is important to use the updated covariate values to dynamically adjust the estimated survival function. We use a simple hypothetical example to fix ideas and illustrate the mechanics of the survival function adjustment. Consider the following simple example, based on the COVID-19 problem noted earlier.
The time-to-event is the time to a positive COVID test, and there is a time-varying covariate X(t) that describes the vaccination status of a subject, assuming a two-shot vaccine (unvaccinated, vaccinated with one dose, vaccinated with two doses, vaccinated with two doses and a booster). Starting with a sample of unvaccinated and COVID-free subjects, the changes of X(t) define distinct subgroups of the sample (and hence the population). Estimation of the population survival function for subjects who have gotten one dose of the vaccine would be based only on those subjects who survived to get that dose, but the overall population survival function estimate would have to be adjusted to account for the probability that unvaccinated subjects did not survive to that time. This same adjustment would be applied at the other times that vaccination status changes, resulting in a survival function estimate that is appropriate for all members of the population who followed the specified temporal vaccination pattern. In the next section we describe how such a survival function estimate can be constructed, and in Appendix A we illustrate the detailed calculations of the dynamically adjusted survival function for this simple example.
The Cox proportional hazards model with time-varying covariates, hereafter referred to as the extended Cox model, relies on restrictive assumptions such as proportional hazards and a log-linear relationship between the hazard function and covariates. Tree-based methods and their ensembles, which are useful non-parametric alternatives to the extended Cox model, also can incorporate time-varying covariates. Recently, two types of survival trees were proposed as extensions of the relative risk tree (LeBlanc and Crowley, 1992) and of the conditional inference tree , respectively, to left-truncated right-censored (LTRC) data, referred to as LTRC trees (Fu and Simonoff, 2017). The proposed LTRC tree algorithms allow for time-varying covariate data after the data structure is reformatted using the counting process approach. Another tree-based method that can handle LTRC survival data and therefore potentially be applied to time-varying covariate data is the novel "transformation tree," and the corresponding ensemble is the "transformation forest" (Hothorn and Zeileis, 2021). These two algorithms are based on a parametric family of distributions characterized by their transformation function and developed to detect distributional alternatives to proportional hazards. None of the above methods have considered the estimation of the survival function. Similarly, recently developed methods for hazard function estimation in the presence of time-varying covariates haven't dealt with survival function estimation in general (Sun et al., 2020;Wongvibulsin et al., 2020). There exist other survival trees and forest methods that can handle time-varying covariate data, but only for discrete-time survival data (Bou-Hamad et al., 2011;Schmid et al., 2020;Kretowska, 2020;Puth et al., 2020;Moradian et al., 2022).
In this paper, we focus on forest algorithms for dynamic estimation of the survival function for continuous-time survival data. Ensemble methods like forest algorithms are known to preserve low bias while reducing variance and therefore can substantially improve prediction accuracy, compared to tree algorithms (Breiman, 2001). The most well-known ensemble methods for survival analysis are perhaps the relative risk forest (Ishwaran et al., 2004), random survival forest (Ishwaran et al., 2008) and conditional inference forest . These forest methods provide estimates of survival functions, but only for right-censored survival data with time-invariant covariates. We propose to extend the relative risk and conditional inference forests, as well as the transformation forest, to allow time-varying (TV) covariates. We refer to them as RRF-TV, CIF-TV, and TSF-TV, respectively.
The proposed methods by design can handle survival data with all combinations of left-truncation and right-censoring in the survival outcome, and with both time-invariant and time-varying covariates. In this paper we focus on survival data with time-varying covariates. Similar analysis for LTRC data with time-invariant covariates is provided in Section S2 in the Supplemental Material.

Proposed forests for time-varying covariate data
Assume p covariates X = (X 1 , X 2 , . . . , X p ) are available, some of which are TV covariates and the others are time-invariant (TI). For example, assume X 1 is the only time-invariant covariate among all p covariates, then at time t, X(t) = (X 1 , X 2 (t), . . . , X p (t)). For ease of exposition, we write X(t) = (X 1 (t), X 2 (t), . . . , X p (t)) with X 1 (t) ≡ X 1 for all t. Observations are obtained from N subjects. Note that the subjects are observed only intermittently, for example, J (i) times for subject i, initially observed at t (i) 0 , and then at t (i) j , j = 1, 2, . . . , J (i) − 1, for followup visits, with corresponding observed values x , the minimum value of the true survival time T (i) and censoring time C (i) . We assume non-informative censoring conditional on the covariates. Denote ∆ = I{T (i) ≤ C (i) }, which indicates whether a subject experienced an event (∆ = 1) or was right-censored (∆ = 0). If t (i) 0 = 0, then we say the survival time is left-truncated. The outcome of interest is the time to the event.
At any time t, let I(t) denote an arbitrary set of time values up to time t, that is, I(t) ⊆ [0, t]. It could be a finite number of time points, a finite number of intervals, or a disjoint set of time intervals and/or time points. Given the historical data for N subjects observed up to the death or censoring time, the goal is to estimate the conditional survival function S(t | X(u) = x(u), u ∈ I(t)) where {x(u), u ∈ I(t)} is a particular stream of covariate values. This true survival function is defined as the population proportion of subjects who have the specified covariate values at the specified times up to either time t or their event time (whichever comes first) that are alive at time t. Note that the population includes all subjects for whom the specified conditions hold up until their time of event if that occurs before the evaluation time t, even if the conditions do not hold after the time of event. This is true if measurement of the covariate is no longer meaningful after the event occurs (as might be the case for a so-called internal covariate, such as blood pressure), or it is meaningful and available but no longer satisfies the conditions (a so-called external covariate, such as pollutant level). The reason is that the influence of the covariate on survival in either case is irrelevant for a subject for whom the event has occurred.
The proposed forest methods provide the survival function estimate by following three steps. First, we adopt the counting process approach to reformat the data structure. This approach assumes that the time-varying covariates are constant between the observed time points, that is, It then splits the i-th subject observation into J (i) pseudo-subject observations: (t , and event indicator δ (i) j = ∆I{j = J (i) − 1}, j = 0, 1, . . . , J (i) − 1. The multiple records from N subjects now become a list of pseudo-subjects, The set of pseudo-subjects is treated as if they were independent in the following form where x l = (x l,1 , . . . , x l,p ) is the vector of the observed values of p covariates from the l-th pseudo-subjects in the reformatted dataset. The second step is to apply the forest algorithms on the reformatted dataset given in (1), to fit a model. Finally, in the third step, given a particular stream of covariate values x * j at the corresponding time values t * j , j = 0, 1, . . ., a survival function estimate is constructed based on the outputs of the proposed forest algorithms. More specifically, at any time t, with X * (t) denoting the covariate information up to time t, we compute the estimated survival probability S(t | X * (t)).

Extending right-censored TI survival forests to the proposed TV forests
The conditional inference forest and the relative risk forest are both tree-based ensemble methods, where B individual trees are grown from B bootstrap samples drawn from the original data. Randomness is induced into each node of each individual tree when selecting a variable to split on. Only a random subset I of the total p covariates is considered for splitting at each node. The node is then split using the candidate covariates based on different criteria for different forest methods. To extend the two forest methods for right-censored survival data with time-invariant covariates to the forests for (left-truncated) right-censored survival data with time-varying covariates, the splitting criteria are modified.

Recursive partitioning in the proposed CIF-TV forest
Consider right-censored survival time data of the form T , ∆, X , with survival/censored timeT , event indicator ∆ (T denotes the survival time if ∆ = 1, or censored time if ∆ = 0), and p time-invariant covariates X = X 1 , . . . , X p . In each node, the recursive partitioning in the conditional inference forest algorithm is based on a test of the global null hypothesis of independence between the response variable in the right censored case V = T , ∆ and any of the covariates in the random subset I. It is formulated in terms of |I| partial hypotheses, where D(V | X k ) denotes the conditional distribution of V given the covariate X k . The independence is measured by linear statistics incorporating the log-rank scores that take censoring into account. In the extension of conditional inference tree to LTRC conditional inference tree, the log-rank score can be modified as follows for LTRC data (Fu and Simonoff, 2017).
Given the list of pseudo-subject observations with LTRC survival times as in (1), the response variable now becomes V = L l , R l , δ l in the test of partial null hypothesis of independence (3) for the l-th observation L l , R l , δ l , x l . The corresponding log-rank score is defined as Note that S is the nonparametric maximum likelihood estimator (NPMLE) of the survival function which takes into account left-truncation. Such an estimator can be constructed using the product-limit estimator, i.e. Kaplan-Meier estimator with pseudo-subjects that fall into the current node (Gross and Lai, 1996;Tsai et al., 1987). We similarly use the log-rank score U l in the proposed extension of conditional inference forest to LTRC conditional inference forest.

Recursive partitioning in the proposed RRF-TV forest
The relative risk forest combines the use of relative risk trees (LeBlanc and Crowley, 1992) with random forest methodology (Breiman, 2001) as a way to reliably estimate relative risk values. The Classification and Regression Tree (CART) paradigm (Breiman et al., 1984) is used to produce a relative risk forest by exploiting an equivalence with Poisson tree likelihoods.
The splitting criterion under the relative risk framework is to maximize the reduction in the one-step deviance between the log-likelihood of the saturated model and the maximized log-likelihood. At node h, let R h denote the set of labels of those observations that fall into the region corresponding to node h, and let λ h (t) and Λ h denote the corresponding hazard and cumulative hazard function, respectively. Under the assumption of proportional hazards, where λ 0 is the baseline hazard and ϕ h is the nonnegative relative risk of the node h. Given the right-censored observations t l , δ l , l ∈ R h , the maximum likelihood estimate of ϕ h is , where the Nelson-Aalen estimator using all of the data at the root node Λ 0 is used for Λ 0 (LeBlanc and Crowley, 1992). The full likelihood deviance residual for node h is defined as For a Poisson regression model, let h denote the event rate, s l and c l be the exposure time and the event count for observation l, respectively, then (5) is equivalent in form to the deviance residual based on the Poisson regression model, with h = l∈R h c l l∈R h s l , by replacing h with ϕ h , s l with Λ 0 (t l ), and c l with δ l (LeBlanc and Crowley, 1992).
To adapt the Poisson regression tree approach for left-truncated right-censored survival observations {(L l , R l , δ l )}, the key is to modify the estimated Λ 0 (t l ) and δ l to replace s l , c l and h in (6). First, compute the estimated cumulative hazard function Λ 0 (·) based on all (pseudo-subject) observations. The exposure time s l and the event count c l for observation l in (6) are then replaced by Λ 0 (R l ) − Λ 0 (L l ) and δ l to obtain the deviance residual appropriate for LTRC data (Fu and Simonoff, 2017).

Implementation of the proposed forests
To implement the CIF-TV and RRF-TV algorithms, we make use of the fast algorithms provided in the packages partykit  and randomForestSRC (Ishwaran and Kogalur, 2020), respectively. The RRF-TV building architecture is based on employing the fast C code from randomForestSRC. The Poisson splitting rule (LeBlanc and Crowley, 1992) is coded in C and is incorporated by exploiting the custom splitting rule feature in the rfsrc function. The CIF-TV is built by extending the survival forest algorithms coded in the cforest function from partykit with the log-rank score adapted for LTRC data.

Bootstrapping subjects vs. bootstrapping pseudo-subjects
In forest-like algorithms, bootstrapped samples are typically used to construct each individual tree to increase independence between these base learners. The nonparametric bootstrap approach is used in all three types of forests being considered here (CIF-TV, RRF-TV, and TSF-TV). It places positive integer weights that sum to the sample size on approximately 63% of the observations in any given bootstrap sample, and the 37% of the data excluded during this procedure is called out-of-bag data (OOB data). As we split each subject into several pseudo-subjects and treat these pseudo-subjects as independent observations on which to build the forests, we have two bootstrapping options: we can bootstrap subjects or bootstrap pseudo-subjects.
Bootstrapping pseudo-subjects is used for some discrete survival forest methods (Bou-Hamad et al., 2011;Schmid et al., 2020). Since all pseudo-subjects are treated as independent observations in the recursive partitioning process (Bacchetti and Segal, 1995;Fu and Simonoff, 2017), bootstrapping pseudo-subjects is just bootstrapping "independent" observations as the first step of any forest algorithm. On the other hand, bootstrapping subjects is a natural approach, as it keeps all of the pseudo-subjects for each subject in the bootstrap sample. In fact, simulations have shown that the two different bootstrapping mechanisms do not result in fundamentally different levels of performance; see Section S1.5 in the Supplemental Material for more details. This paper will focus on forests based on bootstrapping subjects.

Regulating the construction of individual trees in the proposed forests
In a forest algorithm, only a random subset of covariates is considered for splitting at each node. The size of this random set is denoted by mtry. In addition to mtry, many other parameters play an important role in establishing a split in the individual tree. In both the cforest  function for the conditional inference forest and the traforest(Hothorn, 2020) function for the transformation forest algorithms, minsplit (the minimum sum of weights in a node in order to be considered for splitting), minprob (the minimum proportion of observations needed to establish a terminal node) and minbucket (the minimum sum of weights in a terminal node) control whether or not to implement a split; in the rfsrc (Breiman et al., 2018) function for the relative risk forest algorithms, nodesize controls the average terminal node size. These tuning parameters thereby regulate the size of the individual trees. The recommended values for these parameters are usually given as defaults to the algorithm. For example, mtry is usually set to be √ p, where p is the total number of covariates Ishwaran et al., 2008), nodesize to be 15 in the rfsrc function, (minsplit, minbucket) to be (20, 7) in the cforest function and the traforest function, which we refer to as the default parameter settings.
The best values for these parameters would be expected to depend on the problem and they should be treated as tuning parameters (Hastie et al., 2001). It has been shown for conditional inference forests for interval-censored data (Yao et al., 2021) that these parameters have a non-negligible effect on the overall performance of the forest algorithm. As we extend the forest framework to allow for left truncation, and from time-invariant covariate data to time-varying covariate data, we should also consider rules for choosing tuning parameters.
The algorithm designed in survival forests for interval-censored data with time-invariant covariates (Yao et al., 2021) tunes the value of mtry on the "out-of-bag observations". To adapt the same idea to survival forests based on bootstrapping subjects on a dataset with time-varying covariates, we define the "out-of-bag observations" for the b-th tree to be the observations from those subjects that are left out of the b-th bootstrap sample and not used in the construction of the b-th tree. The survival curve can be estimated by using each of the B trees in which that subject was "out-of-bag," denoted as S OOB . To evaluate the fit of the out-of-bag estimate S OOB with a specific value of mtry, we compute the estimation error defined as the integrated Brier score designed for time-invariant covariate data (Graf et al., 1999), adapted here for time-varying covariate data as follows.
For a given dataset D, define the integrated Brier score IBS S; D for the estimated survival function S as where τ (i) determines the length of difference evaluation time span for subject i,Ỹ (i) (t) = I{T (i) > t} is the observed status (T (i) is the survival/censored time), W (i) (t) is the inverse probability of censoring weights, with G the Kaplan-Meier estimate of the censoring distribution based on {(T (i) , 1 − ∆ (i) )} i∈D (Graf et al., 1999). The corresponding Brier score BS(t, S; D) at time t is defined as The resulting estimation error for the ensemble method with a specific value of mtry can then be computed by setting S = S OOB in (7). An appropriate value of mtry is the one that minimizes the "out-of-bag" estimation error.
Regarding the values of other tuning parameters, the optimal values that determine the split vary from case to case. As fixed numbers, the default values may not affect the splitting at all when the sample size is large, while having a noticeable effect in smaller data sets. This inconsistency can potentially result in good performance in some data sets and poor performance in others. In the simulations, we set minsplit, minbucket and nodesize to be the maximum of the default value and the square root of the number of pseudo-subject observations n. This set of values can automatically adjust to the change in size of the data set. We refer to the above choice of mtry, minsplit, minbucket and nodesize as the proposed parameter settings, as opposed to the default settings.

Constructing a survival function estimate for time-varying covariate data
Consider a particular stream of covariate values x * j at time t * j , for j = 0, 1, . . . , J − 1. Denote X * (u) = {x * j : ∀j, t * j ≤ u} the set of the covariate values up to time u. At time t ∈ t * j , t * j+1 , we derive a recursive computation as follows where S A,j (t) P(T > t | x * j ) denotes the output of the algorithm for the input with time-invariant covariate value x * j . See Appendix B for details of derivation. Further, by expansion, (9) is equivalent to for t ∈ t * j , t * j+1 . The formulation in (10) provides another perspective to view the resulting survival curve estimate -it is constructed by combining the pseudo-subject-specific ensemble estimates of the survival function with multiplicative correction factors. These correction factors ensure monotonicity of the overall curve.
Note that the construction in (10) coincides with what the function survfit in the R package survival (Therneau et al., 2020) uses to give a subject's survival function estimate from a coxph fit using the same counting process approach. Figure 1 gives an illustration of the estimated survival functions with or without "updating" the covariate values at time t * 1 , t * 2 , . . . , t * 5 . For each j, the "update" in the estimated survival probability starting from time t * j for all the future time t > t * j reflects the difference in the estimated surviving proportions of two subpopulations of subjects with their covariate trajectories diverging from the shared past before t * j . One can see that the change in covariate information at each time point of change can make a huge impact on the future path of the estimated survival function. In the simulation study, observations from N subjects are generated independently with p covariates X = X 1 , . . . , X p . We set p = 20. Eight of these covariates are time-invariant: X 1 , X 11 ∼ Bern(0.5), X 2 , X 7 , X 10 ∼ Unif(0, 1), X 8 ∼ Unif(1, 2), X 9 follows a categorical distribution with possible outcomes {1, 2, 3, 4, 5} with equal probability, X 12 follows a categorical distribution with possible outcomes {0, 1, 2} with equal probability. The others are timevarying, whose values are obtained at m randomly generated time points, different for each subject. In the simulations, we set m = 11. At each of these preset time points, for some time-varying covariates, the value is randomly resimulated from its distribution: X 3 , X 19 ∼ Bern(0.5), X 4 , X 15 , X 17 ∼ Unif(0, 1), X 5 and X 14 both follow a categorical distribution with possible outcomes {1, 2, 3, 4, 5} with equal probability; for other time-varying covariates, the value is resimulated following particular patterns: X 6 , whose initial value is randomly generated from {0, 1, 2} with equal probability, which will choose to stay at the original value or move one level up but the largest value can only be 2; the changing pattern of X 13 is always 0 → 1; the changing pattern of X 16 is either 0 → 1 or 1 → 2; the changing pattern of X 18 is 0 → 1 → 2; value of X 20 is a linear function of the left-truncated time point of the interval with slope and intercept follows Unif(0, 1). Further details of the changing pattern of X 6 , X 13 , X 16 , X 18 and X 20 can be found in Section S1.2 in the Supplemental Material.
After the time-varying covariates' values are generated at each of those m preset time points, the true survival time T is then computed under different model setups and the right-censoring time C is generated independently.

Model setup
We consider the following factors for different variations of data generating models: Scenario We consider two different proportions of time-varying covariates in the true model: 2TI + 1TV, and 2TI + 4TV; see Table 1. Only the first six covariates are given in the table since X 7 to X 20 are never involved in the true DGP.
Survival relationships Given a survival relationship, the survival time T depends on ϑ(t) = ϑ(X(t)). Here we use scenario 2TI + 4TV for illustration. For a linear survival relationship, with some constants φ 1 , φ 2 , φ 3 , and ψ k , k = 0, . . . , 6. For an interaction model, ϑ is determined by the value of time-varying covariate X 4 and the value of time-varying covariate X 5 . Figure 2 gives an example of the structure of the covariates driving the interaction survival relationship, where T 1 , T 2 , T 3 , T 4 correspond to Survival distributions under the PH and the non-PH setting Given a survival relationship model, the survival time T depends on ϑ via a Weibull distribution.
For proportional hazards models, a closed-form solution can be derived to generate survival times with time-varying covariates for the Weibull distribution (Austin, 2012). For non-proportional hazards models, a closed-form solution exists for the Weibull distribution, with its nonconstant shape term a function of the covariates (note that the proportional hazards relationship is on the scale parameter for the Weibull distribution).
To be more specific, for the proportional hazards setting, we consider the underlying hazard function where the baseline hazard function is given by h 0 (t) = λνt ν−1 with λ > 0 and ν > 0. For the non-proportional hazards setting, the hazard function is set to be where λ > 0. Values of ϑ(t) have been scaled to be between −3 and 3. Note that, compared with the Weibull distribution under the PH setting, now the time-varying effects appear in the shape term instead of the scale term. The survival function is then given by S(t) = exp(− t 0 h 0 (s) exp(ϑ(s))ds). Further details of simulating the survival time T can be found in Section S1.1 in the Supplemental Material.
Histograms of survival times for typical samples with the number of subjects N = 500 in each scenario are provided in Section S1.3 in the Supplemental Material to illustrate the data generating processes. The parameters set in the simulation study can be found in Section S1.4 in the Supplemental Material.
Knowledge of history of changes in covariates' values In practice, it is likely that not all of the changes in the covariates' values are known to the data analyst. For example, suppose that a patient's blood pressure is to be measured at regularly scheduled examination times. If a patient obeys the schedule then, from the doctor's point of view, all changes in blood pressure are known. However, if a patient skips some scheduled examination times, then not all changes in the blood pressure are known. In the latter case, this means that whatever modeling method is used to estimate the survival curves, it is operating with incorrect values as inputs and therefore its performance would be expected to deteriorate. Of course, the fact that blood pressure is actually changing continuously is an extreme example of this phenomenon; in these simulations we limit ourselves to changes at a finite number of time points. The simulations are designed to investigate the performance of different modeling methods in this situation under the following two circumstances: -Case I. When all changes in covariate values are known; -Case II. When only half of the changes in covariate values are known.
The missing changes are selected completely at random. To generate a dataset under Case II, one can start with the dataset generated under Case I. The following example is given to illustrate how to construct such datasets. Suppose the baseline covariates' values of the subject is X(t 0 ) = x 0 and the covariates values X change J − 1 times at time t 1 , . . . , t J−1 , before the subject is censored or the event occurs at t J =T . For J = 3, X(t 1 ) = x 1 and X(t 2 ) = x 2 . The counting process approach assumes The information of the subject under Case I, displays exactly as in (13). For a dataset under Case II when only half of the changes are known, only one of {t 1 , t 2 } is known. If only the change at t k (k = 1, 2) is known, the observed information for the same subject is then Note that for both (13) and (14), the true survival curve is constructed using the information as in (13), when all history of changes in values are known.

Evaluation measures
Since the goal is to estimate the survival function, we evaluate estimation performance using the average integrated L 2 difference between the true and the estimated survival curves S. Given a dataset D, containing N subjects, each with pseudo-subject information up to the survival/censored timeT (i) , X (i) (T (i) ), i = 1, 2, . . . , N , Note that we evaluate the integrated L 2 difference only up toT (i) , the last time point where the survival status is known.
In the simulations, as we generate the true survival time T (i) , we have the trajectory of covariate values up to time T (i) for any given subject i even when it is censored at time C (i) < T (i) . However, here we intend to match the scenario in real world applications where the covariate information is usually no longer recorded after the event occurs (e.g. the patient dies) or the subject is censored (e.g. lost contact). Thus, we define the best modeling method to be the one that gives us the lowest integrated L 2 difference, which is an average value from all subjects; for each subject, the difference between an estimated survival curve and the true survival curve up to its last observed time is measured.

Simulation results
The extended Cox model is included as a benchmark method, since it is one of the most commonly used methods in practice. Another benchmark method used in this paper is the Kaplan-Meier method, which uses no covariates' values given results for the linear, nonlinear and interaction survival relationship; the first three columns show results under the PH setting for the number of subjects N = 100, 300, 500, respectively, and the last three columns for results under the non-PH setting. In each plot, 1-CIF-TV with mtry = 1; 2-CIF-TV with mtry = 2; 3-CIF-TV with mtry = 3; 5-CIF-TV with mtry = 5; 10-CIF-TV with mtry = 10; 20-CIF-TV with mtry = 20; Opt-CIF-TV with value of mtry that gives the smallest Integrated L 2 difference in each round; Tuned-CIF-TV with the value of mtry tuned by the "out-of-bag" tuning procedure. The default value in conditional inference forest is mtry = 5.
to construct the survival function estimate; this helps illustrate the improvement in estimation from incorporating the covariate information.
In this section, we present simulation results based on 500 simulation trials. The number of trees for bootstrap samples is set to be 100 for all forest methods. We also only focus on the Weibull-Increasing distribution, and omit discussion of the Weibull-Decreasing distribution, since results for the latter distribution are similar. Detailed results are given in Section S1.7 of the Supplemental Material. Figure 3 gives an example of how CIF-TV performs with different values of mtry in the scenario 2TI + 4TV, when the censoring rate is 20%, and the signal-to-noise ratio is low. The mtry values are tuned based on the "out-of-bag observations". Similar results for RRF-TV and TSF-TV can be found in Section S1.6 in the Supplemental Material.

Regulating the construction of trees in forests
In these examples, one can see that the forests using the "out-of-bag" tuning procedure give relatively good performance overall. In fact, results from other model setups are broadly similar, in the sense that this tuning procedure provides a relatively reliable choice of mtry and it gains in accuracy as the number of subjects N increases. In contrast, the default value of mtry does not always perform well, and choosing a different value can sometimes significantly improve performance. Table 2: Comparison between forests with default (D) and proposed parameter settings (P) across different numbers of subjects N in the scenario 2TI + 4TV, when the censoring rate is 20%, and the signal-to-noise ratio is low. based on all simulations. * For similar results under other model setups, please refer to Section S1.6 in the Supplemental Material. Table 2 gives examples under the PH setting to show the performance comparison between each forest with its default parameter settings and with the proposed parameter settings.
In Table 2, positive numbers indicate a decrease in integrated L 2 difference compared to a Kaplan-Meier fit on the dataset, while negative numbers indicate an increase. The absolute value of the numbers represents the size of the difference between the integrated L 2 difference of the candidate and that of a Kaplan-Meier fit. The table shows that forests with the proposed parameter settings can provide improved performance over those with default parameter settings across all different numbers of subjects N by a substantial amount. Note that, under the non-PH setting, for datasets with all of the changes in covariates' history known, the negative numbers indicate the poor performance of forests with default parameter settings even compared to a simple Kaplan-Meier curve, showing that the default methods can fail miserably. In contrast, for all forests with the proposed parameter settings, as N increases, the change in sign and in the absolute value of the numbers indicates better and better performance in general. Overall, the performance of the proposed parameter setting is relatively stable and better than that of the default values.
In the following discussion, we therefore only focus on the forest methods with the proposed parameter settings.

Properties of the proposed forest methods
Using factorial designs, we study the difference between each of the proposed forest methods and a simple Kaplan-Meier fit under the effects of the following factors: censoring rate, amount of knowledge, survival relationship, training sample size, scenario, hazards setting and SNR. The effects are estimated based on an analysis of variance model fit with these factors as main effects. Figure 4 provides main effects plots for the integrated L 2 difference improvement from the proposed forest methods over a simple Kaplan-Meier fit.
In Figure 4, the overall center of location is negative, highlighting that both of the proposed forest methods perform better than a simple Kaplan-Meier fit. The overall mean integrated L 2 difference of CIF-TV is slightly smaller than that of RRF-TV. The relative performance of the proposed forest methods can vary with changes in factors. Note that the p-values of the hypothesis testing on the simple main effect of SNR is insignificant at a 0.10 level, which suggests that the impact of the change of its level on the performance of the proposed forest methods over a Kaplan-Meier fit is negligible (more details can be found in Section S1.8 in the Supplemental Material).
For the other factors, the fewer the number of changes in values of covariates that are known, higher censoring rate, smaller training sample size, larger portion of covariates being time-varying, more relaxed assumption on the hazard setting, and more complicated structure of the survival relationship (all reflecting more difficult estimation tasks), the less the proposed forest methods improve over a simple Kaplan-Meier fit. Conversely, in the opposite situations, the stronger the ability of the proposed forest methods to estimate the underlying survival relationship and therefore bring a greater improvement.
In particular, the two proposed forest methods win by a much larger margin under the proportional hazards setting, which is expected since the log-rank-type splitting procedures used in the proposed forest methods rely to some extent on the proportional hazards assumption.
It is also clear that the difference between the number of time-invariant (TI) and the number of time-varying (TV) covariates is driving the scenario effect. When #TV − #TI increases, the relative performance of the proposed forest methods deteriorates. Presumably, this is because the increasing level of local time-varying effects makes the underlying relationship more difficult to estimate.
Note that the improvement from the proposed forest methods over a Kaplan-Meier fit remains relatively stable to the change of levels in both the censoring rate and the scenario factors compared to the change in other factors. In the following discussion, we mainly focus on the factors that are more influential based on our previous study: the number of changes in values of covariates that are known, the underlying survival relationship, the sample size, and the hazard function setting. The simulation results presented are based on the datasets generated under the scenario 2TI + 4TV, the lower signal-to-noise ratio, with 20% censoring rate.  Datasets are generated with survival times following a Weibull distribution, light right-censoring rate (20%). The three rows show results for the number of subjects N = 100, 300, 500, respectively; the first three columns show results under the PH setting for survival relationship linear, nonlinear and interaction, respectively, and the last three columns for results under the non-PH setting. The horizontal dashed line shows the median integrated L 2 difference of a Kaplan-Meier fit on the datasets. In each of the plots, the set of boxplots lightly-shaded shows the performance of different methods on datasets with history of changes in covariates' values known; the set heavily-shaded shows the performance on datasets with half of the changes in covariates' values unknown.

Estimation performance comparison
performs the best. This is expected as the extended Cox model relies exactly on the assumption of proportional hazards and a log-linear relationship between the hazard function and covariates. For nonlinear and interaction survival relationships, all forests outperform the extended Cox model, showing their advantage in dealing with a relatively complex situation. More specifically, for nonlinear setups, CIF-TV performs the best and RRF-TV the second, while for interaction model setups, RRF-TV performs the best and CIF-TV the second. In addition, CIF-TV and RRF-TV outperform TSF-TV across all different numbers of subjects and survival relationships. Under the non-PH setting, the extended Cox model cannot even outperform a simple Kaplan-Meier fit on the dataset, whether all changes in values of covariates are known or not. As discussed, the presence of non-proportional hazards settings poses great challenges to modeling methods that assume proportional hazards; not just Cox, but also the survival forests like CIF-TV and RRF-TV that use a log-rank splitting rule. On the other hand, TSF-TV, which is specifically designed to detect non-proportional hazards deviations, performs the best across all different setups under the non-PH setting.
It is not surprising that having all changes in values known gives increasingly better performance compared to only having half of the changes known as the sample size increases. As N increases, false information due to unknown changes has a negative effect on performance of all modeling methods. In particular, this affects the extended Cox model more than the forest methods when the underlying survival relationship is linear under the PH setting, while it affects the extended Cox model less for all other cases. This is simply because the extended Cox model already performs poorly in nonlinear and non-PH situations, so the misleading information from incorrect knowledge of covariates' values cannot hurt performance very much.
Generally, if the true underlying model setup is known, one should choose CIF-TV or RRF-TV under the PH setting, and TSF-TV under the non-PH setting. However, none of the forest methods can perform well all of the time. In the next section, we provide guidance on how to choose among these forest methods.

Guidance for choosing the modeling method
Cross-validation methods have been used in the past for the error estimation of survival models (Gerds and Schumacher, 2007). We propose to use one of the most common methods, K-fold cross-validation, implemented with integrated Brier scores for survival data, to select the "best" modeling method for a given dataset D as follows.
For a given survival curve estimate S, i. Split the dataset into K non-overlapping subsets D k (k = 1, 2, . . . , K), each containing (roughly) equal number of subjects; ii. For each k = 1, 2, . . . , K (a) Modeling methods S k are then trained with the data D\D k where the k-th subset is removed; (b) Test S k on data in the k-th test set D k and compute the corresponding integrated Brier score IBS S k ; D k as given in (7); iii. Average over all K subsets and obtain We then choose the modeling method that gives the smallest IBSCVErr( S) in (16).
For the simulated data sets, we use 10-fold cross-validation to choose between modeling methods. The measures p B , r B and r W are used to evaluate the performance, where n rep denotes the number of simulations (n rep = 500), x CV denotes the integrated L 2 difference of the method chosen by cross-validation, and x min and x max denote the lowest and highest integrated L 2 difference among all modeling methods, respectively. In each round of simulation, we call the method that gives x min the best modeling method and the method that gives x max the worst modeling method. By definition, p B provides the proportion of the times IBS-based 10-fold CV selects the best modeling method, and r B and r W compute the relative errors from the best and the worst modeling method, respectively. The smaller r B is, or the larger r W is, the better IBS-based 10-fold CV works. Table 3 presents the summary of the performance of the IBS-based 10-fold CV Rule. It is not surprising that IBS-based 10-fold CV works better under Case I where all changes in covariate values are known in general, with larger values of p B , smaller values of r B and larger values of r W . That is, the incorrect knowledge of covariate values also hurts the performance of the selection procedure. In general, as the number of subjects N increases, the value of p B gets larger for most of the scenarios, indicating IBS-based 10-fold CV is able to pick up the best modeling method at a higher frequency; even under those scenarios where p B is lower than 50%, the relative error from the best modeling method, r B , remains within 10% for most of the cases. Note that more than half of the cases for N = 100 have the relative error from the best modeling method, r B , less than 10%, and almost all of the cases for N = 500 have r B under 5%. That means that even when the IBS-based 10-fold CV does not pick the best modeling method, it is still able to pick a method that works reasonably well, resulting in the integrated L 2 difference being not much higher than that of the best method.  (17), with mean value over all simulations.
2 rB and rW are computed as in (18) and (19), respectively, with mean value ± one standard deviation over all simulations.

Proposed forest methods for time-invariant covariate data
We have focused on ensemble methods for survival data with time-varying covariates, as we feel that this is a very common and important situation that has been understudied in the past. Having said that, there are certainly many situations in which only time-invariant (baseline) covariate information is available, and understanding the properties of different methods in that situation is important. Section S2 in the Supplemental Material describes the results of simulations related to this question. In those simulations, datasets with left-truncated and right-censored survival times are generated based on time-invariant covariates.
In fact, the simulation results of all comparative estimation performance in the case of time-invariant covariates are broadly similar to those in the time-varying covariates cases. That is, 1. The "out-of-bag" tuning procedure can provide a reliable choice of mtry that gives relatively good performance in general. One should also consider adjusting other tuning parameters such as minsplit, minbucket in the conditional inference forest and the transformation forest, and nodesize in random survival forest, as the size of dataset grows; 2. Taking into an account all other factors, under the PH setting, the best method is always one of the two proposed forests, while under the non-PH setting, it is the transformation forest method. 3. The IBS-based CV rule is a good option for choosing among the various methods, as the comparative performance of methods appears to be different from setting to setting.

Real data application
We now illustrate application of the proposed time-varying covariates forests to a real data set, the Mayo Clinic Primary Biliary Cirrhosis Data, available in the R package survival. To study the effectiveness of using D-penicillamine as treatment, 312 patients with primary biliary cirrhosis (PBC) were enrolled in a randomized medical trial at the Mayo Clinic from January in 1974 to May in 1984 . The outcome of interest is the time to death for these patients. In this study, medical measurements and other patient information were recorded as covariates' values at entry and at yearly intervals. The Cox model was used to estimate the survival function for patients with primary biliary cirrhosis based on twelve noninvasive, easily collected covariates that require only a blood sample and clinical evaluation ). These twelve covariates include age at entry, alkaline phosphotase (U/liter), logarithm of serum albumin (g/dl), presence of ascites, aspartate aminotransferase (U/ml), logarithm of serum bilirubin (mg/dl), serum cholesterol (mg/dl), condition of edema, presence of hepatomegaly or enlarged liver, platelet count, logarithm of prothrombin time and presence or absence of spiders. As the study was extended for another four years, a total of 1945 visits were generated. All of the twelve covariates except age become time-varying covariates in the follow-up data. At the end of the follow-up study, 169 of the 312 patients were still alive, 140 had died, and three had been lost contact with. The extended study allows the researchers to study the effects of the changes in the prognostic variables, as time-varying covariates (Murtaugh et al., 1989). We therefore fit the proposed forest methods as well as the extended Cox model (Murtaugh et al., 1989) on the dataset with time-varying covariates from the extended study. Note that estimates of this kind for different sets of covariate values over time can be useful in providing guidance from a public policy point of view, as they highlight the different average survival experiences of different subpopulations with different covariate paths. To better illustrate the effects of the time-varying covariates, we also consider the corresponding time-invariant dataset where the covariate values are never updated after the initial observation. For performance comparison, the Brier score-based 10-fold cross-validation error at the t-th recording day is computed for each method on both the time-varying and time-invariant covariate datasets, shown in Figure 6. The corresponding integrated Brier score cross-validation results are given in Table 4.  Figure 6 shows that the cross-validation errors from all methods on the time-varying covariate dataset are lower than the corresponding ones in the time-invariant covariate dataset after t > 300, which suggests that the updated covariate information can significantly improve performance. Given the results on the time-varying covariate data in Figure 6, one can see that the differences in performance are relatively small. The extended Cox model slightly outperforms the others before t = 3000, where around 30% of the subjects are still at risk. CIF-TV and TSF-TV give the best performance between t = 3000 and t = 4200. After t = 4200, since the results are based on very few failures (given only 10% of subjects are still at risk), the Brier scores are highly variable and not as reliable. Note that the Brier score results also suggest that the survival relationship between the hazards and covariates is relatively (log-)linear for t < 3000, but not for t > 3000. Based on these Brier score-based cross-validation results, we would recommend an extended Cox model for 0 < t < 3000, and CIF-TV for t > 3000. Table 4 shows that CIF-TV gives the lowest corresponding integrated Brier score cross-validation errors, suggesting that CIF-TV provides the best estimated accuracy overall. In practice, where the underlying true survival distribution may possess a complex structure with time-varying features, a plot of the cross-validation-based Brier scores as in Figure 6 can potentially help data analysts make better choices for survival estimation at different time points, compared to a universal choice of method for all time points.

Discussion
The estimation of a population-level survival probabilities for time-varying covariate data are useful in many survival analysis settings, causal inference research being an important example. Additional challenges arise when considering time-varying treatment regimens. To draw real-world evidence about the effectiveness of such regimens on patient survival, the key is to account for the time-varying confounding effects and one way to address this issue is by using the inverse probability of treatment weighting (Barber et al., 2004;Robins et al., 2000). For survival data, the estimation of the time-varying weights can potentially be improved by using flexible tree-based methods allowing time-varying covariates (confounders). A recent work (Hu et al., 2022) that uses LTRC forests shows that the use of more flexible models for the estimation of time-varying weights can lead to more accurate treatment effect estimation.
In this paper, we propose extending the relative risk, conditional inference, and transformation forests to provide survival function estimation with time-varying covariates. The extension of random survival forests (Ishwaran et al., 2008) for cumulative hazard function estimation with time-varying covariates based on log-rank-type splitting rules requires reformulation of the log-rank test statistics (LeBlanc and Crowley, 1993;Segal, 1988) and/or the log-rank scores (Hothorn and Lausen, 2003). For a log-rank score rule, one can make use of the log-rank score specified in (4). However, in this case one cannot construct the extension by modifying the custom splitting rule feature and employing the fast C code from randomForestSRC, since its basic structure only allows for right-censored survival. A fast and efficient algorithm to construct random survival forests for left-truncated right-censored survival data is needed, and is a goal of future work.

Conclusion
In this paper, we have proposed two new ensemble algorithms, CIF-TV and RRF-TV, and adapted the transformation algorithm, TSF-TV. These three forest algorithms can handle (left-truncated) right-censored survival data with timevarying covariates and provide dynamic estimation. The tuning parameters in the proposed forest methods for survival data with time-varying covariates affect their overall performance. Guidance on how to choose those parameters is provided to improve on the potentially poor performance of forests with the default parameter settings.
The estimation performance of the proposed forest methods is investigated to understand how the improvement over a Kaplan-Meier fit is related to changes in different factors. Focusing on the more influential factors, the estimation performance comparison against other methods shows that the proposed forest methods outperform others under certain circumstances, while no method can dominate in all cases. We then provide guidance for choosing the modeling method in practice, showing that cross-validation is able to pick the best modeling method most of the time, or at least select a method that performs not much worse than the best method.
The estimation of a population-level survival probabilities for time-varying covariate data are useful in many survival analysis settings, causal inference research being an important example. Additional challenges arise when considering time-varying treatment regimens. To draw real-world evidence about the effectiveness of such regimens on patient survival, the key is to account for the time-varying confounding effects and one way to address this issue is by using the inverse probability of treatment weighting (Barber et al., 2004;Robins et al., 2000). For survival data, the estimation of the time-varying weights can potentially be improved by using flexible tree-based methods allowing time-varying covariates (confounders). A recent work (Hu et al., 2022) that uses the proposed survival forests shows that the use of more flexible models for the estimation of time-varying weights can lead to more accurate treatment effect estimation.
Our developed methodology and algorithms allow for estimation using the proposed forests for (left-truncated) rightcensored data with time-invariant covariates. The same data-driven guidance for tuning the parameters or selecting a modeling method also applies to the time-invariant covariates case (for both left-truncated right-censored survival data and right-censored survival data), which implies its broad effectiveness regardless of additional left-truncation and regardless of the presence of time-varying effects.
6 Code availability R scripts for reproducibility of the simulations and real dataset illustrative example analysis are available from the github repository, https://github.com/WeichiYao/TimeVaryingData_LTRCforests.
An R package, LTRCforests, which implements CIF-TV and RRF-TV for LTRC data with application to time-varying data, is available on CRAN.

Appendices A Dynamically Adjusted Survival Function for the Hypothetical COVID-19 Example
The time-to-event is the time to a positive COVID test. There is a time varying covariate, X(t), which describes the vaccination status of a subject with values in {0, 1, 2, 3} where 0 = unvaccinated, 1 = vaccinated with a single dose and 2 = vaccinated with two doses, and 3 = vaccinated with two doses and a booster. The covariate may change values at discrete time points; in this example at times 1 and 2.
We construct the hypothetical tree for COVID progression, see Figure 7. At time t = 0, there is a sample of n = 100 unvaccinated and COVID-free subjects at the root of the tree, i.e., every subject has X(0) = 0. Thus, the estimated survival function at time t = 0, S(0) = 1. Each of the branches on the tree are labeled b 0 , b 00 , etc. The initial branch is b 0 ; it represents a subject not vaccinated in the time interval [0, 1). Similarly, branch b 00 represents a subject unvaccinated over the time interval [0, 2); branch b 01 , a subject vaccinated with at single dose at time t = 1; b 011 , a subject vaccinated with a single dose at time t = 1 and not receiving second dose at time t = 2, etc. For simplicity, we will assume that there is no censoring.
At time t = 1, there are five subjects infected with COVID and 95 COVID-free. Hence, the estimated survival function at time t = 1 is For the group of three subjects with no second vaccine at time t = 2, the estimated survival probability at time t = 3 is lower: At time t = 2, of the 20 unvaccinated COVID-free subjects, 12 receive their first dose of vaccine and eight do not. Of those 12 subjects, nine are COVID-free at time t = 3. For this group, which gets a single vaccine at time t = 2, the estimated survival probability at t = 3 is S(3 | b 001 ) = S(2 | b 00 )(9/12) = 0.76(0.75) = 0.57, which happens in this example to be the same as the estimated survival probability of the group that got the single vaccine at time t = 1.
The subjects who did or did not receive a booster at the next time period would be handled in a similar way.

B Derivation of the Survival Estimate
Recall that by definition of survival functions, S(t|X * (t)) = P(T > t|X * (t)). At given time t ∈ [t * j , t * j+1 ), note that P(T > t|X * (t)) = P(T > t, T > t * j | X * (t)), we apply the conditional probability and obtain In constructing the survival function estimate, we assume that the hazard at time t is a function only of the current covariate values at time t (but these covariates can include lagged values of some covariates). This allows us to construct the estimate at time t using any subjects in the population with the specified value at that precise time point; that is, we estimate P(T > t | T > t * j , X * (t)) by computing where both the numerator and the denominator are the values of the estimated survival function in the hypothetical case with the covariate x * j at t and t * j , respectively. The risk sets that are used to compute these two quantities consider all subjects with covariate values x * j at t * j , regardless of their covariate paths before t * j . Note that this hypothetical estimated survival function is in fact the output of the algorithm for the input with covariate value x * j .
Therefore, by substituting S A,j (t) P(T > t | x * j ), we can then approximate (20) as which gives us the formula in (9).

Supplemental Material for "Ensemble Methods for
Contents S1 Proposed forests for time-varying covariate survival data 3 S1. Suppose the covariates X are initially observed at t 0 = 0, and change values at t 1 , . . . , t m−1 , before the subject is censored or the event occurs at t m . That means there are in total m intervals, within which the values of all covariates remain the same, The survival times are generated as follows: 1. Compute the values of cumulative hazard function H(t) at t 1 , t 2 , . . . , t m−1 from

Divide time lines into the following regions
3. Randomly generate u ∈ U(0, 1), the survival time is determined by which region − log(u) falls in:

The corresponding survival probability function is
where t falls into mutually exclusive intervals

S1.2 Changing pattern of the time-varying covariates in the simulations
This section provide further details on the changing pattern of the time-varying covariates, specifically, on how to generate the observed values of time-varying covariates X 6 , X 13 , X 16 , X 18 , X 20 at t 0 , . . . , t m−1 : x j,k , for j = 0, . . . , m − 1 and k = 6, 13, 16, 18, 20.

S1.3 Histogram of simulated survival times
Histograms of survival times for typical samples with the number of subjects N = 500 in each scenario are provided in Figure S1.1.

S1.4 Parameters set in the simulation study
In the simulation study, we set the parameters for the basic scenario "2TI + 4TV" as follows.
One can see that for each type of forest, the results when bootstrapping pseudo-subjects are very similar to those when bootstrapping the subjects. That is, the two different bootstrapping mechanisms do not result in fundamentally different levels of performance. S1.6 Regulating the construction of trees in forests Datasets are generated with light right-censoring rate (20%), survival times generated from a Weibull-Increasing distribution. From the top row to the bottom, it gives results for the linear, nonlinear and the interaction survival relationship. From the first column to the last, it gives results for the number of subjects N = 100, 300, 500. In each plot, 1-TSF-TV with mtry = 1; 2-TSF-TV with mtry = 2; 3-TSF-TV with mtry = 3; 5-TSF-TV with mtry = 5; 10-TSF-TV with mtry = 10; 20-TSF-TV with mtry = 20; Opt-TSF-TV with value of mtry that gives the smallest Integrated L 2 difference in each round; Tuned-TSF-TV with the value of mtry tuned by the "out-of-bag" tuning procedure. Table S1.1 gives the performance comparison between each forest method with its default parameter settings and with the proposed parameter settings for the nonlinear survival relationship under the PH setting and the non-PH setting, respectively. Table S1.2 gives the performance comparison between each forest method with its default parameter settings and with the proposed parameter settings for the interaction survival relationship under the PH setting and the non-PH setting, respectively.  Figure S1.9: Boxplots of integrated L 2 difference for performance comparison under the PH setting. Datasets are generated with survival times following a Weibull-Decreasing distribution, light rightcensoring rate 20%. The first row shows results for the number of subjects N = 100, second row for N = 300, bottom row for N = 500; the first column shows results for linear survival relationship, second column for nonlinear, the third column for interaction. The horizontal dashed line shows the median integrated L 2 difference of a Kaplan-Meier fit on the datasets. In each of the plots, the set of boxplots lightly shaded shows the performance of different methods on datasets with history of changes in covariates' values known; the set heavily shaded shows the performance on datasets with part of history of changes in covariates' values unknown.

S2 Proposed forests for time-invariant covariate survival data
The section is organized as follows. Section S2.1 describes the model setup, Section S2.2 provides the performance evaluation of the "out-of-bag" tuning procedure for mtry and gives comparative performance between the proposed parameter settings and the default settings. The overall performance comparison among the Cox model, the three forest methods with the proposed parameter settings, as well as the best method and the method selected by the IBSbased 10-fold CV rule are given in Section S2.3. All results with integrated L 2 difference are computed with τ i = T i . The comparative performance of the different methods was broadly similar when using τ i = τ = max j ‹ T j .
Here we use the proposed forest methods for left-truncated right-censored survival data with time-invariant covariates, the proposed forest methods and the transformation forest method are referred to as LTRC CIF, LTRC RRF and LTRC TSF, respectively.

S2.1 Model setup
We generate left-truncated right-censored survival time data with time-invariant covariates as follows. The left-truncation time T 0 is generated as a U(0, L) random variable with some constant L > 0. The event time T is randomly generated with a Weibull distribution. If the generated T < T 0 , i.e. the event time is less than the left-truncation time, then this observation is discarded. Otherwise, the observation is retained, with censoring time C = T 0 + D, where D ∼ weibull(shape = 2, scale = λ D ) has an weibull distribution. The parameter λ D is selected to ensure 20% censoring rate. If C < T , then this observation is censored (∆ = 0), otherwise the survival time T is observed (∆ = 1). Note that D and T 0 are both independently generated from T and from each other. The observed response for each observation is a triplet (T 0 , T , ∆), where Y = min(T, C).

S2.2 Regulating the construction of trees in forests
In the simulations the value of mtry is tuned based on the "out-of-bag observations," and the values of minprob and minbucket are set to be the maximum of the default value and the square root of the number of pseudo-subject observations n. In each plot, 1-LTRC CIF with mtry = 1; 2-LTRC CIF with mtry = 2; 3-LTRC CIF with mtry = 3; 5-LTRC CIF with mtry = 5; 10-LTRC CIF with mtry = 10; 20-LTRC CIF with mtry = 20; Opt-LTRC CIF with value of mtry that gives the smallest Integrated L 2 difference in each round; Tuned-LTRC CIF with the value of mtry tuned by the "out-of-bag" tuning procedure. The default value in LTRC CIF is mtry = 5. From the first column to the last, are given results for the number of subjects N = 100, 300, 500, 1000.
In each plot, 1-LTRC CIF with mtry = 1; 2-LTRC CIF with mtry = 2; 3-LTRC CIF with mtry = 3; 5-LTRC CIF with mtry = 5; 10-LTRC CIF with mtry = 10; 20-LTRC CIF with mtry = 20; Opt-LTRC CIF with value of mtry that gives the smallest Integrated L 2 difference in each round; Tuned-LTRC CIF with the value of mtry tuned by the "out-of-bag" tuning procedure. The default value in LTRC CIF is mtry = 5. Opt-LTRC TSF with value of mtry that gives the smallest Integrated L 2 difference in each round; Tuned-LTRC TSF with the value of mtry tuned by the "out-of-bag" tuning procedure. The default value in LTRC TSF is mtry = 5. Table S2.1 gives performance comparison between each forest method with its default parameter settings and with the proposed parameter settings, using datasets generated with survival times following a Weibull-Increasing distribution, nonlinear relationship, with 20% right-censoring rate, as an example.
In Table S2.1, positive numbers indicate a decrease in integrated L 2 difference compared to a Cox model on the dataset, while negative numbers indicate an increase. The absolute value of the numbers represents the size of the difference between the integrated L 2 difference of the candidate and that of a Cox model. The two tables show that forests with the proposed parameter settings can provide improved performance over those with default parameter settings across almost all different numbers of subjects N by a substantial amount, under both PH and non-PH settings.

S2.3 Performance comparison and IBS-based CV model selection
Figures S2.7 and S2.8 give side-by-side integrated L 2 difference boxplots on datasets with light (right-)censoring rate, survival times generated following a Weibull-I distribution under the PH setting and the non-PH setting, respectively. These figures compare the performance of the four methods, the best method and the method chosen by IBS-based 10-fold CV rule with τ i = 1.5 ‹ T i .