Neural Frailty Machine: Beyond proportional hazard assumption in neural survival regressions
Authors: Ruofan Wu
Jiawei Qiao
Mingzhe Wu
Wen Yu
Ming Zheng
Tengfei LIU
Tianyi Zhang
Weiqiang Wang
Venue: NeurIPS 2023 poster (2023)
Keywords: Survival Analysis
Theory
Semiparametric statistics

Abstract: We present neural frailty machine (NFM), a powerful and flexible neural modeling framework for survival regressions. The NFM framework utilizes the classical idea of multiplicative frailty in survival analysis as a principled way of extending the proportional hazard assumption, at the same time being able to leverage the strong approximation power of neural architectures for handling nonlinear covariate dependence. Two concrete models are derived under the framework that extends neural proportional hazard models and nonparametric hazard regression models. Both models allow efficient training under the likelihood objective. Theoretically, for both proposed models, we establish statistical guarantees of neural function approximation with respect to nonparametric components via characterizing their rate of convergence. Empirically, we provide synthetic experiments that verify our theoretical statements. We also conduct experimental evaluations over $6$ benchmark datasets of different scales, showing that the proposed NFM models achieve predictive performance comparable to or sometimes surpassing state-of-the-art survival models. Our code is publicly availabel at https://github.com/Rorschach1989/nfm

Neural Frailty Machine: Beyond proportional hazard assumption in neural survival regressions
Ruofan Wu*t, Jiawei Qiao" Mingzhe Wu
Wen Yu', Ming Zheng', Tengfei Liui, Tianyi Zhang?, and Weiqiang Wang
Ant Group
Fudan University
Coupang
{ruofan. wuff aaron. 1tf, zty113091, veiqiang.wwq/ wwt.haatggru..pccoooomo com j jeremyqjw0163.c com, wumingzhe dawcyygmail..ccomoomoommomom.... com, (wenyu,
edu..
Abstract
We present neural frailty machine (NFM), a powerful and flexible neural modeling framework for survival regressions. The NFM framework utilizes the classical idea of multiplicative frailty in survival analysis as a principled way of extend- ing the proportional hazard assumption, the same time being able to leverage the strong approximation power of neural architectures for handling nonlinear covariate dependence. Two concrete models are derived under the framework that extends neural proportional hazard models and nonparametric hazard regres- sion models. Both models allow efficient training under the likelihood objective. Theoretically, for both proposed models, we establish statistical guarantees of neural function approximation with respect to nonparametric components via characterizing their rate of convergence. Empirically, we provide synthetic ex- periments that verify our theoretical statements. We also conduct experimental evaluations over 6 benchmark datasets of different scales, showing that the pro- posed NFM models achieve predictive performance comparable or sometimes surpassing state-of-the-ee art survival models. Our code is publicly availabel at stttps://gite www./oom////eeeeee Nomm/oorscacch198//ee
Introduction
Regression analysis of time-to-event data 40 has been among the most important modeling tools for clinical studies and has witnessed growing interest in areas like corporate finance [22]. rec- ommendation systems (38), and computational advertising [71]. The key feature that differentiates time-to-event data from other types of data is that they are often incompletely observed, with the most prevailing form incompleteness being the right censoring mechanism 40) In the right censoring mechanism, the duration time of sampled subject is (sometimes) only known to be larger than the observation time instead of being recorded precisely. is well known the community of survival analysis that even in the case of linear regression, naively discarding the censored observations produces estimation results that are statistically biased at the same time losses sample efficiency if the censoring proportion is high.
using the convex objective of negative partial likeli- hood [15) is the de facto choice in modeling right censored time-to-event data (hereafter abbreviated as censored data without misunderstandings). The model semiparametrie 4 in the sense that the baseline hazard function needs no parametric assumptions. The original formulation of CoxPH model
Cox's proportional hazard (CoxPH) model [4
'Equal contribution
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

assumes a linear form and therefore has limited flexibility since the truth is not necessarily linear. Subsequent studies extended CoxPH model to nonlinear variants using ideas from nonparametric
regression 36 g 9
ensemble learning [37] and neural networks [25 42 While such extensions allowed a more flexible nonlinear dependence structure with the covariates, the learning objectives were still derived under the proportional hazards (PH) assumption, which was shown to be inadequate in many real-world scenarios [31]. The most notable case was the failure of modeling the phenomenon of crossing hazards [61. It is thus of significant interest to explore extensions of CoxPH that both allow nonlinear dependence over covariates and relaxations of the PH assumption.
Frailty models [69] 21 are among the most important research topics in modern survival analysis, in that they provide a principled way of extending CoxPH model via incorporating a multiplicative random effect to capture unobserved heterogeneity. The resulting parameterization contains many useful variants of CoxPH like the proportional odds model [3], under specific choices of frailty families. While the theory of frailty models has been well-established [48
49 53 45], most of them focused on the linear case. Recent developments on applying neural approaches to survival analysis 42 46 64 56 have shown promising results in terms of empirical predictive performance, with most of them lacking theoretical discussions. Therefore, it is of significant interest to build more powerful frailty models via adopting techniques in modern deep learning [29] with provable statistical guarantees.
In this paper, we present general framework for neural extensions of frailty models called the neural frailty machine (NFM). Two concrete neural architectures are derived under the framework: The first one adopts the proportional frailty assumption, allowing intuitive interpretation of the neural CoxPH model with a multiplicative random effect. The second one further relaxes the proportional frailty assumption and could be viewed as an extension of nonparametric hazard regression (NHR)
13
44), sometimes referred to as "fully neural" models under the context neural survival analysis We summarize our contributions as follows.
$2
We propose the neural frailty machine (NFM) framework as a principled way of incorporating unobserved heterogeneity into neural survival regression models. The framework includes many commonly used survival regression models as special cases.
We derive two model architectures based on the NFM framework that extend neural CoxPH models and neural NHR models. Both models allow stochastic training and scale large datasets.
Theoretically, we show statistical correctness of the two proposed models via characterizing the rates of convergence of the proposed nonparametric function estimators. The proof technique is different from previous theoretical studies on neural survival analysis and is applicable to many other types of neural survival models.
Empirically, we verify the empirical efficacy of the proposed framework via conducting extensive studies on various benchmark datasets at different scales. Under standard performance metrics, both models are empirically shown to perform competitively, matching or sometimes outperforming tate-offtheeart neural survival models.
2
Related works
2.1 Nonlinear extensions of CoxPH
Most nonlinear extensions of CoxPH model stem from the equivalence of partial likelihood and semiparametric profile likelihood (50) of CoxPH model, resulting in nonlinear variants that essentially replaces the linear term in partial likelihood with nonlinear variants: [36] used smoothing splines, [8]9] used local polynomial regression [24). The empirical success of tree-based models inspired subsequent developments like (37] that equip tree-based models such as gradient boosting trees and random forests with losses in the form of negative log partial likelihood. Early developments of neural survival analysis (25] adopted similar extension strategies and obtained neural versions of partial likelihood. Later attempts 42) suggested using the successful practice of stochastic training which is believed to be at the heart of the empirical success of modern neural methods (4) However, stochastic training under the partial likelihood objective is highly non-trivial, as mini-batch versions of log partial likelihood 42 are no longer valid stochastic gradients of the full-sample log partial likelihood [64].

---

2.2 Beyond CoxPH in survival analysis
In linear survival modeling, there are standard alternatives CoxPH such as the accelerated failure time (AFT) model [ [73]. the extended hazard regression model [23], and the family of linear transformation models (74]. While these models allow certain types of nonlinear extensions, the resulting form of (conditional) hazard function is still restricted to be of specific form. The idea of nonparametric hazard regression (NHR) 13
[X] further improves the flexibility of nonparametric survival analysis via directly modeling the conditional hazard function by nonparametric regression techniques such as spline approximation. Neural versions of NHR have been developed lately such as the Cox CoxTime model [46]. [56] used neural network to approximate the conditional survival function and could be thus viewed as another trivial extension of NHR.
Aside from developments in NHR, 47 proposed a discrete- time model with its objective being a mix of the discrete likelihood and rank-based score; [75] proposed neural version of the extended hazard model, unifying both neural CoxPH and neural AFT model; [64] used an ODE approach to model the hazard and cumulative hazard functions.
2.3 Theoretical justification of neural survival models
Despite the abundance of neural survival models, assessment of their theoretical properties remains nascent. In (76), the authors developed minimax theories of partially linear model using neural networks as the functional approximator. [75] provided convergence guarantees of neural estimates under the extended hazard model. The theoretical developments therein rely specific forms of objective function (partial likelihood and kernel pseudo-likelihood) and are not directly applicable to the standard ikelihhoddbbsee objective which is frequently used in survival analysis.
3 Methodology
3.1 The neural frailty machine framework
Let T 2 0 be the interested event time with survival function denoted by S(t) P(T > t) associated with feature(covariate) vector Z E IRd. Suppose that T is continuous random variable and let f(t) be its density function. Then ((4) = (L)) is the hazard function and A()) fo X(s)ds is the cumulative hazard function. Aside from the covariate Z. we use a positive scalar random variable
E IR to express the unobserved heterogeneity corresponding to individuals, or frailty. In this
paper we will assume the following generating scheme of via specifying its conditional hazard function:
X(t Z,w) wāit, Z).
(1)
Here i/ is an unspecified non-negative function, and we let the distribution of w be parameterized by a one-dimensional parameter 8 IR. [ The formulation (T) quite general and contains several important models in both traditional and neural survival analysis:
1. When w follows parametric distributional assumptions, and (t. 2) = X(t)e z [ reduces to the standard proportional frailty model 45). special case when w is degenerate, i.e., has no randomness, then the model corresponds to the classic CoxPH model.
When w is degenerate and T is arbitrary, the model becomes equivalent to nonparametric hazard regression (NHR) [13 44). In NHR, the function parameter of interest is usually the logarithm of the (conditional) hazard function.
In this paper we construct neural approximations to the logarithm of i. i.e., v(t, Z) = log itt, Z). The resulting models are called Neural Frailty Machines (NFM). Depending on the prior knowledge of the function we propose two function approximation schemes:
FFor example in medical biology, it was observed that genetically identical animals kept in as similar an environment as possible will typically not behave the same upon exposure to environmental carcinogens 6 'The choice of one- doeddmmnniioal frailty family is mostly for simplicity and clearness of theoretical derivations. Note that there exist multi sittii-dmmnsionall frailty families like the PVF family (69). Generalizing our theoretical results to such kinds of families would require additional sets of regularity conditions, and will be left to future explorations.
3

---

