# Neural Frailty Machine: Beyond proportional hazard assumption in neural
survival regressions

Ruofan Wu∗ Ant Group  Jiawei Qiao Equal contribution Fudan University  Mingzhe
Wu  Coupang
Wen Yu  Fudan University  Ming Zheng  Fudan University  Tengfei Liu  Ant Group
Tianyi Zhang  Ant Group  Weiqiang Wang  Ant Group

###### Abstract

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

##  1 Introduction

Regression analysis of time-to-event data Kalbfleisch and Prentice (2002) has
been among the most important modeling tools for clinical studies and has
witnessed a growing interest in areas like corporate finance Duffie et al.
(2009), recommendation systems Jing and Smola (2017), and computational
advertising Wu et al. (2015). The key feature that differentiates time-to-
event data from other types of data is that they are often _incompletely
observed_ , with the most prevailing form of incompleteness being the _right
censoring_ mechanism Kalbfleisch and Prentice (2002). In the right censoring
mechanism, the duration time of a sampled subject is (sometimes) only known to
be larger than the observation time instead of being recorded precisely. It is
well known in the community of survival analysis that even in the case of
linear regression, naively discarding the censored observations produces
estimation results that are statistically biased Buckley and James (1979), at
the same time losses sample efficiency if the censoring proportion is high.

Cox’s proportional hazard (CoxPH ) model Cox (1972) using the convex objective
of negative partial likelihood Cox (1975) is the _de facto_ choice in modeling
right censored time-to-event data (hereafter abbreviated as censored data
without misunderstandings). The model is _semiparametric_ Bickel et al. (1993)
in the sense that the baseline hazard function needs no parametric
assumptions. The original formulation of CoxPH model assumes a linear form and
therefore has limited flexibility since the truth is not necessarily linear.
Subsequent studies extended CoxPH model to nonlinear variants using ideas from
nonparametric regression Huang (1999); Cai et al. (2007, 2008), ensemble
learning Ishwaran et al. (2008), and neural networks Faraggi and Simon (1995);
Katzman et al. (2018). While such extensions allowed a more flexible nonlinear
dependence structure with the covariates, the learning objectives were still
derived under the proportional hazards (PH) assumption, which was shown to be
inadequate in many real-world scenarios Gray (2000). The most notable case was
the failure of modeling the phenomenon of crossing hazards Stablein and
Koutrouvelis (1985). It is thus of significant interest to explore extensions
of CoxPH that both allow nonlinear dependence over covariates and relaxations
of the PH assumption.

Frailty models Wienke (2010); Duchateau and Janssen (2007) are among the most
important research topics in modern survival analysis, in that they provide a
principled way of extending CoxPH model via incorporating a multiplicative
random effect to capture unobserved heterogeneity. The resulting
parameterization contains many useful variants of CoxPH like the proportional
odds model Bennett (1983), under specific choices of frailty families. While
the theory of frailty models has been well-established Murphy (1994, 1995);
Parner (1998); Kosorok et al. (2004), most of them focused on the linear case.
Recent developments on applying neural approaches to survival analysis Katzman
et al. (2018); Kvamme et al. (2019); Tang et al. (2022); Rindt et al. (2022)
have shown promising results in terms of empirical predictive performance,
with most of them lacking theoretical discussions. Therefore, it is of
significant interest to build more powerful frailty models via adopting
techniques in modern deep learning Goodfellow et al. (2016) with provable
statistical guarantees.

In this paper, we present a general framework for neural extensions of frailty
models called the neural frailty machine (NFM). Two concrete neural
architectures are derived under the framework: The first one adopts the
proportional frailty assumption, allowing an intuitive interpretation of the
neural CoxPH model with a multiplicative random effect. The second one further
relaxes the proportional frailty assumption and could be viewed as an
extension of nonparametric hazard regression (NHR) Cox and O’Sullivan (1990);
Kooperberg et al. (1995), sometimes referred to as "fully neural" models under
the context of neural survival analysis Omi et al. (2019). We summarize our
contributions as follows.

  * •

We propose the neural frailty machine (NFM) framework as a principled way of
incorporating unobserved heterogeneity into neural survival regression models.
The framework includes many commonly used survival regression models as
special cases.

  * •

We derive two model architectures based on the NFM framework that extend
neural CoxPH models and neural NHR models. Both models allow stochastic
training and scale to large datasets.

  * •

We show theoretical guarantees for the two proposed models via characterizing
the rates of convergence of the proposed nonparametric function estimators.
The proof technique is different from previous theoretical studies on neural
survival analysis and is applicable to many other types of neural survival
models.

  * •

We conduct extensive studies on various benchmark datasets at different
scales. Under standard performance metrics, both models are empirically shown
to perform competitively, matching or outperforming state-of-the-art neural
survival models.

##  2 Related works

###  2.1 Nonlinear extensions of CoxPH

Most nonlinear extensions of CoxPH model stem from the equivalence of partial
likelihood and semiparametric profile likelihood Murphy and Van der Vaart
(2000) of CoxPH model, resulting in nonlinear variants that essentially
replaces the linear term in partial likelihood with nonlinear variants: Huang
(1999) used smoothing splines, Cai et al. (2007, 2008) used local polynomial
regression Fan and Gijbels (1996). The empirical success of tree-based models
inspired subsequent developments like Ishwaran et al. (2008) that equip tree-
based models such as gradient boosting trees and random forests with losses in
the form of negative log partial likelihood. Early developments of neural
survival analysis Faraggi and Simon (1995) adopted similar extension
strategies and obtained neural versions of partial likelihood. Later attempts
Katzman et al. (2018) suggest using the successful practice of stochastic
training which is believed to be at the heart of the empirical success of
modern neural methods Hardt et al. (2016). However, stochastic training under
the partial likelihood objective is highly non-trivial, as mini-batch versions
of log partial likelihood Katzman et al. (2018) are no longer valid stochastic
gradients of the full-sample log partial likelihood Tang et al. (2022).

###  2.2 Beyond CoxPH in survival analysis

In linear survival modeling, there are standard alternatives to CoxPH such as
the accelerated failure time (AFT) model Buckley and James (1979); Ying
(1993), the extended hazard regression model Etezadi-Amoli and Ciampi (1987),
and the family of linear transformation models Zeng and Lin (2006). While
these models allow certain types of nonlinear extensions, the resulting form
of (conditional) hazard function is still restricted to be of a specific form.
The idea of nonparametric hazard regression (NHR) Cox and O’Sullivan (1990);
Kooperberg et al. (1995); Strawderman and Tsiatis (1996) further improves the
flexibility of nonparametric survival analysis via directly modeling the
conditional hazard function by nonparametric regression techniques such as
spline approximation. Neural versions of NHR have been developed lately such
as the CoxTime model Kvamme et al. (2019). Rindt et al. (2022) used a neural
network to approximate the conditional survival function and could be thus
viewed as another trivial extension of NHR.

Aside from developments in NHR, Lee et al. (2018) proposed a discrete-time
model with its objective being a mix of the discrete likelihood and a rank-
based score; Zhong et al. (2021a) proposed a neural version of the extended
hazard model, unifying both neural CoxPH and neural AFT model; Tang et al.
(2022) used an ODE approach to model the hazard and cumulative hazard
functions.

###  2.3 Theoretical justification of neural survival models

Despite the abundance of neural survival models, assessment of their
theoretical properties remains nascent. In Zhong et al. (2021b), the authors
developed minimax theories of partially linear cox model using neural networks
as the functional approximator. Zhong et al. (2021a) provided convergence
guarantees of neural estimates under the extended hazard model. The
theoretical developments therein rely on specific forms of objective function
(partial likelihood and kernel pseudo-likelihood) and are not directly
applicable to the standard likelihood-based objective which is frequently used
in survival analysis.

##  3 Methodology

###  3.1 The neural frailty machine framework

Let $\tilde{T}\geq 0$ be the interested event time with survival function
denoted by $S(t)=\mathbb{P}(\tilde{T}>t)$ associated with a feature(covariate)
vector $Z\in\mathbb{R}^{d}$. Suppose that $\tilde{T}$ is a continuous random
variable and let $f(t)$ be its density function. Then $\lambda(t)=f(t)/S(t)$
is the hazard function and $\Lambda(t)=\int_{0}^{t}\lambda(s)ds$ is the
cumulative hazard function. Aside from the covariate $Z$, we use a positive
scalar random variable $\omega\in\mathbb{R}^{+}$ to express the unobserved
heterogeneity corresponding to individuals, or _frailty_. 111For example in
medical biology, it was observed that genetically identical animals kept in as
similar an environment as possible will typically not behave the same upon
exposure to environmental carcinogens Brennan (2002). In this paper we will
assume the following generating scheme of $\tilde{T}$ via specifying its
conditional hazard function:

$\displaystyle\lambda(t|Z,\omega)=\omega\widetilde{\nu}(t,Z).$ (1)

Here $\widetilde{\nu}$ is an unspecified non-negative function, and we let the
distribution of $\omega$ be parameterized by a one-dimensional parameter
$\theta\in\mathbb{R}$. 222The choice of one-dimensional frailty family is
mostly for simplicity and clearness of theoretical derivations. Note that
there exist multi-dimensional frailty families like the PVF family Wienke
(2010). Generalizing our theoretical results to such kinds of families would
require additional sets of regularity conditions, and will be left to future
explorations. The formulation (1) is quite general and contains several
important models in both traditional and neural survival analysis:

  1. 1.

When $\omega$ follows parametric distributional assumptions, and
$\widetilde{\nu}(t,Z)=\lambda(t)e^{\beta^{\top}Z}$, (1) reduces to the
standard proportional frailty model Kosorok et al. (2004). A special case is
when $\omega$ is degenerate, i.e., it has no randomness, then the model
corresponds to the classic CoxPH model.

  2. 2.

When $\omega$ is degenerate and $\widetilde{\nu}$ is arbitrary, the model
becomes equivalent to nonparametric hazard regression (NHR) Cox and O’Sullivan
(1990); Kooperberg et al. (1995). In NHR, the function parameter of interest
is usually the logarithm of the (conditional) hazard function.

In this paper we construct neural approximations to the logarithm of
$\widetilde{\nu}$, i.e., $\nu(t,Z)=\log\widetilde{\nu}(t,Z)$. The resulting
models are called Neural Frailty Machines (NFM). Depending on the prior
knowledge of the function $\nu$, we propose two function approximation
schemes:

The proportional frailty (PF) scheme assumes the dependence of $\nu$ on event
time and covariates to be completely _decoupled_ , i.e.,

$\displaystyle\nu(t,Z)=h(t)+m(Z).$ (2)

Proportional-style assumption over hazard functions has been shown to be a
useful inductive bias in survival analysis. We will treat both $h$ and $m$ in
(2) as function parameters, and device two multi-layer perceptrons (MLP) to
approximate them separately.

The fully neural (FN) scheme imposes no a priori assumptions over $\nu$ and is
the most general version of NFM. It is straightforward to see that the most
commonly used survival models, such as CoxPH , AFTYing (1993), EHZhong et al.
(2021a), or PF models are included in the proposed model space as special
cases. We treat $\nu=\nu(t,Z)$ as the function parameter with input dimension
$d+1$ and use a multi-layer perceptron (MLP) as the function approximator to
$\nu$. Similar approximation schemes with respect to the hazard function have
been proposed in some recent works Omi et al. (2019); Rindt et al. (2022),
referred to as "fully neural approaches" without theoretical
characterizations.

The choice of frailty family There are many commonly used families of frailty
distributions Kosorok et al. (2004); Duchateau and Janssen (2007); Wienke
(2010), among which the most popular one is the _gamma frailty_ , where
$\omega$ follows a gamma distribution with mean $1$ and variance $\theta$. We
briefly introduce some other types of frailty families in appendix A.

###  3.2 Parameter learning under censored observations

In time-to-event modeling scenarios, the event times are typically observed
under right censoring. Let $C$ be the right censoring time which is assumed to
be conditionally independent of the event time $\tilde{T}$ given $Z$, i.e.,
$\tilde{T}\perp\\!\\!\\!\\!\perp C|Z$. In data collection, one can observe the
minimum of the survival time and the censoring time, that is, observe
$T=\tilde{T}\wedge C$ as well as the censoring indicator
$\delta=I(\tilde{T}\leqslant C)$, where $a\wedge b=\min(a,b)$ for constants
$a$ and $b$ and $I(\cdot)$ stands for the indicator function. We assume $n$
independent and identically distributed (i.i.d.) copies of $(T,\delta,Z)$ are
used as the training sample $(T_{i},\delta_{i},Z_{i}),i\in[n]$, where we use
$[n]$ to denote the set $\\{1,2,\ldots,n\\}$. Additionally, we assume the
unobserved frailties are independent and identically distributed, i.e.,
$\omega_{i}\overset{\text{i.i.d.}}{\sim}f_{\theta}(\omega),i\in[n]$. Next, we
derive the learning procedure based on the observed log-likelihood (OLL)
objective under both PF and FN scheme. To obtain the observed likelihood, we
first integrate the conditional survival function given the frailty:

$\displaystyle\begin{aligned} S(t|Z)&=\mathbb{E}_{\omega\sim
f_{\theta}}\left[e^{-\omega\int_{0}^{t}e^{\nu(s,Z)}ds}\right]\\\
&=:e^{-G_{\theta}\left(\int_{0}^{t}e^{\nu(s,Z)}ds\right)}.\end{aligned}$ (3)

Here the _frailty transform_ $G_{\theta}(x)=-\log\left(\mathbb{E}_{\omega\sim
f_{\theta}}\left[e^{-\omega x}\right]\right)$ is defined as the negative of
the logarithm of the Laplace transform of the frailty distribution. The
conditional cumulative hazard function is thus
$\Lambda(t|Z)=G_{\theta}(\int_{0}^{t}e^{\nu(s,Z)}ds)$. For the PF scheme of
NFM, we use two MLPs
$\widehat{h}=\widehat{h}(t;\mathbf{W}^{h},\mathbf{b}^{h})$ and
$\widehat{m}=\widehat{m}(Z;\mathbf{W}^{m},\mathbf{b}^{m})$ as function
approximators to $\nu$ and $m$, parameterized by
$(\mathbf{W}^{h},\mathbf{b}^{h})$ and $(\mathbf{W}^{m},\mathbf{b}^{m})$,
respectively. 333Here we adopt the conventional notation that $\mathbf{W}$ is
the collection of the weight matrices of the MLP in all layers, and
$\mathbf{b}$ corresponds to the collection of the bias vectors in all layers.
According to standard results on censored data likelihood Kalbfleisch and
Prentice (2002), we write the learning objective under the PF scheme as:

$\displaystyle\begin{aligned}
&\mathcal{L}(\mathbf{W}^{h},\mathbf{b}^{h},\mathbf{W}^{m},\mathbf{b}^{m},\theta)\\\
=&\frac{1}{n}\left[\sum_{i\in[n]}\delta_{i}\log
g_{\theta}\left(e^{\widehat{m}(Z_{i})}\int_{0}^{T_{i}}e^{\widehat{h}(s)}ds\right)+\delta_{i}\widehat{h}(T_{i})+\delta_{i}\widehat{m}(Z_{i})-G_{\theta}\left(e^{\widehat{m}(Z_{i})}\int_{0}^{T_{i}}e^{\widehat{h}(s)}ds\right)\right].\end{aligned}$
(4)

Here we define $g_{\theta}(x)=\frac{\partial}{\partial x}G_{\theta}(x)$. Let
$(\widehat{\mathbf{W}}^{h}_{n},\widehat{\mathbf{b}}^{h}_{n},\widehat{\mathbf{W}}^{m}_{n},\widehat{\mathbf{b}}^{m}_{n},\widehat{\theta}_{n})$
be the maximizer of (4) and further denote
$\widehat{h}_{n}(t)=\widehat{h}(t;\widehat{\mathbf{W}}^{h}_{n},\widehat{\mathbf{b}}^{h}_{n})$
and
$\widehat{m}_{n}(Z)=\widehat{m}(Z;\widehat{\mathbf{W}}^{m}_{n},\widehat{\mathbf{b}}^{m}_{n})$.
The resulting estimators for conditional cumulative hazard and survival
functions are:

$\displaystyle\begin{aligned}
&\widehat{\Lambda}_{{\text{{PF}}}}(t|Z)=G_{\widehat{\theta}_{n}}\left(\int_{0}^{t}e^{\widehat{h}_{n}(s)+\widehat{m}_{n}(Z)}ds\right),\\\
&\widehat{S}_{{\text{{PF}}}}(t|Z)=e^{-\widehat{\Lambda}_{{\text{{PF}}}}(t|Z)},\end{aligned}$
(5)

For the FN scheme, we use
$\widehat{\nu}=\widehat{\nu}(t,Z;\mathbf{W}^{\nu},\mathbf{b}^{\nu})$ to
approximate $\nu(t,Z)$ parameterized by $(\mathbf{W}^{\nu},\mathbf{b}^{\nu})$.
The OLL objective is written as:

$\displaystyle\begin{aligned}
&\mathcal{L}(\mathbf{W}^{\nu},\mathbf{b}^{\nu},\theta)\\\
=&\frac{1}{n}\left[\sum_{i\in[n]}\delta_{i}\log
g_{\theta}\left(\int_{0}^{T_{i}}e^{\widehat{\nu}(s,Z_{i};\mathbf{W}^{\nu},\mathbf{b}^{\nu})}ds\right)+\delta_{i}\widehat{\nu}(T_{i},Z_{i};\mathbf{W}^{\nu},\mathbf{b}^{\nu})-G_{\theta}\left(\int_{0}^{T_{i}}e^{\widehat{\nu}(s,Z_{i};\mathbf{W}^{\nu},\mathbf{b}^{\nu})}ds\right)\right].\end{aligned}$
(6)

Let
$(\widehat{\mathbf{W}}^{\nu}_{n},\widehat{\mathbf{b}}^{\nu}_{n},\widehat{\theta}_{n})$
be the maximizer of (6), and further denote
$\widehat{\nu}_{n}(t,Z)=\widehat{\nu}(t,Z;\widehat{\mathbf{W}}^{\nu}_{n},\widehat{\mathbf{b}}^{\nu}_{n})$.
The conditional cumulative hazard and survival functions are therefore
estimated as:

$\displaystyle\begin{aligned}
&\widehat{\Lambda}_{{\text{{FN}}}}(t|Z)=G_{\widehat{\theta}_{n}}\left(\int_{0}^{t}e^{\widehat{\nu}_{n}(s,Z)}ds\right),\\\
&\widehat{S}_{{\text{{FN}}}}(t|Z)=e^{-\widehat{\Lambda}_{{\text{{FN}}}}(t|Z)}.\end{aligned}$
(7)

The evaluation of objectives like (6) and its gradient requires computing a
definite integral of an exponentially transformed MLP function. Instead of
using exact computations that are available for only a restricted type of
activation functions and network structures, we use numerical integration for
such kinds of evaluations, using the method of Clenshaw-Curtis quadrature Boyd
(2001), which has shown competitive performance and efficiency in recent
applications to monotonic neural networks Wehenkel and Louppe (2019).

######  Remark 3.1.

