FastCPH: Efficient Survival Analysis for Neural Networks

With the rise of deep learning applications in many scientific fields, some studies have tried to combine CoxPH functions with deep neural networks for better time-to-event predictions for larger datasets (Katzman et al. 2018; Spooner et al. 2020; Ching, Zhu, and Garmire 2018). The inclusion of a neural network can simplify the a priori covariates selection and make the model adaptively learn them while preserving the effectiveness of CoxPH functions in survival analysis. Our work focuses on the deep learning method to model the survival hazards using CoxPH.

CoxPH is developed under the assumption of continuous data, but there are often tied event times in real datasets. Several ways have been offered to deal with instances where there are ties, including the exact partial likelihood, the Breslow approximation, and the Efron approximation, the latter two being well-tested in theory and widely used practice (Breslow 1975; Efron 1977; Therneau and Grambsch 2000). The most commonly used approach is the Breslow method which simply uses the number of tied events as the exponent in the denominator of the relative risk. The Efron method is thought to produce superior outcomes, yet the formulation is more complicated to implement efficiently. The difference with the Breslow approach is minor when the number of ties is not large (Therneau 1997). Our work use with both Breslow and Efron methods to break ties.

When most studies focus on the application of linear and nonlinear CoxPH to survival datasets, few of them emphasize the implementation of the Cox partial likelihood function itself, which is the foundation of CoxPH-based methods. We review the four most popular implementations of the CoxPH model in neural networks as in Table 1:

• •

  DeepSurv (Katzman et al. 2018) implements a deep learning generalization of the Cox proportional hazards model using Theano and Lasagne. It supports tensor operation, runs in $\mathcal{O}(n)$ by vectorized cumulative sums over the entire input, but assumes there are no tied events.

• •

  Pycox (Kvamme, Borgan, and Scheel 2019) provides a Python package for survival analysis and time-to-event prediction with PyTorch. It computes in $\mathcal{O}(n)$ a cumulative sum of all input hazards but not the true risk sets of the CoxPH function.

• •

  PySurvival (Fotso et al. 2019–) is an open-source Python package for Survival Analysis modeling published in 2019 and compatible with PyTorch. The implementation of the nonlinear CoxPH model is based on Deepsurv (Katzman et al. 2018), while adding an $\mathcal{O}(n^{2})$ index matrix for Efron’s method to handle ties.

• •

  Scikit-survival (Pölsterl 2020a) (sksurv) is a Python library for survival analysis compatible with scikit-learn (Pedregosa et al. 2011), published in 2020. It implements Brewslow and Efron approximation of the CoxPH model in $\mathcal{O}(n)$ using an inner for-loop at each distinct event time when going through all events. However, it does not support deep learning.

It is important for CoxPH-based methods to properly handle ties, support efficient computation, and use the correct approximation at the same time. Simply assuming the absence of ties or ignoring all tied cases is statistically inappropriate. Failure to handle ties and oversimplifications may cause unexpected consequences, especially when the behaviors of neural networks possess randomness and the results can be hard to interpret (Therneau and Grambsch 2000; Borucka 2014). However, none of the existing popular survival analysis packages have CoxPH implemented in both an efficient and correct way for neural networks.
In recent years, people use the CoxPH model on survival datasets with larger scales and data complexity (Spooner et al. 2020). Existing methods sacrifice the correct formula for computational simplicity or the other way around, but why not use a method that strictly follows the well-tested approximations and possesses efficiency at the same time? Here we present the Fast Cox Proportional Hazard model (FastCPH), a computationally efficient and statistically correct method for survival analysis using neural networks.

To analyze the computational efficiency of FastCPH, we compare its runtime with 4 other popular vectorized implementations including PyCox, DeepSurv, PySurvival. PySurvival use risk / fail matrices to compute Efron method. PyCox and DeepSurv are for deep learning purposes and do not support tie approximations. Detailed descriptions of these methods are in Section 1.1.

We assume all events are uncensored and the data is sorted by duration. We first randomly generate data from size 1 to $10^{3}$. We note down the runtime as calculating the negative log likelihood once and report the mean of 5 random trails.

Among the baselines, FastCPH is the most computationally efficient CoxPH implementation that supports Breslow and Efron approximations. As shown in Fig 1, the curve of FastCPH with Breslow method is very close to the ones of PyCox and DeepSurv, which don’t have tie awareness. As the size of data increases, the difference among FastCPH with Breslow, PyCox, and DeepSurv are getting smaller. For the two baselines both using Efron approximation, FastCPH with Efron method is significantly faster than PySurvival. The gap between FastCPH with Efron and Breslow can be caused by the computation of the weighted terms in denominator. The experimental results align with our claims on the linearity of FastCPH’s runtime.

For the breast cancer dataset, WHAS500, and veteran’s lung cancer dataset, we retrieve the data from Scikit-survival package (Pölsterl 2020b) and obtain one-hot encodings to quantify entries such as treatment received, cell types, prior therapy, etc. The number of final entries is shown in the last row of Table 3. The HNSCC dataset is publicly available via the Cancer Imaging Archive with TCIA Limited Access License (Grossberg et al. 2020). The DICOM imaging data is processed using the Med-ImageTools pipeline (Sejin 2022) to extract the computed tomography (CT) images and gross tumor volume (GTV) segmentation masks with uniform voxel spacing for consistent feature extraction. These images and masks are processed into the NIfTI file format, which is a common standard for 3D medical images, and compatible with PyRadiomics. The processed image and GTV masks into PyRadiomics to extract shape, texture, and statistics features (Van Griethuysen et al. 2017).

