High-Dimensional Survival Analysis: Methods and Applications

3.1. Regularized Cox Models

The approach that dominates survival analysis in the biomedical literature is the Cox (1972) proportional hazards model, famed for presenting both a novel hazard model and a novel concept of partial likelihood. The model links the conditional hazard function (1) to X_i via

where the baseline hazard, λ₀(t), is unspecified, and β = (β₁, …, β_p)^⊤ is the coefficient vector of X_i to be estimated, with a fixed p < n, by maximizing the partial likelihood, i.e.,

with PL_i(β) being the contribution for subject i who is observed to die:

where ℛ(Y_i) = {j: Y_j ≥ Y_i}. In high-dimensional settings, that is, p > n and asymptotically p and n may both go to infinity (Zhao & Yu 2006), directly optimizing the partial likelihood is not feasible because of over-parameterization. Instead, regularized regression adds a penalty term to the negative log partial likelihood, ℓ(β), and optimizes a penalized version of objective function:

where the penalty Pen(β) is controlled by a positive tuning parameter, η, to be selected through cross-validation. A widely recognized family of penalties is based on the l_q-norm,

Regularization approaches with $\mathit{\text{Pen}}(\beta )=\parallel \beta {\parallel }_2^2$, known as ridge regression (Hoerl & Kennard 1970), are applied to the Cox model, and return unique and shrunk estimates (Verweij & Van Houwelingen 1994). However, ridge regression does not promote sparsity, as it cannot shrink individual coefficients to zero. The least absolute shrinkage and selection operator (LASSO) (Tibshirani 1996), with Pen(β) = ∥β∥₁, penalizes the absolute sum of the coefficient estimates and has been routinely used for producing sparse models. Its application (Tibshirani 1997) to survival settings, namely, Cox LASSO, has become a widely used approach for high-dimensional survival analysis by performing feature selection and estimation simultaneously (Figure 2).

LASSO has several notable statistical properties. It possesses model selection consistency under certain regularity conditions, in particular, the strong irrepresentable condition when p grows much faster than n (i.e., that the absolute sum of coefficients for the regression of any noise variable on signal variables must be strictly smaller than 1) (Zhao & Yu 2006). It has a Bayesian interpretation by viewing β as having a “double exponential” prior (Tibshirani 2009). However, as the LASSO penalty term is linear in the size of the coefficients, it leads to biased estimates, especially for the coefficients with large absolute values. To remedy this, Zhang & Lu (2007) proposed the adaptive Cox LASSO by utilizing Pen(β) = ∑_j w_j|β_j|, with smaller weights, w_j, assigned to larger coefficients and vice versa.

The estimates are $\sqrt{n}$ consistent if $\sqrt{n}\eta =O(1)$ and have oracle properties if $\sqrt{n}\eta \to 0$ and nη → ∞. When p > n, it was suggested to use robust estimates such as ridge regression estimates to determine the w_j’s.

Fan & Li (2002) proposed a smoothly clipped absolute deviation (SCAD) penalty, which is a quadratic spline function of |β| with knots at η and αη. Its derivative w.r.t |β|, i.e.,

may more clearly show the role of the penalty in regularizing estimating equations (Fan & Li 2002). While the SCAD penalty retains the penalization rate of LASSO for small coefficients, it relaxes the rate of penalization smoothly as the absolute value of the coefficient increases. Asymptotically, the SCAD penalty yields $\sqrt{n}$-consistent estimates (with a proper rate of η) and possesses oracle properties (if $\sqrt{n}\eta \to 0$ and nη → ∞). Strong oracle properties for LASSO and SCAD were established in Bradic et al. (2011), which further proposed a class of nonconvex penalization procedures for the Cox model. Nonconvex regularization, including SCAD, is appealing as it obtains support recovery properties under much weaker assumptions than for l₁ penalization (Loh & Wainwright 2017).