The interpretation of frailty terms differs in the two schemes. In the PF
scheme, introducing the frailty effect strictly increases the modeling
capability (i.e., the capability of modeling crossing hazard) in comparison to
CoxPH or neural variants of CoxPH Kosorok et al. (2004). In the FN scheme, it
is arguable that in the i.i.d. case, the marginal hazard function is a
reparameterization of the hazard function in the context of NHR. Therefore, we
view the incorporation of frailty effect as injecting a domain-specific
inductive bias that has proven to be useful in survival analysis and time-to-
event regression modeling and verify this claim empirically in section 5.2.
Moreover, frailty becomes especially helpful when handling correlated or
clustered data where the frailty term is assumed to be shared among certain
groups of individuals Parner (1998). Extending NFM to such scenarios is
valuable and we left it to future explorations.

##  4 Theoretical results

In this section, we present theoretical properties of both NFM estimates by
characterizing their rates of convergence when the underlying event data
follows corresponding model assumptions. The proof technique is based on the
method of sieves Shen and Wong (1994); Shen (1997); Chen (2007) that views
neural networks as a special kind of nonlinear sieve Chen (2007) that
satisfies desirable approximation properties Yarotsky (2017). Since both
models produce estimates of function parameters, we need to specify a suitable
function space to work with. Here we choose the following Hölder ball as was
also used in previous works on nonparametric estimation using neural networks
Schmidt-Hieber (2020); Farrell et al. (2021); Zhong et al. (2021b)

$\displaystyle\mathcal{W}^{\beta}_{M}(\mathcal{X})=\left\\{f:\max_{\mathbf{\alpha}:\left|\mathbf{\alpha}\right|\leq\beta}\underset{x\in\mathcal{X}}{\text{esssup}}\left|D^{\mathbf{\alpha}}(f(x))\right|\leq
M\right\\},$ (8)

where the domain $\mathcal{X}$ is assumed to be a subset of $d$-dimensional
euclidean space. $\mathbf{\alpha}=(\alpha_{1},\ldots,\alpha_{d})$ is a
$d$-dimensional tuple of nonnegative integers satisfying
$\left|\mathbf{\alpha}\right|=\alpha_{1}+\cdots+\alpha_{d}$ and
$D^{\mathbf{\alpha}}f=\frac{\partial^{\left|\mathbf{\alpha}\right|}f}{\partial
x_{1}^{\alpha_{1}}\cdots x_{d}^{\alpha_{d}}}$ is the weak derivative of $f$.
Now assume that $M$ is a reasonably large constant, and let $\Theta$ be a
closed interval over the real line. We make the following assumptions for the
_true parameters_ under both schemes:

######  Condition 4.1 (True parameter, PF scheme).

The euclidean parameter $\theta_{0}\in\Theta\subset\mathbb{R}$, and the two
function parameters
$m_{0}\in\mathcal{W}^{\beta}_{M}([-1,1]^{d}),h_{0}\in\mathcal{W}^{\beta}_{M}([0,\tau])$,
and $\tau>0$ is the ending time of the study duration, which is usually
adopted in the theoretical studies in survival analysis Van der Vaart (2000).

######  Condition 4.2 (True parameter, FN scheme).

The euclidean parameter $\theta_{0}\in\Theta\subset\mathbb{R}$, and the
function parameter
$\nu_{0}\in\mathcal{W}^{\beta}_{M}([0,\tau]\times[-1,1]^{d})$,

Next, we construct sieve spaces for function parameter approximation via
restricting the complexity of the MLPs to "scale" with the sample size $n$.

######  Condition 4.3 (Sieve space, PF scheme).

The sieve space $\mathcal{H}_{n}$ is constructed as a set of MLPs satisfying
$\widehat{h}\in\mathcal{W}^{\beta}_{M_{h}}([0,\tau])$, with depth of order
$O(\log n)$ and total number of parameters of order
$O(n^{\frac{1}{\beta+d}}\log n)$. The sieve space $\mathcal{M}_{n}$ is
constructed as a set of MLPs satisfying
$\widehat{m}\in\mathcal{W}^{\beta}_{M_{m}}([-1,1]^{d})$, with depth of order
$O(\log n)$ and total number of parameters of order
$O(n^{\frac{d}{\beta+d}}\log n)$. Here $M_{h}$ and $M_{m}$ are sufficiently
large constants such that every function in
$\mathcal{W}^{\beta}_{M}([-1,1]^{d})$ and $\mathcal{W}^{\beta}_{M}([0,\tau])$
could be accurately approximated by functions inside $\mathcal{H}_{n}$ and
$\mathcal{M}_{n}$, according to (Yarotsky, 2017, Theorem 1).

######  Condition 4.4 (Sieve space, FN scheme).

The sieve space $\mathcal{V}_{n}$ is constructed as a set of MLPs satisfying
$\widehat{\nu}\in\mathcal{W}^{\beta}_{M_{\nu}}([0,\tau])$, with depth of order
$O(\log n)$ and total number of parameters of order
$O(n^{\frac{d+1}{\beta+d+1}}\log n)$. Here $M_{\nu}$ is a sufficiently large
constant such that $\mathcal{V}_{n}$ satisfies approximation properties,
analogous to condition 4.3.

For technical reasons, we will assume the nonparametric function estimators
are constrained to fall inside the corresponding sieve spaces, i.e.,
$\widehat{h}_{n}\in\mathcal{H}_{n}$, $\widehat{m}_{n}\in\mathcal{M}_{n}$ and
$\widehat{\nu}\in\mathcal{V}_{n}$. This will not affect the implementation of
optimization routines as was discussed in Farrell et al. (2021). Furthermore,
we restrict the estimate $\widehat{\theta}_{n}\in\Theta$ in both PF and FN
schemes.

Additionally, we need the following regularity condition on the function
$G_{\theta}(x)$:

######  Condition 4.5.

$G_{\theta}(x)$ is viewed as a bivariate function
$G:\Theta\times\mathcal{B}\mapsto\mathbb{R}$, where $\mathcal{B}$ is a compact
set on $\mathbb{R}$. The functions
$G_{\theta}(x)$,$\frac{\partial}{\partial\theta}G_{\theta}(x)$,$\frac{\partial}{\partial
x}G_{\theta}(x)$,$\log g_{\theta}(x)$,$\frac{\partial}{\partial\theta}\log
g_{\theta}(x)$, $\frac{\partial}{\partial x}\log g_{\theta}(x)$ are bounded on
$\Theta\times\mathcal{B}$.

We define two metrics that measures convergence of parameter estimates: For
the PF scheme, let $\phi_{0}=(h_{0},m_{0},\theta_{0})$ be the true parameters
and
$\widehat{\phi}_{n}=(\widehat{h}_{n},\widehat{m}_{n},\widehat{\theta}_{n})$ be
the estimates. We abbreviate $\mathbb{P}_{\phi_{0},Z=z}$ as the conditional
probability distribution of $(T,\delta)$ given $Z=z$ under the true parameter,
and $\mathbb{P}_{\widehat{\phi}_{n},Z=z}$ as the conditional probability
distribution of $(T,\delta)$ given $Z=z$ under the estimates. Define the
following metric

$\displaystyle
d_{{\text{{PF}}}}\left(\widehat{\phi}_{n},\phi_{0}\right)=\sqrt{\mathbb{E}_{z\sim\mathbb{P}_{Z}}\left[H^{2}(\mathbb{P}_{\widehat{\phi}_{n},Z=z}\parallel\mathbb{P}_{\phi_{0},Z=z})\right]},$
(9)

where
$H^{2}(\mathbb{P}\parallel\mathbb{Q})=\int\left(\sqrt{d\mathbb{P}}-\sqrt{d\mathbb{Q}}\right)^{2}$
is the squared Hellinger distance between probability distributions
$\mathbb{P}$ and $\mathbb{Q}$. The case for the FN scheme is similar: Let
$\psi_{0}=(\nu_{0},\theta_{0})$ be the parameters and
$\widehat{\nu}_{n}=(\widehat{\nu}_{n},\widehat{\theta}_{n})$ be the estimates.
Analogous to the definitions above, we define $\mathbb{P}_{\psi_{0},Z=z}$ as
the true conditional distribution given $Z=z$, and
$\mathbb{P}_{\widehat{\psi}_{n},Z=z}$ be the estimated conditional
distribution, we will use the following metric in the FN scheme:

$\displaystyle
d_{{\text{{FN}}}}\left(\widehat{\psi}_{n},\psi_{0}\right)=\sqrt{\mathbb{E}_{z\sim\mathbb{P}_{Z}}\left[H^{2}(\mathbb{P}_{\widehat{\psi}_{n},Z=z}\parallel\mathbb{P}_{\psi_{0},Z=z})\right]}.$
(10)

Now we state our main theorems. We denote $\mathbb{P}$ as the data generating
distribution and use $\widetilde{O}$ to hide poly-logarithmic factors in the
big-O notation.

######  Theorem 4.6 (Rate of convergence, PF scheme).

In the PF scheme, under condition 4.1, 4.3, 4.5, we have that
$d_{{\text{{PF}}}}\left(\widehat{\phi}_{n},\phi_{0}\right)=\widetilde{O}_{\mathbb{P}}\left(n^{-\frac{\beta}{2\beta+2d}}\right)$.

######  Theorem 4.7 (Rate of convergence, FN scheme).

In the FN scheme, under condition 4.2, 4.4, 4.5, we have that
$d_{{\text{{FN}}}}\left(\widehat{\psi}_{n},\psi_{0}\right)=\widetilde{O}_{\mathbb{P}}\left(n^{-\frac{\beta}{2\beta+2d+2}}\right)$.

######  Remark 4.8.

The idea of using Hellinger distance to measure the convergence rate of sieve
MLEs was proposed in Wong and Shen (1995). Obtaining rates under a stronger
topology such as $L_{2}$ is possible if the likelihood function satisfies
certain conditions such as the curvature condition Farrell et al. (2021).
However, such kind of conditions are in general too stringent for likelihood-
based objectives, instead, we use Hellinger convergence that has minimal
requirements. Consequently, our proof strategy is applicable to many other
survival models that rely on neural function approximation such as Rindt et
al. (2022), with some modification to the regularity conditions. For proper
choices of metrics in sieve theory, see also the discussion in (Chen, 2007,
Chapter 2).

##  5 Experiments

In this section, we assess the empirical performance of NFM. We first conduct
synthetic experiments for verifying the theoretical convergence guarantees
developed in section 4. To further illustrate the empirical efficacy of NFM,
we evaluate the predictive performance of NFM over $6$ benchmark datasets
ranging from small scale to large scale, against state-of-the-art baselines.

###  5.1 Synthetic experiments

We conduct synthetic experiments to validate our proposed theory. The
underlying data generating scheme is as follows: First, we generate a
$5$-dimensional feature $Z$ that is independently sampled from the uniform
distribution over the interval $[0,1]$. The (true) conditional hazard function
of the event time takes the form of the proportional frailty model (2), with
$h(t)=t$ and $m(Z)=\sin(\langle Z,\beta\rangle)+\langle\sin(Z),\beta\rangle$,
where $\beta=(0.1,0.2,0.3,0.4,0.5)$. The frailty $\omega$ is generated
according to a gamma distribution with mean and variance equal to 1. We use
this generating model to assess the recovery guarantee of both NFM modeling
schemes via inspecting the empirical recovery of $\nu(t,Z)$. For the PF
scheme, we have more underlying information about the generating model, and we
present an additional assessment regarding the recovery of $m(Z)$ in appendix
D.1. We generate three training datasets of different scales, with
$n\in\\{1000,5000,10000\\}$. A censoring mechanism is applied such that the
censoring ratio is around $40\%$ for each dataset. The assessment will be made
on a fixed test sample of $100$ hold-out points that are independently drawn
from the generating scheme of the event time. We report a more detailed
description of the implementation of the data generating scheme and model
architectures in appendix C.2. We present the results of our synthetic data
experiments in figure 1. The evaluation results suggest that both NFM schemes
are capable of approximating complicated nonlinear functions using a moderate
amount of data, i.e., $n\geq 1000$.

$0$$50$$100$$0$$2$$4$index$\nu(t,Z)$$N=1000\ \text{(PF)}$ground truthestimated

$0$$50$$100$$0$$2$$4$index$\nu(t,Z)$$N=5000\ \text{(PF)}$ground truthestimated

$0$$50$$100$$0$$2$$4$index$\nu(t,Z)$$N=10000\ \text{(PF)}$ground
truthestimated

$0$$50$$100$$0$$2$$4$index$\nu(t,Z)$$N=1000\ \text{(FN)}$ground truthestimated

$0$$50$$100$$0$$2$$4$index$\nu(t,Z)$$N=5000\ \text{(FN)}$ground truthestimated

$0$$50$$100$$0$$2$$4$index$\nu(t,Z)$$N=10000\ \text{(FN)}$ground
truthestimated

Figure 1: Visualizations of synthetic data results under the NFM framework.
The plots in the first row compare the empirical estimates of the
nonparametric component $\nu(t,Z)$ against its true value evaluated on $100$
hold-out points, under the PF scheme. The plots in the second row are obtained
using the FN scheme, with analogous semantics to the first row.

###  5.2 Real-world data experiments

Datasets We use five survival datasets and one non-survival dataset for
evaluation. The survival datasets include the Molecular Taxonomy of Breast
Cancer International Consortium (METABRIC) Curtis et al. (2012), the Rotterdam
tumor bank and German Breast Cancer Study Group (RotGBSG)Knaus et al. (1995),
the Assay Of Serum Free Light Chain (FLCHAIN) Dispenzieri et al. (2012), the
Study to Understand Prognoses Preferences Outcomes and Risks of Treatment
(SUPPORT) Knaus et al. (1995), and the Medical Information Mart for Intensive
Care (MIMIC-III) Johnson et al. (2016). For all the survival datasets, the
event of interest is defined as the mortality after admission. In our
experiments, we view METABRIC, RotGBSG, FLCHAIN, and SUPPORT as small-scale
datasets and MIMIC-III as a moderate-scale dataset. We additionally use the
KKBOX dataset Kvamme et al. (2019) as a large-scale evaluation. In this
dataset, an event time is observed if a customer churns from the KKBOX
platform. We summarize the basic statistics of all the datasets in table 3.

Baselines We compare NFM with $12$ baselines. The first one is the linear
CoxPH model Cox (1972). Gradient Boosting Machine (GBM) Friedman (2001); Chen
and Guestrin (2016) and Random Survival Forests (RSF) Ishwaran et al. (2008)
are two tree-based nonparametric survival regression methods. DeepSurv Katzman
et al. (2018) and CoxTime Kvamme et al. (2019) are two models that adopt
neural variants of partial likelihood as objectives. SuMo-net Rindt et al.
(2022) is a neural variant of NHR. We additionally chose six latest state-of-
the-art neural survival models: DeepHit Lee et al. (2018), SurvNode Groha et
al. (2020), DeepEH Zhong et al. (2021a), DCM Nagpal et al. (2021), DeSurv
Danks and Yau (2022) and SODEN Tang et al. (2022). Among the chosen baselines,
DeepSurv and SuMo-net are viewed as implementations of neural CoxPH and neural
NHR and are therefore of particular interest for the empirical verification of
the efficacy of frailty.

Table 1: Survival prediction results measured in IBS and INBLL metric (%) on
four small-scale survival datasets. In each column, the boldfaced score
denotes the best result and the underlined score represents the second-best
result.

Model | METABRIC | RotGBSG | FLCHAIN | SUPPORT
---|---|---|---|---
| IBS | INBLL | IBS | INBLL | IBS | INBLL | IBS | INBLL
CoxPH | ${16.46}_{\pm 0.90}$ | ${49.57}_{\pm 2.66}$ | ${18.25}_{\pm 0.44}$ | ${53.76}_{\pm 1.11}$ | ${10.05}_{\pm 0.38}$ | ${33.18}_{\pm 1.16}$ | ${20.54}_{\pm 0.38}$ | ${59.58}_{\pm 0.86}$
GBM | ${16.61}_{\pm 0.82}$ | ${49.87}_{\pm 2.44}$ | ${17.83}_{\pm 0.44}$ | ${52.78}_{\pm 1.11}$ | $\underline{9.98}_{\pm 0.37}$ | $\underline{32.88}_{\pm 1.05}$ | ${19.18}_{\pm 0.39}$ | ${56.46}_{\pm 0.10}$
RSF | ${16.62}_{\pm 0.64}$ | ${49.61}_{\pm 1.54}$ | ${17.89}_{\pm 0.42}$ | ${52.77}_{\pm 1.01}$ | $\mathbf{9.96}_{\pm 0.37}$ | ${32.92}_{\pm 1.05}$ | $\underline{19.11}_{\pm 0.40}$ | $\underline{56.28}_{\pm 1.00}$
DeepSurv | ${16.55}_{\pm 0.93}$ | ${49.85}_{\pm 3.02}$ | ${17.80}_{\pm 0.49}$ | ${52.62}_{\pm 1.25}$ | ${10.09}_{\pm 0.38}$ | ${33.28}_{\pm 1.15}$ | ${19.20}_{\pm 0.41}$ | ${56.48}_{\pm 1.08}$
CoxTime | ${16.54}_{\pm 0.83}$ | ${49.67}_{\pm 2.67}$ | ${17.80}_{\pm 0.58}$ | ${52.56}_{\pm 1.47}$ | ${10.28}_{\pm 0.45}$ | ${34.18}_{\pm 1.53}$ | ${19.17}_{\pm 0.40}$ | ${56.45}_{\pm 1.10}$
DeepHit | ${17.50}_{\pm 0.83}$ | ${52.10}_{\pm 2.16}$ | ${19.61}_{\pm 0.38}$ | ${56.67}_{\pm 1.10}$ | ${11.83}_{\pm 0.39}$ | ${37.72}_{\pm 1.02}$ | ${20.66}_{\pm 0.32}$ | ${60.06}_{\pm 0.72}$
DeepEH | ${16.56}_{\pm 0.65}$ | ${49.42}_{\pm 1.53}$ | ${17.62}_{\pm 0.52}$ | $\underline{52.08}_{\pm 1.27}$ | ${10.11}_{\pm 0.37}$ | ${33.30}_{\pm 1.10}$ | ${19.30}_{\pm 0.39}$ | ${56.67}_{\pm 0.94}$
SuMo-net | ${16.49}_{\pm 0.83}$ | ${49.74}_{\pm 2.21}$ | ${17.77}_{\pm 0.47}$ | ${52.62}_{\pm 1.11}$ | ${10.07}_{\pm 0.40}$ | ${33.20}_{\pm 1.10}$ | ${19.40}_{\pm 0.38}$ | ${56.87}_{\pm 0.96}$
SODEN | ${16.52}_{\pm 0.63}$ | ${49.39}_{\pm 1.97}$ | $\mathbf{17.05}_{\pm 0.63}$ | $\mathbf{50.45}_{\pm 1.97}$ | ${10.13}_{\pm 0.24}$ | ${33.37}_{\pm 0.57}$ | ${19.07}_{\pm 0.50}$ | ${56.15}_{\pm 1.35}$
SurvNode | ${16.67}_{\pm 1.32}$ | ${49.73}_{\pm 3.89}$ | ${17.42}_{\pm 0.53}$ | ${51.70}_{\pm 1.16}$ | ${10.40}_{\pm 0.29}$ | ${34.37}_{\pm 1.03}$ | ${19.58}_{\pm 0.34}$ | ${57.49}_{\pm 0.84}$
DCM | ${16.58}_{\pm 0.87}$ | ${49.48}_{\pm 2.23}$ | ${17.66}_{\pm 0.54}$ | ${52.26}_{\pm 1.23}$ | ${10.13}_{\pm 0.50}$ | ${33.40}_{\pm 1.38}$ | ${19.29}_{\pm 0.42}$ | ${56.68}_{\pm 1.09}$
DeSurv | ${16.71}_{\pm 0.75}$ | ${49.61}_{\pm 2.15}$ | ${17.98}_{\pm 0.46}$ | ${53.23}_{\pm 1.15}$ | ${10.06}_{\pm 0.62}$ | ${33.18}_{\pm 1.93}$ | ${19.50}_{\pm 0.40}$ | ${57.28}_{\pm 0.89}$
NFM-PF | $\underline{16.33}_{\pm 0.75}$ | $\underline{49.07}_{\pm 1.96}$ | $\underline{17.60}_{\pm 0.55}$ | ${52.12}_{\pm 1.34}$ | $\mathbf{9.96}_{\pm 0.39}$ | $\mathbf{32.84}_{\pm 1.15}$ | ${19.14}_{\pm 0.39}$ | ${56.35}_{\pm 1.00}$
NFM-FN | $\mathbf{16.11}_{\pm 0.81}$ | $\mathbf{48.21}_{\pm 2.04}$ | ${17.66}_{\pm 0.52}$ | ${52.41}_{\pm 1.22}$ | ${10.05}_{\pm 0.39}$ | ${33.11}_{\pm 1.10}$ | $\mathbf{18.97}_{\pm 0.60}$ | $\mathbf{55.87}_{\pm 1.50}$

