Survival Analysis via Density Estimation
Authors: Hiroki Yanagisawa
Shunta Akiyama
Venue: ICML 2025 poster (2025)
Keywords: survival analysis
time-to-event analysis
competing risks
density estimation
Primary area: general_machine_learning->supervised_learning

Abstract: This paper introduces a novel framework for survival analysis by reinterpreting it as a form of density estimation. Our algorithm post-processes density estimation outputs to derive survival functions, enabling the application of any density estimation model to effectively estimate survival functions. This approach broadens the toolkit for survival analysis and enhances the flexibility and applicability of existing techniques. Our framework is versatile enough to handle various survival analysis scenarios, including competing risk models for multiple event types. It can also address dependent censoring when prior knowledge of the dependency between event time and censoring time is available in the form of a copula. In the absence of such information, our framework can estimate the upper and lower bounds of survival functions, accounting for the associated uncertainty.

Survival Analysis via Density Estimation
Hiroki Yanagisawa Shunta Akiyama'
Abstract
target variable. Given this conceptual similarity, numerous methodologies originally developed for density estimation have been adapted for survival analysis, particularly in sce- narios with K = 2 under the conditional independence assumption. For instance, random forests (Breiman, 2001) for density estimation have been adapted into random sur- vival forests (Ishwaran et al., 2008) for survival analysis, and modern neural network models for density estimation have been extended into models such as DeepHit (Lee et al., 2018). Furthermore, strictly proper scoring rules for density estimation (Gneiting & Raftery, 2007) have been adapted for survival analysis (Rindt et al., 2022; Yanagisawa, 2023). Calibration metrics, such as the expected calibration er- ror (Naeini et al., 2015; Guo et al., 2017), have similarly been extended to D-calibration for survival analysis (Haider et al., 2020).
This paper introduces a framework for survival analysis by reinterpreting it as a form of density estimation. Our algorithm post-processes density estimation outputs to derive survival functions, enabling the application of any density estimation model to effectively estimate survival functions. This approach broadens the toolkit for survival analysis and enhances the flexibility and applica- bility of existing techniques for density estimation. Our framework versatile enough to handle var- ious survival analysis scenarios, including com- peting risk models for multiple event types. It can also address dependent censoring when prior knowledge of the dependency between event time and censoring time is available in the form of a copula. In the absence of such information, our framework can estimate the upper and lower bounds of survival functions, accounting for the associated uncertainty.
Despite the numerous extensions of density estimation methodologies for survival analysis, these adaptations en- counter several limitations. First, these adaptations are typi- cally tailored to specific methodologies on a case-by case-by-case by-c basis, necessitating the development of a new customized extension for each novel density estimation method. Sec- ond, most survival analysis models rely the conditional independence assumption (or even stronger assumptions such as the proportional hazards assumption (Cox, 1972)), which may not hold in various real-world applications. This issue highlights the need for survival models that operate under weaker assumptions. Third, while Tsiatis (1975) demonstrates that the survival function cannot be identified without making any assumption, the estimation of the up- per and lower bounds of the survival function has not been fully investigated under the condition where no assump- tions can be made. Fourth, while the existence of a strictly proper scoring rule for survival analysis with K = under the conditional independence assumption has been estab- lished (Rindt et al., 2022) and can be used as an evaluation metric, strictly proper scoring rule has been established for K 2.
1. Introduction
Survival analysis constitutes statistical methodology tensively employed across diverse domains, including medicine, engineering, and social sciences, to analyze and predict the probability distribution of the time until the oc- currence of an event of interest. In numerous real-world scenarios, the phenomenon of competing risks emerges when an individual is subject multiple mutually exclu- sive events, where the occurrence of one event precludes the observation alters the probability of the other events. The field of survival analysis with competing risks has a rich historical foundation, as extensively reviewed in survey papers (Wang et al., 2019; Wiegrebe et al., 2024).
Survival analysis with K competing risks can be conceptual- ized as variant of density estimation, a subfield of machine learning aimed at estimating the probability distribution of
In this paper. we propose novel framework that reinterprets survival analysis through the lens of density estimation. By post- post-processing the outputs of density estimation models in the form of cumulative incidence functions, our algo- rithm derives survival functions, thus enabling any density estimation model to be effectively utilized for survival anal-
'CyberAgent, Tokyo, Japan. Correspondenee to: Hiroki Yanagi- sawa <yyanggsawaa_hrokii @ sybbraagent.co.pp...
Proceedings of the nnnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 20255 by the author(s).

---

Survival Analysis via Density Estimation
Cumulative incidence function V2(([z)
Output (1(()))
k =
==
> b
k= = 2 :
k ==2
Step 1
Step 2 postprocess
Input x
t
estimate
:
Figure 1. Two-step (TS) algorithm for survival analysis with K competing risks: it first estimates the cumulative incidence functions (t]z)) via density estimation and then postprocesses them to obtain the outputs Fu (t)z).
ysis, as illustrated in Fig. 1. This reinterpretation not only broadens the toolkit available for survival analysis but also enhances the flexibility and applicability of existing density estimation techniques. Our framework's versatility is ev- ident in its capacity to handle diverse range of survival analysis scenarios, including competing risks and dependent censoring.
age survival function, the Kaplan-Meier estimator (1958) is widely used when the conditional independence assumption is valid. This estimator has been extended as the copula- graphic (CG) estimator (Zheng & Klein, 1995; Carrière, 1995), which operates under the same assumption as ours, that copula is provided as prior knowledge. Therefore, two-step algorithm can be seen as extension of the CG estimators to estimate the individual survival function.
The contributions of this paper are summarized as follows:
Upper and lower bounds estimation. Even in scenarios where prior knowledge of the dependency is absent, our framework capable of estimating the upper and lower bounds of individual survival functions, accounting for the uncertainty arising from the lack of knowledge about the copula. Regarding the estimation of the average survival function, Peterson (1976) presents the upper and lower bounds estimation. Our method can be seen as an extension of this approach for estimating individual survival functions.
Model-agnostic framework. Our two-step algorithm can be integrated with any density estimation model, allowing for the selection of density model based on various cri- teria such as prediction performance, training time, and interpretability. To demonstrate the effectiveness of our al- gorithm, we show that our algorithm, when combined with the LightGBM model (Ke et al., 2017), outperforms base- line models on real datasets. Furthermore, we prove that if the density estimation can achieve an arbitrarily small error c (as the number of data points increases), our algorithm can estimate survival functions with small error for K 2 under plausible assumptions (see Appendices C and D).
Strictly proper scoring rule for competing risks. Given a probabilistic output, a strictly proper scoring rule usually evaluation metric. Several proper ing rules exist for survival analysis with K 2 (Rindt et al., 2022; Yanagisawa, 2023), but no scoring rule has been proven to be proper for K > 2. In this paper, we intro- duce a new strictly proper scoring rule, termed NLL-SC, for K 2 2. NLL-SC based on the copula and can be used as an evaluation metric for survival analysis under dependent censoring, addressing the difficulties discussed by Gharari et al. (2023) in defining an appropriate evaluation metric under dependent censoring. Additionally, by utilizing NLL- SC, we construct a new monotone neural network model for K competing risks based on the copula-based assump- tion, named the survival copula network (SC-Net), which can be seen as an extension of the monotone neural net- work model for K under the conditional independence assumption (Rindt al., 2022).
employed
Dependent censoring. Our approach addresses the chal- lenge of dependent censoring, a prevalent issue in survival data where the censoring mechanism is not independent of the event time. To handle dependent censoring, a widely used assumption is that prior knowledge regarding the de- pendency between event time and censoring time is available in the form of copula (Emura & Chen, 2018). Our algo- rithm is based on this copula-based assumption and, there- fore, relies on weaker assumptions than the conditional inde- pendence assumption. Additionally, our two-step algorithm is capable of handling situations where the dependency in- formation is provided as a parameterized copula, along with appropriate assumptions on the form of the survival function to guarantee the identifiability of the parameter.
While our framework is primarily designed to estimate the individual survival function, for the estimation of the aver-
2

---

