Survival Analysis on Rare Events Using Group-Regularized

9 We propose a Sparse-Group regularized Cox regression method to improve the prediction 10 performance of large-scale and high-dimensional survival data with few observed events. Our 11 approach is applicable when there is one or more other survival responses that 1. has a large 12 number of observed events; 2. share a common set of associated predictors with the rare event 13 response. This scenario is common in the UK Biobank (Sudlow et al. 2015) dataset where 14 records for a large number of common and rare diseases of the same set of individuals are avail15 able. By analyzing these responses together, we hope to achieve higher prediction performance 16 than when they are analyzed individually. To make this approach practical for large-scale data, 17 we developed an accelerated proximal gradient optimization algorithm as well as a screening 18 procedure inspired by Qian et al. (2019). We provide a software implementation of the proposed 19 method and demonstrate its efficacy through simulations and applications to UK Biobank data. 20 21 Cox proportional hazard model; Sparse-Group Lasso; Multi-response regression. 22

where β ∈ R d is the coefficients vector that measures the strength of association between X and the response. This hazard function is equivalent to the cumulative distribution function: exp(X T i β k ) j:Tj ≥Ti exp(X T j β) When the number of observed events is small relative to n, estimating β becomes challenging. This 29 could happen, for example, when the time-to-event response is the age of diagnosis of a rare disease. 30 In particular, if O i are i.i.d Bernoulli random variables with probability p, then the information 31 matrix is proportional to p and thus the asymptotic variance of the maximum partial likelihood 32 estimate is inversely proportional to p.

34
We evaluate a fitted survival model using the concordance index, or the C-index. For a param-35 eter estimateβ, the C-index is defined as (1) The Lasso method (Tibshirani 1996) makes the assumption that only a small subset of predictors are 39 associated with the response. In other words, it assumes that β has only a small number of non-zero 40 entries. A sparse solution can be obtained by optimizing an L 1 -regularized objective function.  In this section we define the notations and the key model assumptions that we will use in the 50 subsequent sections. For an integer n, define [n] = {1, 2, · · · , n} and define x + = max{x, 0} for all 51 x ∈ R.

52
We analyze K ≥ 1 time-to-event responses on n individuals. For example, the responses could be time from birth to K different diseases. The data we observed are in the format: R n×d . For k = 1, · · · , K, O k i = 1 implies that T k i is the true time until the event k and for the ith 54 individual, and O k i = 0 implies that the true time until the event k is right-censored at T k i . We 55 assume each response follows a Cox proportional hazard model: where h 0,k : R + → R + is the baseline hazard function of the kth response. Let be the number of observed event k. 57 We make the assumption that not only β k is sparse for all k ∈ [K], but they also have a small } is a small set relative to d. In human 59 genetics applications, the first assumption means that each response is associated with a small set 60 of genetic variants, and the second assumption implies that there are significant overlap among the 61 genetic variants that are associated with each response. This belief is the main driver for prediction 62 performance improvements on rare diseases, and it translates to the following regularized partial 63 likelihood objective function: Here the first term is the sum of the K negative log-partial-likelihood, normalized by the number of observed events. The second term is the regularization. β j 1 , β j 2 are the 1-norm and 2norm of the coefficients for the jth variable. That is, if we put all the parameters into a matrix B = [β 1 , β 2 , · · · , β K ] ∈ R d×K , then β j is the jth row of B and β k is the kth column. Note that when α = 0, the objective function decouples for each β k and they can be optimized separately. In our implementation we solve a slightly more general problem: where {w j } d j=1 are user provided penalty factor for each variables, which may be useful in the set-  To simplify the notation, for j ∈ [d] denote Here g j ∈ R K is the partial derivative of the smooth part of (4) with respect to the coefficients of the variable j. Finally, let S 1 (·; λ) : R K → R K be the element-wise soft-thresholding function: particularly suitable when the objective function is the sum of a smooth function and a simple 75 function. In our case the smooth function is the sum of the negative log-partial-likelihood functions, 76 and the simple function is the regularization term. The algorithm alternates between two steps until 77 convergence criteria is met: 78 1. A gradient step that decreases the smooth part of the objective: where t is the step size that we determine using backtracking line search. 79 2. A proximal step that keeps the regularization term small: Here the proximal operator prox t : To simplify the notation we omit the dependency of prox t (x) on λ, α. Simple calculation shows that where S 2 (·; tαλ) : R K → R K is a group soft-thresholding operator: We describe the optimization algorithm in pseudocode, including details about Nesterov acceleration  For each We use the following results to get an optimality condition for the solutions β 1 , · · · , β K . The proofs 96 are given in section 6.

97
Algorithm 1: Proximal Gradient Method for (4) Set line search parameter η > 1; Denote the parameter matrix B ∈ R d×K . Initialize B (0) ; Write the smooth part of the objective as: Set iteration count i = 0; Set initial step-size t = 1; Set Nesterov weights w 0 , w 1 = 1; while B has not converged do Nesterov acceleration: Check convergence based on objective value change or parameter change.
This result motivates us to first fit a model (solving (4)) using a small number of variables whose 100 gradient has the largest α * norm, assuming the coefficients for the rest of the variables are all zero.