Table 2: Survival prediction results measured in IBS and INBLL metric (%) on two larger datasets. In each column, the boldfaced score denotes the best result and the underlined score represents the second-best result. Two models are not reported, namely SODEN and DeepEH, as we found empirically that their computational/memory cost is significantly worse than the rest, and we fail to obtain reasonable performances over the two datasets for these two models. Model | MIMIC-III | KKBOX
---|---|---
| IBS | INBLL | IBS | INBLL
CoxPH | ${20.40}_{\pm 0.00}$ | ${60.02}_{\pm 0.00}$ | ${12.60}_{\pm 0.00}$ | ${39.40}_{\pm 0.00}$
GBM | ${17.70}_{\pm 0.00}$ | ${52.30}_{\pm 0.00}$ | ${11.81}_{\pm 0.00}$ | ${38.15}_{\pm 0.00}$
RSF | ${17.79}_{\pm 0.19}$ | ${53.34}_{\pm 0.41}$ | ${14.46}_{\pm 0.00}$ | ${44.39}_{\pm 0.00}$
DeepSurv | ${18.58}_{\pm 0.92}$ | ${55.98}_{\pm 2.43}$ | ${11.31}_{\pm 0.05}$ | ${35.28}_{\pm 0.15}$
CoxTime | ${17.68}_{\pm 1.36}$ | ${52.08}_{\pm 3.06}$ | $\underline{10.70}_{\pm 0.06}$ | $\underline{33.10}_{\pm 0.21}$
DeepHit | ${19.80}_{\pm 1.31}$ | ${59.03}_{\pm 4.20}$ | ${16.00}_{\pm 0.34}$ | ${48.64}_{\pm 1.04}$
SuMo-net | ${18.62}_{\pm 1.23}$ | ${54.51}_{\pm 2.97}$ | ${11.58}_{\pm 0.11}$ | ${36.61}_{\pm 0.28}$
DCM | ${18.02}_{\pm 0.49}$ | ${52.83}_{\pm 0.94}$ | ${10.71}_{\pm 0.11}$ | ${33.24}_{\pm 0.06}$
DeSurv | ${18.19}_{\pm 0.65}$ | ${54.69}_{\pm 2.83}$ | ${10.77}_{\pm 0.21}$ | ${33.22}_{\pm 0.10}$
NFM-PF | $\mathbf{16.28}_{\pm 0.36}$ | $\mathbf{49.18}_{\pm 0.92}$ | ${11.02}_{\pm 0.11}$ | ${35.10}_{\pm 0.22}$
NFM-FN | $\underline{17.47}_{\pm 0.45}$ | $\underline{51.48}_{\pm 1.23}$ | $\mathbf{10.63}_{\pm 0.08}$ | $\mathbf{32.81}_{\pm 0.14}$

Evaluation strategy We use two standard metrics in survival predictions for
evaluating model performance: integrated Brier score (IBS) and integrated
negative binomial log-likelihood (INBLL). Both metrics are derived from the
following:

$\displaystyle\begin{aligned}
\mathcal{S}(\ell,t_{1},t_{2})=\int_{t_{2}}^{t_{1}}\dfrac{1}{n}\sum_{i=1}^{n}\left[\dfrac{\ell(0,\widehat{S}(t|Z_{i}))I(T_{i}\leq
t,\delta_{i}=1)}{\widehat{S}_{C}(T_{i})}+\dfrac{\ell(1,\widehat{S}(t|Z_{i}))I(T_{i}>t)}{\widehat{S}_{C}(t)}\right]dt.\end{aligned}$
(11)

Where $\widehat{S}_{C}(t)$ is an estimate of the survival function $S_{C}(t)$
of the censoring variable, obtained by the Kaplan-Meier estimate Kaplan and
Meier (1958) of the censored observations on the test data.
$\ell:\\{0,1\\}\times[0,1]\mapsto\mathbb{R}^{+}$ is some proper loss function
for binary classification Gneiting and Raftery (2007). The IBS metric
corresponds to $\ell$ being the square loss, and the INBLL metric corresponds
to $\ell$ being the negative binomial (Bernoulli) log-likelihood Graf et al.
(1999). Both IBS and INBLL are proper scoring rules if the censoring times and
survival times are independent. 444Otherwise, one may pose a covariate-
dependent model on the censoring time and use $\widehat{S}_{C}(t|Z)$ instead
of $\widehat{S}_{C}(t)$. We adopt the Kaplan-Meier approach since it’s still
the prevailing practice in evaluations of survival predictions. We
additionally report the result of another widely used metric, the concordance
index (C-index), in appendix D. Since all the survival datasets do not have
standard train/test splits, we follow previous practice Zhong et al. (2021a)
that uses $5$-fold cross-validation (CV): $1$ fold is for testing, and $20\%$
of the rest is held out for validation. In our experiments, we observed that a
single random split into $5$ folds does not produce stable results for most
survival datasets. Therefore we perform $10$ different CV runs for each
survival dataset and report average metrics as well as their standard
deviations. For the KKBOX dataset, we use the standard train/valid/test splits
that are available via the pycox package Kvamme et al. (2019) and report
results based on $10$ trial runs.

Experimental setup We follow standard preprocessing strategies Katzman et al.
(2018); Kvamme et al. (2019); Zhong et al. (2021a) that standardize continuous
features into zero mean and unit variance, and do one-hot encodings for all
categorical features. We adopt MLP with ReLU activation for all function
approximators, including $\widehat{h}$, $\widehat{m}$ in PF scheme, and
$\widehat{\nu}$ in FN scheme, across all datasets, with the number of layers
(depth) and the number of hidden units (width) within each layer being
tunable. We tune the frailty transform over several standard choices detailed
in appendix C.3. We find that the gamma frailty configuration performs
reasonably well across all tasks and is recommended to be the default choice.
A more detailed description of the tuning procedure, as well as training
configurations for baseline models, are reported in appendix C.3.

Results we report experimental results of small-scale datasets in table 1, and
results of two larger datasets in table 2. The proposed NFM framework achieves
the best performance on $5$ of the $6$ datasets. The improvement over
baselines is particularly evident in METABRIC, SUPPORT, and MIMIC-III
datasets.

Benefits of frailty to better understand the additional benefits of
introducing the frailty formulation, we compute the (relative) performance
gain of NFM-PF and NFM-FN, against their non-frailty counterparts, namely
DeepSurv Katzman et al. (2018) and SuMo-net Rindt et al. (2022). The
evaluation is conducted for all three metrics mentioned in this paper. The
results are shown in table 5. The results suggest a solid improvement in
incorporating frailty, as the relative increase in performance could be over
$10\%$ for both NFM models. A more detailed discussion is presented in section
D.4.

##  6 Conclusion

In this paper, we make principled explorations on applying the idea of frailty
models in modern survival analysis to neural survival regressions. A flexible
and scalable framework called NFM is proposed that includes many useful
survival models as special cases. Under the framework, we study two derived
model architectures both theoretically and empirically. Theoretically, we
obtain the rates of convergences of the nonparametric function estimators
based on neural function approximation. Empirically, we demonstrate the
superior predictive performance of the proposed models by evaluating several
benchmark datasets.

## References

  * Kalbfleisch and Prentice (2002) John D Kalbfleisch and Ross L Prentice.  _The Statistical Analysis of Failure Time Data_ , volume 360.  John Wiley & Sons, 2002.
  * Duffie et al. (2009) Darrell Duffie, Andreas Eckner, Guillaume Horel, and Leandro Saita.  Frailty correlated default.  _The Journal of Finance_ , 64(5):2089–2123, 2009\.
  * Jing and Smola (2017) How Jing and Alexander J Smola.  Neural survival recommender.  In _Proceedings of the Tenth ACM International Conference on Web Search and Data Mining_ , pages 515–524, 2017.
  * Wu et al. (2015) Wush Chi-Hsuan Wu, Mi-Yen Yeh, and Ming-Syan Chen.  Predicting winning price in real time bidding with censored data.  In _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 1305–1314, 2015.
  * Buckley and James (1979) Jonathan Buckley and Ian James.  Linear regression with censored data.  _Biometrika_ , 66(3):429–436, 1979.
  * Cox (1972) David R Cox.  Regression models and life-tables.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 34(2):187–202, 1972.
  * Cox (1975) David R Cox.  Partial likelihood.  _Biometrika_ , 62(2):269–276, 1975.
  * Bickel et al. (1993) Peter J Bickel, Chris AJ Klaassen, Ya’acov Ritov, J Klaassen, Jon A Wellner, and YA’Acov Ritov.  _Efficient and adaptive estimation for semiparametric models_ , volume 4.  Springer, 1993.
  * Huang (1999) Jian Huang.  Efficient estimation of the partly linear additive cox model.  _The annals of Statistics_ , 27(5):1536–1563, 1999.
  * Cai et al. (2007) Jianwen Cai, Jianqing Fan, Jiancheng Jiang, and Haibo Zhou.  Partially linear hazard regression for multivariate survival data.  _Journal of the American Statistical Association_ , 102(478):538–551, 2007.
  * Cai et al. (2008) Jianwen Cai, Jianqing Fan, Jiancheng Jiang, and Haibo Zhou.  Partially linear hazard regression with varying coefficients for multivariate survival data.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 70(1):141–158, 2008.
  * Ishwaran et al. (2008) Hemant Ishwaran, Udaya B Kogalur, Eugene H Blackstone, and Michael S Lauer.  Random survival forests.  _The annals of applied statistics_ , 2(3):841–860, 2008.
  * Faraggi and Simon (1995) David Faraggi and Richard Simon.  A neural network model for survival data.  _Statistics in medicine_ , 14(1):73–82, 1995\.
  * Katzman et al. (2018) Jared L Katzman, Uri Shaham, Alexander Cloninger, Jonathan Bates, Tingting Jiang, and Yuval Kluger.  Deepsurv: personalized treatment recommender system using a cox proportional hazards deep neural network.  _BMC medical research methodology_ , 18(1):1–12, 2018.
  * Gray (2000) Robert J Gray.  Estimation of regression parameters and the hazard function in transformed linear survival models.  _Biometrics_ , 56(2):571–576, 2000.
  * Stablein and Koutrouvelis (1985) Donald M Stablein and IA Koutrouvelis.  A two-sample test sensitive to crossing hazards in uncensored and singly censored data.  _Biometrics_ , pages 643–652, 1985.
  * Wienke (2010) Andreas Wienke.  _Frailty models in survival analysis_.  Chapman and Hall/CRC, 2010.
  * Duchateau and Janssen (2007) Luc Duchateau and Paul Janssen.  _The frailty model_.  Springer Science & Business Media, 2007.
  * Bennett (1983) Steve Bennett.  Analysis of survival data by the proportional odds model.  _Statistics in medicine_ , 2(2):273–277, 1983\.
  * Murphy (1994) Susan A Murphy.  Consistency in a proportional hazards model incorporating a random effect.  _The Annals of Statistics_ , 22(2):712–731, 1994\.
  * Murphy (1995) Susan A Murphy.  Asymptotic theory for the frailty model.  _The annals of statistics_ , pages 182–198, 1995.
  * Parner (1998) Erik Parner.  Asymptotic theory for the correlated gamma-frailty model.  _The Annals of Statistics_ , 26(1):183–214, 1998\.
  * Kosorok et al. (2004) Michael R Kosorok, Bee Leng Lee, and Jason P Fine.  Robust inference for univariate proportional hazards frailty regression models.  _The Annals of Statistics_ , 32(4):1448–1491, 2004.
  * Kvamme et al. (2019) Hvard Kvamme, Ørnulf Borgan, and Ida Scheel.  Time-to-event prediction with neural networks and cox regression.  _arXiv preprint arXiv:1907.00825_ , 2019.
  * Tang et al. (2022) Weijing Tang, Jiaqi Ma, Qiaozhu Mei, and Ji Zhu.  Soden: A scalable continuous-time survival model through ordinary differential equation networks.  _Journal of Machine Learning Research_ , 23(34):1–29, 2022.  URL http://jmlr.org/papers/v23/20-900.html.
  * Rindt et al. (2022) D Rindt, R Hu, D Steinsaltz, and D Sejdinovic.  Survival regression with proper scoring rules and monotonic neural networks.  Proceedings of Machine Learning Research. Journal of Machine Learning Research, 2022.
  * Goodfellow et al. (2016) Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  _Deep learning_.  2016\.
  * Cox and O’Sullivan (1990) Dennis D Cox and Finbarr O’Sullivan.  Asymptotic analysis of penalized likelihood and related estimators.  _The Annals of Statistics_ , pages 1676–1695, 1990.
  * Kooperberg et al. (1995) Charles Kooperberg, Charles J Stone, and Young K Truong.  Hazard regression.  _Journal of the American Statistical Association_ , 90(429):78–94, 1995.
  * Omi et al. (2019) Takahiro Omi, naonori ueda, and Kazuyuki Aihara.  Fully neural network based model for general temporal point processes.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper/2019/file/39e4973ba3321b80f37d9b55f63ed8b8-Paper.pdf.
  * Murphy and Van der Vaart (2000) Susan A Murphy and Aad W Van der Vaart.  On profile likelihood.  _Journal of the American Statistical Association_ , 95(450):449–465, 2000.
  * Fan and Gijbels (1996) Jianqing Fan and Irene Gijbels.  _Local Polynomial Modelling and Its Applications: Monographs on Statistics and Applied Probability 66_ , volume 66.  CRC Press, 1996.
  * Hardt et al. (2016) Moritz Hardt, Ben Recht, and Yoram Singer.  Train faster, generalize better: Stability of stochastic gradient descent.  In _International conference on machine learning_ , pages 1225–1234. PMLR, 2016.
  * Ying (1993) Zhiliang Ying.  A large sample study of rank estimation for censored regression data.  _The Annals of Statistics_ , pages 76–99, 1993.
  * Etezadi-Amoli and Ciampi (1987) Jamshid Etezadi-Amoli and Antonio Ciampi.  Extended hazard regression for censored survival data with covariates: a spline approximation for the baseline hazard function.  _Biometrics_ , pages 181–192, 1987.
  * Zeng and Lin (2006) Donglin Zeng and DY Lin.  Efficient estimation of semiparametric transformation models for counting processes.  _Biometrika_ , 93(3):627–640, 2006.
  * Strawderman and Tsiatis (1996) Robert L Strawderman and Anastasios A Tsiatis.  On the asymptotic properties of a flexible hazard estimator.  _The Annals of Statistics_ , 24(1):41–63, 1996\.
  * Lee et al. (2018) Changhee Lee, William Zame, Jinsung Yoon, and Mihaela Van Der Schaar.  Deephit: A deep learning approach to survival analysis with competing risks.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 32, 2018.
  * Zhong et al. (2021a) Qixian Zhong, Jonas W Mueller, and Jane-Ling Wang.  Deep extended hazard models for survival analysis.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, _Advances in Neural Information Processing Systems_ , volume 34, pages 15111–15124. Curran Associates, Inc., 2021a.  URL https://proceedings.neurips.cc/paper/2021/file/7f6caf1f0ba788cd7953d817724c2b6e-Paper.pdf.
  * Zhong et al. (2021b) Qixian Zhong, Jonas W Mueller, and Jane-Ling Wang.  Deep learning for the partially linear cox model.  _The Annals of Statistics_ , 2021b.
  * Brennan (2002) Paul Brennan.  Gene–environment interaction and aetiology of cancer: what does it mean and how can we measure it?  _Carcinogenesis_ , 23(3):381–387, 2002.
  * Boyd (2001) John P Boyd.  _Chebyshev and Fourier spectral methods_.  Courier Corporation, 2001.
  * Wehenkel and Louppe (2019) Antoine Wehenkel and Gilles Louppe.  Unconstrained monotonic neural networks.  _Advances in neural information processing systems_ , 32, 2019.
  * Shen and Wong (1994) Xiaotong Shen and Wing Hung Wong.  Convergence rate of sieve estimates.  _The Annals of Statistics_ , pages 580–615, 1994.
  * Shen (1997) Xiaotong Shen.  On methods of sieves and penalization.  _The Annals of Statistics_ , 25(6):2555–2591, 1997.
  * Chen (2007) Xiaohong Chen.  Large sample sieve estimation of semi-nonparametric models.  _Handbook of econometrics_ , 6:5549–5632, 2007.
  * Yarotsky (2017) Dmitry Yarotsky.  Error bounds for approximations with deep relu networks.  _Neural Networks_ , 94:103–114, 2017.
  * Schmidt-Hieber (2020) Johannes Schmidt-Hieber.  Nonparametric regression using deep neural networks with relu activation function.  _The Annals of Statistics_ , 48(4):1875–1897, 2020.
  * Farrell et al. (2021) Max H Farrell, Tengyuan Liang, and Sanjog Misra.  Deep neural networks for estimation and inference.  _Econometrica_ , 89(1):181–213, 2021.
  * Van der Vaart (2000) Aad W Van der Vaart.  _Asymptotic statistics_ , volume 3.  Cambridge university press, 2000.
  * Wong and Shen (1995) Wing Hung Wong and Xiaotong Shen.  Probability inequalities for likelihood ratios and convergence rates of sieve mles.  _The Annals of Statistics_ , pages 339–362, 1995.
  * Curtis et al. (2012) Christina Curtis, Sohrab P Shah, Suet-Feung Chin, Gulisa Turashvili, Oscar M Rueda, Mark J Dunning, Doug Speed, Andy G Lynch, Shamith Samarajiwa, Yinyin Yuan, et al.  The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups.  _Nature_ , 486(7403):346–352, 2012.
  * Knaus et al. (1995) William A Knaus, Frank E Harrell, Joanne Lynn, Lee Goldman, Russell S Phillips, Alfred F Connors, Neal V Dawson, William J Fulkerson, Robert M Califf, Norman Desbiens, et al.  The support prognostic model: Objective estimates of survival for seriously ill hospitalized adults.  _Annals of internal medicine_ , 122(3):191–203, 1995.
  * Dispenzieri et al. (2012) Angela Dispenzieri, Jerry A Katzmann, Robert A Kyle, Dirk R Larson, Terry M Therneau, Colin L Colby, Raynell J Clark, Graham P Mead, Shaji Kumar, L Joseph Melton III, et al.  Use of nonclonal serum immunoglobulin free light chains to predict overall survival in the general population.  In _Mayo Clinic Proceedings_ , volume 87, pages 517–523. Elsevier, 2012.
  * Johnson et al. (2016) Alistair EW Johnson, Tom J Pollard, Lu Shen, Li-wei H Lehman, Mengling Feng, Mohammad 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.
  * Friedman (2001) Jerome H Friedman.  Greedy function approximation: a gradient boosting machine.  _Annals of statistics_ , pages 1189–1232, 2001.
  * Chen and Guestrin (2016) Tianqi Chen and Carlos Guestrin.  Xgboost: A scalable tree boosting system.  In _Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining_ , pages 785–794, 2016.
  * Groha et al. (2020) Stefan Groha, Sebastian M Schmon, and Alexander Gusev.  A general framework for survival analysis and multi-state modelling.  _arXiv preprint arXiv:2006.04893_ , 2020.
  * Nagpal et al. (2021) Chirag Nagpal, Steve Yadlowsky, Negar Rostamzadeh, and Katherine Heller.  Deep cox mixtures for survival regression.  In _Machine Learning for Healthcare Conference_ , pages 674–708. PMLR, 2021.
  * Danks and Yau (2022) Dominic Danks and Christopher Yau.  Derivative-based neural modelling of cumulative distribution functions for survival analysis.  In _International Conference on Artificial Intelligence and Statistics_ , pages 7240–7256. PMLR, 2022.
  * Kaplan and Meier (1958) Edward L Kaplan and Paul Meier.  Nonparametric estimation from incomplete observations.  _Journal of the American statistical association_ , 53(282):457–481, 1958.
  * Gneiting and Raftery (2007) Tilmann Gneiting and Adrian E Raftery.  Strictly proper scoring rules, prediction, and estimation.  _Journal of the American statistical Association_ , 102(477):359–378, 2007.
  * Graf et al. (1999) Erika Graf, Claudia Schmoor, Willi Sauerbrei, and Martin Schumacher.  Assessment and comparison of prognostic classification schemes for survival data.  _Statistics in medicine_ , 18(17-18):2529–2545, 1999.
  * Hougaard (1984) Philip Hougaard.  Life table methods for heterogeneous populations: distributions describing the heterogeneity.  _Biometrika_ , 71(1):75–83, 1984.
  * van der Vaart et al. (1996) AW van der Vaart, A.W. van der Vaart, A. van der Vaart, and J. Wellner.  _Weak Convergence and Empirical Processes: With Applications to Statistics_.  Springer Series in Statistics. Springer, 1996.  ISBN 9780387946405.  URL https://books.google.com/books?id=OCenCW9qmp4C.
  * Bartlett et al. (2019) Peter L Bartlett, Nick Harvey, Christopher Liaw, and Abbas Mehrabian.  Nearly-tight vc-dimension and pseudodimension bounds for piecewise linear neural networks.  _The Journal of Machine Learning Research_ , 20(1):2285–2301, 2019.
  * Chen and Shen (1998) Xiaohong Chen and Xiaotong Shen.  Sieve extremum estimates for weakly dependent data.  _Econometrica_ , pages 289–314, 1998.
  * Dabrowska and Doksum (1988) Dorota M Dabrowska and Kjell A Doksum.  Partial likelihood in transformation models with censored data.  _Scandinavian journal of statistics_ , pages 1–23, 1988.
  * Cuzick (1988) Jack Cuzick.  Rank regression.  _The Annals of Statistics_ , pages 1369–1389, 1988.
  * Purushotham et al. (2018) Sanjay Purushotham, Chuizheng Meng, Zhengping Che, and Yan Liu.  Benchmarking deep learning models on large healthcare datasets.  _Journal of biomedical informatics_ , 83:112–134, 2018\.
  * Pölsterl (2020) Sebastian Pölsterl.  scikit-survival: A library for time-to-event analysis built on top of scikit-learn.  _Journal of Machine Learning Research_ , 21(212):1–6, 2020.  URL http://jmlr.org/papers/v21/20-729.html.
  * Antolini et al. (2005) Laura Antolini, Patrizia Boracchi, and Elia Biganzoli.  A time-dependent discrimination index for survival data.  _Statistics in medicine_ , 24(24):3927–3944, 2005\.
  * Steck et al. (2007) Harald Steck, Balaji Krishnapuram, Cary Dehing-oberije, Philippe Lambin, and Vikas C Raykar.  On ranking in survival analysis: Bounds on the concordance index.  In J. Platt, D. Koller, Y. Singer, and S. Roweis, editors, _Advances in Neural Information Processing Systems_ , volume 20. Curran Associates, Inc., 2007.  URL https://proceedings.neurips.cc/paper/2007/file/33e8075e9970de0cfea955afd4644bb2-Paper.pdf.