Another extension is the elastic net penalty for Cox models (Wu 2012), which combines the LASSO and ridge penalties; but, unlike LASSO, is capable of selecting more predictors than the sample size (Zou & Hastie 2005). This notion was generalized by Vinzamuri & Reddy (2013) with the kernel elastic net Cox regression model, which replaces the ridge penalty with β^⊤Σβ. Here, Σ is a p × p radial basis function kernel matrix of predictors, which measures pairwise similarity between predictors. This penalty is meant to encourage correlated predictors to have similar strengths on survival prediction. Other regularized Cox methods include the group Cox LASSO (Kim et al. 2012), which selects groups of related covariates as a whole, and the fused LASSO (Chaturvedi et al. 2014), which penalizes both the coefficient estimates and their successive differences for ordered features (Table 1).

3.2. The Dantzig Selector for Survival Data

Candès & Tao (2007) proposed another type of regularized estimator known as the Dantzig selector for linear regression:

where Y, X, β, and ϵ are an n × 1 vector of responses, an n × p covariate matrix, a p × 1 vector of coefficients and an n × 1 vector of zero-mean residual errors, respectively. It estimates β by solving

where η_Q > 0 is a tuning parameter. Empowered by linear programming, the Dantzig selector offers a useful alternative as a regularized estimating equation approach. As a dual problem of LASSO, it often produces the same solution path (Candès & Tao 2007).

On the other hand, accelerated failure time (AFT) models have become a useful alternative to Cox models due to the ease of interpretation (Saikia & Barman 2017). An AFT model links the (log transformed) survival time to covariates via a linear model

where the log transformation ensures the parameter space of β is unconstrained, and the distribution of the errors, e_i, induces a distribution for T_i (Table 2). For parametric AFT models, maximum likelihood estimation (MLE) can be used for inference. When e_i’s distribution is unspecified, the models are semi-parametric and the MLE is difficult to obtain, as the likelihood involves infinite dimensional parameters. With a fixed p < n, Buckley & James (1979) proposed an estimating equation approach by imputing the censored outcomes, and solving a least squares problem.

For AFT models with high-dimensional predictors, one cannot directly apply LASSO estimation, as the objective function again involves infinite dimensional parameters. Motivated by the work of Candès & Tao (2007) for regularized least squares estimation, Li et al. (2014) used Buckley-James imputation to express AFT estimation as a least squares problem and then applied the Dantzig selector:

where $T_i^{\ast }(\beta )=\text{log}\left(Y_i\right)+\left(1-{\delta }_i\right)\frac{{\int }_{e_i(\mathit{\beta })}^{\infty }\widehat{S}(s,\beta )ds}{\widehat{S}\left(e_i(\beta ),\beta \right)}$ are imputed outcomes with the Kaplan-Meier estimate, Ŝ(·, β), based on {e_i(β) = log(Y_i) – β^⊤X_i, δ_i}, i = 1, …, n, and ${\mathit{T}}^{\ast }(\beta )={\left(T_1^{\ast }(\beta ),\dots ,T_n^{\ast }(\beta )\right)}^{\top }$. Here, the projection matrix P_n = I_n – 11^⊤/n, where I_n and 1 are an n × n identity matrix and an n × 1 vector of 1’s, is for centering covariates to avoid estimation of the intercept or the expection of e_i in model (3). An iterative approach is necessary because this is not a linear programming problem, and, like the Dantzig selector for linear models, estimates may be biased and may not possess the oracle property.

To address this, Li et al. (2014) considered an adaptive version of Dantzig selector with data-driven weights that vary inversely with the magnitude of coefficients. They showed that the weighted Dantzig selector has model selection consistency and oracle properties. On the other hand, a Dantzig selector for the Cox model was proposed by Antoniadis et al. (2010) based on partial likelihood score equations.

Note that, in ultra high-dimensional settings where p ≫ n, penalized variable selection methods such as those described so far may incur high computational costs, numeric instability, and poor reproducibility (Fan & Lv 2008). As such, variable screening is a crucial first step in identifying predictive biomarkers and reducing the dimensionality of the feature space before applying regularized methods. Feature screening methods such as sure independence screen (SIS) (Fan & Lv 2008) fit marginal regression models for each covariate one at a time, choose a threshold, and retain those covariates with magnitudes of marginal effects above the threshold. In the ultra high-dimensional survival settings, additional censoring issues need to be addressed. Recent advancements in survival feature screening have included sure independence screening (Fan et al. 2010), principled sure independent screening (Zhao & Li 2012), score test screening (Zhao & Li 2014), concordance-measure based screening (Ma et al. 2017), Buckley-James assisted sure screening (Liu et al. 2020), conditional screening (Kang et al. 2017, Hong et al. 2018b), integrated power density screening (Hong et al. 2018a), L_q norm screening (Hong et al. 2020), and forward regression (Hong et al. 2019, Pijyan et al. 2020). See a focused review of survival feature screening in Hong & Li (2017).