We use Harrell’s concordance index (C-index) as a metric to evaluate the performance of our model and other baselines. It is a generalization of Area under ROC curve (AUC) regarding censored data, reflecting the accuracy of pairwise orders of the risk function as the output of the model (Harrell et al. 1982; Uno et al. 2011; James et al. 2013). For input $\textbf{x}_{i}$, duration $t_{i}$ and events $\delta_{i}$ (1 if uncensored, 0 otherwise), the C-index is computed as:

We also provide a $\mathcal{O}(n\log n)$ implementation of the C-index using ordered data structures.

We first compare FastCPH-LassoNet with other CoxPH-based methods using linear structures. We use CoxPH linear model (Cox 1972), CoxNet (Simon et al. 2011), GlmNet (Friedman, Hastie, and Tibshirani 2010) and DeepSurv (Katzman et al. 2018) as baselines. The first three are classical CoxPH-based models with different regularizations. DeepSurv is considered the most advanced deep learning CoxPH-based method (Katzman et al. 2018; Spooner et al. 2020).

The implementations of the Cox linear model and CoxNet are from Scikit-survival (Pölsterl 2020a). For the CoxPH linear model, we set $\alpha=10^{-6}$ as the regularization parameter in the ridge regression penalty. CoxNet is the CoxPH model with an elastic net penalty. We use cross-validation for choosing the best $\alpha$ of the regularization path from $10^{-1}$ to $10^{-5}$. For GlmNet, we use the R built-in cross-validation cv.glmnet with the Breslow method to select the model for testing. The number of folds in cross-validation (if used) is 5 for breast cancer and veteran’s lung cancer dataset and 10 for WHAS500 and HNSCC dataset. For FastCPH-LassoNet and DeepSurv, we fix the number of hidden dimensions to 1 and the number of hidden layers to 1. They share the same architecture and setting (ReLU, Adam, lr=$10^{-3}$). For FastCPH-LassoNet, we set $M=10$ and start at $\lambda=10^{-6}$. The prox method of LassoNet is called on the dense model following on a geometric path until the model becomes sparse. That value is divided by 10 to give lambda_start. 5-fold cross validation is used to select the best $\lambda$ value from multiple runs. We use Efron’s method to break ties.

We use stratified sampling w.r.t uncensored/censored events to split the training set (80%) and test set (20%) for each of the datasets. For models with randomness in training, we run 5 trails for each set of hyperparameters and obtain the average performance.

Firstly, FastCPH-LassoNet obtains promising results in different datasets and generally outperform existing CoxPH-based methods as in Table 2. Its C-index is ranked 1, 1, 2, 4 for the breast cancer, WHAS500, Veteran’s lung cancer, and HNSCC datasets, rsp, showing its discrimination ability to provide reliable ranking of survival times based on risk scores. Comparing FastCPH-LassoNet to DeepSurv, FastCPH-LassoNet performs better in 4 datasets using the same achitecture.
Moreover, FastCPH-LassoNet is able to attain an effective recovery of signals given a subset of features. As shown in Table 3 and 2, the model gives excellent performance with sparsity in covariates (C-index=0.70 for 10% feature selected and C-index=0.67 for 18% feature in the HNSCC and breast cancer dataset, rsp). It achieves highest validation accuracy with a low ratio of the number of total feature over the number of total feature. However, when the number of total features is small, the model may not be able to obtain sparsity over covariates, as shown by its result of WHAS500. Fig 2 gives an example of the training curve of FastCPH-LassoNet on the breast cancer dataset. We can see LassoNet optimized properly with FastCPH as the loss function. The sparsity demonstrated in the experiments implies its potential on large-scale, more complicated real world datasets.

In the context of NN methods, we further analyze the performance of FastCPH-LassoNet using more complex architectures. We use DeepSurv as the baseline because it is commonly recognized as the most advanced CoxPH-based deep learning method.

We let FastCPH-LassoNet and DeepSurv share the same architecture (as indicated in parentheses in Table 4) and setting (ReLU, Adam, lr=$10^{-3}$). For both methods, we run 15 trails to give 95% CI. The implementation of FastCPH-LassoNet is the same as in Experiment 2.

FastCPH-LassoNet outperforms Deepsurv with the same architecture on all datasets in our experiments. It is a more robust deep learning method with promising results in survival analysis. Looking at the results in Table 4 and 2 together, FastCPH-LassoNet is the best CoxPH-based method on 3 out of 4 survival datasets we use. Notice that for the breast cancer and WHAS500 datasets, FastCPH-LassoNet yields a significantly better C-index than DeepSurv and other existing methods in Table 2. These results reflect that FastCPH-LassoNet has best overall performance than other CoxPH plus regularization models. To summarize, FastCPH-LassoNet is an efficient and robust way to conduct survival analysis based on the FastCPH and $\textit{L}^{1}$ penalty (Tibshirani 1996).

are common in the datasets we use, listed as in Table 6. Breslow and Efron methods give the same log likelihood when no ties are present in the dataset. Our method in Table 2 using Breslow approximation is capable of handling datasets with and without ties.

is a challenging prediction dataset when adpating it for survival analysis. Fig 4 is a visualization of the correlation matrix of clinical covariates in HNSCC dataset generated by PySurvival (Fotso et al. 2019–). As we can see in the figure, many of the covariates are heavily correlated with each other, posing a need of selecting useful features for constructing an efficient solution. According to a single value decomposition computation, the matrix is rank deficient. Despite that, FastCPH-LassoNet can successfully select a subset of covariates (11.04 out of 107) and attain a good and stable performance as shown in Table 2.