##  Appendix A Examples of frailty specifications

We list several commonly used frailty models, and specify their corresponding
characteristics via their frailty transform $G_{\theta}$:

Gamma frailty:

Arguably the gamma frailty is the most widely used frailty model Murphy [1994,
1995], Parner [1998], Wienke [2010], Duchateau and Janssen [2007], with

$\displaystyle G_{\theta}(x)=\frac{1}{\theta}\log(1+\theta x),\theta\geq 0.$
(12)

When $\theta=0$, $G_{0}(x)=\lim_{\theta\rightarrow 0}G_{\theta}(x)$ is defined
as the (pointwise) limit. A notable fact of the gamma frailty specification is
that when the proportional frailty (PF) assumption (2) is met, if $\theta=0$,
the model degenerates to CoxPH . Otherwise if $\theta=1$, the model
corresponds to the proportional odds (PO) model Bennett [1983].

Box-Cox transformation frailty:

Under this specification, we have

$\displaystyle G_{\theta}(x)=\dfrac{(1+x)^{\theta}-1}{\theta},\theta\geq 0.$
(13)

The case of $\theta=0$ is defined analogously to that of gamma frailty, which
corresponds to the PO model under the PF assumption. When $\theta=1$, the
model reduces to CoxPH under the PF assumption.

$\text{IGG}(\alpha)$ frailty:

This is an extension of gamma frailty Kosorok et al. [2004] and includes other
types of frailty specifications like the inverse gaussian frailty Hougaard
[1984], with

$\displaystyle
G_{\theta}(x)=\frac{1-\alpha}{\alpha\theta}\left[\left(1+\frac{\theta
x}{1-\alpha}\right)^{\alpha}-1\right],\theta\geq 0,\alpha\in[0,1).$ (14)

In the one-dimensional parameter paradigm, the parameter $\alpha$ is assumed
known instead of being learnable. When $\alpha=1/2$, we obtain the gamma
frailty model. When $\alpha\rightarrow 0$, the limit corresponds to the
inverse Gaussian frailty.

Satistiability of regularity condition 4.5 In [Kosorok et al., 2004,
Proposition 1], the authors verified the regularity condition of gamma and
$\text{IGG}(\alpha)$ frailties. Using a similar argument, it is
straightforward to verify the regularity of Box-Cox transformation frailty.

##  Appendix B Proofs of theorems

###  B.1 Preliminary

#### Additional definitions

The theory of empirical processes [van der Vaart et al., 1996] will be
involved heavily in the proof. Therefore we briefly introduce some common
notations: For a function class $\mathcal{F}$, define
$N\left(\epsilon,\mathcal{F},\|\cdot\|\right)$ to be the covering number of
$\mathcal{F}$ with respect to norm $\|\cdot\|$ under radius $\epsilon$, and
define $N_{[]}\left(\epsilon,\mathcal{F},\|\cdot\|\right)$ to be the
bracketing number of $\mathcal{F}$ with respect to norm $\|\cdot\|$ under
radius $\epsilon$. We use $\text{VC}\left(\mathcal{F}\right)$ to denote the
VC-dimension of $\mathcal{F}$. Moreover, we use the notation $a\lesssim b$ to
denote $a\leq Cb$ for some positive constant $C$.

Before proving theorem 4.6 and 4.7, we introduce some additional notations
that will be useful throughout the proof process.

In the PF scheme, define

$\displaystyle l(T,\delta,Z;h,m,\theta)=$ $\displaystyle\delta\log
g_{\theta}\left(e^{m(Z)}\int_{0}^{T}e^{h(s)}ds\right)+\delta h(T)+\delta m(Z)$
$\displaystyle-G_{\theta}\left(e^{m(Z)}\int_{0}^{T}e^{h(s)}ds\right),$

where we denote $g_{\theta}=G^{\prime}(\theta)$. Under the definition of the
sieve space stated in condition 4.3, we restate the parameter estimates as

$\displaystyle\left(\widehat{h}_{n},\widehat{m}_{n},\widehat{\theta}_{n}\right)=\mathop{\mathrm{argmax}}\limits_{\widehat{h}\in\mathcal{H}_{n},\widehat{m}\in\mathcal{M}_{n},\theta\in\Theta}\frac{1}{n}\sum_{i\in[n]}l(T_{i},\delta_{i},Z_{i};\widehat{h},\widehat{m},\theta).$

Similarly, in the FN scheme, we define

$\displaystyle l(T,\delta,Z;\nu,\theta)=\delta\log
g_{\theta}\left(\int_{0}^{T}e^{\nu(s,Z)}ds\right)+\delta\nu(T,Z)-G_{\theta}\left(\int_{0}^{T}e^{\nu(s,Z)}ds\right)$

Under the definition of the sieve space stated in condition 4.4, we restate
the parameter estimates as

$\displaystyle\left(\widehat{\nu}_{n}(t,z),\widehat{\theta}_{n}\right)=\mathop{\mathrm{argmax}}\limits_{\widehat{\nu}\in\mathcal{V}_{n},\theta\in\Theta}\frac{1}{n}\sum_{i\in[n]}l(T_{i},\delta_{i},Z_{i};\widehat{\nu},\theta).$

We denote the conditional density function and survival function of the event
time $\tilde{T}$ given $Z$ by $f_{\tilde{T}\mid Z}(t)$ and $S_{\tilde{T}\mid
Z}(t)$, respectively. Similarly, we denote the conditional density function
and survival function of the censoring time $C$ given $Z$ by $f_{C\mid Z}(t)$
and $S_{C\mid Z}(t)$. Under the assumption that
$\tilde{T}\perp\\!\\!\\!\\!\perp C\mid Z$, the joint conditional density of
the observed time $T$ and the censoring indicator $\delta$ given $Z$ can be
expressed as the following:

$\displaystyle p(T,\delta\mid Z)$ $\displaystyle=$ $\displaystyle
f_{\tilde{T}\mid Z}(T)^{\delta}S_{\tilde{T}\mid Z}(T)^{1-\delta}f_{C\mid
Z}(T)^{1-\delta}S_{C\mid Z}(T)^{\delta}$ $\displaystyle=$
$\displaystyle\lambda_{\tilde{T}\mid Z}(T)^{\delta}S_{\tilde{T}\mid
Z}(T)f_{C\mid Z}(T)^{1-\delta}S_{C\mid Z}(T)^{\delta},$

where $\lambda_{\tilde{T}\mid Z}(T)$ is the conditional hazard function of the
survival time $\tilde{T}$ given $Z$.

Under the model assumption of PF scheme, $p(T,\delta\mid Z)$ can be expressed
by

$\displaystyle p(T,\delta\mid
Z;h,m,\theta)=\exp\left(l(T,\delta,Z;h,m,\theta)\right)f_{C\mid
Z}(T)^{1-\delta}S_{C\mid Z}(T)^{\delta}.$

For $\phi_{0}=(h_{0},m_{0},\theta_{0})$ and an estimator
$\widehat{\phi}=(\widehat{h},\widehat{m},\widehat{\theta})$, the defined
distance $d_{\textsf{PF}}\left(\widehat{\phi},\phi_{0}\right)$ can be
explicitly expresses by

$\displaystyle
d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)=\sqrt{\mathbb{E}_{Z}\left[\int\left|\sqrt{p(T,\delta\mid
Z;\widehat{h},\widehat{m},\widehat{\theta})}-\sqrt{p(T,\delta\mid
Z;h_{0},m_{0},\theta_{0})}\right|^{2}\mu(dT\times d\delta)\right]}.$

Here the dominating measure $\mu$ is defined such that for any (measurable)
function $r(T,\delta)$

$\displaystyle\int r(T,\delta)\mu(dT\times
d\delta)=\int_{0}^{\tau}r(T,\delta=1)dT+\int_{0}^{\tau}r(T,\delta=0)dT$

Under the model assumption of FN scheme, $p(T,\delta\mid Z)$ can be expressed
by

$\displaystyle p(T,\delta\mid
Z;\nu,\theta)=\exp\left(l(T,\delta,Z;\nu,\theta)\right)f_{C\mid
Z}(T)^{1-\delta}S_{C\mid Z}(T)^{\delta}.$

For $\psi_{0}=(\nu_{0},\theta_{0})$ and an estimator
$\widehat{\psi}=(\widehat{\nu},\widehat{\theta})$, the defined distance
$d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)$ can be explicitly
expresses by

$\displaystyle
d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)=\sqrt{\mathbb{E}_{Z}\left[\int\left|\sqrt{p(T,\delta\mid
Z;\widehat{\nu},\widehat{\theta})}-\sqrt{p(T,\delta\mid
Z;\nu_{0},\theta_{0})}\right|^{2}\mu(dT\times d\delta)\right]}.$

###  B.2 Technical lemmas

The following lemmas are needed for the proof of Theorem 4.6 and 4.7.
Hereafter for notational convenience, we will use $\widehat{h},\widehat{m}$
for arbitrary elements in the corresponding sieve space listed in condition
4.3, $\widehat{\nu}$ for an arbitrary element in the sieve space listed in
condition 4.4, and $\widehat{\theta}$ for an arbitrary element in $\Theta$.

######  Lemma B.1.

Under condition 4.1, 4.3, 4.5, for
$(T,\delta,Z)\in[0,\tau]\times\\{0,1\\}\times[-1,1]^{d}$, the following terms
are bounded:

  1. 1.

$l(T,\delta,Z;h_{0},m_{0},\theta_{0})$ with true parameter
$(h_{0},m_{0},\theta_{0})$

  2. 2.

$l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})$ with parameter
estimates $(\widehat{h},\widehat{m},\widehat{\theta})$ in any sieve space
listed in condition 4.3.

######  Lemma B.2.

Under condition 4.2, 4.4, 4.5, for
$(T,\delta,Z)\in[0,\tau]\times\\{0,1\\}\times[-1,1]^{d}$, the following terms
are bounded:

  1. 1.

$l(T,\delta,Z;\nu_{0},\theta_{0})$ with true parameter $(\nu_{0},\theta_{0})$

  2. 2.

$l(T,\delta,Z;\widehat{\nu},\widehat{\theta})$ with parameter estimates
$(\widehat{\nu},\widehat{\theta})$ in any sieve space listed in condition 4.4.

######  Lemma B.3.

Under condition 4.1, 4.3, 4.5, let
$(\widehat{h},\widehat{m},\widehat{\theta})$,
$(\widehat{h}_{1},\widehat{m}_{1},\widehat{\theta}_{1})$, and
$(\widehat{h}_{2},\widehat{m}_{2},\widehat{\theta}_{2})$ be arbitrary three
parameter triples inside the sieve space defined in condition 4.3, the
following two inequalities hold.

$\displaystyle\|l(T,\delta,Z;h_{0},m_{0},\theta_{0})-l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})\|_{\infty}\lesssim|\theta_{0}-\widehat{\theta}|+\|h_{0}-\widehat{h}\|_{\infty}+\|m_{0}-\widehat{m}\|_{\infty}$
$\displaystyle\|l(T,\delta,Z;\widehat{h}_{1},\widehat{m}_{1},\widehat{\theta}_{1})-l(T,\delta,Z;\widehat{h}_{2},\widehat{m}_{2},\widehat{\theta}_{2})\|_{\infty}\lesssim|\widehat{\theta}_{1}-\widehat{\theta}_{2}|+\|\widehat{h}_{1}-\widehat{h}_{2}\|_{\infty}+\|\widehat{m}_{1}-\widehat{m}_{2}\|_{\infty}.$

######  Lemma B.4.

Under condition 4.2, 4.4, 4.5, let $(\widehat{\nu},\widehat{\theta})$,
$(\widehat{\nu}_{1},\widehat{\theta}_{1})$, and
$(\widehat{\nu}_{2},\widehat{\theta}_{2})$ be arbitrary three parameter tuples
inside the sieve space defined in condition 4.4,, the following inequalities
hold.

$\displaystyle\|l(T,\delta,Z;\nu_{0},\theta_{0})-l(T,\delta,Z;\widehat{\nu},\widehat{\theta})\|_{\infty}\lesssim|\theta_{0}-\widehat{\theta}|+\|\nu_{0}-\widehat{\nu}\|_{\infty}$
$\displaystyle\|l(T,\delta,Z;\widehat{\nu}_{1},\widehat{\theta}_{1})-l(T,\delta,Z;\widehat{\nu}_{2},\widehat{\theta}_{2})\|_{\infty}\lesssim|\widehat{\theta}_{1}-\widehat{\theta}_{2}|+\|\widehat{\nu}_{1}-\widehat{\nu}_{2}\|_{\infty}.$

######  Lemma B.5 (Approximating error of PF scheme).

In the PF scheme, for any $n$, there exists an element in the corresponding
sieve space $\pi_{n}\phi_{0}=(\pi_{n}h_{0},\pi_{n}m_{0},\pi_{n}\theta_{0})$,
satisfying
$d_{\textsf{PF}}\left(\pi_{n}\phi_{0},\phi_{0}\right)=O\left(n^{-\frac{\beta}{\beta+d}}\right)$.

######  Lemma B.6 (Approximating error of FN scheme).

In the FN scheme, for any $n$, there exists an element in the corresponding
sieve space $\pi_{n}\psi=(\pi_{n}\nu_{0},\pi_{n}\theta_{0})$ satisfying
$d_{\textsf{FN}}\left(\pi_{n}\psi_{0},\psi_{0}\right)=O\left(n^{-\frac{\beta}{\beta+d+1}}\right)$.

######  Lemma B.7.

Suppose $\mathcal{F}$ is a class of functions satisfying that
$N(\varepsilon,\mathcal{F},\|\cdot\|)<\infty$ for $\forall\varepsilon>0$. We
define $\widetilde{N}(\varepsilon,\mathcal{F},\|\cdot\|)$ to be the minimal
number of $\varepsilon$-balls $B(f,\varepsilon)=\\{g:\|g-f\|<\varepsilon\\}$
needed to cover $\mathcal{F}$ with radius $\varepsilon$ and further constrain
that $f\in\mathcal{F}$. Then we have