3.3. Inference with High-Dimensional Covariates

As simultaneous estimation and inference is challenging within the high-dimensional survival framework, we review limited methods available for drawing inference in this area. More broadly, high-dimensional regression inference methods largely fall under post-selection inference and debiased LASSO estimation. Further challenges arise in that post-selection inference is conditional on the selected subset and does not account for variation in model selection. Several authors used debiased LASSO (Van de Geer et al. 2014, Javanmard & Montanari 2014, Yu et al. 2018) to draw inference; however, these methods require estimation of the inverse of a p × p information matrix, which is a daunting task, especially when p > n (Xia et al. 2021a,b).

4.1. Support Vector Machines

Support vector machines (SVMs) fall under the supervised learning family (Vapnik et al. 1995, Noble 2006) and seek to find a hyperplane that provides maximal separation between groups (Figure 4). Specifically, consider a binary outcome Y_i ∈ {−1, 1} for each individual i with a corresponding p-dimensional covariate vector, X_i. The goal of SVM is to identify a hyperplane, $H(\psi ,a)=\left\{v\in {ℝ}^p\mid \langle \psi ,v\rangle +a=0\right\}$, separating these two groups so that the margin, 2/∥ψ∥, can be maximized, where $\psi \in {ℝ}^p$ is the slope vector, and 〈·,·〉 denotes the inner product (Figure 4). Often, the two classes may not be separable in the original feature space within ${ℝ}^p$, and we use F(·) to map the original predictors to a higher dimensional space where the outcomes can be distinguished, in which case, the hyperplane to deal with is $H(\psi ,a)=\left\{v\in {ℝ}^p\mid \langle \psi ,F(v)\rangle +a=0\right\}$ and, with slight overuse of notation, the dimension of ψ is the same as that of F(v). In practice, F(·) does not have to be obtained explicitly and 〈ψ, F(v)〉 can be calculated by using a reproducing kernel (Wahba et al. 1999). We further introduce a slack variable, ξ_i = [1 − Y_i{〈ψ, F(X_i)〉 + a}]₊, to dictate the degree to which the ith data point is misclassified, as illustrated in Figure 4.

SVMs have been extended to model continuous time-to-event data, which are prone to censoring, by predicting the survival time to be 〈ψ, F(X_i)〉 + a. Van Belle et al. (2007) formulated the survival SVM based on the rank concordance between the prediction and observed survival time, Y_i, among comparable individuals in the presence of censoring. Specifically, they introduced a comparability indicator, $v_{ij}={\delta }_i\mathbb{I}(Y_i<Y_j)$, such that the ordering of the observed survival times for subjects i, j can only be determined when v_ij = 1. For a comparable pair with v_ij = 1, a concordance in rank is reached if and only if 〈ψ, F(X_j)〉 − 〈ψ, F(X_i)〉 > 0. Allowing varying degrees of pairwise slacks, i.e., when 〈ψ, F(X_j)〉 − 〈ψ, F(X_i)〉 ≤ 0 with v_ij = 1, across comparable pairs, Van Belle et al. (2007) proposed to solve

where ξ_ij’s are pair-specific slacks, whose summation is to be minimized, and γ > 0 is a regularization parameter controlling the maximal margin and misclassification penalties. This formulation can be shown to maximize the Harrell rank-based concordance index (C-index) (Harrell et al. 1982). Hence, it is termed the rank-based SVM approach for survival data and does not estimate the “intercept” a. An alternative regression approach (Shivaswamy et al. 2007) aimed to find a prediction, 〈ψ, F(X_i)〉+a, for continuous survival times, by identifying a hyperplane that best fit the data that are subject to censoring (Smola & Schölkopf 2004):

