Survival Permanental Processes for Survival Analysis with Time-Varying Covariates
Authors: Hideaki Kim
Venue: NeurIPS 2023 poster (2023)
Keywords: survival analysis
temporal point process
Bayesian estimation
permanental process
representer theorem
kernel method

Abstract: Survival or time-to-event data with time-varying covariates are common in practice, and exploring the non-stationarity in covariates is essential to accurately analyzing the nonlinear dependence of time-to-event outcomes on covariates. Traditional survival analysis methods such as Cox proportional hazards model have been extended to address the time-varying covariates through a counting process formulation, although sophisticated machine learning methods that can accommodate time-varying covariates have been limited. In this paper, we propose a non-parametric Bayesian survival model to analyze the nonlinear dependence of time-to-event outcomes on time-varying covariates. We focus on a computationally feasible Cox process called permanental process, which assumes the square root of hazard function to be generated from a Gaussian process, and tailor it for survival data with time-varying covariates. We verify that the proposed model holds with the representer theorem, a beneficial property for functional analysis, which offers us a fast Bayesian estimation algorithm that scales linearly with the number of observed events without relying on Markov Chain Monte Carlo computation. We evaluate our algorithm on synthetic and real-world data, and show that it achieves comparable predictive accuracy while being tens to hundreds of times faster than state-of-the-art methods.