$\displaystyle
N(\varepsilon,\mathcal{F},\|\cdot\|)\leq\widetilde{N}(\varepsilon,\mathcal{F},\|\cdot\|)\leq
N(\frac{\varepsilon}{2},\mathcal{F},\|\cdot\|).$

######  Lemma B.8.

Suppose $\mathcal{F}$ is a class of functions satisfying that
$N_{[]}(\varepsilon,\mathcal{F},\|\cdot\|_{\infty})<\infty$ for
$\forall\varepsilon>0$. We define
$\widetilde{N}_{[]}(\varepsilon,\mathcal{F},\|\cdot\|_{\infty})$ to be the
minimal number of brackets $[l,u]$ needed to cover $\mathcal{F}$ with
$\|l-u\|_{\infty}<\varepsilon$ and further constrain that $f\in\mathcal{F}$,
$l=f-\frac{\varepsilon}{2}$ and $u=f+\frac{\varepsilon}{2}$. Then we have

$\displaystyle
N_{[]}(\varepsilon,\mathcal{F},\|\cdot\|_{\infty})\leq\widetilde{N}_{[]}(\varepsilon,\mathcal{F},\|\cdot\|_{\infty})\leq
N_{[]}(\frac{\varepsilon}{2},\mathcal{F},\|\cdot\|_{\infty})$

Furthermore, we have
$\widetilde{N}_{[]}(\varepsilon,\mathcal{F},\|\cdot\|_{\infty})=\widetilde{N}(\frac{\varepsilon}{2},\mathcal{F},\|\cdot\|_{\infty})$.

######  Lemma B.9 (Model capacity of PF scheme).

Let
$\mathcal{F}_{n}=\\{l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta}):\widehat{h}\in\mathcal{H}_{n},\widehat{m}\in\mathcal{M}_{n},\widehat{\theta}\in\Theta\\}$.
Under condition 4.5, with $s_{h}=\frac{2\beta}{2\beta+1}$ and
$s_{m}=\frac{2\beta}{2\beta+d}$, there exist constants $c_{h}$ and $c_{m}>0$
such that

$\displaystyle
N_{[]}(\varepsilon,\mathcal{F}_{n},\|\cdot\|_{\infty})\lesssim\frac{1}{\varepsilon}N(c_{h}\varepsilon^{1/s_{h}},\mathcal{H}_{n},\|\cdot\|_{2})\times
N(c_{m}\varepsilon^{1/s_{m}},\mathcal{M}_{n},\|\cdot\|_{2}).$

######  Lemma B.10 (Model capacity of FN scheme).

Let
$\mathcal{G}_{n}=\\{l(T,\delta,Z;\widehat{\nu},\widehat{\theta}):\widehat{\nu}\in\mathcal{V}_{n},\widehat{\theta}\in\Theta\\}$.
Under condition 4.5, with $s_{\nu}=\frac{2\beta}{2\beta+d+1}$, there exists a
constant $c_{\nu}>0$ such that

$\displaystyle
N_{[]}(\varepsilon,\mathcal{G}_{n},\|\cdot\|_{\infty})\lesssim\frac{1}{\varepsilon}N(c_{\nu}\varepsilon^{1/s_{\nu}},\mathcal{V}_{n},\|\cdot\|_{2}).$

###  B.3 Proofs of theorem 4.6 and 4.7

###### Proof of theorem 4.6.

The proof is divided into four steps.

#### Step $1$

We denote $\phi_{0}=(h_{0},m_{0},\theta_{0})$ and
$\widehat{\phi}=(\widehat{h},\widehat{m},\widehat{\theta})$, where
$\widehat{h}\in\mathcal{H}_{n}$,$\widehat{m}\in\mathcal{M}_{n}$ and
$\widehat{\theta}\in\Theta$. For arbitrary small $\varepsilon>0$, we have that

$\displaystyle\inf_{d_{\textsf{PF}}\left(\widehat{\phi},\phi_{0}\right)\geq\varepsilon}\mathbb{E}\left[l(T,\delta,Z;h_{0},m_{0},\theta_{0})-l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})\right]$
$\displaystyle=\inf_{d_{\textsf{PF}}\left(\widehat{\phi},\phi_{0}\right)\geq\varepsilon}\mathbb{E}_{Z}\left[\mathbb{E}_{T,\delta\mid
Z}\left[\log p(T,\delta\mid Z;h_{0},m_{0},\theta_{0})-\log p(T,\delta\mid
Z;\widehat{h},\widehat{m},\widehat{\theta})\right]\right]$
$\displaystyle=\inf_{d_{\textsf{PF}}\left(\widehat{\phi},\phi_{0}\right)\geq\varepsilon}\mathbb{E}_{Z}\left[\text{D}_{\text{KL}}\left(\mathbb{P}_{\widehat{\phi},Z}\parallel\mathbb{P}_{\phi_{0},Z}\right)\right]$

Using the fact that
$\text{D}_{\text{KL}}\left(\mathbb{P}_{\widehat{\phi},Z}\parallel\mathbb{P}_{\phi_{0},Z}\right)\geq
2H^{2}(\mathbb{P}_{\widehat{\phi},Z}\parallel\mathbb{P}_{\phi_{0},Z})$. Thus,
we further obtain that

$\displaystyle\inf_{d_{\textsf{PF}}\left(\widehat{\phi},\phi_{0}\right)\geq\varepsilon}\mathbb{E}\left[l(T,\delta,Z;h_{0},m_{0},\theta_{0})-l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})\right]$
$\displaystyle\geq\inf_{d_{\textsf{PF}}\left(\widehat{\phi},\phi_{0}\right)\geq\varepsilon}2\mathbb{E}_{Z}\left[H^{2}(\mathbb{P}_{\widehat{\phi},Z}\parallel\mathbb{P}_{\phi_{0},Z})\right]$
$\displaystyle=2\inf_{d_{\textsf{PF}}\left(\widehat{\phi},\phi_{0}\right)\geq\varepsilon}d^{2}_{\textsf{PF}}\left(\widehat{\phi},\phi_{0}\right)$
$\displaystyle\geq 2\varepsilon^{2}.$

#### Step $2$

Consider the following derivation.

$\displaystyle\sup_{d_{\textsf{PF}}\left(\widehat{\phi},\phi_{0}\right)\leq\varepsilon}\text{Var}\left[l(T,\delta,Z;h_{0},m_{0},\theta_{0})-l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})\right]$
$\displaystyle\leq\sup_{d_{\textsf{PF}}\left(\widehat{\phi},\phi_{0}\right)\leq\varepsilon}\mathbb{E}\left[\left(l(T,\delta,Z;h_{0},m_{0},\theta_{0})-l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})\right)^{2}\right]$
$\displaystyle=\sup_{d_{\textsf{PF}}\left(\widehat{\phi},\phi_{0}\right)\leq\varepsilon}\mathbb{E}_{Z}\mathbb{E}_{T,\delta\mid
Z}\left[\log p(T,\delta,Z;h_{0},m_{0},\theta_{0})-\log
p(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})\right]^{2}$
$\displaystyle=4\sup_{d_{\textsf{PF}}\left(\widehat{\phi},\phi_{0}\right)\leq\varepsilon}\mathbb{E}_{Z}\left[\int\left(p(T,\delta,Z;h_{0},m_{0},\theta_{0})\left(\log\sqrt{\dfrac{p(T,\delta,Z;h_{0},m_{0},\theta_{0})}{p(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})}}\right)^{2}\right)\mu(dT\times
d\delta)\right]$

By Taylor’s expansion on $\log x$, there exists $\xi(T,\delta,Z)$ between
$p^{\frac{1}{2}}(T,\delta,Z;h_{0},m_{0},\theta_{0})$ and
$p^{\frac{1}{2}}(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})$
pointwisely such that

$\displaystyle
p(T,\delta,Z;h_{0},m_{0},\theta_{0})\left(\log\sqrt{\dfrac{p(T,\delta,Z;h_{0},m_{0},\theta_{0})}{p(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})}}\right)^{2}$
$\displaystyle=p(T,\delta,Z;h_{0},m_{0},\theta_{0})\left(\log\sqrt{p(T,\delta,Z;h_{0},m_{0},\theta_{0})}-\log\sqrt{p(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})}\right)^{2}$
$\displaystyle=\frac{p(T,\delta,Z;h_{0},m_{0},\theta_{0})}{\xi(T,\delta,Z)^{2}}\left(\sqrt{p(T,\delta,Z;h_{0},m_{0},\theta_{0})}-\sqrt{p(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})}\right)^{2}$

Since

$\displaystyle\dfrac{p(T,\delta,Z;h_{0},m_{0},\theta_{0})}{p(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})}=e^{l(T,\delta,Z;h_{0},m_{0},\theta_{0})-l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})}$

by lemma B.1, $l(T,\delta,Z;h_{0},m_{0},\theta_{0})$ and
$l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})$ are bounded
among$[0,\tau]\times\\{0,1\\}\times[-1,1]^{d}$ uniformly on all
$\widehat{\phi}=(\widehat{h},\widehat{m},\widehat{\theta})$. Thus, there exist
constants $C_{1}$ and $C_{2}$ such that $0<C_{1}\leq
p(T,\delta,Z;h_{0},m_{0},\theta_{0})/p(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})\leq
C_{2}$. This leads to the fact that
$p(T,\delta,Z;h_{0},m_{0},\theta_{0})\frac{1}{\xi(T,\delta,Z)^{2}}$ is
bounded. We further obtained that

$\displaystyle
p(T,\delta,Z;h_{0},m_{0},\theta_{0})\left(\log\sqrt{p(T,\delta,Z;h_{0},m_{0},\theta_{0})}-\log\sqrt{p(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})}\right)^{2}$
$\displaystyle\lesssim\left|\sqrt{p(T,\delta,Z;h_{0},m_{0},\theta_{0})}-\sqrt{p(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})}\right|^{2}.$

Thus, we have that

$\displaystyle\sup_{d_{\textsf{PF}}\left[\widehat{\phi},\phi_{0}\right]\leq\varepsilon}\text{Var}(l(T,\delta,Z;h_{0},m_{0},\theta_{0})-l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta}))$
$\displaystyle\lesssim\sup_{d_{\textsf{PF}}\left(\widehat{\phi},\phi_{0}\right)\leq\varepsilon}\mathbb{E}_{Z}\left[\int\left|\sqrt{p(T,\delta,Z;h_{0},m_{0},\theta_{0})}-\sqrt{p(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})}\right|^{2}\mu(dT\times
d\delta)\right]$
$\displaystyle=\sup_{d_{\textsf{PF}}\left(\widehat{\phi},\phi_{0}\right)\leq\varepsilon}d^{2}_{\textsf{PF}}\left(\widehat{\phi},\phi_{0}\right)$
$\displaystyle\leq\varepsilon^{2}.$

#### Step $3$

We define that
$\widetilde{\mathcal{F}}_{n}=\\{l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})-l(T,\delta,Z;\pi_{n}h_{0},\pi_{n}m_{0},\pi_{n}\theta_{0}):\widehat{h}\in\mathcal{H}_{n},\widehat{m}\in\mathcal{M}_{n},\widehat{\theta}\in\Theta\\}$.
Here $(\pi_{n}h_{0},\pi_{n}m_{0},\pi_{n}\theta_{0})$ have been defined in
lemma B.5. Obviously, we have that $\log
N_{[]}(\varepsilon,\widetilde{\mathcal{F}}_{n},\|\cdot\|_{\infty})=\log
N_{[]}(\varepsilon,\mathcal{F}_{n},\|\cdot\|_{\infty})$, where $\mathcal{F}$
is defined in lemma B.9. By lemma B.9, we further have that

$\displaystyle\log
N_{[]}(\varepsilon,\mathcal{F}_{n},\|\cdot\|_{\infty})\lesssim\log\frac{1}{\varepsilon}+\log
N(c_{h}\varepsilon^{1/s_{h}},\mathcal{H}_{n},\|\cdot\|_{2})+\log
N(c_{m}\varepsilon^{1/s_{m}},\mathcal{M}_{n},\|\cdot\|_{2}).$

According to Bartlett et al. [2019, Theorem 7], under condition 4.3, we have
that the VC-dimension of $\mathcal{H}_{n}$ and $\mathcal{M}_{n}$ satisfy that
$\text{VC}\left(\mathcal{H}_{n}\right)\lesssim n^{\frac{1}{\beta+d}}\log^{3}n$
and $\text{VC}\left(\mathcal{M}_{n}\right)\lesssim
n^{\frac{d}{\beta+d}}\log^{3}n$. Thus, we obtain that

$\displaystyle\log
N(c_{h}\varepsilon^{1/s_{h}},\mathcal{H}_{n},\|\cdot\|_{2})\lesssim\frac{\text{VC}\left(\mathcal{H}_{n}\right)}{s_{h}}\log\frac{1}{\varepsilon}\lesssim
n^{\frac{1}{\beta+d}}\log^{3}n\log\frac{1}{\varepsilon},$

and

$\displaystyle\log
N(c_{m}\varepsilon^{1/s_{m}},\mathcal{M}_{n},\|\cdot\|_{2})\lesssim\frac{\text{VC}\left(\mathcal{M}_{n}\right)}{s_{\nu}}\log\frac{1}{\varepsilon}\lesssim
n^{\frac{d}{\beta+d}}\log^{3}n\log\frac{1}{\varepsilon}.$

Thus, we obtain that $\log
N_{[]}(\varepsilon,\widetilde{\mathcal{F}}_{n},\|\cdot\|_{\infty})\lesssim
n^{\frac{d}{\beta+d}}\log^{3}n\log\frac{1}{\varepsilon}$.

#### Step $4$

By the Cauchy-Schwartz inequality, we have that

$\displaystyle\sqrt{\mathbb{E}\left[l(T,\delta,Z;\pi_{n}h_{0},\pi_{n}m_{0},\pi_{n}\theta_{0})-l(T,\delta,Z;h_{0},m_{0},\theta_{0})\right]}$
$\displaystyle\leq\left[\mathbb{E}(l(T,\delta,Z;\pi_{n}h_{0},\pi_{n}m_{0},\pi_{n}\theta_{0})-l(T,\delta,Z;h_{0},m_{0},\theta_{0}))^{2}\right]^{\frac{1}{4}}.$

Similar to the second part and by lemma B.5, we further have that

$\displaystyle\sqrt{\mathbb{E}\left[l(T,\delta,Z;\pi_{n}h_{0},\pi_{n}m_{0},\pi_{n}\theta_{0})-l(T,\delta,Z;h_{0},m_{0},\theta_{0})\right]}\lesssim\sqrt{d_{\textsf{PF}}\left(\pi_{n}\phi_{0},\phi_{0}\right)}\lesssim
n^{-\frac{\beta}{2\beta+2d}}.$

Now let

$\displaystyle\tau=\frac{\beta}{2\beta+2d}-2\frac{\log\log n}{\log n}$

By Step 1,2,3 and Shen and Wong [1994, Theorem 1], we have

$\displaystyle d_{\textsf{PF}}\left(\widehat{\phi}_{n},\phi_{0}\right)$
$\displaystyle=\max\left(n^{-\tau},d_{\textsf{PF}}\left(\pi_{n}\phi_{0},\phi_{0}\right),\right.$
$\displaystyle\left.\sqrt{\mathbb{E}\left[l(T,\delta,Z;\pi_{n}h_{0},\pi_{n}m_{0},\pi_{n}\theta_{0})-l(T,\delta,Z;h_{0},m_{0},\theta_{0})\right]}\right)$

By lemma B.5,
$d_{\textsf{PF}}\left(\pi_{n}\phi_{0},\phi_{0}\right)=O(n^{-\frac{\beta}{\beta+d}})$.

By Step 4,
$\sqrt{\mathbb{E}\left[l(T,\delta,Z;\pi_{n}h_{0},\pi_{n}m_{0},\pi_{n}\theta_{0})-l(T,\delta,Z;h_{0},m_{0},\theta_{0})\right]}=O\left(n^{-\frac{\beta}{2\beta+2d}}\right)$.
Thus, we have
$d_{\textsf{PF}}\left(\widehat{\phi}_{n},\phi_{0}\right)=O(n^{-\frac{\beta}{2\beta+2d}}\log^{2}n)=\widetilde{O}(n^{-\frac{\beta}{2\beta+2d}})$.
∎

###### Proof of theorem 4.7.

The proof is divided into four steps.

#### Step $1$

We denote $\psi_{0}=(\nu_{0},\theta_{0})$ and
$\widehat{\psi}=(\widehat{\nu},\widehat{\theta})$, where
$\widehat{\nu}\in\mathcal{V}_{n}$ and $\widehat{\theta}\in\Theta$. For
arbitrary $0<\varepsilon\leq 1$, we have that

$\displaystyle\inf_{d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)\geq\varepsilon}\mathbb{E}\left[l(T,\delta,Z;\nu_{0},\theta_{0})-l(T,\delta,Z;\widehat{\nu},\widehat{\theta})\right]$
$\displaystyle=\inf_{d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)\geq\varepsilon}\mathbb{E}_{Z}\left[\mathbb{E}_{T,\delta\mid
Z}\left[\log p(T,\delta\mid Z;\nu_{0},\theta_{0})-\log p(T,\delta\mid
Z;\widehat{\nu},\widehat{\theta})\right]\right]$
$\displaystyle=\inf_{d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)\geq\varepsilon}\mathbb{E}_{Z}\left[\text{D}_{\text{KL}}\left(\mathbb{P}_{\widehat{\psi},Z}\parallel\parallel\mathbb{P}_{\psi_{0},Z}\right)\right]$

Using the fact that
$KL(\mathbb{P}_{\widehat{\psi},Z}\parallel\mathbb{P}_{\psi_{0},Z})\geq
2H^{2}(\mathbb{P}_{\widehat{\psi},Z}\parallel\mathbb{P}_{\psi_{0},Z})$. Thus,
we further obtain that

$\displaystyle\inf_{d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)\geq\varepsilon}\mathbb{E}\left[l(T,\delta,Z;\nu_{0},\theta_{0})-l(T,\delta,Z;\widehat{\nu},\widehat{\theta})\right]$
$\displaystyle\geq\inf_{d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)\geq\varepsilon}2\mathbb{E}_{Z}\left[H^{2}(\mathbb{P}_{\widehat{\psi},Z}\parallel\mathbb{P}_{\psi_{0},Z})\right]$
$\displaystyle=2\inf_{d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)\geq\varepsilon}d^{2}_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)$
$\displaystyle\geq 2\varepsilon^{2}.$

#### Step $2$

We consider the following derivation.

