Hazard Gradient Penalty for Survival Analysis
Authors: Seungjae Jung
Kyung-Min Kim
Venue: Submitted to ICLR 2023 (2023)
Keywords: survival analysis
gradient penalty
KL divergence

Abstract: Survival analysis appears in various fields such as medicine, economics, engineering, and business.
Recent studies showed that the Ordinary Differential Equation (ODE) modeling framework integrates many existing survival models while the framework is flexible and widely applicable.
However, naively applying the ODE framework to survival analysis problems may model fiercely changing density function with respect to covariates which may worsen the model’s performance.
Though we can apply L1 or L2 regularizers to the ODE model, their effect on the ODE modeling framework is barely known.
In this paper, we propose hazard gradient penalty (HGP) to enhance the performance of a survival analysis model.
Our method imposes constraints on local data points by regularizing the gradient of hazard function with respect to the data point.
Our method applies to any survival analysis model including the ODE modeling framework and is easy to implement.
We theoretically show that our method is related to minimizing the KL divergence between the density function at a data point and that of the neighborhood points. 
Experimental results on three public benchmarks show that our approach outperforms other regularization methods. 

Under review as a conference paper at ICLR 2023
HAZARD GRADIENT PENALTY FOR SURVIVAL ANALYSIS
Anonymous authors Paper under double-blind review
ABSTRACT
Survival analysis appears in various fields such as medicine, economics, engineer- ing, and business. Recent studies showed that the Ordinary Differential Equation (ODE) modeling framework integrates many existing survival models while the framework is flexible and widely applicable. However, naively applying the ODE framework survival analysis problems may model fiercely changing density function with respect to covariates which may worsen the model's 's performance. Though we can apply 1.1 or 1.2 regularizers to the ODE model, their effect on the ODE modeling framework is barely known. In this paper, we propose hazard gradient penalty (HGP) to enhance the performance of a survival analysis model. Our method imposes constraints ent of hazard function with respect to the data point. Our method applies to any survival analysis model including the ODE modeling framework and is easy implement. We theoretically show that our method is related to minimizing the KL. divergence between the density function at a data point and that of the neigh- borhood points. Experimental results on three public benchmarks show that our approach outperforms other regularization methods.
local data points by regularizing the gradi-
INTRODUCTION
Survival analysis (a.k.a time-to-event modeling) is a branch of statistics that predicts the duration of time until an event occurs (Kleinbaum & Klein 2012). Survival analysis appears in various fields such as medicine (Schwab et 2001), economics (Meyer 1988), engineering (O'Connor & Kleyner 201T), and business (Jing & Smola 2017, Li et al. 2021). Due to the presence of right- censored data, which is data whose event has not occurred yet, survival analysis models require special considerations. Cox proportional hazard model (CoxPH) (Cox 1972; Katzman et al.) 2018) and accelerated time failure model (AFT) (Wei 1992) are widely used to handle right-censored data. Yet the assumptions made by these models are frequently violated in the real world (Lee et al. 2018 Tang et al. 2022aa.. Recent studies showed that the Ordinary Differential Equation (ODE) modeling framework integrates many existing survival analysis models including CoxPH and AFT (Groha
al. 2020 Tang et al., 2022a b). They also showed that the ODE modeling framework is flexible and widely applicable.
However, naively applying the ODE framework to survival analysis problems may result in wildly oscillating density function that may worsen the model's performance. Regularization techniques that can regularize this undesirable behavior are understudied. Though applying L1 or L.2 regulariz- ers to the ODE model is one option, their effects on the ODE modeling framework are barely known. The cluster assumption from semi-supervised learning states that the decision boundaries should not cross high-density regions (Chapelle et al. 2006). Likewise, survival analysis models need hazard functions that slowly change in high-density regions.
Suppose we attempt predict the time to death of three individuals A, B, and C. Assume the traits of A and B are similar and the traits of B and C are dissimilar. is natural to expect that the probability distribution of time-to-death of A should be close to that of B while far from that of C. The expectation aligns with the cluster assumption. Explicitly modeling the assumption enhances the performance long it holds.
In this paper, we propose hazard gradient penalty to make a slowly changing (with respect to co- variates) survival analysis model in high-density regions. In nutshell, the hazard gradient penalty

---