Survival Permanental Processes for Survival Analysis with Time- Varying Covariates
Hideaki Kim
NTT Human Informatics Laboratories NTT Corporation
hideaki wwikntt..com com
Abstract
Survival or time-to-event data with time-varying covariates are common in prac- tice, and exploring the non-stationarity in covariates is essential to accurately an- alyzing the nonlinear dependence of time-to-event outcomes covariates. Tra- ditional survival analysis methods such as Cox proportional hazards model have been extended to address the time-varying covariates through counting process formulation, although sophisticated machine learning methods that can accom- modate time-varying covariates have been limited. In this paper, we propose non-parametric Bayesian survival model to analyze the nonlinear dependence time-to-event outcomes on time-varying covariates. We focus on a computation- ally feasible Cox process called permanental process, which assumes the square root of hazard function to be generated from Gaussian process, and tailor it for survival data with time-varying covariates. We verify that the proposed model holds with the representer theorem, a beneficial property for functional analysis, which offers fast Bayesian estimation algorithm that scales linearly with the number of observed events without relying on Markov Chain Monte Carlo com- putation. We evaluate our algorithm on synthetic and real-works world data, and show that it achieves comparable predictive accuracy while being tens to hundreds of times faster than state-of-the-arr methods.
Introduction
Survival or time-to-event data analysis has been widely applied for analyzing the dependence of survival time, the time until an event occurs, on covariates. They have a wide ranging list of appli- cations in reliability engineering [[7]], finance [[E]], marketing [20]], and especially, clinical research [HI, 5 28).
An essential part of the survival data analysis is to estimate a hazard function, that is, the instanta- neous probability of events occurring at any particular time, regression function of covariates. Given a hazard function, we can evaluate the impact of covariates on survival times, and assess the degree of risk of experiencing an event in the future, through the survival function, i.e., the prob- ability of no events occurring during a specified interval. The literature contains a vast number of studies that have proposed survival models estimate hazard functions, most of which assume that covariates are time- invarinvariant (e.g., gender and age at the time of diagnosis in clinical research): they range from the classical Cox proportional hazards model (5] to modern machine learning models based generalized boosting machines [[8]] random forests [0Z)], Gaussian processes [[ 15), and deep neural networks [ICC, 22, [9] However, time-varying covariates are common in survival data, and the rapid advances in data collection allow to access long-term and high temporal res- olution covariate data. This is encouraging researchers to explore the non-stationarity in covariates for the accurate estimation/prediction of time-to-event outcomes. An important application of sur-
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

vival analysis with vary-ng covariates is what-if analysis: through simulations of future events under various covariate functions over time, we can find the optimal covariate function or policy that would yield desirable future state. For example, in reliability engineering applications where the events are the failures of machines, time-varying covariates could be the maintenance schedule and the room temperaturehnumidity controlled by air conditioners. The maintenance manager can opti- mize the schedules of maintenance and air conditioning by balancing the risk of failure and the costs of maintenance and air conditioning, If survival models for static covariates were applied to non- stationary data, then the static survival models would fail to estimate the underlying dependence of the hazard function covariates, resulting in unreliable decision decision-making, Therefore, survival mod- els that accommodate time-vary time-varying covariates are needed. This paper assesses survival data with time-varying covariates to estimate a hazard function that takes covariates into account.
The standard survival model the Cox proportional hazards model [5], B. 20I (CoxPH), which forms the logarithm of the hazard function as linear combination of covariates. The original CoxPH as- sumes that covariates are ime-innariant and have constant log-linear effect over time on the hazard function. It was extended to address time-varying covariates through the technique of counting pro- cess formulation [T]. Although the extended CoxPH is very efficient to compute, it fails to address the nonlinear dependence of survival times on covariates, which is often the case in real-world data. To rectify this limitation, two non non- con-aammet machine learning models have been proposed: gen- eralized boosted models [[T]] and random forest- based models [37, ESI. both of which adopt the extended CoxPH as learner/trce, and construct non- con-paamettric predictors in an ensemble manner. They succeeded in addressing the nonlinear dependence of survival times timevaryying covari- ates, but the large number of trees needed for accurate predictions tend to make the algorithm too slow and ineffective for real-time applications.
In this paper, we propose a novel Bayesian survival model to estimate a hazard function as a nonlin- ear regression function of covariates from survival data with time-varying covariates. To construct a scalable algorithm, we utilize permanental process, a doubly-stochastic point process which as- sumes the square root of hazard function to be generated from a Gaussian process; its computational advantages have been highlighted by recent studies [ I& 24 60. 95].. We tailor it for survival data in the counting process format, and show that the tailored process, which call the survival perma- nental process, holds the representer theorem [A4]] the maximum a posteriori estimator of the latent function can be represented as finite linear combination of a transformed kernel product evaluated on the event time points. The representer theorem offers us a feasible Bayesian estimation algorithm that scales linearly with the number of events observed without relying on Markov Chain Monte Carlo computation. Furthermore, we derive the predictive distribution and the marginal likelihood in one feasible form, enabling us implement the uncertainty evaluation and the hyper-parameter optimization in a Bayesian manner. To the best of our knowledge, it is the first time to exploit the representer theorem for survival data analysis. We provide python codes to reproduce the results in this paper"
In Section 2 we introduce the survival permanental process (SurvPP) and construct a scalable Bayesian estimation algorithm for it. In Section E. we outline related work. In Section EL we compare SurvPP with reference survival models synthetic and real-world data, and confirm that our algorithm achieves comparable predictive accuracy while being tens to hundreds of times faster than -tate-of-the-art survival models. Finally, Section E states our conclusions.
2 Methods
2.1 Survival Permanental Processes
We assume that right-censored dataset with time-varying covariates, D (t) Tu. nu) Pu))== is observed, where ya(t)= : (0, T.] + y Rdd T, E RR and 1) E {0, 1} are the map of d- dimensional covariate, the end time of observation, and the indicator that represents whether an individual experienced an event (7), = 1) or was right-censored (n, 0) at T., for individual u {1, 2, ... U'. respectively. We consider tailoring a permanental process for the survival data, where latent function, 2z)): y + I, is generated from Gaussian process (GP) on covariate space y and an event time for each individual u is generated from a point process with hazard function, X (t),
"Code and data are provided at https /tte::///ttttt com/H...Kim/SurwPPP
2

---

such that
((D(=ly) (P(x(y))
p(z(y)]D)
(1)
(r(yy) P (D)x(y)) GP(x(y))k)
1f0 Th >m (f)dtt Aa(t)= x2 (y,(t))
log logp(D(z(()) = Σ nu log Au (TT)
(2)
where (P(x(y)) [k) represents GP with kernel function k(y. ((D) (())((())))))))))))))))))) (()) is the likelihood func- tion of the point process, while (x(yy) in the denominator represents the integral over the function or the infinie-dimmesionall variable xzy). We call the model defined by (I-2) the survival perma- nental process (SurvPP). Hazard function Au (t) represents an instantancous probability of events occurring at each point in time, and our goal is to estimate the functional form of the hazard func- tion over covariate domain, x² (3)) from the survival data D.
Given a hazard function over covariate domain, x² (y), and covariate map for individual (y),), we can evaluate the survival function and the probability density distribution of event time at arbi- trary time point 1, denoted by S(t) and Pa (t), respectively, as follows:
S(t) exp /. (33(((()))
(t(t) = (t)) exp (x( 22 (')(dss))
By using (E), we can perform survival analyses that include the analysis of the expected duration of event time, and the assessment of the degree of risk that individual will experience the event before a specified period.
2.2 Counting Process Format of Data
Traditionally, survival data with time-varying covariates have taken the counting process format [[]]
The format assumes that, for each individual u, covariates measured Jn finite representative time points and can be regarded as constant between successive time points,
... Ju where (s.. sy--) = (0, T.). then splits individual u's record, (t), Tu, nu), into Ju pseudo- individual records in
(l) yd. I LE (s/, s
{ $ E/. o =
na- I(j = Ju -1),
(5)
where I(.) represents the indicator. Joining the pseudo-individual records of U individuals together results in the counting process format of data:
D = (TT Tj. Eg. D
Ju
(6)
In this paper, we employ the counting process format of (E), and rewrite the likelihood function of SurvPP (2) as
(Dgr(D)z()))
Eg log x (v,)
-TT -TT--T))=² Σ log ( -
ΣS (7)
where A, and {#n J represent the durations of observation and the covariates observed at event times, respectively:
V
A, = T 17,
( - N
&
The
(
J=1
Note that the number of events, N, is usually much smaller than the number of pscudo-individual records: N < J.
In practice, there are two ways of defining the representative time points { one is as the points that were observed, which might be sparse over time due to measurement constraints; the other is as the denser points obtained by using interpolation methods (e.g., splines). This paper assumes the latter, and thus various values of Jn were considered in the synthetic data experiment (see Section H) Although the counting process format of input assumes piecewise stationarity in covariates, this assumption not so strong because we can adopt the representative time points in any density if necessary.
3

---

2.3 Maximum A Posteriori Estimator
We consider the problem of obtaining the maximum a posteriori (MAP) estimator of xly), denoted by kyy), to maximize the posterior probability (II). We derive f(y) through the approach of using the path integral representation of GP [23]].
1 |K]|
[- I//12xa k" (v. y ())rr())(())dddl Dx(y),
(P(z(y)(k) Dx(y)
V
exp
(9)
where K. and J (y, y dy' are the integral operator with kernel function k(y, y/ and its inverse operator, respectively: x(y) = Jy k0*) (y.y') )')a(/))dd'''
(y.y)) = 6(y - y²). Here [|] represents the function determinant [ISS] of K. defined by the product of its eigenvalues [23]]. Using the representation of (H), we write the posterior of SurvPP (II)2) in the following functional form,
p(r(y)D))xx=
1 P(D) exp - -(((())).(())
log Ki] Øx.
(10)
where 5(x(y),z(y) is the action integral, defined by
5(x(y),r(u)) {
Σ x² (y)), yz) - Σ log (r()))o - (m) dy, (11)
and x(y) In (2, 3 (yy''dy'' Then we apply calculus of variations to the action integral, where the functional derivative of S(x (y). on the MAP estimator ((y) should be equal to zero: 69 ,2(1, + (ā(y] x (ily) O, which results in the exact MAP estimator.
68
N
a (ty) 2 hly Un)un, Un = =(ūū))
(12)
where hyy, (() is a transformed kernel function that solves discretized Fredholm integral equation of the second kind [36]].
My. + Σ k(y (,,,)),,,(y, v') k(y.
(13)
j=1
See Appendix A for the detailed derivations. Equation (122) shows that the MAP estimator of SurvPP involves the representer theorem under the transformed kernel function h(y, y'), and thus the Bayesian estimation reduces finite-dimensional optimization problem. Here, we call h y') the equivalent kernel following studies by Flaxman et al. [I2]) Walder & Bishop [B5],, and Kim et al. [24]]. Given the equivalent kernel, the unknown coefficients v, in (l2) solve the simultaneous quadratic equations derived from (I2)..
A U
Σ h(Tn: Un.) 1= 0, n= 121,2.... N. ===1
(14)
In this paper, we estimate set of coefficients, {vn, )N by solving minimization problem of the mean of the squared residuals, Σn rr,l/N, with popular gradient descent algorithm, Adam [26]-
2.4 Equivalent Kernels
The equivalent kernel My. y solves the discretized Fredholm integral equation (.) When the kernel function of GP has degenerate form with rank M(< oo) such that
M
k(g.y))
E
On (1) (3/) (b (dyy",
(15)
m=1

---

it is easily shown that the discretized Fredholm integral equation P 24 as,
can be solved analyticllly
My. = (b4))(())))))) (TM (TM
2 '(y'), A = Σ A (,
(16)
where IM the M x MM identity matrix, and ( = (+) (y): (y u (v)) Note that A defined by sum of outer products, which can be rewritten through matrix- multiplication as
ABB' B -V[Vii(())).... V A
(17)
Empirically, implementation by matrix-man -matrix multiplication of outer products when J > 1.
is substantially faster than the sum
When kyy, y²) has degenerate form with rank M, the relation (LG shows that the equivalent kernel hy yo) also has degenerate form obtained through Cholesky decomposition:
hyy = (Lo(v)) (ooov)))(LL(())) (Lo(y')) L = (IM + 2A)
(18)
The degenerate equivalent kernel (LE)) offers fast Bayesian estimation that scales linearly with N of a Gaussian kernel
(see Section 2.G). In this paper, we used the random feature map [67] 42 33
to obtain degenerate form of kernel (M=500).
If the covariate map, yu,(t,, (t), can assumed to be smooth enough over time, a more accurate evalu- ation of A than equation (ITG) is possible (see Appendix D.I., but this was not exploited in the main experiments.
2.5 Predictive Distribution and Marginal Likelihood
One of the advantages of GP models over non- non-Bayead nan-Bayeiian approaches is that they can provide pre- dictive distributions and marginal likelihoods, which enable us to perform uncertainty evaluations, hyper-parameter optimization, and model selection in Bayesian manner. Following the methodology with the path integral representation GP [23] 244), SurvPP (I-2) adopts a Laplace approximation in the functional space, and finds the approximate form of the predictive distribution and the marginal likelihood. We only show the results due to space limitations. For details, see Appendix
The marginal likelihood, p(D), is given as
1 log IN Z H
1 log [IM + 2A + (log - 1)N,
P(D) log Z
(19)
2
2
where IN is the identity matrix with size N, and
nnn' 'ymn', Han' h(gn, In)), h(y) (h(y. Ø1)..... ..... h TN
(20)
The predictive distribution of intensity function on a covariate value, p() = (z))) is given as
21 to
(8² + 0
(A)(()) = T(v) X exp((AAp))) μ == 20
V/ 20 +a)
(21)
2 +
where & is the abbreviation of f(y). and the predictive variance, o, is defined by h(y, y) h(y) (Z- HH - (LII) represents the predictive distribution of the hazard function is worth noting that the posterior process of SurvPP is given
regression function of covariates. by permanental process, and we can evaluate the distribution of the survival function and perform a risk-aware survival analysis by generating random samples of the estimated hazard function of covariates, A(y).
2.6 Computational Complexity
The computational complexity of evaluating the equivalent kernel (LG for covariate pair (y. y is 0(M3 +dJM+MM where O(dJM) and (JJM) come from the computation of feature maps, on (yy), and the sum of outer products (I7)), respectively.
When the equivalent kernel is given in degenerate form with rank M (< N) such that hly, y') = M
D pm (v) (m (v'), the objective function to be minimized in MAP estimation, EN Ir, 1r,11//0 I2/N,
5

---

incurs a linear computation with the number of observed events N, that is, O(NM), for each evalu- ation in gradient descent algorithms: The vector of residual, r = (TI,...,'NN) T can be expressed as [R RR v)) @ (2v) 1, where R is the N matrix defined by Rum = Pm (. and O represents the Hadamard product. R costs (d)
Given an equivalent kernel with degenerate form, the evaluation of the predictive variance, o' 1y.. h (Z + H needs the computation of O(M + N 2M)) NxN matrix H can be decomposed into a product of N x M matrix R and its transpose as H = RR The matrix inversion is transformed as (Z + RR ) Z Z Z-IR((I + R Z R) -R) Z 11 through the Woodbury matrix identity, which costs O(M + NM2).. Note that Z- R and R Z-! cost M( because Z is a diagonal matrix.
In computing the marginal likelihood (II), the complexities of the first, the second, and the third terms are O(N), (MS + N 2 and O(M* + JM,,, respectively, where the matrix determinant lemma is used in computing the second term.
In total, the computational complexity of SurvPP is O(N + (N J) M)M+ M2). where Q is the number of gradient descent iterations. This feasible computation, which scales linearly with the number of observed events, N, the data size, J, and the covariate dimensionality, d. is achieved by exploiting the representer theorem. It should be emphasized here that data size J is not part of the gradient descent iteration term. This is clear advantage over conventional survival models which naively incur ((dN + dJ + JN))) computation: Typically, J is larger than N in time-varying covariate scenario (see the counting process format in Section 22)), and the discrepancy between and N becomes more substantial when the measurements of covariates are made more frequently or/and over a longer period of time.
3 Related Work
Survival Models for Time-invariant Covariates: The most popular survival model is the Cox proportional hazards model [5]] (CoxPH), which is a semi-parametric model that forms the hazard function as
(() h(() exp(31 + 9282 + +...,
(22)
where the base hazard function, f((t. is obtained by the non-parametric Breslow/Fleming- Harrington estimator [G] [I3] and the log-linear regression coefficient ($,....) is estimated by maximizing the partial likelihood function. CoxPH is very efficient to compute and scales linearly with the number of events [3%], but cannot address the nonlinear dependence of survival times on covariates. To overcome this limitation, a vast number of survival models have been proposed that replace the log-linear parametric function (22) with non-linear one, such as generalized boosted models []8], random survival forests [17]], Gaussian process models [[]]. [2]], and deep neural net-
work models | 22 39
For comprehensive review, see Wang et al. [B6]]]
Survival Models for Time-varying Covariates: The counting process format of input [M] plays a central role in extending static survival models into those suitable for time-varying covariates. For each individual, the counting process format splits her/his observation period into multiple short intervals, assigns a constant value of covariate to each interval, and marks the intervals as being censored when no event present, resulting a set of pseudo-individus
right/left-censonee servations with a constant covariate (see Section 272). CoxPH was extended to accommodate the counting process format, and the extended CoxPH can estimate the hazard function based on sur- vival data with time-varying covariates [[]]. To alleviate the simplest assumption in the extended CoxPH, generalized boosted models [[6]] and random survival forests [27 [N]. both of which are ensemble approaches, have been also extended to accommodate the counting process format, where the extended CoxPH is employed as a weak learner or a tree. While they can address the nonlinear dependence of survival times time-varying covariates, the large number of learners/trees needed for accurate predictions could make the algorithms too slow and ineffective for real-time applica- tions. Also, Cygu et al. reported that random survival forests for time-varying covariates required substantial amounts of computer memory for large datasets [[]].
ob-
Exogeneity and Endogeneity of Covariates: When considering survival analysis in the presence of time-varying covariates, we need to distinguish between exogenous and endogenous covariates. Kalbfleisch and Prentice [[T] define an endogenous (internal) covariate as the output of a stochastic

---

A
TLL: Jo =
TLL: Ju =0
B
CPU [sec]
IRFM
GBM
0.2
GBM
0.0
10
0.4
0.6
10P
RFM
AUC:Ju
AUC.J,
SurvPP
0.8
10
0.7
10
CoxPH
0.3
0.7 Evaluation Point
0.5
0.9
0.3
0.7 Evaluation Point
0.5
0.9
1.5 10
20
50
Ju
Figure 1: Performance on the dataset of log-linear hazard function Atsn (A) Box plot of TLL and AUC as functions of evaluation point: the higher, the better. (B) The CPU times demanded for estimating hazard function. The error bars represent the standard deviations. For GBM and SurvPP, the average epu times over 9-points grid search of the hyperparameter are displayed.
process that is generated by the individual under study. In contrast, exogenous (external) variate is not influenced by the individual under study. Diggle et al. [ suggest using similar but slightly different definitions. The primary purpose of survival analysis with exogenous covariates is to estimate a hazard function as an explicit function of covariates, and to make predictions of fail- times under various possible future covariate functions (i.e., what-if analysis). In contrast, the primary purpose of survival analysis with endogenous covariates is to make predictions of failure times by using the past observations of covariates, where joint modeling approaches with recurrent neural networks, which jointly model the stochastic process of covariates and failure times, have been developed intensively [29,, [3]] In this paper, we consider exogenous covariates, and survival models for exogenous covariates cannot be compared directly with those for endogenous covariates because the tasks are different. Note that a joint modeling approach for exogenous time-varying covariates has been proposed very recently[32], which was not included in the benchmark models because was published a month after this paper's submission.
Permanental Process:
The permanental process is a variant of Gaussian Cox process that assumes Its computational
the square root of hazard function be generated from Gaussian process [51]]
advantages have recently been highlighted in machine learning research [T2
I8
24
60, ESS] In
particular, the representer theorem in the permanental process has been exploited through the RKHS
theory [[2], the Mercer's 's theorem [45], and the path integral formulation [23, 24j]. Although the key derivation of SurvPP is based on the path integral methodology used for augmented permanental process (APP) [244] SurvPP is non con-tiivial extension of APP: (i) SurvPP can accommodate multiple trials of event sequence data, while APP assumes one trial of evet sequence;
(ii) in SurvPP, the end time of observation is stochastic variable that depends on the time of event occurrence, while APP assumes that the end time of observation is given. We discovered, for the first time, that the representer theorem holds for such complicated point process or a survival model.
4 Experiments
We examined the validity of our proposed model by comparing it with conventional survival mod- els on synthetic and real-world data. As benchmark models, we adopted Cox proportional haz- ards model (CoxPH), generalized boosted model (GBM), and random forest-based model (RFM): We implemented CoxPH and GBM with the established algorithms provided in the packages survival. coxph [II]] and gbm3. gbmt [IG], respectively; We implemented RFM by using class ltrcrrf in the open R code provided by Yao et al. [[S]. The benchmark models can accommodate time-varying covariates via counting process format of input (Surv (t0, t1, event) in survival). For our proposal, we implemented SurvPP by using TensorFlow-22
A MacBook Pro with 12- core CPU (Apple M2 Max) was used, where GPU set off (tf device (* (c/cpu:0)) for

---

A
TLL: Ju =
TLL: Ju =0
B
CPU [sec]
0.2
CoPH GBM
RFM
0.0
SurPPP
GBM
101
10P
-0.6
RFM
A
AUC:Jy 20
1.0
0.8
10"
0.4
10
CoxPH
0.3
0.7 Evaluation Point
0.5
0.9
0.3
0.7 Evaluation Point
0.5
0.9
1.5 10
20
50
Ju
Figure 2: Performance on the dataset of non-linear hazard function Anon (t). (A) Box plot of TLL and AUC as functions of evaluation point: the higher, the better. (B) The CPU times demanded for estimating hazard function. The error bars represent the standard deviations. For GBM and SurvPP, the average epu times over 9-point grid search of the hyperparameter are displayed.
fair benchmark comparisons. For GBM and SurvPP, the hyper-paramen parameters optimized through 9 9-point grid search: the number of trees and the shrinkage for GB.M, and the kernel parameters for SurvPP. For details of the model configurations, see Appendix Q
4.1 Synthetic Data
We created two survival datasets with vire-narying 2-D covariate. One was generated from a hazard function with log-linear dependence covariates, and the other with nonlinear dependence covariates as follows:
Auto (t) = h(t) exp 1.5y1 (t) (t), Anon(t) exp 2 5(y + v²())
(23)
where htt) is the base hazard function defined by Weibull hazard function, h(t) 2 13/2 For each of the datasets, we considered U individuals, each of which had 2-D covariate function of time:
(t)
(a, cwscooossaaac + xY)), a', cos(2 + mTaa), u= u= U, (24)
where an((a(s) and 7u()()) were sampled uniformly over [0, 1], Wa was sampled uniformly over [5, 10]. and w/ was sampled uniformly over [20, 30). We set the censoring time as for all individu- als, and event time generated from denoted by tx, was censored (77, when was over 1: T. = min(t... 1).
The predictive performances were evaluated based on the dynamic area under the ROC curve (AUC) [34] and the test log-likelihood (TLL),
ΣV w1= I(T ()T(TT (S(b,(((())))))) ((1l)w
AUC(t)
(25)
(E I(T (E I(T, f)w.)
TLL(() = T Σ 1(1 Tw) na lo..oooglllllloooo...... log (minn(,, Tw) yn ()){
(26)
where W. is the inverse probability of censoring weight, and (()) and S(tlyu())) are the estimated hazard and survival functions given covariate map yu(-), respectively. The dynamic AUC essentially estimates the C-index at each time, and is commonly used the literature of survival analysis with varying covariates.
The evaluation point, 0 < la < 1, was set as E {0.3, 0.5, 0.7, 0.9). For each dataset, we randomly split the U individuals into 10 subgroups, assigned one to test and the others to training data, and

---

A: PBC
TLL
AUC
CPU [sed]
1.0
10
-2
0.9
-4
10
-6-
0.8
10-1
-8
CoxPH GBM
RFM
SurvPP
0.7
3
4
3
x10²
CoxPH GBM
RFM SurvPP
Evaluation Time Point
Evaluation Time Point
B: SANAD
[sec]
1.0
10²
0.8
10"
0.6
0.4
17
18
19 Point
20
x10²
17
18
19 Point
20
x10?
CoxPH GBM
RFM SurVPP
Evaluation
Figure 3: Performances on real-world datasets. TLL (the higher, the better), AUC (the higher, the better), and CPU times (the lower, the better) on PBC dataset (A) and SANAD dataset (B).
conducted 10-fold cross evaluation the predictive performances. We reformatted the training data into counting process format by equally splitting individual u's observation periods into Ja sub- region and assuming the covariates to be constant within sub-region; the survival models were applied to the training data. In this experiment, considered U = 10x and Ju E {1, 5, 10, 20, 50]: the number of observed events N was 809 for Atin and 818 for Anon respectively.
It should be noted here that when time-varying covariates are considered, the time duration from the entry, t, is itself time-varying covariate. SurvPP considered a 3-D covariate map, (t) (e' (t), y² (t), t), and performed hazard function estimation, while the other models assumed the hazard function to be (t))), and (() and f(y were estimated separately.
Figure MI displays the predictive performance on the dataset of Atiz (t). Figure IIA shows that, except for RFM, there were significant differences in the performances between the compared models, but CoxPH slightly outperformed the other models in TLL at some evaluation periods. The result is plausible because the underlying generative process was perfectly consistent with the log-linear and proportional assumptions of CoxPH. The comparison in TLL between Ju = and J. = 20 suggests that more frequently measured covariates would yield better estimations/prediicionn in scenarios of time-varying covariates. Therefore, the feasible computational complexity of SurvPP regarding Ju a clear advantage over the reference models, which is confirmed by Figure DIB
Figure 2 displays the predictive performance on the dataset of Anon (t), showing that the non- parametric models, SurvPP and GBM, significantly outperformed CoxPH, while the performance gaps between SuryPP and GBM were marginal and no significance. However, Figure ZB shows that SurvPP could estimate the hazard function hundreds of times faster than GBM when covariates were measured frequently (J,, = 50), and was even faster than the simplest CoxPH. The preferable scalability of SurvPP's algorithm regarding Ju comes from the fact that J, (or J=Σ., J,,) is not part of the iterative optimization (Section 2.G). which is a fruit of SurvPP's representer theorem. RFM was comparable with SurvPP and GBM in AUC, but was worse than either of them in TLL. RFM might need more careful tuning of the
but this was not fully investigated this time.
We conducted additional experiments independent validation datasets (see Appendix (2.2), and on larger synthetic datasets (U < 110)) to examine the model's 's computation scalability regarding the event number N (see Appendix D.
Real-world Data
We examined the validity of SurvPP against the benchmark models on two real-world survival data sets, Mayo Clinic Primary Biliary Cholangitis data (PBC) and Standard And New Antiepileptic

---

Drugs study data (SANAD), provided by R packages survival (LGPL-3) [ and joineR (GPL- 3) |لاا| respectively. In PBC, 312 patients with primary biliary cirrhosis were enrolled in a random- ized medical trial at the Mayo Clinic between 1974 and 1984 [8], where events were the time of death; one static and eleven time-varying covariates were measured at entry and at yearly intervals, which include age at entry, alkaline phosphotase. logarithm of serum albumin, presence of ascites, aspartate aminotransferasee, logarithm of serum bilirubin, serum cholesterol, condition of edema, presence of hepatomegaly or enlarged liver, platelet count, logarithm of prothrombin time, and pres- ence or absence of spiders. SANAD was an unblind randomized trial that recruited patients with epilepsy for whom carbamazepine (CBZ) was considered to be standard treatment and they were randomized to CBZ or the newer drug lamotrigine (LTG), where 605 patients were included and event was the time to treatment failure; adopted calibrated dose time-varying covariate, and three static covariates including age of patient at randomization, gender, and randomized treatment (CBZ. or LTT);; calibrated dose was measured at 166-day intervals on average, and used linear interpolation to obtain the values of calibrated dose at which the measurement intervals were reg- ularly trisected. We randomly split the individuals into 10 subgroups, assigned one to test and the others to training data, and conducted 10-fold cross evaluation of the predictive performances.
Figure DA displays the predictive performance on PBC. It shows that CoxPH achieved very high AUCs (> 0.9), suggesting that the underlying generative process could be consistent with CoxPH's simple assumption of log- linearity and proportionality. Thus the semi-parametric model, CoxPH, is likely to achieve equal or slightly better performances than the nonparametric models, which is consistent with the result. Figure BB plots the predictive performance SANAD, where the not so high AUCs of CoxPH suggests that the underlying dependence of intensity on covariates is nonlinear. The figure shows that SurvPP achieved better performance than CoxPH, and achieved comparable performance while being substantially faster than GBM.
5 Conclusions
We have proposed con-parametric Bayesian survival model to address survival data with time- varying covariates. We tailored a permanental process such that the latent hazard function is defined on covariate space and right-censored observations in a counting process format can be handled, which we call the survival permanental process (SurvPP). Through the path integral formulation of Gaussian process, we showed that SurvPP encompasses representer theorem, and derived a feasible estimation algorithm that scales linearly with the number of observed events. We evaluated SurvPP on synthetic data, confirming that it achieved comparable predictive accuracy while being tens to hundreds of times faster than state-of-the-art methods.
Limitations & future work: We examined SurvPP for relatively low-dimensional covariates, and its potential suitability for high-dimensional covariates remains be clarified. Because SurvPP is based on normal Gaussian process, the (cquivalent) kernel function might be too simple to discover meaningful representations in nigh-dimensional covariate data. A promising direction is to apply deep kernel learning to SurvPP, where high-dimensional covariates are transformed by non- linear mapping with a deep architecture. A technical issue that we could not search an appropriate set of (a lot of) parameters of neural networks, because our proposed scheme performs the hyper- parameter optimization by grid search, not by gradient descent. Also, Gaussian kernels, which were used in the paper, naively require scale parameter for each dimension of data, resulting in a high dimensional kernel parameter for high-dimensional covaritate scenarios. The grid search becomes prohibitively costly with high dimensional parameters, but we addressed the problem by a well- known approach that normalizes (c.g., centering and scaling) each dimension of the data and puts a common scale parameter across all dimensions of normalized data. This approach empirically works robustly, but a more sophisticated approach could improve the performance of SurvPP. Variational Bayesian approximations with inducing points might address the technical issues of hhper- parameer optimization, and thus is an important next step in our study.
As in ordinary permanental processes, the nodal line problem could arise in SurvPP: the posterior distribution the latent variable zr(-) has many local modes since =-((.) can lead to similar hazard functions X =-)) and artificial zero crossings of zz-) could happen, especially on locations where
the hazard function low. John and Hensman [[S] have proposed to extend the quadratic link
function to include an offset parameter 8, so that A = (z(y) + 3) which is valid for SurvPP.
10

---

References
[1] Per Kragh Andersen and Richard D. Gill. Cox's 's regression model for counting processes: A large sample study. The Annals of Statistics, pages 1100- 120, 1982.
[2] Kendall Atkinson. A personal perspective on the history of the numerical analysis of Fredholm integral equations of the second kind. In The Birth of Numerical Analysis, pages 53-72. World Scientific, 2010.
[3] Norman E. Breslow. Contribution to discussion of paper by DR Cox. Journal of the Royal Statistical Society, Series B, 34:216-217, 1972.
[4] Taane G. Clark, Michael Bradburn, Sharon B. Love, and Douglas G. Altman. Survival Analysis Part I: Basic concepts and first analyses. British Journal of Cancer, 89(2):232-238, 2003.
[5] David R. Cox. Regression models and life-tables. Journal of the Royal Statistical Society: Series B (Methodological), 34(2):177-22,, 1972.
[6] Steve Cygu, Jonathan Dushoff, and Benjamin M. Bolker. pcoxtime: Penalized Cox propor- tional hazard model for time-dependent covariates. arXiv preprint arXiv:2102.0229) 2021.
[7] Steve Cygu, Hsien Seow, Jonathan Dushoff, and Benjamin M. Bolker. Comparing machine learning approaches to incorporate time-varying covariates in predicting cancer survival time. Scientific Reports, 13(1):1370, 2023.
[8] E. Rolland Dickson, Patricia M. Grambsch, Thomas R. Fleming, Lloyd D. Fisher, and Alice Langworthy. Prognosis in primary biliary cirrhosis: Model for decision making, Hepatology, 10(1):1-7, 1989.
[9] Peter J. Diggle. Analysis of Longitudinal Data. Oxford University Press, 2002.
[10] David Faraggi and Richard Simon. A neural network model for survival data. Statistics in Medicine, 14(1):73-82, 1995.
[11] Tamara Fernández, Nicolás Rivera, and Yee Whye Teh. Gaussian processes for survival anal- ysis. In Advances in Neural Information Processing Systems 29, 2016.
[12] Seth Flaxman, Yee Whye Teh, and Dino Sejdinovic. Poisson intensity estimation with repro- ducing kernels. In Artificial Intelligence and Statistics, pages 270-279. PMLR, 2017.
[13] Thomas R. Fleming and David P. Harrington. Nonparametric estimation of the survival dis- tribution in censored data. Communications in PattisticsTheony and Methods, 13(20):2469- 2486, 1984.
[14] Adrian Gepp and Kuldeep Kumar, The role of survival analysis in financial distress prediction. International Research Journal of Finance and Economics, 16(16):13-3 34, 2008
[15] Israel Gohberg, Seymour Goldberg, and Nahum Krupnik. Traces and Determinants of Linear Operators, volume 16. Birkhäuser, 2012.
[16] James Hickey and Greg Ridgeway. Generalized Boosted Regression Models, 2023. URL https://gt thub. com/ www.ccooooooccecccccccceceeee www-ddveeoppereeeeoc gmb33
[17) Hemant Ishwaran, Udaya B. Kogalur, Eugene H. Blackstone, and Michael S. Lauer. Random survival forests. The Annals of Applied Statistics, 2(3):841 2008
[18) S. T. John and James Hensman. Large-scale Cox process inference using variational Fourier features. In International Conference on Machine Learning, volume 80, pages 2362-2370. PMLR, 2018.
[19] Alistair E.W. Johnson, Tom J. Pollard, Lu Shen, Li-wei H. Lehman, Mengling Feng, Moham- mad Ghassemi, Benjamin Moody, Peter Szolovits, Leo Anthony Celi, and Roger G. Mark. MIMIC-III, a freely accessible critical care database. Scientific Data, 3(1):1-9, 2016.
11

---

[20] Euy- Suy-Young Jung, Chulwoo Back, and Jeong-Dong Lee. Product survival analysis for the app store. Marketing Letters, 23:929- 9419 2012.
[21]
John D. Kalbfleisch and Ross L. Prentice. The Statistical Analysis of Failure Time Data. John Wiley & Sons, 2011.
[22] Jared L. Katzman, Uri Shaham, Alexander Cloninger, Jonathan Bates, Tingting Jiang, and Yu- val Kluger. DeepSurv: personalized treatment recommender system using Cox proportional hazards deep neural network. BMC Medical Research Methodology, 18(1):1-12, 2018.
[23] Hideaki Kim. Fast Bayesian inference for Gaussian Cox processes via path integral formula- tion. In Advances in Neural Information Prucessing Systems 34, 2021.
[24] Hideaki Kim, Taichi Asami, and Hiroyuki Toda. Fast Bayesian estimation of point process intensity function of covariates. In Advances Neural Information Processing Systems 35, 2022.
[25] Minyoung Kim and Vladimir Pavlovic. Variational inference for Gaussian process models for survival analysis. In Uncertainty in Artificial Intelligence, pages 435 445, 2018.
[26] Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint 1rXiv: 2014.
[27] Chin Diew Lai and Min Xie. Stochastic Ageing and Dependence for Reliability, Springer Science Business Media, 2006.
[28] András Lánczky and Balázs Gyōrffy. Web-based survival analysis tool tailored for medical research (kmplot): development and implementation. Journal of Medical Internet Research, 23(7))e27633,, 2021.
[29] Changhee Lee, Jinsung Yoon, and Mihaela Van Der Schaar. Dynamic-DeepHiit A deep learn- ing approach for dynamic survival analysis with competing risks based on longitudinal data. IEEE Transactions on Biomedical Engineering, 67(1):122-13) 133, 2019.
[30] Chris Lloyd, Tom Gunter, Michael Osborne, and Stephen Roberts. Variational inference for Gaussian process modulated Poisson processes. In International Conference on Machine Learning, volume 37, pages 1814-1822. PMLR, 2015.
[31] Peter McCullagh and Jesper Moller. The permanental process. Advances in Applied Probabil- iry, 38(4))873-888 2006.
[32] Victor Medina- Medina-Olivares, Raffaella Calabrese, Jonathan Crook, and Finn Lindgren. Joint mod- els for longitudinal and discrete survival data in credit scoring. European Journal of Opera- tional Research, 107731145511444, 2023.
[33] Chirag Nagpal, Vincent Jeanselme, and Artur Dubrawski. Deep parametric time-to-event re- gression with varying covariates. In Survival Frediction-Algoriihen Challenges and Ap- plications, pages 184-193. PMLR, 2021.
[34] Norberto Pantora- Pantoja-Galicia, Olivia I. Okereke, Deborah Blacker, and Rebecca A. Betensky. Con- cordance measures and time-dependent ROC methods. Biostatistics & Epidemiology, 5(2): 232-249. 2021.
[35] Pete Philipson, Ines Sousa, Peter J. Diggle, Paula Williamson, Ruwanthi Kolamunnage-Dona, Robin Henderson, and Graeme L. Hickey, joineR: Joint Modelling of Repeated Measurements and Time-to-Eveny Data, 2018. URL https: /ettp:///itthbb com/ graemeeleehicctee//oonieeit package version 1.2.8.
[36] Andrei D. Polyanin and Alexander V. Manzhirov. Handbook of Integral Equations. CRC press, 1998.
[37] Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In Ad- vances in Neural Information Processing Systems 20, 2007.
12

---

[38] Greg Ridgeway. Generalized Boosted Models: A guide to the gbm package. Update, 1(1): 2007, 2007.
[39] Noah Simon, Jerome Friedman, Trevor Hastie, and Rob Tibshirani. Regularization paths for Cox's proportional hazards model via coordinate descent. Journal of Statistical Software, 39 15):1, 2011.
[40] Terry Therneau, Cindy Crowson, and Elizabeth Atkinson. Using time dependent covariates and time dependent coefficients in the Cox model. Survival Vignettes, 2:3, 2017.
[41] Terry M. Thernea,. A Package for Survival Analysis in R, 2023. URL https://CRAM:
R- projeeect..o org/ package R package version 3.5-7.
[42] Andrea Vedaldi and Andrew Zisserman. Efficient additive kernels via explicit feature maps.
IEEE Transactions on Puttern Analysis and Machine Intelligence, 34(3):480-45) 2012.
[43] Sreekanth Vempati, Andrea Vedaldi, Andrew Zisserman, and C. Jawahar. Generalized RBF feature maps for efficient detection. In BMVC, pages 1-11, 2010.
[44] Grace Wahba. Spline Models for Observational Data, volume 59. SIAM, 1990.
[45] Christian J. Walder and Adrian N. Bishop. Fast Bayesian intensity estimation for the permanen- tal process. In International Conference on Machine Learning, volume 70, pages 3579-3588. PMLR, 2017.
[46] Ping Wang, Yan Li, and Chandan Reddy. Machine learning for survival analysis: A survey. ACM Computing Surveys (CSUR), 11(6):1-3,, 2019.
[47] Shannon Wongvibulsin, Katherine C. Wu, and Scott L. Zeger. Clinical risk prediction with random forests for survival, longitudinal, and multivariate (RF-SLAM) data analysis. BMC Medical Research Methodology, 20(1):1-14, 2020.
[48] Weichi Yao, Halina Frydman, Denis Larocque, and Jeffrey S. Simonoff. Ensemble methods for survival function estimation with time-varying covariates. Statistical Methods in Medical Research, 31 (1(11))2217-2236 2022.
[49] Qixian Zhong, Jonas W. Mueller, and Jane-Ling Wang. Deep extended hazard models for survival analysis. In Advances in Neural Information Processing Systems 34, 2021.
13

---

A Derivation of MAP Estimator
We detail the derivation of the MAP estimator should be zero MAP estimator B(y):
The functional derivative of ((x(y))
/
8S
5S
85 (S(a f(y))
8x (y)
55(y)
dy + 0((xx)))
(5
NV
2
2
A
(
5(y (m.) + 2
-2((1) oody
n=1
S
1
2
- -0 43) Σ
2
8(y (n.) + f(y)
body 0,
((v)
where the following relation used,
J 6(y (gdy
J k(v) / k" v'
/ dy dy'dz )z(())))) / k- k" v' (yyydd
/
(u. = ())
Thus the following equation is derived,
B(v) + 2 Σ A Gl - y,)) = Σ
8(y Un), VEW y.
(A1)
(()))
n=1
By applying operator IC to IA),, we obtain a linear integral equation that derives the MAP estimator f(y) as follows,
N Σ
( + 2 Σ A v,)
(y.
y EEV.
(A2)
3=1
n=1
The linearity of the integral equation permits a representation of the form
( 2 h(y (n))(()))'
where 1y y²) is a positive semi-definite kernel that solves integral equation Note that the
derivations follow Kim [23]
B Derivation of Predictive Distribution and Marginal Likelihood
We now know the mode the posterior, h(y), and consider a Taylor expansion of functional action potential ((z(y). x(²)) centered on the mode such that
s(z(g).(()) ~ ((((()(())))))))))(())) (t() + / o o'
Adyyy' (B1)
yxy
where o² oy 3¹) = ă br(y)az(yy) 1== is the second derivative of S. The first term in the Taylor ex- pansion vanishes due to the stationary condition. The quadratic approximation of the action integral corresponds the approximation of the posterior process by GP, and the predictive covariance or the kernel function for the posterior GP, denoted by a (3,, y can be obtained by the functional inversion of o² (')))) which results
52 S(
oly. v') hly, y') ( H) --(2'))',
(B2)
14

---

[II].. The calls coxph to where the definitions of Z and H are taken from (2)). The full derivation of (E2) is given in [24|] When the latent function zyy) follows a posterior GP with mean of f(y) and kernel o y is easily verified that the value of the squared function, À = x² (v), on each point of the covariate domain y E y follows a Gamma distribution defined by (2II).
Furthermore, under Laplace approximation (B.I,, we can obtain the marginal likelihood, p(D), in (IT)) by performing the path integral as,
log (DD log f exp(-S(x(y), x())) - log XC] 2r -S(((y) x(y) + log [
(B3)
where | the functional determinant of integral operator 2' Sy ay y dyy' We can rewrite the result in more tractable form by substituting (III, B2) into (B33)
(() log |Z|
log IN + Z- H
log MM 2A -1 + (log 1)N.
Full derivations are provided [244]
C Model Configuration
C.1 Synthetic Data
Survival Permanental Process (SurvPP)
We set the number of features for Random feature map [37]] (M), learning parameter (Ir), and stop condition (G) for Adam [26]] follows:
M 500, 1-= 0.05, G 10
We applied to SurvPP multiplicative Gaussian kernel
k(y.y') y¹) = III 10(000))))) d=1
(1, v1.
where hyper-parameter was optimized for each data by maximizing the marginal likelihood through grid search. In the experiments on synthetic data, we selected a set of nine values for 0 as the grid points,
0 E {0.1, 0.2, 0.5, 0.7, 1.0, 2.0, 5.0, 7.0, 10.0).
We implemented SurvPP by using TensorFlow-2.10. A MacBook Pro with 12-core CPU (Apple M2 Max) was used, with the GPU inactivated (tf device (' (/cpu:0)) for fair comparison with the benchmarks.
Cox Proportional Hazards Model (CoxPH))
We implemented CoxPH through package survival. coxph (I.GPI--3) fit model and compute a base hazard function were
cfit = coxpt Stop, ~ cov1 + cov2, df)
sfit survfit list cov2 0)),
where df was the survival data in counting process format.
Generalized Boosted Model (GBM)
We implemented GBM through package
[G] (GPL). The call in gbmt to model
was
gfit gbmt SunttSurrrStarrt,ttoe Stop. Event) ~ cov1 + cov2, data = df,
distribution glg_gggm_______ dist("CoxPHH))
cv cv_folds 10, trrinnppaaame = params,
par details gimParallel(num threads = 12)),
15

---

where params represents the hyperparameter. We selected set of nine hyperparameters for the grid search, num trees E (500, 1000, 2000) shrinkage E {0.001, 0.005, 0.01),
and found one that minimized the cross validation error (gfittSaalid.error) where num trees and shrinkage represent the number of trees and the shrinkage/learning rate, respec- tively.
Random Forest-based Model (RFM)
We implemented RFM through package LTRCforests (GPL.) [88]. The call to fit model was rfit Ltrcrrf(Surv(Starrt) Event ) cov1 + cov2, data df, id = ID, mtry = ceiling(10), ntree = 100).
C.2 Real-world Data: PBC
Survival Permanental Process (SurvPP)
We set the number of features for Random feature map [6>]]] (MM), learning parameter (Ir), and stop condition (G) for Adam [26] as follows:
M = 500, Ir 50,
G < 10-5.
We applied to SurvPP multiplicative Gaussian kernel
13
k ky y') = II 10(ppa v(iia
y = (t. V1.... 912),
where the hyper-parameter was optimized for each data by maximizing the marginal likelihood through grid search. In the experiments on synthetic data, we selected a set of nine values for 0 as the grid points,
{0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09}.
Here, we normalized the 13 covariates so that - 0.1 mean(yal)/steee
As in the experiments on synthetic data, GPU off (tf devicev (f/cpu:²)))) for a fair com- parison with the benchmarks.
Cox Proportional Hazards Model (CoxPH)
The calls in coxph to fit model and compute base hazard function were
cfit coxph(Surv(Start,, Stop, Event) ~ age + edema + alk.phos chol + ast
+ spidees hepato ascites + albumin + bili + protime, df)
sfit survfit(cfit, list 0, edema = 0, alk 0, chol 0, ast = 0,
platelet = 0, spiders = 0, hepato 0, ascites = 0, al almumin O, bili protime 0)),
where df was the survival data in counting process format.
Generalized Boosted Model (GBM)
The call in gbmt to fit model was
gfit=
gbmt (SrttSuurSStartt,S.. Stop, Event) age + edema + alk.phos + chol 4 ast +platelet + spiders + hepato + ascites + al albumin bili +protime, data =f, df, distribution gbm bim__dis((( Ailst(Cossii "CooxH"'))
cv cv_folds 10, tettttrrarareeeee = params,
par details = gbmParallel(num threads = 12)),
16

---

where params represents the hyperparameter. We selected set of nine hyperparameters for the grid search, num trees E (500, 1000, 2000) shrinkage E {0.001, 0.005, 0.01).
and found the one that minimized the cross validation error (gfittvalid.eerorr) where
num trees and shrinkage represent the number of trees and the shrinkage/learming rate, respec- tively.
Random Forest-based Model (RFM)
The call to fit a model was
rfit = 1trcrrf (SrrvSttu Stop, Event) ~ age + edema + alk phos + chol + ast +platelet + spiders + hepato + ascites albumin bili tprotime, data = df, id ID, stepFactor 1.5).
C.3 Real-world Data: SANAD
Survival Permanental Process (SurvPP)
We set the number of features for Random feature map [37]] (M), learning parameter (Ir), and stop condition (G) for Adam [26]] as follows:
M = 500,
Ir 10, G 10-5,
We applied to SurvPP a multiplicative Gaussian kernel
ky 3/) 0.1
II --(ouuu-val) y (t, .t.
where hyper-parameter was optimized for each data by maximizing the marginal likelihood through grid search. In the experiments on synthetic data, we selected set of nine values for as the grid points,
0 e {0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4). Here, we normalized the covariates so that ya - 0. 0.1 (ya - mean(yal)
As in the experiments on synthetic data, GPU was set off (tffdeeiie(( (</cpu:0)) for a fair com- parison with the benchmarks.
Cox Proportional Hazards Model (CoxPH)
The calls in coxph to fit a model and compute base hazard function were
cfit coxph(SurvSSaail) Stop, Event) ~ age + gender + treat + dose, df)
sfit survfit(cee cfit, list(age gender 0, treat 0, dose 0)),
where df was the survival data in counting process format.
Generalized Boosted Model (GBM)
The call in gbmt to fit a model was
gfit- = gbat Stop, Event) ~ age + gender + treat + dose,
data df, distribution gbm dist("CoxPP)) )
cv cv_folds = 10, train params params.
par details gomParallel (num threads 12)).
where params represents the hyperparameter. We selected a set of nine hyperparameters for the grid search,
num trees E (500, 1000, 2000) x shrinkage {0.001, 0.0005,0.011
and found the one that minimized the cross validation error (gfitSvalid.error) where num trees and shrinkage represent the number of trees and the shrinkage/learning rate, respec- tively.
17

---

TLL
AUC
0.1
0.9
SurPPP
Sim10
SurvP_ Rie20
0.0
0.8
0.1
0.7
0.2
0.3
0.5
0.3
0.5
0.7
0.3
0.9
Evaluation Point
Evaluation Point
Figure D1: Performances on the dataset of non-linear hazard function Anon (f) for different numer- ical integration rules. SurvPP_Riis and SurvPP_Simg represent SurvPPs with Riemann sum and Simpson's s rule approximations of integration, respectively. is the number of representative time points for covariate, Ju.
Random Forest-based Model (RFM)
The call to fit a model was
rfit ltrerrf (SureSttar,t Stop, Event) ~ age + gender treat + dose, data df, id ID, stepFactor 1.5).
D Additional Experiments
D.1 More Accurate Evaluation of Equivalent Kernel
The functional form of covariate map, yt), involves the sum of outer products of feature maps, A, defined by Equation (IG), and A can be considered as Riemann sum approximation of the time integral follows:
Σ L To ((..()))(
A
Σ ((g,)
(D1)
Thus, when the covariate map, (t), and the feature map, d(y), are both smooth, the Riemann sum approximation can be replaced by more accurate approximation,
ply (l)(())))) di Σ (u 1=1
(D22
where the weights, Wy. depend on the approximation rule. This time, adopt Simpson's 's 1/3 rule, and add an experiment on the synthetic nonlinear data, Anon (4), to check for an improvement in accuracy. Note that Simpson's 1/3 rule corresponds to quadratic interpolation. Figure DD shows that adopting Simpson's 's rule achieved better predictive performance on TLL with a smaller discretization number (Ju).
D.2
Experiment on Independent Validation Data Sets
In Section H. we evaluated the performances on synthetic data with the cross-validation approach because it is common in the machine learning literature. But synthetic data offers the unique ad- vantage of being able to generate as much data as needed, and we here conducted the performance assessment with the independent validation approach. Figure D displays the result on independent validation sets of synthetic data (U 900 for each of 10 training data, and 100 for each of 10 test data), showing the result similar to the rros-validation approach (Figure (IH2).
D.3 Experiment on Larger Data Sets
To examine the computation scalability of the compared models versus the number observed events N, created data sets with user size U e {103, 10 5-10", 100 according the nonlinear
18

---

02-
0.0
COM
Sundo
-02
-04
-04
0.6
0.9
AUC A +
AdD-u 20
1.0
AURC.A
AUC J -00
0.8
0.4
02
0.3
0.5
0/7
D-7 0.9
05
0.5
D-7
0.9
0-7
0.
Evoluation Point
Evaluation Point
Evalusion Point
Evaluation Paint
Figure D2: Performance on the independent validation synthetic datasets. (A) Log-linear hazard
function Ain (t). (B) Non-linear hazard function Ano (t). The results are consistent with Figure [I-2
10a
RFM
102
GBM
CPU [sec]
101
SurvPP
100
CoxPH
10a1
1
10a
104
10
Number of Observed Events
Figure D3: The CPU times demanded for estimating a hazard function served events. The error bars represent the standard deviations across 10 trials. The dashed line represents line of CPU x N reference. For GBM and SurvPP, the average cpu times over 9- point grid search of the hyperparameter are displayed. The CPU times of GBM and RFM exceeded 103 seconds with N 10 and the estimations were given up.
the number of ob-
scenario (see Section
Xnon (t) h(t) exp ((v + )Z())) h(t) 2 23/
which resulted in the data sets with N E (818, 8082, 40660, 81066). Here, we set J, to 10 for the counting process format of data. For each dataset, we randomly split the U individuals into 10 training data, and conducted trials of evaluations of the CPU times demanded for estimating hazard function. Figure DΣ displays the CPU times as function
subgroups, repeated assigning subgroups
of training data. shows that the CoxPH computation clearly scaled linearly with N. while that of SurvPP seems to be a little more than linear. This is because that each iteration of gradient descent algorithm scaled linearly with N, but the number of iterations to meet stop condition G < 10-5 increased moderately with N. Among the non non-parametrie approaches (SurvPP, GBM, and RFM), SurvPP achieved the fastest computation at hundreds of times faster than the others, regardless of the number of observed events N.
D.4 Sensitivity Hyper-Parameters
Kernel methods/GP models are generally sensitive to kernel parameter values, and a effective way of optimizing the kernel parameter is essential. Figure D plots AUC performance on PBC dataset for SurvPPs with different search ranges of kernel parameters: the range of SurvPP, is E (..1, 0.5, 0.7, 1.0, 2.0, 5.0, 7.0, 10.0), which was used in the supplement; the range of SurvPP, is 0 E {0.01, 0.02, 0.03... ... 0.09). Figure DE shows that SurvPP; achieved substantially
19

---

AUC
AUC
1.0
1.0
0.9
0.9
0.8
0.8
CoxPH GBM
CoxPH IBBM
0.7
SurvPP_ 0
0.7
SurVP
2
x10
2
4
x10²
Evaluation Time Point
Evaluation Time Point
Figure D4: Performances on PBC dataset for SurvPPs with different search ranges of kernel param-
eters. (A) The search range of SurvPPo is 0 {0.1, 0.2, 0.2,,0.2,0 0.5, 0.7, 1.0, 2.0, 5.0, 7.0, 10.0). (B) The search range of 1 is E {0.01, 0.02, 0.03, 0.09}.
A
TLL
AUC
CPU [sec]
-1.4
-1.6
0.8
10'
0.6
1..
2.0
0.4
10
2.2
0.2
GBM
SurvPP
18
19
20
21
x10²
18
19 Evaluation
20
21 x10
GBM SurVPP
Evaluation Time Point
Figure D5: Performances on MIMIC III dataset. (A) Box plot of TLL, and AUC as functions of evaluation point: the higher, the better. (B) The CPU times demanded for estimating hazard function. The error bars represent the standard deviations. The average cpu times over 9-point grid search of the hyperparameter are displayed.
better AUC performance than SurvPPo. where SurvPP, displayed in Figure more sophisti- cated algorithm for hyperparameter optimization could enhance the performance of SurvPP, which the next step our study.
D.5 Experiment on ICU Data Set
We examined the validity of SurvPP against GBM on MIMIC-III Clinical Database (MIMIC III,, the large publicly available dataset over 50,000 ICU admissions from the Beth Israel Deaconess Medical Center [I1%], where events were deaths. We extracted admissions and measurements from Vue (ITEMID) such that admissions shared 14 measured covariates with each other, which sulted in 133 admissions. We adopted sex static covariate and the 14 measurements as time- varying covariates. We randomly split the individuals into subgroups, assigned one to test and others to training data, and conducted 3-fold cross evaluation of the predictive performances. The model configuration follows the experiment on PBC data set (see Appendix I))
Figure D displays the predictive performance on MIMIC III. It shows that SurvPP achieved better TLL performance than GBM, and SurvPP achieved comparative AUC performance with GBM. The smaller variance of SurvPP's AUC implies that SurvPP works more robustly than GBM, but further experiments are needed to investigate
20


Reviews:
[{"review_id": "2rxqgViCNs", "reviewer": "Reviewer_4Stf", "summary": "This paper introduces a new hazard function in survival analysis. The hazard is a squared Gaussian process, where the index set is the covariates at a given time. They assume that the covariates are piecewise constant, propose a MAP estimator and derive the predictive distribution and marginal likelihood. They use ideas from kernel methods and achieve computational complexities that are lower than alternate approaches, scaling linearly with the number of observed events.", "questions": "Are your benchmarks also assuming piecewise constant time-varying covariates?", "limitations": "While the limitation they mention are important, some of the limitations that I mentioned in the weaknesses are not addressed.", "rating": 3, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The idea of connecting this to kernel methods to reduce computational complexity is quite clever, and the speedups over GBM and RFM are impressive.", "weaknesses": "The piecewise constant time-varying covariate assumption is very strong. You should mention this.\n\nThere also needs to be some discussion of exogeneous vs endogeneous time-varying covariates. It looks like you're implicitly assuming that your covariates are exogeneous, but you should make that clear.\n\nYou need a real data experiment in the main paper. In the introduction, you mention 'confirm that our algorithm achieves comparable pre67 dictive accuracy while being hundreds of times faster than state-of-the-art survival models' but saying nothing about real data. You have something in the supplement, but again, there's no reference to how good the results are in the main paper. Given that, one has to give less weight to the supplement.\n\nIn your benchmarks, I assume that all of them are assuming piecewise constant time-varying covariates. You should include something using interpolation: e.g. linear mixed models with splines. They do this a lot in the joint modeling literature.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for the comment, to which we provide a detailed response below. Especially, the critical concern about the piecewise constant assumption of covariates arises from our confusing explanations or introductions, and we believe that we can fully dispel the concern. We hope that our explanation is satisfactory. \n\n**The piecewise constant time-varying covariate assumption is very strong. You should mention this.**\n**You should include something using interpolation: e.g. linear mixed models with splines. They do this a lot in the joint modeling literature.**\n**Are your benchmarks also assuming piecewise constant time-varying covariates?**\n\nThe piecewise constant time-varying covariate assumption, which is referred as the *counting process format of input* in the paper, approximates the continuous covariate function over time by a set of representative points (RPs) at finite time points. There are two ways of defining the RPs: one is as the points that were observed, which might be sparse over time due to measurement constraints ; the other is as the denser points obtained by using interpolation methods (e.g., splines). We think that the reviewer assumes the former definition and concludes that the covariate function should be approximated by denser points. But our model (SurvPP) and the benchmarks (CoxPH, GBM, RFM) assume the latter definition in the paper, and thus various numbers of splitting ($J_u$) were considered in the synthetic data experiments. Therefore, the counting process format of input is not so strong because we can adopt RPs in any density if necessary. Note that our model achieves an efficient computation regarding $J_u$ (see Figure 1-2). \n\nAlso, we consider the possibility that the reviewer suggests us to utilize the functional form of interpolated covariate function, without discretizing the interpolated covariate function. The suggestion makes sense, but unfortunately, we cannot find such approaches in the joint modeling literature: covariate function needs to be discretized over time in RNN-based approaches [a,b]; in Cox proportional hazards-based models, the non-parametric estimator of baseline intensity is zero except at observed event times, which results in the discretization of covariate function for likelihood function evaluation [c]. However, thanks to the reviewer's suggestion, we found that our proposed algorithm could be improved as follows. The functional form of (interpolated) covariate function involves the sum of outer products of feature maps, denoted by $\\mathbf{A}$ in Eq. (16), and it can be considered as a Riemann sum approximation of the time integral,\n\n$$\n\\mathbf{A} = \\int_0^T \\phi(\\mathbf{y}(t)) \\phi(\\mathbf{y}(t))^{\\top} dt \\sim \\sum_{j=1}^J \\Delta_j \\mathbf{\\phi}(\\bf{y}_j)\\mathbf{\\phi}(\\mathbf{y}_j)^{\\top}.\n$$\n\nGiven a smooth covariate function, the Riemann sum approximation can be replaced by a more accurate approximation,\n\n$$\n\\mathbf{A} = \\int_0^T \\phi(\\mathbf{y}(t)) \\phi(\\mathbf{y}(t))^{\\top} dt \\sim \\sum_{j=1}^J w_j \\mathbf{\\phi}(\\mathbf{y}_j)\\mathbf{\\phi}(\\mathbf{y}_j)^{\\top},\n$$\n\nwhere the weight, $w_j$, depends on the approximation rule. This time, we adopt the Simpson's 1/3 rule, and add an experiment with the synthetic nonlinear data to examine the improvement in accuracy. Note that the Simpson's 1/3 rule corresponds to a quadratic interpolation. Figure R1 (see a pdf file submitted) shows that adopting the Simpson's rule achieved a better predictive performance in TLL with a smaller discretization number ($J_u$). We will add the above discussions and results in the manuscript. We appreciate the reviewer's valuable comment.\n\n**There also needs to be some discussion of exogeneous vs endogenous time-varying covariates. It looks like you're implicitly assuming that your covariates are exogeneous, but you should make that clear.**\n\nAs the reviewer pointed out, we consider exogenous covariates in the paper. Although we discussed about the topic in Related Work (Two Scenarios in ...), we found that the discussion was inadequate. According to the reviewer's suggestion, we will revise the paragraph to explain the following points: the definitions of exogenous and endogenous covariates [d]; the primary purpose of survival analysis with exogenous covariates is to estimate a hazard function as an explicit function of covariates, and to make predictions of failure times under various possible *future* covariate functions (what-if analysis); the primary purpose of survival analysis with endogenous covariates is to make predictions of failure times by using the *past* observations of covariates, where joint modeling approaches, which jointly model the stochastic process of covariates and failure times, have been developed intensively; we consider exogenous covariates, and survival models for exogenous covariates cannot be compared directly with those for endogenous covariates because the tasks are different with each other.\n\n**You need a real data experiment in the main paper.**\n\nAccording to the reviewer's suggestion, we will put in the main paper the real data experiment (now in the supplement). Furthermore, we added one more experiment on real-world data (for details, see Figure R2 and our response to Reviewer LBpm).  We will use an extra page given for accepted papers. \n\n[a] Lee et al. Dynamic-DeepHit: A deep learning approach for dynamic survival analysis with competing risks based on longitudinal data. IEEE Transactions on Biomedical Engineering, 67(1):122–133, 2019.  \n[b] Nagpal et al. Deep parametric time-to-event regression with time-varying covariates. In Survival Prediction-Algorithms, Challenges and Applications, 184–193, 2021.  \n[c] Jarrett et al. Dynamic prediction in clinical survival analysis using temporal convolutional networks. IEEE Journal of Biomedical and Health Informatics, 24(2):424–436, 2019.  \n[d] Kalbfleisch and Prentice. The Statistical Analysis of Failure Time Data. John Wiley & Sons, 2011.  "}, {"review_id": "QgzK2KEDmA", "reviewer": "Reviewer_LBpm", "summary": "This paper addresses a survival analysis problem with time-varying covariates by proposing a novel Bayesian approach to estimate a hazard function as a nonlinear function of covariates using the kernel trick. MAP, predictive distribution, and computational complexity of the proposed approach were presented. Experiments were done on both synthetic and real-world benchmark datasets.", "questions": "- The authors might not have enough space to put the experimental results on the real-world benchmark dataset. Would the authors plan to move it to the manuscript?\n\n- Is there any specific reason that the authors used AUC rather than C-index, which is more commonly used in the survival analysis?\n\n- It is interesting that SurvPP, which is a non-parametric model, showed as good TLL as CoxPH, which is a (semi) parametric model, unlike the other nonparametric tree-based approaches. Could authors give some explanation why this happened? Is this something expected?\n\n- I am also curious about what each method's survival curve looks like. Might be interesting if the authors could compare the average prediction curve against the KM survival curve.\n\n- As the proposed method is a Bayesian approach, it can provide more than just better prediction performance. Maybe the authors could add some discussion on that.", "limitations": "Did not find any other limitations than the authors discussed in the last section.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- Well-written and easy to follow. The relationship with the related work is clearly presented.\n\n- Proposed a novel Bayesian approach to tackle the survival analysis with time-varying covariates that allows modeling of a nonlinear relationship between hazard and covariates. The authors presented MAP, predictive distribution, and computational complexity of the proposed approach.\n\n- Experimental results clearly support that SurvPP provides good predictive performance as well as is computationally efficient.", "weaknesses": "- Just as other GP or kernel-based methods, SurvPP might not be ideal for complex high-dimensional data. But the authors have pointed this out, and a potential way to handle this has been clearly discussed.\n\n- In the experimental results on the real-world benchmark dataset, SurvPP showed slightly worse AUC than the other models, unlike in synthetic data. Could the authors kindly give some explanation why? Interestingly, in terms of TLL, SurvPP is doing quite well.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for the highly positive comments and suggestions, by which we are strongly encouraged. Below we provide a detailed response to each of the comments.\n\n**The authors might not have enough space to put the experimental results on the real-world benchmark dataset. Would the authors plan to move it to the manuscript?**\n\nYes, we will put in the main paper the real data experiment. Furthermore, we added one more experiment on real-world data (for details, see Figure R2 and our response to Reviewer LBpm).  We will use an extra page given for accepted papers. Sorry for the inconvenience. \n\n**Is there any specific reason that the authors used AUC rather than C-index, which is more commonly used in the survival analysis?**\n\nC-index considers the concordance between the ranking of the test results and that of the observed times-to-events, where the ranking of the test results is assumed to be static. It is not appropriate as an evaluation criterion when covariates are time-varying because the ranking of the test results can also be time-varying. As is clear from the definition in the paper (Eq.(28)), the dynamic AUC essentially estimates the C-index at each time, and is commonly used in the literature of survival analysis with time-varying covariates. Note that it is known that a weighted average of the dynamic AUC equals C-index [e]. We found that our introduction was inadequate from the reviewer comment, and will add the above explanation in the main paper.\n\n[e] Pantoja-Galicia et al. Concordance measures and time-dependent ROC methods. Biostatistics & Epidemiology 5(2): 232-249, 2021.\n\n**In the experimental results on the real-world benchmark dataset, SurvPP showed slightly worse AUC than the other models, unlike in synthetic data. Could the authors kindly give some explanation why?**\n\nThe result that SurvPP showed slightly worse AUC than the other models could be mainly due to the inadequate optimization of kernel parameters. Figure R3 (see a pdf file submitted) plots the AUC performances on PBC dataset for SurvPPs with different search ranges of kernel parameters: the range of SurvPP_0 is $\\theta \\in \\\\{ 0.1, 0.2, 0.5, 0.7, 1.0, 2.0, 5.0, 7.0, 10.0 \\\\}$, which was used in the supplement; the range of SurvPP_1 is $\\theta \\in \\\\{ 0.01, 0.02, 0.03, \\dots, 0.09 \\\\}$. Figure R3 shows that SurvPP_1 achieved substantially better AUC than SurvPP_0, and SurvPP_1's AUC is now comparable to the other models on average. More sophisticated algorithm for hyperparameter optimization could enhance the performance of SurvPP, which is the next step of our study. We will add the above discussion in the supplement.\n\n**It is interesting that SurvPP, which is a non-parametric model, showed as good TLL as CoxPH, which is a (semi) parametric model, unlike the other nonparametric tree-based approaches. Could authors give some explanation why this happened? Is this something expected?**\n\nThank the reviewer for the good point of view. SurvPP treats the duration from the entry, $t$, as a time-varying covariate, and estimates the hazard function of $t$ with the smoothness specified by the kernel function. But the tree-based approaches, GBM and RFM, estimate the base hazard function $h(t)$ (Eq.(25)) by the Breslow/Fleming-Harrington estimator, where $h(t)$ is treated as piecewise constant between uncensored failure times [f]. Therefore, if the underlying base hazard function is smooth and the number of observed failure times is small, then GBM and RFM could estimate the hazard function regarding $t$ poorly compared to SurvPP. TLL evaluates the base hazard function directly, while AUC disregards the contribution of the base hazard function common to all individuals, which is a possible reason of the obtained results. \n\nAlso, CoxPH's very high AUC ($\\gtrsim 0.9$) in PBC dataset suggests that the underlying generative process could be consistent with the log-linear and proportional assumptions of CoxPH, in which the semi-parametric model, CoxPH, is likely to achieve equal or better performances than the nonparametric model, SurvPP. As a response to another reviewer, we added one more experiment on real-world data, showing that SurvPP performed better than CoxPH when the log-linear assumption is not valid (for details, see Figure R2 and our response to Reviewer 3BPs). \n\nWe will add the detailed discussions about the experimental results in the main paper.\n\n[f] Lin. On the Breslow estimator. Lifetime data analysis 13:471-480, 2007.  \n\n**As the proposed method is a Bayesian approach, it can provide more than just better prediction performance. Maybe the authors could add some discussion on that.**\n\nThank the reviewer for a good suggestion. We will add the following discussion. \n\nOne of the advantages of Bayesian approaches over non-Bayesian ones is that they enable us to perform uncertainty evaluations. More concretely, the posterior process of SurvPP is also given by a permanental process, and thus we can generate sample paths of the estimated hazard function of covariates, $\\lambda (\\mathbf{y})$. Using the samples, we can evaluate the distribution of survival function and perform risk-aware survival analyses."}, {"review_id": "gfNU69iiYI", "reviewer": "Reviewer_3BPs", "summary": "The paper proposes a non-parametric Bayesian survival model called Survival Permanental Process (SurvPP) for analyzing survival data with time-varying covariates. The model is based on a permanental process that defines the latent hazard function on covariate space and can handle right-censored observations in a counting process format. The paper shows that SurvPP holds a representer theorem, which offers a fast Bayesian estimation algorithm that scales linearly with the number of observed events without relying on Markov Chain Monte Carlo computation. The paper evaluates SurvPP on synthetic data and shows that it achieves comparable predictive accuracy while being hundreds of times faster than state-of-the-art methods (specifically, CoxPH and GBMs).", "questions": "Although this paper is well written and technically sound, the reviewer believes that the paper has very limited practical utility. Below are the points that seem serious limitations, which the authors need to address:\n\n**1:** While the paper presents a novel approach to survival analysis using permanental processes and provides a comparison with conventional models (Cox-PH and GBM), the experimental validation appears to be limited to synthetic datasets with only two-dimensional covariates. The reviewer would like to see how the proposed model performs on real-world datasets, which often have higher dimensionality and more complex structures. Specifically, it would be beneficial to test the model on widely recognized and publicly available datasets such as the MIMIC III (Medical Information Mart for Intensive Care III) or the U.K. Cystic Fibrosis Registry. These datasets contain a rich variety of covariates (that do vary over time) and present real-world challenges that synthetic datasets may not fully capture. Using such datasets would not only provide a more rigorous test of the model's performance but also make the results more relatable to practitioners in the field. It would also allow for a more direct comparison with other models in the literature that have been tested on these datasets.\n\nThe authors should explain why they chose to use these synthetic datasets for their experiments, by motivating the particular covariate functions they used–Equations 26, 27, etc.--explaining why they match a realistic situation. Additionally, they should extend the experimental validation to include real-world datasets, to further demonstrate the applicability and robustness of the proposed model.  This would address some other concerns–eg, how to determine the appropriate kernel to use, etc.\n\n**2:** Additionally, the scalability of the proposed model in the context of high-dimensional, low-sample size problems remains an area that warrants further exploration. While we see the authors have acknowledged this limitation (Section 5, lines 308 onwards), the current treatment of this issue does not fully address the potential implications for the practical application of the proposed model. In real-world scenarios, datasets often exhibit high dimensionality, and it is common to encounter situations where the sample size is relatively small. These conditions pose significant challenges for many machine learning models, and it is crucial to understand how the proposed model would perform under such circumstances.The Bayesian estimation framework proposed in this paper adds another layer of complexity to this issue.. Thus, a more detailed discussion of how the proposed Bayesian estimation framework would scale with data dimensionality would be highly beneficial. Therefore, the authors are advised to expand their discussion on this topic. Specifically, it would be valuable to include theoretical considerations, potential strategies for handling high-dimensional data, and, if possible, empirical results from preliminary experiments with (maybe synthetic) high-dimensional, low-sample size datasets.\n\nIn addition to a discussion about the sample efficiency, the paper should also discuss the learner’s computational efficiency.\n\n**3:** The performance of the proposed model is compared with two benchmark methods, namely the Cox proportional hazards model (CoxPH) and the generalized boosted regression model (GBM). While these comparisons provide valuable insights into the model's performance, the choice of benchmark methods raises some questions. Specifically, the CoxPH and GBM models, although well-established in the field, do not represent the most recent advancements in survival analysis. It would be interesting to see how the proposed model compares with more recent methods that leverage advanced machine learning techniques. Specifically, the authors should compare their approach with the following methods:\n\n*Nagpal, Chirag, Vincent Jeanselme, and Artur Dubrawski. \"Deep parametric time-to-event regression with time-varying covariates.\" Survival Prediction-Algorithms, Challenges and Applications. PMLR, 2021.*\n\n*C. Lee, J. Yoon and M. v. d. Schaar, \"Dynamic-DeepHit: A Deep Learning Approach for Dynamic Survival Analysis With Competing Risks Based on Longitudinal Data,\" in IEEE Transactions on Biomedical Engineering, vol. 67, no. 1, pp. 122-133, Jan. 2020, doi: 10.1109/TBME.2019.2909027.*\n\n*Tomašev, Nenad, et al. \"Use of deep learning to develop continuous-risk models for adverse event prediction from electronic health records.\" Nature Protocols 16.6 (2021): 2765-2787.*\n\nComparing the proposed model with these more recent methods would provide a more comprehensive evaluation of its performance and relevance in the current research landscape. It would also allow for a more direct comparison with state-of-the-art models, enhancing the paper's contribution to the field.\n\n**4:** The authors have presented their experimental results in a cross-validation setting for synthetic datasets. While cross-validation is a widely accepted method for model evaluation, particularly when dealing with limited data, it is not immediately clear why this approach was chosen for synthetic data, where the volume of data can be controlled. Synthetic data offers the unique advantage of being able to generate as much data as needed, which could potentially allow for the creation of independent validation sets. This approach would provide a more stringent test of the model's generalization ability, as it would involve testing on data that the model has not seen during training, even indirectly through cross-validation folds.\n\nCould the authors provide some insights into their decision to use cross-validation for synthetic data? Furthermore, would it be possible to conduct additional experiments using independent validation sets with synthetic data to further assess the model performance?", "limitations": "N/A", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 1, "strengths": "The main strengths of this paper are listed below:\n\n1. The paper introduces a novel non-parametric model, called SurvPP, for estimating hazard functions based on survival data with time-varying covariates. This model addresses the limitation of the Cox proportional hazards model (CoxPH) by allowing for non-linear dependence between survival times and covariates, and also where these covariates can change, after the time when the predictions are made.\n\n2. The paper shows that SurvPP holds a representer theorem, which means it can be trained efficiently using a fast Bayesian estimation algorithm that scales linearly with the number of observed events without relying on Markov Chain Monte Carlo computation.\n\n3. The authors have conducted experiments on two synthetic datasets to demonstrate that their proposed SurvPP outperforms CoxPH and another non-parametric model called Generalized Boosted Regression Models (GBM) in terms of predictive performance. SurvPP achieves significantly better performance than CoxPH, but the performance gaps between SurvPP and GBM are marginal and not significant.\n\n4. The paper shows that SurvPP is computationally efficient, particularly when covariates are measured frequently. SurvPP can be trained faster than GBM and CoxPH in estimating hazard functions when covariates are measured frequently.", "weaknesses": "(See next entry: Questions )", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for carefully reading our paper and giving valuable comments, to which we provide a detailed response below. We believe that we can fully address all the concerns. \n\n**1: the experimental validation appears to be limited to synthetic datasets with only two-dimensional covariates. The reviewer would like to see how the proposed model performs on real-world datasets, which often have higher dimensionality and more complex structures....**\n\nWe totally agree with the reviewer about the importance of validation on real-world datasets. Actually, we performed an experiment on a real-world dataset (PBC dataset), but put the result in the supplement (Section S4) due to space limitation. We apologize for the inconvenience. According to the reviewer's suggestion, we added one more experiment on another real-world data, and will put in the main paper the following two real data experiments.\n\nPBC dataset has twelve covariates, which we believe alleviates the reviewer's concern about the dimensionality of covariates. But Figure S2 in the supplement shows that CoxPH achieved very high AUCs ($\\gtrsim 0.9$), suggesting that the underlying generative process could be consistent with CoxPH's simple assumption of log-linearity and proportionality. Then, the semi-parametric model, CoxPH, is likely to achieve equal or better performances than the nonparametric model, SurvPP, being consistent with the result. \n\nThis time, we ran an additional experiment on a real-world dataset, *Standard And New Antiepileptic Drugs study data* (SANAD), provided by R package joineR (GPL-3). SANAD was an unblinded randomized trial recruiting patients with epilepsy for whom carbamazepine (CBZ) was considered to be standard treatment and they were randomized to CBZ or the newer drug lamotrigine (LTG), where 605 patients were included and event was the time to treatment failure. We adopted a time-varying covariate (calibrated dose), and three static covariates (age of patient, gender, and randomized treatment (CBZ or LTG)). The procedure of 10-fold cross evaluation of performance follows that of the experiment on PBC data. Figure R2 plots the result, showing that SurvPP achieved better performance than CoxPH, and achieved comparable performance while being substantially faster than GBM.\n\nMIMIC III and U.K. Cystic Fibrosis Registry, kindly suggested by the reviewer, seem to require more than a week for access approval. Thus, we gave up testing them in the rebuttal period. We hope that it is acceptable to the reviewer.  \n\n**The authors should explain why they chose to use these synthetic datasets for their experiments, by motivating the particular covariate functions they used–Equations 26, 27, etc...**\n\nWe prepared the synthetic datasets (Eq.(26-27)) for evaluating the methods' predictive and time performances in the condition that time-varying covariates have linear and nonlinear effects on the hazard function. Synthetic data with similar nonlinearity has been adopted in the literature [g], and we believe that the synthetic datasets are in line with the purpose of this paper.\n\n[g] Katzman et al. DeepSurv: personalized treatment recommender system using a Cox proportional hazards deep neural network. BMC medical research methodology, 18(1):1-12, 2018.\n\n**2: the scalability of the proposed model in the context of high-dimensional, low-sample size problems remains an area that warrants further exploration...Thus, a more detailed discussion of how the proposed Bayesian estimation framework would scale with data dimensionality would be highly beneficial.**\n\nThank the reviewer for an important point of view! One of the biggest advantages of kernel method/GP is the good scalability with data dimensionality, $d$, from which SurvPP have benefitted fully: the complexity of SurvPP regarding $d$ is $\\mathcal{O}(NMQ+(N+J)M^2+M^3+dJM)$, where $\\mathcal{O}(dJM)$ stems from the computation of feature maps, $\\phi_m(\\mathbf{y})$. As is $J$, $d$ is not part of the iterative optimization ($Q$), resulting in a very feasible computation regarding the data dimensionality. We will add the discussion in the main paper, which highlights the merit of our model.\n\n**3: the choice of benchmark methods raises some questions. Specifically, the CoxPH and GBM models, although well-established in the field, do not represent the most recent advancements in survival analysis...**\n\nAll the references raised by the reviewer are joint modeling approaches for endogenous covariates, while our SurvPP is for exogeneous covariates. Please see our 2nd response to Reviewer 4Stf about the difference between endogenous and exogenous covariates: SurvPP cannot be compared directly with the approaches for endogenous covariates because the tasks considered are different with each other. Also, it is notable that as well as CoxPH and GBM, we adopted RFM as a benchmark, which is one of the most recent advancements in survival analysis for exogenous covariates published in 2020 (arXiv) and in 2021 (journal).\n\n**4: it is not immediately clear why this approach (cross-validation) was chosen for synthetic data...Could the authors provide some insights into their decision to use cross-validation for synthetic data? Furthermore, would it be possible to conduct additional experiments using independent validation sets with synthetic data to further assess the model performance?**\n\nWe chose the cross-validation approach just because it is common in the machine learning literature. We are pleased to follow the reviewer's suggestion. Figure R4 displays the result on independent validation sets of synthetic data ($U$ = 900 for each of 10 training data, and $U$ = 100 for each of 10 test data), showing the result similar to the cross-validation approach (Figure 1-2). We will add the result in the supplement. If the reviewer is interested in more large samples or trials, we promise to add experiments in the final manuscript."}, {"review_id": "xdxE6TR8iA", "reviewer": "Reviewer_SVLR", "summary": "This manuscript proposes a Gaussian process based method for survival modeling (or time-to-event modeling) assuming time-varying covariates.  Authors assume covariates are time-varying and they are treated using the counting process format which assumes covariates remain fixed until to the next measurement time point.  The hazard function (or square root of it) is assumed GP prior. Authors rely on previous work on path integral representation of GPs and derive a maximum a posterior estimator of the latent function that serves the hazard function for the likelihood/point process. The method is evaluated on a simulated data set with the log-linear as well as non-linear dependence on covariates.  On simulated data the proposed method performs on par with previous methods but has some beneficial performance in terms of computational efficiency. ", "questions": "Authors note in Conclusions section (and in Appendix) that the kernel hyper-parameters are optimised using the grid-search. Is the proposed method sensitive to kernel parameter values?  Is it due to the grid search that you need to assume the same length scales for all input features, and is that limiting the performance?", "limitations": "I believe authors have adequately described limitations.", "rating": 4, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "Authors study the important question of time-to-event prediction in the case of time-varying covariates. Model formulation that assumes GP prior for the (square of the) hazard function and development of the methodology using path integration representation of GPs in survival modeling appear novel. Quality and clarity of the work appears very good. Problem setting is important. ", "weaknesses": "The development of the maximum a posterior estimate of the latent function as such relies strongly on the previous work (e.g. ref. [19]) and it is unclear what part of the method development is actually novel. The proposed method is tested only on a single simulated data setting with linear and non-linear dependency on time-varying covariates.  On simulated data, performance of the model (in terms of accuracy) is only comparable relative to previous methods.  The proposed method is tested primarily on a very low-dimensional setting (2-dimensions).  Additional experiments presented in Appendix show that the performance of the method in a real data set is slightly weaker than that of the standard CoxPH model.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for carefully reading our paper and giving valuable comments, to which we provide a detailed response below. We believe that our explanation is satisfactory.\n\n**The development of the maximum a posterior estimate of the latent function as such relies strongly on the previous work (e.g. ref. [19]) and it is unclear what part of the method development is actually novel.**\n\nThank the reviewer for pointing out our inadequate explanation about the novelty. Although the key derivation of algorithms is based on the methodology proposed in ref. [18,19], our proposed model (SurvPP) is a non-trivial extension of the point process model proposed by ref.[19]: (i) SurvPP can accommodate multiple trials of event sequence data, while ref.[19] assumes one trial of evet sequence; (ii) in SurvPP, the end time of observation is a stochastic variable which depends on the time of event occurring, while ref.[19] assumes that the end time of observation is given. We discovered for the first time that the representer theorem, a beneficial property for functional analysis, holds for such a complicated point process, that is, survival analysis models, which is a remarkable technical contribution of our paper. We believe that the result makes a non-incremental impact on the survival analysis community.\n  \n**Is the proposed method sensitive to kernel parameter values? Is it due to the grid search that you need to assume the same length scales for all input features, and is that limiting the performance?**\n\nKernel methods/GP models are generally sensitive to kernel parameter values, and a clever way of optimizing the kernel parameter is essential. Gaussian kernels, which was used in the paper, naively require a scale parameter for each dimension of data, resulting in a high dimensional kernel parameter. The grid search becomes prohibitively costly with high dimensional parameter, but a well-known approach to the problem is to normalize (e.g., centering and scaling) each dimension of the data and put a common scale parameter across all dimensions of normalized data, which we adopted in the paper. This approach empirically works robustly, but a more sophisticated approach could improve the performance of SurvPP, which is an important next step of our study. We will add the discussion in the manuscript.\n\n**The proposed method is tested primarily on a very low-dimensional setting (2-dimensions).**\n\nIn the real-world data experiment on PBC dataset, the proposed method was tested on 12-dimensional covariates, which we believe alleviates the reviewer's concern about the dimensionality of covariates. We will move it to the main paper. \n\n**Additional experiments presented in Appendix show that the performance of the method in a real data set is slightly weaker than that of the standard CoxPH model.**\n \nThe result in Figure S2 that SurvPP showed slightly worse AUC than the other models could be mainly due to the inadequate optimization of kernel parameters. Figure R3 (see a pdf file submitted) plots the AUC performances on PBC dataset for SurvPPs with different search ranges of kernel parameters: the range of SurvPP_0 is $\\theta \\in \\\\{ 0.1, 0.2, 0.5, 0.7, 1.0, 2.0, 5.0, 7.0, 10.0 \\\\}$, which was used in the supplement; the range of SurvPP_1 is $\\theta \\in \\\\{ 0.01, 0.02, 0.03, \\dots, 0.09 \\\\}$. Figure R3 shows that SurvPP_1 achieved substantially better AUC than SurvPP_0, and SurvPP_1's AUC is now comparable to the other models on average. More sophisticated algorithm for hyperparameter optimization could enhance the performance of SurvPP, which is the next step of our study. \n\nAlso, Figure S2 shows that CoxPH achieved very high AUCs ($\\gtrsim 0.9$), suggesting that the underlying generative process could be consistent with CoxPH's simple assumption of log-linearity and proportionality. Then, the semi-parametric model, CoxPH, is likely to achieve equal or better performances than the nonparametric model, SurvPP, which is consistent with the result. \n\nTo test the proposed method on a more non-linear structure, we ran an additional experiment on a real-world dataset, *Standard And New Antiepileptic Drugs study data* (SANAD), provided by R package joineR (GPL-3). SANAD was an unblinded randomized trial recruiting patients with epilepsy for whom carbamazepine (CBZ) was considered to be standard treatment and they were randomized to CBZ or the newer drug lamotrigine (LTG), where 605 patients were included and event was the time to treatment failure. We adopted a time-varying covariate or calibrated dose, and three static covariates including age of patient at randomization, gender, and randomized treatment (CBZ or LTG). The way of 10-fold cross evaluation of performance follows the experiment on PBC data. Figure R2 plots the result, showing that SurvPP achieved better performance than CoxPH, and achieved comparable performance while being substantially faster than GBM.\n\nWe will add the two real-world experiments in the main paper."}, {"review_id": "lrjjEsQh79", "reviewer": "Reviewer_feqp", "summary": "The paper proposes the first coherent probabilistic formulation for survival analysis in scenarios where the covariates depend on time, present all practical details for learning MAP estimates and marginal likelihoods and predictive distributions, and provide a range of synthetic and real-data experiments that demonstrate the approach is superior to earlier heuristic solutions and standard survival models. The limitations are adequately addressed.", "questions": "1) What would be the most interesting motivational example of importance of time-varying context that you could describe in 4-5 lines, with detailed examples of what the time-varying covariates could be and what goes wrong if using standard survival models? Could you fit that in your response and then in the Introduction?", "limitations": "None", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The paper is overall strong, providing comprehensive solution to an open problem and delivering it in a clear manner with sufficient empirical validation. I cannot immediately think of any specific way the paper could easily be improved, except perhaps by adding a clear motivational example for readers who cannot immediately understand the setting. It is by no means a perfect paper, but nevertheless represents solid scientific work in the field.\n\nThe problem is highly relevant and rarely studied. Even though some previous attempts have been made, we have been missing a proper formulation and I think that this paper provides one in a form that other can build on. The solution proposed here is definitely sound and covered well in the paper. The technical solution is not extremely novel in the sense that the hazard function is simply a transformation of a GP and standard MAP estimation is used, but this is appropriate and perhaps even ideal for a paper that proposes the formal definition of the process. It makes the paper easier to read for a broad audience, compared to one that would attempt to also use novel methodological elements within the process. The paper as written not encourages further work in improved modelling of the hazard functions and other inference techniques, and some ideas are already outlined in the conclusion.", "weaknesses": "The main empirical result is that we do not really observe an improvement in accuracy over the previous methods, which is naturally a bit of a disappointment. However, the proposed method is faster.\n\nAnother minor weakness is evaluation only on a fairly simple artificial data. A brief demonstration on the value of properly addressing time-varying covariates in some easy-to-understand application task would have made the paper more interesting also for audiences outside the core NeurIPS readers. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for the highly positive comments and suggestions, by which we are strongly encouraged!\n\n**What would be the most interesting motivational example of importance of time-varying context that you could describe in 4-5 lines, with detailed examples of what the time-varying covariates could be and what goes wrong if using standard survival models? Could you fit that in your response and then in the Introduction?**\n\nAccording to the reviewer's suggestion, we will fit the following explanation in the Introduction. We hope that our explanation is satisfactory. \n\nWhen exogenous and controllable covariates are considered, survival models enable us to estimate a hazard function as an explicit function of the covariates, by which we can perform *what-if analysis*, that is, simulations of the future failure events under various covariate function over time. Then we can find the optimal covariate function or policy that would yield a desirable future state, which is an important application of survival analysis with controllable covariates. Time-varying covariates are common in survival data, and if survival models for static covariates are applied to such a non-stationary data, then the static survival models fail to estimate the underlying dependence of hazard function on covariates, resulting in unreliable decision-makings. Therefore, survival models that accommodate time-varying covariates are needed. For example, in reliability engineering applications where the events are the failure times of machines, time-varying covariates could be the room temperature/humidity, which is controlled by air conditioners, and the scheduled maintenance. The maintenance manager can optimize the schedules of maintenance and air conditioning by balancing the risk of failure and the costs of maintenance and air conditioning. 　\n"}]

Author rebuttal:
We would like to thank all reviewers for valuable comments. We provided responses in as much detail as possible. We believe that we can fully dispel the concerns raised by the reviewers. Please see the submitted pdf file which includes newly added figures (Figure R1, R2, R3, R4) as responses to the reviewers' suggestions.

Meta-review: This well-written paper has been assessed by five knowledgeable reviewers who varied in their opinions: two recommended straight acceptance, one borderline acceptance, one borderline rejection, and one straight rejection. The work proposes new nonparametric hazard estimator for survival analysis, supported by empirical studies and theoretical discussion. Most of the concerns raised by the reviewers have been addressed by the authors' rebuttals, however the two reviewers who issued the lowest scores did not participate in the discussion. Overall assessment is that this work is sufficiently novel to attract attention of the NeurIPS audience and should be accepted. The authors are urged to reflect the key points raised in the discussion with the reviewers in the final version of their manuscript.
Decision: Accept (poster)