With censoring indicators incorporated into the constraints, the formulation utilizes available information from both censored and non-censored observations. To make full use of the strengths of both approaches, Van Belle et al. (2011) and Pölsterl et al. (2015) further proposed hybrid approaches, combining the penalties imposed by both methods.

4.2. Tree-Based Methods

While SVMs are adept at estimating non-linear relationships, they do not scale well for large datasets and often under-perform when the outcomes are noisy. Also there may be no clear interpretations for classifying data points above or below the estimated hyperplane (Somvanshi et al. 2016). Decision trees are an alternative for classifying patients that provide an intuitive interpretation of the hierarchical relationships between predictors. Broadly, classification and regression trees (CART) is an umbrella term for a set of recursive partitioning algorithms, which predict the group membership (classification) or target value (regression) for an observation based on a set of binary decision rules (Figure 5).

Gordon & Olshen (1985) first presented survival trees, and Ciampi et al. (1986, 1987) solidified the notion and established splitting criteria based on the log-rank and likelihood ratio test statistics, respectively, gaining predictive accuracy and interpretability. A recursive partitioning algorithm for generating a survival tree is given as follows.

1. Discretize each covariate to be a binary variable (categorical variables with m levels are expressed as m − 1 dummy variables).

2. For every binary covariate, X_j, j = 1, …, p, compute the log-rank statistic to test the difference between the survival curves for the two groups defined by X_j.

3. Choose the covariate, say, X_j∗, with the largest significant test statistic and partition the full sample (i.e., the root node) into two groups (child nodes) based on X_j∗.

4. Repeat steps 2–3 for each subset (child node) until reaching the terminal nodes, that is, no covariates produce a significant test statistic and there are enough events (exceeding a prespecified number) in each terminal node.

The resulting terminal nodes split the original sample into distinct groups, who are deemed more homogeneous within each group, and will output survival estimates via Kaplan-Meier estimation in each group. Further variations in splitting are based on metrics that accommodate censored data and by either minimizing within-node homogeneity or maximizing between-node heterogeneity. For example, these metrics can be Martingale residuals (Therneau et al. 1990) or deviance residuals (LeBlanc & Crowley 1992). With an established splitting criterion, to select a final tree, either a full survival tree is ‘grown’ and ‘pruned’ or a stopping rule is applied in backward or forward selection (Bou-Hamad et al. 2011).

4.3. Ensemble Learners

While survival trees provide a fast and intuitive means of studying hierarchical relationships of predictors with outcomes, they are prone to over-fitting and high variability (Hu & Steingrimsson 2018, Steingrimsson et al. 2016). Ensemble learners overcome the instability issues by using techniques such as bagging, boosting, and random forests.

4.4. Deep Learning and Artificial Neural Networks

Deep learning has emerged as a powerful tool for risk prediction. This work stems from artificial neural networks that tried to mirror how the human brain functions (Rosenblatt 1958), wherein nodes (or neurons) are connected in a network as a weighted sum of inputs through a series of affine transformations and non-linear activations.

A fully-connected, feed-forward artificial neural network is made up of L layers, with k_l neurons in the lth layer (l = 1, …, L) (Figure 7). With an input, network predictions are made based on an L-fold composite function, f_L ◦ f_L−1 ◦ ⋯ ◦ f₁(·) with (g ◦ f)(·) = g(f(·)). At the lh layer, f_l(·), is defined as

where v is a k_l−1 × 1 input vector fed from the (l − 1)th layer, ${\sigma }_l(\cdot ):{ℝ}^{k_l}\to {ℝ}^{k_l}$ is an activation function, W_l is a k_l × k_l−1 weight matrix, b_l is a k_l × 1 bias vector, and the 0th layer is the input layer. Typical choices of σ_l(·) include the sigmoid function or the rectified linear unit activation function (ReLU), that is, σ_l(b) = max(0, b), where $b\in {ℝ}^{k_l}$ and max(0, ·) operates component-wise.

For survival prediction, several deep learning approaches have emerged, beginning with the seminal work of Faraggi & Simon (1995), which adopted a fully-connected, feed-forward neural network to extend the Cox model to nonlinear predictions. Other feed-forward neural networks (Liestbl et al. 1994, Brown et al. 1997, Biganzoli et al. 1998, Eleuteri et al. 2003) used the survival status as a training label, and output predicted survival probabilities. Further developments have been made in Bayesian networks (Bellazzi & Zupan 2008, Lisboa et al. 2003, Fard et al. 2016), convolutional neural networks (Yao et al. 2017, Katzman et al. 2016, 2018, Ranganath et al. 2016), and recurrent neural networks (Yang et al. 2018).

