{'fraction_non_alphanumeric': 0.0305811468805829, 'fraction_numerical': 0.011559883700563982, 'mean_word_length': 4.886185567010309, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Landmark Proportional Subdistribution Hazards Models for Dynamic Prediction of Cumulative Incidence Functions

An individualized risk prediction model that dynamically updates the probability of a clinical event from a specific cause is valuable for physicians to be able to optimize personalized treatment strategies in real-time by incorporating all available information collected over the follow-up. However, this is more complex and challenging when competing risks are present, because it requires simultaneously updating the overall survival and the cumulative incidence functions (CIFs) while adjusting for the time-dependent covariates and time-varying covariate effects. In this study, we developed a landmark proportional subdistribution hazards (PSH) model and a more comprehensive supermodel by extending the landmark method to the Fine-Gray model. The performance of our models was assessed via simulations and through analysis of data from a multicenter clinical trial for breast cancer patients. Our proposed models have appealing advantages over other dynamic prediction models for data with competing risks. First, our models are robust against violations of the PSH assumption and can directly predict the conditional CIFs bypassing the estimation of overall survival and greatly simplify the prediction procedure. Second, our landmark PSH supermodel enables researchers to make predictions at a set of landmark points in one step. Third, the proposed models can easily incorporate various types of time-dependent information using existing standard software without computational burden.


Introduction
An individualized risk prediction model that estimates the probability of a clinical event dynamically using time-dependent information can enable physicians to optimize personalized treatment strategies in real-time. For example, for an early-stage breast cancer patient who received surgery, the physicians can adjust the approach to the subsequent chemotherapies if knowing the risk of a distant metastasis within the next few years given the timing of the patient's occurrence of locoregional recurrence. Note that the occurrence of locoregional recurrence is a time-dependent intermediate event. Other types of time-dependent information include longitudinal measurements of biomarkers and time-varying covariate effects, which are also crucial to the success of a risk prediction model. Two popular approaches in dynamic prediction that received significant attention are joint modeling (Proust-Lima and Taylor 2009, Rizopoulos 2011, Mauguen et al. 2013, Rizopoulos et al. 2014) and landmarking (van Houwelingen 2007, van Houwelingen and Putter 2008, Parast et al. 2012. To incorporate the time-dependent information, joint modeling captures the processes of longitudinal and time-to-event data simultaneously whereas the landmark method applies time-to-event modeling on patients who are still at risk at the updated prediction baseline, i.e., the landmark time, using the information collected up to the landmark time. Although the joint models provide additional information from the joint distributions, the implementation is complex and the correlation structure is difficult to specify correctly thus often results in inaccurate prediction. By contrast, the landmark Cox model and supermodels proposed by van Houwelingen (2007) showed significant advantages in simplicity, flexibility, and accuracy in estimating the risk. Moreover, they become robust against model misspecifications by avoiding the full specification of the joint distribution. These models can easily incorporate various types of time-dependent information without using computationally intensive estimation procedures, and the estimation can be implemented in everyday clinical practice by using standard statistical software. Recently, landmark modeling was used in many clinical applications (Zamboni et al. 2010, Fontein et al. 2015. To account for data with competing events, Cortese and Andersen (2009) and Nicolaie et al. (2013a) applied the landmark method onto the cause-specific hazards model. The disadvantage of this approach is that to estimate the risk from one cause we need to construct the landmark models or supermodels based on the cause-specific hazards for all causes. For direct modeling cause-specific risk, the Fine-Gray model (Fine and Gray 1999), pseudo-observation approach (Klein and Andersen 2005), and direct binomial regression model (Scheike et al. 2008) are the commonly used methods, and the Fine-Gray model is the most widely used one because the estimated risk can be easily calculated. Cortese et al. (2013) applied the landmark method to the Fine-Gray model at a small set of fixed landmark time points, however, their method did not assess robustness of their proposed model against violations of the proportional subdistribution hazards (PSH) assumption, nor did they construct the comprehensive landmark supermodel to estimate the risks of a sequence of prediction time points simultaneously. The major challenge encountered in developing a landmark model/supermodel based on the Fine-Gray model was in constructing the landmark subset at each landmark time to properly account for competing events in the setting of subdistribution hazard (Nicolaie et al. 2013a). Nicolaie et al. (2013b) developed a landmark supermodel based on the pseudoobservations created at each landmark time point, however, additional steps are required to calculate the risks.
In this study, we overcame the aforementioned challenges, extended the landmark method to the Fine-Gray model, and proposed landmark PSH model and landmark PSH supermodel for dynamic prediction in the presence of competing risks. Our proposed models can directly predict the risk for a specific cause and provide accurate prediction even when the PSH assumption fails to hold. The landmark PSH supermodel uses smoothing techniques to aggregate a set of simple landmark PSH models, allowing users to make dynamic predictions at arbitrary landmark points using only one model, also provides a simple and explicit estimation form to incorporate time-dependent covariates and/or time-varying covariate effects. Compared to the currently available landmark dynamic prediction models for competing risks, our proposed models are more straightforward and can be effortlessly implemented with standard statistical software. This paper is organized as follows. In Section 2, we first introduce notations and dynamic predictive quantities for data with competing risks, and then review the Fine-Gray PSH model and present our proposed landmark PSH model and supermodel. In Section 3, we assess the performance in prediction for our proposed models and compare with the existing methods using simulations. In Section 4, we apply the proposed models to predict the dynamic cumulative incidences of distant metastasis based on the locoregional recurrence status and a set of prognostic factors of breast cancer patients in a multicenter clinical trial. The discussions are provided in Section 5.

Notations and Conditional CIF
Let T and C be the failure and censoring time, respectively; ε ∈ {1, . . . , k} be the cause of failure; and Z(·) be a p-dimensional vector of covariates, which could be time-fixed covariates measured at the baseline or be time-dependent covariates. Here, for each subject i (i = 1, . . . , n), we assume C i is independent of T i and Z i (·), and refer to it as random censoring.
For right censored data, we observe an independently and identically distributed quadruplet Suppose cause 1 is the primary event of interest, the subdistribution function or CIF of cause 1 is defined as F 1 (t) = Pr(T ≤ t, ε = 1).
In dynamic prediction, the aim is to predict the conditional CIF that is the cumulative incidence of occurring the primary event within the next w time units (also named as a prediction window with a fixed width of w) for a subject who has not failed from any cause at a landmark time (denoted by s) by incorporating the time-dependent information available up to the landmark time s. The width of prediction window w depends on the clinical relevance that defines the time interval within which the risk of having the primary event is of interest.
Then, the conditional CIF is defined as where Z(s) is the covariates Z(·) whose values are measured or available at s. Using the definition of conditional probabilities, we can rewrite the conditional CIF as where F j {s|Z(s)} is the CIF for cause j and S{s|Z(s)} is the overall survival at time s.  (2007) to the competing risks setting and propose a landmark PSH model which is also robust against the misspecification of the proportional subdistribution hazards.

Fine-Gray Model
Fine and Gray (1999) proposed a regression model with respect to the subdistribution hazard With time-dependent covariates Z(t), the Fine-Gray model takes the form λ 1 {t|Z(t)} = λ 10 (t) exp{Z(t) T β}; and The regression coefficients β are estimated through a partial likelihood approach with modified risk sets which are defined as who have not failed from the event of interest by time T i .

Landmark PSH Model
To directly estimate the conditional CIF F 1,LM {s + w|Z(s), s} at a given landmark time s, we define a conditional subdistribution hazard for all t ≥ s, which is the subdistribution hazard λ 1 {t|Z(s)} conditional on assuming no event occurred from any cause prior to s. It can be shown that (1) is the hazard function of the conditional CIF: where Λ 1 {·|Z(s), s} is a cumulative conditional subdistribution hazard at landmark time s.
Therefore, we propose a landmark PSH model with the form for s ≤ t ≤ s + w.
To estimate the regression parameters β LM and the baseline conditional subdistribution hazard λ 10 (t|s), we can simply apply the Fine-Gray model to subjects who have not yet failed from any cause at s and ignoring the events occurring after s + w by adding an administrative censoring at s+w. As discussed in Liu et al. (2016), under non-PSH, the Fine-Gray model with an administrative censoring at the horizon time provides an accurate prediction of CIF at this horizon time (but not other time points) without constructing complex procedures to estimate time-varying effects; and the partial likelihood estimatorβ LM of β LM in model (2) (2012), this provides a dynamic prediction framework with a sliding prediction window over the specified interval. Theoretically we can fit L + 1 landmark PSH models (2) to obtain L + 1 conditional CIFs. When the number of landmark points increased dramatically, it is less practical and computationally inefficient to fit many models separately.
Instead, we adopt the smoothing strategy of van Houwelingen (2007) and propose a landmark PSH supermodel with the form The parameters (θ, η) can be consistently estimated by maximizing a Breslow pseudo partial log-likelihood for tied events, because that in the stacked dataset, one subject with event time T i has n i = #{s : s ≤ T i ≤ s + w, s ∈ [s 0 , s L ]} repeated observations; andΛ 10 (t) can also be estimated by a Breslow estimator for the cumulative subdistribution hazard (see Web Appendix A for the details). To obtain the standard errors for the estimated parameters, a robust sandwich estimator is required to adjust for the correlation between the risk sets which exists because the same subject is repeatedly used when we estimate the parameters based on the stacked dataset. Thus, the target dynamic prediction probabilities F 1,LM {s + w|Z(s), s} have a simple and explicit estimation form, which is given bŷ for all s ∈ [s 0 , s L ]. Therefore, the landmark PSH supermodel can provide the prediction of For implementation, before stacking all landmarking subsets into a super dataset, we need to transform each subset into the counting process style and include time-varying IPCW weights for the subjects experienced competing events to adjust for the random right censoring (Geskus 2011). In practice, fitting the landmark PSH supermodel in the stacked dataset requires a software that allows for delayed entry or left truncation at s. The R function coxph() can be used to fit model (3) and it can also provide robust sandwich covariance matrix for (θ, η) which can be used in the significance test of the estimated regression coefficients. As discussed, the landmark effect on λ 10 (t|s) is through β LM (s), so that there is a correlation between θ and η. It is suggested to center the covariates before fitting the model (van Houwelingen and Putter 2012).
For simplicity, we assumed that C is independent of T and Z, but, it can be generalized to a conditional independence between C and T given Z, as discussed in the PSH model (Fine and Gray, 1999). The dependence between C and Z can be handled by modeling C as a function of Z in IPCW. Following the same strategy, the properties of landmark PSH model and landmark PSH supermodel will be retained under conditional independent censoring.

Measure of Predictive Performance
To evaluate dynamic predictive performance of the proposed procedures, we adapted the timedependent Brier score, O/E ratio, and the area under the receiver operating characteristics curve (AUC) in terms of predictive accuracy, calibration, and discrimination, respectively.
The time-dependent Brier score is an estimate of the mean-squared prediction errors of the predicted event probabilities at s+w over the observed event status for subjects who are event free at landmark time s (Schoop et al. 2011, Cortese et al. 2013. For data with competing risks, the expected time-dependent Brier score at landmark time s for the prediction at s + w is defined as When random right censoring exists, we utilized a pseudovalue-based consistent estimator 1,LM (s + w|s) is a jackknife pseudovalue for the ith subject who is still at risk at time s, whereF 1,LM (s + w|s) is the nonparametric estimate of the marginal cumulative incidence Pr(T ≤ s + w|T > s), and F (i) 1,LM (s + w|s) is the same estimate but is based on the data where the ith subject has been removed.
To assess calibration, we adapted the time-dependent O/E ratio which is defined as the observed number of main events from the landmark time point divided by the prediction horizon time to the expected number of events estimated from the predictive model, As discussed in Pfeiffer and Gail (2017), to address the right censoring, the jackknife pseu- To measure discrimination capability of the landmark PSH supermodel, we followed the work of Blanche et al. (2015) and Huang et al. (2016) to calculate the time-dependent AUC which is given by

Simulation Studies
We evaluated the performance of the proposed dynamic prediction models using simulated data under two different settings of the nonproportional subdistribution hazards (non-PSH).
Because that nonparametric method is model-free, we used it as the reference to compare with other model-based methods in scenarios where only categorical variables are included.
For simplicity, only two failure types were considered: type 1 failure is the primary event of interest; type 2 failure indicates competing events. We also simulated a competing risks data under non-PSH with a time-dependent covariate.
In the first non-PSH setting, we generated the type 1 failure times from a two-parameter Weibull mixture distribution with the subdistribution F 1,1s (t; where the rate parameter depends on the covariate Z i . In the second non-PSH setting, we let the coefficient of Z i be a function of time; and the subdistribution of the primary event became We generated the type 2 failure times from an exponential distribution Pr(T i ≤ t|ε i = 2, Z i ) = 1 − exp{− exp(Z i β c t)} by taking Pr(ε i = 2|Z i ) = 1 − Pr(ε i = 1|Z i ), where β c = 0.5. Sample size of n = 1, 000 was chosen and the data were simulated repeatedly for N = 1, 000 times. The censoring times were generated independently from a uniform distribution.  We also provided the estimates of the conditional failure probability from the landmark Cox supermodel (van Houwelingen 2007) by treating the competing risks as random censoring. We found that in two different non-PSH scenarios, the performance of the landmark PSH model and the landmark PSH supermodel are as good as that of the nonparametric methods. But the landmark Cox supermodel overestimated the risk of failure for main event in the presence of competing risks.
For each approach we evaluated the prediction errors in the dynamic conditional CIFs by estimating time-dependent Brier scores, and we used a 3-fold cross-validation to correct for possible overfitting. Web Table 1 presents averaged estimates of the cross-validated timedependent Brier score and its empirical standard deviation. To quantify the improvement of predictive accuracy for the proposed landmark PSH model and supermodel to the standard PSH model under nonproportional hazards, we utilized a relative increment (or reduction) of prediction errors by treating the nonparametric estimates as the reference. The relative increment of prediction errors are presented in Figure 2. As expected, in both non-PSH settings the predictive accuracy of the landmark PSH model was almost the same as that obtained from the nonparametric method. As compared with the landmark PSH model or the nonparametric method, the landmark PSH supermodel has slightly lower accuracy, yet the differences in prediction errors are negligible.
To further evaluate predictive performance of the proposed landmark PSH supermodel for data with time-dependent covariates, we added Z i (t) into the second non-PSH setting using the simulation strategy introduced in Huang et al. (2016) (see Web Appendix B for details).
We chose a prediction window of w = 0.4 and sample size of n = 5, 000 to calculate the 3-fold cross-validated time-dependent O/E ratio, Brier score, and AUC for the landmark PSH supermodel and compared it with the landmark Cox supermodel. To fit the supermodels, we selected landmark time points from 0 to 4 with a step of 0.1, and used the same quadratic basis function f (s) and g(s) for Z i (t) because that only the effect of Z i (t) depends on landmark s in the Wald test. Table 1 provides averaged estimates of the O/E ratio, Brier score, and AUC with the corresponding empirical standard errors from 1, 000 simulated data sets. The landmark PSH supermodel showed a good predictive performance in terms of calibration, predictive accuracy, and discrimination.

Application
To illustrate the application of our proposed landmark PSH supermodel in prediction of the conditional cumulative incidence probabilities for a moving (or dynamic) time interval with a fixed width, we used the data from a multicenter phase III clinical trial for breast cancer patients with estrogen receptor positive and historically nodes negative (Fisher et al. 1997). In this trial, 2,363 patients were randomly assigned to receive one of the following three regimens: tamoxifen 10mg daily for 5 years, tamoxifen 10mg daily for 5 years plus metrotrexate (M) and Fluorouracil (F), and tamoxifen 10mg daily for 5 years plus M, F, and cyclophosphamide (C); denoted TAM, TAM+MF, and TAM+CMF, respectively for simplicity. The median follow-up time was 11.2 years.
Among the 2,272 clinically eligible patients who were followed, 241 developed distant metastasis, 127 died due to other causes before distant metastasis could occur, and the remaining 1,904 were censored due to withdrawal of the study, lost to follow up, or event-free on or before the analysis cutoff date. In analyses, the censoring cases were treated as random noninformative censoring where the justification is described in Web Appendix C. For early stage breast cancer patients after surgery, development of locoregional recurrence (LRR) is an important prognostic clinical event affecting the risk of distant metastasis. In this data, 15.1% of the patients developed LRR before progressing to distant metastasis; but only 6.6% patients experienced LRR before death. Our main interest in this application is to dynamically predict the risk of distant metastasis within the subsequent 3-years for a breast cancer patient, based on her LRR status measured during follow-up and other prognostic covariates measured at baseline, including the treatment (TAM, TAM+MF, or TAM+CMF), surgery type (lumpectomy plus radiation therapy [L+XRT] vs. mastectomy), age at the study entry (< 50 vs. ≥ 50 years old), clinical tumor size (≤ 2 vs. > 2 cm), and tumor grade (well, moderate, and poor). We also compared the dynamic 3-year fixed width probabilities of distant metastasis and death based on a patient's LRR history.
Web Figure 2(a) shows that the estimated CIFs for both distant metastasis and death in which the maximum event time around 13 years. Web Figure 2(b) depicts the estimated CIF of LRR. There were only a few random censoring events occurred during the first 10 years of follow-up (figure not shown). Thus, we chose the range of landmark points from 0 to 10 years and prediction window with a fixed width of 3-years. To fit a landmark PSH supermodel to this dataset, we took 51 equally spaced landmark points s (0 ≤ s ≤ 10) and set the basis functions for β LM (s) and γ(s) as β LM (s) = θ 1 + θ 2 s + θ 3 s 2 and γ(s) = η 1 s + η 2 s 2 , respectively.
The frequencies of distant metastasis and death in each of the landmark sub-dataset for s = 0, 1, . . . , 10 years are shown in Web Figure 2(c).
We began our analysis using the backward selection procedure to select those covariates of which the effects were dependent on the landmark points. We tested the landmark-covariate interactions for each covariate via the Wald test based on the robust covariance matrix of the estimated coefficients from the landmark PSH supermodel. We found that for distant metastasis, the effects of TAM+MR and age were significantly dependent on the landmark points, whereas for death, only the effect of LRR status was dependent on the landmark points. The estimated coefficient and the corresponding robust standard error for a given prognostic factor are summarized in Table 2. Multivariate Wald tests for the baseline parameters (η 1 , η 2 ) were significant for both distant metastasis and death, indicating that the baseline subdistribution hazard also depended on the choice of landmark points.
Given different LRR status ( including no LRR developed over the course of study, with the first LRR occurred at 3, 5, and 7 years) and different treatment regimens (TAM, TAM+MF, or TAM+CMF), Figure 3 depicts the predicted dynamic 3-year fixed width cumulative incidences of distant metastasis and death with the associated bootstrap 95% confidence intervals for patients younger than 50 years old with poorly differentiated tumors, tumor size larger than 2cm and treated by L+XRT. If a patient did not develop LRR, the risk of having distant metastasis within the subsequent 3 years is very close to the risk of death in any treatment group. However, if the patient occurred LRR, she would have a much higher risk of developing distant metastasis compared to the risk of death, especially for TAM+MF treatment group. Similar trends were obtained for patients if they experience the first LRR at 3, 5, and 7 years, respectively.
We compared the risk prediction performance of our proposed landmark PSH supermodel against the landmark Cox supermodel. The fitted landmark Cox supermodel is summarized in Web Table 2. In addition to measure the time-dependent O/E ratio, Brier score, and AUC at a selected set of landmark time points with bootstrap standard errors with B = 500, we also calculated the squared bias (Bias2) for both models. As results shown in Table 3, the landmark PSH supermodel is well-calibrated compared to the landmark Cox supermodel which tends to underestimate the risk of event with larger O/E ratios. Although two models provide similar Brier score results, the Cox supermodel has notably larger squared bias in predictive accuracy, which indicates that the Brier score is dominated by the Bernoulli variation of the outcome not by the squared bias in this application data, so leads to a limited variation in Brier score between two models. Because of the two models have the similar model form (Table 2 and   Web Table 2), the difference of AUCs are limited in this application example.
To demonstrate an additional value that dynamic prediction provides, we compared the performance of the landmark PSH supermodel that dynamically incorporated LRR status at each landmark time to the PSH model that simply used LRR status as a predictor to predict the CIFs at the horizon time. Results are summarized in Web Table 3. Apparently, the landmark PSH supermodel outperformed the PSH model in terms of both calibration and predictive accuracy.

Discussion
In this study, we developed dynamic prediction models for data containing competing risks by extending the landmark approach to the Fine-Gray model. The resulting landmark PSH model and the supermodel can be used to directly predict the dynamic cumulative incidences for the occurrence of a specific event within a given prediction window of a fixed width through incorporating all available information up to the landmark time under the condition that the patient has not failed at the landmark time.
Our proposed models have several advantageous features over the currently available methods in predicting the dynamic cumulative inferences. First, our landmark PSH models can provide accurate estimates even when covariate effects are potentially time-varying and violate the PSH assumption whereas the model developed by Cortese et al. (2013)  Furthermore, the proposed landmark PSH supermodel is more straightforward and simpler to implement using the existing standard software packages. In contrast with the landmark supermodel based on the cause-specific hazards (Nicolaie et al. 2013a), our landmark PSH supermodel predicts the conditional CIFs only in one step and provides a direct interpretation of the landmark-specific effects on the predictive probabilities. Compared with the landmark supermodel based on the pseudo-observations (Nicolaie et al. 2013b), our supermodel also shows simplicity in computation whereas the GEE-based method used in pseudo-observations would have convergence issues for large sample size especially when dealing with many landmark points of interest.
We evaluated the prediction performance of our proposed models through simulations and compared to other existing methods. We determined how closely the estimated conditional CIFs would approximate the true probabilities for our models, which is the first of its kind done in this area of study. We also evaluated the landmark supermodel in competing risks data from the perspectives of calibration, predictive accuracy, and discrimination capabilities by utilizing the time-dependent O/E ratios, Brier scores and AUCs. The simulation results showed that our models performed well in prediction after incorporating time-fixed and timedependent covariates even when the PSH assumption was violated.
In the landmark PSH supermodel, we assume that the effect of s on baseline hazards λ 10 (t|s) was an additive effect. As discussed in van Houwelingen (2007), this assumption will hold if the follow-up is not too long or the effect of covariates is not too large. If we choose an optimal width w for the prediction window and a suitable range [s 0 , s L ] for the landmark time points, the landmark PSH supermodel can provide a correct approximation for the conditional CIF at time s + w for any s ∈ [s 0 , s L ]. An alternative method is to fit a stratified landmark PSH supermodel by applying a Fine-Gray model with the landmark-covariate interactions Z(s) * f j (s) to the stacked dataset with stratification on s, which estimates the baseline conditional subdistribution hazard for each s.
Note that in the presence of time-dependent covariates, if the covariate values vary too often and too much over time, the estimated covariate effects could be attenuated in the proposed models and subsequently lead to biased estimation of the conditional CIF. This is also shown in our simulation result for data with a time-dependent covariate. To adjust for this issue, which was also discussed by van Houwelingen and Putter (2012), we will explore a set of suitable landmark time points and use an additional monotonically decreasing function to model the attenuation process of the covariate effects as future work.