101
Then to verify the validity of the assumption we check that g j α * ≤ λ for variables assumed to 102 have zero coefficients. We refer to this step as KKT checking. Note that based on its definition (6), 103 it's not clear how we can compute v α * . Here we give a more explicit characterization.

105
Using the above we can check if v α * ≤ λ quite easily and compute v α * using binary search 107 108 Now we are ready to state the overall structure of our algorithm with variable screening. Sup-109 pose valid solutions for λ 1 , · · · , λ l have been obtained. Next we follow these steps:

113
• Top M variables with the largest g j α * that are also not ever-active. 114 We denote the set of variables used to fit (4) as the strong set S ⊆ [d]. assuming β j = 0 for all j ∈ S. This is done using proximal gradient descent (algorithm 1).

117
To speed-up the computation we initialize the variables at the previous valid solution (warm 118 start).   We implemented the proximal gradient method in section 2.1 as an R package, available at https: 132 //github.com/rivas-lab/multisnpnet-Cox. In this package we also implement the screening pro-133 cedure described in this section for genetics data in Plink2 format (Chang et al. 2014).
M . This will be the variables used in the fitting step; Fitting: for λ ∈ Λ (i) do Obtain parameter estimates by solving (4) using only variables in S (i) ; The optimization algorithm is described in algorithm 2; Coefficients for variables not in S (i) are set to be zero.; end Checking: Compute the full gradients g 1 (λ), · · · , g d (λ) at solutions obtained with regularization parameters λ ∈ Λ (i) ; Find the smallest λ ∈ Λ (k) such that the KKT condition (7) is satisfied: which satisfies proportional hazard. The parameters β 1 , β 2 have a common support of size 35. We 148 use a large (5000 samples) and uncensored validation set to select the optimal λ. α is fixed at √ 2.

149
The parameter estimate corresponding to the best λ is then evaluated at a large, uncensored test The assumption that the coefficients for different responses have the same support can come from 156 domain knowledge (such as biology). In practice this is usually not exactly satisfied. Here we use 157 simulation to examine the robustness of our approach. The setup is the same as the previous sim-  2} using a small validation set of size 200 with 5% uncensored second event. The whiskers indicates 95% confidence intervals. 161 We can see that, when the support overlap percentage is less than around 40% the prediction 162 performance actually becomes worse when the two responses are trained together. One solution to 163 this problem is to also treat α as a hyperparameter and use the validation set to determine it. For In this section we apply the proposed method to UK Biobank data. We focus on thyroiditis, which 172 has 808 observed events in the study population (337, 129 white British participants). We randomly 173 assign 70% of the samples to the training set, 10% to the validation set, and 20% to the test set. The when we use the validation set to determine both the optimal λ and whether to set α = 0 or α = √ K.  Table 1: Test C-index of thyroiditis when this response is paired with other ones. The second column is obtained when α is fixed at √ K (which is 1 in the first row, √ 7 in the last row, and √ 2 in the rest). The third column is obtained when we use the validation set to choose α between 0 (single-response) and √ K. Improved C-indices are given in boldface. Baseline test C-index of thyroiditis, when only the 11 covariates are used for fitting, is 0.649.

188
We developed a regularized regression method for multiple survival responses. This method is par-189 ticularly suitable when we can pair a response with few observed events with others having a larger 190 number of events, such that these responses have a same set of useful predictors. We demonstrate 191 the improvements in the prediction accuracy for the rare events responses through simulation and 192 real data applications. We also provide efficient implementation for the proposed method.

194
Here are two directions for future studies. When there are more than two responses, the relationship 195 between the prediction performance and the degree of overlapping in the coefficients support is yet 196 to be understood. On the practical side, it is reasonable to have different types of models (not 197 just Cox model), or even different response types (such as binary or count) to boost the accuracy 198 of survival analysis on rare events. In principle the same type of regularization could still be applied.

200
For the GPU implementation, one challenge is the limited amount of memory available on GPUs.

201
In modern computer clusters it is common to have machines with hundreds of Gigabytes of system 202 memory, but most GPUs have less than 32GB of memory. In our implementation we use single pre-6 Proof of Propositions 211 6.1 Proof of Proposition 1 212 We show a slightly more general result from convex analysis. Let f : R K → R be continuously 213 differentiable, convex, and bounded from below. Let · be a norm on R K , and · * be its 214 corresponding dual norm. Let λ > 0, and set 215 x * := arg min We will show that It is clear that (8) is equivalent to the constrained optimization problem 217 min x,y f (x) + λ y such that x = y.
The Lagrangian dual is 219 g(z) := inf x,y L(x, y, z) = inf y λ y + z T y + inf Using the definition of dual norm, when z * > λ, the infimum of the first term above is −∞, and 220 when z * ≤ λ, the infimum is 0. Therefore Therefore the dual solution z * := arg max z g(z) must satisfy z * ≤ λ. Now we go back to the 222 Lagrangian. Since the primal objective is convex and Slaters condition holds, the solution to the 223 primal problem can be obtained through minimizing 224 L(x, y, z * ) = f (x) + λ y + z * T (y − x).
That is v α * ≤ 1 if and only if v ∈ B * α , which by definition means that v = v 1 + v 2 for some