$\displaystyle\sup_{d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)\leq\varepsilon}\text{Var}\left[l(T,\delta,Z;\nu_{0},\theta_{0})-l(T,\delta,Z;\widehat{\nu},\widehat{\theta})\right]$
$\displaystyle\leq\sup_{d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)\leq\varepsilon}\mathbb{E}\left[\left(l(T,\delta,Z;\nu_{0},\theta_{0})-l(T,\delta,Z;\widehat{\nu},\widehat{\theta})\right)^{2}\right]$
$\displaystyle=\sup_{d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)\leq\varepsilon}\mathbb{E}_{Z}\left[\mathbb{E}_{T,\delta\mid
Z}\left[\left(\log p(T,\delta,Z;\nu_{0},\theta_{0})-\log
p(T,\delta,Z;\widehat{\nu},\widehat{\theta})\right)^{2}\right]\right]$
$\displaystyle=4\sup_{d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)\leq\varepsilon}\mathbb{E}_{Z}\left[\int\left(p(T,\delta,Z;\nu_{0},\theta_{0})(\log\sqrt{\dfrac{p(T,\delta,Z;\nu_{0},\theta_{0})}{p(T,\delta,Z;\widehat{\nu},\widehat{\theta})}})^{2}\right)\mu(dT\times
d\delta)\right]$

By Taylor’s expansion on $\log x$, there exists $\eta(T,\delta,Z)$ between
$\sqrt{p(T,\delta,Z;\nu_{0},\theta_{0})}$ and
$\sqrt{p(T,\delta,Z;\widehat{\nu},\widehat{\theta})}$ pointwisely such that

$\displaystyle
p(T,\delta,Z;\nu_{0},\theta_{0})(\log\sqrt{\dfrac{p(T,\delta,Z;\nu_{0},\theta_{0})}{p(T,\delta,Z;\widehat{\nu},\widehat{\theta})}})^{2}$
$\displaystyle=p(T,\delta,Z;\nu_{0},\theta_{0})\left(\log\sqrt{p(T,\delta,Z;\nu_{0},\theta_{0})}-\log\sqrt{p(T,\delta,Z;\widehat{\nu},\widehat{\theta})}\right)^{2}$
$\displaystyle=\dfrac{p(T,\delta,Z;\nu_{0},\theta_{0})}{\eta(T,\delta,Z)^{2}}\left(\sqrt{p(T,\delta,Z;\nu_{0},\theta_{0}})-\sqrt{p(T,\delta,Z;\widehat{\nu},\widehat{\theta})}\right)^{2}$

Since
$p(T,\delta,Z;\nu_{0},\theta_{0})/p(T,\delta,Z;\widehat{\nu},\widehat{\theta})=e^{l(T,\delta,Z;\nu_{0},\theta_{0})-l(T,\delta,Z;\widehat{\nu},\widehat{\theta})}$,
by lemma B.2, $l(T,\delta,Z;\nu_{0},\theta_{0})$ and
$l(T,\delta,Z;\widehat{\nu},\widehat{\theta})$ are bounded on
$[0,\tau]\times\\{0,1\\}\times[-1,1]^{d}$ uniformly for all
$\widehat{\psi}=(\widehat{\nu},\widehat{\theta})$. Thus there exist constants
$C_{3}$ and $C_{4}$ such that $0<C_{3}\leq
p(T,\delta,Z;\nu_{0},\theta_{0})/p(T,\delta,Z;\widehat{\nu},\widehat{\theta})\leq
C_{4}$. This leads to the fact that
$p(T,\delta,Z;\nu_{0},\theta_{0})\frac{1}{\eta(T,\delta,Z)^{2}}$ is bounded.
We further have that

$\displaystyle
p(T,\delta,Z;\nu_{0},\theta_{0})\left(\log\sqrt{p(T,\delta,Z;\nu_{0},\theta_{0})}-\log\sqrt{p(T,\delta,Z;\widehat{\nu},\widehat{\theta})}\right)^{2}$
$\displaystyle\lesssim\left|\sqrt{p(T,\delta,Z;\nu_{0},\theta_{0})}-\sqrt{p(T,\delta,Z;\widehat{\nu},\widehat{\theta})}\right|^{2}.$

Thus, we have that

$\displaystyle\sup_{d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)\leq\varepsilon}\text{Var}\left[l(T,\delta,Z;\nu_{0},\theta_{0})-l(T,\delta,Z;\widehat{\nu},\widehat{\theta})\right]$
$\displaystyle\lesssim\sup_{d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)\leq\varepsilon}\mathbb{E}_{Z}\left[\int\left|\sqrt{p(T,\delta,Z;\nu_{0},\theta_{0})}-\sqrt{p(T,\delta,Z;\widehat{\nu},\widehat{\theta})}\right|^{2}\mu(dT\times
d\delta)\right]$
$\displaystyle=\sup_{d_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)\leq\varepsilon}d^{2}_{\textsf{FN}}\left(\widehat{\psi},\psi_{0}\right)$
$\displaystyle\leq\varepsilon^{2}.$

#### Step $3$

We define that
$\widetilde{\mathcal{G}}_{n}=\\{l(T,\delta,Z;\widehat{\nu},\widehat{\theta})-l(T,\delta,Z;\pi_{n}\nu_{0},\pi_{n}\theta_{0}):\widehat{\nu}\in\mathcal{V}_{n},\theta\in\Theta\\}$.
Here $(\pi_{n}\nu_{0},\pi_{n}\theta_{0})$ have been defined in lemma B.6.
Obviously, we have that $\log
N_{[]}(\varepsilon,\widetilde{\mathcal{G}}_{n},\|\cdot\|_{\infty})=\log
N_{[]}(\varepsilon,\mathcal{G}_{n},\|\cdot\|_{\infty})$, where $\mathcal{G}$
is defined in lemma B.10. By lemma B.10, we further obtain that

$\displaystyle\log
N_{[]}(\varepsilon,\mathcal{G}_{n},\|\cdot\|_{\infty})\lesssim\log\frac{1}{\varepsilon}+\log
N(c_{\nu}\varepsilon^{1/s_{\nu}},\mathcal{V}_{n},\|\cdot\|_{2}).$

According to Bartlett et al. [2019, Theorem 7], under condition 4.4, we have
that the VC-dimension of $\mathcal{V}_{n}$ satisfies that
$\text{VC}\left(\mathcal{V}_{n}\right)\lesssim
n^{\frac{d+1}{\beta+d+1}}\log^{3}n$. Thus, we obtain that

$\displaystyle\log
N(c_{h}\varepsilon^{1/s_{\nu}},\mathcal{V}_{n},\|\cdot\|_{2})\lesssim\frac{\text{VC}\left(\mathcal{V}_{n}\right)}{s_{\nu}}\log\frac{1}{\varepsilon}\lesssim
n^{\frac{d+1}{\beta+d+1}}\log^{3}n\log\frac{1}{\varepsilon}.$

Furthermore, we get that $\log
N_{[]}(\varepsilon,\widetilde{\mathcal{G}}_{n},\|\cdot\|_{\infty})\lesssim
n^{\frac{d+1}{\beta+d+1}}\log^{3}n\log\frac{1}{\varepsilon}$.

#### Step $4$

By the Cauchy-Schwartz inequality, we have that

$\displaystyle\sqrt{\mathbb{E}[l(T,\delta,Z;\pi_{n}\nu_{0},\pi_{n}\theta_{0})-l(T,\delta,Z;\nu_{0},\theta_{0})]}\leq\left[\mathbb{E}\left(l(T,\delta,Z;\pi_{n}\nu_{0},\pi_{n}\theta_{0})-l(T,\delta,Z;\nu_{0},\theta_{0})\right)^{2}\right]^{\frac{1}{4}}.$

Similar to the second part and by lemma B.6, we further obtain that

$\displaystyle\sqrt{\mathbb{E}[l(T,\delta,Z;\pi_{n}\nu_{0},\pi_{n}\theta_{0})-l(T,\delta,Z;\nu_{0},\theta_{0})]}\lesssim\sqrt{d_{\textsf{FN}}\left(\pi_{n}\psi_{0},\psi_{0}\right)}\lesssim
n^{-\frac{\beta}{2\beta+2d+2}}$

Now let

$\displaystyle\tau=\frac{\beta}{2\beta+2d+2}-2\frac{\log\log n}{\log n}.$

By step 1,2,3 and Step 1,2,3 and Shen and Wong [1994, Theorem 1],

$\displaystyle d_{\textsf{FN}}\left(\widehat{\psi}_{n},\psi_{0}\right)$
$\displaystyle=\max\left(n^{-\tau},d_{\textsf{FN}}\left(\pi_{n}\psi_{0},\psi_{0}\right),\right.$
$\displaystyle\left.\sqrt{\mathbb{E}[l(T,\delta,Z;\pi_{n}\nu_{0},\pi_{n}\theta_{0})-l(T,\delta,Z;\nu_{0},\theta_{0})]}\right)$

By lemmaB.6,
$d_{\textsf{FN}}\left(\pi_{n}\psi_{0},\psi_{0}\right)=O(n^{-\frac{\beta}{\beta+d+1}})$

By Step 4,
$\sqrt{\mathbb{E}[l(T,\delta,Z;\pi_{n}\nu_{0},\pi_{n}\theta_{0})-l(T,\delta,Z;\nu_{0},\theta_{0})]}=O(n^{-\frac{\beta}{2\beta+2d+2}})$.
Thus, we have
$d_{\textsf{FN}}\left(\widehat{\psi}_{n},\psi_{0}\right)=O(n^{-\frac{\beta}{2\beta+2d+2}}\log^{2}n)=\widetilde{O}(n^{-\frac{\beta}{2\beta+2d+2}})$.
∎

###  B.4 Proofs of technical lemmas

###### Proof of lemma B.1.

Since $h_{0}(T)\in\mathcal{W}^{\beta}_{M}([0,\tau])$ and
$m_{0}(Z)\in\mathcal{W}^{\beta}_{M}([-1,1]^{d})$, we have that $h_{0}(T)\leq
M$, $m_{0}(Z)\leq M$ and $e^{m_{0}(Z)}\int_{0}^{T}h_{0}(s)ds\leq\tau e^{2M}$.
Let $\mathcal{B}=[0,\tau e^{2M}]$, we have that

$\displaystyle|l(T,\delta,Z;h_{0},m_{0},\theta_{0})|$
$\displaystyle\leq\left|\log
g_{\theta_{0}}\left(e^{m_{0}(Z)}\int_{0}^{T}e^{h_{0}(s)}ds\right)\right|+|h_{0}(T)|+|m_{0}(Z)|+\left|G_{\theta_{0}}\left(e^{m_{0}(Z)}\int_{0}^{T}e^{h_{0}(s)}ds\right)\right|$
$\displaystyle\leq 2M+\sup_{x\in\mathcal{B}}\left|\log
g_{\theta_{0}}(x)\right|+\sup_{x\in\mathcal{B}}\left|G_{\theta_{0}}(x)\right|$

By condition 5, we have that $l(T,\delta,Z;h_{0},m_{0},\theta_{0})$ is bounded
for $(T,\delta,Z)\in[0,\tau]\times\\{0,1\\}\times[-1,1]^{d}$. The proof of the
boundness of $l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})$ is
similar. ∎

###### Proof of lemma B.2.

Since $\nu_{0}(T,Z)\in\mathcal{W}^{\beta}_{M}([0,\tau]\times[-1,1]^{d})$, we
have $\nu_{0}(T,Z)\leq M$ and $\int_{0}^{T}e^{\nu(s,Z)}ds\leq\tau e^{M}$. Let
$\mathcal{B}=[0,\tau e^{M}]$, we have that

$\displaystyle|l(T,\delta,Z;\nu_{0},\theta_{0})|$ $\displaystyle\leq\left|\log
G^{\prime}_{\theta_{0}}\left(\int_{0}^{T}e^{\nu_{0}(s,Z)}ds\right)\right|+|\nu_{0}(T,Z)|+\left|G_{\theta_{0}}\left(\int_{0}^{T}e^{\nu_{0}(s,Z)}ds\right)\right|$
$\displaystyle\leq M+\sup_{x\in\mathcal{B}}\left|\log
G^{\prime}_{\theta_{0}}(x)\right|+\sup_{x\in\mathcal{B}}\left|G_{\theta_{0}}(x)\right|.$

By condition 4.5, we have that $l(T,\delta,Z;\nu_{0},\theta_{0})$ is bounded
among $(T,\delta,Z)\in[0,\tau]\times\\{0,1\\}\times[-1,1]^{d}$ The proof of
the boundness of $l(T,\delta,Z;\widehat{\nu},\widehat{\theta})$ is similar. ∎

###### Proof of lemma B.3.

By definition, we have that

$\displaystyle|l(T,\delta,Z;h_{0},m_{0},\theta_{0})-l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})|$
$\displaystyle\leq\left|\log
g_{\theta_{0}}\left(e^{m_{0}(Z)}\int_{0}^{T}e^{h_{0}(s)}ds\right)-\log
g_{\widehat{\theta}}\left(e^{\widehat{m}(Z)}\int_{0}^{T}e^{\widehat{h}(s)}ds\right)\right|+\left|h_{0}(T)-\widehat{h}(T)\right|$
$\displaystyle+|m_{0}(Z)-\widehat{m}(Z)|+\left|G_{\theta_{0}}\left(e^{m_{0}(Z)}\int_{0}^{T}e^{h_{0}(s)}ds\right)-G_{\widehat{\theta}}\left(e^{\widehat{m}(Z)}\int_{0}^{T}e^{\widehat{h}(s)}ds\right)\right|.$

Let $\mathcal{B}=[0,\tau\max(e^{2M},e^{M_{h}+M_{m}})]$. By Taylor’s expansion,
we can further show that

$\displaystyle|l(T,\delta,Z;h_{0},m_{0},\theta_{0})-l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})|$
$\displaystyle\leq\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial\log
g_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{\theta}}}\right|\cdot\left|\theta_{0}-\widehat{\theta}\right|$
$\displaystyle+\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial\log
g_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{x}}}\right|\cdot\left|e^{m_{0}(Z)}\int_{0}^{T}e^{h_{0}(s)}ds-e^{\widehat{m}(Z)}\int_{0}^{T}e^{\widehat{h}(s)}ds\right|$
$\displaystyle+|h_{0}(T)-\widehat{h}(T)|+|m_{0}(Z)-\widehat{m}(Z)|+\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial
G_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{\theta}}}\right|\cdot\left|\theta_{0}-\widehat{\theta}\right|$
$\displaystyle+\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial
G_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{x}}}\right|\cdot\left|e^{m_{0}(Z)}\int_{0}^{T}e^{h_{0}(s)}ds-e^{\widehat{m}(Z)}\int_{0}^{T}e^{\widehat{h}(s)}ds\right|.$

Again, by Taylor’s expansion, we have that

$\displaystyle\left|e^{m_{0}(Z)}\int_{0}^{T}e^{h_{0}(s)}ds-e^{\widehat{m}(Z)}\int_{0}^{T}e^{\widehat{h}(s)}ds\right|$
$\displaystyle\leq\left|e^{m_{0}(Z)}\int_{0}^{T}(e^{h_{0}(s)}-e^{\widehat{h}(s)})ds\right|+\left|(e^{m_{0}(Z)}-e^{\widehat{m}(Z)})\int_{0}^{T}e^{\widehat{h}(s)}ds\right|$
$\displaystyle\leq e^{M}\cdot\tau
e^{\max(M,M_{h})}\left\|h_{0}-\widehat{h}\right\|_{\infty}+\tau e^{M_{h}}\cdot
e^{\max(M,M_{m})}\|m_{0}-\widehat{m}\|_{\infty}.$

Finally, we obtain that

$\displaystyle|l(T,\delta,Z;h_{0},m_{0},\theta_{0})-l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})|$
$\displaystyle\leq\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial\log
g_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{x}}}\right|\cdot\left[e^{M}\cdot\tau
e^{\max(M,M_{h})}\|h_{0}-\widehat{h}\|_{\infty}+\tau e^{M_{h}}\cdot
e^{\max(M,M_{m})}\|m_{0}-\widehat{m}\|_{\infty}\right]$
$\displaystyle+\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial\log
g_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{\theta}}}\right|\cdot\left|\theta_{0}-\widehat{\theta}\right|+\left|h_{0}(T)-\widehat{h}(T)\right|+|m_{0}(Z)-\widehat{m}(Z)|$
$\displaystyle+\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial
G_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{\theta}}}\right|\cdot\left|\theta_{0}-\widehat{\theta}\right|$
$\displaystyle+\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial
G_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{x}}}\right|\cdot\left[e^{M}\cdot\tau
e^{\max(M,M_{h})}\|h_{0}-\widehat{h}\|_{\infty}+\tau e^{M_{h}}\cdot
e^{\max(M,M_{m})}\|m_{0}-\widehat{m}\|_{\infty}\right].$

Taking supremum on both sides, we conclude that

$\|l(T,\delta,Z;h_{0},m_{0},\theta_{0})-l(T,\delta,Z;\widehat{h},\widehat{m},\widehat{\theta})\|_{\infty}\lesssim|\theta_{0}-\widehat{\theta}|+\|h_{0}-\widehat{h}\|_{\infty}+\|m_{0}-\widehat{m}\|_{\infty}.$

The proof of the second inequality is similar. ∎

###### Proof of lemma B.4.

By definition, we have that

$\displaystyle|l(T,\delta,Z;\nu_{0},\theta_{0})-l(T,\delta,Z;\widehat{\nu},\widehat{\theta})|$
$\displaystyle\leq\left|\log
g_{\theta_{0}}\left(\int_{0}^{T}e^{\nu_{0}(s,Z)}ds\right)-\log
g_{\widehat{\theta}}\left(\int_{0}^{T}e^{\widehat{\nu}(s,Z)}ds\right)\right|+\left|\nu_{0}(T,Z)-\widehat{\nu}(T,Z)\right|$
$\displaystyle+\left|G_{\theta_{0}}\left(\int_{0}^{T}e^{\nu_{0}(s,Z)}ds\right)-G_{\widehat{\theta}}\left(\int_{0}^{T}e^{\widehat{\nu}(s,Z)}ds\right)\right|.$

Let $\mathcal{B}=[0,\tau\max(e^{M},e^{M_{\nu}})]$. By Taylor’s expansion, we
can further show that

$\displaystyle|l(T,\delta,Z;\nu_{0},\theta_{0})-l(T,\delta,Z;\widehat{\nu},\widehat{\theta})|$
$\displaystyle\leq\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial\log
g_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{\theta}}}\right|\cdot\left|\theta_{0}-\widehat{\theta}\right|+\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial\log
g_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{x}}}\right|\cdot\left|\int_{0}^{T}e^{\nu_{0}(s,Z)}ds-\int_{0}^{T}e^{\widehat{\nu}(s,Z)}ds\right|$
$\displaystyle+|\nu_{0}(T,Z)-\widehat{\nu}(T,Z)|+\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial
G_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{\theta}}}\right|\cdot|\theta_{0}-\widehat{\theta}|$
$\displaystyle+\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial
G_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{x}}}\right|\cdot\left|\int_{0}^{T}e^{\nu_{0}(s,Z)}ds-\int_{0}^{T}e^{\widehat{\nu}(s,Z)}ds\right|.$

Again, by Taylor’s expansion,

$\displaystyle\left|\int_{0}^{T}e^{\nu_{0}(s,Z)}ds-\int_{0}^{T}e^{\widehat{\nu}(s,Z)}ds\right|\leq\tau
e^{\max(M,M_{\nu})}\|\nu_{0}-\widehat{\nu}\|_{\infty},$

Finally, we obtain that