5.1. Competing Risks

In a competing risk setting, observing an event type, labeled by c ∈ {1, …, K}, effectively eliminates the chance of observing other event types happening to the same individual (Young et al. 2020). For example, when studying the survival of patients with cancer, competing events can be cancer-related death (c = 1) or death by cardiac disease (c = 2) (Figure 8); an individual cannot die of cardiac disease once they have died of cancer, and vice versa. For characterizing the risk of competing events, there are two commonly used statistical metrics, namely, the cause-specific hazard and the subdistribution hazard, which target different counterfactual scenarios. The former describes the risk under hypothetical elimination of competing events, while the latter is about the observable risk without elimination of any competing events (Rudolph et al. 2020).

Several authors (Lau et al. 2009, Koller et al. 2012) have stated that the subdistribution hazard is useful for predicting the probability of having an event of a type of interest by a given time, termed the cumulative incidence function (CIF), which reflects an individual’s actual risks and prognosis. In the following, we focus on the subdistribution hazard, which is derived from CIF, i.e., F_c(t) = Pr(T_i < t, 𝒞_i = c), where 𝒞_i marks the event type for subject i. Specifically, for each event type c = 1, …, K, it is defined as

which denotes the instantaneous risk of failure from event type c among those who have not experienced this type of event. That is, the risk set at t includes those who are event free as well as those who have experienced a competing event (other than type c) by t. The subdistribution hazard model (Fine & Gray 1999) links a subdistribution hazard function to covariates via

where λ_0c(t) is the baseline subdistribution hazard function for event type c, and β specifies the effect of X_i on the probability of event c occurring over time. In fact, model (7) implies that $1-F_c\left(t\mid X_i\right)={\left\{1-F_{0c}(t)\right\}}^{\text{exp}\left(X_i^{\top }\beta \right)}$, where F_c(t|X_i) and F_0c(t) are the CIF given X_i and the baseline CIF, respectively.

With high-dimensional predictors, several authors (Kawaguchi et al. 2019, Ha et al. 2014, Ahn et al. 2018) proposed regularized subdistribution hazard models for variable selection, and Hou et al. (2019) further performed inference using a one-step debiased LASSO estimator. For prediction, several deep learning works for competing risks have been proposed based on CIFs. For example, DeepHit (Lee et al. 2018) developed a multi-task network to nonparametrically estimate F_c(t|X_i) for c = 1, …, K. The network is trained to minimize a loss function, which is constructed based on the joint distribution of the first hitting time for competing events of each subject, while ensuring the concordance of estimates across subjects (Harrell et al. 1982), that is, a patient who died at a given time should have a higher risk at that time than a patient who survived longer. Dynamic DeepHit (Lee et al. 2019) further incorporated longitudinal information for dynamic predictions. Other approaches have included DeepCompete (Aastha & Liu 2020), as well as a hierarchical approach to multi-event survival analysis (Tjandra et al. 2021).

5.2. Semi-Competing Risks

Semi-competing risk problems, a variant of competing risk problems, have commonly been encountered in clinical studies. By semi-competing, we mean that the occurrence of one event, i.e., a non-terminal event, is subject to the occurrence of another terminal event, but not vice versa (Figure 9). As the non-terminal event (e.g., cancer progression) is often a strong precursor to the terminal event (death), semi-competing events are often related and, hence, the terminal event may informatively censor the non-terminal event (Jazić et al. 2016). To overcome such informative censoring, researchers either consider only the terminal event (i.e., mortality) or a composite outcome such as progression-free survival, that is, time to progression or death, whichever comes first.

What is lacking here is how to model a predictor’s potentially different roles in disease progress and death, while utilizing the crucial information about the sojourn time between progression and death. Even in settings where the non-terminal and terminal event times are only modestly correlated, failing to acknowledge this sojourn time may lead to incorrect inference or biased predictions (Crilly et al. 2021).