The proportional frailty (PF) scheme assumes the dependence of v on event time and covariates to be completely decoupled, i.e.
- - (Z).
(2)
Proportional-style assumption over hazard functions has been shown to be a useful inductive bias in survival analysis. We will treat both h and in (2) as function parameters, and device two multi-layer
perceptrons (MLP) to approximate them separately.
The fully neural (FN) scheme imposes no a priori assumptions over y and is the most general version of NFM. It straightforward to see that the most commonly used survival models, such as CoxPH, AFT 773, EH[75), or PF models are included in the proposed model space as special cases. We treat V/ rtt, Z) as the function parameter with input dimension + and use multi-layer perceptron (MLP) as the function approximator to v. Similar approximation schemes with respect to the hazard function have been proposed in some recent works [52] [6], referred as "fully neural approaches" without theoretical characterizations.
The choice of frailty family There are many commonly used families of frailty distributions 45 21
69). among which the most popular one is the gamma frailty, where w follows a gamma distribution with mean and variance 0. We briefly introduce some other types of frailty families in appendix A
3.2 Parameter learning under censored observations
In time-to-event modeling scenarios, the event times are typically observed under right censoring. Let C' be the right censoring time which is assumed to be conditionally independent of the event time T given Z. i.e., T C In data collection, one can observe the minimum of the survival time and the censoring time, that is, observe T T as well as the censoring indicator å I(T < C), where A b = min(a, b) for constants a and b and (()) stands for the indicator function. We assume n independent and identically distributed (i.i.d)) copies of (T, 6, Z) are used as the training sample ...2,................. .... Additionally, we assume the (Ti, 8,, , E [n]]. where we use [n] to denote the set .1.
unobserved frailties are independent and identically distributed, i.e., Wi ~ fo(w), E [n]. Next, we derive the learning procedure based the observed log-likelihood (OLL) objective under both PF and FN scheme. To obtain the observed likelihood, we first integrate the conditional survival function given the frailty:
sid.
2( = Eu (o c 2i
(e- = GG 21d)))
Here the frailty transform Gz(z) = log (Eew [e ()) is defined as the negative of the logarithm of the Laplace transform of the frailty distribution. The conditional cumulative hazard function is thus (4)2) = Go(So (ds) For the PF scheme NFM, we use two MLPs π = Ti(1; w* and m. = m Z: W bmd as function approximators to and m, parameterized by (Wh ba) and
(W"m b²), respectively. F According to standard results on censored data likelihood 40 we write
the learning objective under the PF scheme as:
L(W", bd, W b²t,0)
Th 8 log go mizz) f" Mollds (.T((())) + d,m((z) Go (L(z)) foh gioodss
(4)
VE[M]
Here we define go(x) a Go(x). Let (Wh b Wm m b. 0m)) be the maximizer of (4) and n further denote (t) = (1( Wh. no bb) and m, m,) = m Wm. n. D The resulting estimators for conditional cumulative hazard and survival functions are:
ApF(t|Z) = Gā. ( (fennoottmnnl
SpF(t) = -
(5)
'Here we adopt the conventional notation that W is the collection of the weight matrices of the MLP in all layers, and b corresponds to the collection of the bias vectors in all layers.

---

For the FN scheme, we use 1 = itt, Z: W b²) to approximate (t, 2) parameterized by (W", b"). The OLL objective is written as:
L(W",b" .))
& log (unnwwwaa & w",b") Go LEIN
(6)
Let (W Dy 0n) be the maximizer of (6. and further denote 2 (t, 2) itt Zi W bg). The conditional cumulative hazard and survival functions are therefore estimated as:
(T. (eulo. (
AFN(4IZ= GGö
SEN (L =
The evaluation of objectives like ( and its gradient requires computing definite integral of an exponentially transformed MLP function. Instead of using exact computations that are available for only restricted type of activation functions and network structures, use numerical integration for such kinds of evaluations, using the method of Clenshaw- competitive performance and efficiency in recent applications to monotonic neural networks (68]
quadrature (5]. which has shown
Remark 3.1. The interpretation of frailty terms differs in the two schemes. the PF scheme, introducing the frailty effect strictly increases the modeling capability (i.e., the capability of modeling crossing hazard) in comparison to CoxPH or neural variants CoxPH [45). In the FN scheme, is arguable that in the i.i.d. case, the marginal hazard function is a reparameterization of the hazard function the context of NHR. Therefore, we view the incorporation of frailty effect as injecting domain- -spacifii inductive bias that has proven to be useful in survival analysis and time-to-even regression modeling and verify this claim empirically in section 5.2, Moreover, frailty becomes especially helpful when handling correlated clustered data where the frailty term is assumed to be shared among certain groups of individuals [53] Extending NFM to such scenarios is valuable and we left it to future explorations.
4 Statistical guarantees
In this section, we present statistical guarantees regarding both NFM estimates in the sense of nonparametric regression [65), where we obtain rates of convergence the ground truth function parameters (which is frequently referred to as the true parameter in statistics literature). The results in this section is interpreted as showing the statistical correctness of our approach.
Proof strategy Our proof technique is based on the method of sieves [60] 59 that views neural networks as special kind of nonlinear sieve [1] that satisfies desirable approximation properties (72). Our strategy is different from previous thcoretical works neural survival models [75 76 where the developments implicitly requires the loss function to be well-controlled by the Le loss and therefore not directly applicable to model due to the flexibility in choosing the frailty transform. Since both models produce estimates of function parameters, we need to specify a suitable function space work with. Here we choose the following Hölder ball was also used in previous works nonparametric estimation using neural networks [57 26 76
W Mr (X)
{5=
max esssup D* = M
(8)
a:|a|g =EX
where the domain X is assumed to be subset of d-dimensional euclidean space. (a) a d-dimensional tuple of nonnegative integers satisfying |a| = 01 +... 0s and D Bo ...z
Ota) is
al
the weak derivative of f- Now assume that M is reasonably large constant, and let t be a closed interval over the real line. We make the following assumptions for the true parameters under both schemes:
Condition 4.1 (True parameter, PF scheme). The euclidean parameter Og E O C TR and the two function parameters mo W (I-1, 1]d), ho ( r)), and t > the ending time of the study duration, which is usually adopted in the theoretical studies in survival analysis 167).
Condition 4.2 (True parameter, FN scheme). The euclidean parameter Og θ C R. and the function parameter Po W ((0, 7] [11,
5

---

Next, we construct sieve spaces for function parameter approximation via restricting the complexity of the MLPs "scale" with the sample size n.
Condition 4.3 (Sieve space, PF scheme). The sieve space Hn is constructed as a set of MLPs satisfyyng h W [O,, with depth of order (logn) and total number of parameters of order On aLa n). The sieve space Mn is constructed as set of MLPs satisfying in - (1-1,
with depth of order (lognn and total number of parameters of order O s=aa lognt). Here Mn and Mn are sufficiently large constants such that every. furcctinn in W (1-1, 1]d) and W ([0, 7]) could be accurately approximated by functions inside Hr and MM, according to [72 Theorem 1].
Condition 4.4 (Sieve space, FN scheme). The sieve space Vn is constructed as a set of MLPs yatisfying E WM, (10, r]), with depth of order Ollo n) and total number of parameters of order (n #+0+I logn). Here M, is sufficiently large constant such that Vn satisfies approximation properties, analogous condition 4.3
For technical reasons, we will assume the nonparametric function estimators are constrained to fall inside the corresponding sieve spaces, i.c., fi E Hn. Mn E Mn and Vn This will not affect the implementation of optimization routines as was discussed in [26]. Furthermore, restrict the estimate 8n E in both PF and FN schemes.
Additionally, we need the following regularity condition on the function Gz(z):
Condition 4.5. Gz( is viewed as bivariate function G : θ x B + R. where B is compact set on R. The functions Go(z). 30 Gz(z), de (o 06 + x B.
& log (z), log ge(z) are bounded on
We define two metrics that measures convergence of parameter estimates: For the PF scheme, let Φ0 = (ho, mo, 0o) be the true parameters and の. (hn, min 0n) be the estimates. We abbreviate F 2= as the conditional probability distribution of (T, 5) given Z = under the true parameter, and P- == -= the conditional probability distribution of (T,5) given Z = 2 under the estimates.
Define the following metric
dpF (fn..h) BEnPZ (2 @n, l Pro..
where (TP Q) = [VvP--
vTT is the squared Hellinger distance between probability distributions and Q. The case for the FN scheme similar: Let s/n (Vo. 0o) be the parameters and i2 (Dz, 0n) be the estimates. Analogous to the definitions above, we define Pwo.Z=z as the
2
true conditional distribution given = 2, and Po.,, be the estimated conditional distribution, we will use the following metric the FN scheme:
dFN (wn, wh)
V BEMPZ (2(Pp || P (oo.
(10)
Now we state our main theorems. We denote IP as the data generating distribution and use to hide poga-togaii factors in the big-O notation.
Theorem 4.6 (Rate of convergence, PF scheme). In the PF scheme, under condition 4.1
4.5 we
have that dpF Pa, Φ6 = Öz n 20122
Theorem 4.7 (Rate of convergence, FN scheme), In the FN scheme, under condition we have that dFN
4.2
4.4
4.5
Wn: 10) = 0p (n #742842
Remark 4.8. The idea of using Hellinger distance measure the convergence rate of sieve MLEs was proposed in (70.. Obtaining rates under stronger topology such as Lg is possible if the likelihood function satisfies certain conditions such as the curvature condition [26]. However, such kind of conditions is in general too stringent for likelihood-base objectives, instead, we use Hellinger distance that has minimal requirements. Consequently, our proof strategy is applicable to many other survival models that rely on neural function approximation such [56), with some modification to the regularity conditions. For proper choices of metrics in sieve theory, see also the discussion in Chapter 2].

---

1000
5000 (P9)
1000
groand trui
groaned trah
graand manh
50
100
100
50
100
1000 JPN)
5000
10000 (PN)
groand moams
100
100
50
100
Figure 1: Visualizations of synthetic data results under the NFM framework. The plots in the first row compare the empirical estimates of the nonparametrie component (tt, Z) against its true value evaluated on 100 hold-out points, under the PF scheme. The plots in the second row are obtained using the FN scheme, with analogous semantics to the first row.
Experiments
In this section, we report the empirical performance of NFM, we will focus on the following two research questions:
of statistical correctness): The results listed in section characterized the conver- gence results in theory, providing crude guide on the number of samples required for an accurate estimate. Nonetheless, theoretical rates often pessimistie, thus want to investigate whether moderate number of sample size suffices for good approximation.
RQ2(Assessment of empirical efficacy): While NFM is theoretically sound in terms of estimation accuracy, the theory we have developed does not necessarily guarantee its empirical efficacy as method of doing prognosis. It is therefore valuable to inspect how useful NFM is regarding real-world predictive tasks in survival analysis.
5.1 Synthetic experiments
To answer RQ1, conduct synthetic experiments check the empirical convergence. Specifically, we investigate the empirical recovery of underlying ground truth parameters under various level of sample size.
Ground truth We set the true underlying model to be nonlinear gamma-frailty model with 5-dimensional feature. We generate three training datasets of different scales, with n (1000, 5000, 10000). The assessment will be made on fixed test sample of 100 hold-out points that are independently drawn from the generating scheme of the event time. A censoring mechanism is applied such that the censoring ratio is around 40% for each dataset. The precise form of the frailty model well the generating distribution of the feature vectors are detailed in appendix C.2 Empirical recovery results We report the empirical recovery of the nonlinear component (t,, Z) regarding the hold-out test set in figure [ We observe from graphical illustrations that under a moderate sample size n 1000, NFM already exhibits satisfactory recovery for (relatively low-dimensional) feature space, which is the prevailing case in most public benchmark datasets. We also present additional assessments about: (i) The recovery of m(() when using PF scheme in appendix [..] (ii) The recovery of survival functions under both and FN scheme in appendix D.2 and (ii) The numerical recovery results of survival function in appendix D.3

---

Table 1: Survival prediction results measured in IBS and INBLL metric (%) on four small survival datasets. In each column, the boldfaced score denotes the best result and the underlined score represents the second-best result (both in mean).
Model
METABRIC
RoIGBSG
FLCHAIN
SUPPORT
IBS
INBLL
IBS
INBLL
IBS
INBLL
IBS
INBLL
CoxPH
16.4620.90
49.57 66
18. 22.5220.4 34
63.76g1.11
10.05
33.18-1.16
20.54.0.35
59.58 s6
GBM
16.61
49.87=2.44
17.8 17.83+0.44
52.78
9.98_...7
32. 32.88
19.18
56, 56.46
RSF
16.6220.64
49.61
17.89
52.77
9.96
32.92=1.05
19.11
56.28
DeepSurv
16.55go.93
49.85.3.00 oz
17.802o.a aM
62.62
10.09
33.28
19.20
56.48 21 os
CoxTime
16.51g0..8
49.67=2.67
17 17.80
52.56
10.28
34.18:11.55
19.17
56, 56,45
DeepHiit
17.50
52. 22.1002..1
19.61 +0.38
56.67
11 11.83
37.72=1100
20.66 20.66:0.332
60.06
DeepEH
16.56gu..6
49.42 53
17.62.0.52
62.08
10.11
33.30_1.10
19.30
56.67
SuMo-net
16.49:0.83
49.74g2.22 21
17.77+0047
62.62+1.11
10.07
33.20=:1.10
19.40
56.87
SODEN
16.52:0.6a
49.39=+1.9
17.05 17.05
50.45=1.97
10.13:20.24
33.37
19.07 #0.50
56.15g1.35
SurvNode
16.67
49.73ga3. s9
17.42.0.55
11.70
10.40
34.37=1.03
19.58
57.49
DCM
16.58
49.48gg2.0 23
17.66.0.50
52:26::125
10.13
33.40-.1.88
19.29
56, 56.68%21
DeSurv
16.71
49.61
17.98
13.23
10.0640.02
33. 33.18=1.0
19.50
57.28
NEM-PF
16.73
19.07.
40.00
52.12
9.96
32 84
19.1420.30
56.35g1
NFM-FN
16.11 16.11
48. +2.04
17.66
52.41+1.22
10.05
33. 33.11g1.10
18.97
55.87
5.2 Real-world data experiments
To answer RQ2, we conduct extensive empirical assessments over 6 benchmark datasets, comprising five survival datasets and one non-survival dataset. The survival datasets include the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) [16), the Rotterdam tumor bank and German Breast Cancer Study Group (RoIGBSG)33 the Assay Of Serum Free Light Chain (FLCHAIN) 20), the Study to Understand Prognoses Preferences Outcomes and Risks of Treatment (SUPPORT) 43 and the Medical Information Mart for Intensive Care (MIMIC-III [39]. For all the survival datasets, the event of interest is defined as the mortality after admission. In our experiments, we view METABRIC, RotGBSG, FI.CHAIN, and SUPPORT as small-scale datasets and MIMIC-III as a moderate-scall dataset. We additionally use the KKBOX dataset [46) as a large-scale evaluation. In this dataset, an event time is observed if customer churns from the KKBOX platform. We summarize the basic statistics of all the datasets in table
Baselines We compare NFM with 12 baselines. The first one is linear CoxPH model [4]. Gradient Boosting Machine (GBM)
2 10 and Random Survival Forests (RSF) 37 are two tree-based
nonparametric survival regression methods. DeepSurv 42 and CoxTime 46 are two models that adopt neural variants of partial likelihood as objectives. SuMo-net [56] is neural variant of NHR. We additionally chose six latest state-of theth-rrt neural survival models: DeepHit [47], SurvNode (2), DeepEH 751, DCM [11. DeSurv [9] and SODEN 64]. Among the chosen baselines, DeepSury and SuMo-net are viewed as implementations of neural CoxPH and neural NHR and are therefore of particular interest for the empirical verification of the efficacy of frailty.
Evaluation strategy We use two standard metrics in survival predictions for evaluating model performance: integrated Brier score (IBS) and integrated negative binomial log-likelihood (INBLL). Both metrics are derived from the following:
1
8(0, (( 11, (1,3(1)))) >
S(e,
dr. (11)
(T.)
So(1)
Where (t) is an estimate of the survival function (t) of the censoring variable, obtained by the Kaplan-Meier estimate [41 of the censored observations on the test data. l: (0, x [0, 1] + R+ is some proper loss function for binary classification (28). The IBS metric corresponds to being the square loss, and the INBLL metric corresponds to being the negative binomial (Bernoulli) log-likelihood E00. Both IBS and INBLL are proper scoring rules if the censoring times and survival times are independent. 5 We additionally report the result of another widely used metric, the concordance index (C-index), in appendix D Since all the survival datasets do not have standard
'Oteerwise, one may pose covariate-dependen model on the censoring time and use (c)() instead of (t()). We adopt the Kaplan-Mes Meier approach since it's still the prevailing practice in evaluations of survival predictions.
8

---

Table Survival prediction results measured in IBS and INBLL metric (%) on two larger datasets. In each column, the boldfaced score denotes the best result and the underlined score represents the second-best result (both mean). Two models are not reported, namely SODEN and DeepEH, as we found empirically that their computational/memory cost is significantly worse than the rest, and we fail to obtain reasonable performances over the two datasets for these two models.
Model
MIMIC-III
KKBOX
IBS
INBLL
IBS
INBLL
CoxPH
20.40+a.u0
60.02+0.00
12.60=0.00
39.40+o.00
GBM
17.70
52.30
11.81
38. 15.15110
RSF
17.7940.11
53.44
14.46=2000
44..99 00
DeepSurv
18.58:1.92
55.98 +2.43
11.31=00.5
15.28100.15
CoxTime
17.68
52. 02.008+3.06
10.70 00.06
33. 10 40.21
DeepHit
19.80+1.31
59.03.14.20
16.00=0.34
18.6141.04
SuMo-net
18.62+1.28
54.5 +2.97
11.58_0.11
36.61
DCM
18.02#4.49
52.83
10.71
33.24+0.06
DeSurv
18.19
54.69
10.77
33.22
NFM-PF
16.28.0.30
49.18 1.092
11.02
35.1010.22
NFM-FN
17.47
51.48
11.23 $1.23
10.63
32.81 10.14
train/test splits, we follow previous practice [75] that uses 5-fold cross-validation (CV): 1 fold is for testing, and 20% of the rest is held out for validation. In our experiments, observed that single random split into folds does not produce stable results for most survival datasets. Therefore we perform 10 different CV runs for each survival dataset and report average metrics well as their standard deviations. For the KKBOX dataset, we use the standard train/valid/test splits that are available via the pyсoх package 46) and report results based on 10 trial runs.
that standardize continuous features into zero mean and unit variance, and do one-hot encodings for all categorical features. We adopt MLP with ReLU activation for all function approximators, including h, m in PF scheme. and e in FN scheme, across all datasets, with the number of layers (depth) and the number of hidden units (width) within each layer being tunable. We tune the frailty transform over several standard choices; gamma frailty, Box-Cox transformation frailty and IGG(() frailty, with their precise forms detailed in appendix C.3 A more detailed description of the tuning procedure, as well as training configurations for baseline models, are reported in appendix C.3
Experimental setup We follow standard preprocessing strategies 42
46
Results We report experimental results of small-scale datasets in table I and results of two larger datasets in table - The proposed NFM framework achieves competitive performance which is comparable the other ttate-of-the-ae models. In partieular, NFM attains best performance in mean on 5 of the 6 datasets, and is statistically significantly better over all the baselines at 0.05 empirical level on the MIMIC-III dataset.
Allaiion on the benefits of frailty To better understand the additional benefits of introducing the frailty formulation, we compute the (relative) performance gain of NFM-PF and NFM-FN, against their non-frailty counterparts, namely DeepSurv [42 and SuMo-net [56]. The evaluation is conducted for all three metrics mentioned in this paper. The results shown in table The results suggest solid improvement in incorporating frailty, as the relative increase in performance could be over 10%
for both NFM models. A more detailed discussion is presented in section D.5
6 Discussion and conclusion
We have introduced NFM as a flexible and powerful neural modeling framework for survival analysis, which is shown to be both statistically correct in theory, and empirically effective in predictive tasks. While our proposed framework provides a theoretically-prinipllee tool of neural survival modeling, few limitations and challenges need be addressed in future works including predictive guarantees and better evaluation protocols, which we elaborate in appendix D.6

---

7
Acknowledgements
We would like to thank professor Zhiliang Ying and professor Guanhua Fang for helpful discus- sions. Wen Yu's research is supported by the National Natural Science Foundation of China Grants (12071088). Ming Zheng's research is supported by the National Natural Science Foundation of China Grants (12271106).
References
[1] L. Antolini, P. Boracchi, and E. Biganzoli, A time-dependent discrimination index for survival data. Statistics medicine, 24(24):322-3394, 2005.
[2] P. L.. Bartlett, N. Harvey, C. Liaw, and A. Mehrabian. Hearly-tight vc-dimension and pseu- dodimension bounds for piecewise linear neural networks. The Journal of Machine Learning Research, 20(1):2285-2333 22001, 2019.
[3] S. Bennett. Analysis of survival data by the proportional odds model. Statistics in medicine, 2(2):273 1983.
[4] J. Bickel, C. A. Klaassen, Y. Ritov, J. Klaassen, J. A. Wellner, and Y. Ritov. Efficient and adaptive estimation for semiparametric models, volume Springer, 1993.
[5] J. P. Boyd. Chebyshev and Fourier spectral methods. Courier Corporation, 2001.
[6] P. Brennan. cene-environment interaction and actiology of cancer: what does it mean and how
can we measure it? Carcinogenesis, 23(3):88-3877 2002.
[7] J. Buckley and I. James. Linear regression with censored data. Biometrika, 66(3):429
1979.
[8] Cai, J. Fan, Jiang, and H. Zhou. Partially linear hazard regression for multivariate survival data. Journal of the American Statistical Association, 102(478):538-551 2007.
[9] Cai, J. Fan, Jiang, and Zhou. Partially linear hazard regression with varying coefficients for multivariate survival data. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 70(1):141-158, 2008
[10] T. Chen and C. Guestrin. Xgboost: A scalable tree boosting system. In Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining, pages 785-794, 2016.
[11] X. Chen. Large sample sieve estimation of semi-nonparametric models. Handbook of econo- metrics, 6:5549-5635 2007.
[12] X. Chen and X. Shen. Sieve extremum estimates for weakly dependent data. Econometrica, pages 289-314, 1998.
[13] D. D. Cox and F. O'Sullivan. Asymptotic analysis of penalized likelihood and related estimators. The Annals of Statistics, pages 1676-1695, 1990.
[14] D. R. Cox. Regression models and life-tables. Journal of the Royal Statistical Society: Series (Methodological), 34(2):187-202, 1972.
[15] D. R. Cox. Partial likelihood. Biometrika, 62(2):66--276, 1975.
[16] C. Curtis, S. P. Shah, S.-F. Chin, G. Turashvili, O. M. Rueda, M. J. Dunning, D. Speed, A. G. Lynch, S. Samarajiwa, Y. Yuan, et al. The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups. Nature, 486(7403) 2012.
[17] J. Cuzick. Rank regression. The Annals of Statistics, pages 1369-1389, 1988.
[18] D. M. Dabrowska and K. A. Doksum. Partial likelihood in transformation models with censored data. Scandinavian journal of statistics, pages 1-23, 1988.
[19] D. Danks and C. Yau. Derivative-based neural modelling of cumulative distribution functions for survival analysis. In International Conference Artificial Intelligence and Statistics, pages
7240-7256. PMLR, 2022.
[20] A. Dispenzieri, J. Katzmann, R. A. Kyle, D. R. Larson, M. Therneau, C. L. Colby, J. Clark, G. P. Mead, S. Kumar, J. Melton III, et al. Use of nonclonal serum immunoglobulin free light chains to predict overall survival the general population. In Mayo Clinie Proceedings, volume 87, pages 517-523. Elsevier, 2012.
10

---

[21] L. Duchateau and P. Janssen. The frailty model. Springer Science & Business Media, 2007.
[22] D. Duffie, A. Eckner, G. Horel, and L. Saita. Frailty correlated default. The Journal of Finance, 64(5):2089-2133, 2009.
[23] Etezadi-Amoli and A. Ciampi. Extended hazard regression for censored survival data with covariates: spline approximation for the baseline hazard function. Biometries, pages 181-192, 1987.
[24] Fan and Gijbels. Local Polynomial Modelling and Its Applications: Monographs on Statistics and Applied Probability 66, volume 66. CRC Press, 1996.
[25] D. Faraggi and R. Simon. A neural network model for survival data. Statistics in medicine, 14(1):73-82, 1995.
[26] M. H. Farrell, T. Liang, and S. Misra. Deep neural networks for estimation and inference. Econometrica, 89(1):8812213, 2021.
[27] H. Friedman. Greedy function approximation: gradient boosting machine. Annals of statistics, pages 1189-1233 2001.
[28] T. Gneiting and A. E. Raftery. Strietly proper scoring rules, prediction, and estimation. Journal of the American statistical Association, 102(477):359-377 378, 2007.
[29]
Goodfellow, Y. Bengio, and A. Courville. Deep learning. 2016.
[30] E. Graf, C. Schmoor, W. Sauerbrei, and M. Schumacher. Assessment and comparison of prognostic classification schemes for survival data. Slutistics in medicine, 18(17-18):259--5555 25445, 1999.
[31] R. J. Gray. Estimation of regression parameters and the hazard function in transformed lincar survival models. Biometrics, 56(2):571-576 2000.
[32] S. Groha, S. M. Schmon, and A. Gusev. A general framework for survival analysis and multi-state modelling. arXiv preprint arXiv:2006.04883 2020.
[33] X. Han, M. Goldstein, A. Puli, T. Wies, A. Perotte, and R. Ranganath. Inverse- weighted survival
games. Advances in neural information processing systems, 34:2160-2172, 2021.
[34] M. Hardt, B. Recht, and Y. Singer. Train faster, generalize better: Stability of stochastic gradient descent. In International conference on machine learning, pages 1225-1234. PMLR, 2016.
[35] P. Hougaard. Life table methods for heterogeneous populations: distributions describing the heterogeneity. Biometrika, 71(1):75-85 83, 1984.
[36] Huang. Efficient estimation of the partly linear additive cox model. The annals of Statistics, 1765):153- 1999.
[37] H. Ishwaran, U. B. Kogalur, E. H. Blackstone, and M. S. Lauer. Random survival forests. The annals of applied statistics, 2(3):841-88) -600, 2008
[38] H. Jing and A. J. Smola. Neural survival recommender. In Proceedings of the Tenth ACM International Conference Web Search and Data Mining, pages 515-524, 2017.
[39] A. E. Johnson, T. J. Pollard, L. Shen, L.-w. H. Lehman, M. Feng. M. Ghassemi, B. Moody, P. Szolovits, I.. Anthony Celi, and R. G. Mark. Mimic-iii, a freely accessible critical care database. Scientific data, 3(1):1-9, 2016.
[40] J. D. Kalbfleisch and R. L.. Prentice. The Statistical Analysis of Failure Time Data, volume 360.
John Wiley & Sons, 2002.
[41] E. L. Kaplan and P. Meier. Nonparametric estimation from incomplete observations. Journal of the American statistical association, 53(282):457-481, 1958.
[42] L. Katzman, U. Shaham, A. Cloninger, Bates, Jiang. and Y. Kluger. Deepsurv: person- alized treatment recommender system using cox proportional hazards deep neural network. BMC medical research methodology, 18(1):112, 2018.
[43] W. A. Knaus, F. E. Harrell, J. Lynn, L. Goldman, R. Phillips, A. F. Connors, V. Dawson, W. J. Fulkerson, R. M. Califf, N. Desbiens, et al. The support prognostic model: Objec- tive estimates of survival for seriously ill hospitalized adults. Annals of internal medicine, 122(3):191-205 1995.
11

---

[44] C. Kooperberg, C. J. Stone, and Y. K. Truong. Hazard regression. Journal of the American Statistical Association, 90(429):78 94, 1995.
[45] M. R. Kosorok, B. L. Lee, and J. Fine. Robust inference for univariate proportional hazards frailty regression models. The Annals of Statistics, (000:)448-199) 2004.
[46] H. Kvamme, O. Borgan, and I. Scheel. Time-to-event prediction with neural networks and cox regression. arXiv preprint arXiv:1907.008: 000200 2019.
[47] C. Lee, W. Zame, J. Yoon, and M. Van Der Schaar. Deephit: A deep learning approach to survival analysis with competing risks. In Proceedings of the AAAI conference on artificial intelligence, volume 32, 2018.
[48] A. Murphy. Consistency in a proportional hazards model incorporating random effect. The Annals of Statistics, 22(2):712-731 1994.
[49] S. A. Murphy. Asymptotic theory for the frailty model. The annals of statistics, pages 182-198, 1995.
[50] A. Murphy and A. W. Van der Vaart. On profile likelihood. Journal of the American Statistical Association, 454450);449-465, 2000.
[51] C. Nagpal, S. Yadlowsky, N. Rostamzadeh, and K. Heller. Deep cox mixtures for survival regression. In Machine Learning for Healthcare Conference, pages 674-708. PMLR, 2021.
[52] T. Omi, ueda, and K. Aihara. Fully neural network based model for general temporal point processes. In H. Wallach, H. Larochelle, A. Beygelzimer, F, d'Alché-Buc; E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
[53] E. Parner. Asymptotic theory for the correlated gamma-frailty model. The Annals of Statistics, 26(1):183-214, 1998.
[54] S. Pölsterl. scikit-survival: A library for time-to-event analysis built on top of seikit-leam. Journal of Machine Learning Research, 21(212):1-6, 2020.
[55] S. Purushotham, C. Meng, Z. Che, and Y. Liu. Benchmarking deep learning models on large healthcare datasets. Journal of biomedical informatics,
2018.
[56] D. Rindt, R. Hu, D. Steinsaltz, and D. Sejdinovic. Survival regression with proper scoring rules and monotonic neural networks. Proceedings of Machine Learning Research. Journal of Machine Learning Research, 2022.
[57] J. Schmidt-Hicber. function. The Annals of Statistics, 48(4): 2020.
Nonparametrie regression using deep neural networks with relu activation
[58] S. Shalev-Shwartz and S. Ben-David. Understanding machine learning: From theory to algorithmss Cambridge university press, 2014.
[59] X. Shen. On methods of sieves and penalization. The Annals of Statistics, 25(6):2555-251 1997.
[60] X. Shen and W. H. Wong. Convergence rate of sieve estimates. The Annals of Statistics, pages 580-615, 1994.
[61] M. Stablein and I. Koutrouvelis. A two-sample test sensitive to crossing hazards in uncensored and singly censored data. Biometrics, pages 643-652, 1985.
[62] H. Steck, B. Krishnapuram, C. Dehing--berije, Lambin, and V. C. Raykar. On ranking in survival analysis: Bounds on the concordance index. In J. Platt, D. Koller, Y. Singer, and S. Roweis, editors, Advances in Neural Information Processing Systems, volume 20. Curran Associates, Inc., 2007.
[63] R. L. Strawderman and A. A. Tsiatis. On the asymptotic properties of a flexible hazard estimator. The Annuls of Statistics, 24(1):41-66 1996.
[64] W. Tang, Ma, Q. Mei, and J. Zhu. Soden: A scalable continuous-time survival model through ordinary differential equation networks. Journal of Machine Learning Research, 23(34):1-29, 2022.
[65] A. Tsybakov. Introduction to Nonparametric Estimation. Springer Series Statisties. Springer New York, 2008
12

---

[66] A. van der Vaart, A. van der Vaart, A. van der Vaart, and J. Wellner. Weak Convergence and Empirical Processes: With Applications Statistics. Springer Series in Statistics. Springer. 1996.
[67] A. W. Van der Vaart. Asympiotic statistics, volume 3. Cambridge university press, 2000.
[68] A. Wehenkel and G. Louppe. Unconstrained monotonic neural networks. Advances in neural information processing systems, 32, 2019.
[69] A. Wienke. Frailty models in survival analysis. Chapman and Hall/CRC, 2010.
[70] W. H. Wong and X. Shen. Probability inequalities for likelihood ratios and convergence rates of sieve mles. The Annals Statistics, pages 339-362, 1995.
[71] W. C.-H. Wu, M.-Y. Yeh, and M.-S. Chen. Predicting winning price in real time bidding with censored data. In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 1305-1314, 2015.
[72] D. Yarotsky. Error bounds for approximations with deep relu networks. Neural Networks, 94:103-114, 2017.
[73] Z. Ying. large sample study of rank estimation for censored regression data. The Annals of Statistics, pages 76-99, 1993.
[74] D. Zeng and D. Lin. Efficient estimation of semiparametric transformation models for counting processes. Biometrika, 93(3))627-644) 2006.
[75] Q. Zhong, J. W. Mueller, and J.-L. Wang. Deep extended hazard models for survival analysis. M. Ranzato, A. Beygelzimer, Dauphin, P. Liang. and W. Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages 151 151114 11-1112.. Curran Associates, Inc., 2021.
[76] Q. Zhong, J. W. Mueller, and J.-L. Wang. Deep learning for the partially linear cox model. The Annals of Statistics, 2021.
13

---

A Examples of frailty specifications
We list several commonly used frailty models, and specify their corresponding characteristics via their frailty transform Gg:
Gamma frailty: Arguably the gamma frailty is the most widely used frailty model 48
49
53
69
21, with
1
Go(z) log(1 0z), 2
(12)
When 0 0, Go(r) lo (z is defined as the (pointwise) limit. notable fact of the
gamma frailty specification is that when the proportional frailty (PF) assumption 2 is met, if 0 0, the model degenerates to CoxPH. Otherwise if 1, the model corresponds to the proportional odds (PO) model [3].
Box-Cox transformation frailty: Under this specification, we have
Go(x) (1 x) - 0 0 0.
(13)
The case 0 is defined analogously to that of gamma frailty, which corresponds to the PO model under the PF assumption. When 0 = 1, the model reduces to CoxPH under the PF assumption. (GG(a) frailty: This an extension of gamma frailty specifications like the inverse gaussian frailty [35,, with
[45 and includes other types of frailty
L-a Ga(z) = af
Ox - (1 11+1==)--11100
0, 10,
(14)
In the one-dimensional parameter paradigm, the parameter a is assumed known instead of being learnable. When = 1/2, we obtain the gamma frailty model. When 0, the limit corresponds to the inverse Gaussian frailty.
Satistiability of regularity condition 4.5 In 45 Proposition 1], the authors verified the regularity condition of gamma and IGG(a) frailties. Using similar argument, straightforward verify the regularity of Box-Cox transformation frailty.
B Proofs of theorems
B.1 Preliminary
Additional definitions The theory of empirical processes [66] will be involved heavily in the proof. Therefore briefly introduce some common notations: For function class F, define N (e, F. to be the covering number of F with respect to norm || under radius e, and define Nμ (e, F. 11 to be the bracketing number of FF with respect to norm || under radius e. We use VC (F) to denote the VC-dimension of F. Moreover, constant C.
the notation S b to denote S Cb for some positive
Before proving theorem 4.6 and 4.7 we introduce some additional notations that will be useful
throughout the proof process.
In the PF scheme, define
mm(z) l" choolds
1(T, 6, Z; 2,h,m, m, 8) =ō loggo
+ +T((T) + 6m(()
m()) (" (cilo)dss
Ge
where we denote ge = (0). Under the definition of the sieve space stated in condition 4.3 we restate the parameter estimates as
(hrn. mn. 0n
argmax
Σ
1(T,, Zg; h, m. 0).
MEH.,,E.MM,
ve(n)
14

---

Similarly, in the FN scheme, we define
I(T, 6, Z; = 5log (L, e(ls.2)ds) + (u( - Go (f elle.2
Under the definition of the sieve space stated in condition
4.4 we restate the parameter estimates as
( (t((t:)).B,)
argmax DEV. iE|n]
[[T 8, 2);
We denote the conditional density function and survival function the event time T given Z by (t) and (t), respectively. Similarly, denote the conditional density function and survival function of the censoring time C given by Sciz(() and Saiz(1), (t). Under the assumption that IL C Z, the joint conditional density of the observed time T and the censoring indicator ó given Z can expressed the following:
p(T, 6 Z) =
z (T))
(T
where >112 (T) is the conditional hazard function the survival time 7 given Z.
Under the model assumption of PF scheme, p(T, Z) can be expressed by
pTT, 8 Z; h,, m, 0) = exp (((T, 6, Z; h, m, 0)) Sc)z)
For Φo (ho, mo, 0o) and an estimator ē (h, m. 0), the defined distance deF (ō. 00) explicitly expresses by
can be
[VVeC2, 2;h, vp(T, K; mo, 0o) u(dT d6)
2
dFN (e. vio)
Ez
Here the dominating measure ş is defined such that for any (measurable) function r(T,
/ ((T,8)y(dt d6) [ rTT,8 8 = 1)oT L. r(T, 0)df
Under the model assumption of FN scheme, PTT,5 Z) can be expressed by p(T 8 2; 0) exp (((T, 5, Z; 0)) (
For vo = (No, 0o) and an estimator T0 (i, 0). the defined distance dFN (8)))) can be explicitly
expresses by
dFN ( th (0) = Ez [ VP(T,8 Z: pldT dō)
B.2
Technical lemmas
The following lemmas are needed for the proof of Theorem 4.6 and 4.7 Hereafter for notational convenience, we will use h, m for arbitrary elements in the corresponding sieve space listed in condition 4.3 & for arbitrary element in the sieve space listed in condition 4.4 and for an arbitrary element e
Lemma B.I. Under condition 4.1 terms bounded:
4.3
4.5 for (T,
(0, (0,1) [-1, 1]d, the following
15

---

1. 1(T, 6, ho, mo, 80) with true parameter (ho, mg, 0u)
1(T, 6, Z;h, m. 8) with parameter estimates
(h, m, 0) in any sieve space listed in condition 4.3
Lemma B.2. Under condition 4.2
4.4
4.5
for (T, d Z) E [0, x {0, [-1, the following
terms are bounded:
1(T, 6, Z; Vo, 06) with true parameter (vo, (o)
2. I(T, 6, Z;0,0) 8) with parameter estimates (i, (D,0) in any sieve space listed in condition 4.4
Lemma B.3. Under condition 4.1
4.3 4.5 let (h, m, 0). (h1. m1 01), and (h2, m2g 02) be arbitrary
three parameter triples inside the sieve space defined in condition 4.3 the following inequalities hold.
[K[T, 6, Z; ho, mo, 0o) -
m, 00 S 100 1 +1++||u -illoo + +1mmo -fellox
[I[T, 6, Z: 01) 1(T, Z: hig m2, 02) X x (0, -11+ + 11 Tialloo + || - m21lo..
Lemma B.4. Under condition 4.2
4.4 4.5 let (D, 0) (2,,0), 01), and (D 02) be arbitrary three
parameter tuples inside the sieve space defined in condition 4.4. the following inequalities hold.
[I, 6, Z: ug. 0o) I(T, Z: 2, ~ S 100 0 + 1110 11
[[ 5, 2;0, 81) I(T, Z; D2, 12) S 101 021 + 10 D21100
Lemma B.S (Approximating error of PF scheme). In the PF scheme, for any n, there exists an element in the corresponding sieve space Tu = ho, Mo, Oo). satisfying dpF (Tn00, Φu)
0 72 PTd
Lemma B.6 (Approximating error of FN scheme). In the FN scheme, for any n, there exists an element in the corresponding sieve space nn (Rn NnDo) satisfying deN (FN(TT (Tn tho) =
O
7+4+I
Lemma B.7. Suppose F is a class of functions satisfying that F1 < x: oo Ve > 0. We define Fe || II) to be the minimal number of e-balls B(f,e) = {g : 11g fll < <e) needed to cover FF with radius e and further constrain that F. Then we have
F, S (E, F, 11 -I) S N F. II).
Lemma B.&. Suppose F is a class of functions satisfying that N1 (E, F, (1o)) < for VA > We define Ny(e,F;11 11,0)) to be the minimal number of brackets [1, needed to cover F with || alloo E and further constrain that f FFF, 1 and
5. Then have
Nn(e, 11 - (n F, (1oo) S N F. F, (1o0)
Furthermore, we have (e, (1~m) = F, 1100).
Lemma B.9 (Model capacity PF scheme). Let Fn {(T, Z; h m 0) h E Hn, m 2B
Mn.
}}. Under condition 4.5 with st.
and 8m 23+4'
23
there exist constants Ch and cm
20+1
such that
Nn 1bx)) S - N(cne1/en Hn, (121)) x N (( 1/em Mn. 1122)
Lemma B.10 (Model capacity of FN scheme). Let Gn (I(T, 8, Z; DE Vn, O 0). Under 28
condition 4.5 with sy = 2044+I. there exists a constant er > 0 such that
(e Gnsll ll. (1 1 N(cvs Vn. 112).
B.3
Prools of theorem 4.6 and 4.7
Proof of theorem 4.6 The proof is divided into four steps.
16

---

We denote de = (ho. mo, 6o) and = (h, m. 0). where E Hnmm E Mn and 1 E θ, For arbitrary small 0. we have that
Step 1
inf dop (0.00) )e
E [[T, 8, Z; ho, mo, 0o) -([T, Z; h,
Ez ET,aiz log Z; ho mu. 00) (T, Z;h, m. ]]]
inf des (0.00) 2=
inf
Ez DKL (P || (Po0.2)
(0.00) 25
Using the fact that DKL.
(Paz ll Pob.2) 2 2H22(P5,2 8, P.o.z). Thus, we further obtain that
inf (0(0000)226
E I(T, 6, Zz ho mo. 0o) 1(T, 6, Z; m. 0)
inf doFl (5.bh))
2Ez
(Pazz
(Pzz || Pdo.z)
dBF (ō,
inf dpr (0.00) )2=
2 2
Step 2
Consider the following derivation.
mo, Var [(T, 6, ho mo, 0o) I[T, 2; m, 22.
sup dos (5e
(IT,,0, 2; ma, ma. 0o) I(T, Z; 7; mr, 0)
< sup dor (0)
E
Ry, [ogg(T, 8, Z; ho, mn, 0o) (og (T, (T, m, 0)]²
sup doF (o)
Ez /
?) uldt dō)
mo p(T, 8, Z3 Mo, 0o) p(T, 5, Z; h mr, 0)
=4 sup dop (0.0) Se
p(T, 0, Z; ho, mo, 0o) log
By Taylor's expansion on logı, there exists E(T, 8, 2)
between på (T, 2; ho, mo, 00) and
pa (T, 8, Z; m. 0) pointwisely such that
p(T, 8, Z; ho, ma, 00) (108 V
p(T, 5, Z;ho, mo, 0n)
p(T, 8, Z; f,
2
pTT, 8, Z; ho, mo 0v)
log Vp(T, 8, Z; hog my, 0o) - log V p(T, 6, Z; m.
P(T,8, ZiZho ho, mo. 00) S(T, 8, 2)²
(Vp(I, Zi ho, mo, 0o) V((T. 8, Z; Z;h, m,
Since
p(T, 5, Z; ho, mp 0o) mo)) p(T, 0, h;
by lemma B.1 I(T, Z: ho ma, 0)) and I(T, 5, 7: hr m 0) are bounded among 7] x {0, 1} x
[-1, 1ld uniformly on all = m, 0). Thus, there exist constants C₁ and C₂ such that 0 < C p(T, 8, Z; ho, mo, (oT/ 6, Z;h, m < C2. This leads to the fact that
17

---

P(T, 6, Z;ho, mo. 80) E(T,8,27 is bounded. We further obtained that
2
p(T, 6, Z: ho, mg 80) log V pTT, hn mp. 80)
-log VP(T, Z; m.
x VP(T, ò. Zi ho. mo, mo,0u) V P(T, Z. m.
Thus, we have that
Var([ Z; ho mo, 00) ((T 8, Z; m. 0)
sup
do. Se.
/ /VEIT VP(T. 8, Z; ho, mo, mo, 0o) v v((T,O, 2: m. 0)
12
sup de (5.d0)
Bz
MLT x d6)
dpF (2.00)
sup doF( (ō.0))
S 5.a
Step We define that Fr {1(T, 8, Z; h. m. I(T, 6, Z; ho. T. T, h Hr, m Mn, 1 E @). Here (o, T ma TT 00) have been defined in lemma B.5 Obviously, we have that (loo) = log N1,
log Fs, ||
lloo), where F is defined in lemma B.9 By lemma B.9
we further have that
log Fn. (loc) S log + log (C Ha, -1l1)2) + log N (cmm2//5m Ma. (112..
According to 12] Theorem 7], under condition 4.3 we have that the VC-dimension of Hn and Mn
satisfy that VC (H.) S nगvत log" n and VC (MM S n u.r log n. Thus, we obtain that
VC (Ha) log S nota log nlog E
log N
2/5b Ho, l12) x
and
log (Cme Mn, (1)1112) (1)) x VM log S nnzaa log n log _
Thus, obtain that log (D, Fr, || (1oo) S n log nlog 1
Step 4 By the Cauchy-Schward inequality, we have that
v E (TT,, Z; n, T. mo, n. 0o) 1(T, 8, Z; ho, mg, 0o) < [B( 8, Z; ho, m, To) - I(T, 5, Z; ho, mo. 0o)) 1 Similar to the second part and by lemma B.5 we further have that
E [(TT, 8, Z; n. by, my On) - I(T, 8, Z; by, mp, 0o)] S Vdpe (un Qu) S n -
Now let
B
log log n
2
2,3 + 2d
log
By Step 1,2,3 and 160 Theorem 1]. we have
dpF
00) = max (n -", ddpp (0o 00),
v EE [[T Ta, Mo, 0o) - -(T,, 6, 2; ho, mo, bo]]
By lemma B.5 dpF (T do, Φ0) = O(n #4d
By Step 4, / [T(T, 8, Z; Thho, nn mo, 00) I(T, 8, Z; ho, mp, 00) = =0 =0(n (n-))()) ). Thus, we have dpF
(ons s0 = O(n حाد log² n) = Ö(n I8+22
Proof of theorem 4.7
The proof is divided into four steps.
18

---

Step 1
We denote vu = (vo. 0b) and eे = (D. -0,00) where v E Vn and O E O. For arbitrary 0 <<< 1,
we have that
inf dru (8,00)
I3 [(T, d, 2; - 11T 2;
ET,aaz logp(T,6) 8 2; logp(T,8 Z: 02 R]]
inf diu(W-No)2z
Ez
inf
Ez DKL (P7.2 III| P0..22
deu (n.(voo)2e
Using the fact that K [L(
(2,.2) 2 2H (T. 3.2 || 22).2). Thus, we further obtain that
322
inf dav (v.u)))
E [(T 6, Z: Mu, Bo) L(T, 6, 2:0,0)
inf dru )2
2 H2 H2 $, P
2 inf dev (v.V))))
dPN (5. 1b)
25a.
Step 2
We consider the following derivation.
sup (E
Var
I(T, (T 6, Z; 0o) ((T,6,
sup dew (.v0) Se
E
[[T,, Z: Bo) I(T, -VTT, (;
Er,oi (logp(T,o, 0, 2; 0o) log (TT 0,, ) I]]
sup den (on( Sc
Ez
Ez // p( 0o)
p 210, Onl
=4 sup deu e0)
uldT x dō)
p(T, 6(T,6,2;37 2:0,0)
By Taylor's expansion on logæ, there exists v(,, Z)
between Vp(T, 6, 2; 0o) and
V PT 5, 2:0,0) 8) pointwisely such that
(T, 8, Z; Z; B00)2
p(T, 5, Z; (lol
V
P(T, 2;2,0)
2
= p(T, 2: 00)
log V p(T, 5, Z;vo, 00) log V p(T 2;0,
p(T, 6, Z; 0o)
VV(TT Z; up, 0o) (TTT,, 2
6
Since p(T, 5, Z; 0o) 0p)/p(T, 8, 2: = J(T.AZ.0000 0o) (TT66 by lemma B.2
1(T, 6, Z; 0o) and 1(T, 5, 2; 0) are bounded on - x {0.1} x 1_ uniformly (D,0). Thus there exist constants C3 and C4 such that 0 C3s
for all 1
p(T, 8, 2o 00)0 6, Z; S C.. This leads to the fact that p(T, 8, Z; vo, 0o) ,ा,,
bounded. We further have that
PTT, 5, 2: 9o)
log v P(T, 8, 8, - log V p(T, 5, 2:0. 0)
$ DOTT,,2:0000 00)
p(T (,T,8,220,0)
19

---

Thus, we have that
Var [(TT, 6, 2: 80) ((T,6,2;000)
sup (o.ows
Ez / 07:00002 2o0 p(
sup den (.v
uldT x d9)
dFN vi)
sup dev dnn(o.o)se
EE²
Step We define that 9n (1(T, 2: I(T, 7: No Tn E Vn, 0 E θ). Here (TnMo, Tnlo) have been defined lemma B.6 Obviously, we have that log Nnfe, Gn. l|x)) =
log (E, lloo), where is defined in lemma B.10 By lemma B.10 we further obtain that
log (ID Gn. |l (l z log
log N (c 1/2, 1/4,
Vn. 1122)
According to 12 Theorem 7J, under condition 4.4 VC (V,) S n n442i1 log" n. Thus, we obtain that
we have that the VC-dimension of Vn, satisfies that
d-1 log (cne! Vn, 112) S (Vn)n) log S n=+ log n log
1
Furthermore, we get that log (e, Gn. || (1oo) S n niili log log 1
Step 4 By the Cauchy-Schwartz, inequality, we have that
(E, 8, Z; No, 0o) ((T, 6, vo 80) S E (I(T, 6, Z; T, T. [(T, Z: 0n))
Similar to the second part and by lemma B.6 further obtain that
v E[[ Z; Mn Bo) I(T, 6, Vo 0o)l S v dFN (un S n 2012
Now let
8
log log n log
2
23 2d +
By step 1,2,3 and Step 1,2,3 and (60 Theorem 1].
dFN @nn+ vb)
max (n "n (T to. vho),
[[ 5, Z; 4o, 0o) 1(T, 2;vo. 0n)]
By lemma B.6 dFN (nn No. vho) = On गтकका
By Step 4, / [[ 8, Z: n vo, 0o) ITT, 8, Z: Vo, vo) O(n 22-20+2). Thus, we have
dEN
Vna Mu = O(n 20424+2 log n) Ö(n 20+22+2).
B.4 Proofs of technical lemmas
Proof of lemma B.I Since ho(T) E W r]) and ma(Z) E (1-11-11))) we have that
ho(T) < M, m 7) M and eml(z) fo ho(o)ds Te2M Let B = [0, Te2M], we have that
[[T, 6, Z: ho, mo. 0o)l
mo) (" cholo)dss ds + +hoo(T) + (2) + Goo emo(2) 1" cho(e)dss
< log gao
< 2. + sup llog go (x)) + sup Goo (x) seB
By condition 5, we have that ((T, 6, Z; ho, mo, 0p) is bounded for (T, E
The proof of the boundness of (TT 8, Z; h m, 0) is similar.
20

---

Proof of lemma B.2 Since 2) E (10, x 1-1, we have vo(T, I
M and
Sor el(,,2) ds TM Let B - [0, Te we have that
1(T,, 6, Z; 80 00)
( (0000022) c00(2:2) ds ds)
[To (T,,TT,) 2)1 + Goo (Jo² coolaz2)ds)
log Co
S M sup =EB
log C'g (3:)
sup =EB
|oo (x))
By condition 4.5
we have that I[T, 6, 2o 0o) is bounded among (T, E [0. x {0,1}
[-,, The proof of the boundness of [[T, 5, Z; 12 0) is similar.
Proof of lemma B.3 By definition, have that
[[ 6, Z; ho, mo. Bu)
1(T, 6, Z: m. 0)
mm)(2) !'"
T
emi2) ''" Moodss ho(T)) - -((T)
log 98.
co(l) ds
logga ga
L'" ch(s) ds -G c S"" (iioids)
(2) - -m((Z) + Goo e
Let B = 00, max(c2M M. +m }]. By Taylor's expansion, we can further show that I(T, 6, Z: f.
[(T,, 6, Z: ho: mo. 8o)
a log ga
0 0
sup
0g
0
i log 97 (7)
mmllz) {o
-
sup 8E8
0g
(Gg(2)
+ (ho(T) R(TT)) + (mo(Z) m(2)) 8
-100 100-0
sup
0z
(Gz(T)
(mol))) 6" (hoss) do (2) 1' Moddd
sup 6EB
Again, by Taylor's expansion, we have that
e (" bo(s) - mi( Sod p (ss)
mmll2) (" (cholo) - -colo) + (emo(Z) (
c2122)) S²" chio)ds
S M max(
Ino - T
+ Mo. (aax(M || mllxo.
Finally, we obtain that
i log (8) [T, 8, Z; ho. mu, 00) - UTT, 5, Z: m.
M re ma ( Tolloo + Te (max(MM,MM) || - melloo
< sup
0z
dEB,iEB
i log 98 (i)
0g Ø +|bu(T) A((T)] + +lmo(Z)- - m((2)
+ sup
0g
(G)(
+ sup
00
0g
d
(Gg(ā)
M Te max MMM Rolloo + Mr. mmax(M.Ma))lng - ml110
+ sup
0g
de0,FEB
21

---

Taking supremum on both sides, we conclude that
[I 6, Z; ho, mo: 0b) I(T, 2: m. 0l 2 10u
0 11/g - -11100 + 1|mo
m/loc.
The proof of the second inequality similar.
Proof of lemma B.4 By definition, we have that
[[T, 6, Mo Bu) 1(T, 6, Z;0,0)
-log 96 (f" cpla2)dss)
< log 9oo
+ (oo(2, -1)((; ()r,
+ Geo (%"" cro(ai)dss ds
Gā (f" pla ds
T
Let B (0, max(e )....)). (.M)) By Taylor's expansion, can further show that
[[, 8, Mo 8o) - 1(T, ((T, 6,
i log 98 (ō)
Olo gg(
l'd culla2)ds 10" ala
0 0
< sup dE0.EEB
+
sup
0z
6E
aG
+ (T, 7) - (T, 2)1
+ sup
100 Ø
Og
iE0,FEB
(Gg(ā)
cu(a. ds
((a.2)
+ sup
Again, by Taylor's expansion,
Lo cola.2) - S¹ dls S max
D1100
Finally, we obtain that
[(T, 6, 2; 0o) I(T, 2;0,0)
3 log 9a (ō)
i log ga (z)
< sup
00
0
+
sup
ma.( 1110 1
0g
0z
8EB
0
(G)
+ +110 (T - (T, 2)
+
sup
0o - 0
0g
8
(G₂(i)
+ sup 8
Te max MM.
Taking supremum on both sides, we conclude that
[M 6, Mo, 0n) I(T, 8, 2:0, 2; S 10p - -0 + 1|v0 D2
The proof of the second inequality is similar.
Proof of lemma B.5 f such that + bolloc
According to 772 Theorem 1], there exist approximating functions πo- and
Y1 and || molloc O (n 3+d Let π. 1o Te'
22

---

Mu f and TT = 80. We have that
dpF (TnDo (b))
/ | Tnho, Tn T,, - Vp(T, 8
Ez
Z;
Z; ho, mo. 00) x dô)
/
Ez
x d6)
-
1
x Ez | (azz(T))-5502(17)) x dō)
By lemma
B.1 and
B.3
we have that
||e 3: - - 1l
2 1 - foll + | - holloc + || - molex
=0 (πス
Since (T) < (T) + and (Tiz(T) S 1, we also have that
Ez V
( x d6)
/ (1+ dō)
Ez
V2 2T.
Thus, obtain that dpF (Tndo (bo) 0 JFu Jin
Proof of lemma B.6 According to
72 Theorem 1], there exists an approximating function D' such
that 1|0" - volloo O (n गहका
Let TT = D and T,00 = 0o- We have that
dFN (TT vio)
Rz [I /IVeTT, 6 2; nI Vo, Tolo) VP(T,ō Z; 7: On) M(dT d6)
Ez [ [edtrrraeemmmmmm -
- 33:0015522,000000) (cru(T)=-Sc2((()uudt x dd)
1 (T Zim.
1 (
|f
fonz/Th)-SSsp2(T99 d6)
1 Bzz
2
By lemma
B.2 and B.4
we have that
e ( ((0,))), 2:7,66.220000)
I || 0g - 0oll + || -olloo
n काम
Since [T) < fc]z((T) + and (T)z(T)
< 1, we also have that
Ez
(tc2TT)-SScp2TTyuuua x do)
Ez
[
[(1+ x dō)
V2+ + 2r.
Thus, we obtain that dEN (T, vho) = =0 - 8+041
23

---

Proof of lemma B.7 The left inequality is trivial according to the definition of covering number. We need show that the correctness of the right inequality.
Suppose that we have [B(gs: 5),;;= 1..., N, where = N F, || [I), are the minimal number of 5-ball that covers F. Then there exists at least one Si E such that Ji B(g,,)). Consider the following E - balls [Bi s)}, = 1. I...,NN. For arbitrary f E F F ((g., 2), we have that Ilf fill = flf gill + [[fe gell S E, Thus [B(fs, =)),, ,..., N forms a e-covering of F. By definition, we have that (N, F. -|- II) S F. I)).
Proof of lemma B.8
The proof of the first two inequalities follows exactly the same steps of lemma
B.7
Here we just need mention the rest of the statement that Ni(e, || l1 NNG F. lbo))
We first choose a set of 5 5-covering balls { [B(ft. 5)), i 1...., NN. where N1 _ F. Il lce).
Now we construct set of brackets {14, u,]}, é = 1 -1.... .N.. where 4 = f. and u4 f. Noting
that the bracket {[4, u|]} is exactly the same as B(f, 5), The set {14, 2.] = l...., N1 covers F.
which leads to (e, F1 ||x) S F. ||. lloo). Likewise, we have that (e, F, 11x)) 2
(5 FI (l Consequently, we have that Ne, ll (1 - (G, F, lloo))
Proof of lemma B.9
By lemma B.8 first we have that NI Fr. 1l (1o0) S Nole, Fr, ll 12). By
lemma B.3, there exists a constant oj such that for arbitrary h1, lr2 E Hn.mi m2 E Mn and
0,, 02 E O, we have that
[I(T, 6, 2;7,, mi, 01) - I(T, 6, m2, 02) S c1 - 022 + 17 - Rial + Mi m2 mool
which indicates that as long as (01 021 S 3e lh Tuzlloo S fet and || mallas S م we have that []<(T, 6, Z: hu, my, 01) I(T, 5, Z; hug m2, 02) S Consequently, we have that
3es
(n Fn, ll )1x))= S No 3c1
0 || 11 Ni - Hn. 11o0)) x [ lloo). 3e1
3e1
Since θ is compact set on R., by lemma
B.8 and traditional volume argument, we have that
( 3ry O, 11oo) S N beg O, l| (00)
For Np Be, Hn. || (l by lemma B.8 we have that N (ST Hn, || |loo)) Ñ ( Hn: || l.
By 12 Lemma 21, there exists a constant > O such that [7 - hzlloo S ez Th - hall@" which
leads to ( Hn, |l (1oo) S N (136)12)) Ho, || l12). By lemma B.7 we further have that
N 21/08 (370 Hn, Il (12) S N el/on
2(36)10) Hnoll l12). Let Cr. 2(3c)
We have that
Ni 3e. Ha, || 11x0) < N (Ccnet/es, Hn, -122)
Similarly, there exists a constant cm > 0 such that ND M 1loo) S N(cmm Mn: ll l12).
Thus, finally we can obtain that
Nn(e,Fmi |ll.-l (lx)) x N(cne Hn, || (12) x N NCcnet/an MMn -12).
11=0) S No Gn, 1l (loo). By lemma B.4 there exists a constant Cg > 0 such that for arbitrary 22, 02 E Vn and 0 O2 E O we
Proof of lemma B.10 By lemma B.8 first we have Nj Gn. |l
have that
| - I(T, 5, 2; 02)1 S <cal - 021 + ||0, - 1pii-b2llad 12
which indicates that long 10, 0₂ S 203 and || D211o0 I
E. we have that
203
[[, 8, 2;2, 02) ((T) L; 0z) )l S Thus, we have:
Nn(e, En, 11 (loo) < No120 20g"
8
O. (1 x NoL, 6
2cg 20s"
Vn. ll 11bo).
24

---

Table 3: Descriptive statistics of benchmark datasets
METABRIC
RotGBSG
FLCHAIN
SUPPORT
MIMIC-III
KKBOX
Size
1904
2232
6524
8873
35953
2646746
Censoring rate
0.423
0.432
0.699
0.320
0.901
0.280
Features
9
7
14
26
15
Since è is a compact set on I. by lemma B.8 and traditional volume argument, we have that
2 20s 0, 11s0) S NI dey O, ll (
For (2 Vn >1 by lemma B.8 we have that N
2ng Vn, [1~0)) (2 Vnl 1120). By
12 Lemma 2), there exists constant 0q >> 0 such that 1 i3 oo S 0 22l12" which leads to
( 20g" Vn, 1111110000 ) S ( (20gea))/= Vn. -1.2))). By lemma B.7 we further have N (20ge4) Vns
e1/"
el/au
ll2) I N 21/00
2220004))1r0 Vns ll [12). Let cu 2(20gg4) we have that Ni = Va, 11,0) S
2cs
(Ccoet/>,, Vn:11 112).
Thus, finally we can obtain that
1 N Vns 112). Vn, (lxo) S E
C Experimental details
C.1 Dataset summary
We report summaries of descriptive statistics 6 benchmark datasets used in section 5.2 in table
C.2
Details of synthetic experiments
Since the true model assumed to be form, we generate event time according to following transformed regression model [18]:
log HTT) (2) +e,
(15)
where (t)) Jee chs with h defined in 2 The error term is generated such that e has cumulative hazard function Gg- The formulation (15) is the equivalent to (2) [88 17 45). In our experiments, the covariates of dimension 5, sampled independently from the uniform distribution over [0, 1]. We set (1)) = and hence H(t) = e'. The function form of m(() is set to be m(Z) sin( (Z, B3) (sin(Z), where 3 (0.1, 0.2, 0.3, 0.4, 0.5). Then censoring time C is generated according to
log HCC) -m( ccs
(16)
which reuses covariate 7, and draws independently noise vector Ec such that the censoring ratio controlled at around 40%. We generate three datasets with E (1000, 5000, 10000) respectively.
Hyperparameter configurations We specify below the network architectures and optimization configurations used in all the tasks:
PF scheme: For both in and 7ú we use 64 hidden units for n = 1000, 128 hidden units for n = 5000 and 256 hidden units for n = 10000. We train each model for 100 epochs with batch size 128, optimized using Adam with learning rate 0.0001, and no weight decay.
FN scheme: For both D we use 64 hidden units for n 1000. 128 hidden units for 5000 and 256 hidden units for n = 10000. We train each model for 100 epochs with batch size 128, optimized using Adam with learning rate 0.0001, and no weight decay.
25

---

C.3 Details of public data experiments
Dataset preprocessing For METABRIC, RotGBSG, FLCHAIN, SUPPORT and KKBOX dataset, we take the version provided in the русoх package [o]. We standardize continuous features into zero mean and unit variance and one-hot encodings for all categorical features. For the MIMIC-III dataset, we follow the preprocessing routines in [55] which extracts 26 features. The event of interest. defined as the mortality after admission, and the censored time defined as the last time of being discharged from the hospital. The definition is similar to that in [64]. But since the dataset is not open sourced, according to our implementation the resulting dataset exhibits a much higher censoring rate (90.2 2% compared to 61.0% reported in the SODEN paper [64)). Since the major purpose of this paper is for the proposal of the NFM framework, We use our own version the processed dataset to further verify the predictive performance of NFM.
Hyperparameter configurations We follow the general training template that uses MLP as all nonparametric function approximators (i.e., m and h in the PF scheme, and in the FN scheme), and train for 100 epochs across all datasets using Adam the optimizer. The tunable parameters and their respective tuning ranges are reported as follows:
Number of layers (network depth) We tune the network depth L {2, 3, 4]. Typically, the perfor- mance of two-layer MLPs is sufficiently satisfactory.
Number of hidden units in each layer (network width) We tune the network width W E (2*, 5 10).
Optional dropout We optionally apply dropout with probability p (0.1, 0.2,0.3,0.5,, 0.7.
Batch size We tune batch size within the range {128, 256, 512]. in the KKBOX dataset, we also tested with larger batch sizes {1024}.
Learning rate and weight decay We tune both the learning rate and weight decay coefficient of Adam within range {0.01, 0.001, 0.0001).
Frailty specification We tested gamma frailty, Box-Cox transformation frailty, and IGG((a) frailty with a E (0, 0.25, 0.75). Here note that (GG(( is equivalent to gamma frailty. We also empirically tried to set to be a learnable parameter and found that this additional flexibility provides little performance improvement regarding the datasets used for evaluation.
C.4 Implementations
We use pytorch implement NEM. The source code is provided in the supplementary material. For the baseline models:
We the implementations of CoxPH, GBM, and RSF from the sksurv package [54] for the KKBOX dataset, we use the XGBoost library [10 to implement GBM and RSF, which might yield some performance degradation.
We use the pyсox package to implement DeepSurv, CoxTime, and DeepHit models.
We use the official code provided in the SODEN paper [64] to implement SODEN.
We obtain results of SuMo and DeepEH based our re-implementttions.
D Additional experiments
D.I Recovery assessment of m2)) in PF scheme
We plot empirical recovery results targeting the function in in figure The result demonstrates
satisfactory recovery with moderate amount of data, i.e., n 2 1000.
D.2 Recovery assessment of survival functions
To assess the recovery performance of NFM with respect survival functions, we consider the following setup: under the same data generation framework in section C.2 we compute the test
feature Z as the sample mean of all the 100 hold-out test points. And plot S Z) against the ground
26

---

N 1000
5000
10000
2
grreured crushs colinnated
gritand Irutt
Eroaand truth
cutimnttdd
eotimnated
Mydee Ma I yp Malhet
"y WWy Mh
1.00
100
1.00
Indox
andex
index
Figure 2: Visualizations of synthetic data results under the PF scheme of NFM framework, regarding empirical recovery of the m function in (2)
N 1000
N 5000
N' 10000
troth retimosd
giound sermoss
troth
0.75
0.75
rermund
801223
0.25
0.25
N 1000
N 5000
10000
ground estinone
ground ectmmed
ground
0.75
cotinned
810123
850123
0.25
0.25
0.25
0
3
Figure 3: Visualizations of synthetic data results under the NFM framework. The plots in the first row compare the empirical estimates of the survival function S 2) against its true value with Z being the average of the features of the 100 hold-out points, under the PF scheme. The plots in the second row are obtained using the FN scheme, with analogous semantics to the first row.
truth S(t|2) regarding both PF and FN schemes. The results are shown in figure 3 The results suggest that both scheme provides accurate estimation of survival functions when the sample size is sufficiently large.
D.3 Numerical results of the synthetic experiments
we report the relative integrated mean squared error (RISE) of the estimated survival function against the ground truth and list the results in table The reults suggest that the goodness of fit becomes better with larger sample size. Moreover, since the true model in the simulation is generated as an PF model, we found PF to perform slightly better than FN, which is reasonable since the inductive bias of PF is more correct in this setup.
Following 755
Table 4: RISE of the estimated survival function in synthetic experiments
N
1000
N
5000
N
10000
NEM-PF
0.0473
0.0145
0.0137
NFM-FN
0.0430
0.0184
0.0165
27

---

Table 5: Survival prediction results measured C-index (%) on all the 6 benchmark datasets. In each column, the boldfaced score denotes the best result and the underlined score represents the second-best result. The average rank of each model reported in the rightmost column. We did not manage obtain reasonable results for DeepEH and SODEN on two larger datasets MIMIC-III III and KKBOX, and we set corresponding ranks to the worst on those datasets.
Model
METABRIC
RotGBSG
FLCHAIN
SUPPORT
MIMIC-I1I
KKBOX
Ave. Rank
CoxPH
1..421
66.14 1.46
79.09
56.89
74.9110.0
83.01 83.0110.00
11.33
GBM
64.02+1.79
67.35-1.16
79.47
61.46+0.80
75.20+0..0
85. 85.8440.00
7.17
RSF
1..4
67.33_1.34
78.75-1.07
61.63
75.47
85.7920.00
8.00
DeepSurv
63.95
67.20_1.22
19.04.11.14
60.91
80.08
85.59 20.06
8.50
CoxTime
66.22+1.69
67.41-1.35
78.95=1.01
61.54+0.87
78. 78.78+0.02
87 31 +0.24
5,00
DeepHit
+6.3
66.38
78.48.11.0
63.20=a
79. 79.1640.50
86.12
6.50
DeepFH
66.59+2..0
67.93
78.71
61.51 +:1 1.00
6.33
SuMo-net
64.82+1.80
67.20
79.28 11.02
68.18 18.78
76.23
84.77
7.00
SODEN
6-1.82+1.05
66.97
79.00
61.10 +0.59
10.17
SurvNode
64..4444.91
67.30 67.30_1.66
76.11 26.11 20.98
55.37 55.77 +0.77
11.50
DCM
15.7611.25
66.75 1.35
78.61 18.79
19.1910.99
76.44110.4
83.4810
8.33
DeSurv
65.88
67.30 67.30_1.45
78.97
61.47
80. 10.37 10.30
86.11
5.67
NFM-PF
64.98
67.77 1.35
79.45
61.33
79.56
86.23.0.00
NFM-FN
11.82
67.73 1.29
10.93
10.41
80.18 10.20
86 86.61
Table 6: Relative improvement of NFM models in comparison to their non-frailty counterparts, measured in IBS, INBLL, and C-index.
Dataset
NFM-PF DeepSury
NFM-FN SuMo-net
IBS
INBLL
C-index
IBS
INBLL
C-index
METABRIC
+1.33%
+1.56%
+1.61%
+2.30%
+3.08%
+2.79%
RotGBSG
+1.11%
+0.95%
+0.84%
+0.62%
+0.40%
+0.79%
FLCHAIN
+1.29%
+1.32%
+0.52%
+0.20%
+0.27%
+0.01%
SUPPORT
+0.31%
+0.23%
+0.69%
+2.22%
+1.76%
+0.05%
MIMIC-III
+12.38%
+12.15%
-0.64%
+6.18%
+5.56%
+5.18%
KKBOX
+2.56%
+0.51%
+0.75%
+8.20%
+10.38%
+2.17%
D.4
Performance evaluations under the concordance index (C-index)
The concordance index (C-index) [ is yet another evaluation metric that is commonly used in survival analysis. The C-index estimates the probability that, for random pair of individuals, the predicted survival times of the two individuals have the same ordering as their true survival times. Formally, C-index is defined as
C S(T, Z.) < S(T, 22) T, =1
(17)
We report performance evaluations based on C-index over all the benchmark datasets in table 5 From table 5] it appears that there's no clear winner regarding the C-index metric across the 6 selected
datasets. We conjecture this phenomenon to be closely related to the loose correlation between the C-index and the likelihood- baseed learning objective,
observed in [56]. Therefore compute
the average rank of each model as an overall assessment of performance, as illustrated in the last The results suggest that the two NFM models perform better on average.
column in table 5
D.5 Benefits frailty
We compute the (relative) performance gain of NFM-PF and NFM-FN, against their non-fran non-frailty counterparts, namely DeepSurv 42 and SuMo-net (56) based on results in table table and table
5. The results are shown in table 6 The results suggest solid improvement in incorporating frailty, especially for IBS and INBLL metrics, as the relative increase performance could over 10% for both NFM models. For the IBS and INBLL metrics, the performance improvement is consistent across all datasets. The only performance degradation appears on the MIMIC-II dataset evaluated under C-index. This phenomenon is also understandable: Since the DeepSurv model utilized variant
28

---

of partial likelihood (PL) for model training, as previous works 62] pointed out that PL type objective closely related the ranking problem. As C-index could be considered certain type of ranking measure, it is possible that DeepSurv obtains better ranking performance than NFM-type models which are trained using scale-sensitive likelihood objective.
D.6 Limitations
In this section we discuss the limitations of this paper form both theoretical and empirical standpoints.
Theoretical limitations As we have established formal statistical guarantees regarding the esti- mation properties of NFM, the guanrantees do not necessarily imply that NFM perform well on prediction tasks under metrics such as IBS and INBLL. Following the spirit of classical learning theory [58], for prediction guarantees it is ideal to directly optimize the underlying metric or its surrogates, which difficult in survival problems since the metrics involve both a model over the event time and working model the censoring time. So far have noticed, the only effort that aims to address this issue is the method of inversely-weighte survival games [33]. However, the authors in [3] did not provide rigorous thhooeec statements, which is promising research direction for future works.
While NFM is shown to be a competitive survival model for prognosis empirically, we have observed from the empirical results that we can hardly obtain statistically significant improvements over the baselines, which is a common problem exhibited in previous works neural survival regressions [75 56 We conjecture that this phenomenon is primarily due to two facts: Firstly, there lacks open-to-public large-scale survival datasets that allows scalable evaluation of neural survival models. Secondly, for most of the current available datasets, there are no authoritative train-test splits. Consequently, most experiments are done using cross- cros-validatiin on moderate scale datasets, causing the resulting variability of the modeling algorithms to be relatively large. Therefore we suggest the survival analysis community to release (in a manner) more large-scale, sanitized datasets equipped with standard train-tess splits, which we believe will greatly benefit the state-of- the-art for neural survival modeling.
Empirical limitations
29


Reviews:
[{"review_id": "X8s6RvDPRb", "reviewer": "Reviewer_U6h7", "summary": "The authors propose a framework for survival regression called neural frailty machine. They have shown that most of the existing methods are a special case of the proposed framework. Also, they have drawn statistical convergence guarantees for the proposed model. The experiments show marginal improvement compared to existing models in terms of some metrics.", "questions": "- What is the intuition behind introducing PF? It is mentioned (in third page, line 127) that \"Proportional-style assumption over hazard functions has been shown to be a useful inductive bias in survival analysis\". If this type of function approximation has been already shown deployed in survival analysis, then it wouldn't be very novel.\n- Checking results, sometimes FN performs better, while PF works better other times. Why? Are there any explanations for each dataset?\n- Is it possible to report numerical results for Figure 1, so we can compare PF and FN?\n - Addition of frailty variable (w) has marginal effect (<~1%) on most of the dataset (Table 5 in Appendix) while having huge effect on one the datasets. Any explanations for this discrepancy?\n\nOverall, I think more discussion around numerical results is needed.  ", "limitations": "Authors discussed some the limitations in the appendix.", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper is well-written and has clear motivation. The overall approach is intuitive and makes sense (although I haven't checked the proofs). The main contributions are introducing frailty variable and the the statistical guarantees drawn for the proposed approaches.", "weaknesses": "The main weakness is lack of discussion/explanation around the results and the two proposed approaches (PF and FN). More on this in the next section. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for providing insightful comments. Below we address your specific points  \n## Q1: The intuition of proportional frailty model\nFirstly, we would like to emphasize that random effect (called frailty), which serves as a principled tool to model *unobserved heterogeneity*, has played an important role in modern survival analysis. There has been significant and active research concerning the addition of random effect to survival models. The random effect can describe risk or frailty for distinct categories, such as individuals or famlies. There eixst literature on providing theoretical and practical motivation for frailty models by discussing the impact of heterogeneity on analyses. Below is a quote from [1]:\n> It is a basic observation of medical statistics that individuals are dissimilar ... Still, there is tendency to regard this variation as a nuisance, and not as something to be considered seriously in its own right. Statisticians are often accused of being more interested in averages, and there is some truth to this.\n\nIn the influential article [1], Aalen showed that with frailty, the population relative risk can go from r to 1/r over time, even though the true relative risk stays at r. Frailty is shown to be very useful to adjust for heterogeneity in real data. \n\nSecondly, while the algorithmic extension of incorporating frailty is straightforward, **it is highly nontrivial to establish statistical guarantees**. In this paper the goal is to introduce a rigorous statistical model with provable guarantees, which was lacked in many previous works. As was pointed out in the paper, the current theoretical developments on neural survival models rely on the fact that the loss function is well-controlled by the L2 loss [2, 3], which is not directly applicable to our model due to the flexibility in choosing the frailty transform. **In this paper, we developed a framework that further extends the proof strategy of previous works to more general losses, which we think is a significant contribution to the theoretical developments of neural survival models**.  \n\n## Q2: The benefits of frailty and the implications from empirical evaluations\nHere we answer the second and fourth question simultaneously. From an approximation point of view, FN is a more general scheme than PF. However it has been shown in recent developments of deep learning that the universality of some model does not necessarily imply its empirical effectiveness, but instead inductive bias matters, which is sometimes hard to formally quantify. At this stage we would like to argue that the additional inductive bias provided by the PF scheme may sometimes show its usefulness on specific datasets.  \nWhile the precise notion of explaining frailty seems to go beyond the scope of our paper, we conducted heuristic experiments to investigate the **strength of frailty effects** in the real-world data used in our paper. The methodology is to use the bootstrap test [4]. In particular, we construct $R=10$ bootstrap samples of each dataset, and then computes the frailty parameter estimates of each bootstrap replicate and report the means and standard deviations of bootstrap estimates. The test allows standard asymptotics if the estimate itself is asymptotically normal which is to be established in our future study. Therefore we treat the method as being heuristic. The results are summarized in the following table:\n|               |metabric |gbsg |flchain|support|mimic-iii|kkbox    |\n|---------------|---------|-----|-------|-------|---------|---------|\n|bootstrap mean |0.650    |0.569|0.678  |1.945  |0.857    |1.391    |\n|bootstrap stdev|0.017    |0.017|0.070  |0.118  |0.023    |0.066    |\n|mean/stdev     |**38.23**|33.47|9.69   |16.48  |**37.26**|21.08    |\nTable 1: bootstrap tests against the presence of frailty\n\nHere we use (bootstrap mean)/(bootstrap stdev) as an intuitive measure of the presence of frailty effect (the rationale is similar to the Z-score in standard hypothesis testing). According to the table, the effect of frailty is mostly evident in metabric and mimic-iii dataset, which is coherent with the gain of frailty reported in Table 5 in appendix D of the paper.  \n\n## Q3: Numerical evaluations of the synthetic experiments\nFollowing [5], we report the relative integrated mean squared error (RISE) of the estimated survival function against the ground truth and list the results in the following table. The reults suggest that the goodness of fit becomes better with a larger sample size. Moreover, since the true model in the simulation is generated as an PF model, we found PF to perform slightly better than FN, which is reasonable since the inductive bias of PF is more correct in this setup.\n|       |N=1000|N=5000|N=10000|\n|-------|------|------|-------|\n|NFM(PF)|0.0473|0.0145|0.0137 |\n|NFM(FN)|0.0430|0.0184|0.0165 |\n\n[1] Aalen, Odd O. \"Heterogeneity in survival analysis.\" Statistics in medicine 7.11 (1988): 1121-1137.  \n[2]. Zhong, Qixian, Jonas W. Mueller, and Jane-Ling Wang. \"Deep extended hazard models for survival analysis.\" Advances in Neural Information Processing Systems 34 (2021): 15111-15124.  \n[3]. Zhong, Qixian, Jonas Mueller, and Jane-Ling Wang. \"Deep learning for the partially linear Cox model.\" The Annals of Statistics 50.3 (2022): 1348-1375.  \n[4]. Davison, Anthony Christopher, and David Victor Hinkley. Bootstrap methods and their application. No. 1. Cambridge university press, 1997.\n[5]. Zhong, Qixian, Jonas W. Mueller, and Jane-Ling Wang. \"Deep extended hazard models for survival analysis.\" Advances in Neural Information Processing Systems 34 (2021): 15111-15124."}, {"review_id": "FASG4RRiSm", "reviewer": "Reviewer_PAZm", "summary": "The frailty model is one of the popular models in survival analysis, and it is an extension of the classical Cox proportional hazard model.   This paper extends the frailty model by using neural networks, and this paper provides the theoretical analysis of the proposed models.", "questions": "I have no question.", "limitations": "I couldn’t find any description on limitations.", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "+ This paper proposes two new models NFM-PF and NFM-PM by combining the frailty model and neural networks.\n+ This paper provides the theoretical analysis of the proposed models to show its correctness.", "weaknesses": "I think that the (original) frailty model is interesting, simply because it uses slightly weaker assumptions than the classical Cox model.  Even though their difference is small, the frailty model can be significantly better than the Cox model in practical prediction performance.  This is because the Cox model is essentially a linear model and therefore a slight change (i.e., using frailty) can be a huge differentiator.  However, the proposed extensions of frailty models by using neural networks are not so interesting. We already know many neural network models for survival analysis such as DeepSurv [42] and DeepEH [75], which are extensions of the Cox and AFT models.  All of the neural network models for survival analysis are flexible enough, and the advantage of using the concept of frailty in neural network models should be minimal.\n\nThe experiments are not good enough with respect to showing the practical advantage of using frailty in neural network models.  Although the proposed models were compared with many existing models, the performance differences seem coming from the differences of the neural network architectures (e.g., the network structures, the number of layers, the number of neurons in each layer, and many other implementation details) rather than the frailty.  Since the goal of the experiments is to show the advantage of using frailty in the neural network models, the authors should compare the two neural network models with and without frailty by using (almost) the same neural network architecture.\n\nMoreover, the authors should have used the right-censored log-likelihood, which is shown to be a strictly proper scoring rule in [56], as the evaluation metric instead of IBS and IBNLL in the experiments.  The statement “Both IBS and INBLL are proper scoring rules” in Line 301 is not true, and it is shown that both IBS and INBLL are not proper in [56].", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for providing insightful comments. Below we address your specific points:\n## Q1: The motivation of frailty and its advantage\nFirstly, we would like to emphasize that random effect (called frailty), which serves as a principled tool to model *unobserved heterogeneity*, has played an important role in modern survival analysis which is beyond a simple relaxation of proportional hazard assumption. There has been significant and active research concerning the addition of random effect to survival models. The random effect can describe risk or frailty for distinct categories, such as individuals or famlies. There eixst literature on providing theoretical and practical motivation for frailty models by discussing the impact of heterogeneity on analyses. In the influential article [1], Aalen showed that with a sufficiently strong frailty the population relative risk can go from r to 1/r over time, even though the true relative risk stays at r.  As you have pointed out, DeepSurv and DeepEH are flexible models compared to CoxPH. However, DeepSurv still encodes the belief of proportional hazard, and frailty models relax this assumption and is strictly more powerful than DeepSurv.  \n\nSecondly, while the algorithmic extension of incorporating frailty is straightforward, **it is highly nontrivial to establish statistical guarantees**. In this paper the goal is to introduce a rigorous statistical model with provable guarantees, which was lacked in many previous works. As was pointed out in the paper, the current theoretical developments on neural survival models rely on the fact that the loss function is well-controlled by the L2 loss [2, 3], which is not directly applicable to our model due to the flexibility in choosing the frailty transform. **In this paper, we developed a framework that further extends the proof strategy of previous works to more general losses, which we think is a significant contribution to the theoretical developments of neural survival models**.  \n\nFinally, regarding your comments \" the authors should compare the two neural network models with and without frailty by using (almost) the same neural network architecture.\", we have provided such comparison in Table 5 in appendix D. The result shows the benefits of frailty in terms of predictive performance. Next we provide a further study on the presence of frailty:  \nWhile the precise notion of explaining frailty seems to go beyond the scope of our paper, we conducted heuristic experiments to investigate the **strength of frailty effects** in the real-world data used in our paper. The methodology is to use the bootstrap test [4]. In particular, we construct $R=10$ bootstrap samples of each dataset, and then computes the frailty parameter estimates of each bootstrap replicate and report the means and standard deviations of bootstrap estimates. The test allows standard asymptotics if the estimate itself is asymptotically normal which is to be established in our future study. We treat the method as being heuristic. The results are summarized in the following table:\n|               |metabric |gbsg |flchain|support|mimic-iii|kkbox    |\n|---------------|---------|-----|-------|-------|---------|---------|\n|bootstrap mean |0.650    |0.569|0.678  |1.945  |0.857    |1.391    |\n|bootstrap stdev|0.017    |0.017|0.070  |0.118  |0.023    |0.066    |\n|mean/stdev     |**38.23**|33.47|9.69   |16.48  |**37.26**|21.08    |\nTable 1: bootstrap tests against the presence of frailty\n\nHere we use (bootstrap mean)/(bootstrap stdev) as a intuitive measure of the presence of frailty effect (the rationale is similar to the Z-score in standard hypothesis testing). According to the table, the effect of frailty is mostly evident in metabric and mimic-iii dataset, which is coherent with the gain of frailty reported in Table 5 in appendix D of the paper.  \n\nWhile there are many off-the-shelf neural survival models that are flexible and performative, only very few of them come with formal statistical guarantees (We believe some of them actually are theoretically sound without proofs). And the existing proof strategy has its own limitations. In this paper we make contributions to the theoretical developments via adopting new proof strategies to establish statistical guarantees, and we believe the strategy could be used to derive results for other approaches.  \n\n## Q2: On evaluation with proper score metrics\nPlease kindly refer to the first part of the general response for a detailed explanation. First we would like to clarify that the statement in the paper was \"Both IBS and INBLL are proper scoring rules if the censoring times and survival times are independent\". In our opinion, while right-censored log-likelihood (CLL) is a proper metric, it is not well defined for many neural models with implicitly uses the method of nonparametric maximum likelihood (NPMLE) like DeepSurv and Coxtime, and the best thing can do is to evaluate some approximated version, with the approximation itself being not theoretically-principled. For a more thorough evaluation, we provide additional experiments comparing NFM and SuMo-net, which we find is the only baseline method that allows straightforward computation of CLL metric, and NFM is shown to perform slightly better in this metric in comparison to SuMo-net. As shown in Table 1 in the first part of the global response.\n\n[1] Aalen, Odd O. \"Heterogeneity in survival analysis.\" Statistics in medicine 7.11 (1988): 1121-1137.  \n[2]. Zhong, Qixian, Jonas W. Mueller, and Jane-Ling Wang. \"Deep extended hazard models for survival analysis.\" Advances in Neural Information Processing Systems 34 (2021): 15111-15124.  \n[3]. Zhong, Qixian, Jonas Mueller, and Jane-Ling Wang. \"Deep learning for the partially linear Cox model.\" The Annals of Statistics 50.3 (2022): 1348-1375.  \n[4]. Davison, Anthony Christopher, and David Victor Hinkley. Bootstrap methods and their application. No. 1. Cambridge university press, 1997.  "}, {"review_id": "hrUpLvHXh0", "reviewer": "Reviewer_82in", "summary": "The authors propose a neural architecture to estimate the survival function for observations of survival times and censored survival times. The authors describe a methodology to include heterogeneity among the population by including the frailty component. The authors then theoretically and empirically illustrate properties of their approach. The theoretical results are used to provide additional justification for why a neural network architecture is appropriate for this kind of problem. Empirical results are added to place into context the performance of the proposed methods with respect to existing methods. \n", "questions": "How does the empirical convergence compare with the theoretical? \n\nThe authors mention using Clenshaw-Curtis as the integrator, how does the parameters of the integrator (e.g. step size) affect the performance of the method? Is there some numerical bottleneck related to this integration?\n\nIn the theoretical results, the MLPs scale with the number of samples, is there a rigorous connection to the number of parameters in the network and the number of samples? Additionally, was this enforced in the experimental setup, I could not tell if this was the case in the hyper parameter section of the appendix?\n\nDoes the method extend to representing the joint distribution of multiple survival times? E.g. something like $P(T_1 > t_1, \\ldots, T_d > t_d \\mid Z)$? \n\nWas a nonparametric model for the frailty component examined, i.e. could $f_\\theta$ be represented by some generative model constrained to be positive? ", "limitations": "The authors did mention limitations in the appendix. Alongside their notes, it may be good to get a review of where the authors think that their method is appropriate and where other methods would be more appropriate. Additionally, it would be nice to see conditions where the method fails to get a better idea of the properties of the method. This could possibly come from some of the theoretical assumptions being violated (e.g. it’s unclear which ones are really necessary in practice versus for proving the theorems).\n", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper provides a thorough analysis on using neural networks for survival estimation beyond the analyses present in previous works. The authors also consider frailty which they naturally include within their modelling framework, removing some of the homogeneity assumptions of previous works. The empirical results suggest that the method performs well in the considered scenarios. I did not check the details of the proofs, but the theoretical component seems useful and original for predicting convergence rates of the estimator and applicable to other similar estimators. Finally, the writing is clear and easy to follow. ", "weaknesses": "The main weakness involves the empirical evaluation of the paper. None of the empirical gains are major when considering the empirical experiments. This is expected with real world data since it's difficult to estimate the counterfactual (i.e. maybe the survival time could have been much longer/shorter). I would suggest to add an empirical result that compares the methods in idealized settings where the true survival distribution is known and that the metrics can be appropriately compared. \n\nAnother potential weakness involves the theoretical guarantees. If I'm understanding correctly, the theory applies in the case where the networks are well-learned, which is still subject to the issues associated with neural network optimization. It may be good to emphasize this in the text.  ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for providing insightful comments. Below we address your specific points:\n### Q1: Evaluation of idealized settings\nAccording to our understanding, if the true surival time is always observed (i.e., no censoring effects), then the problem boils down to an ordinary regression problem and there are many appropriate metrics and learning objectives that allow a formal study of generalization. In our synthetic experiments, the true survival is known and we provided a pictorial illustration in figure 3 at appendix D that assesses the recovery of the true surival function. We also computed the relative integrated mean squared error (RISE) corresponding to different sample sizes. The results are listed in the following table, showing that the goodness of fit becomes better with a larger sample size.\n|       |N=1000|N=5000|N=10000|\n|-------|------|------|-------|\n|NFM(PF)|0.0473|0.0145|0.0137 |\n|NFM(FN)|0.0430|0.0184|0.0165 |\nTable 1. RISE for the synthetic experiments\n\n### Q2: Theoretical guarantees and optimization issues\nThis is a very good point. Throughout this paper, we rely on the empirical observation that state-of-the-art stochastic optimization algorithms like Adam produce satisfactory results for neural models, and set aside the convergence issue (in the language of optimization). A careful analysis of the optimization landscape of the NFM objective is beyond the scope of the paper, but is a valuable topic and worth future explorations. We will add a discussion in the camera-ready version of the paper.  \nBesides, regarding your comments \"..networks are well-learned..\", according to our understanding, to place ourselves in a machine learning context the issue is somewhat related to some notion of *generalization* in survival analysis, which is to our knowledge an open problem. In some recently proposed frameworks like [6] it might offer opportunities to study formal generalization in survival analysis, which is a promising future direction.\n\n### Q3: Comparison of empirical and theoretical convergence\nIn this paper, we establish statistical guarantees using the Hellinger distance which is hard to empirically evaluate [1]. Therefore we instead compute a intuitive metric RISE and reported in Table 1. We can see from the results that the goodness of fit becomes better with a larger sample size.\n\n### Q4: On the complexity of using Clenshaw-Curtis quadrature\nPlease kindly refer to the second part of the general response for a detailed explanation. We illustrate in Table 2 and Table 3 in the second part of the general response that using $10$ discretization steps in the quadrature method suffices for competitive performance, with the running time comparable to efficient methods like SuMo-net and being much faster than kernel approximation methods like DeepEH. \n\n### Q5: On the scale of the network versus sample size\nIn our theoretical analysis, we construct sieve spaces as a set of MLPs with depth and number of parameters that grow with the sample size at certain speed. This is actually the definition of sieve method [2] that construct a series of parameter spaces that eventually becomes dense in the target function space (the Holder ball defined in (8) in the paper). The construction is mostly driven by approximation-theoretic type arguments [3] and serves as a guide instead of enforcement for empirical hyperparameter choice. As you have pointed out, there are multiple additional factors that may affect the empirical results including optimization issues and training strategy. In our experiments, the final models are selected using early stopping on validation datasets.\n\n### Q6: On extension to multiple survival times\nThis is another very interesting question. Actually the frailty model is even more appropriate for scenarios where multiple surivival times are present, and we may use the powerful tool of frailty to describe the correlation structure among the individual survival times. Such extension requires carefully modeling the dependence structure [5] and will be left to further studies. See also remark 3.1 in the paper.\n\n### Q7: Is it possible for a nonparametric model of frailty?\nThis is yet another interesting and intriguing question. In general, if we allow the frailty transform to be nonparametric, it becomes statistically hard to separate the effect of the frailty transform and its input argument, as they are both infinite-dimensional and extremely expressive (see definition (3) in the paper). Therefore to study such kind of extensions we may need some extra restrictions on the generating process [5]. That said, we did try this option during our empirical evaluations and find this to perform pretty well empirically (comparable but not significantly surpassing parametric NFM). However as they are very subtle in theory, we chose not to report this line of results.\n\n[1] Sreekumar, Sreejith, and Ziv Goldfeld. \"Neural estimation of statistical divergences.\" The Journal of Machine Learning Research 23.1 (2022): 5460-5534.  \n[2] Chen, Xiaohong. \"Large sample sieve estimation of semi-nonparametric models.\" Handbook of econometrics 6 (2007): 5549-5632.  \n[3] Yarotsky, Dmitry. \"Error bounds for approximations with deep ReLU networks.\" Neural Networks 94 (2017): 103-114.  \n[4] Parner, Erik. \"Asymptotic theory for the correlated gamma-frailty model.\" The Annals of Statistics 26.1 (1998): 183-214.  \n[5] Tang, Weijing, et al. \"Survival analysis via ordinary differential equations.\" Journal of the American Statistical Association (2022): 1-16.  \n[6] Han, Xintian, et al. \"Inverse-weighted survival games.\" Advances in neural information processing systems 34 (2021): 2160-2172.\n"}, {"review_id": "rzwOQI2JZ4", "reviewer": "Reviewer_nR17", "summary": "The paper proposes to parametrize frailty models with neural networks and provide results on learning rates for these models.\n\nExperiments demonstrate the efficacy of the proposed models. ", "questions": "1. Everything is evaluated using IBS and IBLL, non-proper scores when censoring is not independent of event-times. Can you prove that the datasets you use indeed exhibit this independence. Can you provide evaluations using the likelihood, a proper score?\n2. The objective seems computationally intensive to compute, can you provide computational comparisons with other methods?", "limitations": "N/A", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. The theoretical analysis is interesting and cool.", "weaknesses": "1. The paper states that: \"[56] used a neural network to approximate the conditional survival function and could be thus viewed as another trivial extension of NHR.\" How is this not exactly the same idea but for frailty models, you're proposing a neural parametrization very much akin to [56]. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for providing insightful comments. Below we address your specific points:\n### Q1: On evaluation with proper score metrics\nPlease kindly refer to the first part of the general response for a detailed explanation. In our opinion, while right-censored log-likelihood (CLL) is a proper metric, it is not well defined for many neural models with implicit uses of the nonparametric maximum likelihood (NPMLE) method like DeepSurv and Coxtime as well as models using pseudo likelihood like DeepEH. To provide more intuitions, think of the Breslow 's estimate for the cumulative hazard used in proportional-hazard type models (CoxPH, DeepSurv, RSF...), the estimate is a step function and therefore there are no corresponding estimates for the hazard function and hinders the evaluation of CLL. Breslow's estimate could be regarded a by-product of the NPMLE construction of partial likelihood. \n \nNow if we still want to evaluate CLL, the best one can do is to evaluate some approximated version as done in [1], with the approximation itself being not theoretically-principled. This is the reason why we did not incorporate CLL into our evaluation in the first place.  \n\nFor a more thorough evaluation, we provide additional experiments comparing NFM and SuMo-net, which we find is the only baseline method that allows straightforward computation of CLL metric, and NFM is shown to perform slightly better in this metric in comparison to SuMo-net. As shown in Table 1 in the second part of the general response. For completeness we paste it below:\n|               | metabric         | gbsg             | flchain          | support          | mimic-iii        | kkbox            |\n|---------------|------------------|------------------|------------------|------------------|------------------|------------------|\n| SuMo-net      | -0.256(0.052)    | 0.367(0.056)     | -1.184(0.023)    | -0.673(0.023)    | -0.052(0.003)    | 0.128(0.004)     |\n| NFM(PF)       | -0.184(0.020)    | 0.378(0.054)     | -1.169(0.022)    | -0.622(0.028)    | **0.125(0.006)** | 0.647(0.003)     |\n| NFM(PF)       | **-0.148(0.027)**| 0.369(0.046)     | -1.166(0.026)    | **-0.588(0.036)**| -0.026(0.001)    | **0.786(0.009)** |\nTable 1: Comparison of NFM with SuMo-net in the CLL metric\n\n## Q2: On the computational complexity of using Clenshaw-Curtis quadrature\nPlease kindly refer to the second part of the general response for a detailed explanation. We illustrate in Table 2 and Table 3 in the second part of the general response that using $10$ discretization steps in the quadrature method suffices for competitive performance, with the running time comparable to efficient methods like SuMo-net and being much faster than kernel approximation methods like DeepEH.  For completeness we paste the results below:\n\n|dataset |          |NFM(PF)   |          |          |NFM(FN)   |          |SuMo-net|DeepEH  | \n|--------|----------|----------|----------|----------|----------|----------|--------|--------|\n|        |steps = 10|steps = 20|steps = 50|steps = 10|steps = 20|steps = 50|        |        |                \n|metabric|9.7s      | 12.5s    |17.7s     |11.4s     | 19.4s    | 36.1s    | 13.58s | 85.6s  |\nTable 2: Running time comparisons of using different number of discretization steps, along with two baselines\n\n|          |           | NFM(PF)   |           |           | NFM(FN)   |           |\n|----------|-----------|-----------|-----------|-----------|-----------|-----------|\n|          |cindex     |ibs        |inbll      |cindex     |ibs        |inbll      |\n|steps = 10|65.16(1.46)|16.28(0.76)|49.02(2.29)|66.82(1.62)|16.03(0.87)|47.96(2.53)|\n|steps = 20|65.16(1.46)|16.28(0.76)|49.02(2.29)|66.63(1.68)|16.07(0.84)|48.04(2.44)|\n|steps = 50|65.13(1.46)|16.28(0.76)|49.03(2.29)|66.79(1.49)|16.10(0.83)|48.10(2.43)|\nTable 3: Performance comparisions of using different number of discretization steps\n\n1] Rindt, David, et al. \"Survival regression with proper scoring rules and monotonic neural networks.\" International Conference on Artificial Intelligence and Statistics. PMLR, 2022."}]

Author rebuttal:
We'd like to thank all the reviewers for providing insightful comments, we will integrate some of the suggestions into the camera-ready version.  
We have found that there are several common issues raised by different reviewers:
- The motivation of frailty and its advantage.
- IBS/IBNLL/Cindex are not proper scoring rules in general.
- The complexity of using Clenshaw-Curtis quadrature to compute the learning objective.

Due to limited space, the first issue is postponed to the response to specific reviewers. Below we provide detailed explanations for the rest issues:

## IBS/IBNLL/Cindex are not proper scoring rules in general.
It was pointed out by several reviewers that IBS and IBNLL are improper scoring rules when the censoring time is not marginally independent with the event time. In [1] the authors advocate using censored log-likelihood (CLL) which is a proper metric. We explain below why CLL is not adopted in our paper, and provide empirical evaluations regarding this metric:
- **CLL might be not well-defined for a lot of survival models**: Many survival models are semiparametric with a nuisance parameter being assumed to lie in an infinite-dimensional functional space. The construction of the renowned partitial likelihood for the CoxPH model is a representative case of nonparametric maximum likelihood (NPMLE) (see [4] for the construction). It is well known in statistical literature that **without proper restriction, nonparametric likelihood can be unbounded** [3, 5]. In the construction of partial likelihood, the estimate of the cumulative hazard function is restricted to be a step function and there are no proper estimate of hazard function under NPMLE, thereby hindering the evaluation of CLL. However, most current survival baselines more or less uses the rationale of partial likelihood, including RSF, DeepSurv, Coxtime and many others. Besides, it is also tricky to apply CLL to models using pseudo likelihood like DeepEH. In [1], the authors used a finite-difference approximation to compute an approximate version of CLL that was not well-defined for models like DeepSurv or Coxtime. Therefore, we think that despite its properness, the applicability of CLL is limited for a fair comparison.
- **A comparison in CLL with SuMo-net**: Although CLL is not a generally applicable metric for survival models as argued above, it is applicable to the NFM framework. Therefore, we provide an additional evaluation of CLL with a comparison to SuMo-net [1], the reults are summarized in table 1 in the attachment pdf. The results suggest that NFM performs on par with SuMo-net, with statistically significantly better results (under $0.05$ significance level) on 4-out-of-6 datasets.   

## The complexity of using Clenshaw-Curtis quadrature to compute the learning objective.
The only hyperparameter involved in the computation of Clenshaw-Curtis quadrature (CCQ) is the number of discretization steps which are common to many numerical integration methods. Here we provide assessments on this hyperparameter on both computational efficiency (measured in wall clock running time) and the trade-off with model performance via comparing the model performance using different number of discretization steps. All experiments are made on the metabric dataset under the optimally tuned architecture, trained for $50$ epochs using a M1 Max MacBook Pro 2021 in its cpu version. We first report the efficiency evaluation in table 2 from the attachment pdf, where we assess three configurations with $\text{steps} \in \{10, 20, 50 \}$ and compare with SuMo-net as well as DeepEH (with a computationally identical configuration of neural function approximator). From the results we find that
- PF is faster than FN and is less affected by increasing the number of integration steps. This stems from the fact that PF decouples the computation of $h$ and $m$, and we only need to evaluate the numerical integral of a one-dimensional function. For the FN scheme, integral of a multi-dimensional function is involved and is thus more time-consuming.
- With a proper choice of steps (i.e., less than $20$), both NFM schemes are comparable in efficiency to SuMo-net, which is regarded as an efficient implementation that involves only a single feedforward neural network call. Besides, NFM are much more efficient then DeepEH, since the computation of DeepEH's objective involves kernel approxiamtions and is quadratic in sample size.
Overall, **the computational cost of both NFM schemes are controllable for a moderate number of steps**. Next we study the effects on model performance. The experiments are conducted on the metabric dataset and summarized in table 3 from the attachment pdf.  We conclude from the table that even using $10$ integration steps suffices for competitive performance, and using a larger number of integration steps seems to have little extra gain. Therefore, we conclude that the additional complexity of using CCQ for evaluation of learning objective is totally affordable, especially when GPUs are deviced for parallel evaluation, the computation time might be further reduced. 


[1] Rindt, David, et al. "Survival regression with proper scoring rules and monotonic neural networks." International Conference on Artificial Intelligence and Statistics. PMLR, 2022.  
[2] Bickel, Peter J., et al. Efficient and adaptive estimation for semiparametric models. Vol. 4. Baltimore: Johns Hopkins University Press, 1993.  
[3] Kiefer, Jack, and Jacob Wolfowitz. "Consistency of the maximum likelihood estimator in the presence of infinitely many incidental parameters." The Annals of Mathematical Statistics (1956): 887-906.  
[4] Murphy, Susan A., and Aad W. Van der Vaart. "On profile likelihood." Journal of the American Statistical Association 95.450 (2000): 449-465.  
[5] Zeng, Donglin, and D. Y. Lin. "Efficient estimation for the accelerated failure time model." Journal of the American Statistical Association 102.480 (2007): 1387-1396.  

Meta-review: The authors provided fairly detailed responses to the reviewers that I think sufficiently address reviewer concerns. Note that of the four reviewers, two leaned toward acceptance. Among the two reviewers who were negative (both voting for "borderline reject"), the main concerns appear to be on experimental results roughly in terms of two aspects: (i) how to know when we would benefit from modeling frailty, and (ii) how the proposed methods compare to non-frailty-based methods. The authors have in their rebuttal as well as their discussions with reviewers addressed these aspects to some extent, where I realize that fully addressing them can be quite challenging (especially for the real datasets considered, I agree with what the authors said in that even for SOTA methods, the performance gaps are often not large). The experimental results are informative even though the proposed methods do not always work the best (they work reasonably well though). That the experimental results (including the newly added ones during the author-reviewer discussion period) are not amazing does not detract from the main strength of this work which is on the theoretical development of modeling frailty with neural nets (including statistical guarantees), which I think is a valuable contribution.
Decision: Accept (poster)