$\displaystyle\left|l(T,\delta,Z;\nu_{0},\theta_{0})-l(T,\delta,Z;\widehat{\nu},\widehat{\theta})\right|$
$\displaystyle\leq\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial\log
g_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{\theta}}}\right|\cdot\left|\theta_{0}-\widehat{\theta}\right|+\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial\log
g_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{x}}}\right|\cdot\tau
e^{\max(M,M_{\nu})}\left\|\nu_{0}-\widehat{\nu}\right\|_{\infty}$
$\displaystyle+|\nu_{0}(T,Z)-\widehat{\nu}(T,Z)|+\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial
G_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{\theta}}}\right|\cdot\left|\theta_{0}-\widehat{\theta}\right|$
$\displaystyle+\sup_{\tilde{\theta}\in\Theta,\tilde{x}\in\mathcal{B}}\left|\frac{\partial
G_{\tilde{\theta}}(\tilde{x})}{\partial_{\tilde{x}}}\right|\cdot\tau
e^{\max(M,M_{\nu})}\left\|\nu_{0}-\widehat{\nu}\right\|_{\infty}.$

Taking supremum on both sides, we conclude that

$\|l(T,\delta,Z;\nu_{0},\theta_{0})-l(T,\delta,Z;\widehat{\nu},\widehat{\theta})\|_{\infty}\lesssim|\theta_{0}-\widehat{\theta}|+\|\nu_{0}-\widehat{\nu}\|_{\infty},$

The proof of the second inequality is similar. ∎

###### Proof of lemma B.5.

According to Yarotsky [2017, Theorem 1], there exist approximating functions
$\widehat{h}^{*}$ and $\widehat{m}^{*}$ such that
$\|\widehat{h}^{*}-h_{0}\|_{\infty}=O\left(n^{-\frac{\beta}{\beta+d}}\right)$
and
$\|\widehat{m}^{*}-m_{0}\|_{\infty}=O\left(n^{-\frac{\beta}{\beta+d}}\right)$.
Let $\pi_{n}h_{0}=\widehat{h}^{*}$, $\pi_{n}m_{0}=\widehat{m}^{*}$, and
$\pi_{n}\theta=\theta_{0}$. We have that

$\displaystyle d_{\textsf{PF}}\left(\pi_{n}\phi_{0},\phi_{0}\right)$
$\displaystyle=\sqrt{\mathbb{E}_{Z}\left[\int|\sqrt{p(T,\delta\mid
Z;\pi_{n}h_{0},\pi_{n}m_{0},\pi_{n}\theta_{0})}-\sqrt{p(T,\delta\mid
Z;h_{0},m_{0},\theta_{0})}|^{2}\mu(dT\times d\delta)\right]}$
$\displaystyle=\sqrt{\mathbb{E}_{Z}\left[\int[e^{\frac{1}{2}l(T,\delta,Z;\pi_{n}h_{0},\pi_{n}m_{0},\pi_{n}\theta_{0})}-e^{\frac{1}{2}l(T,\delta,Z;h_{0},m_{0},\theta_{0})}]^{2}f_{C\mid
Z}(T)^{1-\delta}S_{C\mid Z}(T)^{\delta}\mu(dT\times d\delta)\right]}$
$\displaystyle\leq\left\|e^{\frac{1}{2}l(T,\delta,Z;\pi_{n}h_{0},\pi_{n}m_{0},\pi_{n}\theta_{0})}-e^{\frac{1}{2}l(T,\delta,Z;h_{0},m_{0},\theta_{0})}\right\|_{\infty}$
$\displaystyle\qquad\qquad\times\sqrt{\mathbb{E}_{Z}\left[\int f_{C\mid
Z}(T)^{1-\delta}S_{C\mid Z}(T)^{\delta}\mu(dT\times d\delta)\right]}.$

By lemma B.1 and B.3, we have that

$\displaystyle\|e^{\frac{1}{2}l(T,\delta,Z;\pi_{n}h_{0},\pi_{n}m_{0},\pi_{n}\theta_{0})}-e^{\frac{1}{2}l(T,\delta,Z;h_{0},m_{0},\theta_{0})}\|_{\infty}$
$\displaystyle\lesssim\|\pi_{n}\theta_{0}-\theta_{0}\|+\|\pi_{n}h_{0}-h_{0}\|_{\infty}+\|\pi_{n}m_{0}-m_{0}\|_{\infty}$
$\displaystyle=O\left(n^{-\frac{\beta}{\beta+d}}\right).$

Since $f_{C\mid Z}(T)^{1-\delta}\leq f_{C\mid Z}(T)+1$ and $S_{C\mid
Z}(T)^{\delta}\leq 1$, we also have that

$\displaystyle\sqrt{\mathbb{E}_{Z}\left[\int f_{C\mid Z}(T)^{1-\delta}S_{C\mid
Z}(T)^{\delta}\mu(dT\times d\delta)\right]}$ $\displaystyle\leq$
$\displaystyle\sqrt{\mathbb{E}_{Z}\left[\int(1+f_{C\mid Z}(T))\mu(dT\times
d\delta)\right]}$ $\displaystyle\leq$ $\displaystyle\sqrt{2+2\tau}.$

Thus, we obtain that
$d_{\textsf{PF}}\left(\pi_{n}\phi_{0},\phi_{0}\right)=O\left(n^{-\frac{\beta}{\beta+d}}\right)$.
∎

###### Proof of lemma B.6.

According to Yarotsky [2017, Theorem 1], there exists an approximating
function $\widehat{\nu}^{*}$ such that
$\|\widehat{\nu}^{*}-\nu_{0}\|_{\infty}=O\left(n^{-\frac{\beta}{\beta+d+1}}\right)$.
Let $\pi_{n}\nu_{0}=\widehat{\nu}^{*}$ and $\pi_{n}\theta_{0}=\theta_{0}$. We
have that

$\displaystyle d_{\textsf{FN}}\left(\pi_{n}\psi_{0},\psi_{0}\right)$
$\displaystyle=\sqrt{\mathbb{E}_{Z}\left[\int\left|\sqrt{p(T,\delta\mid
Z;\pi_{n}\nu_{0},\pi_{n}\theta_{0})}-\sqrt{p(T,\delta\mid
Z;\nu_{0},\theta_{0})}\right|^{2}\mu(dT\times d\delta)\right]}$
$\displaystyle=\sqrt{\mathbb{E}_{Z}\left[\int\left[e^{\frac{1}{2}l(T,\delta,Z;\pi_{n}\nu_{0},\pi_{n}\theta_{0})}-e^{\frac{1}{2}l(T,\delta,Z;\nu_{0},\theta_{0})}\right]^{2}f_{C\mid
Z}(T)^{1-\delta}S_{C\mid Z}(T)^{\delta}\mu(dT\times d\delta)\right]}$
$\displaystyle\leq\left\|\frac{1}{2}e^{l(T,\delta,Z;\pi_{n}\nu_{0},\pi_{n}\theta_{0})}-\frac{1}{2}e^{l(T,\delta,Z;\nu_{0},\theta_{0})}\right\|_{\infty}\sqrt{\mathbb{E}_{Z}\left[\int
f_{C\mid Z}(T)^{1-\delta}S_{C\mid Z}(T)^{\delta}\mu(dT\times
d\delta)\right]}.$

By lemma B.2 and B.4, we have that

$\displaystyle\left\|e^{\frac{1}{2}l(T,\delta,Z;\pi_{n}\nu_{0},\pi_{n}\theta_{0})}-e^{\frac{1}{2}l(T,\delta,Z;\nu_{0},\theta_{0})}\right\|_{\infty}$
$\displaystyle\lesssim\|\pi_{n}\theta_{0}-\theta_{0}\|+\|\pi_{n}\nu_{0}-\nu_{0}\|_{\infty}$
$\displaystyle=O\left(n^{-\frac{\beta}{\beta+d+1}}\right).$

Since $f_{C\mid Z}(T)^{1-\delta}\leq f_{C\mid Z}(T)+1$ and $S_{C\mid
Z}(T)^{\delta}\leq 1$, we also have that

$\displaystyle\sqrt{\mathbb{E}_{Z}\left[\int f_{C\mid Z}(T)^{1-\delta}S_{C\mid
Z}(T)^{\delta}\mu(dT\times d\delta)\right]}$ $\displaystyle\leq$
$\displaystyle\sqrt{\mathbb{E}_{Z}\left[\int(1+f_{C\mid Z}(T))\mu(dT\times
d\delta)\right]}$ $\displaystyle\leq$ $\displaystyle\sqrt{2+2\tau}.$

Thus, we obtain that
$d_{\textsf{FN}}\left(\pi_{n}\psi_{0},\psi_{0}\right)=O\left(n^{-\frac{\beta}{\beta+d+1}}\right)$.
∎

###### Proof of lemma B.7.

The left inequality is trivial according to the definition of covering number.
We need to show that the correctness of the right inequality.

Suppose that we have $\\{B(g_{i},\frac{\varepsilon}{2})\\},i=1\ldots,N$, where
$N=N(\frac{\varepsilon}{2},\mathcal{F},\|\cdot\|)$, are the minimal number of
$\frac{\varepsilon}{2}$-ball that covers $\mathcal{F}$. Then there exists at
least one $f_{i}\in\mathcal{F}$ such that $f_{i}\in B(g_{i},\varepsilon)$.
Consider the following $\varepsilon-balls$
$\\{B(f_{i},\varepsilon)\\},i=1\ldots,N$. For arbitrary $f\in\mathcal{F}\cap
B(g_{i},\frac{\varepsilon}{2}),$ we have that
$\|f-f_{i}\|\leq\|f-g_{i}\|+\|f_{i}-g_{i}\|\leq\varepsilon$. Thus
$\\{B(f_{i},\varepsilon)\\},i=1\ldots,N$ forms a $\varepsilon$-covering of
$\mathcal{F}$. By definition, we have that
$\widetilde{N}(\varepsilon,\mathcal{F},\|\cdot\|)\leq
N(\frac{\varepsilon}{2},\mathcal{F},\|\cdot\|)$. ∎

###### Proof of lemma B.8.

The proof of the first two inequalities follows exactly the same steps of
lemma B.7. Here we just need to mention the rest of the statement that
$\widetilde{N}_{[]}(\varepsilon,\mathcal{F},\|\cdot\|_{\infty})=\widetilde{N}(\frac{\varepsilon}{2},\mathcal{F},\|\cdot\|_{\infty})$.
We first choose a set of $\frac{\varepsilon}{2}$-covering balls
$\\{B(f_{i},\frac{\varepsilon}{2})\\},i=1,\ldots,N_{1}$, where
$N_{1}=\widetilde{N}(\frac{\varepsilon}{2},\mathcal{F},\|\cdot\|_{\infty})$.
Now we construct a set of brackets $\\{[l_{i},u_{i}]\\},i=1\ldots,N_{1}$,
where $l_{i}=f_{i}-\frac{\varepsilon}{2}$ and
$u_{i}=f_{i}+\frac{\varepsilon}{2}$. Noting that the bracket
$\\{[l_{i},u_{i}]\\}$ is exactly the same as $B(f_{i},\frac{\varepsilon}{2})$,
The set $\\{[l_{i},u_{i}]\\},i=1,\ldots,N_{1}$ covers $\mathcal{F}$, which
leads to
$\widetilde{N}_{[]}(\varepsilon,\mathcal{F},\|\cdot\|_{\infty})\leq\widetilde{N}(\frac{\varepsilon}{2},\mathcal{F},\|\cdot\|_{\infty})$.
Likewise, we have that
$\widetilde{N}_{[]}(\varepsilon,\mathcal{F},\|\cdot\|_{\infty})\geq\widetilde{N}(\frac{\varepsilon}{2},\mathcal{F},\|\cdot\|_{\infty})$.
Consequently, we have that
$\widetilde{N}_{[]}(\varepsilon,\mathcal{F},\|\cdot\|_{\infty})=\widetilde{N}(\frac{\varepsilon}{2},\mathcal{F},\|\cdot\|_{\infty})$.
∎

###### Proof of lemma B.9.

By lemma B.8, first we have that
$N_{[]}(\varepsilon,\mathcal{F}_{n},\|\cdot\|_{\infty})\leq\widetilde{N}_{[]}(\varepsilon,\mathcal{F}_{n},\|\cdot\|_{\infty})$.
By lemma B.3, there exists a constant $c_{1}>0$ such that for arbitrary
$\widehat{h}_{1},\widehat{h}_{2}\in\mathcal{H}_{n}$,$\widehat{m}_{1},\widehat{m}_{2}\in\mathcal{M}_{n}$
and $\widehat{\theta}_{1},\widehat{\theta}_{2}\in\Theta$, we have that

$\displaystyle\|l(T,\delta,Z;\widehat{h}_{1},\widehat{m}_{1},\theta_{1})-l(T,\delta,Z;\widehat{h}_{2},\widehat{m}_{2},\theta_{2})\|_{\infty}\leq
c_{1}[|\widehat{\theta}_{1}-\widehat{\theta}_{2}|+\|\widehat{h}_{1}-\widehat{h}_{2}\|_{\infty}+\|\widehat{m}_{1}-\widehat{m}_{2}\|_{\infty}],$

which indicates that as long as
$|\widehat{\theta}_{1}-\widehat{\theta}_{2}|\leq\frac{\varepsilon}{3c_{1}}$,
$\|\widehat{h}_{1}-\widehat{h}_{2}\|_{\infty}\leq\frac{\varepsilon}{3c_{1}}$
and
$\|\widehat{m}_{1}-\widehat{m}_{2}\|_{\infty}\leq\frac{\varepsilon}{3c_{1}}$,
we have that
$\|l(T,\delta,Z;\widehat{h}_{1},\widehat{m}_{1},\theta_{1})-l(T,\delta,Z;\widehat{h}_{2},\widehat{m}_{2},\theta_{2})\|_{\infty}\leq\varepsilon$.
Consequently, we have that

$\displaystyle\widetilde{N}_{[]}(\varepsilon,\mathcal{F}_{n},\|\cdot\|_{\infty})\leq\widetilde{N}_{[]}(\frac{\varepsilon}{3c_{1}},\Theta,\|\cdot\|_{\infty})\times\widetilde{N}_{[]}(\frac{\varepsilon}{3c_{1}},\mathcal{H}_{n},\|\cdot\|_{\infty})\times\widetilde{N}_{[]}(\frac{\varepsilon}{3c_{1}},\mathcal{M}_{n},\|\cdot\|_{\infty}).$

Since $\Theta$ is a compact set on $\mathbb{R}$, by lemma B.8 and traditional
volume argument, we have that
$\widetilde{N}_{[]}(\frac{\varepsilon}{3c_{1}},\Theta,\|\cdot\|_{\infty})\leq
N_{[]}(\frac{\varepsilon}{6c_{1}},\Theta,\|\cdot\|_{\infty})\lesssim\frac{1}{\varepsilon}$.

For
$\widetilde{N}_{[]}(\frac{\varepsilon}{3c_{1}},\mathcal{H}_{n},\|\cdot\|_{\infty})$,
by lemma B.8, we have that
$\widetilde{N}_{[]}(\frac{\varepsilon}{3c_{1}},\mathcal{H}_{n},\|\cdot\|_{\infty})=\widetilde{N}(\frac{\varepsilon}{3c_{1}},\mathcal{H}_{n},\|\cdot\|_{\infty})$.
By Chen and Shen [1998, Lemma 2], there exists a constant $c_{2}>0$ such that
$\|\widehat{h}_{1}-\widehat{h}_{2}\|_{\infty}\leq
c_{2}\|\widehat{h}_{1}-\widehat{h}_{2}\|_{2}^{s_{h}}$, which leads to
$\widetilde{N}(\frac{\varepsilon}{3c_{1}},\mathcal{H}_{n},\|\cdot\|_{\infty})\leq\widetilde{N}(\frac{\varepsilon^{1/s_{h}}}{(3c_{1}c_{2})^{1/s_{h}}},\mathcal{H}_{n},\|\cdot\|_{2})$.
By lemma B.7 we further have that
$\widetilde{N}(\frac{\varepsilon^{1/s_{h}}}{(3c_{1}c_{2})^{1/s_{h}}},\mathcal{H}_{n},\|\cdot\|_{2})\leq
N(\frac{\varepsilon^{1/s_{h}}}{2(3c_{1}c_{2})^{1/s_{h}}},\mathcal{H}_{n},\|\cdot\|_{2})$.
Let $c_{h}=\frac{1}{2(3c_{1}c_{2})^{1/s_{h}}}$. We have that
$\widetilde{N}_{[]}(\frac{\varepsilon}{3c_{1}},\mathcal{H}_{n},\|\cdot\|_{\infty})\leq
N(c_{h}\varepsilon^{1/s_{h}},\mathcal{H}_{n},\|\cdot\|_{2})$.

Similarly, there exists a constant $c_{m}>0$ such that
$\widetilde{N}_{[]}(\frac{\varepsilon}{3c_{1}},\mathcal{M}_{n},\|\cdot\|_{\infty})\leq
N(c_{m}\varepsilon^{1/s_{m}},\mathcal{M}_{n},\|\cdot\|_{2})$.

Thus, finally we can obtain that

$\displaystyle
N_{[]}(\varepsilon,\mathcal{F}_{n},\|\cdot\|_{\infty})\lesssim\frac{1}{\varepsilon}N(c_{h}\varepsilon^{1/s_{h}},\mathcal{H}_{n},\|\cdot\|_{2})\times
N(c_{m}\varepsilon^{1/s_{m}},\mathcal{M}_{n},\|\cdot\|_{2}).$

∎

###### Proof of lemma B.10.

By lemma B.8, first we have
$N_{[]}(\varepsilon,\mathcal{G}_{n},\|\cdot\|_{\infty})\leq\widetilde{N}_{[]}(\varepsilon,\mathcal{G}_{n},\|\cdot\|_{\infty})$.
By lemma B.4, there exists a constant $c_{3}>0$ such that for arbitrary
$\widehat{\nu}_{1},\widehat{\nu}_{2}\in\mathcal{V}_{n}$ and
$\widehat{\theta}_{1},\widehat{\theta}_{2}\in\Theta$, we have that

$\displaystyle\|l(T,\delta,Z;\widehat{\nu}_{1},\widehat{\theta}_{1})-l(T,\delta,Z;\widehat{\nu}_{2},\widehat{\theta}_{2})\|_{\infty}\leq
c_{3}[|\widehat{\theta}_{1}-\widehat{\theta}_{2}|+\|\widehat{\nu}_{1}-\widehat{\nu}_{2}\|_{\infty}],$

which indicates that as long as
$|\widehat{\theta}_{1}-\widehat{\theta}_{2}|\leq\frac{\varepsilon}{2c_{3}}$and
$\|\widehat{\nu}_{1}-\widehat{\nu}_{2}\|_{\infty}\leq\frac{\varepsilon}{2c_{3}}$,
we have that
$\|l(T,\delta,Z;\widehat{\nu}_{1},\widehat{\theta}_{1})-l(T,\delta,Z;\widehat{\nu}_{2},\widehat{\theta}_{2})\|_{\infty}\leq\varepsilon$.
Thus, we have:

$\displaystyle\widetilde{N}_{[]}(\varepsilon,\mathcal{G}_{n},\|\cdot\|_{\infty})\leq\widetilde{N}_{[]}(\frac{\varepsilon}{2c_{3}},\Theta,\|\cdot\|_{\infty})\times\widetilde{N}_{[]}(\frac{\varepsilon}{2c_{3}},\mathcal{V}_{n},\|\cdot\|_{\infty}).$

Since $\Theta$ is a compact set on $\mathbb{R}$, by lemma B.8 and traditional
volume argument, we have that
$\widetilde{N}_{[]}(\frac{\varepsilon}{2c_{3}},\Theta,\|\cdot\|_{\infty})\leq
N_{[]}(\frac{\varepsilon}{4c_{3}},\Theta,\|\cdot\|_{\infty})\lesssim\frac{1}{\varepsilon}$.

For
$\widetilde{N}_{[]}(\frac{\varepsilon}{2c_{3}},\mathcal{V}_{n},\|\cdot\|_{\infty})$,
by lemma B.8, we have that
$\widetilde{N}_{[]}(\frac{\varepsilon}{2c_{3}},\mathcal{V}_{n},\|\cdot\|_{\infty})=\widetilde{N}(\frac{\varepsilon}{2c_{3}},\mathcal{V}_{n},\|\cdot\|_{\infty})$.
By Chen and Shen [1998, Lemma 2], there exists a constant $c_{4}>0$ such that
$\|\widehat{\nu}_{1}-\widehat{\nu}_{2}\|_{\infty}\leq
c_{4}\|\widehat{\nu}_{1}-\widehat{\nu}_{2}\|_{2}^{s_{h}}$, which leads to
$\widetilde{N}(\frac{\varepsilon}{2c_{3}},\mathcal{V}_{n},\|\cdot\|_{\infty})\leq\widetilde{N}(\frac{\varepsilon^{1/s_{\nu}}}{(2c_{3}c_{4})^{1/s_{\nu}}},\mathcal{V}_{n},\|\cdot\|_{2})$.
By lemma B.7 we further have
$\widetilde{N}(\frac{\varepsilon^{1/s_{\nu}}}{(2c_{3}c_{4})^{1/s_{\nu}}},\mathcal{V}_{n},\|\cdot\|_{2})\leq
N(\frac{\varepsilon^{1/s_{\nu}}}{2(2c_{3}c_{4})^{1/s_{\nu}}},\mathcal{V}_{n},\|\cdot\|_{2})$.
Let $c_{\nu}=\frac{1}{2(2c_{3}c_{4})^{1/s_{\nu}}}$, we have that
$\widetilde{N}_{[]}(\frac{\varepsilon}{2c_{3}},\mathcal{V}_{n},\|\cdot\|_{\infty})\leq
N(c_{\nu}\varepsilon^{1/s_{\nu}},\mathcal{V}_{n},\|\cdot\|_{2})$.

Thus, finally we can obtain that

$\displaystyle
N_{[]}(\varepsilon,\mathcal{G}_{n},\|\cdot\|_{\infty})\lesssim\frac{1}{\varepsilon}N(c_{\nu}\varepsilon^{1/s_{\nu}},\mathcal{V}_{n},\|\cdot\|_{2}).$

∎

##  Appendix C Experimental details

###  C.1 Dataset summary

We report summaries of descriptive statistics of the $6$ benchmark datasets
used in section 5.2 in table 3.

Table 3: Descriptive statistics of benchmark datasets | METABRIC | RotGBSG | FLCHAIN | SUPPORT | MIMIC-III | KKBOX
---|---|---|---|---|---|---
Size | $1904$ | $2232$ | $6524$ | $8873$ | $35953$ | $2646746$
Censoring rate | $0.423$ | $0.432$ | $0.699$ | $0.320$ | $0.901$ | $0.280$
Features | $9$ | $7$ | $8$ | $14$ | $26$ | $15$

###  C.2 Details of synthetic experiments

Since the true model is assumed to be of PF form, we generate event time
according to the following transformed regression model Dabrowska and Doksum
[1988]:

$\displaystyle\log H(\tilde{T})=-m(Z)+\epsilon,$ (15)

where $H(t)=\int_{0}^{t}e^{h(s)}ds$ with $h$ defined in (2). The error term
$\epsilon$ is generated such that $e^{\epsilon}$ has cumulative hazard
function $G_{\theta}$. The formulation (15) is the equivalent to (2) Dabrowska
and Doksum [1988], Cuzick [1988], Kosorok et al. [2004]. In our experiments,
the covariates are of dimension $5$, sampled independently from the uniform
distribution over $[0,1]$. We set $h(t)=t$ and hence $H(t)=e^{t}$. The
function form of $m(Z)$ is set to be $m(Z)=\sin(\langle
Z,\beta\rangle)+\langle\sin(Z),\beta\rangle$, where
$\beta=(0.1,0.2,0.3,0.4,0.5)$. Then censoring time $C$ is generated according
to

$\displaystyle\log H(C)=-m(Z)+\epsilon_{C},$ (16)

which reuses covariate $Z$, and draws independently a noise vector
$\epsilon_{C}$ such that the censoring ratio is controlled at around $40\%$.
We generate three datasets with $n\in\\{1000,5000,10000\\}$ respectively.

Hyperparameter configurations We specify below the network architectures and
optimization configurations used in all the tasks:

PF scheme:

For both $\widehat{m}$ and $\widehat{h}$, we use $64$ hidden units for
$n=1000$, $128$ hidden units for $n=5000$ and $256$ hidden units for
$n=10000$. We train each model for $100$ epochs with batch size $128$,
optimized using Adam with learning rate $0.0001$, and no weight decay.

FN scheme:

For both $\widehat{\nu}$, we use $64$ hidden units for $n=1000$, $128$ hidden
units for $n=5000$ and $256$ hidden units for $n=10000$. We train each model
for $100$ epochs with batch size $128$, optimized using Adam with learning
rate $0.0001$, and no weight decay.

###  C.3 Details of public data experiments

Dataset preprocessing For METABRIC, RotGBSG, FLCHAIN, SUPPORT and KKBOX
dataset, we take the version provided in the pycox package Kvamme et al.
[2019]. We standardize continuous features into zero mean and unit variance
and do one-hot encodings for all categorical features. For the MIMIC-III
dataset, we follow the preprocessing routines in Purushotham et al. [2018]
which extracts $26$ features. The event of interest is defined as the
mortality after admission, and the censored time is defined as the last time
of being discharged from the hospital. The definition is similar to that in
Tang et al. [2022]. But since the dataset is not open sourced, according to
our implementation the resulting dataset exhibits a much higher censoring rate
($90.2\%$ as compared to $61.0\%$ as reported in the SODEN paper Tang et al.
[2022]). Since the major purpose of this paper is for the proposal of the NFM
framework, We use our own version of the processed dataset to further verify
the predictive performance of NFM.

Hyperparameter configurations We follow the general training template that
uses MLP as all nonparametric function approximators (i.e., $\widehat{m}$ and
$\widehat{h}$ in the PF scheme, and $\widehat{\nu}$ in the FN scheme), and
train for $100$ epochs across all datasets using Adam as the optimizer. The
tunable parameters and their respective tuning ranges are reported as follows:

Number of layers (network depth)

We tune the network depth $L\in\\{2,3,4\\}$. Typically, the performance of
two-layer MLPs is sufficiently satisfactory.

Number of hidden units in each layer (network width)

We tune the network width $W\in\\{2^{k},5\leq k\leq 10\\}$.

Optional dropout

We optionally apply dropout with probability $p\in\\{0.1,0.2,0.3,0.5,0.7\\}$.

Batch size

We tune batch size within the range $\\{128,256,512\\}$, in the KKBOX dataset,
we also tested with larger batch sizes $\\{1024\\}$.

Learning rate and weight decay

We tune both the learning rate and weight decay coefficient of Adam within
range $\\{0.01,0.001,0.0001\\}$.

Frailty specification

We tested gamma frailty, Box-Cox transformation frailty, and
$\text{IGG}(\alpha)$ frailty with $\alpha\in\\{0,0.25,0.75\\}$. Here note that
$\text{IGG}(0.5)$ is equivalent to gamma frailty. We also empirically tried to
set $\alpha$ to be a learnable parameter and found that this additional
flexibility provides little performance improvement regarding the datasets
used for evaluation.

###  C.4 Implementations

We use pytorch to implement NFM. The source code is provided in the
supplementary material. For the baseline models:

  * •

We use the implementations of CoxPH , GBM, and RSF from the sksurv package
Pölsterl [2020], for the KKBOX dataset, we use the XGBoost library Chen and
Guestrin [2016] to implement GBM and RSF, which might yield some performance
degradation.

  * •

We use the pycox package to implement DeepSurv, CoxTime, and DeepHit models.

  * •

We use the official code provided in the SODEN paper Tang et al. [2022] to
implement SODEN.

  * •

We obtain results of SuMo and DeepEH based on our re-implementations.

##  Appendix D Additional experiments

###  D.1 Recovery assessment of $m(Z)$ in PF scheme

We plot empirical recovery results targeting the $m$ function in (2) in figure
2. The result demonstrates satisfactory recovery with a moderate amount of
data, i.e., $n\geq 1000$.

$0$$50$$100$$0$$1$$2$index$m(Z)$$N=1000$ground truthestimated

$0$$50$$100$$0$$1$$2$index$m(Z)$$N=5000$ground truthestimated

$0$$50$$100$$0$$1$$2$index$m(Z)$$N=10000$ground truthestimated

Figure 2: Visualizations of synthetic data results under the PF scheme of NFM
framework, regarding empirical recovery of the $m$ function in (2)

###  D.2 Recovery assessment of survival functions

To assess the recovery performance of NFM with respect to survival functions,
we consider the following setup: under the same data generation framework as
in section C.2, we compute the test feature $\bar{Z}$ as the sample mean of
all the $100$ hold-out test points. And plot $\widehat{S}(t|\bar{Z})$ against
the ground truth $S(t|\bar{Z})$ regarding both PF and FN schemes. The results
are shown in figure 3. The results suggest that both scheme provides accurate
estimation of survival functions when the sample size is sufficiently large.

$0$$1$$2$$3$$0$$0.25$$0.75$$1$$t$$S(t|Z)$$N=1000$ground truthestimated

$0$$1$$2$$3$$0$$0.25$$0.75$$1$$t$$S(t|Z)$$N=5000$ground truthestimated

$0$$1$$2$$3$$0$$0.25$$0.75$$1$$t$$S(t|Z)$$N=10000$ground truthestimated

$0$$1$$2$$3$$0$$0.25$$0.75$$1$$t$$S(t|Z)$$N=1000$ground truthestimated

$0$$1$$2$$3$$0$$0.25$$0.75$$1$$t$$S(t|Z)$$N=5000$ground truthestimated

$0$$1$$2$$3$$0$$0.25$$0.75$$1$$t$$S(t|Z)$$N=10000$ground truthestimated

Figure 3: Visualizations of synthetic data results under the NFM framework.
The plots in the first row compare the empirical estimates of the survival
function $S(t|\bar{Z})$ against its true value with $\bar{Z}$ being the
average of the features of the $100$ hold-out points, under the PF scheme. The
plots in the second row are obtained using the FN scheme, with analogous
semantics to the first row.

###  D.3 Performance evaluations under the concordance index (C-index)

The concordance index (C-index) Antolini et al. [2005] is yet another
evaluation metric that is commonly used in survival analysis. The C-index
estimates the probability that, for a random pair of individuals, the
predicted survival times of the two individuals have the same ordering as
their true survival times. Formally, C-index is defined as

$\displaystyle\mathcal{C}=\mathbb{P}\left[\widehat{S}(T_{i}\mid
Z_{i})<\widehat{S}(T_{j}\mid Z_{j})\mid T_{i}<T_{j},\delta_{i}=1\right].$ (17)

We report performance evaluations based on C-index over all the $6$ benchmark
datasets in table 4.

Table 4: Survival prediction results measured in C-index (%) on all the $6$ benchmark datasets. In each column, the boldfaced score denotes the best result and the underlined score represents the second-best result. The average rank of each model is reported in the rightmost column. We did not manage to obtain reasonable results for DeepEH and SODEN on two larger datasets MIMIC-III and KKBOX, and we set corresponding ranks to be the worst on those datasets. Model | METABRIC | RotGBSG | FLCHAIN | SUPPORT | MIMIC-III | KKBOX | Ave. Rank
---|---|---|---|---|---|---|---
CoxPH | ${63.42}_{\pm 1.81}$ | ${66.14}_{\pm 1.46}$ | ${79.09}_{\pm 1.11}$ | ${56.89}_{\pm 0.91}$ | ${74.91}_{\pm 0.00}$ | ${83.01}_{\pm 0.00}$ | $11.33$
GBM | ${64.02}_{\pm 1.79}$ | ${67.35}_{\pm 1.16}$ | $\mathbf{79.47}_{\pm 1.08}$ | ${61.46}_{\pm 0.80}$ | ${75.20}_{\pm 0.00}$ | ${85.84}_{\pm 0.00}$ | $7.17$
RSF | ${64.47}_{\pm 1.82}$ | ${67.33}_{\pm 1.34}$ | ${78.75}_{\pm 1.07}$ | ${61.63}_{\pm 0.84}$ | ${75.47}_{\pm 0.17}$ | ${85.79}_{\pm 0.00}$ | $8.00$
DeepSurv | ${63.95}_{\pm 2.12}$ | ${67.20}_{\pm 1.22}$ | ${79.04}_{\pm 1.14}$ | ${60.91}_{\pm 0.85}$ | ${80.08}_{\pm 0.44}$ | ${85.59}_{\pm 0.08}$ | $8.50$
CoxTime | ${66.22}_{\pm 1.69}$ | ${67.41}_{\pm 1.35}$ | ${78.95}_{\pm 1.01}$ | ${61.54}_{\pm 0.87}$ | ${78.78}_{\pm 0.62}$ | $\mathbf{87.31}_{\pm 0.24}$ | $5.00$
DeepHit | ${66.33}_{\pm 1.61}$ | ${66.38}_{\pm 1.07}$ | ${78.48}_{\pm 1.09}$ | $\mathbf{63.20}_{\pm 0.85}$ | ${79.16}_{\pm 0.59}$ | ${86.12}_{\pm 0.26}$ | $6.50$
DeepEH | ${66.59}_{\pm 2.00}$ | $\mathbf{67.93}_{\pm 1.28}$ | ${78.71}_{\pm 1.44}$ | ${61.51}_{\pm 1.04}$ | $-$ | $-$ | $6.33$
SuMo-net | ${64.82}_{\pm 1.80}$ | ${67.20}_{\pm 1.31}$ | ${79.28}_{\pm 1.02}$ | ${62.18}_{\pm 0.78}$ | ${76.23}_{\pm 1.06}$ | ${84.77}_{\pm 0.02}$ | $7.00$
SODEN | ${64.82}_{\pm 1.05}$ | ${66.97}_{\pm 0.75}$ | ${79.00}_{\pm 0.96}$ | ${61.10}_{\pm 0.59}$ | $-$ | $-$ | $10.17$
SurvNode | ${64.64}_{\pm 4.91}$ | ${67.30}_{\pm 1.65}$ | ${76.11}_{\pm 0.98}$ | ${55.37}_{\pm 0.77}$ | $-$ | $-$ | $11.50$
DCM | ${65.76}_{\pm 1.25}$ | ${66.75}_{\pm 1.35}$ | ${78.61}_{\pm 0.79}$ | ${62.19}_{\pm 0.95}$ | ${76.45}_{\pm 0.34}$ | ${83.48}_{\pm 0.07}$ | $8.33$
DeSurv | ${65.88}_{\pm 2.02}$ | ${67.30}_{\pm 1.45}$ | ${78.97}_{\pm 1.64}$ | ${61.47}_{\pm 0.97}$ | $\mathbf{80.97}_{\pm 0.30}$ | ${86.11}_{\pm 0.05}$ | $5.67$
NFM-PF | ${64.98}_{\pm 1.87}$ | $\underline{67.77}_{\pm 1.35}$ | $\underline{79.45}_{\pm 1.03}$ | ${61.33}_{\pm 0.83}$ | ${79.56}_{\pm 0.15}$ | ${86.23}_{\pm 0.01}$ | $\underline{4.67}$
NFM-FN | $\mathbf{66.63}_{\pm 1.82}$ | ${67.73}_{\pm 1.29}$ | ${79.29}_{\pm 0.93}$ | $\underline{62.21}_{\pm 0.41}$ | $\underline{80.18}_{\pm 0.20}$ | $\underline{86.61}_{\pm 0.01}$ | $\mathbf{2.16}$

From table 4, it appears that there’s no clear winner regarding the C-index
metric across the $6$ selected datasets. We conjecture this phenomenon to be
closely related to the loose correlation between the C-index and the
likelihood-based learning objective, as was observed in Rindt et al. [2022].
Therefore we compute the average rank of each model as an overall assessment
of performance, as illustrated in the last column in table 4. The results
suggest that the two NFM models perform better on average.

###  D.4 Benefits of frailty

We compute the (relative) performance gain of NFM-PF and NFM-FN, against their
non-frailty counterparts, namely DeepSurv Katzman et al. [2018] and SuMo-net
Rindt et al. [2022] based on results in table 1, table 2 and table 4. The
results are shown in table 5

Table 5: Relative improvement of NFM models in comparison to their non-frailty counterparts, measured in IBS, INBLL, and C-index. Dataset | NFM-PF vs DeepSurv | NFM-FN vs SuMo-net
---|---|---
| IBS | INBLL | C-index | IBS | INBLL | C-index
METABRIC | $+1.33\%$ | $+1.56\%$ | $+1.61\%$ | $+2.30\%$ | $+3.08\%$ | $+2.79\%$
RotGBSG | $+1.11\%$ | $+0.95\%$ | $+0.84\%$ | $+0.62\%$ | $+0.40\%$ | $+0.79\%$
FLCHAIN | $+1.29\%$ | $+1.32\%$ | $+0.52\%$ | $+0.20\%$ | $+0.27\%$ | $+0.01\%$
SUPPORT | $+0.31\%$ | $+0.23\%$ | $+0.69\%$ | $+2.22\%$ | $+1.76\%$ | $+0.05\%$
MIMIC-III | $+12.38\%$ | $+12.15\%$ | $-0.64\%$ | $+6.18\%$ | $+5.56\%$ | $+5.18\%$
KKBOX | $+2.56\%$ | $+0.51\%$ | $+0.75\%$ | $+8.20\%$ | $+10.38\%$ | $+2.17\%$

The results suggest solid improvement in incorporating frailty, especially for
IBS and INBLL metrics, as the relative increase in performance could be over
$10\%$ for both NFM models. For the IBS and INBLL metrics, the performance
improvement is consistent across all datasets. The only performance
degradation appears on the MIMIC-III dataset evaluated under C-index. This
phenomenon is also understandable: Since the DeepSurv model utilized a variant
of partial likelihood (PL) for model training, as previous works Steck et al.
[2007] pointed out that PL type objective is closely related to the ranking
problem. As C-index could be considered a certain type of ranking measure, it
is possible that DeepSurv obtains better ranking performance than NFM-type
models which are trained using scale-sensitive likelihood objective.