Survival Analysis via Density Estimation
2. Preliminaries
R/13
Rix
Rai
Rax
Ra.1
In this study, we investigate survival analysis with K com- peting risks. Let X denote the random variable representing feature vector whose support is X. Let Ti, .,, TKk be the K random variables corresponding to the event times of K distinct types of events, with each T3 supported on Rzo. Due to censoring, direct observation of samples (f1, .() ~ (2.,72..... Tx) is not feasible. How- ever, we can observe their minimum T = m | and the corresponding index A = arg min; (Tk) We note that in much of the existing literature on survival analy- sis, the index A starts at 0 (i.e.,& A E {0. 1) for K = 2); how- ever, in this paper, we assume A starts at (i.e., A E {1,2) for K = 2).
R2a
6-
Rao
R1o
(a) Observations in 2D space
(b) CJD divided by Rsa
Figure 2. (a) Two observations ((x(1, 20, 1), (2x()),35, ,35, 2)} are illustrated as vertical and horizontal line segments in the twa- dimensional space. (b) The CJD space with B 6 and K 2, which is divided into subregions Rs.k-
The primary objective of survival analysis is to estimate the cumulative distribution function (CDF) (t/x) = Pr(T): S t/æ>) of Tk conditioned on x X for each k [K] from samples (z) iid. (X, A), where [K] = (1,2...., ,,). In this study,
In this study, we refer to this partitioned region as the Cen- sored Joint Distribution (CJD) representation.
frame the prob- lem to estimate (e(()))) for each / in a finite set of times {o such that 0 = So < C1 < .... < {B. where {B is sufficiently large to ensure 0 S < SB for any observed time l ~ T. We assume that the true (4(())) a continu- ous and monotonically increasing function of t, satisfying Fi (Colz))) 0 and (<)(( This discretization approach is commonly adopted in nu- merous survival models (e.g.. (Leee et al., 2018; Yanagisawa, 2023; Hickey et al., 2024)).
Copula and survival copula. Many models for survival analysis with K operate under the conditional indepen- dence assumption, denoted as Ti IL T2 or even stronger assumptions such as the proportional hazards assumption, exemplified by the Cox model (1972). In this work, we adopt more flexible assumption that the dependence struc- , Tx can be modeled using a copula,
for all & E [K] and x E X.
ture among Ti, .......777 a widely recognized method for capturing dependencies in survival analysis (Emura & Chen, 2018; Gharari et al., 2023; Zhang et al., 2024).
While survival analysis often focuses on estimating the sur- vival function, defined as 5 (f)x: = F& (f( this study aims to estimate the CDF (t of Tk- Additionally, we consider the estimation of the average CDF (t) of (tlx)) over X and the average survival function (t = 1 - (t).
In probability theory and statistics, copula is defined as multivariate cumulative distribution function where each univariate marginal distribution is uniformly distributed over the interval [o, 1]. Copulas are particularly useful for charac- terizing the dependence structure among random variables. Formally, copula is defined as follows (Nelsen, 2006; Gharari et al., 2023).
Censored Joint Distribution (CJD) representation. Given an observation (,) (X,T, we interpret the pair (t, 6) in a K -idimmnnionnal space. For instance, in the case where K = 2, the pair (t, 8) can be represented as a line segment in a two-dimensional plane, as illustrated in Fig. 2(a). In this figure, the pair 51) = (20, 1) observed for (1() is depicted as a vertical line segment, indicating that Event is observed at time t1- = 20 and to 2 f1. Similarly, the pair (t, 6) (35, 2) observed for (2) is represented as horizontal line segment, indicating that Event 2 is observed at time f2- = 35 and t 2 f2.
Definition 2.1. (Copula.) A copula is K-dimensional function C: [0, 1|" + [0,1] that satisfies the following conditions:
(i) C(u1, ((u1,u2,... Ug,..., Uk- O, 02+11-... uk) = 0 for ev- ery u E [0,1] 1]
(i) C(I,..... 1, u. L,..,...,11 1) for every [0, 1],
(ii) Given u, E [0, 1] such that us < U2 is valid for all k E [K], the following condition is satisfied:
Given that the time horizon is discretized by the boundaries Hz the K-dimensional space is partitioned as illus- trated in Fig. 2(b) by defining the set Ro, of observations (t, 8) ~ (T, A) as follows:
Σ (-1)1t1at-m+t [E(0,1)
-C(uft 11 1-h 12 12121112
u uk 1--k] 2 0,
R = (tt, : C6 < C441, 5 k).
where = (41 ... ......, ...,, ) a length- K binary vector.
3

---

Survival Analysis via Density Estimation
notable example of copula is the independence copula, defined
Subsequently, utilizing the estimated (C(=),, we estimate the probability
K
Cina(ur, u2.-... uk)
II
uk-
(1)
bb Pr((t, 8) E Ro,s(2)
1=1
within the context of the CJD representation via the follow- ing equation:
Another example is the Frank copula for the bivariate case with non-zero parameter 0:
fb,A/= = (<6ll)) -aa--VVisssa Vi
(5)
CFrant(Uu, w2) = log (1+ (e (e-0ur - () Ou2
1)
c 6-
Thereafter, the output distribution (C((C/o) is derived from the estimation fo,k|z assuming that we have prior knowl- edge regarding the dependencies among the random vari- ables T1, T2..... TK in the form of copula C. It is noteworthy that the conditional independence assumption (Ti 1 T2|X) used in many survival models for K = 2 is equivalent to using the independence copula (1) for the copula C.
A significant attribute of copulas is that the joint distribu- tion Pr(Ti < 4, T2 S 1,,..., Tk S (x) can be uniquely represented using a copula, as stated in Sklar's theorem.
Theorem 2.2. (Sklar's Theorem (1959).) There exists a copula C such that for all t,, ,,,
Pr(Ti S t1, T2 S .,, .......,T Tk S tt)
(b)) Fa(12)....., Fx((x))
3.1. Step 1: Cumulative Incidence Function Estimation
If the marginal distribution Fg is continuous for all k; then C is unique.
A straightforward methodology for estimating the CIF (4) is the application of a distribution regression model, which is specifically engineered to estimate conditional CDF Flylz) Pr(Y S y|æx) given samples (zz, v) (X,Y), where the random variables X and Y represent feature vec- tors and target values from IR, respectively. Exemplary dis- tribution regression models include those based on mono- tone neural networks (Chilinski & Silva, 2020) and random forests (Schlosser et al., 2019; Hothorn & Zeileis, 2021; Čevid et al., 2022), and NGBoost (Duan al., 2020), which is founded on gradient boosting.
Copulas are instrumental in computing joint probabilities. For instance, the joint probability Pr(śi < Ti S Sa. T s (3) can be determined using copula as follows:
Pr(Si < T, S S4. T2 S Ca) = (C), C(E((()) (3)
where (x(()= = Pr(Ts 1).
In the context of survival analysis, a survival copula C is frequently employed, which satisfies the following equation:
An alternative approach involves directly estimating *s, (defined in (5)) through the utilization of density estima- tion model. Here, density estimation refers to estimating Pr(Y y|æ) from samples (x, iid. (X,Y). where the random variables X and Y represent feature vectors and tar- get values from a discrete set, respectively. Most multiclass classification models can be used for density estimation, in- cluding random forests (Breiman, 2001), 1, gradient boosting (e.g., LightGBM (Ke et al., 2017)), and neural networks. Recent advancements in density estimation techniques are found in (Filho al., 2023).
Pr(Ti > 11, T2 >
Tk > tx) -
= - F2(2)..... 1 (t)).
is well-established that any survival copula can be repre- sented using its corresponding copula C' (see, e.g., (Georges al., 2001)). For instance, the case where K 2, the survival copula C can be expressed
C(u1, u2) 01 + U2 - C(1 U1, 1 U2).
Additionally, is important to note that if = Cind (de- fined in (I)), then C Cind-
While broader array of models for density estimation exists compared to distribution regression models, an advantage of employing a distribution regression model is that, if we wish to adjust the hyperparameter B, there is no need to retrain the predictive model to estimate (t(1));; we only need to recompute Eq. (5) obtain fb
3. Two-Step Algorithm
We introduce a two-step algorithm designed to estimate the CDF F& (G) > for survival analysis with K competing risks. As illustrated in Algorithm 1, the initial step involves esti- mating (t the conditional k--th cumulative incidence function (CIF), which formally defined as follows:
3.2. Step 2: Postprocessing
The second step of our algorithm involves computing Fi (G utilizing the estimates bz obtained in the first
V (4)x)) Pr(T s t, k))z).
(4)
4

---

Survival Analysis via Density Estimation
step and specified copula C. Let Fbjæ E [0, 1] denote the vector of length K whose k-thh element is rb,k/zz and let F Falz E IRK denote the vector of length K whose k-th element is (C)
We begin by representing Tbj2 as a function of F 1/1e Fojz. and the copula C. For simplicity, we consider the case where K 2 this section, with generalizations for K > 2 detailed in Appendix B.1. By the definition of rb,k|za we have the following representations:
51
65
1b = Pr((s-1 < Ti S Sb. T - Tzl2)
Figure Illustration ((1)22 and 2(1.2)22
= 911),b/z - W1 ((122),b
(6)
2/ = Pr(Sp-1 < T2 E So S 21
Algorithm Two-Step (TS) Algorithm
1: Estimate V4 (t)xr)
012),b/z
W2 ((12),bb
2: Let fb = (& (x) (5e-6=1/2))
where
3: Let F
4: for b 0| 0
1.2,..., B do
Pr(Gs-1 < Ti S Cb. 66- S T2/x),
(8)
5: Calculate Fotz by solving Eq. (13)
Pr(Gi-1 < <T2 S Sb G-1 S Ti (x))
6: end for
7: return {F B-1
= Pr(Go-1 < Ti S Sb G-1 < T3 S Sb/z)
(9)
and W1, w2 2 0 are unknown weight parameters such that wi + W2 = 1. See Fig. 3 for an illustration of quantities (8) and (9) with b = 2. Unless otherwise stated, we set W1 = U² = 1/2. Note that, we use a sufficiently large B. the correction term (9) should be small, and therefore the choices of the weight parameters W1 and w2 should have minimal impact on these equations. Then, recalling be computed using copula
by using function gc that depends C.
Having established (13), we can obtain the estimation Fu of F3 for all b by solving this equation based on the esti- mation f3 of Taja: as outlined in Steps 3-7 of Algorithm 1. We leverage the initial condition F, Fojz - for all 2, where O is the K-dimensional vector of zeros. At the initial step for 6 = 1, we obtain F by solving (13). This equation is solvable because it contains K equality con- straints and the unknown value is the length- K vector Faa-- We repeat this procedure for = 2, ..., F, for all b. See Appendix B.2 for more details of the algorithm.
that any joint distribution (as demonstrated in equation (3)), it is straightforward to represent (8)-9) using Fi Fi(G6-1 (r), Fi (G.lz), (2)5)_1 F2(G)l)), and the copula C:
to obtain
9/1).6la = C(Pi(G))=)) C(F(())-))
Note that this second step of our algorithm is similar to the algorithm based on a bisection root-finding algo- rithm (Zheng & Klein, 1995), but their algorithm is valid only for K = 2 and its extension for K > 2 unknown. In contrast, our algorithm is extendable for K > as shown in Appendix B.I. We also note that as B - oc, Carrière (1995) shows another method to estimate F Fojz by solving equation similar to (13). We discuss this method further in Appendix B.3.
+ (3
(10)
4/2).6/a
+ (2
(11)
4(1.2).6|z
= C(Fi(G))Z), F2(6/zz))
-C(FF(G))))) F2(G6-1(z))
-C(FF(G--111) (2(G)lr))
Simplified implementation of our algorithm. In our implementation, adopted simpler approach to solve Eq. (13) for all b simultaneously, rather than sequentially solving Eq. (13) for each b as outlined in Lines 3-7 of Algorithm Specifically, by utilizing an automatic differ- entiation library for Python (e.g., PyTorch and Tensorflow), we estimate Fojz by minimizing the following objective
+ ( (x)), F (2 [z)).
(12)
Combining (6)- ing Fi (G)- Fi (Galar), F2(Go-122)) l F2 (2((),), and the copula C. This means that we can write this relationship as
represent F6, 1j= and 2
Tbl = 9C Fo F 11
(13)
5

---

Survival Analysis via Density Estimation
6
64
6
41
6
6
6
6
50
04
41
T.
7
6
6
41
4
6
44
6
62
G 42
Ti
7
71
Censored joint distribution
column-wise
column-wise sum
column-wise sum
Lower bound estimation
Copula- based estimation
Upper bound estimation
Figure 4. Illustration of the upper and lower bounds estimation. Here, region Rs. is divided into grids, with denser color indicating higher probability that a data point is contained in the corresponding region. The lower bound estimation achieved by assigning the probability mass Fa, the last time slot within region Ro, and then calculating the column-wise sum. Conversely, the upper bound estimation is obtained by assigning the probability mass fa,6 to the earliest time slot within region Rs, and calculating the column-wise sum.
function:
distribution (t((())) is completely known, Schwarz et al. (2013) discuss the identifiability of Archimedean copulas and the non-identifiability symmetric copulas.
gc (/oo/)) Fbjor)
(14)
It is important note that our two-step algorithm is suffi- ciently flexible to incorporate these identifiability results. Suppose that the distribution Fk.,(t]a) and the copula C are parameterized by v and 0, respectively. If the identi- fiability of the parameters and 0 is established, then the objective function (14) should have unique solution, where Foj and are replaced with the parameterized F b.nyla and Co, respectively.
for simultaneously. Note that Eq. (14) minimized if Eq. (13) holds for all b.
3.3. Non-Identifiablity
limitation our two-step algorithm is the necessity for prior knowledge of the copula C, which may not be read- ily available in practical scenarios. One potential approach to address this issue is to employ parameterized copula Cg. where f represents some parameter (e.g., the Frank cop- ula (2) with parameter f), and to extend the survival analysis addition to
4. Upper and Lower Bounds Estimation
As we discuss in Appendix A, there exist scenarios where it is infeasible to verify dependency in any manner and we cannot have any copula as a prior knowledge on the dependency. Even in such instances, our algorithm can still employed estimate the upper and lower bounds of Fa(())). By definition, the following inequalities can be readily verified:
framework to estimate the copula parameter the CDFs (t(z). However, Tsiatis (1975) demonstrated that this approach is not feasible because is impossible to identify (4(x)) without making some assumptions about (f)ar))
Consequently, many researchers have explored the introduc- tion of additional assumptions on the distribution (t]x)) to ensure the identifiability of both the copula parameter 0 and (t)x). For instance, Gharari et al. (2023) assume that F& (t]x)) follows Weibull distribution and present an algo- rithm estimate the copula parameter 0. Other examples include the study of copula identifiability under the propor- tional hazards assumption in (Heckman & Honoré, 1989; Deresa & Keilegom, 2024), and the investigation of copula identifiability for other restricted classes of distributions F& (t ( in (Czado & Keilegom, 2022; Wang, 2023; Zhang al., 2024). Under the strongest assumption that the true
Fx((//)) 2 Pr (t, E U Roc x
Fx(G//)) S Pr
(t,
Ro Ne I
which are equivalent
ro S FR(Golz)
Σ
(15)
6

---

Survival Analysis via Density Estimation
Based these inequalities, we estimate the upper and lower bounds of (2) by substituting Fb,A|z with the output fb,A|z from Step 1 of our algorithm. It is important to note that these inequalities are derived without utilizing the pa- rameters wi and W2 in our two-step algorithm. As illustrated in Fig. 4, the estimation of upper and lower bounds and the second step of our algorithm to estimate F& Fe(Gi (z) can be interpreted as redistributing the probability mass within the CJD representation into fine-grate mne-grained grid cells.
Subsequently, demonstrate the existence of strictly proper scoring rule for any K. According to Tsiatis (1975), the derivative un (t/zr) of the cumulative incidence function (CIF) (i(()) (as defined in Eq. (4)) can be represented using F& (4)))) and a survival copula C:
d
un (x(t)x) = di
Vi((l(z)
0(1 (1(1)))))
dt.
While these bounds for (1]r are computed only for dis- crete values t E {Go) it is possible to compute these bounds for any f using a distribution regression model to estimate the CIF (t)z)) (as defined in (4)). This approach explained in Appendix E. These bounds can be viewed as variants of the bounds established for the average survival function in (Peterson, 1976).
1 F2(t222),
1 - Fx(tk(z)
b==2=---
By using this equation, we propose a scoring rule NLL- SC, which stands for Negative Log-Likelihood based on Survival Copula.
It is crucial to distinguish that our upper and lower bounds differ from the concept of a confidence interval, which quan- tifies the epistemic uncertainty inherent in the prediction model (Bengs et al., 2022). Our bounds quantify uncertainty arising from the absence of prior knowledge about the true copula C. As discuss in Appendix E, these two types of bounds can be combined to account for both sources of uncertainty.
Assumption 5.2. The Kullback-Leibler (KL) diver- gence between ex (t(xr and its estimate (t)xz) satisfies DKL (xz(u=ll0 (0) < oo for all k, where ez and 0r denote the probability distribution over (ki, f) of ux (t/z)) and 0n (t|x), respectively.
Theorem 5.3. (A Strictly Proper Scoring Rule for Compet- ing Risks.) If Assumption 5.2 holds, the following scoring rule SNLI--SC is strictly proper:
5. Strictly Proper Scoring Rule for Competing Risks
SNLL-sc((FF((t)))) (t, = -1-= 8y (t)x).
Strictly proper scoring rules hold significant importance in the domain of statistics, particularly for the evaluation of probabilistic estimates (Gneiting& & Raftery, 2007). These scoring rules ensure that the expected score is minimized when the estimated probabilities accurately reflect the true distribution of outcomes.
Proof. Since Assumption 5.2 ensures the existence of the KL divergence, we have
(t)
-[SSuL-s(((PFA()))) (t,
In this section, establish the existence of a strictly proper scoring rule for any number of competing risks, K, given that the dependency structure is defined by survival copula. While Rindt et al. (2022) demonstrated the existence of a strictly proper scoring rule for K = under the assump- tion of conditional independence, no such rule has been established for K > 2 for cases where the conditional independence assumption does not hold.
K
(t(()) log (t(x)))
2
Equality in the last inequality holds and only if (t/xr) = (x(()) holds for all i and t, which is equivalent to that (e(()2)) = F& (f/xr) holds for all ls and t. Hence the scoring rule SNLL-SC is strictly proper.
We begin by defining proper and strictly proper scoring rules as follows:
Definition 5.1. (Proper and Strietly Proper Scoring Rules.)
A scoring rule S for the estimation (t(rr) of (t)x) is proper if the following inequality is satisfied:
Neural network model based on NLL-SC. Utilizing the scoring rule NLL-SC, we propose new monotone neu- ral network model, the survival copula network (SC-Net), Similar to the monotone neural network model for K = 2 proposed by Rindt et al. (2022), our SC-Net employs a monotone neural network to represent a CDF (Chilinski & Silva, 2020), and we use the NLL-SC as its loss function.
((S(E))()))) (1)))) (I((EE (t, 8))
(16)
A scoring rule is strictly proper if equality in (16) holds if and only if (t)()) is equal to (l(() for all and

---

Survival Analysis via Density Estimation
Table 1. Evaluation metrics
Task
Metric Name
Assumption
Metric Type
Density estimation
CJD-Brier
Discrimination
CUD-Logarithmic
Discrimination
CJD-KS
Calibration
Survival analysis
NLL-SC
Copula
Discrimination
Cen-log
Independence
Discrimination
D-calibration
Independence
Calibration
K.M-calibration
Independence
Calibration
We utilized the smooth min-max neural network (Igel, 2024) as the monotone neural network, though other monotone neural networks, such as those proposed in (Yanagisawa al., 2022; Kim & Lee, 2024), could also be used. Note can be viewed as incorporating
a = 0 to ensure its loss function was a proper scoring rule as recommended by Yanagisawa (2023).
Evaluation metrics. We used three metrics to evaluate the estimates of the CJD representation fb,A|a and four metries to assess the estimated distribution (t(()))) summarized in Table 1. For the CJD representation, we used the Brier and Logarithmic scores (Gneiting & Raftery, 2007) as dis- crimination metrics. Additionally, we applied the sum of the Kolmogonov-Smirno calibration error (Gupta al., 2021) as calibration metric. This metric, based the Kolmogorov-Smirno test (Kolmogorov, 1933; Smirnov, 1939), is defined as follows:
that our SC-Net with K prior knowledge of the ground truth copula into the DC- Survival model (Zhang et al., 2024), although the primary objective of DCSurvival appears to be the identification of parameters within an Archimedean copula.
6. Experiments
We conducted series of experimental evaluations assess the performance of our proposed two-step algorithm. We conducted the experimental procedures on a virtual machine possessing single CPU devoid of any GPU, equipped with a memory of 64 GB, and operating on CentOS Stream 9. The software implementation was achieved using Python 3.11.6 and PyTorch 2.2 The datasets employed for this purpose were the Dialysis and oldmort datasets, sourced from the Python package SurvSet (Drysdale, 2022).
hRo
max US
where
hk
and
N Σ J (e(
T
Models. In the first step of our algorithm, we utilized five different models. Specifically, the TS-Bries Brier and TS- Log models employed neural networks for density estima- tion, utilizing the Brier and Logarithmic scores (Gneiting & Raftery, 2007) as their respective loss functions. The TS-LGB, TS-RF, and TS-DRF models utilized the Light- GBM (LGB) model (Ke al., 2017), the random forest (RF) model available in the Python package sklearn, and the distribution regression forest (DRF) (Cevid et al., 2022), respectively. The prefix TS denotes the Two-Step algorithm, with each model implementing a distinct first-step model but sharing common second-step algorithm (specifically, the simplified implementation as detailed in Sec. 3.2). The hyperparameter set to B- = 32. Additionally, em- ployed our SC-Net model proposed in Sec. 5.
Here, fa (z) represents the probability of the feature vec- tor (**) being classified in class k, and each fx( is equal to its corresponding ft, For models that output only the distribution Fx (4))), the CJD representation fo was esti- mated using Eq. (13) with the output and the independence copula Cind (as defined in (1)).
For the estimated distribution (1)x)) we utilized the sim- plified variant of the censored logarithmic score (referred to as cen-log) (Yanagisawa, 2023) and our strictly proper scor- ing rule, NLL-SC, as evaluation metries. Additionally, D- calibration (Haider al., 2020) and KM-calibration (Yanagi- sawa, 2023) were employed as calibration metrics.
For comparative purposes, we included the Cox model (Cox, 1972), the random survival forest (RSF) (Ishwaran et al., 2008), and the DeepHit model (Lee et al., 2018), a neural network model. the DeepHit model, we the parameter
Results. Figure 5 illustrates the results for the Dialysis and oldmort datasets. The results indicate that while the DeepHit model demonstrated superior performance com- pared the Cox and RSF models, our TS-LGB model
8

---

Survival Analysis via Density Estimation
Dialysis
Dlayyus
Dalysis
Con
Goc
Con
Dialysis
15:08
15
15.089
15:0
15:008
15:099
15-nner
T5
Sc
SC
0.50
0.25
0.50
0.75
1.00
1:25 3.50
0.007
0.055
0.000
0,615
0.0
07
0.1
06
Cat-lng (mance)) 1)
Dicalomiiin Jeseet
dl bostoon (evert
Dislysis
Dialysis
Dia
Dialysis
Con
Coc
Doush4
Despht
Desphi
Do0/16
15-408
15:08
15-488
15:088
13.01P
15.0mP
15-00-
13:09
Log
15.0gg
15:brul
15.010
15.000
Sc-Nc
S0/4cc
1:5
2.0
30
25
0%
0.8
1
2.0
2.5
3.0
35
4.0
0.0
62
0/6
66
03
oldmort
oldmoot
aldmertt
Coo
RSF
050
Drepet
Deesit
oldmort
TScar
15
Tear-
rea
T5-4
TG
SC-lirt
SC
0.00
0.5
0.23
1.35 3.30
0.000
0.002
0.006
0.0
Con (overk 3)
D-caltration Jovor.
oVdmmrt
oldmort
alamart
Coe
Con
R0F
15p
Dogeit
Desse9
De0H41
Use
15:mm
15:0
15
15-M9
15
15-45
15:0g
SC
Schet
50 AM
Sc
MIL-CC
ryp-tree
CP-s
Figure 5. Performance comparison the Dialysis and oldmort datasets with B = 32 (lower is better). While Deepllit performs the best among the baseline models, our TS-LGB model shows comparable or better performances than DeepHit.
exhibited comparable or superior performance relative to DeepHit. It is also noteworthy that neural network-bass -ased models, including DeepHit, typically require longer train- ing times compared to tree-based models such as RSF and TS-LGB. Consequently, the TS-LGB model emerged as the most effective model for the Dialysis and oldmort datasets. Additional evaluation results using other datasets for K = 2 and scenarios involving competing risks with K- = 3, in comparison with models such as Deep Survival Machines (DSM) (Nagpal et al., 2021), DeSurv (Danks Yau, 2022), and Neural Gray (NeuralFG) (Jeanselme al., 2023), are detailed in Appendix F. Furthermore, the appendix in- cludes additional evaluation results encompassing the esti- mation of upper and lower bounds, as well as an ablation study on the hyperparameter B.
erates under the assumption of having prior knowledge of the copula C, which is a weaker assumption compared to the conditional independence assumption. This approach is consistent with Tsiatis's 's non-identifiability result (1975), and we have shown that algorithm also be utilized to estimate the upper and lower bounds of individual survival functions.
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
References
Source Codes. The implementationss of our models are ac- cessible http:://ggthbbbooom//hheeeeegett
Bengs, V., Hüllermeier, E.,, and Wacgeman, W. Pitfalls of epistemic uncertainty quantification through loss minimi- sation. In Advances in Neural Information Processing Systems, pp. 39205-29216 2022.
7. Conclusion
We demonstrated a reduction from survival analysis to den- sity estimation, which allows the application any density estimation model to survival analysis. This algorithm op-
Bilodeau, B., Foster, J., and Roy, D. M. Minimax rates for conditional density estimation via empirical entropy. The Annals of Statistics, 51(2):762-790 2023.

---

Survival Analysis via Density Estimation
Breiman, L. Random forests. Machine Learning, 45:5-32, 2001.
Fine, P. and Gray, R. J. A proportional hazards model for the subdistribution of competing risk. Journal of the American Statistical Association, 94(446):496-099 1999.
Carrière, J. F. Removing cancer when it correlated with other causes of death. Biometrical Journal, 37(3):339 350, 1995.
Georges, P., Lamy, A.-G., Nicolas, E., Quibel, G., and Roncalli, T. Multivariate survival modelling: A unified approach with copulas. Technical report, SSRN, 2001.
Čevid, D., Michel, L., Nåf, J., Bühlmann, P., and Mein- shausen, N. Distributional random forests: Heterogene- ity adjustment and multivariate distributional regression. Journal of Machine Learning Research, 23(333):1-79, 2022.
Gharari, A. H. F., Cooper, M., Greiner, R., and Krishnan, G. Copula-based deep survival models for dependent censoring. Proceedings of the Thirty-Ninth Conference on Uncertainty in Artificial Intelligence, pp. 669-680, 2023.
Chilinski, P. and Silva, R. Neural likelihoods via cumulative distribution functions. In Proceedings of the 36th Confer- ence on Uncertainty in Artificial Intelligence (UAI), pp. 420-429. 2020.
Gneiting, T. and Raftery, A. E. Strictly proper scoring the American Statistical Association, 102(477):359-377 2007.
rules, prediction, and estimation. Journal the
Cox. D. R. Regression models and life-tables. Journal of the Royal Statistical Society. Series B (Methodological), 34(2)187-220, 1972.
Greenwood, M. A report on the natural duration of can- cer. In Reports on Public Health and Medical Subjects. Ministry of Health. London: H.M.S.O., 1926.
Czado, C. and Keilegom, I. V. Dependent censoring based on parametric copulas. Biometrika, 110(3):721-738, 2022.
Guo, C., Pleiss, G., Sun, Y., and Weinberger, K. Q. On calibration of modern neural networks. In Proceedings of the 34th International Conference on Machine Learning. pp. 1221-1330, 2017.
Danks, D. and Yau, C. Derivative-based neural modelling of cumulative distribution functions for survival analysis. In Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, pp. 7240-7200555 2022.
Gupta, K., Rahimi, A., Ajanthan, T., Mensink, T., Sminchis- escu, C.,, and Hartley, R. Calibration of neural networks using splines. In ICLR 2021, 2021.
Defazio, A., Yang. X. A., Mehta, H., Mishchenko, K., Khaled, A., and Cutkosky, A. The road less scheduled. In Advances in Neural Information Processing Systems, 2024.
Haider, H., Hoehn, B., Davis, S., and Greiner, R. Effective ways to build and evaluate individual survival distribu- tions. Journal of Machine Learning Research, 21(1), 2020.
Deresa, N. W. and Keilegom, V. Copula based Cox propor- tional hazards models for dependent censoring. Journal of the American Statistical Association, 1 119(546): 1044 1054, 2024.
Heckman. J. J. and Honoré, B. E. The identifiability of the competing risks model. Biometrika, 76(2):325-330, 1989.
Drysdale, E. SurvSet: An open-source time-to-event dataset repository. Technical report, arXiv:2203.030%) 2022.
Hickey, J., Henao, R., Wojdyla, D., Pencina, M., and En- gelhard, M. Adaptive discretization for event predic- tion (ADEPT). In Proceedings of The 27th International Conference on Artificial Intelligence and Statistics, pp. 1351-1352 2024.
Duan, T., Avati, A., Ding, D. Y., Thai, K. K., Basu, S., Ng. A., and Schuler, A. NGBoost: Natural gradient boosting for probabilistic prediction. In Proceedings of the 37th International Conference on Machine Learning, pp. 1690-2700, 2020.
Hothorn, T. and Zeileis, A. Predictive distribution modeling using transformation forests. Journal of Computational and Graphical Statistics, 10(4):118-1119) 2021.
Emura, T. and Chen, Y.H. Analysis of Survival Data with Dependent Censoring: Copula-Bast Based Approaches. Springer, 2018
Igel, C. Smooth min-max monotonic networks. Proceed- ings of the 41st International Conference on Machine Learning, pp. 20908 20923, 2024.
Filho, T. S., Song, H., Perello-Nieto, M., Santos- Rodriig, R., Kull, M., and Flach, P. Classifier calibration: a survey on how to assess and improve predicted class probabili- ties. Machine Learning, 11::2211 1-3600, 2023.
Ishwaran, H., Kogalur, U. B., Blackstone, E. H.,, and Lauer, M. S. Random survival forests. The Annals of Applied Statistics, 2(3):841 1600, 2008
10

---

Survival Analysis via Density Estimation
Jeanselme, V., Yoon, C. H., Tom, B., and Barrett, J. Neural Fine-Gray: Monotonic neural networks for competing risks. In Proceedings of the Conference on Health, Infer- ence, and Learning, pp. 379-392, 2023.
Sart, M. Estimating the conditional density by histogram type estimators and model selection. ESAIM: Probability and Statistics, 21:34-55, 2017.
Schlosser, L.,, Hothorn, T., Stauffer, R., and Zeileis, A. Dis- tributional regression forests for probabilistic precipita- tion forecasting in complex terrain. The Annals of Applied Statistics, 13(3):1564-1588) 2019.
Kannel, W. B. and McGee, D. L. Diabetes and cardiovas- cular disease: the Framingham study, JAMA, 241(;9): 2035-2038, 1979.
Kaplan, E. L. and Meier, P. Nonparametric estimation from incomplete observations. Journal of the American Statis- tical Association, 53(282))457-481, 1958.
Schwarz, M., Jongbloed, G., and Van Keilegom, I. On the identifiability of copulas in bivariate competing risks models. Canadian Journal of Statistics, 41(2):291-303 2013.
Ke, G., Meng, Q., Finley, T., Wang. T., Chen, W., Ma, W., Ye, Q., and Liu, Y. LightGBM: highly efficient gra- dient boosting decision tree. In Proceedings of the 31st International Conference on Neural Information Process- ing Systems, pp. 3149-3157, 2017.
Sklar, A. Fonctions de répartition à n dimensions et leurs marges. Publications de l'Institut Statistique de l'Université de Paris, 8:229 231, 1959.
Smirnov, N. On the estimation of the discrepancy between empirical curves of distribution for two independent sam- ples, 1939.
Kim, H. and Lee, J.-S. Scalable monotonic neural networks.
In The Twelfth International Conference on Learning Representations, 2024.
Therneau, T. M. and Grambsch, P. M. Modeling Survival Data: Extending the Cox Model. Springer Science & Business Media, 2000.
Kolmogorov, A. Sulla determinazione empírica di uma legge di distribuzione, 1933.
Tsiatis, A. A nomidentifiability aspect of the problem of competing risks. Proc. Natl Acad. Sci. USA, 72(1):20-22 1975.
Lee, C., Zame, W., Yoon, J., and van der Schaar, M. Deep- Hit: A deep learning approach to survival analysis with competing risks. Proceedings of the AAAI Conference on Artificial Intelligence, 32(1). 2018.
Vallender, S. Calculation of the Wasserstein distance between probability distributions on the line. Theory of Probability & Its Applications, 18(4):784-786, 1974.
Nacini, M. P.. Cooper, G., and Hauskrecht, M. Obtain- ing well calibrated probabilities using Bayesian binning. Proceedings of the AAAI Conference on Artificial Intelli gence, 29),), 2015.
Wang, A. The identifiability of copula models for depen- dent competing risks data with exponentially distributed margins. Statistica Sinica, 33(2):983-1011, 2023.
Nagpal, C., Li, X., and Dubrawski, A. Deep survival ma- chines: Fully parametric survival regression and repre- sentation learning for censored data with competing risks. IEEE Journal of Biomedical and Health Informatics, 25 (8):3163-3175, 2021.
Wang. P., Li, Y., and Reddy, C. K. Machine learning for survival analysis: A survey. ACM Comput. Surv., 51()), 2019.
Wiegrebe, S., Kopper, P., Sonabend, R., Bischl, B., and Bender, A. Deep learning for survival analysis: A review. Artificial Intelligence Review, 2024.
Nelsen, R. B. An Introduction Copulas. Springer New York, NY, 2006.
Yanagisawa, H. Proper scoring rules for survival analysis. In ICML 2023, 2023.
Weed, J. and Berthet, Q. Minimax estimation of smooth densities in Wasserstein distance. The Annals of Statistics, 50(3):1519-1555
2022.
Yanagisawa, H., Miyaguchi, K., and Katsuki, T. Hicrarchi- cal lattice layer for partially monotone neural networks. In Advances in Neural Information Processing Systems, pp. 11092-11100 1103, 2022.
Peterson, A. V. Bounds for a joint distribution function with fixed sub-distribution functions: Application to com- peting risks. Proceedings of the National Academy of Sciences, 73(1):11-13, 1976.
Zhang, W., Ling, C. K., and Zhang, X. Deep copula-based survival analysis for dependent censoring with identifia- bility guarantees. In AAAI 2024, 2024.
Rindt, D., Hu, R., Steinsaltz, D., and Sejdinovic, D. Survival regression with proper scoring rules and monotonic neu- ral networks. In Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, pp. 1190-1205, 2022.
Zheng, M. and Klein, J. P. Estimates of marginal survival for dependent competing risks based on an assumed copula. Biometrika, 127 1995.

---

Survival Analysis via Density Estimation
A. On Applications of Proposed Methods
Survival analysis a crucial tool in various fields such as medicine, engineering, and social sciences, where the time until an event interest occurs is studied. The applicability our proposed methods can be classified into three distinct types based on the nature of the dependency between the event time and the censoring time.
(Conditional Independence.) The first class of applications encompasses scenarios where the conditional independence assumption between the event time and the censoring time valid or highly likely to hold. This situation typically arises in cases of administrative censoring. where data points are censored solely due to the limited window observation times. For instance, in clinical trials, patients might be censored the end of the study period regardless of whether the event of interest has occurred. In such cases, the conditional independence assumption simplifies the analysis, allowing the use of survival models based on the conditional independence assumption including two-step algorithm.
(Verifiable Dependency.) The second class of applications includes scenarios where there is a dependency between the event time and the censoring time, but this dependency can be verified, albeit at a significant cost, for a small subset of data points. example of this situation is found medical studies where the primary event of interest is patient mortality, and censoring occurs when patients are discharged from the hospital. In such cases, it might be feasible investigate the true event time for small fraction of patients, thereby assessing the dependency between the event and censoring times form of copula. This estimated copula can then be used to model the dependency in the entire dataset, allowing for more accurate estimation of survival functions.
(Unverifiable Dependency.) The third class applications involves situations where the dependency between the event time and the censoring time cannot be determined, even with extensive resources. such scenarios, identifying the survival function inherently challenging due to the unknown nature of the dependency. We note that, even without any knowledge of the dependency, method can estimate the upper and lower bounds to account for the uncertainty in the dependency structure. Additionally, if we have some confidence that the dependency can be represented by parameterized copula with some range of the parameters (e.g., the Frank copula with parameters -5 <. 0 S 5), we can also estimate the survival functions using the copula information to narrow down the bounds.
B. Step 2 of Proposed Algorithm
B.1. Generalization for competing risks
We generalize the second step of our algorithm presented in Sec. for K > 2. For notational simplicity, omit in this section.
For subset C [K]. let Q1, = {(t1,ty .tk., Awer(GG-1 < tw < (b) and Augr (G--1 S (w)). Then we can
compute the probability 9r.6 = Pr((t1, .,, QI,)) by using the inclusion-exclusion principle:
@T Σ
Σ K
J(I,J,b)
where b, C(PI.....P.J... .P..,,,,... ).) and
(1
ifke inkeJnI,
PT,
Fe(C)) Fi(Go-1)
if
Note that, if K = 2, the quantity (T.b here is equal to the quantities computed in equations (10), (11), and (12).
As generalizations of equations (6) and (7) for K 2, we represent rb, using qr
K
Tb.k = {(8),b
1-1)
Σ
WH HH
1=2
HAKECC
Correction term
12

---

Survival Analysis via Density Estimation
Algorithm 2 Algorithm to solve equations
1: Initialize F/ 1-11
2: while F/ is not converged do
3: for & E 21. . K do
4: Increase the k-th element of b so that the k-th equation of Fblu= 9c:(Folr Fo-1le) is satisfied (while other
k'-th element (k k') of F b/ is fixed)
5: end for
6: end while
7: return F
where WH 2 0 is weight parameter and we assume that wH = H/ Note that, if K = 2, this equation is equal to equations (6) and (7).
Since 2 are functions of F,, Fb-1, and copula C, we can write this relationship as
rb = Ic Fb-1):
by using a function 9c that depends on C. Having established this equation, we can estimate F5 by using Algorithm 1.
B.2. Solving Equations
In this section, present algorithm to solve Eq. (13), Algorithm 2 shows pseudo-code to solve this equation, and it capitalizes on the following property:
Property B.1. Assuming that Fo-1|= is fixed:
(1) The k-th element of the length-K vector (C In Fo-1)p)) 1|x) is monotonically increasing with respect to the k-th element of F M)..
(2) The (+ k))th element of the length- K vector ge (F Fo-1w) monotonically decreasing with respect to the k'-th element of May
First, we demonstrate that the following inequality always holds during the execution Algorithm 2:
Tblz 2 gc [F 1/2 Fo-11))
(17)
At line of Algorithm 2, Fo is initialized with 1/-1/-r Hence, by the definition of Eq. (13), we have go (Fo F
which means that inequality (17) holds. At line 4 of this algorithm, we can increase the k-th element of F Falz to satisfy Paja 9c b/. FF, )11))) due Property (.1 and this increment does not violate inequality (17) due to Property B.1
Since each element in the length-K K vector Film does not decrease during the execution of this algorithm, we can find the solution to equation (13) by repeating the while-loop (lines 2-6) until the convergence of F biu-
Parameterization.
When formulate equation (13), we assume that F& (Gola))) represented by using some parameters. In our implementation, we used the parameterization by using the softmax function
Eb exp(av.k.t
Fa(Co/2)
ΣB= exp(av,k..
with parameters kg, for each and æ. The softmax function was chosen to ensure that any function F& (c can be expressed. Note that, any other parameterization can be used to represent F& (Cul)))) in our algorithm.
13

---

Survival Analysis via Density Estimation
B.3. Alternative Algorithm for Step
We briefly explain the algorithm presented in (Carrière, 1995). This algorithm exploits the following equation:
(t) = 01k
...., 1- Fx((K))
① ((S) Sk(tk)) 11
0tk
8C(ui, u2, uk)
Ss((x)
Oug
(()====((0)....xx=kkk dt
14=1
where C the survival copula corresponding to copula C. Assuming that B is sufficiently large, his algorithm solves this equation with respect to (t) on & {G))= by using estimate Vi)(() of (t) and these approximations for all k:
S&(Co) + Sk((6+1)
uk ~
2 Sa(Go+1 - Sa((s)
d
(1)
dt
66+1
t=
d
Va(Go+1) - (Cn)
(t)
df
S0+1
1=
This algorithm is designed to estimate the average survival function (t)(t) and thereby (t),,, but we can modify the algorithm to obtain Fx (t) conditional on x if a conditional estimate (t is available.
C. Theoretical Analysis
In this section, we theoretically verify that the two-step algorithm outputs solutions with sufficiently small errors. We consider the case K 2 for simplicity, and assume So BB. As discussed in the preceding section, various models can be implemented to estimate the CIF in the first step of our algorithm. Therefore, we evaluate errors affected by Step 2, solving (13), under the assumption that the models employed in the first step accurately approximate the true probabilities such that
A Tb,kl S.
(18)
holds for all b = 1..... B and 1, 2. Note that while how small a value we can take as in (18) depends on the choice of the model in Step we can apply the results exhibited in this section. We provide examples of achieving (18) in Appendix D. To formally state our theoretical results, introduce the following assumption:
Assumption C.1. We assume the following conditions:
(1) (True probability is not biased.) There exists a global constant co > 0 such that for every b = 1,..., B and k = 1, 2,
(G(()()) - FA (66-1/2) = Pr((6-1 < Tp I (6) < B
holds.
(2) (Copula.) We assume that the copula of class C² and satisfies
a²
i = inf
C(u, > 0,
(0.1) dudv
8² v),
8 dudo u(u,v, A(u, m)) < +00.
8²
L = sup max (v.r)e[o,11
Ou
(3) (All tu are equally observed.) There exist constants C1 > 0 depend on e, L, and T > such that the following condition holds: Let 6, > 0 be a constant determined by l and I. and by = max b Vk, (Cola))) S min Fg (x)) 2 Tog B
8n
Then,
b
+ log 1
d1
k
We make some remarks on the assumption. The first condition required to bound the error by the choice of W1 and wz: the probability concentrates on a squared region partitioned by suboptimal wh and u2, significant errors are inevitable.
14

---

Survival Analysis via Density Estimation
This condition is satisfied if Pg is Lipschitz-contimuuus with coLB serving the Lipschitz constant. The second condition manages the sensitivity of the estimation relative to the true distribution and noise. For example, if &CC < 1, indicating Budv
that C' exhibits minimal variation as u and v change, substantial adjustments to the estimation are necessary to accommodate for noise and achieve (13). This condition is typically met for the independence copula C(u, v) = we with any < and <>1 1. The third condition appears to be technical. As will be demonstrated in subsequent analyses, errors between Fi: and FR can only be effectively bounded for b S bo- As approaches B, and consequently r diminishes, the impact of intensifies. Condition (3) excludes scenarios where a part of f&s is concentrated in the region b > bo- In other words, all tys are equally observed in the region b < bo-
Let W1(..) -) be the Wasserstein distance', Then, we provide the statement about the W1 distance between the estimated and true probabilities. We consider the extension of Fs (Gala) to a CDF on [0, (B) by F (z) = (x(GG( where (o < S Cb+1-
Theorem C.2. Suppose that Assumption C.1 holds. Then, there exists constant Ce > 0 depending co, l and L such that if c < ce the following inequality holds:
B.
B B-
W1 (PMpes S (B B B+
(19)
B Blog B
where Pi and Pklz are probability measures whose CDFs are given by (x(()x) and Fx respectively.
The proof is deffered to the following subsection. Suppose that the condition = oB + holds. Then, we obtain an upper bound as Wi (ix|e= MMtee) = Ça o(1) Thus, we can ensure that as B - +0o and the sample size increases as we can take sufficiently small e, the output of Step 2 converges to the ground truth distribution in terms of the W1 distance.
We provide some comments on Theorem C.2. We can observe trade-off in (19) based on the choice of B: while the second term decreases as B increases, B++T and in the first term should increase. Consequently, an optimal choice of B should be considered under appropriate assumptions that determine and bu, such as the model utilized in Step 1 and the properties of F- It is also significant to examine that the derived bound achieves a statistical min-max lower bound exhibited in (Niles-Wer & Berthet, 2022; Bilodeau et al., 2023), for example. We reserve these considerations for future research endeavors.
C.1. Proof of Theorem C.2
This subsection provides proofs of Theorem C.2. Here, for notational simplicity, we abbreviate the conditional variable I. For example, we denote F& (C6) instead of Fe(Gsl))). We assume w1 = U12 = just for simplicity. We can extend our analysis to arbitrary w1, w2 E (0, 1). Moreover, we assume that Step exactly solves (13), i.e., Fi. exactly satisfies the equation (13).
First, we evaluate the error between the outputs of Step 2 and the true probability.
Proposition C.3 (Estimation error when solving (13)). Suppose that Assumption C.I holds. Let Fb be the output of Algorithm 1. Then, for every b = I....,B B and k = 1, 2, there exists a constant C2 0 depending Co. l and L such that under the condition c < -- for sufficiently large B and every b satisfying max Fs(C.) < 1 - 5 with a positive constant &
kE(1.2)
&egL satisfying Fog
S ô with T E (0, 1)
(R(GG)) Fr(Gsa) 4
16L8-2 12
(e
(e+
I.
}
(
B 4coL.
)
(1
ScoL
bb+1
(20)
BL5
Proof. Let us denote Ao, = Fa(Co) FF(Co). Instead of (20), we aim to obtain tighter bound
161.5-
(c (C+L. (.1..g))
])9 [
(1- )
))
8c0 T.
-0+1
A S c+ l2
(21)
4cgL)
Béd
4cpL
"W,(p,))= inf Sealt yldn (2. y). where 11(H, v)) denotes the set of all couplings of two probability measures μ and y on rell(urr)
R.
15

---

Survival Analysis via Density Estimation
We give its proof by induction on b. For the case have
(G), 1) = F1, + wi ((()
and
1) = wi, C(F((G),
with wf, [0. By taking the difference of the both sides, obtain
1) - C((L()))) = 61,1 -1,1 wi C(A(G)) wi,1 C(A()))
(22)
First, we evaluate the term (S)), Rearranging terms gives
1) - (C) P((1)) = -1, - -I,1 + 1(F(()) "in ((R)
6 ((F((())))
co B
sct L
where the first inequality follows from f1,1 -1,1 S c (18) and wi,1 C(FI(G)), F2(C.) 2 0, and the last inequality is derived from ((1) < da by Assumption C.1-(1) and B
C(R(G)) Fi(C.) J.'
82
dudv C c( dedu s L. F(() < L. ca
SL
Since ((F((S), & C(F)(G1), 1), we obtain C(R(GG) -
- w1 C(C())) 5 L.
Co B
By using 1 - w1 1.2 we obtain
( CD
B
Moreover, have
8²
C(F((G))
Clu, v) dodu 2 (F)()))
Oudr
2
and hence,
2 c c0) B
Fi(S)) < 1
A similar argument gives
2 ( -00 B
(a(C1) < 7
By combining these bounds, we obtain
A(5i)
{
Az(S1)
82
4L
C 2
A(G)))
C(u, v) dudu L. (LG))BG)
<
c+
L
(23)
Oudo
B
1t
Thus we get the bound on the term (
16

---

Survival Analysis via Density Estimation
Moreover, have
((P(G)))) Fi(())
FalCi) !
0²
v), dedu S ((()))((())) SL (L.())
)
(24)
dudr
B
ST.
where we use Assumption C.1-(1) for the last inequality.
Then, by taking the absolute value of the both sides of (22), we have
1)
= A.. 10.--,, F1,1 + wi C(FI(G)..B((II) wi, ((F((S),
- f11,1 - 1111-1,11 w1 (2R) - wf,1 C(Fi(((),, F((1))
S 1++mmx(((E B(LG))
44 E L.
60 2
C0
21
S s max
L
1
B
B
Set 4L
C0) 0o)
2
L
where the triangle inequality for the first inequality, 0 S W1, w S for the second one, and (23), (24) for the third one. Since (C(A())). - C(Fi(C)).
1) 2 / A(( Fi((1) A111 we obtain
A - cl 4L
27 S c E (E+L )
16L8 2
(e E L.
90 2
B
where the last term coincides to the right hand side of (21) with b 1.
We can obtain the bound for A: by utilizing the same argument. Thus we get (21) for b 1. Assume that (21) holds for = 1/, i.e.,
A S + 16L5-2
(E+ 10 3)
B AcgL
)
8c0L
-b'+1
B 4coL
BPd
holds for l 1. We consider the case b = b' + 1. This bound gives that by taking
e with sufficiently small c.,
-6'+1
A 5
16L8-2
T,
Cu B
B
ScyL
12
8
4.oL.
B88
16L8-2
Cu
1
B
8c
5
c L.
+
B
8
4cgT.
BLd
1615-2
4
(1
B
8eg T
(+ L,
+ 4coL
exp
1
B
85
elog B
CD 27
]]
(1
B
(c+
B
(25)
2cB
B
4coL
2 BI-T"
(26)
where we use - x < e in the third inequality and the definition of 5 in the fourth inequality. This and Fs((w) < 1 - 8
gives Fx (Cw) < 2 for sufficiently large B.
We remind
f4411 = -(11,N w1 Q(12),.
(27)
TV41,1 -(1).8 - ((12),
(28)
17

---

Survival Analysis via Density Estimation
Now, we consider integral representation of 411,./ and (1122,00
11)000 1())(()) C111111
0² dudr
(29)
Q(1.2).8 C(Fi(GV+1) Falle)) -CAA://HHHHGGNAA/NCCCCCCCCCCAAA- (ilGort))
(altor)))
8²
(AIGW) Palcs)
Clu, vydudu.
dudv
The same expression holds for 9(1),0 and ((123,1 by replacing Fn and R2 with Fi and Fg. By taking the difference between (27) and (28), we obtain
FB+11 - 11+ = ((1)v Q1). wi (1.2).0 + i ((12),
(30)
Similar the case b first evaluate the term ((1.23.0 (note that ((12).1 (((G))) By rearranging terms, we have
9(1), w1 -(1,2),1 fv++111 - + 9(1)- wi,1 ((1.)),
s 9(1).8'
= [ (m(Caw Fa(())
(
8² Clu, v) dudu
Iudu
SL
Set L Set L (F((S+41) Fi (L(F(())))
Co B'
(31)
where the first inequality follows from fM 11+1111 S by (18) and 0(1.2),0 2 O, and the last inequality, follows from Assumption C.I-(I) and 1 F26()) < 1. Moreover, the left hand side is lower bounded by
((1). wi ((12).. 2 (1 - w1)6(19M
1 (ilGuee)
0² dudu
2 Jecco
Clu, v) dodu
2t
2 2 (AA(GG+4))) - A(Go)) F2(C))),
where we use 1 - w1 = and (29) for the equality. Combining this with (31), we have
F - Fi(Sx) -
2
E L. S0
(32)
e(1- F2(Gv))
This and the triangle inequality give
A =
(P1(G++1) - P((())) + (Fi(Ge) - Fi(G))) + (Fi((C)) - Fi (S(+411)
< -1(5)+1)- - F(GG +|Ar(G)) - -FI(SE)) + +[F((Go) - Fi(G++1))
< 2
( L c0) + | B'
cu
8(1 Pa(Gx))
which we will use in the latter of this proof.
The same bound as (32) holds for = 2, which gives
(Ri(Gr+)) (G())) (F2(Gor+)) f(((Gr) S
( ( L. C 20)
1 (1(1-f((5)) Fi(Gv)) (--F)((()))
B
18

---

Survival Analysis via Density Estimation
Thus, obtain
8
11.23.0
AllW))
B2(33)
Oudo
C(u, dudu
(33)
SL
L1 -
(F2(5041) - F2(S5r)
4L
( (et
-0
(2(1 11 - (1-Ff(ss)) (1-F((sw))
B
Moreover, we have
11.2),001 = (
(
8² Oudo
JPASS
C(u, v) dedu
(34)
St.
< L. (R((Gr)))
- -F((e))))))
$ Co 2
where we use Assumption C.1-(1) for the last inequality.
Then, by taking the absolute value of the both sides of (30), we have
(0),81 9(1),111
| TV+11 W1 (112), 16+11 ((12),112
< |fv41,1 - 1841111 + [w ((1,),,61 - w ((1.2),11
< + max (axx((@11100 (Q(1.2) (1.2),1
4L
E+ L. Co 00
Cu 2
< + max
L.
12
Filser) (1 (e(1-Filsw))
Fo(sar)
B
B
(I-R( -FalGr) (c+ L. 20)2
4I.
set
(35)
where we use the triangle inequality for the first inequality, S W1. 464111 S for the second one, and (33) and (34) for the third one.
Then, evaluate the left hand side. We have
|AAcson Jeucs Fudr Clu, vidudu
82
FilGo++1) p1
82
[9(1)--- 900
FFlGV) Jatc dudr Clu, ddddd
19

---

Survival Analysis via Density Estimation
and
Fillo41)) Fi(So))
S
0² Clu, v)dudu
&² JESCS audr Clu, v)dvdu
JES(C)) dudv
(Pllor CCu, 1)du
FilGor 0
= -J
C(d,
du -FIlCN) ( Fa,, Falss))
(((CC)
du
Bi(G/++) a F2ls))d du
du
P((c))
Folloete) á ((u,1)du- Filla)
-JAICOL) Ōu
J Su
C(u, 1)du
Ai(GN) 0 (aGr))
(filbw) Clu, F2(Se) d)
JAIGGAai du
+
P. (Sr44)) du
0
Clu, (2(x)) du
Ou
Ou
FilGs)
C(u, 1)du
C(u,, 1)du
Du.
Ou
FilGu+t) 0
FilGy) θ C(u, F2
C(u, ((u,F2((G)))d
(36)
Ou
PAlG du
AllGu+1)
0
C(u, F2(Gr))
Clu, Fa((sr)) du.
Ou
Du
By the mean value theorem integrals, there exist constants
min (i((s). max max(AA(s))) Fi(())
uy 1. 2 E
and
min (1) max (G)
W4111.44+12. E
such that
A C(uv411;1) Av,1 Du
C(u,, A841.1 2 C(u Fe(Gx+1))
(36) A64111 du
Du
0
Fille)s)) C(u,
0 Clu, Fe(sx)) du.
A2+1,1 du
C(ur,2 C(uv))
+
(((CV)
Ou
du
We denote
(I) = A84111 C(u - AN. 2
C(w,11
du
In.
AV-11 C(u6441,2 F2(Gv+1)) + AUt11 a
du C(uy,2,
du
A((r441)
(II) =
C(u, (2lGr) - C(u, (((.F)(()))) du.
Aisw)
Ou
Du
20

---

Survival Analysis via Density Estimation
We evaluate the each term. First, have
θ
i C(u
(I) =(A++11 1v
C(w 1. 1) A8+11
c(uy, ((w...1)
du
du
du
a Cluv 2, F2(Gu))) AM+1,1
O
8 Fels)))
(Av+11 Av11)
(( F2((o++1)) du
du
du
θ C(wy,,1)
0 C(urt1a,1) 0 ((uy,,1))
=(A++1, 1v.1)
c( F2(Sor)) +10+11
du
Ou
du
0 C(ur++12: F2((v+1))
0
Cluy F2(Su))
Av+111
Du
du
Thus, we obtain
((D) 2 [ - Av.1)
② C(usr,,1) Clu 2, Fa(sc)))
du
du
((w
Cuu,,1 1)
+1+1,1
C(u F2(Gr+1)) Ou
C(uy, Fe(Ser))
du
where we use the triangle inequality.
Moreover, we can bound the terms in the above inequality by
C(1 C(ux,1.1) =
u6641,3 02
Clu, du
du I L(w0+111 -,,,) s <I(Iav=aal+ + )
tu
du
Ou
and
2 F2(S+41))
Cluy 2, F2(Sx.))
Ou
Fu
<
C( F2((G(+1)
CCC Fe(Gor))
+
C(wx+11,22 F2(Cor))
CCCC F2(Gs)) ))
Ou
Ou
Ou
Eu
(Pilloss) 8²
+ f
ÿ2
FelGGw
C(u du du
Olu Fe(Ce) ) du
O²2
1o Ou
SL
SL.
Co
Co
+Avtaal)= = I 10---11AA (1Autaal+
200
S L.
+
+
+
B
B
B
where we use Assumption C.1-(1) for the first term and |usr2 uv4122 S A + co for the second term. Moreover, we have
(mI)
( dude Clu.v) dudu
82
<
L (c L C0 Av.2l S
4c
B&8 Ar
2o
ZL
where we use e < B in the last inequality.
Then, (35) gives
1 - Aval-
C(uy,1, 1) 1) CC(uu 2, F2(()))))
Du
Ou
4L
2 2
4c
3co (1w+tial=
5
(e+ b.
+
Ar.2 + M/
(2 (P1--i((w) () F2((v))
B
Bed
B
21

---

Survival Analysis via Density Estimation
Moreover, by the triangle inequality, for sufficiently large B satisfying L- B
S
cb
2
we have
0 Cuuy 1, 1)
0 C( F2((())
Iu
du
i Cluy,2,
0
C(uy 2, F2(Su))
((ur1,1)
C(uy 22,
du
du
Ou
Iu
8 Cluy 2: du du
0² Cu, 1) du
Ou2 (FalGu))
du
28
27
2 E(1 F2(Sr)) T/ur,, uv 2 88
I.
C0
2
f6
B
2
where the third inequality follows from F2(G6) S 8 and
1u.2 20,1 I max (i (i((v))
A $
Co
by (26).
Then, have
Av+1.1 Avall
2
41.
20)2
4cgL.²
(Aassal+
300
c+
(CLL L.
+
A + L|
88
8(1- (l-f((cr) (1 Fe(Ger)
B
B&8
B
4L
00)²
Co 2
4c
2L F/b
+ 1
3c0
<
E +
E
+
BBB A
+
85
82 (11-
(-(G) (1-))11-11 R2(Grr)
B
e5
B
We consider two cases (i)
A84111 S o and (ii) A841,1 >
B If (i) holds, we have
2 to
4L
C0 2
4caL²
8cpL
A - Ar A S
c @2 (2(1 (2(1- Fi (Gx))
( L.
+
Av.22
+
A Asttal
-1
(2(r))
B
BPd
B18
By using the triangle inequality again, we obtain
2
4L
-
4co
ScoL
A AI. + 4o
[ 0 e011- Rr (11-A(G))) -1-())))))
c+
Pe(Gw))
c [ L.
B
+
Bed
Avv
+
Bed
Avt1al
Finally, by rearranging the above inequality and using the induction hypothesis, we have
(1-
11-
| &
c+
16L5
[-L.
1)
c 21
B [(t+
1
80q
-1+1
B
BE8
B
B
4coL
BE8
4c
2
4L
Co
+
6
(e+ L.
85
(1- (1 Fi(())) (1-
12
B
16L8-2 2
Co 27
B {(t
4CoL )
B
(1
ScgL
-b++1
B
<
+
( L.
BEd
4
2
16L8-2
(E+L) L. 200 }
c+
88
22

---

Survival Analysis via Density Estimation
where we Ft(() < 1
_/2 for the last inequality, and hence,
16L0-2
2
1coL ) (
)
B
8co
-1+1
B
(1-
8cp L
A041n1 S c+
(++
L
e2
Blö
4coL
BL
(
800
-1
2
16L6-2
Co cu)21
+
c
6+ T..
BLd
f2
B
+
16L6
-2
( c L
Cu
27
+
B
1
8coL
-1'4+1)+1
B
(2
B
4cgL
BLd
L
which ensures (21) for b- = b' 1 and k 1. case (i).
Since the same argument holds for k = 2, obtain the conclusion for the
If (ii) holds, we have
A - Aw.il
2 E8
6
AL
E c L.
C0 2
+
4coL2
+
8L A
A
P(1- Fi(Go)) (1 (1- F2(Cv))
B
BE5
Then, by using (26), we have |B+1,11 1- B (2-((1-7), This and triangle inequality gives
2
AL
00)
3coL²
8L
A < Av +
c L.
+
Ax22
+
A Av
E8
12
(1-M(CH))
(1 P2(Cr))
B
Bld
26
< AAv +
(1-
8coL
-1
2
+
16L
(e+
(e+L.
L.
CD
40,L2
Avr
B8d
td
@ (2 (i (1-
B
Biö
with taking sufficiently large B. This gives the conclusion for the case (ii).
To obtain the bound with respect to W1, we utilize the following lemma:
Lemma C.4 (Vallender (1974)). Let and v be probability measures on R whose CDFs are defined by F and G. Then,
[F(t) - ((t))- W2(,vv).
Then, we move to the proof of Theorem C.2.
Proof of Theorem C.2. Set dy in Assumption C.1-(3) by 6, = ScgL
By using Lemma C.3, we obtain
+oc
(F FR) =
(t) - FF(() dt
[a(()) (t( dt
b0
Re(1) Fic1))at+ ar+
(x(1)) Fa(1))ar.
(37)
(m)
where bn is defined by Assumption C.1-(3). Then, we bound the each term.
We denote AD.E = F((Gs) - Fk(G)) again. By using the bound (25), we have
A S [e log B
Bl+r
B
B2
23

---

Survival Analysis via Density Estimation
Then, the first term can be bounded by
bn
6
[Pc(t) P Foccs))
(1)
5
max
b=1
<h
max maxx|Ab,ll Ab, dt
b=1
b0
Sh
(AaAll +
C0
dt
B
b=1
log B
SB
B1+T.
c
132
B
b=1
log B
SB
log B
E
B²a
B²+T
B
B²
BL
(38)
For the second term, since 1 - (G;)))
2 log 1 B by Assumption C.1-(3), we have [Pa(1) FF(()) x log-11 B for t 2 Sbo
with sufficiently large B. This implies
B bo
(LI) x
log
Bdt -
B Blog B CB.
(39)
By substituting the bounds (38) and (39) into (37), we obtain the conclusion.
D. Examples of Bounds on Step 1 Error
In this section, we present an example that satisfies (18). As mentioned in Sections 3 and C, any estimation method for the probability distribution can be employed, and theoretical results pertaining to those models can be leveraged to guarantee (18). Among the various methods, we introduce results derived from the Distributional Random Forest (DRF) (Cevid et al., 2022) and histogram type estimators (Sart, 2017).
D.1. Distributional Random Forest
DRF constructs random forests designed to estimate the conditional distribution of multivariate responses. It achieves this by splitting the data using distributional metric, specifically the maximal mean discrepancy (MMD), with the goal of maximizing the differences in distributions between child nodes. DRF then estimates targets, such as the CIF in this study, by employing weight function that reflects how frequently the training data points end up in the same leaf as the test point across different trees.
Čevid et al. (2022) impose the following assumptions:
(P1) (Data Sampling.) Instead of the traditional bootstrap sampling with replacement, commonly used in forest- methods, a subsampling approach is employed. For each tree, random subset of size 8n is selected from n training data points. It assumed that sn approaches infinity as n increases, with the rate specified below.
(P2) (Honesty.) The data used to construct each tree split into two parts: one part is used for determining the splits, and the other is used for populating the leaves and thus for estimating the response.
(P3) (a-Regularity.) Each split leaves at least fraction < a S 0.2 of the available training sample on each side.
Additionally, trees are grown until each leaf contains between k and 2x 1 observations, where E is a fixed tuning parameter.
(P4) (Symmetry.) The (randomized) output of tree does not depend on the ordering of the training samples.
(P5) (Raiton-Spiit.) At every split point, the probability that the split occurs along the feature Xy is bounded below by π T/P. for some æ > 0 and for all I.....p.
24

---

Survival Analysis via Density Estimation
Note that each of these conditions can be verified by inspecting the constructed forest. The following proposition is direct consequence from Corollary 5 of (Cevid al., 2022).
Proposition D.1. Under the assumptions (PI)-(P5,, it holds that
fb, 4 Tb,k|z
for any b and l as the sample size goes to infinity².
The proposition above ensures the probabilistic convergence of the estimation. Specifically, for arbitrary values of e >> 0 and 5 >> O, there exists sample size threshold ne,5 > O such that if the sample size exceeds n6,6, then (18) holds with a probability of at least
D.2. Histogram Type Estimators
Sart (2017) proposes histogram type conditional density estimators on X' x y by
Σ 1x v1)
Pr(yl=) 8(x,1)
1x (x,
Σ (6 ()(K)
KEm
where m is a partition of X' x V/ is a reference measure of the conditional density, and óz is the Dirac measure at a E X. In the context of survival analysis, we can choose y = [0, 1] equipped with the Lebesgue measure and m as a set of regions defined by (So-1 < tx S 6,, = k) for E [B], k E [].
Let be measure defined on and
= SX-D (VFle.u) -Vol,,) xxY
h(f.g) =
Vo(z,v) v) dv(x) du(y)
be the Hellinger distance. Then, Sart (2017) provides the following result: Proposition D.2 (Proposition 2.6 of (Sart, 2017)). Let be true conditional density. Then, there exists global constants C1. C2 > 0 such that for any 0.
[m|
Pr h2 (s. inf (2(s, v) +1 vEVn
n + +C2 2 n
where
Vm Σ axlk, AK E a, 20
See (Sart, 2017) for specific examples of deriving the term inf 1i² v) under the conditions on and X. The bound on the VEVn
Hellinger distance implies the bound on the total variation distance TV(.,)) g) SxxylS(s, g(x, v)) du(x) du(y)
TV(f.9 S h(f.g),
which follows from the inequality between the Lnmmorm and the L2-norm. Thus, by utilizing Proposition D.2, we obtain (18) by taking e as the total variation distance between f and r.
E. Upper and Lower Bounds Based on Cumulative Incidence Function
If the CIF (4(ar) (as defined in (4)) is available, the upper and lower bounds of F& (4[ar) can be easily computed by
Pr(T s t, 5 = k|z) = Fx(t)z) S Pr(T $ t|æ)
[i((|>>) < (1(() S Σ Vv (t
2.P denotes the probability convergence.
25

---

Survival Analysis via Density Estimation
Table 2. Real datasets used in our experiments
Name
K
N
# categorical
# numuerical
censored
time
dataDIVATI
5943
83.6%
6225
oldmort
6495
69.7%
20
Dialysis
6805
2
76.4%
44
flchain
7874
4
72.5%
5215
support2
9105
11
24
31.9%
2029
prostateSurvival
2
14294
3
94.4%
119
PBC
3
312
12
45.8%
15
Framingham
4434
10
56.2%
8767
By averaging over ~ X, can derive the same upper and lower bounds as in (Peterson, 1976):
Σ
valeeo))
E [[((t(())]] Fi()) < x~X
(40)
x~x
Given dataset D = {(x(0) (00,80)) ) we can compute the empirical estimates of the expectations in (40) as follows:
E [+((|z)]
~
Σ
1,0
x~x
N
(aft) ,06)) (260)
(80)ED
E X
Vi(((z)
(af1),8t)) (269) (800)DD
Note that these values are equivalent to empirical CDFs, and they may not correspond to the actual bounds if the number of data points is insufficient. In such cases, the confidence intervals of these empirical CDFs should also be computed using methods such as Greenwood's method (1926).
F. A Additional Experiments
Datasets. We used eight datasets, summarized in Table 2, where N denotes the number of data points, and the fourth and fifth columns indicate the numbers of categorical and numerical features in the feature vectors, respectively. The six datasets with K- were obtained from the Python package SurvSet. (Drysdale, 2022). The Framingham (Kannel & McGee, 1979) and PBC (Therneau & Grambsch, 2000) datasets with K 3 were ones used in (Jeanselme al., 2023).
All datasets were randomly split into training (65%), validation (15%), and testing (20%) sets. The results reported in this section are the mean and standard deviation over five random splits. We divided the time horizon [0, {max] into B - 1 evenly spaced boundaries and added an additional time slot to represent times greater than fmax.. where tmax is the maximum observed time within the dataset. This setup means that we used time slots divided by { We set B = 32 unless otherwise stated.
Models and hyperparameters. We used a multi-layer perceptron (MLP) with three hidden layers as a neural network model. The dropout layer was employed with dropout rate of 0.5, and the ReLU function was utilized as the activation layer. The softmax function served as the output layer. The neural network was trained using the AdamWScheduleff optimizer (Defazio al., 2024) with early stopping. For each dataset, we performed a hyperparameter search determine the number of neurons in the hidden layers and the learning rate of the optimizer: the number of neurons chosen from the set (4, 8, 16, 32, 64, 128, 256), and the learning rate was chosen from the set (0.00001, 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0).
For TS-RF and TS-DRF models, hyperparameter searches were performed on these parameters: n n.est immt was chosen from the integers between 100 and 1000, max_depte was chosen from the integers between 10 and 50, min.samples.s camplesssppitt was chosen from the integers between 2 and 64, wain_samples.l integers between and 32, max features was chosen from sqrt or log2, critert on was chosen from log loss,
was chosen from the
26

---

Survival Analysis via Density Estimation
gini, or entropy, splitting.rule was chosen from CART or Fourie num_feature was chosen from the integers between and 100, sample fract on was chosen from the numbers between 0.1 and 0.5, min node s ze was chosen from the integers between I and 10, and alpha was chosen from the numbers between 0.01 and 0.3.
For TS-LGB model, hyperparameter search was performed on these parameters: integers between 30 and 300, max heeth was chosen from the integers between 1 and 64, min n.chiiiiiiii 1 samples was chosen from the integers between and 32, num_leaves was chosen from the integers between 2 and 64, criterion was chosen from log loss, gini, or entropy, learning rate was chosen from the numbers between 0.001 and 0.01, lambda.11 was chosen from the numbers between 10-8 and 1.0, and 1 lambda.1. 12 was chosen from the numbers between 10
was chosen from the
Prediction performance on the other datasets. To complement the results in Sec. 6 for the Dialysis and oldmort datasets, we conducted experiments using the dataDIVATI, flchain, prostateSurvival, and support datasets with K = 2. The outcomes of these experiments are presented in Fig. 6. Unlike the Dialysis and oldmort datasets, no single model consistently outperformed the others across different metrics and datasets.
We also evaluated the prediction performance using the Framingham and PBC datasets with K 3. As baseline methods, we compared our models with those utilized in Jeanselme al. (2023). Specifically, we compared against: DeepHit (Lee et al., 2018), Deep Survival Machines (DSM) (Nagpal et al., 2021), DeSurv (Danks & Yau, 2022), and Neural Fine-Gray (Neu- ralFG) (Jeanselme et al., 2023), which is a neural network model extending the Fine-Gray model (Fine & Gray, 1999). For these models, we used the implementations available at https: / ////tttthb//////// www//ooa/nsseeeeec Ime/ /eeeeeeaieeerrrreeeeeeeee under MIT license, and performed hyperparameter searches based on the guidelines provided in the source code. We compared our models with DeepHit, DSM, DeSurv, and NeuralFG models using the independence copula, and the results are displayed Fig. 7. These results demonstrate that our two-step algorithm is competitive with these baseline models.
Ablation study on hyperparameter B. We conducted an ablation study on the hyperparameter B in our two-step algorithm. This study aimed to evaluate the prediction performance on the cen-log metric using the TS-LGB model with the parameters Wa E {0, 0.5, 1} and B E (4, 16, 32, 64), where wi is the parameter for the primary event of interest. In this study, we also implemented the method proposed in (Carrière, 1995), which is labeled as 'middle'.
Figure 8 shows the results, where the prediction performances are normalized relative to those with w1 = 0.5 for each B. The results indicate that prediction performances varied significantly depending on the choice of the parameter wi on several datasets when B was small. However, these differences diminished for B Z 32. Regarding the method proposed in (Carrière, 1995), while it is valid only if B + oc in theory, somewhat surprisingly, the results demonstrate that it performed comparable to our algorithm even for small B in practice.
Upper and lower bounds estimation. We illustrate survival functions along with bounds for the six datasets with K = 2 in Fig. 9. The six graphs on the left depict average survival functions. In these graphs, we employed the Kaplan-Meier (KM) estimator (1958) and the copula-graphic (CG) estimator (Zheng& & Klein, 1995) using the Frank copula (2) with parameters 0 = 5 and 0 = 5. Recall that the CG estimator is generalization of the KM estimator. The shaded regions indicate the bounds enclosed by the upper and lower limits (Peterson, 1976). The twelve graphs in the center and right display individual survival functions. To estimate these functions, we used the TS-LGB model combined with the independence copula and the Frank copula with parameters = 5 and 0 = 5. The shaded regions represent the bounds enclosed by the upper and lower limits given in Inequality (15).
The figures displaying the average survival functions indicate that uncertainty due to the unknown copula increases as time progresses. The figures displaying the individual survival functions also show similar uncertainty, but the degree of uncertainty varies by individual. In particular, the right figures for the flchain and support2 datasets showed small uncertainty, even without prior knowledge of the copula. We also see that the estimated uncertainty due to the unknown copula can be narrowed by using the estimated survival functions based on the Frank copula.
27

---

Survival Analysis via Density Estimation
detaDNATT
CataltivAT1
dataDNATI
Con
Coc
Con
15:00
15.519
15:008
dataDIVATT
15
15:0
T5 TS
SC
S0
0.7
0.0
0.0
1/5
0.007
5.000
0.001
0.000
0.0
0.5
1.5
200
2.5
Can-ing (munne
D-cstomiiin Jnsert
CH Incniiotttonn (mvert 13
dataDIvAT1
dataDWVtT
cataCivAAT
dataDMAT1
Con
Coc
Doush4
Desght
Desphi
Do0/16
124408
15:08
15:08
15:008
13.019
15.01P
15:00-
13:09
Log
15.0g
15:bric
15
15 15.0
Sc-Nc
50/4cc
6
0%
0.8
10
12
200
25
30
0.0
622
0'6
06
0.
CD
fichain
fcchann
fichain
Coo
Cnc
Deepli
Seepl
Deasllt-
To-c0
T5-100
TS1G2
fichain
T5.00r
15-687
15-05-
TSlls
Tr-dnie
T5
T5 Drer
Sc-eet
6G
1b
0.003 0.655 0.010 adis 6.50 0,05
0.0 65 05
03
0.4
85
0%
07
Jever
fichain
fichain
fichain
fichain
852
DeegH
Drepih
Dee:44
T54664
15: 06-
YGA
7s4r
7647
1545
so
Tolas
TS
Totcg
Tsdeer
Tidree
50
55
6.0
25
7.
75
00
0.5
0.0
3,0
122
25
3.0
02
0.0
0.
C4
Logarti
prosiateSurvieal
prodtateS.rvvval
prostateSuriival
Con
Coc
837
157
057
Creprts
Despl
prostateSurvival
T5.0n
15-001
15-77
"5.1
TS-17
T5
Sc-Net
Sc
Can-lng (euece 3)
D-calonaiin (ever 1)
(M brstoon (evert
prostateSunvivl
grostateSurivell
prostatee.rvvval
prostateSunstval
Cor
Coo
Eec
C00
T5- 15-
15-02
75400
15 nar
T5:07
15-587
T5-r-
75-87
T5ddeer
TSd
Schet
Sc
SCAe
Sc
SL-SC
qp-srit
CD-55
support2
Suppor12
support2
Deeplt
15: oo
support2
Tsar
TEA
15-47
Toins
TS
Toico
TS-nner
T5
0.
0.5
1.0
1.5
0.000
0.005
0.000
0,15
6.000 2.555
0.0
0.2
84
0.6
0.
3.0
Supgort2
support2
support2
support2
R51
357
Despite
Diep111
Dasphi::
De
15-02
15-000
15-t7
T5
75-0
15-17
T5-ing
15
T5-kop
T5-bniu
15
6.0
055
7/2
0.0
1.5
335
0.0
0.6
0.3
op-teer
qptogettmm
CD45
Figure 6. Prediction performance comparison on VatTDIVAT1, flchain, prostateSurvivall, and support2 datascts with various metrics (lower is better).
28

---

Survival Analysis via Density Estimation
Framingham
Framingham
Framingham
Framingham
MeunltS
Meureff
Wesnir
6S9
6S4
Deapeio
Duspr9i
Despet
Dusphti
15-660
75 282
15
75-17-0
T5-00g
Iserien
15
15
0.20 0.23
c.3
2.25
1.50
0.0000 0.0235 0030
0130 0.0125
0:50
0.02
0.03
0.04
6:25
21
lsven)
lig (
(even)
Framingham
Framingham
Framingham
Framingham
Meurat0
Meuni
Meura
De5urv
De5ure
Deppee
Deepp4
15-DRT
15-DR7
T5.087
15
15:ung
15:009
Mer
SC
SC Am
NLUSSC
QJ enerr
CID-55
PBC
FBC
PBC
Neurahn
Neuring
ESM
Defar
Geepeen
Brep-e:
Deeppes
15-DR7
75.087
15-DR7
15.097
15
15-
15:009
Men
SC-Ner
SC
SC
0.77
0.00 ados
Con log revent
D-caltration levenn
Con-log rovert 21
D-colibariin tevent 21
Meurai7
MeurafG
Meaning
Merrifc
015arv
U:5urv
T5Lco
754G8
Te408
T5400
15-687
75588
T5-87
7S.87
15g
15:uy
15.0gg
T5-ugg
CC Ape
SC Me
SC Mae
SC-
NLUSC
QU Stief
CEAKS
Figure 7. Prediction performance comparison on Framingham and PBC datasets with various metries (lower better).
dataDivAT1
Dialysis
Relstive Con-loo
2.0
Rolstive Con (onilo
2.0
w1=0
wl=0
w1=0.5
w1=0.5
1.5
w1=1
3.5
w1=1
middle
middle
1.0
0g
2 2.0%
2088
0003
3.0
ooow
-00%
0004
166
64
26
64
fichain
oldmort
2.0
2.0
8
wlm0
Rolstive Con log
wlm0
Rolstvv Con
w1=0.5
w1=0.5
1.5
w1=1
2.5
wl=1
middle
middle
1.0
2692
2.0
2g
-.
166
64
36
64
prostateSurvival
support2
2.0
2.0
Ressive Can 100
8
w1=0
Reasive Con log
8
w1=0
w1m0.5
w1m0.5
1.5
wl=1
2.5.
wl-1
middle
middle
2.0
432
3044
3604
3554
2.0
16
64
16
64
B
B
Figure Ablation study on hyperparameter B on the six datasets (lower is better).
29

---

Survival Analysis via Density Estimation
detAD1 Average survival function
dataDivAT1: Individual survival Function
dataDIVAT1: Individual survival function
Copura
(theta=55)
Frank (t:sta=a5))
Surwhel preah
Kapaan Melrr
Sunvhal prob.
Indescncence
Sunviol prot.
0.8
Indopandence
0.5-
0.8
0.6
0.4
0.0
0.4
0.2
0.2
0.2
0.9
0.0
0.0
1000
2000
3000
5000
6000
1000
0000
Time
Dialysis: Average survival function
Dialysis: Inciviidul survival function
Dialysis: Individual survival function
Copulo Graphic
2.0
Frank [thota=25) Indesendence
1.0
Fank (theta-15) Indopendonce
Kaplan Mecer
0.8-
0.8
Survical
0.0
30
60
10
20
40
Time
Richain: Average survival function
fichain: Individual survival function
fichain: Individual survival function
1.0
2.0
1.0
Frank (theta=45)
0.5
0.0
0.8
Independence
0.6
0.6
0.4
0.4
0.4
0.2
Copuls Grophic
Frank (chcta=25)
Keplan-Mear
Independence
1000
2009
3000
4000
5000
1000
2000
2000
4000
5000
0
1000
2000
3000
4000
5000
Time
Time
Time
oldmort: Average survival function
oldmort: Individual survival function
aldmort: Individual survival function
Copuls Graphic
3.0
1.0
Frank (thota=25)
0.8
Keplen-Meer
0.0
0.8
Indepaendence
Survoal
0.6
0.6
Sunvvall
0.6
0.4-
0.0
0.4
0.2
0.2
Fronk (th)ta=25)
0.2
Independence
0.3
0.0
0.0
15
10
10
Time
Time
Time
prostateSurvival: Average survival function
prostateSurvival: Individual survival function
prostatesurvival: Individual survival function
0.6
Suoweel
0.4
0.4
Copula-Grappic
0.7
(thelam+3)
Keplan-Mwer
Independence
Indeppendence
0.0
0.0
0.0
:00
80
support2: Average survival function
Support2: Individual survival function
support2: Individual survival function
Copute-Graphic
2.0
Frank (5hetame5)
1.0
Frank (thetam:5)
0.8
Kepler-Reeeer
0.0
Independence
0.8
Independence
0.5
0.6
0.6
Sunkhol
Supvtal
0.0
500
1000
1500
2000
500
1000
1500
2000
500
1000
1500
2000
Time
Time
Time
Figure 9. Estimated survival functions with upper and lower bounds for the six datasets with K = 2. The left graphs show average survival functions, while the graphs in the center and the right show arbitrary chosen individual survival functions. In these graphs, the shaded region corresponds to the upper and lower bounds of survival functions accounting for the uncertainty arising from the lack of knowledge about the copula.
30


Reviews:
[{"review_id": "u91zAyeprz", "reviewer": "Reviewer_gMtj", "summary": "This paper presents a framework that reframes survival analysis as a density estimation problem. By post-processing density estimates to derive survival functions, the approach enables the use of any density estimation model for survival analysis, including handling competing risks and dependent censoring.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "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": "> The experimental evaluation involves only a few datasets, raising concerns about the generalizability of the approach.\n\nAs we state in the last paragraph of Section 6, we include additional evaluation results in the appendix.  Specifically, we present our experimental results (Fig. 6) on four datasets with $K=2$ in Section G, in addition to the experimental results (Fig. 5) on two datasets with $K=2$ in the main body of our paper.  Furthermore, Section G includes experimental results (Fig. 7) on two datasets with $K=3$.   In summary, we used eight datasets in total, as summarized in Table 2  (page 27).\n\n> The proposed method appears to perform poorly with certain classes of density estimators and exhibits large variations across different methods. This issue is not thoroughly discussed in the paper.\n\nThank you for your insightful comments. Compared to standard regression analysis, we need to deal with two types of difficulties in survival analysis: (i) estimating a probability distribution and (ii) handling censored data. The key to our approach is to decouple these two difficulties into Step 1 (which handles (i)) and Step 2 (which handles (ii)). The large variations coming from (i) in our experimental results indicate that survival analysis is difficult even if the censored data are absent (i.e., (ii) is removed) and our future work can be seeking for the best density estimator to improve the performance of survival analysis models.", "claims_and_evidence": "The paper highlights several limitations of existing survival analysis frameworks, including reliance on the conditional independence assumption, issues with identifiability, and the absence of strictly proper scoring rules for cases with K > 2. It provides empirical evidence and theoretical arguments to support these claims.", "methods_and_evaluation_criteria": "The proposed method is primarily a two-step process. In the first step, the authors estimate the cumulative incidence function, and in the second step, they compute the cumulative distribution function using the estimates obtained in the first step along with a specified copula. The paper presents the general case of a sequential approach to solving the proposed function and also provides a simplified implementation of the algorithm. Additionally, the method can be used to estimate upper and lower bounds and is flexible enough to incorporate identifiability results. The proposed approach is evaluated on two datasets and across five different models using a standard survival analysis setup.", "theoretical_claims": "No", "experimental_designs_or_analyses": "The experimental design appears standard; however, the primary concern lies with the choice and number of datasets used.", "relation_to_broader_scientific_literature": "The paper is primarily related to the survival analysis literature.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "Strengths:\n- Model-agnostic framework.\n\nWeaknesses:\n- The experimental evaluation involves only a few datasets, raising concerns about the generalizability of the approach.\n- The proposed method appears to perform poorly with certain classes of density estimators and exhibits large variations across different methods. This issue is not thoroughly discussed in the paper.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "N/A", "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": "", "extra_scores": {}}, {"review_id": "DTU4XkUSiv", "reviewer": "Reviewer_hJPD", "summary": "**Summary:**\nThe authors introduce an algorithm that can post-process density estimators to perform survival analysis. Additionally, a relaxed assumption to the typical conditional independence between event times is used, by modelling the joint distribution using copulas. Several generalisations and scenarios involving more than two competing events are analysed.\n\n## update after rebuttal\nThanks for your response. You've satisfactorily addressed all my comments. It seems I was also the most negative reviewer, and the other reviewers have somewhat swayed my decision. Given these two points, I increase my score.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "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": "Thank you for your comments. We hope our answers resolve your concerns regarding the clarity of our paper.\n\n> Only two datasets are evaluated --- more datasets would make sense for the problem at hand.\n\nAs we state in the last paragraph of Section 6, we include additional evaluation results in the appendix.  Specifically, we present our experimental results (Fig. 6) on four datasets with $K=2$ in Section G, in addition to the experimental results (Fig. 5) on two datasets with $K=2$ in the main body of our paper.  Furthermore, Section G includes experimental results (Fig. 7) on two datasets with $K=3$.   In summary, we used eight datasets in total, as summarized in Table 2  (page 27).\n\n> Preliminaries. Up until the last sentence, it seems like we are dealing with general K. Then in the last sentence, we seem to be dealing exclusively with K=2 - \"however, in this paper, we assume blah takes values from \". Whether the index starts at 0 or 1 is separate to my confusion. Are we dealing with K=2 or general K?\n\nOur paper deals with general K. The sentence is meant to note that the start index is different from most of the existing literature (with K=2). We will rewrite this sentence to avoid confusion.\n\n> I am confused as to how we can interpret the observation (x, t, \\delta) as belonging to a \n-dimensional space. Doesn't x belong to ${\\cal X}$, an arbitrary feature space? Can't this be of any dimension? I can understand how (t, \\delta) can be interpeted in a $K$-dimensional space.\n\nYou are right. We mean that $(t, \\delta)$ can be interpreted in a $K$-dimensional space. We will fix the description in our future revision.\n\n> I am confused about \"This implies the existence of a function $g_C$ such that\". Don't we actually know exactly what this function is, described by (6-13)? That is, not only does it imply an existence, but also that we have this function. I am assuming this is the function that is actually used in algorithm 1? The wording seems to hint that somehow this function exists but we might have to find it, whereas I don't think this is actually the case.It is essentially just the definition of the CJD, CIF and copula.\n\nYou are right.  The function $g_C$ is exactly the combination of (6)-(13).  We will fix the description in our future revision.\n\n> I am confused about equation 15. The symbol $\\hat{r}_{b|x}$  is never defined or explained. What are you minimising (15) with respect to? What is the reason for the gap between algorithm 1 and the \"simplified implementation\"?\n\nWe admit that the definition of $\\hat{r_{b|x}}$ is missing, and we will fix this in our future revision. The definition of $\\hat{r_{b|x}}$ is analogous to the definition of $r_{b|x}$, where each element $r_{b,k|x}$ of $r_{b|x}$ is replaced with $\\hat{r}_{b,k|x}$.\n\nWe minimize (15) with respect to $F_{b|x}$ (and $F_{b-1|x}$) for a given $\\hat{r}_{b|x}$, which is equivalent to solving (14) for all $b$ simultaneously (because the value of (15) is equal to zero if the equality (14) holds for all $b$). Whereas we need to implement a specialized algorithm to solve (14) as we explain in Section B.3, it is much easier to implement an algorithm to minimize (15) because we can use a PyTorch library to minimize (15). We will clarify this fact in our future revision.\n\n> The abstract and intro and several other places talk about survival analysis being performed as a post-processing of density estimation procedures. However, I am not seeing this in algorithm 1, equation 14 or equation 15. Could the authors be more explicit about exactly how this is a postprocessing of density estimation? It is clear that the claim is true, as the authors later show different density estimation techniques in the experiments.\n\nGiven an estimate $r_{b,k|x}$ (possibly together with an estimate in the form of $V_{k}(t|x)$) as an output of Step 1, Step 2 just solves equation (14) (or equivalently minimizes (15)) to obtain $F_{b,k|x}$. In other words, Step 2 is simply a transformation from the representation $r_{b,k|x}$ to another equivalent representation $F_{b,k|x}$. Therefore, we use the term “postprocessing” for Step 2.\n\n> The first sentence of the experiments section talks about evaluating the two-step algorithm. Does this mean you don't use the previously mentioned \"simplified implementation\"? It is not clear which algorithm you are using.\n\nWe used the two-step algorithm, and we used the simplified implementation in Step 2 of this algorithm. We will clarify this fact in our future revision.", "claims_and_evidence": "The claims appear to mostly be supported by evidence, however I have some confusions. See detailed comments below.", "methods_and_evaluation_criteria": "- Only two datasets are evaluated --- more datasets would make sense for the problem at hand.", "theoretical_claims": "- Preliminaries. Up until the last sentence, it seems like we are dealing with general K. Then in the last sentence, we seem to be dealing exclusively with K=2 - \"however, in this paper, we assume blah takes values from ${1,2}$\". Whether the index starts at 0 or 1 is separate to my confusion. Are we dealing with K=2 or general K?\n- I am confused as to how we can interpret the observation (x, t, \\delta) as belonging to a $K$-dimensional space. Doesn't x belong to $\\mathcal{X}$, an arbitrary feature space? Can't this be of any dimension? I can understand how (t, \\delta) can be interpeted in a $K$-dimensional space.\n- I am confused about \"This implies the existence of a function $g_C$ such that\". Don't we actually know exactly what this function is, described by (6-13)? That is, not only does it imply an existence, but also that we have this function. I am assuming this is the function that is actually used in algorithm 1? The wording seems to hint that somehow this function exists but we might have to find it, whereas I don't think this is actually the case.It is essentially just the definition of the CJD, CIF and copula.\n- I am confused about equation 15. The symbol $\\hat{r}_{b \\mid x}$ is never defined or explained. What are you minimising (15) with respect to? What is the reason for the gap between algorithm 1 and the \"simplified implementation\"?\n- The abstract and intro and several other places talk about survival analysis being performed as a post-processing of density estimation procedures. However, I am not seeing this in algorithm 1, equation 14 or equation 15. Could the authors be more explicit about exactly how this is a postprocessing of density estimation? It is clear that the claim is true, as the authors later show different density estimation techniques in the experiments.", "experimental_designs_or_analyses": "- The first sentence of the experiments section talks about evaluating the two-step algorithm. Does this mean you don't use the previously mentioned \"simplified implementation\"? It is not clear which algorithm you are using. \n- Only two datasets are evaluated.", "relation_to_broader_scientific_literature": "I believe the results are related to the broader scientific literature, as discussed in the introduction. In particular, prior works that use density estimation to perform survival analysis are discussed. Previous works have used stronger assumptions.", "essential_references_not_discussed": "None to the best of my knowledge.", "other_strengths_and_weaknesses": "**Detailed Comments/Questions/Weaknesses:**\n- Preliminaries. Up until the last sentence, it seems like we are dealing with general K. Then in the last sentence, we seem to be dealing exclusively with K=2 - \"however, in this paper, we assume blah takes values from ${1,2}$\". Whether the index starts at 0 or 1 is separate to my confusion. Are we dealing with K=2 or general K?\n- I am confused as to how we can interpret the observation (x, t, \\delta) as belonging to a $K$-dimensional space. Doesn't x belong to $\\mathcal{X}$, an arbitrary feature space? Can't this be of any dimension? I can understand how (t, \\delta) can be interpeted in a $K$-dimensional space.\n- I am confused about \"This implies the existence of a function $g_C$ such that\". Don't we actually know exactly what this function is, described by (6-13)? That is, not only does it imply an existence, but also that we have this function. I am assuming this is the function that is actually used in algorithm 1? The wording seems to hint that somehow this function exists but we might have to find it, whereas I don't think this is actually the case.It is essentially just the definition of the CJD, CIF and copula.\n- I am confused about equation 15. The symbol $\\hat{r}_{b \\mid x}$ is never defined or explained. What are you minimising (15) with respect to? What is the reason for the gap between algorithm 1 and the \"simplified implementation\"?\n- The abstract and intro and several other places talk about survival analysis being performed as a post-processing of density estimation procedures. However, I am not seeing this in algorithm 1, equation 14 or equation 15. Could the authors be more explicit about exactly how this is a postprocessing of density estimation? It is clear that the claim is true, as the authors later show different density estimation techniques in the experiments.", "other_comments_or_suggestions": "See above.", "questions_for_authors": "See above.", "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": "", "extra_scores": {}}, {"review_id": "HKvvqDJijG", "reviewer": "Reviewer_qovr", "summary": "The authors consider the problem of survival analysis with competing and potentially dependent risks. The paper has two main contributions. First, the authors propose a two-step plug-and-play method which uses the output of a generic density estimator and transforms it into an estimate of the joint survival function. Their method relies on knowledge of a copula describing the dependence structure between the competing risks. When the independence copula is used, this reduces to the standard assumption of conditional independence among the competing risks. When the copula is unknown, they also provide upper and lower bounds on the resulting survival function. Second, the authors introduce a proper scoring rule for competing risks which is valid for any number of competing risks. They use this proper scoring rule both as an evaluation metric, and as a training objective for a neural network-based model of the joint survival function. Finally, they compare several instantiations of their two-step method, their neural network approach, and several standard survival analysis baselines on two real datasets. They obtain consistent improvements in terms of both accuracy and calibration compared to the baselines.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "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": "Thank you for your comments. We appreciate your remark, \"The two-step method proposed by the paper is flexible and innovative,\" and we are grateful for the suggestions to improve the presentation of our paper.\n\n> The suite of evaluation metrics are extensive and well-explained in the paper. The datasets used are relevant, though somewhat limited (there are only two).\n\n> The only obvious shortcoming of the paper is the lack of baseline methods in the experiments.\n\nAs we state in the last paragraph of Section 6, we include additional evaluation results in the appendix.  Specifically, we present our experimental results (Fig. 6) on four datasets with $K=2$ in Section G, in addition to the experimental results (Fig. 5) on two datasets with $K=2$ in the main body of our paper.  Furthermore, Section G includes experimental results (Fig. 7) on two datasets with $K=3$.   In summary, we used eight datasets in total, as summarized in Table 2  (page 27).\n\nWe would also like to note that we used different sets of baseline methods for $K=2$ and $K=3$, and we used six baseline methods (i.e., not only three) in total.\n\n> I checked the proof of Theorem 5.3 in Appendix F and did not find any errors. A minor comment: I assume $v$ and $\\hat{v}$ are the distributions over $(t,k)$  induced by each of the $v_k$ and $\\hat{v}_k$. This may be apparent but I did not see it explicitly stated anywhere.\n\nWe will update the description based on your comments if our paper is accepted.\n\n> Please also provide some references for prior works which used bespoke extensions of specific density estimators to survival analysis (these are mentioned in the related work without references).\n\nWe would like to remind you that we show some examples of the bespoke extensions with references in the second paragraph of Section 1.\n\n> Is there an intuitive explanation for condition (iii) in the definition of a copula (pg. 3, just above equation (1)?\n\nWhen $K=2$, this condition is $C(v_1,v_2) - C(u_1,v_2) - C(v_1,u_2) + C(u_1,u_2) \\geq 0$. This is a necessary condition to represent a probability $\\Pr(\\zeta_1 < T_1 \\leq \\zeta_2, \\zeta_3 < T_2 \\leq \\zeta_4)$ as $\\Pr(\\zeta_1 < T_1 \\leq \\zeta_2, \\zeta_3 < T_2 \\leq \\zeta_4) = C(v_1,v_2) - C(u_1,v_2) - C(v_1,u_2) + C(u_1,u_2)$, where $v_1 = F_1(\\zeta_2)$, $v_2 = F_2(\\zeta_4)$, $u_1 = F_1(\\zeta_1)$, and $u_2 = F_2(\\zeta_3)$, because the probability must be non-negative.\n\n> Can the authors provide the derivation for the equation of the derivative $v_k(t|x)$ of the cumulative incidence function (pg. 7, just above Assumption 5.2)?\n\nThe proof is immediate from Theorem 1 of (Tsiatis, 1975). Note that $H^{(k)}(t_1,t_2,...,t_k)$ in (Tsiatis, 1975) corresponds to $\\overline{C}(1-F_{1}(t_{1}),1-F_{2}(t_{2}),\\ldots,F_{K}(t_{K}))$ in our paper and $Q_i(t)$ in (Tsiatis, 1975) corresponds to $V_k(\\infty) - V_k(t)$ in our paper. We will cite (Tsiatis, 1975) in our future revision.\n\n> How are the class probabilities $f_k(x^{(i)})$ which are used to compute the KS calibration error (bottom of pg. 7) recovered from the survival function estimates?\n\nThey can be computed by using $r_{k}(t|x)$, and we can compute each $r_{k}(t|x)$ by using equation (14) (equivalently, equations (6)-(13)) since we have an estimate $F_{k}(t|x)$.  We will add explanation on the computation in our future revision.", "claims_and_evidence": "The only obvious shortcoming of the paper is the lack of baseline methods in the experiments. The authors compare against a Cox model, a random survival forest, and a neural network-based method (DeepHit). However, there are many other methods from the literature which are applicable to this problem. It would greatly strengthen the paper to compare against these as well. Some examples include:\n\nKatzman, Jared L., et al. \"DeepSurv: personalized treatment recommender system using a Cox proportional hazards deep neural network.\" BMC medical research methodology 18 (2018): 1-12.\n\nKvamme, Håvard, Ørnulf Borgan, and Ida Scheel. \"Time-to-event prediction with neural networks and Cox regression.\" Journal of machine learning research 20.129 (2019): 1-30.\n\nHu, Bingqing, and Bin Nan. \"Conditional distribution function estimation using neural networks for censored and uncensored data.\" Journal of Machine Learning Research 24.223 (2023): 1-26.\n\nIt seems it would also be reasonable to include results for some of the methods which are already discussed in the related work, in particular the monotone neural network model of Rindt et al. (2022).", "methods_and_evaluation_criteria": "The two-step method proposed by the paper is flexible and innovative. The ability to turn any density estimator into a method for survival analysis with competing risks greatly expands the toolkit available to practitioners, and the experimental results so far indicate that it obtains consistent improvement over existing methods. While the assumption of a known copula may seem restrictive at first, the authors clearly explain that the more common assumption of conditional independence of the competing risks is actually a special case of this copula assumption, so it is never a more severe modeling restriction than the independence assumption. The derived upper and lower bounds on the survival function even go a step further in alleviating this modeling assumption.\n\nWhile the introduced strictly proper scoring rule is helpful as an evaluation metric in the setting of this paper, especially in the $K>2$ setting. However, the introduction of neural network model based on this score feels a bit out of place in the paper, as relatively little space is devoted to it and its performance is not up to par with the other methods.\n\nThe suite of evaluation metrics are extensive and well-explained in the paper. The datasets used are relevant, though somewhat limited (there are only two).", "theoretical_claims": "I checked the proof of Theorem 5.3 in Appendix F and did not find any errors. A minor comment: I assume $v$ and $\\hat{v}$ are the distributions over $(t, k)$ induced by each of the $v_k$ and $\\hat{v}_k$. This may be apparent but I did not see it explicitly stated anywhere.", "experimental_designs_or_analyses": "The experimental design is sound. As mentioned in the Claims and Evidence section, while the structure of the experiments is sound, it is missing some relevant materials and could benefit from experiments on more datasets.", "relation_to_broader_scientific_literature": "The authors give a very clear description of the relationship between their work and the prior literature. Specifically:\n- Previous extensions of density estimation to survival analysis are bespoke for that particular density estimator, whereas their method works for any density estimation method. (A remark on this point: it would be helpful for the authors to provide some references for the prior works they're talking about.)\n- Previous works rely on stronger assumptions on the competing risks, such as conditional independence or proportional hazards.\n- Previous work has shown that the survival function is not identifiable without additional assumptions, but it has not examined whether or not the survival function can be *bounded* even if it is not exactly identifiable.\n- Previous works established strictly proper scoring rules in more restricted settings ($K=2$ competing risks and conditional independence), but a strictly proper scoring rule has not been established in the more general setting of this paper.", "essential_references_not_discussed": "See the Claims and Evidence section for a collection of baselines which should be added.  Please also provide some references for prior works which used bespoke extensions of specific density estimators to survival analysis (these are mentioned in the related work without references).", "other_strengths_and_weaknesses": "This paper is exceptionally well-written. I was especially impressed by the color-coded visual explanation of the core formula for their two-step method, represented in equations (6)-(10) and Figure 3.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "- Is there an intuitive explanation for condition (iii) in the definition of a copula (pg. 3, just above equation (1)?\n- Can the authors provide the derivation for the equation of the derivative $v_k(t|x)$ of the cumulative incidence function (pg. 7, just above Assumption 5.2)?\n- How are the class probabilities $f_k(x^{(i)})$ which are used to compute the KS calibration error (bottom of pg. 7) recovered from the survival function estimates?\n---\nAfter the rebuttal phase, the authors have addressed my main concerns. I looked at the other reviews and author responses, and it seems that these concerns have been addressed as well. I have raised my score accordingly.", "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": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: After the discussion period, the reviewers unanimously favored accepting the paper (including two reviewers increasing their scores, with one going from leaning on rejection to leaning on acceptance). Looking over the discussion, I find that the authors did a great job addressing reviewer concerns with detailed explanations. I am recommending this paper to be accepted.
Decision: Accept (poster)