Under review as a conference paper at ICLR 2023
regularizes the gradient of the hazard function with respect to the data point from the real data dis- tribution. Our method has several advantages. 1) The method is computationally efficient. 2) The method is theoretically sound. 3) The method is applicable to any survival analysis model including the ODE modeling framework as long as it models hazard function. 4) It easy to implement. We theoretically show that method is related to minimizing the KI. divergence between the density function at a data point and that of the neighborhood points of the data point.
Experimental results on three public benchmarks show that our approach outperforms other regular- ization methods.
PRELIMINARIES
Survival analysis data comprises of an observed covariate æ, failure event time and an event indicator e. If event is observed, corresponds to the duration time from the beginning of the follow-up of an individual until the event occurs. In this case, the event indicator e = 1. If an event unobserved, corresponds to the duration time from the beginning of follow-up of an individual until the last follow-up. In this case, we cannot know the exact time of the event occur and event indicator 0. An individual is said be right-censored if 0. The presence of right-censored data differentiates survival analysis from regression problems. In this paper, we only focus on the single-risk problem where event is binary-valued variable.
Given a set of triplet D = { 44, e;)] the goal of survival analysis is to predict the likelihood of an event occur p(t x) or the survival probability S(t x). The likelihood and the survival probability have the following relationship:
S(t x) = [ p(T x)
Modeling p(t or S(t |)) should satisfy the following constraints:
p(t x)) > 0, p(r x)d =1
S(0 x) = 1,
lim S(t x) = 0, S(t x) z S(12 x)) i S t2
t1000
Previous works instead modeled the hazard function (a.k.a conditional failure rate) la(t 1972
x) (Cox
Katzman et al. 2018: Wei. 1992 Zhong et al. 2021).
P At T2t,2) S(t a))
p(t x)
h(( x)) = lim
(2)
A1=0
At
As the hazard function is a probability per unit time, it is unbounded upwards. Hence, the only constraint of the hazard function is that the function is non-negative: hit a>) 2 0
2.1 THE ODE MODELING FRAMEWORK
We can obtain an ODE which explains the relationship between the hazard function and the survival function by putting derivative of equation Tinto equation (Kleinbaum Klein, 2012).
h(t x)) = p(t x)
(- dS(t x))
log S(t | x)) dt
(3)
5(t x
S(t
x))
dt
Starting from initial value log S(0 x) = O, we can define log S(t a) as the solution of the ODE initial value problem where the ODE is defined equation 3
log S(t x)) = log 5(0 x:) + -h((T w))d 1=)d===f -n-((T x:)
Tang et al. 2202bbss 's formulation is slightly different in that their hazard function also depends on the cumulative hazard. To our understanding, depending cumulative hazard is redundant so we conduct experi- ments without it.
2

---

Under review as a conference paper at ICLR 2023
We can train the ODE model by minimizing the negative oog-likeihhod.
- - log Po x) - e) log St { x)
(4)
= - (log ht x) iz) + log So =)) (1 - e) log St x)
Following Groha et al. (020)), update the model parameters using Neural ODEs (Chen al.
2018). The hazard function (t (t] æ) is modeled using a neural network followed by the softplus activation function ensure that the output is always non-negative.
2.2 NEURAL ODES
Neural ODEs model the continuous dynamics of variables (Chen et al. 2018). Starting from =(0,, we can define the output =TT) to be the solution of the following ordinary differential equation (ODE) initial value problem.
dz(() f()()),t,)) 2(T) =(0) + L/' f(z((),t,
Naively applying an ODE solver to an ODE initial value problem leads to practical difficulties. An ODE solver builds big computation graph which incurs high memory cost and additional numerical errors may occur in backpropagation steps. Chen et al. (2018) showed that we can obtain the gradients of rcalar-valued loss w.r.t all inputs of any ODE solver with constant memory cost. We can calculate the gradients without backpropagating through the operations of the solver but with another call to ODE solver.
3 METHODS
In this section, we introduce the hazard gradient penalty and show that is related to minimizing the KI. divergence between the density function at data point and that of its neighbours. See Figure for the graphical overview of our method.
Poft (x)
Polt (x)
DRL. (patt x))pp (x])
ZL
The cluster assumption from semi-supervised learn- ing states that the decision boundaries should not cross high-density regions (Chapelle et al.) (2006). In similar vein, hazard functions of survival analy- sis models should change slowly in high-density re- gions.
Em-
Valligtt x
M(x++(
Adt hitt |x)
Consider a case where two data points σ1, D2 E IRd failed 1,,42 (t1 12) each. Under the cluster assumption, a point æy E B( should fail at t/ fg. p(t <<) is skewed for some reason and puts high density at ti < f1, it worsens the model's performance. To evade such situation, $(t a²)) should not deviate too much from p(t 22). To achieve this, we propose the following regularizer
Figure 1: Graphical overview of proposed method. Our method min- imize the hazard gradient penalty (s,(t(a) h/ x)]12 and the nega- tive ligeliiel in equation 4 at the same time. Intuitively speaking, we regularize the model so that the hazard function does not vary much when a small noise is added or subtracted to the data point æ. section 3.2 show that minimizing the hazard gradient (5) penalty connected to minimizing the KL divergence between the density and the density at x' E B(x, e).
Rs = EE-SS(() hoa 2)
Bar, e) e-ball centered at 'In practice, we implement ho(t x)) using neural network whose input is a combination (concatenation, addition or both) of t and z. Hence we can write ho | æ.) and ho (t, æ)) interchangeably. The gradient V. (t, an) is naturally defined and is Vaho x).
3

---

Under review as a conference paper at ICLR 2023
3.1 EFFICIENT
SAMPLING FROM THE SURVIVAL DENSITY
The sampling operation ~ S(L a: in equation may induce computational overhead. To boost the sampling operation, use log So x:) which was computed during the negative log-likelihood calculation in equation Let [f...... tk] k be the union of the time points minibatch. The time points are sorted in increasing order. The adaptive time stepping in ODE solvers are sensitive to the time interval lx - t rather than the number of time points (Rubanova al. 2019). We can access log (t x) with negligible overhead long as t < to <tx.
We sample ts from a categorical distribution whose k-th weight defined as Sp(te x) exp(-_o' h(T | x))r). We finalize the sampling process by sampling t from the uniform dis- tribution (((t tx+1l). In this way, we don't have to calculate S(t x) again for sampling t. See Algorithm in Appendix for the pseudo code of ODE based survival analysis with the hazard gra- dient penalty.
3.2 CONNECTION TO KL DIVERGENCE
We now show that the hazard gradient penalty in equation S is equivalent to minimizing the ap- proximation of the upper bound of the KI. divergence between the density function at a data point and that of the neighborhood points of the data point. Henceforth, we denote X' by the subset of d-dimensional real space Idd
Theorem Suppose the hazard function is strictly positive function for all data point x, x' X. The KL divergence
Ep(t(a)) [log
(>) - log,((t x')]
is upper bounded by
Ep((()) |log h(t x)) - log h(t | x')ll2 + Es(ela)) |[h(t] (x) - / x')1
(6)
To prove Theorem I we need the following lemma.
Lemma The expectation of survival densities difference under the density is the negative of the expectation of hazard functions difference under the survival density. In other words,
Ep(t(a) [log S(t x)) - log S(t x')] = -Es(((a)) [h(t] x - hht xx)]
for all x, x' EX
Proof) We use the fact that Es(s))) (log S(t x)) - log S(t x²)) is constant with respect
d Es((a)) (log S(t) w) log S(1 ax)))
dt
dt / / S(t x)) (log S(t x)) log S(t x
d
= - / p(l (x) (log S(1 () - log S(1 x')) di + / S(t (-)((( () ((1(=()) di
Hence,
Ep(((e) [log S(L æ) log 5(1 xx)] - -Es(( [M(E a:) -1(1 xx))
*S(t ..) not valid prohability distribution as we cannot guarantee SSt | a=)dt = Rigorously, we
sample f a(t. m:) where s(t æ) = S(t <>/ J S(t) at We use S(t a)) for notational simplicity. See Appendix for the existence of so x))
4

---

Under review as a conference paper at ICLR 2023
We now go back to Theorem | and prove the theorem.
Pp((a)) [logr( m²) log (t] az)]]
[pp((()) [lt p( æ) - log x x')]]|2 (:: DxL 2 0)
[Bp(((()) [log htt x)) log h(t [z']] Bp(t(ae) [log S(t a)) log S(t a²)) lla (: equation 2
|[Ep((|=) [log (tt x)) log h(t [z']] + Bs(((a)) [h(t x) - h(t x²)] [l| (: Lemma Dp
< [ [log (t x)) - logh(t htt x x2)]ll2 + [| [h(t x) - h(t x a')lll2 (t. triangle inequality)
< Ep((()) [log (tt x) - log ( + (ss()))) [] m) - h(t x')ll2
Theorem An approximation the the upper bound of the KL divergence given in equation 6 is upper bounded by ((())( V| x)| l2 x' is in the epsilon ball centered at æ, i.e. x' E [B(x.c). c).
To prove the theorem, we first find the approximation.
Lemma 2 2e Vl hit a=)"(x² - 12)l2 2 is an approximation of the upper bound of the KL. divergence which is given equation
Ep(())) [llog h(t =) - log h(t x 1l1 Es(( [l 2:) - h(t x 22)122
Ep((/m) || V log h(t x)) 2)||2 + Es( mn) || Vih(t x)) (x' - 2>)12
(:: log h(t
x²) ~ log htt æ)) + V lot x² (a' - x)
and h(t x²) h(t x)) + vah(t -ex - x²)
V att( x²
Ep(((=)
li(t x)) (x' x))
+ Es(((a) Vah(t TT²(( =>)12
p(t
V/ah((() w²)² (a) - 20)122 + Es( V =0)((a" -x>===))-> x))12(( (: h(4 [x) > 0)
ht
x))
())())))
Vat((t) x) (x' - x))12
Obviously, 2[s(()))) Val(( " - x: 2 S ())(()))) max V ax)² )=) =))12- As
we assumed x' E B(x,,) in Theorem exx Vah(t x) 21)12 is achieved when
x' cVah( x) V h( 2) Hence,
sB)(( ha "x) - 2)112 S 2c Vah(( 2)112
and this concludes the proof.
Theoremp 2 shows that regularizating the hazard gradient penalty in equation 5 is equivalent to min-
imizing the approximation of the upper bound of the KL divergence Ep(t]a) [log p(t x) - log p(t a²)]- To incorporate the regularizer into the negative log-likelihood loss, we minimize the Lagrange multiplier defined as the sum of the negative log- ligeeiidl and the hazard gradient penalty regular- izer.
C E(a.t.))< [e + XR
(7)
Here, is a coefficient that balances the negative log- likelihoddo and the regularizer. See Appendix B for the code snippet of our JAX implementation (Bradbury et al. 2018).
Minimizing the hazard gradient penalty in equation 5 has two advantages over minimizing the KL divergence directly: a) computational efficiency and b) reduced burden of hyperparameter tuning. To compute the KI. divergence, we first sample xx' E B(xr,e) We then need to compute four values: (((]z) ((t]))),(((t)) and S(t)='). In this case, we have to compute hazard values of every t ~ S(f)2)) Further, we need one more hazard function integration S(t x') = exp(- f hi(t x')). On the other hand, regularizing the hazard gradient penalty only need to calculate the gradient of the hazard function.
When it comes to regularizing the KL divergence, we have to set the appropriate value of the reg- ularizing coefficient X' and the size of the ball e. On the other hand, if we regularize the hazard gradient penalty, we don't need to tune as À in equation incorporates e.
5

---

Under review as a conference paper at ICLR 2023
4 EXPERIMENTS
In this section, we experimentally show that the hazard gradient penalty outperforms other regular- izers. Further, check the hyperparameter sensitivity of hazard gradient penalty. Throughout the experiments, we use three public datasets: Study to Understand Prognoses Preferences Outcomes and Risks of Treatment (SUPPORT) the Molecular Taxonomy of Breast Cancer International Con- sortium (METABRIC) and the Rotterdam tumor bank and German Breast Cancer Study Group (RotGBSG) Table summarizes the statisties of the datasets. See Appendix A for evaluation metrics and experimental details.
4.1 METHODS COMPARED
We compare ODE HGP with four methods: vanilla ODE, ODE + LI, ODE + L2 ODE LCI. Vanilla ODE minimizes the expectation of the negative log- likeliiod in equation imizes the Lagrange multiplier defined as the sum of the expectation of the negative log- ligliiel and the LI penalty term: Eit,,,)))Do + a Ep-
ODE LI min-
|ump|
Here, wpS are model parameters and a is a coefficient that balances the negative log-likelihood and the I.1 penalty term. This is extension of Lasso-Cox [Tibshirani] 1997) to the ODE modeling framework. ODE + 1.2 minimizes the Lagrange multiplier defined as the sum of the expectation of the negative log- likeliikld and the L2 penalty term: E(a,,e)-DCa + >p ",
Here, ups are model parameters and a is a coefficient that balances the negative log-likelihood and the L2 penalty term. This is an extension of Ridge-Cox (Verweij & Van Houwelingen 1994) the ODE modeling framework. ODE + LCI minimizes the Lagrange multiplier defined as the sum of the expectation of the negative log- liglikklihood and the negative of the lower bound of a simplified version of time-dependent C-index. The regularizer is defined as
N Σ Ei (T) <; < () + (log o (t x.) < So x3))\ log 2)
M Σソ I < T;, <)
This is equivalent to time dependent concordance index in Section[A.1.I. if we don't 't take the Kaplan- Meier estimator into account. The regularizer is reminiscent of the lower bound of C-index (Steck al. 2007]. Although the lower bound of C-index was originally proposed as a substitute of the negative log-likelihood, [Chaffuwa et al. (2018) used the lower bound (Steck et al. (007) as a regu- larizer of the AFT model (Wei 1992).
4.2 RESULTS
Table shows the mC mAUC, and iNLLL scores The hazard gradient penalty outperforms other methods across almost all metrics and datasets. The interesting point is that both LI and L2 penalties do not affect the ODE model's performance in most cases. We speculate that regularizing the weight norm is effective in CoxPH as the model is simple and has a strong assumption that the hazard rate is constant. On the contrary, regularizing the norm of the weight may not be able to affect the ODE model's performance as ODE models are much more complex than CoxPH. Also, the experimental results highlight the possibility that the performance of the survival analysis models is more related to the local information such as the gradient at each data point rather than the global information such as the weight norm of the model. Figure shoos that the ODE HGP effectively regularized the variation of the density with respect to the input while other methods could not.
Table T also shows that regularizing the lower bound of the C-index not effective in many cases. We conjecture that the method is ineffective as the ODE modeling framework is flexible and optimiz- ing the negative log- likellhoooo can discriminate each data point's rank. Furthermore, regularizing the lower bound of the C-index does not harness the information of neighbors of data points. We also compare HGP against a Neural ODEs specific regularizer. See Appendix D for details.
"httss://ittuu.com/// thtpss//github.coom
See Appendix A.1 for the details of mC MUUC, and iN BL.L.

---

Under review as a conference paper at ICLR 2023
Table 1: Experimental Results on three datasets. Averages and standard deviations of different random seeds for each setting are shown. See Appendix A.2 for the details of experimental setups. The dagger mark indicates that the result is statistically significant (p 0.05) compared the result of vanilla ODE.
(a) mC (1)
Method
SUPPORT
METABRIC
RotGBSG
CoxPH
0.672 + 0.008
0.649 + 0.012
0.710 + 0.007
DeepHit
0.762 + 0.004
0.688 + 0.010
0.711 + 0.009
ODE
0.771 + 0.003
0.695 + 0.008
0.718 + 0.005
ODE + I.1
0.771 + 0.002
0.697 + 0.008
0.715 + 0.006
ODE L2
0.770 + 0.006
0.695 = 0.006
0.716 + 0.005
ODE + LCI
0.771 + 0.003
0.698 + 0.002
0.716 + 0.005
ODE HGP
0.775 0.004
0.702 + 0.009
0.723 + 0.006
(b) AUC((()
Method
SUPPORT
METABRIC
RotGBSG
CoxPH
0.706 + 0.010
0.685 0.008
0.739 1 0.009
DeepHit
0.800 = 0.004
0.720
0.012
0.741 = 0.010
ODE
0.810 + 0.002
0.729
0.005
0.746 + 0.006
ODE + L1
0.810 0.002
0.729
0.005
0.742 = 0.006
ODE L2
0.809 = 0.005
0.728
0.005
0.742 + 0.006
ODE LCI
0.810 + 0.002
0.003 0.731
0.743 - 0.006
ODE HGP
0.814 + 0.002
0.732 0.005
0.753 + 0.005
(c) IN (L)(L
Method
SUPPORT
METABRIC
RotGBSG
CoxPH
0,564 + 0.025
0,474 + 0.005
0.530 + 0.005
DeepHit
0.519 0.004
0.515 = 0.008
0.531 + 0.006
ODE
0.516 0.015
0.472 = 0.005
0.530 + 0.012
ODE + L1
0.515 + 0.015
0,470 = 0.005
0.533 + 0.013
ODE L2
0.517 0.017
0.470 0.004
0.537 1 0.010
ODE + LCI
0.516 + 0.015
0.469 + 0.002
0.528 + 0.010
ODE HGP
0.506 0.011
0.479 + 0.003 0.530 + 0.003
Table 2: Experimental Results on SUPPORT ablated in terms of sample size M. We set À = 10.
Method
mCC1(t)
mAUC(t)
iN
No reg.
0.771 +
0.810 0.002
0.516 0.015
M = 1
0.775 + 0.004
0.814 + 0.002
0.505 0.010
M
0.775 0.004
0.814 0.002
0.506 0.011
M 10
0.775 + 0.004
0.815 + 0.002
0.505 + 0.009
Tsows the results by varying the number of samples M in the sampling process ( S(L x) in equation 5
As long as the regularizer applied, the number of samples M does not affect the performance. Even when M 1. the regularizer works well. Figure 3 shows the results on SUPPORT and RotGBSG datasets by varying the cocfficient in equation]? 7 Since the performance variation by À is stable, the hyperparameter À can be tuned without much difficulty in practical setups.
Table 3 shows the time taken for training/vvaluation step for each methods. The time it takes to train ODE and to train ODE + L2 is slightly faster than vanilla ODE training time. It is straightforward to that the I.1 and 1.2 penalty makes smooth dynamics as the number of function

---

Under review as a conference paper at ICLR 2023
Variation (410-1
Veriaion (x10-1)
7.5
Veriatoon
6.5
5..
CDE
L1
12
LCI 1-1
1-1
000
101 A-1
o'
17
IC1
(a) SUPPORT
(b) METABRIC
(c) RotGBSG
Figure 2: A boxplot of the log-likelihood variation with respect the input perturbation ((a.,))De==1ll log pe x) - log pa x2)l2 on three datasets. We denote D==1 by the set of uncensored data. We choose Po 0) log po æ')lla is maximized under x' B(x,,) constraint. We set le across all experiments. The hazard gradient penalty effectively regularizes the variation. As the density of probability distribution p(- x) should be concentrated / for uncensored data (æ, the figures show that the hazard gradient penalty effectively regularizes the KL divergence between the density function at point x and that of neighborhood points x' B(x,e). e).
= + cg/llgl2 where g V log po x) so that
0'387
0.880
0.14
0776
053
0774
852
800
0768
0.66
1-
A=0
A-5
(a) SUPPORT mt)
(b) SUPPORT mAUC(()
(c) SUPPORT IN BBL(L)
0.445
07550
0335
07225
03522
8530
07200
03300
07176
8526
07400
8315
10
(d) RotGBSG mC (t)
(e) RoiGBSG mAUC(()
(1) RoIGBSG iNBLL())
Violin plots of experimental results on SUPPORT and RoIGBSG by varying A. Red, green, and blue denote mC'd, mAUC, and BLL. The thickness of a plot denotes the probability density of the results. The hazard gradient penalty may conflict with the negative log-likelin likelihood if we set high X. The that achieves the best scores across all metrics on the RotGBSG dataset could have been acquired between RotGBSG dataset in Table Ic for consistency.
Figure 3:
and À 10. However, report the result at À the
evaluations (NFE) in the training time of ODE LI and ODE + L2 are lower than that vanilla ODE. The decrease in NFE compensates for the overheads of calculating LI and L2 penalties, which makes the training time of ODE same applies to ODE HGP. Despite the additional sampling process and gradient calculation, the training time of ODE + HGP is on par with vanilla ODE. The decreased NFE thanks to smooth dynamics made by the hazard gradient penalty compensates for the additional computations. The smooth dynamics made by LI. L2, and the hazard gradiens penalty can also be observed in the evaluation phase. The NEE of ODE + LI, ODE + L2, and ODE HGP is much smaller than that of vanilla ODE. Not all regularizers make smooth dynamics. The NFE ODE + LCI is comparable to the NEE of vanilla ODE. Clearly, there's 's no point in regularizing L.CI in terms of inference speed.
and ODE 1.2 faster than that of vanilla ODE. The
8

---

Under review as a conference paper at ICLR 2023
Table 3: The taken time (in seconds) and the number of function evaluations (NFE) for each step in the training/evaluation time on RotGBSG. The numbers in parenthesis indicate relative performance against vanilla ODE. HGP incurs negligible overhead (1% slowdown) on training time while gives rise to 10% speedup on evaluation time. See Table for the taken time and the NFE for each step of regularizers on SUPPORT.
ODE
ODE + L1
ODE + L2
ODE + LCI
ODE + HGP
train time
0.1 (1)
0.1043 (0.90)
0.1125 (0.97)
0.1173 (1.01)
0.1179 (1.01) 0.0016 (0.80)
eval time
0.0020 (1)
0.0016 (0.80)
0.0014 (0.71)
0.0019 (0.97)
train NFE
13.207 (1)
11.519 (0.87)
10.637 (0.80)
13.207 (1.00)
12.129 (0.91) 9.929 (0.85)
eval NFE
11.893 (1)
9.946 (0.83)
8.411 (0.70)
11.643 (0.97)
RELATED WORKS
A line of research integrated deep neural networks to CoxPH (Faraggi & Simon 1995 Katzman
let al. 2018) and Extended Hazards (Zhong et al. 2021) for more model flexibility. Another line of research proposed distribution-free survival analysis models via the time domain discretization
[Le et al. 2018), adversarial learning approach (Chapfuwa et al. [018 or derivative-basse models Danks & Yau 2022). Previous works (Goldstein et al. 2020: Han et al.) 2021) proposed new objectives to optimize Brier score (Graf et al., 1999), Binomial log-likelihood, or distributional calibration directly. Yet to the best of our knowledge, none of the previous works focused on the
effect of gradient penalty on survival analysis models.
Previous works proposed 1.1 and 1.2 regularization in the survival analysis literature Tiisshirani 1997 Verweij & Van Houwelingen] 1994). Those methods regularize the survival analysis models so that the L.1 or 1.2 norm the model parameters does not increase so much. Our method is different from those methods in that penalize the norm of the gradient on each local data point.
Our method is closely related to semi-supervised learning (Chapelle et al. 2006). Among many semi-supervised learning methods, our method is germane to virtual adversarial training (Miyato) jet al., 2018) in that it regularizes function variation between a local data point and its neighbours. However, virtual adversarial training is different from ours in that the method was demonstrated in the classification setting and the output a discrete distribution.
In Generative Adversarial Nets (GANs) literature (Goodfellow al. 2014), the gradient penalty had been studied actively. Gulrajani et al. 2017) proposed the gradient penalty to satisfy the 1-Lipschitz function constraint in Santrovich-Rubinstenn duality. Mescheder al. (2018) proposed the gradient penalty to penalize the discriminator for deviating from the Nash equilibrium. Ours is different from these works in that we propose gradient penalty so that the density at x does not deviate much from that of x's 's neighborhood points.
6 CONCLUSION
In this paper, we introduced a novel regularizer for survival analysis. Unlike previous methods, we focus on individual local data point rather than global information. We theoretically showed that regularizing the norm of the gradient of hazard function with respect to the data point is related to minimizing the KL divergence between the data point and that of its neighbours. Empirically, we showed that the proposed regularizer outperforms other regularizers and it is not sensitive to hyperparameters. Furthermore, the proposed regularizer is computationally efficient and incurs an ignorable overhead. Nonetheless, as minimizing the proposed regularizer may conflict with opti- mizing the negative likeliiel,hood, practitioners should tune the balancing coefficient for each dataset. The paper highlights the new possibility that the recent advancements in semi-supervised learning could enhance the performance of survival analysis models.
9

---

Under review as a conference paper at ICLR 2023
REFERENCES
Lei Jimmy Ba, Jamie Ryan Kiros, and Geoffrey E. Hinton. Layer normalization. CoRR, Abs//607.06455) 2016. http:://a/a/ wwwpaiiww...ooaaaooooooo httg//bb/aaee 1607. 06450
James Bradbury. Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal
Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL http: //github. cowp///gg////og//aa/a/a/
Olivier Chapelle. Bernhard Schölkopf, and Alexander Zien (eds.). Semi-Supervised Learning, The MIT Press, 2006. ISBN 9780262033589.
Paidamoyo Chapfuwa, Chenyang Tao, Chunyuan Li, Courtney Page, Benjamin Goldstein, Lawrence Carin Duke, and Ricardo Henao. Adversarial time-to-event modeling. In Jennifer Dy and Andreas Krause (eds.), Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 735-744. PMLR, 10-15 Jul 2018. URL
mwr.iiinnsmmm...m...mmmmmmmmm..m press / / v80/ cwwp fwwallwwwll.. 8a. html
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary dif- ferential equations. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 31. Curran As- sociates, Inc., 2018. URL / hww.:iinnnsn.euuieeeeeeee neurips...... wt/ppaper//1///EEL....
pdf
David R Cox. Regression models and life-tables. Journal of the Royal Statistical Society: Series (Methodological), 34(2):187-202, 1972.
Dominic Danks and Christopher Yau. Derivative-bass neural modelling of cumulative distribution functions for survival analysis. In Gustau Camps- Valls, Francisco J. R. Ruiz, and Isabel Valera (eds.), Proceedings The 25th International Conference on Artificial Intelligence and Statistics, volume 151 of Proceedings of Machine Learning Research, pp. 7240-7256. PMLR, 28-30 Mar 2022. URL https::/proceedgggs mlr http://asaa........eeeee///////a
David Faraggi and Richard Simon. A neural network model for survival data. Statistics in medicine, 14(1):73-2, 1995.
Arnab Ghosh, Harkirat Behl, Emilien Dupont, Philip Torr, and Vinay Namboodiri. Steer Simple temporal regularization for neural ode. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Bal- can, and H. Lin (eds.), Advances in Neural Information Processing Systems, volume 33, pp- 14831-14843. Curran Associates, Inc., 2020. URL
neur Paper.
Mark Goldstein, Xintian Han, Aahlad Puli, Adler Perotte, and Rajesh Ranganath. X-cal: Ex- plicit calibration for survival analysis. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Bal- can, and H. Lin (eds.), Advances in Neural Information Processing Systems, volume 33, pp. 18296-18307, Curran Associates, Inc., 2020. URL https::/proceeddng/ neurips. Papeee
cc
lan Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sher- jil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K.Q. Weinberger (eds.), Advances in Neural Information Processing Systems, volume 27. Curran Associates. Inc., 2014. URL https:: https::/ neurips................ sc/papper/214//iie// 5a22999.222611333099977222e1fff3
Erika Graf, Claudia Schmoor, Willi Sauerbrei, and Martin Schumacher. Assessment and comparison of prognostic classification schemes for survival data. Statistics in medicine, 18( 18(17-18):2529 2545, 1999.
Stefan Groha, Sebastian M Schmon, and Alexander Gusev. Neural odes for multi-state survival analysis. https:/aarii..or/aaa/200..... 2020.
10

---

Under review as a conference paper at ICLR 2023
Ishaan Guirajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron C Courville. Im- proved training of wasserstein gans. In I. Guyon, U. Von Luxburg. S. Bengio, H. Wallach, R. Fer- gus, S. Vishwanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017. URL bt http:/ttttssoooeeiiiiig. ngs. neur ps ec/ pdf
Xintian Han, Mark Goldstein, Aahlad Puli, Thomas Wies, Adler Perotte, and Rajesh Ranganath. Inverse-weipe weighted survival games. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan (eds.), Advances in Neural Information P Processing Systems, volume 34, pp. newri cc/ 2160-2172. Curran Associates, Inc., 2021. /tR//ttttss::/preeeeeiii Paperr pdf
Hung Hung and Chin-Tsang Chiang. Estimation methods for time-dependent models with survival data. The Canadian Journal of Statistics / La Revue Canadienne de Statistique, 38(1): 8-26, 2010. ISSN 03195724, URL http:ttp:////
jstor www.orgg//teee wttg/staable/2285522.
How Jing and Alexander J. Smola. Neural survival recommender. In Proceedings of the Tenth ACM International Conference on Web Search and Data Mining, WSDM '17, pp. 515-524, New York, NY, USA, 2017. Association for Computing Machinery. ISBN 9781450346757. doi: 10.1145/
301 5018661.3018711 .3018719. URL 1BR.httpps///oo..rrr//00111111//
3018719
Edward L Kaplan and Paul Meier. Nonparametric estimation from incomplete observations. Journal of the American statistical association, 53(282):457-481, 1958.
Jared I. Katzman, Uri Shaham, Alexander Cloninger, Jonathan Bates, Tingting Jiang, and Yuval Kluger. Deepsurv: personalized treatment recommender system using a cox proportional hazards deep neural network. BMC medical research methodology, 18(1):24, 2018.
Patrick Kidger. On Neural Differential Equations. PhD thesis, University of Oxford, 2021.
David G. Kleinbaum and Mitchel Klein. Survival Analysis: A Self-Learning Text. Springer-Verles New York, 2012. ISBN 9781 178144199 101114196666 44196666, doi: 10. 1007/97 1-4419--64465
Håvard Kvamme, Ørnulf Borgan, and Ida Scheel. Time-to-event prediction with neural networks and cox regression. arXiv preprint 1rXiv:1907.00825, 2019.
Changhee Lee, William R Zame, Jinsung Yoon, and Mihaela van der Schaar. Deephit: A deep learning approach to survival analysis with competing risks. In Thirty-Second AAAI Conference on Artificial Intelligence, 2018
Jiayu Li, Hongyu Lu, Chenyang Wang. Weizhi Ma, Min Zhang, Xiangyu Zhao, Wei Qi, Yiqun Liu, and Shaoping Ma. A difficulty-aware framework for churn prediction and intervention in games. In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, KDD '21, pp. 943-952, New York, NY, USA, 2021, Association for Com- puting Machinery. ISBN 9781450383325.
doi: 10.1 10.1445/34475443447777. URL https:
//doi. 10. 1et..1114//44411 3447548 3467277
Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization. In International Confer-
ence on Learning Representations, 2019. URL
Bkg6RiCqY7
Lars Mescheder, Andreas Geiger, and Sebastian Nowozin. Which training methods for GANs do ac- tually converge? In Jennifer Dy and Andreas Krause (eds.), Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 3481-3490 PMLR, 10-15 Jul 2018. URL https://pproeeeinggg...... wtt.sm....peeee/////
mescheder 18a. ntml
Bruce D Meyer. Unemployment insurance and unemployment spells. Working Paper 2546, Na- tional Bureau of Economic Research, March 1988. URL/http: http:///wwwwwbbbbbo//ppperrr.r
w2546
Takeru Miyato, Shin-ichi Maeda, Masanori Koyama, and Shin Ishii. Virtual adversarial training: regularization method for supervised and supervise learning. IEEE transactions on pattern analysis and machine intelligence, 1999 2018.

---

Under review as a conference paper at ICLR 2023
Patrick D. T. O'Connor and Andre Kleyner. Introduction Reliability Engineering, chap- ter 1, pp. 1-18. John Wiley & Sons. Ltd, 2011. ISBN 9781119961200. doi: 10.1002/ 9781 19961260. URL 99811 9961
www.ew.wi....................... com/ http///ooaa/////aeee 1052
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang. Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, d d'Alché-Buc, E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems 32, pp. 6024-8035. Curran Associates, Inc., 2019. URL http:: httpt:::/ppperrs.eerppp.......
www.nneuii.... wtw.:/s..niiiss.ccc.ccooooo cctppprrr
pdf
Sebastian Pölsterl. scikit-survival: A library for time-to-event analysis built on top of scikit- learn. Journal of Machine Learning Research, 21(212):1 2020. URL httt:t/tt//mmm org/
papers/v21/20 729. 9html
Yulia Rubanova, Ricky Q Chen, and David K Duvenaud. Latent ordinary differential equations for iregularly-sampled time series. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché- Buc, E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems, vol- ume 32. Curran Associates, Inc., 2019. URL https://proceediing neurips.comm 7be
pdf.
Patrick Schwab, Arash Mehrjou, Sonali Parbhoo, Leo Anthony Celi, Jürgen Hetzel, Markus Hofer, Bernhard Schölkopf, and Stefan Bauer. Real-time prediction of covid-19 related mortality using electronic health records. Nature commumications, 12(1):1-16, 2021.
Harald Steck, Balaji Krishnapuram, Cary Dehing-oberije Philippe Lambin, and Vikas C Raykar. On ranking in survival analysis: Bounds on the concordance index. In J. Platt, D. Koller, Y. Singer, and S. Roweis Advances in Neural Information Processing Systems, volume 20, Curran As- sociates, Inc., 2007. httppp/pps//ppppooeeeee......nrr
Weijing Tang, Kevin He, Gongjun Xu, and Ji Zhu. Survival analysis via ordinary differential equa- tions. Journal of the American Statistical Association, (just-accepted):1-41 2022a.
Weijing Tang, Jiaqi Ma, Qiaozhu Mei, and Ji Zhu. Soden: A scalable continuous time survival model through ordinary differential equation networks. Journal of Machine Learning Research, 23(34):1-29, 2022b.
hmtp://mmmlrrr http:////mmmlllrrr//pppee/////// 900 htm
Robert Tibshirani. The lasso method for variable selection in the cox model. Statistics in medicine, 16(4))385-39,, 1997.
Hajime Uno, Tianxi Cai, Michael J Pencina, Ralph B D' Agostino, and Lee-Jen Wei. On the c- statistics for evaluating overall adequacy of risk prediction procedures with censored survival data. Statistics in medicine, 30(10):1005-111 2011.
Pierre JM Verweij and Hans C Van Houwelingen. Penalized likelihood in cox regression. Statistics in medicine, 13/23-24):247--433, 1994.
Lee-Jen Wei. The accelerated failure time model: a useful alternative to the cox regression model in survival analysis. Statisties in medicine, 11(14-15)117111179, 1992.
Qixian Zhong, Jonas W Mueller, and Jane-Ling Wang- Deep extended hazard models for survival analysis. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, volume 34, pp. 15111-15124 Curran As- www..nneripp...e... wttpu//sscccppppee///2/////////.
sociates, Inc., 2021. URL http::tttt://atto //pp:
12

---

Under review as a conference paper at ICLR 2023
Table Summary statisties of the datasets used in our experiments. N denotes the number of data points and denotes the dimension of each data points.
Durations
Event Quantiles
Dataset
N
Censoring (%)
# unique
domain
f
25%
50%
75%
SUPPORT
9105
43
31.89%
1724
N+
14
58
252
METABRIC
1904
42.06%
1686
R²
42.68
85.86
145.33
RotGBSG
2232
43.23%
1230
R-
13.61
24.01
40.32
Algorithm Hazard Gradient Penalty
Require: hg. learning rate Y repeat
sample (ar, e) D
Retrieve unique times f1,..., fx from minibatch.
integrate hh (t æ) from f and store log Stiti æ)...., log So(tk x)
log (t x t choose from log (t1 x).. log So 2) that corresponds to (æ, t, e)
log (t xx) log ht x + log So x))
L. <- - Po =) (1 e) log St xx)
Negative log eigliiklihood
sample i1, ..,...; iM ~ Categorical(S x).... Sol1x x))
So x
I'm
4m -1 Uniform(44.,
D 4
RS
m= V/ (f/ x2)12
M
Hazard gradient penalty
+ XRx)
until Convergence
A EVALUATION METRICS AND EXPERIMENTAL DETAILS
A.1 EVALUATION METRICS
Throughout this subsection, we denote S(1 w) as the estimate of S(t æ). I(-) the indicator function, (z:, T;, e;) the ith covariate, time, event indicator of the dataset, G(t) the Kaplan- Meier estimator for censoring distribution (Kaplan Meier 1958 and W; as 1/G(T.).
A.1.1 TIME DEPENDENT CONCORDANCE INDEX (Ctd)
The concordance index, or C-index is defined as the proportion of correctly ordered pairs among all comparable pairs. We use time dependent variant of C-index that truncates pairs within the prespecified time point Uno et al. (2011). The time dependent concordance index at Ord(t), is defined as
MM & - < <T, < ()(((()) w.) < Š(L ax))
& - <T T <)
To evaluate Cdd at [t,....,. [ the same time, take its mean mC DEL (td(tt).
A.1.2 TIME DEPENDENT AREA UNDER CURVE (AUC)
an extension of the ROC-AUC to survival data Hung & Chiang (010). It measures how well
model can distinguish individuals who fail before the given time (7, < and who fail after the given time (T) > t).
The AUC at time f, AU)( is defined
MN Σ I(T, > ()T(T, S t(,,I(((tt) xx) S S(t xx)))
EN ITT + I(T. S t)w:)
To evaluate AUC [I,,..., . lr.] at the same time, take its mACC Σt AUC(L).
13

---

Under review as a conference paper at ICLR 2023
A.1.3 NEGATIVE BINOMIAL LOG-LIKELIHOGE
We can evaluate the negative binomial log-likelihood (NBLL) measure both discrimination and calibration performance K Kvamme et al. (2019). The negative binomial log-likelihood at é measures how close the survival probability to the given data survived and how close the survival probability is to O if the given data failed before t. The NBLL at t, NBLL((), is defined as
log(1 S(t aT) st, = =) + log S(t 2 xT) > t)
((T))
G(1)
For the convenience of evaluation, we integrate the NBLL, IN BLL = b2- S NBL((t)dt.
A.2 EXPERIMENTAL DETAILS
Across all datasets, we split training set, validation set and test set into 70%, 10% and 20% each using PyTorch's 's random_splis function Paszke et al. (2019). We set seed ting.
when split-
Across all experiments, we use an MLP with two hidden layers where each layer has 64 hidden units. Across all layers, we apply Layer normalization Ba et al. (2016). Instead of naively feeding time t into the neural network, we feed scaled time = (t (2) f/) f1) where f1, f2. and t3 first, second, and third quartile of failure event times. We found that this strategy enhances the ODE model's performance and boosts training time. To incorporate time into the survival analysis model, project the time into eight dimensional vector using single layer MLP and then concatenate to the input data. The time is also specified by adding projected output into each layer output. We use the AdamW optimizer norm so that it does not exceed 1. We set the learning rate to 0.001. We have implemented the code using JAX Bradbury et al. (2018) and Diffrax Kidger 2021
Loshchilov ootimizer & Hutter (2019) and clipped the gradient
To find the best À in equation 7 we run experiments with À = 1,5, 10, 50 and report the results A 10 as it shows decent performance across all metrics and datasets. We also have set the number of samples M from the time sampling process & ~ Po(t x) in equation 5 We set M across all the hazard gradient penalty experiments. To find the best coefficient a ODE LI, ODE + L2, and ODE + LCI experiments, we set a le - 1,1e 2, le and run the experiments. We report the best a terms mAUC. To report mC and mUC,, we calcu- late Ctd and AUC at 10%, 20%, .... 90% event quantiles and average them. To report BLLL we integrate from the minimum time the test set to the maximum time of the test set. We use Pölsterl (020 to report mCdd and mAUC. We use pycox Kvamme et al. (2019) to report iN B..L. Across all experiments, we run experiments with different seeds and report their mean and the standard deviation.
B CODE SNIPPET OF HGP
import jax
import jax. random jrandom
import jax.nn nn jnn
import jax. as jnp
odeint ode_nt (hazard func, args.rtol, args.atol)
@jax.jit
def calc_dh_dx(
10
params: optax Params, t: Inarray,, X; ndarray
12
def sum hun_h_func (params, t, X):
13
zercs jnp. zeros shapeape 1))
14
states_t0 jnp. concatenate (X, zeros), axis=-1) return jnp. sum (hazard_ func (t, states_t0. to, args=parans) [i, -1])
15
The code will he made publicly available in the near future.
14

---

Under review as a conference paper at ICLR 2023
17
return (rum argnums=2) (params, t, x)
18
19
@jax. jit
def loss_func(
params: optax. Param,, X: jnp.ndarray, event: inp.nayrray,
22
t: jnp. timestamp: jnp. ndarray, unique_ids idx: jmp.ndarray. key: jap. ndarray
23
24
batch size X. shape [0] time_size jnp. size (timestamp)
27
28
timestamp jnp. ins.insert (timestamp, 0, 0)
int hazara t0 inp. zeros ( (X. shape(01, 11)
states t0 inp. concatenate (XX int hazard t0), axis=-1) ys odeint (timestamp, states_tO, params)
31
32
int_hazarz ysl1: ;, -1]
33
log_surv int_hazard
log_surv_t log_surv [unique_idx, inp. arange (x.sh [0])]
36
hazard_t hazard_func states_to, args-params) [:, -1]
37
38
probtt jnp. log (hazard_t 1e-6) sur_surv_t
40
assert log surv t. shape log prob shape assert event. shape -- log_prcb_t, shape
41
42
nl1 (event log_prob_t + (1. event) log_surv_t) mean ()
45
now [batch_size, time_ time_size]
46
logits transooee (jax. waw..la...tto. stx.stoo_gradient (log_surv), (1, 0))
47
= jrandom, categorical
key, jnp. tile (logitsl;, ;, Nonel, (1, 1, args. sample_sipeeeeeeeeeeeeeeee_eemee axis=1
19
50
51
to timestamp [0:] [time_idx]
tl timestamp [1:] [time_ióx]
54
t_prime + (tl-t0) jrandem. uniform(key shape=t0. shape)
55
assert shape == (batch size, args. size)
(t x) II \nabla_x h(t x) dh_ox jax. vmap (calc_dh_áx, in_axes= (None, 1, None) (params, t_prime X)
58
59
assert dh dx shape == (args.sample_size batch_sizae size, shape dh cx__normmm jnp. linalg norm ord=2, axig=-1) mean
62
loss all args. lamoda dh dx norm
63
return Loss, (nil, dh norm)
C THE EXISTENCE OF THE PROBABILITY DISTRIBUTION so(t x))
In general, we cannot guarantee the existence of probability distribution
So(t x))
so(1 x)) J Sa x)
as the integration f (t x)dt may not exist.
To ensure the existence of Se(t dt, we simply add a constraint: he(t x)) 2 c. Here, é is very small constant (c.g. le 8). As is very small constant, has negligible impact on the algorithm. The constraint can be achieved easily by adding to the softplus output of the hazard function.
15

---

Under review as a conference paper at ICLR 2023
Table 5: Comparison against STEER.
(1)
Method
SUPPORT
METABRIC
RotGBSG
ODE STEER
0.771 + 0.003
0.695 + 0.008
0.718 + 0.005
0.772 + 0.002
0.699 = 0.010
0.718 + 0.006
HGP
0.775 + 0.004
0.702 + 0.009
0.723 + 0.006
(b) mAUC(T)
Method
SUPPORT
METABRIC
RotGBSG
ODE
0.810 + 0.002
0.729 + 0.005
0.746 + 0.006
STEER
0,810 = 0.001
0.730 = 0.006
0.745 + 0.006
HGP
0.814 + 0.002
0.732 = 0.005
0.753 = 0.005
(c) N. BTLT(L)
Method
SUPPORT
METABRIC
RotGBSG
ODE
0.516 + 0.015
0.472 + 0.005
0.530 + 0.012
STEER
0.521 + 0.012
0.475 + 0.011
0.523 + 0.006
HGP
0.506 + 0.011
0.479 + 0.003
0.530 + 0.003
Table 6: The taken time (in seconds) and the number of function evaluations (NFE) for each step in the training/evaluation time on SUPPORT. The numbers in parenthesis indicate relative perfor- mance against vanilla ODE. HGP boosts training time (8% faster) and gives rise 9% speedup evaluation time.
ODE
ODE + LI
ODE + 1.2
ODE + LCI
ODE HGP
train time
0.3392 (1)
0.2688 (0.79)
0.2758 (0.81)
0.3434 (1.01)
0.3135 (0.92) 0.0018 (0.91)
eval time
0.0019 (1)
0.0015 (0.77)
0.0015 (0.78)
0.0018 (0.95)
train NFE
10.980 (1)
9.214 (0.83)
9.151 (0.83)
10.973 (0.99)
10.556 (0.96)
eval NFE
8.580 (1)
7.073 (0.81)
7.073 (0.82)
8.446 (0.98)
7.988 (0.93)
Under the constraint, Sp(t x>) is bounded by exponential function:
Sp x) exp (-/ h(T z)dr)
exp(-at)
Also, f St æ)dt bounded:
f So x) - <f fexp(-et)dt
D COMPARISON AGAINST NEURAL ODES SPECIFIC REGULARIZER
In this section, we compare HGP against Neural ODEs specific regularizer: STEER (Ghosh et al. 2020). STEER regularizes Neural ODEs by perturbing the final time of the integration. Table 5 compares the performance of HGP against STEER. Overall, HGP outperforms STEER on various setups.
E THE TAKEN TIME AND THE NUMBER FUNCTION EVALUATIONS
In this section, compare HGP against competing regularizers in terms of the taken time and the number of function evaluations on SUPPORT dataset. See Tor the details. Overall, the result aligns with that of Table 3
16


Reviews:
[{"review_id": "wHAnYCZFS2", "reviewer": "Reviewer_CgBz", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper proposes a hazard gradient penalty algorithm for regularising neural ode-based survival models. Additionally, the paper formulates a theoretical connection of the proposed hazard gradient penalty to KL divergence  minimization of  KL ($p(t|x) || p(t| x^{\\prime})$), where $x^{\\prime}$ are neighbouring data points of $x$. Experimental results on three real-world datasets show performance improvement over baselines per metrics C-Index, negative binomial log-likelihood (NBLL), and AUC. ", "summary_of_the_review": "Some researchers may find the theoretical analysis of the proposed gradient penalty interesting. However, the paper is lacking due to  several technical issues and underwhelming experiments.", "strength_and_weaknesses": "**Strengths**:\n- The paper focuses on survival analysis, an important problem across many domains\n- Theoretical analysis of connections of the hazard gradient penalty to KL divergence  minimization of  KL ($p(t|x) || p(t| x^{\\prime})$), where $x^{\\prime}$, may be of interest to some researchers\n- Paper provides a snippet code implementation in JAX\n\n**Weaknesses**:\n\n*Overall Clarity*\n-  The title of the paper needs to be more focused since the proposed regularizer  is mostly applicable to neural ODE survival models\n- The paper should include a figure, illustrating (i) model parameters $\\theta$; (ii) how $h_{\\theta}(t|x)$ are learned, and (iii) model outputs at test time.\n\n*Technical Quality*\n- The paper claims modeling $p(t|x)$ or $S(t|x)$ is challenging:  This is not necessarily true, since there exist several approaches modeling $p(t|x)$ or $S(t|x)$ including AFT, Lee et al. (2018),  Chapfuwa et al. (2018), *etc.*\n- It's unclear how the paper imposes the survival  constraints for learning valid distributions: Survival function is a monotonically decreasing function and $\\lim_{t \\rightarrow \\infty } S(t|x) =0$\n- It's unclear what the model outputs at test time.  E.g., Is the model predicting hazard, event times, or survival probabilities?\n- Eqn 2: Should be $P(t < T < t + Δt | T> t, x)$\n- Eqn 3: What are the inputs to the ODE solver?\n- Paper claims modeling the hazard function is better that Tang et al. (2022b)’s cumulative hazard function without providing any justification, *e.g.*, experimental results\n- Eqn 5: Should be $t \\sim p_{\\theta}(t|x)$, it's unclear how the paper samples from the survival function $t \\sim S_{\\theta}(t|x)$?\n- Hazard gradient penalty: It's unclear why the ODESolver is evaluated at varying batch-specific times, instead of shared range $[0, t_{\\rm max}]$\n- Eqn 7: How is $\\lambda$ selected?\n\n*Underwhelming Experiments*\n- The performance improvements over the vanilla neural ODE baseline seem marginal:\n*  How significant are the claimed improvements\n* Table 3: Computational complexity analysis is only provided for one dataset. How does the proposed approach scale with larger datasets (*e.g.*, Support), batch sizes, *etc.*? \n-  Paper should also include qualitative results comparing predicted $h(t|x), S(t|x), p(t|x)$ against ground truth times and baselines\n- For completeness, the paper should consider comparisons to other non-ODE baselines, e.g., AFT, Lee et al. (2018),  Chapfuwa et al. (2018), etc.\n- Table 2: It seems the model performance is not sensitive to increasing from $M=1$ to $M=10$. Is this expected?\n\n*Novelty*\n- Gradient penalty regularisation has been proposed before. An expanded discussion highlighting how the proposed approach is better (or different) than previous approaches, *e.g.*, 1-Lipschitz function constraint in W-GAN is crucial.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "w3bvYAeyibp", "reviewer": "Reviewer_T7HL", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "ODEs have recently found their applications in survival analysis where in the case of right-censoring the log-likelihood can be optimized quite straightforwardly using an ODE solver. This work is inspired by an assumption known the clustering literature according to which the decision boundaries do not cross high-density regions in the data space. They use the same intuition in their hazard estimation. They propose a regularizer on the hazard gradient which is the expectation of the hazard gradient with respect to a random variable whose PDF is the normalized survival function. They show that it is equivalent to minimizing the KL divergence between the PDFs of two input vectors which are \"close\" in the input space.", "summary_of_the_review": "I find this paper interesting but at the same time I am not fully convinced of its applicability and superiority over baselines. Not sure about applicability because it only incorporates time-static covariates and not sure about superiority given the missing comparisons mentioned above.", "strength_and_weaknesses": "* A major short-coming of this work is being restricted to time-static covariates. There has been recent work on time-varying covariates in survival analysis so considering time-static covariates only comes at the price of usability.\n\n\n* Since you are proposing a hazard estimator, it would have been beneficial to also provide RMSE results on recovering the true hazard in a dataset. Of course for that you needed to create a synthetic dataset where you simulate survival data using a known hazard function but it would have been beneficial.\n\n\n* The only non-ODE-solver baseline was CoxPH, whose dependence on input covariates is linear. Comparing against non-linear hazard estimators or even CoxPH extensions whose covariate effect is non-linear would have made the a stronger case of the developed algorithm. Without the mentioned comparisons, only ODE solvers are compared with one another.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "W9BclWFvDtB", "reviewer": "Reviewer_pnuq", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors propose a hazard gradient penalty on ODE models to ensure the smoothness on the density in survival analysis. The main contribution is this hazard gradient penalty.", "summary_of_the_review": "It's unclear whether we should use the proposed method. The motivation from ssl is vague and there is obviously another way to penalize the function based on the theoretical claim. So I give \"reject\".", "strength_and_weaknesses": "Strength:\n1. The authors propose a novel hazard gradient penalty specifically on survival analysis problems using ODE methods. \n\nWeakness:\n1. The use of the hazard gradient penalty is not well-supported. \na. The authors claim that \n\"The cluster assumption from semi-supervised learning states that the decision boundaries should not cross high-density regions (Chapelle et al., 2006). In a similar vein, hazard functions of survival analysis models should change slowly in high-density regions.\" I am not sure how the idea of cluster assumption applies here in survival analysis. It does not sound similar to me.\nb. The authors somehow prove that their gradient penalty is an upper bound on KL-divergence between p(t|x) and p(t|x'). There are many ways to penalize that. I am not sure why the authors care about the gradient of the hazard. You can actually directly penalize the gradient of p(t|x). You can also model the CDF as an ODE[1] and do not need to model the cumulative hazard. You probably also want to cite [1] since they also talk about using ODE in survival analysis and is very similar to SODEN.\n2. Empirical results are not very convincing. The authors do not tune much for the baseline ODE model. But it's actually sensitive to many hyperparams. See the SODEN repo. They random search 100 groups of hyperparams, select the best on the validation and report on the test. If the authors can show that their method is better than the numbers reported in SODEN, it will be more convincing empirically. Some large real world experiments may also be needed (for example, MIMIC, SEER). The current datasets are not very large and do not contain many features which may not be the right test-bed for deep learning.\n\n[1]Derivative-based neural modelling of cumulative distribution functions for survival analysis.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}, {"review_id": "DiaHegXV6u", "reviewer": "Reviewer_ReJF", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper applied the neural ODE method to log trasformed survival rate model expressed by integration of hazard funciton. The penalty proposed is hazrd gradient penalty. The equivalence to KL divergence is proved, and empiral comparison to four other penalties under three criterions demonstrated the benifit of the proposed hazard gradient penalty.   \n", "summary_of_the_review": "This paper proposed hazard gradient penalty to address survial analyses using ODE framework facusing on data point rather than global information. The general presentation is good. The experimental result (if the author can clarify/defend the result with more than only ten quantiles selected) compared to existing method based on the paper showed advantage of implementing this regularizer. \n", "strength_and_weaknesses": "This paper started from a very general form of survival analyses, making the proposed method applicable to broad survival models shown in section 2.1. The alogorithm proposed is computationally efficient than KL convergence by requiring less derivation of intermediant values. Theoretical proof of the equivalence of KL divergence and the hazard gradient method is proved. The experiments performed on three public datasets showed the proposed hazard gradient penalty method perform good under the three adopted criterions. The reviewer like the idea in A2 of feeing scale time to boost training time. \n\nIn A2, it clarifies \"calculate C, and AUC at 10%, 20%, ..., 90% event quantiles and average them\". This reviewer is wondering the comparision of the proposed method and existing one would be like what if closer quantiles (say each 1% apart), or if higher weight to cetral quantiles. Or say, would it possible the proposed method with estimated hazard change slow in boundary to make calculated mAUC etc look good while not better performance than existing penalties at the central high-density regions.      \n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject