# A presmoothing approach for estimation in semiparametric mixture cure models

Eni Mustalabel=e1 [<EMAIL_ADDRESS>Valentin
<EMAIL_ADDRESS>[ Ingrid Van Keilegomlabel=e3 [
<EMAIL_ADDRESS>ORSTAT, KU Leuven, Belgium,  CREST,
Ensai, France,

###### Abstract

A challenge when dealing with survival analysis data is accounting for a cure
fraction, meaning that some subjects will never experience the event of
interest. Mixture cure models have been frequently used to estimate both the
probability of being cured and the time to event for the susceptible subjects,
by usually assuming a parametric (logistic) form of the incidence. We propose
a new estimation procedure for a parametric cure rate that relies on a
preliminary smooth estimator and is independent of the model assumed for the
latency. We investigate the theoretical properties of the estimators and show
through simulations that, in the logistic/Cox model, presmoothing leads to
more accurate results compared to the maximum likelihood estimator. To
illustrate the practical use, we apply the new estimation procedure to two
studies of melanoma survival data.

cure models,

kernel smoothing,

logistic model,

survival analysis,

###### keywords:

, and

##  1 Introduction

There are many situations in survival analysis problems where some of the
subjects will never experience the event of interest. For instance, as
significant progress is being made for treatment of different types of
cancers, many of the patients get cured of the disease and do not experience
recurrence or cancer-related death. Other examples include study of time to
natural conception, time to default in finance and risk management, time to
early failure of integrated circuits in engineering, time to find a job after
a layoff. However, because of the finite duration of the studies and
censoring, the cured subjects (for which the event never takes place) cannot
be distinguished from the ‘susceptible’ ones. We can just get an indication of
the presence of a cure fraction from the context of the study and a long
plateau (containing many censored observations) with height greater than zero
in the Kaplan-Meier estimator of the survival function. Predicting the
probability of being cured given a set of characteristics is often of
particular interest in order to make better decisions in terms of treatment,
management strategies or public policies. This lead to the development of
mixture cure models.

Mixture cure models were first proposed by [5] and [4]. They assume that the
population is a mixture of two groups: the cured and the susceptible subjects.
Within this very wide class of models, various approaches have been considered
in the literature for modelling and estimating the incidence (probability of
being uncured) and the latency (survival function of the uncured subjects).
Initially, fully parametric models with a logistic regression form of the
incidence and various parametric distributions for the latency were used in
[11, 31, 15]. Later on, more flexible semi-parametric approaches were proposed
for the latency based on the Cox proportional hazards model [24, 21] or
accelerated failure time models [17, 32]. However, they still maintain the
logistic regression model for the incidence. More recently, nonparametric
methods have been developed for both or one of the model components in [30,
20, 2]. In this wide range of models, probably the most commonly used one in
practice is the logistic/Cox mixture cure model [23, 29, 16].

There have been different proposals for estimation in the logistic/Cox mixture
cure model. The presence of a latent variable (the unknown cure status), does
not allow for a ‘direct’ approach as in the classical Cox proportional hazards
model. [15] adapted a marginal likelihood approach computed through Monte
Carlo approximations, whereas [21] and [24] computed the maximum likelihood
estimator via the Expectation-Maximization algorithm. Asymptotic properties of
the latter estimators are investigated in [18], while the procedure is
implemented in the package smcure [7]. One concern about the previous
estimators is that they are obtained by iterative procedures which could be
unstable in practice. In particular, when the sample size is small there are
situations in which the EM algorithm fails to converge (even though the smcure
package can still provide without error the estimates obtained when the
maximum number of iterations is reached). Such problems are for example
reported in [14]. In addition, the maximum likelihood estimator for the
incidence component depends on which variables are included in the latency
model (see for example the illustration in Section 7) and this instability
might in practice lead to unobserved effects (when the effect is not very
strong). In particular, if the latency model is misspecified, even the
estimators of the incidence parameters suffer from induced bias (see for
example [6]).

In this paper, we introduce an alternative estimation method which applies
very broadly and, in particular, for the logistic/Cox mixture cure model. Our
approach focuses on direct estimation of the cure probability without using
distributional assumptions on the latency and iterative algorithms. It relies
on a preliminary nonparametric estimator for the incidence which is then
‘projected’ on a parametric class of functions (like logistic functions). The
idea of constructing a parametric estimator by nonparametric estimation has
been previously proposed for the classical linear regression by [9]. Later on
it was shown to be effective also in the context of variable selection and
functional linear regression [1, 12]. However, its extension to nonlinear
setups has been very little investigated. Here we show that in the context of
mixture cure models, even when a parametric form is assumed for the incidence,
the use of a presmoothed estimator as an intermediate step for obtaining the
parameter estimates often leads to more accurate results. Once the cure
fraction is estimated, we estimate the survival distribution of the uncured
subjects. In the case of the logistic/Cox cure model, this is done by
maximizing the Cox component of the likelihood. In this step, an iterative
algorithm is used to compute the estimators of the baseline cumulative hazard
and the regression parameters. This new approach is of practical relevance
given the popularity of the semiparametric logistic/Cox mixture cure model.
However, the method can be applied more in general to a mixture cure model
with a parametric form of the incidence and other type of models for the
uncured subjects, such as the semiparametric proportional odds model or the
semiparametric AFT model. Our findings suggest that presmoothing has potential
to improve parameter estimation for small and moderate sample size.

The paper is organized as follows. In Sections 2 and 3 we describe the model
and the estimation procedure. Section 4 focuses on the estimation method in
the case of the logistic/Cox mixture cure model. Consistency and asymptotic
normality of the estimators are shown in Section 5. Thanks to the
presmoothing, we are able to present theoretical results under more reasonable
assumptions and thus we contribute to fill a gap between unrealistic technical
conditions and applications. The finite sample performance of the method is
investigated through a simulation study and results are reported in Section 6.
For practical purposes, we propose to make simple and commonly used choices
for the bandwidth and the kernel function in the presmoothing step, and we
show that these choices provide satisfactory results. The proposed estimation
procedure is applied to two medical datasets about studies of patients with
melanoma cancer (see Section 7). We conclude in Section 8 with some discussion
and ideas for further research. Finally, some of the proofs can be found in
Section 9, while the remaining proofs and additional simulation results are
collected in the online Supplementary Material.

##  2 Model description

In the mixture cure model the survival time $T$ can be decomposed as

$T=BT_{0}+(1-B)\infty,$

where $T_{0}$ represents the finite survival time for an uncured individual
and $B$ is an unobserved $0$-$1$ random variable giving the uncured status:
$B=1$ for uncured individuals and $B=0$ otherwise. By convention
$0\cdot\infty=0$. Let $C$ be the censoring time and
$(X^{\prime},Z^{\prime})^{\prime}$ a $(p+q)$-dimensional vector of covariates,
where $x^{\prime}$ denotes the transpose of the vector $x$. Let $\mathcal{X}$
and $\mathcal{Z}$ be the supports of $X$ and $Z$ respectively. Observations
consist of $n$ i.i.d. realizations of $(Y,\Delta,X,Z)$, where $Y=\min(T,C)$ is
the finite follow-up time and $\Delta=\mathds{1}_{\\{T\leq C\\}}$ is the
censoring indicator. Since $Y$ is finite, then necessarily
$\mathbb{P}(C<\infty)=1$, that means the censoring times are finite (which
makes sense given the limited duration of the studies). As a result, censored
survival times of the uncured subjects cannot be distinguished from the cured
ones.

The covariates included in $X$ are those used to model the cure rate, while
the ones in $Z$ affect the survival conditional on the uncured status. This
allows in general to use different variables for modelling the incidence and
the latency but does not exclude situations in which the two vectors $X$ and
$Z$ share some components or are exactly the same. Apart from the standard
assumption in survival analysis that $T_{0}\perp(C,X)|Z$, here we also need

$B\perp(C,T_{0},Z)|X.$ (1)

This implies in particular that

$T\perp C|(X,Z)$ (2)

(see Lemma 1 in the Suplementary Material). Moreover, (1) implies

$\mathbb{P}(T=\infty|X,Z)=\mathbb{P}(T=\infty|X).$ (3)

In addition, in the cure model context we need that the event time $T_{0}$ has
support $[0,\tau_{0}]$, i.e. $\\{T>\tau_{0}\\}=\\{T=\infty\\}$, such that

${\inf_{x}\mathbb{P}(C>\tau_{0}|X=x)>0.}$ (4)

(If the support of $T_{0}$ given $Z=z$ depends on $z$, then we let
$\tau_{0}=\sup\tau_{0}(z)$, where $\tau_{0}(z)$ is the right endpoint of this
support.) This condition tells us that all the observations with $Y>\tau_{0}$
are cured. Even if it might seem restrictive, it is reasonable when a cure
model is justified by a ‘good’ follow-up beyond the time when most of the
events occur and it is commonly accepted in the cure model literature in order
for the mixture cure model to be identifiable and not to overestimate the cure
rate. Since $T_{0}\perp X|Z$, we have

$\mathbb{P}(T_{0}\leq t|X,Z)=\mathbb{P}(T_{0}\leq t|Z),\quad\forall
t\in[0,\tau_{0}].$

We assume a parametric model for the cure rate and we denote by $\pi_{0}(x)$
the cure probability of a subject with covariate $x$, i.e

$\pi_{0}(x)=\mathbb{P}(T=\infty|X=x)=1-\phi(\gamma_{0},x),$

for some parametric model $\\{\phi(\gamma,x):\gamma\in G\\}$ and
$\gamma_{0}\in G$. The first component of $X$ is equal to one and the first
component of $\gamma$ corresponds to the intercept. In order for $\gamma$ to
be identifiable we need the following condition

$\mathbb{P}\left(\phi(\gamma,X)=\phi(\tilde{\gamma},X)\right)=1\qquad\text{
implies that }\qquad\gamma=\tilde{\gamma}.$ (5)

Choosing a parametric model for the incidence seems quite standard in the
literature of mixture cure models ([20, 6, 22]) because of its simplicity and
ease of interpretability (particularly for multiple covariates). To check the
fit of this model in practice, one can compare the prediction error with that
of a more flexible single-index model as done in [2] and for our real data
application in Section 7. It is also possible to test whether this assumption
is reasonable using the test proposed in [19], but this is currently developed
only for one covariate. Among the parametric models for the incidence
component, the most common example is the logistic model, where

$\phi(\gamma,x)=1/(1+{\exp(-\gamma^{\prime}x)}).$ (6)

We state the results in Section 5 for a general parametric model for the
incidence, but then we focus on the logistic function in the simulation study
in Section 6 since it is more of interest in practice. For the uncured
subjects, we can consider a general semiparametric model defined through the
survival function

$S_{u}(t|z)=S_{u}(t|z;\beta,\Lambda)=\mathbb{P}(T_{0}>t|{Z=z},B=1)\quad\text{and}\quad
S_{u}(\tau_{0}|z)=0,$ (7)

where the conditional survival function $S_{u}$ is allowed to depend on a
finite-dimensional parameter, denoted by $\beta\in\mathcal{B}$, and/or an
infinite-dimensional parameter, denoted by $\Lambda\in\mathcal{H}$, with
$\mathcal{B}$ and $\mathcal{H}$ the respective parameter sets. Let
$\beta_{0}\in\mathcal{B}$ and $\Lambda_{0}\in\mathcal{H}$ be the true values
of these parameters. As a result, the conditional survival function
corresponding to $T$ is then

$S(t|x,z)=\mathbb{P}(T>t|{X=x,Z=z})=1-\phi(\gamma_{0},x)+\phi(\gamma_{0},x)S_{u}(t|z).$

The main example we keep in mind is the Cox proportional hazards (PH) model
where $\Lambda_{0}$ is the baseline cumulative hazard. In this case

$S_{u}(t|z)=S_{0}(t)^{\exp(\beta^{\prime}_{0}z)}=\exp(-\Lambda_{0}(t)\exp(\beta^{\prime}_{0}z)),$
(8)

where $S_{0}$ is the baseline survival and $\beta_{0}$ does not contain an
intercept.

##  3 Presmoothing estimation approach

The estimation method we propose is based on a two step procedure. We first
estimate nonparametrically the cure probability for each observation and then
compute an estimator of $\gamma$ as the maximizer of the logistic likelihood,
ignoring the model for the uncured subjects. In the second step, we plug-in
this estimator of $\gamma$ in the full likelihood of the mixture cure model
and fit the latency model using maximum likelihood estimation. In what
follows, we describe in more details these two steps.

Step 1. Even though a parametric model is assumed for the incidence, we start
by computing a nonparametric estimator of the cure probability for each
subject. One possibility is to use the method followed by [20] (see also
[30]), but other estimators are possible as well, as long as the conditions
given in Section 5 are satisfied. The estimator of [20] is defined as follows:

$\hat{\pi}(x)=\prod_{t\in\mathbb{R}}\left(1-\frac{\hat{H}_{1}(dt|x)}{\hat{H}([t,\infty)|x)}\right){,}$
(9)

where
$\hat{H}([t,\infty)|x)=\hat{H}_{1}([t,\infty)|x)+\hat{H}_{0}([t,\infty)|x)$,
$\hat{H}_{1}(dt|x)=\hat{H}_{1}((t-dt,t]|x)$ for small $dt$ and

$\hat{H}_{k}([t,\infty)|x)=\sum_{i=1}^{n}\frac{\tilde{K}_{b}(X_{i}-x)}{\sum_{j=1}^{n}\tilde{K}_{b}(X_{j}-x)}\mathds{1}_{\\{Y_{i}\geq
t,\Delta_{i}=k\\}},\quad k=0,1{,}$

are estimators of

$H_{k}([t,\infty)|x)=\mathbb{P}\left(Y\geq t,\Delta=k|X=x\right),$

$H([t,\infty)|x)=H_{1}([t,\infty)|x)+H_{0}([t,\infty)|x)$. Here
$\tilde{K}_{{b}}$ is a multidimensional kernel function defined in the
following way. If $X$ is composed of continuous and discrete components,
$X=(X_{c},X_{d})\in\mathcal{X}_{c}\times\mathcal{X}_{d}\subset\mathbb{R}^{p_{c}}\times\mathbb{R}^{p_{d}}$
with $p_{c}+p_{d}=p$, then

$\tilde{K}_{b}(X_{i}-x)=K_{b}(X_{c,i}-x_{c})\mathds{1}_{\\{X_{d,i}=x_{d}\\}},$

where $b=b_{n}$ is a bandwidth sequence, $K_{b}(\cdot)=K(\cdot/b)/b^{p_{c}}$
and $K(u)=\prod_{j=1}^{p_{c}}k(u_{j})$, with $k$ a kernel. Note that, one can
compute this estimator with any covariate but here we only use $X$ because of
our assumption (3). The estimator $\hat{\pi}(x)$ coincides with the Beran
estimator of the conditional survival function $S$ at the largest observed
event time $Y_{(m)}$ and does not require any specification of $\tau_{0}$.
Since $\hat{H}_{1}(dt|x)$ is different from zero only at the observed event
times, computation of $\hat{\pi}(x)$ requires only a product over $t$ in the
set of the observed event times. Afterwards, we consider the logistic
likelihood

$\hat{L}_{n,1}(\gamma)=\prod_{i=1}^{n}\phi(\gamma,X_{i})^{1-\hat{\pi}(X_{i})}(1-\phi(\gamma,X_{i}))^{\hat{\pi}(X_{i})}{,}$

and define $\hat{\gamma}_{n}$ as the maximizer of

$\log\hat{L}_{n,1}(\gamma)=\sum_{i=1}^{n}\Big{\\{}\left[1-\hat{\pi}(X_{i})\right]\log\phi(\gamma,X_{i})+\hat{\pi}(X_{i})\log\left[1-\phi(\gamma,X_{i})\right]\Big{\\}}.$
(10)

Existence and uniqueness of $\hat{\gamma}_{n}$ holds under the same conditions
as for the maximum likelihood estimator in the binary outcome regression model
where $1-\hat{\pi}(X_{i})$ is replaced by the outcome $B_{i}$. For example, in
the logistic model, it is required that $p<n$ and the matrix of the variables
$X$ has full rank.

Step 2. Now we consider the likelihood of the mixture cure model. Let

$f_{u}(t|z;\beta,\Lambda)=-(\partial/\partial t)S_{u}(t|z;\beta,\Lambda)$

with $S_{u}(t|z;\beta,\Lambda)$ as defined in (7), denote an element in the
model for the conditional density of $T_{0}$ given $Z=z$, which is supposed to
exist and belong to the model. Assuming non informative censoring and that the
distribution of the covariates does not carry information on the parameters
$\beta$, $\Lambda$, the likelihood criterion is then

$L_{n,2}(\beta,\Lambda,\gamma)=\prod_{i=1}^{n}\left\\{\phi(\gamma,X_{i})f_{u}(Y_{i}|Z_{i};\beta,\Lambda)\right\\}^{\Delta_{i}}\\!\left\\{1\\!-\phi(\gamma,X_{i})+\phi(\gamma,X_{i})S_{u}(Y_{i}|Z_{i};\beta,\Lambda)\right\\}^{1-\Delta_{i}},$
(11)

and we maximize it w.r.t. $\beta$ and $\Lambda$ for $\gamma=\hat{\gamma}_{n}$,
i.e. $(\hat{\beta}_{n},\hat{\Lambda}_{n})$ are the maximizers of

$\hat{l}_{n}(\beta,\Lambda,\hat{\gamma}_{n})=\frac{1}{n}\sum_{i=1}^{n}\ell(Y_{i},\Delta_{i},X_{i},Z_{i};\beta,\Lambda,\hat{\gamma}_{n}),$
(12)

over a set of possible values for $\beta$ and $\Lambda$, where

$\begin{split}\ell(Y_{i},\Delta_{i},X_{i},Z_{i};\beta,\Lambda,\gamma)&=\Delta_{i}\log
f_{u}(Y_{i}|Z_{i};\beta,\Lambda)\\\
&\quad+(1-\Delta_{i})\log\left\\{1-\phi(\gamma,X_{i})+\phi(\gamma,X_{i})S_{u}(Y_{i}|Z_{i};\beta,\Lambda)\right\\}.\end{split}$
(13)

##  4 Presmoothing estimation for the parametric/Cox mixture cure model

In the sequel we focus on the case of a Cox PH model defined in (8) for the
conditional law of $T_{0}$. The criterion defined in (12) becomes

$\hat{l}_{n}(\beta,\Lambda,\hat{\gamma}_{n})=\frac{1}{n}\sum_{i=1}^{n}\Delta_{i}\left\\{{\mathds{1}_{\\{Y_{i}<\tau_{0}\\}}}[\log\Delta\Lambda(Y_{i})+\beta^{\prime}Z_{i}]-\Lambda(Y_{i})e^{\beta^{\prime}Z_{i}}\right\\}\\\
+\frac{1}{n}\sum_{i=1}^{n}(1-\Delta_{i})\log\left\\{1-\phi(\hat{\gamma}_{n},X_{i})+\phi(\hat{\gamma}_{n},X_{i})\exp\left(-\Lambda(Y_{i})e^{\beta^{\prime}Z_{i}}\right)\right\\},\,$
(14)

and has to be maximized with respect to $\beta$ and $\Lambda$ in the class of
step functions $\Lambda$ defined on $[0,\tau_{0}]$ (thus by definition
$\Lambda(t)={\infty}$ if $t>\tau_{0}$), with jumps of size $\Delta\Lambda$ at
the event times. The indicator of the event $\\{Y_{i}<\tau_{0}\\}$ in the
first term is needed in case the distribution of the event times has a jump at
$\tau_{0}$ meaning that $\mathbb{P}(T_{0}=\tau_{0}|Z)>0$. In such a case
$f_{u}(\tau_{0}|Z;\beta,\Lambda)=\exp(-\Lambda(\tau_{0})e^{\beta^{\prime}Z})$
where $\Lambda(\tau_{0})=\lim_{t\uparrow\tau_{0}}\Lambda(t)$. Otherwise, if
$\mathbb{P}(T_{0}=\tau_{0}|Z)=0$, then for all uncensored observations we have
$\mathds{1}_{\\{Y<\tau_{0}\\}}=1$ with probability one. Thus, the presence of
the indicator function can be neglected. As in [18], it can be shown that

$(\hat{\beta}_{n},\hat{\Lambda}_{n})=\arg\max_{\beta,\Lambda}\hat{l}_{n}(\beta,\Lambda,\hat{\gamma}_{n})$
(15)

exists and it is finite. Moreover, for any given $\beta$ and $\gamma$, the
${\Lambda}_{n,\beta,\gamma}$ which maximizes
$\hat{l}_{n}(\beta,\Lambda,\gamma)$ in (14), with respect to $\Lambda$ with
jumps at the event times, can be characterized as

${\Lambda}_{n,\beta,\gamma}(t)=\frac{1}{n}\sum_{i=1}^{n}\frac{\Delta_{i}\mathds{1}_{\\{Y_{i}\leq{t{,Y_{i}<\tau_{0}}}\\}}}{\frac{1}{n}\sum_{j=1}^{n}\mathds{1}_{\\{{{Y_{i}\leq
Y_{j}\leq\tau_{0}}}\\}}\exp(\beta^{\prime}Z_{j})\left\\{{\Delta_{j}}+(1-{\Delta_{j}})g_{j}(Y_{j},{\Lambda}_{n,\beta},\beta,\gamma)\right\\}},$
(16)

where

$g_{j}(t,\Lambda,\beta,\gamma)=\frac{\phi(\gamma,X_{j})\exp\left(-\Lambda(t)\exp\left(\beta^{\prime}Z_{j}\right)\right)}{1-\phi(\gamma,X_{j})+\phi(\gamma,X_{j})\exp\left(-\Lambda(t)\exp\left(\beta^{\prime}Z_{j}\right)\right)}.$
(17)

Next, we could define

$\hat{\beta}_{n}=\arg\max_{\beta}\hat{l}_{n}(\beta,{\Lambda}_{n,\beta,\hat{\gamma}_{n}},\hat{\gamma}_{n})\quad\text{and}\quad\hat{\Lambda}_{n}={\Lambda}_{n,\hat{\beta}_{n},\hat{\gamma}_{n}}.$

To compute $(\hat{\beta}_{n},\hat{\Lambda}_{n})$ we use an iterative algorithm
based on profiling. To be precise, we start with initial values which are the
maximum partial likelihood estimator and the Breslow estimator (as if there
was no cure fraction) and we iterate between the next two steps until
convergence:

  * a)

Compute the weights

$w_{j}^{(m)}={\Delta_{j}}+(1-{\Delta_{j}})\frac{\phi(\hat{\gamma}_{n},X_{j})\hat{S}^{(m)}_{u}(Y_{j}|Z_{j})}{1-\phi(\hat{\gamma}_{n},X_{j})+\phi(\hat{\gamma}_{n},X_{j})\hat{S}^{(m)}_{u}(Y_{j}|Z_{j})},$

where

$\hat{S}^{(m)}_{u}(Y_{j}|Z_{j})=\exp\left(-\hat{\Lambda}^{(m)}_{n}(Y_{j})\exp\left(\hat{\beta}^{(m)^{\prime}}_{n}Z_{j}\right)\right),$

using the estimators $\hat{\Lambda}^{(m)}_{n}$, $\hat{\beta}^{(m)}_{n}$ of the
previous step.

  * b)

Using the previous weights, update the estimators for $\Lambda$ and $\beta$,
i.e. $\hat{\beta}^{(m+1)}_{n}$ is the maximizer of

$\prod_{i=1}^{n}\left\\{\frac{e^{\beta^{\prime}Z_{i}}}{\sum_{{Y_{k}\geq
Y_{i}}}w_{k}^{(m)}e^{\beta^{\prime}Z_{k}}}\right\\}^{\Delta_{i}}$

and

$\hat{\Lambda}_{n}^{(m+1)}(t)=\sum_{i=1}^{n}\frac{\Delta_{i}\mathds{1}_{\\{{Y_{i}\leq
t{,Y_{i}<\tau_{0}}}\\}}}{\sum_{j=1}^{n}\mathds{1}_{\\{{{Y_{i}\leq
Y_{j}\leq\tau_{0}}}\\}}w_{j}^{(m)}\exp\left(\hat{\beta}_{n}^{(m+1)^{\prime}}Z_{j}\right)}.$
(18)

The update of $\Lambda$ an $\beta$ in Step (b) coincides with the maximization
step of the EM algorithm and the weights $w^{(m)}$ correspond to the
expectation of the latent variable $B$ given the observed data and the current
parameter values. However, unlike the maximum likelihood estimation [24], we
are keeping $\hat{\gamma}_{n}$ fixed while performing this iterative
algorithm. The estimator $\hat{\Lambda}_{n}$ seems to depend on the unknown
$\tau_{0}$. However, with data at hand, one could easily proceed without
knowing $\tau_{0}$. Indeed, if there are ties at the last uncensored
observation, then $\tau_{0}$ is revealed by the data. On the other hand, if
there are no ties, all uncensored observations will be smaller than
$\tau_{0}$, hence no need to know $\tau_{0}$.

As suggested in [25, 24], we impose the zero-tail constraint, meaning that
$\hat{S}_{u}^{(m)}$ is forced to be equal to zero beyond the last event. In
this way, all censored observations in the plateau are assigned to the cured
group.

##  5 Asymptotic results

We first explain why presmoothing allows for more realistic asymptotic results
in semiparametric mixture cure models. Next, we show consistency and
asymptotic normality of the proposed estimators $\hat{\gamma}_{n}$,
$\hat{\beta}_{n}$ and $\hat{\Lambda}_{n}$ for the parametric/Cox mixture cure
model when, in Step 1, we use a general nonparametric estimator $\hat{\pi}$ of
$\pi_{0}$ that satisfies certain assumptions. Afterwards, we verify these
conditions for the particular estimator $\hat{\pi}$ in (9). Some of the proofs
can be found in Section 9 and the rest in the online Supplementary Material.
The assumptions mentioned in Section 2 are assumed to be satisfied throughout
this section. In addition $Var(Z)$ is supposed to have full rank.

###  5.1 A challenge with mixture cure models

To derive asymptotic results, in most of the existing literature it has been
assumed that

$\inf_{z}\mathbb{P}(T_{0}\geq\tau_{0}|Z=z)>0,$ (19)

[18, 20]. In nonparametric approaches such a condition keeps the denominators
away from zero. In the parametric/Cox mixture cure model, it guarantees that
the baseline distribution stays bounded on the compact support $[0,\tau_{0}]$.
However, condition (19) implies that
$\inf_{z}\mathbb{P}(Y=\tau_{0}{,\Delta=1}|Z=z)>0$, a condition which is not
frequently satisfied in real-data applications.

One could imagine that, instead of imposing condition (19), it could be
possible to proceed as follows: first restrict to events on $[0,\tau^{*}]$ for
some $\tau^{*}<\tau_{0}$ such that

$\inf_{z}\mathbb{P}(Y\geq\tau^{*}{,\Delta=1}|Z=z)>0,$ (20)

next derive the asymptotics, and finally let $\tau^{*}$ tend to $\tau_{0}$.
This idea is used, for instance, in Cox PH model, see [13] chapter 8, or [3].
However, this idea does not seem to work for mixture cure models without
suitable adaptation. This is because it implicitly requires that $\beta_{0}$
and $\Lambda_{0}$ are identifiable from the restricted data. Here,
identifiability means that the true values $\beta_{0}$ and $\Lambda_{0}$ of
the parameters maximize the expectation of the criterion maximized to obtain
the estimators. Two aspects have to be taken into account when analyzing this
identifiability. The first aspect is related to the parameter identifiability
in the semiparametric model for $T_{0}$ when the events are restricted to
$[0,\tau^{*}]$. This property is satisfied in the common models, in particular
it holds true in the Cox PH model as soon as $Var(Z)$ has full rank. The
second aspect is the additional complexity induced by the mixture with a cure
fraction. If the cure fraction is unknown and one decides to restrict to
events on $[0,\tau^{*}]$, the parameter identifiability is likely lost because
the events $\\{T_{0}\in(\tau^{*},\tau_{0}]\\}$ and $\\{T=\infty\\}$ are not
distinguishable. The usual remedy for this is to impose (19), so that
$\tau^{*}$ could be taken equal to $\tau_{0}$.

Presmoothing allows to avoid condition (19) and thus to fill the gap between
the technical conditions and the reality of the data. This is possible
because, when using the presmoothing, the conditional probability of the event
$\\{T=\infty\\}$ is identified by other means. We are thus able to prove the
consistency of $\hat{\beta}$ and $\hat{\Lambda}$ without imposing (19).
Deriving the asymptotic normality without (19) remains an open problem which
will be addressed elsewhere.

###  5.2 Consistency

We first prove consistency of $\hat{\gamma}_{n}$ and then use that result to
obtain consistency of $\hat{\Lambda}_{n}$ and $\hat{\beta}_{n}$. In order to
proceed with our results, the following conditions will be used.

  * (AC1)

$\sup_{x\in\mathcal{X}}\left|\hat{\pi}(x)-\pi_{0}(x)\right|\to 0$ almost
surely.

  * (AC2)

The parameters $\beta_{0}$ and $\gamma_{0}$ lie in the interior of compact
sets $B\subset\mathbb{R}^{q}$, $G\subset\mathbb{R}^{p}$.

  * (AC3)

There exist some constants $a>0$, $c>0$ such that

$\left|\phi(\gamma_{1},x)-\phi(\gamma_{2},x)\right|\leq
c\|\gamma_{1}-\gamma_{2}\|^{a},\qquad\forall\gamma_{1},\gamma_{2}\in
G,\,\forall x\in\mathcal{X},$

where $\|\cdot\|$ denotes the Euclidean distance.

  * (AC4)

$\inf_{\gamma\in G}\inf_{x\in\mathcal{X}}\phi(\gamma,x)>0$ and
$\inf_{\gamma\in G}\inf_{x\in\mathcal{X}}\phi(\gamma,x)<1$.

  * (AC5)

The covariates are bounded: $\mathbb{P}\left(\|Z\|<{m}\text{ and
}\|X\|<m\right)=1$ for some $m>0$.

  * (AC6)

The baseline hazard function $\lambda_{0}(t)=\Lambda^{\prime}_{0}(t)$ is
strictly positive and continuous on $[0,{\tau_{0})}$.

  * (AC7)

With probability one, the conditional distribution function of the censoring
times $F_{C}(t|x,z)$ is continuous in $t$ on $[0,{\tau_{0}}]$ and there exists
a constant $C>0$ such that

$\inf_{0\leq
t_{1}<t_{2}\leq\tau_{0}}\inf_{x,z}\frac{F_{C}(t_{2}|x,z)-F_{C}(t_{1}|x,z)}{t_{2}-t_{1}}>C.$

(AC1) is a minimal assumption given that we want to match $\phi(\gamma,\cdot)$
to $\hat{\pi}(\cdot)$. (AC2) to (AC4) are mild conditions satisfied by usual
binary regression models, like for instance the logistic one, and (AC5) is
always satisfied in practice for large $m$.

######  Theorem 1.

Let the estimator $\hat{\gamma}_{n}$ be defined as in (10). Assume that
(AC1)-(AC4) hold. Then, $\hat{\gamma}_{n}\to\gamma_{0}$ almost surely.

######  Theorem 2.

Let the estimators $\hat{\beta}_{n}$ and $\hat{\Lambda}_{n}$ be defined as in
Section 4. Assume that (AC1)-(AC7) hold. Then, with probability one,
$\|\hat{\beta}_{n}-\beta_{0}\|\to 0$, where $\|\cdot\|$ denotes the Euclidean
distance. Moreover, for any $\tau^{*}\leq\tau_{0}$ satisfying (20), with
probability one,

$\sup_{t\in[0,{\tau^{*}}]}\left|\hat{\Lambda}_{n}(t)-\Lambda_{0}(t)\right|\to
0.$

When condition (19) is satisfied and $\tau^{*}=\tau_{0}$ in the previous
Theorem, we are referring to the continuous version of $\Lambda_{0}$, i.e.
$\Lambda_{0}(\tau_{0})=\lim_{t\uparrow\tau_{0}}\Lambda_{0}(t)$. Note that, by
definition, we also have
$\hat{\Lambda}_{n}(\tau_{0})=\lim_{t\uparrow\tau_{0}}\hat{\Lambda}_{n}(t)$.

###  5.3 Asymptotic normality

We first derive asymptotic normality of $\hat{\gamma}_{n}$ following the
approach in [8]. Theorem 2 in that paper provides sufficient conditions for
the $\sqrt{n}$ normality of parametric estimators obtained by minimizing an
objective function that depends on a preliminary infinite dimensional
estimator $\hat{\pi}$. In our case, since $\hat{\gamma}_{n}$ solves

$\frac{1}{n}\nabla_{\gamma}\log\hat{L}_{n,1}(\gamma)=0,$

where $\nabla_{\gamma}$ denotes the vector-valued partial differentiation
operator with respect to the components of $\gamma$, it follows that
$\hat{\gamma}_{n}$ minimizes the function

$\left\|\frac{1}{n}\nabla_{\gamma}\log\hat{L}_{n,1}(\gamma)\right\|=\left\|\frac{1}{n}\sum_{i=1}^{n}m(X_{i};\gamma,\hat{\pi})\right\|,$

where

$m(x;\gamma,\pi)=\left[\frac{1-\pi(x)}{\phi(\gamma,x)}-\frac{\pi(x)}{1-\phi(\gamma,x)}\right]\nabla_{\gamma}\phi(\gamma,x).$
(21)

Hence, we only need to check that the conditions of Theorem 2 in [8] are
satisfied. To do that, we need the following assumptions which are stronger
than the previous (AC1)-(AC4).

  * (AN1)

The parameter $\gamma_{0}$ lies in the interior of a compact set
$G\subset\mathbb{R}^{p}$ and, for each $x\in\mathcal{X}$, the function
$\gamma\mapsto\phi(\gamma,x)$ is twice continuously differentiable with
uniformly bounded derivatives in $G\times\mathcal{X}$ and satisfies (AC4).

  * (AN2)

$\pi_{0}(\cdot)$ belongs to a class of functions $\Pi$ such that

$\int_{0}^{\infty}\sqrt{\log
N(\epsilon,\Pi,\|\cdot\|_{\infty})}\,\mathrm{d}\epsilon<\infty,$

where $N(\epsilon,\Pi,\|\cdot\|_{\infty})$ denotes the $\epsilon$-covering
number of the space $\Pi$ with respect to
$\|\pi\|_{{\infty}}=\sup_{x\in\mathcal{X}}|\pi(x)|$.

  * (AN3)

The matrix
$\mathbb{E}\left[\nabla_{\gamma}\phi(\gamma_{0},X)\nabla_{\gamma}\phi(\gamma_{0},X)^{\prime}\right]$
is positive definite.

  * (AN4)

The estimator $\hat{\pi}(\cdot)$ satisfies the following properties:

    * (i)

$\mathbb{P}\left(\hat{\pi}(\cdot)\in\Pi\right)\to 1$.

    * (ii)

$\left\|\hat{\pi}(x)-\pi_{0}(x)\right\|_{{\infty}}=o_{P}(n^{-1/4})$.

    * (iii)

There exists a function $\Psi$ such that

$\begin{split}&\mathbb{E}^{*}\left[\left(\hat{\pi}(X)-\pi_{0}(X)\right)\left(\frac{1}{\phi(\gamma_{0},X)}+\frac{1}{1-\phi(\gamma_{0},X)}\right)\nabla_{\gamma}\phi(\gamma_{0},X)\right]\qquad\qquad\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad=\frac{1}{n}\sum_{i=1}^{n}\Psi(Y_{i},\Delta_{i},X_{i})+R_{n},\end{split}$

where $\mathbb{E}^{*}$ denotes the conditional expectation given the sample,
taken with respect to the generic variable $X$. Moreover,
$\mathbb{E}[\Psi(Y,\Delta,X)]=0$ and $\|R_{n}\|=o_{P}(n^{-1/2})$.

######  Theorem 3.

Let the estimator $\hat{\gamma}_{n}$ be defined as in (10). Assume that
(AN1)-(AN4) hold. Then,

$n^{1/2}\left(\hat{\gamma}_{n}-\gamma_{0}\right)\xrightarrow{d}N(0,\Sigma_{\gamma})$

with covariance matrix $\Sigma_{\gamma}$ defined in (A28).

For deriving the asymptotic distribution of $\hat{\beta}_{n}$ and
$\hat{\Lambda}_{n}$ we assume, for simplicity, that condition (19) is
satisfied. In such case, in Theorem 2 we can take $\tau^{*}=\tau_{0}$ and
obtain uniform strong consistency of $\hat{\Lambda}_{n}$ on the whole support
$[0,\tau_{0}]$. We believe that, at the price of additional technicalities,
asymptotic distributional theory can be obtained also without imposing (19),
as we did for the consistency in Theorem 2. This conjecture is supported by
simulations but we leave the problem to be addressed by future research.

######  Theorem 4.

Let the estimators $\hat{\beta}_{n}$ and $\hat{\Lambda}_{n}$ be defined as in
Section 4. Assume that condition (19), (AN1)-(AN4) and (AC2), (AC5)-(AC7)
hold. Then,

$\left\langle\sqrt{n}\left(\hat{\Lambda}_{n}-\Lambda_{0}\right),\sqrt{n}\left(\hat{\beta}_{n}-\beta_{0}\right)\right\rangle\to
G$

weakly in $l^{\infty}(\mathcal{H}_{\mathfrak{m}})$, where
$\mathcal{H}_{\mathfrak{m}}$ is a functional space defined in Section 9.1,
$l^{\infty}(\mathcal{H}_{\mathfrak{m}})$ denotes the space of bounded real-
valued functions on $\mathcal{H}_{\mathfrak{m}}$, $G$ is a tight Gaussian
process in $l^{\infty}(\mathcal{H}_{\mathfrak{m}})$ with mean zero and
covariance process given in (A39) and for
$h=(h_{1},h_{2})\in\mathcal{H}_{\mathfrak{m}}$

$\langle\Lambda,\beta\rangle(h)=\int_{0}^{\tau_{0}}h_{1}(t)\,\mathrm{d}\Lambda(t)+h_{2}^{\prime}\beta.$

The asymptotic variances of each component of $\hat{\beta}_{n}$ and of
$\hat{\Lambda}_{n}(t)$ can be obtained from the covariance process in (A39) by
taking $h_{1}(t)=0$ for all $t$ and $h_{2}=e_{i}$ (the $i$th unit vector) or
$h_{2}=0$ and $h_{1}(s)=\mathds{1}_{\\{s\leq t\\}}$. We leave the details
about these covariance matrices in the Supplementary Material because they
have quite complicated expressions that require definitions of several other
quantities. Even though it could be possible in principle to estimate the
asymptotic standard errors through plug-in estimators and numerical inverse,
we think that this is not feasible in practice and we do not intend to exploit
it further. Instead, we use a bootstrap procedure for estimation of the
standard errors in the application discussed in Section 7. However, the
maximum likelihood estimators are not more favorable in this regard. For
example, in the logistic/Cox model, the proposed estimators of the asymptotic
variance in [18] also involve solving numerically complicated nonlinear
equalitons. For this reason, bootstrap is used in practice to estimate the
standard errors even for the maximum likelihood estimators.

By considering a two-step procedure, where estimation of the incidence
parameters is performed independently of the latency model, we expect to loose
efficiency of the estimators. However, this does not cause major concern
because our purpose is to provide an alternative estimation method that
performs better than the maximum likelihood estimation with sample sizes
usually encountered in practice. Efficiency is a key concept for the
asymptotics of the estimators, and in general there is no particular need for
another method since the MLE would be the best choice. However, in many
nonlinear models, like the mixture cure models, the asymptotic approximation
is poor and the efficiency becomes a less relevant purpose for real data
sample sizes. Hence, we choose to trade efficiency for better performance in a
wider range of applications.

###  5.4 Verification of assumptions for $\hat{\pi}$

Next we show that our assumptions (AN1)-(AN4) of the asymptotic theory are
satisfied for the nonparametric estimator $\hat{\pi}$ defined in (9) and the
logistic model in (6). For reasons of simplicity, since we use results
available in the literature only for a one-dimensional covariate, we consider
only cases with one continuous covariate. In order for assumption (AN4) to be
satisfied we need the following conditions:

  * (C1)

The bandwidth $b$ is such that $nb^{4}\to 0$ and $nb^{3+\xi}/(\log
b^{-1})\to\infty$ for some $\xi>0$.

  * (C2)

The support $\mathcal{X}$ of $X$ is a compact subset of $\mathbb{R}$. The
density $f_{X}(\cdot)$ of $X$ is bounded away from zero and twice
differentiable with bounded second derivative.

  * (C3)

The kernel $k$ is a twice continuously differentiable, symmetric probability
density function with compact support and $\int uk(u)\,\mathrm{d}u=0$.

  * (C4)

(i) The functions $H([0,t]|x)$, $H_{1}([0,t]|x)$ are twice differentiable with
respect to $x$, with uniformly bounded derivatives for all $t\leq\tau_{0}$,
$x\in\mathcal{X}$. Moreover, there exist continuous nondecreasing functions
$L_{1}$, $L_{2}$, $L_{3}$ such that $L_{i}(0)=0$, $L_{i}(\tau_{0})<\infty$ and
for all $t,s\in[0,\tau_{0}]$, $x\in\mathcal{X}$,

$\begin{split}\left|H_{c}(t|x)-H_{c}(s|x)\right|\leq\left|L_{1}(t)-L_{1}(s)\right|,&\quad\left|H_{1c}(t|x)-H_{1c}(s|x)\right|\leq\left|L_{1}(t)-L_{1}(s)\right|\\\
\left|\frac{\partial H_{c}(t|x)}{\partial x}-\frac{\partial
H_{c}(s|x)}{\partial x}\right|&\leq\left|L_{2}(t)-L_{2}(s)\right|\\\
\left|\frac{\partial H_{1c}(t|x)}{\partial x}-\frac{\partial
H_{1c}(s|x)}{\partial
x}\right|&\leq\left|L_{3}(t)-L_{3}(s)\right|,\end{split}$

where the subscript c denotes the continuous part of a function.

(ii) The jump points for the distribution function $G(t|x)$ of the censoring
times given the covariate, are finite and the same for all $x$. The partial
derivative of $G(t|x)$ with respect to $x$ exists and is uniformly bounded for
all $t\leq\tau_{0}$, $x\in\mathcal{X}$. Moreover, the partial derivative with
respect to $x$ of $F(t|x)$ (distribution function of the survival times $T$
given $X=x$) exists and is uniformly bounded for all $t\leq\tau_{0}$,
$x\in\mathcal{X}$.

  * (C5)

The survival time $T$ and the censoring time $C$ are independent given $X$.

(C1) to (C5) are conditions guaranteeing the rates of convergence and the
i.i.d. representation [10]. In case of discrete covariates we also need to
have only a finite number of atoms. Assumption (C5) is needed because we are
dealing with the distribution of $T$ conditional only on the covariate $X$
(since the cure rate depends only on $X$).

######  Theorem 5.

Under the conditions (C1)-(C5), the assumptions (AN1)-(AN4) hold true for the
logistic model and the estimator $\hat{\pi}(x)$ defined in (9).

##  6 Simulation study

In this section we focus on the logistic/Cox mixture cure model and evaluate
the finite sample performance of the proposed method. Comparison is made with
the maximum likelihood estimator implemented in the package smcure.

We first illustrate through a brief example the convergence problems of the
smcure estimator. We consider a model where the incidence depends on four
independent covariates: $X_{1}\sim N(0,2)$, $X_{2}\sim\text{Uniform}(-1,1)$,
$X_{3}\sim\text{Bernoulli}(0.8)$, $X_{4}\sim\text{Bernoulli}(0.2)$. The
latency depends on $Z_{1}=X_{1}$, $Z_{2}=X_{3}$ and $Z_{3}=X_{4}$. We generate
the cure status $B$ as a Bernoulli random variable with success probability
$\phi(\gamma,X)$ where $\phi$ is the logistic function and
$\gamma=(0.6,-1,1,2.5,1.2)$. The survival times for the uncured observations
are generated according to a Weibull proportional hazards model

$S_{u}(t|z)=\exp\left(-\mu t^{\rho}\exp(\beta^{\prime}z)\right),$

and are truncated at $\tau_{0}=14$ for $\rho=1.75$, $\mu=1.5$,
$\beta=(-0.8,0.9,0.5)$. The censoring times are independent from $X$ and $T$.
They are generated from the exponential distribution with parameter
$\lambda_{C}=0.22$ and are truncated at $\tau=16$. We generate $1000$ datasets
according to this model with sample size $n=100$, and we observe that smcure
fails to converge in $43\%$ of the cases. Convergence fails mainly in the
$\gamma$ parameter, with only $17\%$ of the cases failing to converge also for
the $\beta$ parameter (because of the unreasonable $\gamma$ estimators). On
the other hand, there was no convergence problem in the second step of the
presmoothing approach. In addition, even among the cases where smcure
converged, the presmoothing approach showed significantly better behavior, as
can be seen in Table 1.

Table 1: Bias, variance and MSE of $\hat{\gamma}$ and $\hat{\beta}$ for smcure
and our approach among the iterations that converged for smcure.

presmoothing smcure Par. Bias Var. MSE Bias Var. MSE $\gamma_{1}$ $-0.113$
$0.620$ $0.633$ $0.200$ $8.318$ $8.358$ $\gamma_{2}$ $-0.073$ $0.156$ $0.162$
$-0.388$ $3.085$ $3.236$ $\gamma_{3}$ $-0.071$ $0.546$ $0.551$ $0.280$ $1.957$
$2.035$ $\gamma_{4}$ $0.037$ $1.326$ $1.327$ $0.704$ $14.395$ $14.891$
$\gamma_{5}$ $-0.250$ $8.398$ $8.461$ $1.621$ $36.450$ $36.945$ $\beta_{1}$
$-0.014$ $0.011$ $0.012$ $-0.017$ $0.012$ $0.012$ $\beta_{2}$ $0.024$ $0.064$
$0.065$ $0.026$ $0.065$ $0.065$ $\beta_{3}$ $-0.053$ $0.165$ $0.168$ $-0.053$
$0.166$ $0.169$

Hence, in the cases in which smcure exhibits very poor behavior, the
presmoothing is obviously superior. Next, we focus on models for which smcure
behaves reasonable (there are convergence problems in less then $3\%$ of the
cases) and show that, even in such scenarios presmoothing can lead to more
accurate results.

We consider four different models and for each of them various choices of the
parameters in order to cover a wide range of scenarios. The models are as
follows.

Model 1. Both incidence and latency depend on one covariate $X$, which is
uniform on $(-1,1)$. We generate the cure status $B$ as a Bernoulli random
variable with success probability $\phi(\gamma,X)$ where $\phi$ is the
logistic function. The survival times for the uncured observations are
generated according to a Weibull proportional hazards model

$S_{u}(t|x)=\exp\left(-\mu t^{\rho}\exp(\beta x)\right),$

and are truncated at $\tau_{0}$ for $\rho=1.75$, $\mu=1.5$, $\beta=1$ and
$\tau_{0}=4$. The censoring times are independent from $X$ and $T$. They are
generated from the exponential distribution with parameter $\lambda_{C}$ and
are truncated at $\tau=6$.

Model 2. Both incidence and latency depend on one covariate $X$ with standard
normal distribution. The cure status and the survival times for the uncured
observations are generated as in Model 1 for $\rho=1.75$, $\mu=1.5$, $\beta=1$
and $\tau_{0}=10$. The censoring times are generated according to a Weibull
proportional hazards model

$S_{C}(t|x)=\exp\left(-\nu\mu t^{\rho}\exp(\beta_{C}x)\right),$

for $\beta_{C}=1$ and various choices of $\nu$ and are truncated at $\tau=15$.

Model 3. For the incidence we consider three independent covariates: $X_{1}$
is normal with mean zero and standard deviation $2$, $X_{2}$ and $X_{3}$ are
Bernoulli random variables with parameters $0.6$ and $0.4$ respectively. The
latency also depends on three covariates: $Z_{1}=X_{1}$, $Z_{2}$ is a uniform
random variable on $(-3,3)$ independent of the previous ones and
$Z_{3}=X_{2}$. The cure status and the survival times for the uncured
observations are generated as in Model 1 for $\rho=1.75$, $\mu=1.5$ and
different choices of the other parameters. The censoring times are generated
independently of the previous variables from an exponential distribution with
parameter $\lambda_{C}$ and are truncated at $\tau$, for given choices of
$\lambda_{C}$ and $\tau$.

Model 4. This setting is obtained by adding an additional continuous covariate
to the incidence component of Model 3. To be precise, $X_{1}$ is normal with
mean zero and standard deviation $2$, $X_{2}$ is uniform on $(-1,1)$
independent of the other variables, $X_{3}$ and $X_{4}$ are Bernoulli random
variables with parameters $0.6$ and $0.4$ respectively. As in Model 3,
$Z_{1}=X_{1}$, $Z_{2}$ is a uniform random variable on $(-3,3)$ independent of
the previous ones and $Z_{3}=X_{3}$. The event and censoring times are
generated as in the previous model.

For the four models we choose the values of the unspecified parameters in such
a way that the cure rate is around $20\%$, $30\%$, $50\%$ (corresponding
respectively to scenarios $1$, $2$ and $3$) and the censoring rate corresponds
to three levels (with a difference of $5\%$ between each other). The
specification of the parameters and the corresponding censoring rate and
percentage of the observations in the plateau are given Table 2. Note that,
within each scenario, the fraction of the observations in the plateau
decreases as the censoring rate increases because more cured observations are
censored earlier and as a result are not observed in the plateau. This makes
the estimation of the cure rate more difficult. The truncation of the survival
and censoring times on $[0,\tau_{0}]$ and $[0,\tau]$ is made in such a way
that $\tau_{0}<\tau$ and condition (19) is satisfied but in practice it is
unlikely to observe event times at $\tau_{0}$. In this way, we try to find a
compromise between theoretical assumptions and real-life scenarios.

Table 2:  Parameter values and model characteristics for each scenario.

Model Parameters Scenario Cens. Cens. Cens. Plateau level parameters rate $1$
$\lambda_{C}=0.1$ $25\%$ $15\%$ $\gamma=(1.75,2)$ $1$ $2$ $\lambda_{C}=0.2$
$30\%$ $11\%$ $3$ $\lambda_{C}=0.3$ $35\%$ $9\%$ $1$ $\lambda_{C}=0.1$ $34\%$
$22\%$ 1 $\gamma=(1,1.5)$ $2$ $2$ $\lambda_{C}=0.25$ $40\%$ $15\%$ $3$
$\lambda_{C}=0.4$ $46\%$ $10\%$ $1$ $\lambda_{C}=0.2$ $54\%$ $32\%$
$\gamma=(0.1,5)$ $3$ $2$ $\lambda_{C}=0.4$ $59\%$ $23\%$ $3$ $\lambda_{C}=0.7$
$65\%$ $15\%$ $1$ $\nu=1/15$ $25\%$ $7\%$ $\gamma=(1.5,0.5)$ $1$ $2$ $\nu=1/7$
$30\%$ $4\%$ $3$ $\nu=1/4$ $35\%$ $2\%$ $1$ $\nu=1/13$ $35\%$ $14\%$ 2
$\gamma=(1,1)$ $2$ $2$ $\nu=1/10$ $40\%$ $9\%$ $3$ $\nu=5/18$ $45\%$ $6\%$ $1$
$\nu=1/9$ $56\%$ $38\%$ $\gamma=(-0.1,5)$ $3$ $2$ $\nu=1/4$ $60\%$ $30\%$ $3$
$\nu=2/5$ $65\%$ $25\%$ $\gamma=(0.5,-1,2.5,1.2)$ $1$ $\lambda_{C}=0.12$
$25\%$ $10\%$ $\beta=(-1,0.5,1.5)$ $1$ $2$ $\lambda_{C}=0.25$ $30\%$ $6\%$
$\tau_{0}=30$, $\tau=35$ $3$ $\lambda_{C}=0.45$ $35\%$ $4\%$
$\gamma=(1,2,1.8,0.5)$ $1$ $\lambda_{C}=0.2$ $35\%$ $16\%$ 3 $\beta=(1,0.5,2)$
$2$ $2$ $\lambda_{C}=0.5$ $40\%$ $9\%$ $\tau_{0}=6$, $\tau=8$ $3$
$\lambda_{C}=0.8$ $45\%$ $6\%$ $\gamma=(-0.8,1.3,1.5,-0.2)$ $1$
$\lambda_{C}=0.3$ $55\%$ $24\%$ $\beta=(1,-0.1,0.8)$ $3$ $2$ $\lambda_{C}=0.7$
$59\%$ $14\%$ $\tau_{0}=5$, $\tau=7$ $3$ $\lambda_{C}=1.3$ $65\%$ $8\%$
$\gamma=(0.6,-1,1,2.5,1.2)$ $1$ $\lambda_{C}=0.1$ $25\%$ $11\%$
$\beta=(-0.8,0.3,0.5)$ $1$ $2$ $\lambda_{C}=0.22$ $30\%$ $7\%$ $\tau_{0}=14$,
$\tau=16$ $3$ $\lambda_{C}=0.35$ $35\%$ $5\%$ $\gamma=(0.45,0.5,2{,}1,0.5)$
$1$ $\lambda_{C}=0.15$ $35\%$ $11\%$ 4 $\beta=(1,0.5,2)$ $2$ $2$
$\lambda_{C}=0.35$ $40\%$ $7\%$ $\tau_{0}=18$, $\tau=20$ $3$ $\lambda_{C}=0.6$
$45\%$ $5\%$ $\gamma=(-0.22,0.3,-0.4,0.5,-0.2)$ $1$ $\lambda_{C}=0.2$ $55\%$
$30\%$ $\beta=(0.4,-0.1,0.5)$ $3$ $2$ $\lambda_{C}=0.4$ $59\%$ $20\%$
$\tau_{0}=6$, $\tau=8$ $3$ $\lambda_{C}=0.7$ $65\%$ $12\%$

Table 3: Bias, variance and MSE of $\hat{\gamma}$ and $\hat{\beta}$ for smcure
(second rows) and our approach (first rows) in Model 1 and 2.

Cens. level 1 Cens. level 2 Cens. level 3 Mod. n scen. Par. Bias Var. MSE Bias
Var. MSE Bias Var. MSE 1 $200$ $1$ $\gamma_{1}$ $0.001$ $0.060$ $0.060$
$0.020$ $0.065$ $0.065$ $0.005$ $0.078$ $0.078$ $0.021$ $0.063$ $0.063$
$0.050$ $0.068$ $0.071$ $0.044$ $0.084$ $0.086$ $\gamma_{2}$ $-0.034$ $0.164$
$0.165$ $-0.014$ $0.202$ $0.202$ $-0.051$ $0.209$ $0.212$ $0.026$ $0.173$
$0.173$ $0.067$ $0.222$ $0.226$ $0.044$ $0.229$ $0.230$ $\beta$ $0.008$
$0.028$ $0.028$ $0.015$ $0.029$ $0.029$ $0.013$ $0.034$ $0.035$ $0.007$
$0.028$ $0.028$ $0.012$ $0.029$ $0.029$ $0.009$ $0.035$ $0.035$ $3$
$\gamma_{1}$ $-0.001$ $0.059$ $0.059$ $0.009$ $0.065$ $0.065$ $-0.014$ $0.091$
$0.092$ $0.010$ $0.064$ $0.064$ $0.029$ $0.074$ $0.075$ $0.037$ $0.113$
$0.115$ $\gamma_{2}$ $-0.034$ $0.536$ $0.537$ $-0.111$ $0.595$ $0.608$
$-0.085$ $0.809$ $0.816$ $0.201$ $0.649$ $0.689$ $0.218$ $0.768$ $0.816$
$0.400$ $1.146$ $1.306$ $\beta$ $0.011$ $0.090$ $0.090$ $0.024$ $0.109$
$0.110$ $0.014$ $0.128$ $0.128$ $0.007$ $0.091$ $0.091$ $0.014$ $0.110$
$0.110$ $-0.001$ $0.129$ $0.129$ $400$ $1$ $\gamma_{1}$ $0.001$ $0.028$
$0.028$ $0.007$ $0.032$ $0.032$ $0.001$ $0.037$ $0.037$ $0.015$ $0.029$
$0.030$ $0.027$ $0.033$ $0.034$ $0.024$ $0.039$ $0.039$ $\gamma_{2}$ $-0.024$
$0.083$ $0.084$ $-0.004$ $0.088$ $0.088$ $-0.018$ $0.107$ $0.107$ $0.021$
$0.087$ $0.087$ $0.049$ $0.093$ $0.095$ $0.041$ $0.111$ $0.113$ $\beta$
$0.003$ $0.013$ $0.013$ $0.007$ $0.015$ $0.015$ $0.002$ $0.016$ $0.016$
$0.002$ $0.013$ $0.013$ $0.005$ $0.015$ $0.015$ $0.000$ $0.016$ $0.016$ $3$
$\gamma_{1}$ $-0.004$ $0.029$ $0.029$ $-0.004$ $0.030$ $0.030$ $-0.007$
$0.048$ $0.048$ $0.002$ $0.030$ $0.030$ $0.009$ $0.033$ $0.033$ $0.015$
$0.053$ $0.053$ $\gamma_{2}$ $-0.050$ $0.237$ $0.239$ $-0.080$ $0.312$ $0.318$
$-0.134$ $0.432$ $0.450$ $0.111$ $0.260$ $0.273$ $0.142$ $0.361$ $0.381$
$0.167$ $0.491$ $0.519$ $\beta$ $-0.003$ $0.039$ $0.039$ $0.024$ $0.051$
$0.052$ $0.013$ $0.071$ $0.071$ $-0.007$ $0.039$ $0.039$ $0.017$ $0.051$
$0.052$ $0.000$ $0.071$ $0.071$ 2 $200$ $1$ $\gamma_{1}$ $0.004$ $0.040$
$0.040$ $0.020$ $0.045$ $0.045$ $-0.016$ $0.060$ $0.060$ $0.017$ $0.040$
$0.040$ $0.058$ $0.047$ $0.050$ $0.083$ $0.079$ $0.086$ $\gamma_{2}$ $0.001$
$0.039$ $0.039$ $-0.022$ $0.042$ $0.043$ $-0.027$ $0.055$ $0.056$ $0.016$
$0.040$ $0.040$ $0.008$ $0.047$ $0.047$ $0.029$ $0.072$ $0.073$ $\beta$
$0.006$ $0.011$ $0.011$ $0.000$ $0.014$ $0.014$ $0.011$ $0.015$ $0.015$
$0.005$ $0.011$ $0.011$ $-0.002$ $0.014$ $0.014$ $0.004$ $0.016$ $0.016$ $3$
$\gamma_{1}$ $-0.016$ $0.071$ $0.071$ $-0.057$ $0.065$ $0.068$ $-0.139$
$0.083$ $0.102$ $0.029$ $0.092$ $0.092$ $0.051$ $0.119$ $0.121$ $0.024$
$0.175$ $0.176$ $\gamma_{2}$ $-0.468$ $0.723$ $0.942$ $-0.943$ $0.823$ $1.713$
$-1.348$ $0.829$ $2.646$ $0.364$ $0.926$ $1.058$ $0.495$ $1.453$ $1.698$
$0.596$ $2.128$ $2.482$ $\beta$ $0.017$ $0.035$ $0.035$ $0.022$ $0.039$
$0.039$ $0.036$ $0.052$ $0.054$ $0.014$ $0.035$ $0.035$ $0.017$ $0.040$
$0.040$ $0.025$ $0.053$ $0.054$ $400$ $1$ $\gamma_{1}$ $0.011$ $0.019$ $0.019$
$0.019$ $0.023$ $0.023$ $0.002$ $0.032$ $0.032$ $0.018$ $0.019$ $0.019$
$0.037$ $0.023$ $0.025$ $0.047$ $0.034$ $0.036$ $\gamma_{2}$ $-0.002$ $0.018$
$0.018$ $-0.010$ $0.023$ $0.023$ $-0.019$ $0.027$ $0.028$ $0.009$ $0.018$
$0.018$ $0.007$ $0.025$ $0.025$ $0.008$ $0.032$ $0.032$ $\beta$ $0.000$
$0.006$ $0.006$ $0.004$ $0.006$ $0.006$ $0.003$ $0.008$ $0.008$ $0.000$
$0.006$ $0.006$ $0.002$ $0.006$ $0.006$ $0.000$ $0.008$ $0.008$ $3$
$\gamma_{1}$ $-0.015$ $0.031$ $0.031$ $-0.071$ $0.034$ $0.039$ $-0.086$
$0.041$ $0.048$ $0.014$ $0.037$ $0.038$ $0.001$ $0.050$ $0.050$ $0.047$
$0.072$ $0.074$ $\gamma_{2}$ $-0.444$ $0.330$ $0.527$ $-0.802$ $0.410$ $1.053$
$-1.191$ $0.463$ $1.881$ $0.149$ $0.364$ $0.386$ $0.244$ $0.557$ $0.616$
$0.325$ $0.739$ $0.845$ $\beta$ $0.007$ $0.016$ $0.016$ $0.015$ $0.019$
$0.020$ $0.017$ $0.024$ $0.024$ $0.004$ $0.016$ $0.016$ $0.010$ $0.019$
$0.020$ $0.010$ $0.024$ $0.024$

For each setting we consider samples of size $n=200,400,1000$. This leads to a
total of $108$ settings ($4$ models, $3$ scenarios for the cure rate, $3$
censoring levels and $3$ sample sizes). In this way, we hope to address a
number of issues such as the effect of the cure proportion, the sample size,
amount and type of censoring, covariates (number, relation between $X$ and $Z$
and their distribution). For each configuration $1020$ datasets were generated
and the estimators of $\beta_{0}$ and $\gamma_{0}$ were computed through
smcure and our method. We report the bias, variance and mean squared error
(MSE) of the estimators, computed after omitting the lowest and the highest
$1\%$ of the estimators (for stability of the reported results) and rounded to
three decimals. Tables 3-5 show some of the results, while the rest can be
found in the online Supplementary Material. We aim to provide a ready-to-use
method that works well in practice without needing to think about how to
choose the kernel function or the bandwidth. Hence, we illustrate the
performance of the method for some standard and commonly used choices. The
kernel function $k$ is taken to be the Epanechnikov kernel
$k(u)=(3/4)(1-u^{2})\mathds{1}_{\\{|u|\leq 1\\}}$. We use the cross-validation
bandwidth (implemented in the R package np) for kernel estimators of
conditional distribution functions, in our case for estimation of
$H=H_{0}+H_{1}$ given the continuous covariates (affecting the incidence). In
addition, we restrict to the interval $[0,Y_{(m)}]$, where $Y_{(m)}$ is the
last observed event time since the estimator of the cure probability
$\hat{\pi}$ in (9) is essentially a product over values of $t$ that are equal
to the observed event times. This means that we use the cross-validation
bandwidth for estimation of the conditional distribution $H(t|x)$ for $t\leq
Y_{(m)}$. This choice of bandwidth improves significantly the performance of
the estimators, compared to the cross-validation bandwidth on the whole
interval $[0,\tau]$, in situations with a large percentage of observations in
the plateau, while it leads to little difference otherwise.

Table 4:  Bias, variance and MSE of $\hat{\gamma}$ for smcure (second rows)
and our approach (first rows) in Model 3.

Cens. level 1 Cens. level 2 Cens. level 3 Mod. n scen. Par. Bias Var. MSE Bias
Var. MSE Bias Var. MSE 3 $200$ $1$ $\gamma_{1}$ $0.025$ $0.147$ $0.147$
$0.010$ $0.192$ $0.192$ $-0.008$ $0.243$ $0.243$ $0.034$ $0.147$ $0.148$
$0.034$ $0.191$ $0.192$ $0.062$ $0.249$ $0.253$ $\gamma_{2}$ $-0.045$ $0.042$
$0.044$ $-0.078$ $0.049$ $0.055$ $-0.085$ $0.059$ $0.066$ $-0.077$ $0.050$
$0.056$ $-0.122$ $0.065$ $0.080$ $-0.148$ $0.092$ $0.144$ $\gamma_{3}$ $0.081$
$0.366$ $0.373$ $0.074$ $0.485$ $0.491$ $0.029$ $0.536$ $0.537$ $0.174$
$0.397$ $0.427$ $0.266$ $0.574$ $0.644$ $0.309$ $0.799$ $0.895$ $\gamma_{4}$
$-0.046$ $0.326$ $0.373$ $-0.160$ $0.412$ $0.437$ $-0.289$ $0.453$ $0.537$
$0.087$ $0.366$ $0.374$ $0.089$ $0.528$ $0.535$ $0.186$ $0.908$ $0.943$ $3$
$\gamma_{1}$ $-0.059$ $0.161$ $0.164$ $-0.091$ $0.258$ $0.266$ $-0.223$
$0.419$ $0.468$ $-0.053$ $0.163$ $0.166$ $-0.071$ $0.261$ $0.266$ $-0.138$
$0.524$ $0.543$ $\gamma_{2}$ $0.018$ $0.046$ $0.046$ $0.026$ $0.063$ $0.064$
$0.086$ $0.088$ $0.096$ $0.080$ $0.052$ $0.058$ $0.121$ $0.080$ $0.095$
$0.252$ $0.170$ $0.233$ $\gamma_{3}$ $0.060$ $0.235$ $0.238$ $0.076$ $0.366$
$0.372$ $0.135$ $0.517$ $0.535$ $0.091$ $0.242$ $0.251$ $0.135$ $0.375$
$0.393$ $0.228$ $0.642$ $0.694$ $\gamma_{4}$ $-0.030$ $0.202$ $0.203$ $-0.040$
$0.292$ $0.293$ $-0.081$ $0.479$ $0.486$ $-0.027$ $0.205$ $0.205$ $-0.017$
$0.277$ $0.277$ $-0.037$ $0.534$ $0.535$ $400$ $1$ $\gamma_{1}$ $0.016$
$0.074$ $0.074$ $0.021$ $0.091$ $0.092$ $0.003$ $0.128$ $0.128$ $0.017$
$0.072$ $0.073$ $0.022$ $0.082$ $0.082$ $0.023$ $0.108$ $0.108$ $\gamma_{2}$
$-0.026$ $0.019$ $0.019$ $-0.039$ $0.023$ $0.025$ $-0.070$ $0.032$ $0.037$
$-0.042$ $0.020$ $0.021$ $-0.049$ $0.025$ $0.027$ $-0.081$ $0.035$ $0.041$
$\gamma_{3}$ $0.039$ $0.194$ $0.195$ $0.028$ $0.219$ $0.220$ $0.026$ $0.298$
$0.298$ $0.093$ $0.190$ $0.198$ $0.097$ $0.206$ $0.215$ $0.158$ $0.297$
$0.322$ $\gamma_{4}$ $0.010$ $0.171$ $0.171$ $-0.091$ $0.193$ $0.201$ $-0.178$
$0.276$ $0.307$ $0.070$ $0.177$ $0.182$ $0.038$ $0.198$ $0.200$ $0.088$
$0.289$ $0.297$ $3$ $\gamma_{1}$ $-0.023$ $0.089$ $0.089$ $-0.051$ $0.118$
$0.121$ $-0.124$ $0.212$ $0.228$ $-0.029$ $0.092$ $0.093$ $-0.032$ $0.112$
$0.113$ $-0.062$ $0.200$ $0.204$ $\gamma_{2}$ $0.003$ $0.023$ $0.023$ $0.006$
$0.033$ $0.033$ $0.042$ $0.048$ $0.050$ $0.042$ $0.023$ $0.025$ $0.057$
$0.034$ $0.037$ $0.104$ $0.055$ $0.066$ $\gamma_{3}$ $0.010$ $0.113$ $0.113$
$0.042$ $0.166$ $0.168$ $0.090$ $0.276$ $0.284$ $0.039$ $0.111$ $0.111$
$0.060$ $0.152$ $0.156$ $0.108$ $0.250$ $0.262$ $\gamma_{4}$ $0.012$ $0.110$
$0.110$ $-0.021$ $0.131$ $0.131$ $-0.047$ $0.220$ $0.223$ $0.014$ $0.111$
$0.111$ $-0.018$ $0.117$ $0.118$ $-0.020$ $0.183$ $0.183$

Simulations show that, for not large sample size, the new method performs
better than smcure for estimation of $\gamma_{0}$, mostly because of a smaller
variance. As the sample size increases, they tend to behave quite similarly.
On the other hand, both methods give almost the same estimates for $\beta_{0}$
and $\Lambda$. The most favorable situation for our method is when there is
little censoring among uncured observations and the censored uncured
observations are in the region of covariates that corresponds to higher cure
rate. This comes from the fact that the nonparametric estimator in (9) takes
larger values when the product has more terms equal to one. This should not be
a problem when we expect that subjects with high probability of being cured
correspond to longer survival times, meaning that it is more probable for them
to be censored compared to those with small cure probability and shorter
survival times.

Table 5:  Bias, variance and MSE of $\hat{\gamma}$ for smcure (second rows)
and our approach (first rows) in Model 4.

Cens. level 1 Cens. level 2 Cens. level 3 Mod. n scen. Par. Bias Var. MSE Bias
Var. MSE Bias Var. MSE 4 $200$ $1$ $\gamma_{1}$ $0.041$ $0.157$ $0.159$
$0.016$ $0.187$ $0.188$ $-0.010$ $0.210$ $0.210$ $0.077$ $0.178$ $0.184$
$0.096$ $0.228$ $0.238$ $0.127$ $0.285$ $0.301$ $\gamma_{2}$ $-0.017$ $0.039$
$0.039$ $-0.019$ $0.042$ $0.042$ $-0.015$ $0.049$ $0.049$ $-0.090$ $0.052$
$0.060$ $-0.125$ $0.069$ $0.085$ $-0.164$ $0.108$ $0.135$ $\gamma_{3}$
$-0.245$ $0.159$ $0.219$ $-0.281$ $0.165$ $0.244$ $-0.355$ $0.179$ $0.305$
$0.064$ $0.244$ $0.249$ $0.084$ $0.304$ $0.311$ $0.114$ $0.395$ $0.408$
$\gamma_{4}$ $-0.068$ $0.331$ $0.336$ $-0.162$ $0.385$ $0.411$ $-0.285$
$0.443$ $0.524$ $0.171$ $0.401$ $0.430$ $0.241$ $0.561$ $0.619$ $0.314$
$0.842$ $0.941$ $\gamma_{5}$ $-0.095$ $0.301$ $0.310$ $-0.234$ $0.349$ $0.404$
$-0.366$ $0.371$ $0.505$ $0.106$ $0.363$ $0.375$ $0.143$ $0.509$ $0.529$
$0.177$ $0.693$ $0.724$ $3$ $\gamma_{1}$ $-0.044$ $0.079$ $0.081$ $-0.079$
$0.095$ $0.101$ $-0.148$ $0.132$ $0.154$ $0.000$ $0.079$ $0.079$ $0.003$
$0.096$ $0.096$ $0.009$ $0.141$ $0.141$ $\gamma_{2}$ $0.018$ $0.007$ $0.008$
$0.024$ $0.008$ $0.009$ $0.041$ $0.010$ $0.012$ $0.015$ $0.008$ $0.008$
$0.017$ $0.009$ $0.010$ $0.028$ $0.013$ $0.014$ $\gamma_{3}$ $0.034$ $0.066$
$0.067$ $0.046$ $0.073$ $0.075$ $0.067$ $0.087$ $0.091$ $-0.025$ $0.080$
$0.080$ $-0.033$ $0.091$ $0.092$ $-0.041$ $0.120$ $0.122$ $\gamma_{4}$ $0.041$
$0.102$ $0.104$ $0.054$ $0.125$ $0.128$ $0.082$ $0.166$ $0.173$ $0.022$
$0.100$ $0.101$ $0.023$ $0.126$ $0.126$ $0.026$ $0.179$ $0.180$ $\gamma_{5}$
$-0.031$ $0.103$ $0.104$ $-0.034$ $0.120$ $0.121$ $-0.054$ $0.159$ $0.162$
$-0.016$ $0.099$ $0.099$ $-0.013$ $0.115$ $0.115$ $-0.018$ $0.150$ $0.151$
$400$ $1$ $\gamma_{1}$ $0.013$ $0.067$ $0.067$ $0.015$ $0.079$ $0.080$ $0.005$
$0.097$ $0.097$ $0.024$ $0.067$ $0.068$ $0.037$ $0.080$ $0.082$ $0.043$
$0.101$ $0.103$ $\gamma_{2}$ $-0.001$ $0.017$ $0.017$ $-0.003$ $0.020$ $0.020$
$-0.007$ $0.023$ $0.023$ $-0.042$ $0.020$ $0.021$ $-0.055$ $0.025$ $0.028$
$-0.079$ $0.034$ $0.041$ $\gamma_{3}$ $-0.229$ $0.089$ $0.141$ $-0.207$
$0.090$ $0.133$ $-0.275$ $0.102$ $0.178$ $0.046$ $0.107$ $0.109$ $0.061$
$0.137$ $0.141$ $0.066$ $0.162$ $0.166$ $\gamma_{4}$ $-0.063$ $0.161$ $0.165$
$-0.143$ $0.175$ $0.196$ $-0.222$ $0.215$ $0.265$ $0.085$ $0.176$ $0.183$
$0.107$ $0.222$ $0.234$ $0.145$ $0.318$ $0.339$ $\gamma_{5}$ $-0.075$ $0.146$
$0.151$ $-0.192$ $0.177$ $0.214$ $-0.299$ $0.199$ $0.289$ $0.043$ $0.157$
$0.159$ $0.049$ $0.194$ $0.196$ $0.060$ $0.253$ $0.257$ $3$ $\gamma_{1}$
$-0.024$ $0.038$ $0.039$ $-0.038$ $0.047$ $0.048$ $-0.092$ $0.064$ $0.073$
$-0.003$ $0.036$ $0.036$ $0.002$ $0.044$ $0.044$ $0.004$ $0.060$ $0.060$
$\gamma_{2}$ $0.006$ $0.003$ $0.003$ $0.010$ $0.004$ $0.004$ $0.019$ $0.005$
$0.006$ $0.005$ $0.004$ $0.004$ $0.005$ $0.004$ $0.004$ $0.006$ $0.006$
$0.006$ $\gamma_{3}$ $0.028$ $0.033$ $0.034$ $0.045$ $0.038$ $0.040$ $0.065$
$0.045$ $0.049$ $-0.010$ $0.036$ $0.036$ $-0.008$ $0.042$ $0.042$ $-0.008$
$0.052$ $0.052$ $\gamma_{4}$ $0.021$ $0.052$ $0.053$ $0.032$ $0.062$ $0.063$
$0.060$ $0.088$ $0.092$ $0.015$ $0.050$ $0.050$ $0.016$ $0.060$ $0.060$
$0.019$ $0.083$ $0.083$ $\gamma_{5}$ $-0.024$ $0.049$ $0.050$ $-0.041$ $0.059$
$0.061$ $-0.048$ $0.077$ $0.079$ $-0.017$ $0.049$ $0.049$ $-0.022$ $0.056$
$0.057$ $-0.022$ $0.069$ $0.069$

This is indeed the case in Model 1 and we observe that our approach
outperforms smcure in all the scenarios. The difference between the two is
more marked when $n$ is small and the absolute value of the $\gamma$
coefficient is larger. In Model 2, the situation is more difficult because
censoring depends on the covariate in such a way that, the non-cured subjects
have the same probability of being censored independently of their cure
probability. However, for the first two scenarios the new method is still
superior. The third scenario is more problematic because the cure probability
drops very fast from almost one to almost zero, resulting in a large fraction
of uncured observations with almost zero cure probability. The presence of
censoring in this region leads to overestimation of the cure rate. If we would
take $\beta_{C}=0.1$ (meaning larger probability of being censored for higher
cure rate), then the new approach is significantly superior (see Table 6 for
$n=400$ and scenario 3). In Model 3, complications arise because of the
presence of different covariates for the incidence and latency. Hence,
subjects with higher cure rate might correspond to shorter survival times. As
a result, the previous problem might still happen and its effects are more
visible for large sample size and large censoring rate. Finally, Model 4
suggests that, even though the assumptions in Section 5 were shown to be
satisfied only for one continuous covariate, the method could be applied in
more general cases. We noticed that, when a continuous covariate affects only
the incidence and not the latency, the bandwidth selected by the np package is
often very large, meaning that it fails to capture the effect of this
covariate on the conditional distribution function. In those cases, we
truncate the selected bandwidth from above at $2$. Note that the bandwidth is
chosen for standardized covariates so the truncation level can be fixed
regardless of the distribution of the covariate. We decided to truncate at $2$
since it seems to be a kind of boundary for a ‘reasonable’ bandwidth with
standardized covariates (we do not want to externally affect chosen bandwidths
smaller than $2$ but we only replace extremely large values by $2$). However,
even when reasonable, the np bandwidth for $X_{2}$ seems to be larger than it
should, resulting in more bias in the estimator of $\gamma_{3}$. Nevertheless,
in terms of mean squared error, the method performs well for not large sample
size. If $X_{2}$ would affect also the latency, the selected bandwidth would
be more adequate and there would be no bias problems.

Table 6: Bias, variance and MSE of $\hat{\gamma}$ and $\hat{\beta}$ for smcure
and our approach in Model 2, scenario 3 when $\beta_{C}=0.1$ and $n=400$.

smcure package Our approach Parameter Bias Var MSE Bias Var MSE $\gamma_{1}$
$0.014$ $0.123$ $0.123$ $-0.058$ $0.103$ $0.106$ $\gamma_{2}$ $0.418$ $1.243$
$1.418$ $-0.535$ $0.652$ $0.937$ $\beta$ $0.001$ $0.025$ $0.025$ $0.001$
$0.027$ $0.027$

To conclude, the new approach seems to perform significantly better than
smcure when the sample size is not large and the fraction of censored
observations is not much higher than the expected cure proportion. In other
situations, both methods are comparable. However, one has to be more careful
when there is no reason to expect that the censored subjects correspond to
higher cure probabilities.

In the previous settings, we truncated the event times at $\tau_{0}$ in such a
way that condition (19) is satisfied but in practice it is unlikely to observe
event times at $\tau_{0}$. Next, we consider one additional model for which
condition (19) is not satisfied. The covariates and the parameters are as in
Model 3 described above, but the event times are generated from a Weibull
distribution on $[0,\tau_{0}]$ with $\tau_{0}=15$, i.e.

$S_{u}(t|z)=\frac{\exp\left\\{-\mu
t^{\rho}\exp(\beta^{\prime}z)\right\\}-\exp\left\\{-\mu\tau_{0}^{\rho}\exp(\beta^{\prime}z)\right\\}}{1-\exp\left\\{-\mu\tau_{0}^{\rho}\exp(\beta^{\prime}z)\right\\}}$

The censoring times are exponentially distributed as in Model 3 and truncated
at $\tau=20$. Results for sample size $n=200$ and three censoring levels are
shown in Table 7. Compared to Model 3 above, we observe that, when condition
(19) is not satisfied, presmoothing is even more superior than the smcure
estimator.

Table 7:  Bias, variance and MSE of $\hat{\gamma}$ for smcure (second rows)
and our approach (first rows) in Model 3 without condition (19).

Cens. level 1 Cens. level 2 Cens. level 3 Par. Bias Var. MSE Bias Var. MSE
Bias Var. MSE $\gamma_{1}$ $0.015$ $0.152$ $0.152$ $0.000$ $0.196$ $0.196$
$-0.032$ $0.246$ $0.247$ $0.017$ $0.150$ $0.151$ $0.027$ $0.193$ $0.194$
$0.035$ $0.260$ $0.262$ $\gamma_{2}$ $-0.054$ $0.044$ $0.047$ $-0.077$ $0.052$
$0.058$ $-0.109$ $0.064$ $0.076$ $-0.085$ $0.050$ $0.057$ $-0.119$ $0.069$
$0.083$ $-0.171$ $0.101$ $0.130$ $\gamma_{3}$ $0.087$ $0.379$ $0.386$ $0.073$
$0.450$ $0.456$ $0.045$ $0.578$ $0.580$ $0.197$ $0.423$ $0.462$ $0.249$
$0.561$ $0.623$ $0.343$ $0.885$ $1.002$ $\gamma_{4}$ $-0.010$ $0.339$ $0.339$
$-0.106$ $0.373$ $0.385$ $-0.228$ $0.498$ $0.550$ $0.125$ $0.364$ $0.380$
$0.156$ $0.523$ $0.548$ $0.260$ $1.513$ $1.581$

Finally we conclude with a remark about the computational aspect. The proposed
approach is computationally more intensive than the MLE mainly because of the
bandwidth selection through a cross-validation procedure. For example, for one
iteration in Model 3 with sample size $200$ and $400$, smcure computes the
estimates in $0.7$ and $0.8$ seconds respectively, while the new approach
requires $4.1$ and $23.5$ seconds (with a Core i7-8665U CPU desktop). However,
this seems still reasonable since the method is not meant for much larger
sample sizes.

##  7 Application: Melanoma study

To illustrate the practical performance, we apply the proposed estimation
procedure to two medical datasets for patients with melanoma and compare the
results with smcure. Melanoma is the third most common skin cancer type with
overall incidence rate 21.8 out of 100,000 people in the US (Cancer statistics
from the Center for Disease Control and Prevention) and according to the
American Cancer Society, $6850$ people are expected to die of melanoma in
$2020$. However, in the recent years, the chances of survival for melanoma
patients have increased due to earlier diagnosis and improvement of treatment
and surgical techniques. The 5-year survival rates based on the stage of the
cancer when it was first diagnosed are $92\%$ for localized, $65\%$ for
regional and $25\%$ for distant stage. It is also known that this disease is
more common among white people and the death rate is higher for men than
women. Even though most melanoma patients are cured by their initial
treatment, it is not possible to distinguish them from the uncured patients.
Hence, accurately estimating the probability of being cured is important in
order to plan further treatment and prevent recurrence of uncured patients.

###  7.1 Eastern Cooperative Oncology Group (ECOG) Data

We use the melanoma data (ECOG phase III clinical trial e1684) from the smcure
package [7] in order to compare our results with those of smcure. The purpose
of this study was to evaluate the effect of treatment (high dose interferon
alpha-2b regimen) as the postoperative adjuvant therapy. The event time is the
time from initial treatment to recurrence of melanoma and three covariates
have been considered: age (continuous variable centered to the mean), gender
(0=male and 1=female) and treatment (0=control and 1=treatment). The data
consists of $284$ observations (after deleting missing data) out of which
$196$ had recurrence of the melanoma cancer (around $30\%$ censoring). The
Kaplan-Meier curve is shown in Figure 1. The parameter estimates, standard
errors and corresponding p-values for the Wald test using our method and the
smcure package are given in Table 8. Standard errors are computed through
$500$ naive bootstrap samples.

Figure 1: Left panel: Kaplan-Meier survival curve for ECOG data. Right panel:
Kaplan-Meier survival curves for the treatment group (solid) and control group
(dotted) in the ECOG data. Table 8: Results for the incidence (logistic
component) and the latency (Cox PH component) from the ECOG data.

smcure package Our approach Covariates Estimates SE p-value Estimates SE
p-value incidence Intercept $1.3649$ $0.3457$ $8\cdot 10^{-5}$ $1.6697$
$0.3415$ $10^{-6}$ Age $0.0203$ $0.0159$ $0.2029$ $0.0220$ $0.0104$ $0.0344$
Gender $-0.0869$ $0.3347$ $0.7949$ $-0.3039$ $0.3448$ $0.3493$ Treatment
$-0.5884$ $0.3706$ $0.1123$ $-0.9345$ $0.3603$ $0.0095$ latency Age $-0.0077$
$0.0069$ $0.2663$ $-0.0079$ $0.0060$ $0.1861$ Gender $0.0994$ $0.1932$
$0.6067$ $0.1240$ $0.1653$ $0.4534$ Treatment $-0.1535$ $0.1715$ $0.3707$
$-0.0947$ $0.1692$ $0.5756$

We observe that, for both methods, the effects of the covariates have the same
direction. Only the intercept was found significant for the incidence with
smcure, while our method concludes that also age and treatment are
significant. In particular, the probability of recurring melanoma is higher
for the control group compared to the treatment group. This seems to be indeed
the case if we look at the Kaplan Meier survival curves for the two groups in
Figure 1. On the other hand, both methods agree that none of the covariates is
significant for the latency.

To illustrate another advantage of the new approach, we also compute the
maximum likelihood estimator with the smcure package for different choices of
the latency model. We see in Table 9 that the estimators of the incidence
component (and their significance) change depending on which variables are
included in the latency. On the other hand, the new method does not suffer
from this problem because it estimates the incidence independently of the
latency.

Table 9: Results for the incidence (logistic component) and the latency (Cox
PH component) from the ECOG data.

Model 1 Model 2 Model 3 Covariates Estimates SE p-value Estimates SE p-value
Estimates SE p-value incidence Intercept $1.3507$ $0.3001$ $7\cdot 10^{-6}$
$1.4148$ $0.3213$ $10^{-5}$ $1.4181$ $0.3073$ $4\cdot 10^{-6}$ Age $0.0164$
$0.0125$ $0.1905$ $0.0205$ $0.0154$ $0.1803$ $0.0209$ $0.0146$ $0.1528$ Gender
$-0.0265$ $0.3113$ $0.9320$ $-0.0673$ $0.3352$ $0.8407$ $-0.0222$ $0.3130$
$0.9432$ Treatment $-0.6060$ $0.3509$ $0.0842$ $-0.6773$ $0.3223$ $0.0415$
$-0.6913$ $0.3439$ $0.0444$ latency Age $-0.0074$ $0.0066$ $0.2568$ $-0.0073$
$0.0064$ $0.2579$ Gender $0.0789$ $0.1863$ $0.6719$ Treatment $-0.1324$
$0.1561$ $0.3963$

###  7.2 Surveillance, Epidemiology and End Results (SEER) database

The SEER database collects cancer incidence data from population-based cancer
registries in US. These data consist of patient demographic characteristics,
primary tumor site, tumor morphology, stage at diagnosis, length of follow up
and vital status. We select the database ‘Incidence - SEER 18 Regs Research
Data’ and extract the melanoma cancer data for the county of San Francisco in
California during the period $2004-2015$. We consider only patients with stage
at diagnosis: localized, regional and distant and exclude those with unknown
or zero follow-up time and restrict the study to white people because of the
very small number of cases from other races. The event time is death because
of melanoma. This cohort consists of $1445$ melanoma cases out of which $596$
are female and $849$ male. The age ranges from $11$ to $101$ years old, the
follow-up from $1$ to $155$ months. For most of the patients the cancer has
been diagnosed at early stage (localized), while for $101$ of them the stage
at diagnosis is ‘regional’ and only for $42$ it is ‘distant’. We aim at
evaluating how age, gender and stage at diagnosis affect the survival of
melanoma patients in this cohort. The use of cure models is justified from the
presence of a long plateau containing around $20\%$ of the observations (see
the Kaplan-Meier curve in Figure 2). Moreover, the Kaplan-Meier curves
depending on gender and stage at diagnosis in Figure 2 confirm that gender and
stage affect the cure rate.

Figure 2: Upper panel: Kaplan-Meier survival curves for SEER data. Left panel:
group division based on gender, females (solid) and males (dotted). Right
panel: group division based on cancer stage at diagnosis, localized (solid),
regional (dashed) and distant (dotted).

We checked the fit of the logistic model by comparing it with the single-index
mixture cure model proposed in [2] through the prediction error of the
incidence. More precisely, as in [2], we divide the data into a training set
and a test set of size $964$ and $481$ respectively. Using the training set,
we estimate the logistic/Cox model and the single-index/Cox model. Afterwards,
we compute the prediction error in the test set given by

$PE=-\sum_{j=1}^{481}\left\\{\hat{W}_{j}\log[1-\hat{\pi}(X_{j}^{\text{test}})]+(1-\hat{W}_{j})\log\hat{\pi}(X_{j}^{\text{test}})\right\\}$

where $\hat{\pi}(X_{j}^{\text{test}})$ and $\hat{W_{j}}$ are the predicted
cure probability and the predicted weight for the $j$th observation in the
test set, computed based on the parameter estimates (and the link function for
the single-index model) in the training set. More precisely, for the
logistic/Cox model we have
$\hat{\pi}(X_{j}^{\text{test}})=\phi(\hat{\gamma}_{n},X_{j}^{\text{test}})$
and

$\hat{W}_{j}=\Delta_{j}^{\text{test}}+(1-\Delta_{j}^{\text{test}})\frac{\hat{\pi}(X_{j}^{\text{test}})\exp\left(-\hat{\Lambda}_{n}(Y_{j}^{\text{test}})e^{\hat{\beta}^{\prime}_{n}Z_{j}^{\text{test}}}\right)}{1-\hat{\pi}(X_{j}^{\text{test}})+\hat{\pi}(X_{j}^{\text{test}})\exp\left(-\hat{\Lambda}_{n}(Y_{j}^{\text{test}})e^{\hat{\beta}^{\prime}_{n}Z_{j}^{\text{test}}}\right)}$

where $\hat{\gamma}_{n}$, $\hat{\beta}_{n}$ and $\hat{\Lambda}_{n}$ are the
estimated parameters and the estimated hazard function in the training set.
For the single-index/Cox model, the only difference is that
$\hat{\pi}(X_{j}^{\text{test}})=\hat{g}_{n}(\hat{\gamma}_{n},X_{j}^{\text{test}})$
where $\hat{g}_{n}$ is the estimated link function as in [2]. The weights
$\hat{W}_{j}$ correspond to the conditional expectation of the cure status $B$
given the observations. We find that the prediction error for the logistic
model is $98.53$, whereas for the single-index model it is $156.55$. This
means that the logistic model performs better.

Table 10:  Results for the incidence (logistic component) and the latency (Cox
PH component) from the SEER data.

smcure package Our approach Covariates Estimates SE p-value Estimates SE
p-value incidence Intercept $-4.2071$ $0.3817$ $0$ $-4.2436$ $0.3980$ $0$ Age
$0.0304$ $0.0122$ $0.0124$ $0.0328$ $0.0172$ $0.0565$ Gender $1.1318$ $0.4211$
$0.0072$ $1.2341$ $0.4792$ $0.010$ $S_{1}$ $2.6738$ $0.3702$ $5\cdot 10^{-13}$
$2.4474$ $0.4247$ $8\cdot 10^{-9}$ $S_{2}$ $4.0763$ $0.5067$ $8\cdot 10^{-16}$
$3.9426$ $0.4536$ $0$ latency Age $-0.0139$ $0.0098$ $0.1577$ $-0.0143$
$0.0106$ $0.1756$ Gender $-0.0549$ $0.4065$ $0.8925$ $-0.0871$ $0.3687$
$0.8131$ $S_{1}$ $0.5176$ $0.3993$ $0.1949$ $0.6130$ $0.3971$ $0.1226$ $S_{2}$
$1.8039$ $0.4529$ $7\cdot 10^{-5}$ $1.8623$ $0.5072$ $0.0002$

The parameter estimates, standard errors and corresponding p-values for the
Wald test using our method and the smcure package are given in Table 10.
Standard errors are computed with $500$ naive bootstrap samples. The covariate
stage is classified using two dummy Bernoulli variables $S_{1}$ and $S_{2}$,
where $S_{1}=1$ indicates the regional stage and $S_{2}=1$ indicates the
distant stage. The gender variable is equal to zero for females and one for
males. We observe that both methods agree that all the considered covariates
are significant for the incidence (with age being a borderline case for our
approach). For the latency, only being in the distant stage is found
significant with both methods. Moreover, again the effects of all the
covariates on the latency and incidence have the same direction for both
methods.

##  8 Discussion

In this paper we proposed a new estimation procedure for the mixture cure
model with a parametric form of the incidence (for example logistic) and any
semiparametric model for the latency. We investigated more in detail the
logistic/Cox model given its practical relevance. Instead of using an
iterative algorithm for dealing with the unknown cure status, this method
relies on a preliminary nonparametric estimator of the cure probabilities. We
showed through simulations that the new approach improves upon the classical
maximum likelihood estimator implemented in the package smcure, mainly for
smaller sample sizes. For the latency, both methods behave similarly. Hence,
it is of particular interest in situations in which the focus is on the
estimation of cure probabilities. The real data application on the ECOG
clinical trial also showed that the improvement in estimation can be
meaningful in practice and help detecting significant effects.

The proposed method has the advantage of direct estimation of the incidence
component, without relying on the latency, which makes it robust to latency
model misspecification. On the contrary, the smcure estimator strongly depends
on the choice of the variables for the latency and could be biased for a
misspecified Cox model. Hence, for practical reason, confronting the
estimators obtained with the two methods is valuable for confirming the
results or obtaining new insights. From the theoretical point of view, unlike
the standard maximum likelihood estimation, presmoothing allows us to obtain
consistency and asymptotic normality without requiring the ‘unrealistic’
assumption that the distribution of uncured subjects has a positive mass at
the end point of the support.

It might be argued that since the proposed method relies on smoothing, it is
more complex and the results can be affected by the choice of the kernel
function or the bandwidth. Our purpose was to show that the user doesn’t have
to think about this because the standard choices proposed in this paper
perform well in practice. In addition, since the final estimator is a
parametric one and the kernel estimator is only a preliminary step of the
procedure, the results would anyway be more stable with respect to these
choices than in a nonparametric setting. The main challenge this method faces
is extension to many continuous covariates for the incidence. We did not
deeply investigate such situations since, in that case, multiple bandwidths
have to be chosen, which can be more problematic and computationally
intensive. However, our approach based on presmoothing allows to efficiently
handle these situations if the estimator $\hat{\pi}$ is constructed in a more
adequate way. One possibility would be to construct the estimator assuming a
single-index model for the latency, which is reasonable since the final goal
is a parametric estimator. With this approach one can avoid the choice of
multiple bandwidths and perform the estimation as in the one dimensional case.
However, this problem will be addressed by future research. In this regard,
even though considering only one continuous covariate might seem restrictive
in practice, the proposed procedure constitutes the basis for further
developments of new estimators for general dimension scenarios that do not
require multidimensional smoothing.

##  9 Appendix

###  9.1 Proof of Theorem 4

We obtain the asymptotic normality of $\hat{\Lambda}_{n}$, $\hat{\beta}_{n}$
following the proof of Theorem 3 in [18]. In order to work with a one-
dimensional submodel, for $d$ in a neighbourhood of the origin, let
$\Lambda_{d}(t)=\int_{0}^{t}\\{1+dh_{1}(s)\\}\mathrm{d}\hat{\Lambda}_{n}(s)$
and $\beta_{d}=dh_{2}+\hat{\beta}_{n}$, where $h_{1}$ is a function of bounded
variation on $[0,{\tau_{0}}]$ and $h_{2}$ is a $q$-dimensional real vector.
Let $\hat{S}_{n}(\hat{\Lambda}_{n},\hat{\beta}_{n})(h_{1},h_{2})$ denote the
derivative of ${\hat{l}_{n}}(\Lambda_{d},\beta_{d})$ (defined in (14)) with
respect to $d$ and evaluated at $d=0$. We have

$\begin{split}\hat{S}_{n}(\hat{\Lambda}_{n},\hat{\beta}_{n})(h_{1},h_{2})&=\frac{1}{n}\sum_{i=1}^{n}\Delta_{i}\mathds{1}_{\\{Y_{i}<\tau_{0}\\}}\left[h_{1}(Y_{i})+h^{\prime}_{2}Z_{i}\right]\\\
&\quad-\frac{1}{n}\sum_{i=1}^{n}\left\\{\Delta_{i}+(1-\Delta_{i}){\mathds{1}_{\\{Y_{i}\leq\tau_{0}\\}}}g_{i}(Y_{i},\hat{\Lambda}_{n},\hat{\beta}_{n},\hat{\gamma}_{n})\right\\}\\\
&\qquad\qquad\qquad\quad\times\left\\{e^{\hat{\beta}^{\prime}_{n}Z_{i}}\int_{0}^{Y_{i}}h_{1}(s)\mathrm{d}\hat{\Lambda}_{n}(s)+e^{\hat{\beta}^{\prime}_{n}Z_{i}}\hat{\Lambda}_{n}(Y_{i})h^{\prime}_{2}Z_{i}\right\\},\end{split}$

where $g_{j}$ is defined in (17) and $\hat{\gamma}_{n}$ is the maximizer of
(10). Let $\Upsilon_{n}=(\hat{\Lambda}_{n},\hat{\beta}_{n})$ and
$\Upsilon_{0}=(\Lambda_{0},\beta_{0})$. Furthermore, denote by $S$ the
asymptotic version of $\hat{S}_{n}$:

$\begin{split}S(\Lambda,\beta)(h_{1},h_{2})&=\mathbb{E}\bigg{[}\Delta\mathds{1}_{\\{Y<\tau_{0}\\}}\\{h_{1}(Y)+h^{\prime}_{2}Z\\}-\left\\{\Delta+(1-\Delta){\mathds{1}_{\\{Y\leq\tau_{0}\\}}}g(Y,\Lambda,\beta,\gamma_{0})\right\\}\\\
&\qquad\qquad\qquad\qquad\left.\times\left\\{e^{\beta^{\prime}Z}\int_{0}^{Y}h_{1}(s)\mathrm{d}\Lambda(s)+e^{\beta^{\prime}Z}\Lambda(Y)h^{\prime}_{2}Z\right\\}\right].\end{split}$

We have $\hat{S}_{n}(\Upsilon_{n})=0$ and $S(\Upsilon_{0})=0$. The score
function $S_{n}$ and $S$ are respectively a random and a deterministic map
from $\Xi$ to $l^{\infty}(\mathcal{H}_{\mathfrak{m}})$ (the space of bounded
real-valued functions on $\mathcal{H}_{\mathfrak{m}}$), where

$\Xi=\left\\{(\Lambda,\beta)\,:\,\sup_{h\in\mathcal{H}_{\mathfrak{m}}}\left|\int_{0}^{{\tau_{0}}}h_{1}(s)\mathrm{d}\Lambda(s)+h^{\prime}_{2}\beta\right|<\infty\right\\}$

and
$\mathcal{H}_{\mathfrak{m}}=\\{h\in\mathcal{H}\,:\,\|h\|_{H}\leq\mathfrak{m}\\}$.
Here $\|h\|_{H}=\|h_{1}\|_{v}+\|h_{2}\|_{L_{1}}$,
$\|h_{2}\|_{L_{1}}=\sum_{j=1}^{q}|h_{2,j}|$,
$\|h_{1}\|_{v}=|h_{1}(0)|+V_{0}^{{\tau_{0}}}(h_{1})$ and
$V_{0}^{{\tau_{0}}}(h_{1})$ denotes the total variation of $h_{1}$ on
$[0,{\tau_{0}}]$. This means that $S_{n}$ is a random variable defined in the
abstract probability space $(\Omega,\mathcal{F},\mathbb{P})$ (where the random
vector $(B,T_{0},C,X,Z)$ is defined) with values in the space of bounded
functions $\Xi\mapsto l^{\infty}(\mathcal{H}_{\mathfrak{m}})$ with respect to
the supremum norm. The latter one is a Banach space equipped with the Borel
$\sigma$-field.

We need to show that conditions 1-4 of Theorem 4 in [18] (or Theorem 3.3.1 in
[26]) are satisfied. The main difference of the function $S$ from the one in
[18] is that here $\gamma=\gamma_{0}$ fixed. We are only considering variation
with respect to $\beta$ and not $\gamma$, so the components of $h$ that
correspond to $\gamma$ are set to zero. However, conditions 2 and 3 of Theorem
4 in [18] for $S$ can be shown in the same way as in [18]. Details about
conditions 1 and 4 can be found in the online Supplementary Material. ∎

###  9.2 Proof of Theorem 5

The logistic model for the cure probability obviously satisfies assumptions
(AN1) and (AN3). Let $\Pi$ be the space of continuously differentiable
functions $f$ from $\mathcal{X}$ to $[0,1]$ such that
$\sup_{x\in\mathcal{X}}|f^{\prime}(x)|\leq M$ and

$\sup_{x_{1},x_{2}\in\mathcal{X}}\frac{|f^{\prime}(x_{1})-f^{\prime}(x_{2})|}{|x_{1}-x_{2}|^{\xi}}\leq
M$

for some $M>0$ and $\xi\in(0,1]$. If such space is equipped with the supremum
norm, the covering numbers satisfy

$\log N\left(\epsilon,\Pi,\|\cdot\|_{\infty}\right)\leq
K\frac{1}{\epsilon^{1/(1+\xi)}}$

for some constant $K>0$ independent of $\epsilon$ (see Theorem 2.7.1 in [26]).
Obviously, for $\epsilon>1$, $\log N(\epsilon,\Pi,\|\cdot\|_{\infty})=0$.
Hence, assumption (AN2) is satisfied. It remains to check (AN4). Recall that
the estimator of the cure probability $\hat{\pi}(x)$ is the value at time
$\tau_{0}$ of the Beran estimator $\hat{S}(t|x)$, while
$\pi_{0}(x)=S(\tau_{0}|x)$. Moreover, by assumption (4), we have
$\inf_{x}H((\tau_{0},\infty)|x)>0$. From Proposition 4.1 and 4.2 in [27] it
follows that

$\displaystyle\sup_{x}\left|\hat{\pi}(x)-\pi_{0}(x)\right|$
$\displaystyle=O\left((nb)^{-1/2}(\log b^{-1})^{1/2}\right)\quad a.s.,$
$\displaystyle\sup_{x}\left|\hat{\pi}^{\prime}(x)-\pi_{0}^{\prime}(x)\right|$
$\displaystyle=O\left((nb^{3})^{-1/2}(\log b^{-1})^{1/2}\right)\quad a.s.$

and

$\sup_{x_{1},x_{2}\in\mathcal{X}}\frac{|\hat{\pi}^{\prime}(x_{1})-\pi^{\prime}_{0}(x_{1})-\hat{\pi}^{\prime}(x_{2})+\pi^{\prime}_{0}(x_{2})|}{|x_{1}-x_{2}|^{\xi/2}}=O\left(\left(nb^{3+\xi}\right)^{-1/2}(\log
b^{-1})^{1/2}\right)\quad a.s.,$

where $\xi$ is as in assumption (C1). Since $\pi_{0}$ is twice continuously
differentiable, from assumption (C1) it follows that $\hat{\pi}$ satisfies
(i,ii) of (AN4). From Theorem 3.2 of [10] (with $T=\tau_{0}$) we have
$\hat{\pi}(x)-\pi_{0}(x)=\frac{1}{n}\sum_{i=1}^{n}A_{i}(x)+R_{n}(x)$, where

$\begin{split}A_{i}(x)=-\frac{1-\phi(\gamma_{0},x)}{f_{X}(x)}\frac{1}{b}k\left(\frac{x-X_{i}}{b}\right)\left\\{\frac{\Delta_{i}\mathds{1}_{\\{Y_{i}\leq\tau_{0}\\}}}{H([Y_{i},\infty)|x)}-\int_{0}^{Y_{i}\wedge\tau_{0}}\frac{H_{1}(ds|x)}{H^{2}([s,\infty)|x)}\right\\}\end{split}$
(22)

and $\sup_{x}|R_{n}(x)|=O\left((nb)^{-3/4}(\log n)^{3/4}\right)$ a.s.. Hence

$\begin{split}&\mathbb{E}^{*}\left[\left(\hat{\pi}(X)-\pi_{0}(X)\right)\left(\frac{1}{\phi(\gamma_{0},X)}+\frac{1}{1-\phi(\gamma_{0},X)}\right)\nabla_{\gamma}\phi(\gamma_{0},X)\right]\\\
&=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}^{*}\left[A_{i}(x)\left(\frac{1}{\phi(\gamma_{0},X)}+\frac{1}{1-\phi(\gamma_{0},X)}\right)\nabla_{\gamma}\phi(\gamma_{0},X)\right]\\\
&\quad+\mathbb{E}^{*}\left[R_{n}(X)\left(\frac{1}{\phi(\gamma_{0},X)}+\frac{1}{1-\phi(\gamma_{0},X)}\right)\nabla_{\gamma}\phi(\gamma_{0},X)\right].\end{split}$

The second term on the right hand side of the previous display is bounded by
$c\sup_{x}|R_{n}(x)|=o(n^{-1/2})$ for some $c>0$ because of assumptions (C1)
and (AN1). Furthermore, from (AN1) and (AC4) and a Taylor expansion, it
follows that the generic element of the sum in the first term is equal to

$\begin{split}&-\int_{\mathcal{X}}\frac{1}{b}k\left(\frac{x-X_{i}}{b}\right)\left\\{\frac{\Delta_{i}\mathds{1}_{\\{Y_{i}\leq\tau_{0}\\}}}{H([Y_{i},\infty)|x)}-\int_{0}^{Y_{i}\wedge\tau_{0}}\frac{H_{1}(ds|x)}{H^{2}([s,\infty)|x)}\right\\}\frac{1}{\phi(\gamma_{0},x)}\nabla_{\gamma}\phi(\gamma_{0},x)\,\mathrm{d}x\\\
&=-\left\\{\frac{\Delta_{i}\mathds{1}_{\\{Y_{i}\leq\tau_{0}\\}}}{H([Y_{i},\infty)|X_{i})}-\int_{0}^{Y_{i}\wedge\tau_{0}}\frac{H_{1}(ds|X_{i})}{H^{2}([s,\infty)|X_{i})}\right\\}\frac{1}{\phi(\gamma_{0},X_{i})}\nabla_{\gamma}\phi(\gamma_{0},X_{i})+O(b^{2}).\end{split}$

Since because of (C1) we have $O(b^{2})=o(n^{-1/2})$, (AN4-iii) holds with

$\Psi(Y,\Delta,X)=-\left\\{\frac{\Delta\mathds{1}_{\\{Y\leq\tau_{0}\\}}}{H([Y,\infty)|X)}-\int_{0}^{Y\wedge\tau_{0}}\frac{H_{1}(ds|X)}{H^{2}([s,\infty)|X)}\right\\}\frac{1}{\phi(\gamma_{0},X)}\nabla_{\gamma}\phi(\gamma_{0},X).$

∎

## Acknowledgements

I. Van Keilegom and E. Musta acknowledge financial support from the European
Research Council (2016-2021, Horizon 2020 and grant agreement 694409). For the
simulations we used the infrastructure of the Flemish Supercomputer Center
(VSC).

## Supplement

Supporting information may be found in the online appendix. It contains the
proofs of Theorems 1, 2 and 3 in Section 5 and additional simulation results.

## References

  * [1] [author] Aerts, MarcM., Hens, NielN. and Simonoff, Jeffrey S.J. S. (2010). Model selection in regression based on pre-smoothing. J. Appl. Stat. 37 1455-1472.
  * [2] [author] Amico, MaïlisM., Van Keilegom, IngridI. and Legrand, CatherineC. (2019). The single-index/Cox mixture cure model. Biometrics 75 452–462.
  * [3] [author] Andersen, P. K.P. K. and Gill, R. D.R. D. (1982). Cox’s Regression Model for Counting Processes: A Large Sample Study. Ann. Stat. 10 1100–1120.
  * [4] [author] Berkson, JosephJ. and Gage, Robert PR. P. (1952). Survival curve for cancer patients following treatment. J. Am. Stat. Assoc. 47 501–515.
  * [5] [author] Boag, John WJ. W. (1949). Maximum likelihood estimates of the proportion of patients cured by cancer therapy. J. R. Stat. Soc. B 11 15–53.
  * [6] [author] Burke, KevinK. and Patilea, ValentinV. (2020). A likelihood-based approach for cure regression models. TEST 1–20.
  * [7] [author] Cai, ChaoC., Zou, YuboY., Peng, YingweiY. and Zhang, JiajiaJ. (2012). smcure: An R-Package for estimating semiparametric mixture cure models. Comput. Meth. Prog. Bio. 108 1255–1260.
  * [8] [author] Chen, XiaohongX., Linton, OliverO. and Van Keilegom, IngridI. (2003). Estimation of semiparametric models when the criterion function is not smooth. Econometrica 71 1591–1608.
  * [9] [author] Cristobal, JA CristobalJ. C., Roca, P FaraldoP. F. and Manteiga, W GonzalezW. G. (1987). A class of linear regression parameter estimators constructed by nonparametric estimation. Ann. Stat. 603–609.
  * [10] [author] Du, YunlingY. and Akritas, MGM. (2002). Uniform strong representation of the conditional Kaplan-Meier process. Math. Methods Stat. 11 152–182.
  * [11] [author] Farewell, Vernon TV. T. (1982). The use of mixture models for the analysis of survival data with long-term survivors. Biometrics 1041–1046.
  * [12] [author] Ferraty, FrédéricF., González-Manteiga, WenceslaoW., Martínez-Calvo, AdelaA. and Vieu, PhilippeP. (2012). Presmoothing in functional linear regression. Stat. Sin. 69–94.
  * [13] [author] Fleming, Thomas RT. R. and Harrington, David PD. P. (2011). Counting Processes and Survival Analysis 169\. John Wiley & Sons.
  * [14] [author] Han, XiaoxiaX. (2017). Statistical Methods for Analysis of Genetic and Survival Data with Latent Heterogeneity, PhD thesis, New York University.
  * [15] [author] Kuk, Anthony YCA. Y. and Chen, Chen-HsinC.-H. (1992). A mixture model combining logistic regression with proportional hazards regression. Biometrika 79 531–541.
  * [16] [author] Lee, Tamsin ET. E., Fisher, Diana OD. O., Blomberg, Simon PS. P. and Wintle, Brendan AB. A. (2017). Extinct or still out there? Disentangling influences on extinction and rediscovery helps to clarify the fate of species on the edge. Global Change Biol. 23 621–634.
  * [17] [author] Li, Chin-ShangC.-S. and Taylor, Jeremy MGJ. M. (2002). A semi-parametric accelerated failure time cure model. Stat. Med. 21 3235–3247.
  * [18] [author] Lu, WenbinW. (2008). Maximum likelihood estimation in the proportional hazards cure model. Ann. I. Stat. Math. 60 545–574.
  * [19] [author] Müller, Ursula UU. U. and Van Keilegom, IngridI. (2019). Goodness-of-fit tests for the cure rate in a mixture cure model. Biometrika 106 211–227.
  * [20] [author] Patilea, ValentinV. and Van Keilegom, IngridI. (2020). A general approach for cure models in survival analysis. Annals of Statistics 48 2323–2346.
  * [21] [author] Peng, YingweiY. and Dear, Keith BGK. B. (2000). A nonparametric mixture model for cure rate estimation. Biometrics 56 237–243.
  * [22] [author] Sposto, RichardR. (2002). Cure model analysis in cancer: an application to data from the Children’s Cancer Group. Statistics in medicine 21 293–312.
  * [23] [author] Stringer, SvenS., Denys, DamiaanD., Kahn, René SR. S. and Derks, Eske ME. M. (2016). What cure models can teach us about genome-wide survival analysis. Behav. Genet. 46 269–280.
  * [24] [author] Sy, Judy PJ. P. and Taylor, Jeremy MGJ. M. (2000). Estimation in a Cox proportional hazards cure model. Biometrics 56 227–236.
  * [25] [author] Taylor, Jeremy MGJ. M. (1995). Semi-parametric estimation in failure time mixture models. Biometrics 899–907.
  * [26] [author] van der Vaart, Aad W.A. W. and Wellner, John A.J. A. (1996). Weak Convergence and Empirical Processes, with Applications to Statistics. Springer Series in Statistics. Springer-Verlag, New York.
  * [27] [author] Van Keilegom, IngridI. and Akritas, Michael GM. G. (1999). Transfer of tail information in censored regression models. Ann. Stat. 27 1745–1784.
  * [28] [author] Wellner, Jon A.J. A. (1978). Limit theorems for the ratio of the empirical distribution function to the true distribution function. Z. Wahrsch. Verw. Gebiete 45 73–88.
  * [29] [author] Wycinka, EwaE. and Jurkiewicz, TomaszT. (2017). Mixture cure models in prediction of time to default: comparison with logit and Cox models. In Contemporary Trends and Challenges in Finance 221–231. Springer.
  * [30] [author] Xu, JianfengJ. and Peng, YingweiY. (2014). Nonparametric cure rate estimation with covariates. Can. J. Stat. 42 1–17.
  * [31] [author] Yamaguchi, KazuoK. (1992). Accelerated failure-time regression models with a regression model of surviving fraction: an application to the analysis of “permanent employment” in Japan. J. Amer. Stat. Assoc. 87 284–292.
  * [32] [author] Zhang, JiajiaJ. and Peng, YingweiY. (2007). A new estimation method for the semiparametric accelerated failure time mixture cure model. Stat. Med. 26 3157–3171.

A presmoothing approach

for estimation in semiparametric mixture cure models

Supplementary Material

Eni Musta∗, Valentin Patilea† and Ingrid Van Keilegom∗

∗KU Leuven, †Ensai

This supplement is organized as follows. Appendix A contains technical lemmas
and proofs. Appendix B collects additional simulation results, that were
omitted from the main paper due to page limits.

##  Appendix A Technical lemmas and proofs

######  Lemma 1.

Let $B$ be a Bernoulli random variable, $T_{0}$ a nonnegative random variable
and let $T=T_{0}$ if $B=1$ and $T=\infty$ if $B=0$. Let $X$ and $Z$ be two
real-valued random vectors. Then

$T_{0}\perp(C,X)\mid Z\quad\text{ and }B\perp(C,T_{0},Z)\mid
X\quad\Longrightarrow\quad T\perp C\mid(X,Z)$

###### Proof.

This lemma is similar to Lemma 8.1 in [2]. We provide the proof for
completeness. By elementary properties of conditional independence we have

$B\perp(C,T_{0},Z)\mid X\quad\Longleftrightarrow\quad B\perp
C\mid(X,Z,T_{0})\quad\text{and}\quad B\perp T_{0}\mid(X,Z)\quad\text{and}\quad
B\perp Z\mid X$

and

$T_{0}\perp(C,X)\mid Z\quad\Longleftrightarrow\quad T_{0}\perp
C\mid(X,Z)\quad\text{and}\quad T_{0}\perp X\mid Z.$

Then,

$C\perp B\mid(X,Z,T_{0})\quad\text{and}\quad C\perp
T_{0}\mid(X,Z)\quad\Longleftrightarrow\quad(B,T_{0})\perp C\mid(X,Z)$

The result follows from the fact that $T$ is completely determined by $B$ and
$T_{0}$. ∎

###  A.1 Identifiability with restricted survival times

For any $0<\tau^{*}\leq\tau_{0}$, let

$T_{0}^{*}=\min(T_{0},\tau^{*}),\qquad
T^{*}=BT_{0}^{*}+(1-B)\infty\qquad\text{and}\qquad C^{*}=\min(C,\tau^{*}).$

Moreover, let

$Y^{*}=\min(T^{*},C^{*})\quad\text{ and }\Delta^{*}=\mathds{1}_{\\{T^{*}\leq
C^{*}\\}}.$

A first aspect to study is the identifiability of the true values of the
parameter when $(Y,\Delta)$ is replaced by $(Y^{*},\Delta^{*})$. Here,
identifiability means that the true values $\beta_{0}$ and $\Lambda_{0}$ of
the parameters maximize the expectation of the criterion maximized to obtain
the estimators. This issue is addressed in Lemma 2. Let us introduce some
additional notation: for any $0<\tau^{*}\leq\tau_{0}$ and
$\Lambda\in\mathcal{H}$, $\Lambda_{|\tau^{*}}$ is defined as

$\Lambda_{|\tau^{*}}(t)=\Lambda(t),\;\;\forall t\in[0,\tau^{*})\quad\text{ and
}\quad\Delta\Lambda_{|\tau^{*}}(\tau^{*})=\Lambda_{|\tau^{*}}(\\{\tau^{*}\\})=1.$
(A1)

The dominating measure for the model of $T_{0}$ changes with such a stopped
cumulative hazard measure to allow for a positive mass at $\tau^{*}$. Then,
$\ell$ defined in (13) becomes

$\ell(y,d,x,z;\beta,\Lambda_{|\tau^{*}},\gamma)=\mathds{1}_{\\{y<\tau^{*}\\}}\left[d\log
f_{u}(y|z;\beta,\Lambda)\right.\\\
\left.+(1-d)\log\left\\{1-\phi(\gamma,x)+\phi(\gamma,x)S_{u}(y|z;\beta,\Lambda)\right\\}\right]\\\
+\mathds{1}_{\\{y\geq\tau^{*}\\}}\left[d\log
S_{u}(\tau^{*}|z;\beta,\Lambda)+(1-d)\log\left\\{1-\phi(\gamma,x)\right\\}\right].$
(A2)

######  Lemma 2.

Let $0<\tau^{*}\leq\tau_{0}$. Assume that for any $\tilde{\beta}\in B$ and
$\tilde{\Lambda}\in\mathcal{H}$,

$S_{u}(t|z;\tilde{\beta},\tilde{\Lambda}_{|\tau^{*}})=S_{u}(t|z;\beta_{0},\Lambda_{0|\tau^{*}}),\;\forall
t\in[0,\tau^{*})\quad\Longrightarrow\quad\tilde{\beta}=\beta_{0}\;\text{ and
}\;\tilde{\Lambda}_{|\tau^{*}}=\Lambda_{0|\tau^{*}}.$ (A3)

Then $(\beta_{0},\Lambda_{0|\tau^{*}})$ is the unique solution of

$\max_{\beta\in
B,\Lambda\in\mathcal{H}}\mathbb{E}\left[\ell(Y^{*},\Delta^{*},X,Z;\beta,\Lambda_{|\tau^{*}},\gamma_{0})\right].$
(A4)

Condition (A3) is a minimal requirement of identification of the true value of
the parameters in the model for the uncured subjects if the variable
$T_{0}\wedge C$ was observed and only the events in a subset of the support of
$T_{0}$ are considered. In the Cox PH model (A3) is guaranteed by the
requirement that $Var(Z)$ has full rank.

###### Proof of Lemma 2.

First, let

$H_{k}([0,t]|x,z)=\mathbb{P}(Y\leq t,{\Delta=k}|X=x,Z=z),\quad
k\in\\{0,1\\},\quad t\in[0,\infty),$

and let $H_{k}(dt|x,z)$ be the associated conditional measures. These
conditional measures characterize the distribution of $(Y,\Delta)$ given $X=x$
and $Z=z$. By the model and independence assumptions, for any $t\geq 0$,

$H_{1}(dt|x,z)=\phi(\gamma_{0},x)F_{C}([t,\infty)|x,z)f_{u}(t|z;\beta_{0},\Lambda_{0})dt,$
(A5)

and

$H_{0}(dt|x,z)=\\{1-\phi(\gamma_{0},x)+\phi(\gamma_{0},x)S_{u}(t|z;\beta_{0},\Lambda_{0})\\}F_{C}(dt|x,z).$
(A6)

Following an usual notation abuse, herein we treat $dt$ not just as the length
of a small interval but also as the name of the interval itself. Note that up
to additive terms which do not depend on the parameters $\beta,\Lambda$,

$(y,d)\mapsto d\log
f_{u}(y|z;\beta_{0},\Lambda_{0})+(1-d)\log\left\\{1-\phi(\gamma_{0},x)+\phi(\gamma_{0},x)S_{u}(y|z;\beta_{0},\Lambda_{0})\right\\},$

is the conditional log-density of $(Y,\Delta)$ given $X=x$ and $Z=z$. From
this and Kullback information inequality one can deduce that the expectation
of $\ell$ defined in (13) is maximized by $\beta_{0},\Lambda_{0}$ and
$\gamma_{0}$.

Let $0<\tau^{*}\leq\tau_{0}$. Note that

$H_{1}([\tau^{*},\tau_{0}]|x,z)=H_{1}([\tau^{*},\infty)|x,z)=\phi(\gamma_{0},x)\int_{[\tau^{*},\tau_{0}]}F_{C}([t,\infty)|x,z)f_{u}(t|z;\beta_{0},\Lambda_{0})dt,$

and

$H_{0}([\tau^{*},\infty)|x,z)=\phi(\gamma_{0},x)\int_{[\tau^{*},\tau_{0}]}S_{u}(t|z;\beta_{0},\Lambda_{0})F_{C}(dt|x,z)\\\
+\\{1-\phi(\gamma_{0},x)\\}F_{C}([\tau^{*},\infty)|x,z).$

Moreover,

$d(x,z;\tau^{*}):=\phi(\gamma_{0},x)\int_{[\tau^{*},\tau_{0}]}F_{C}([t,\infty)|x,z)f_{u}(t|z;\beta_{0},\Lambda_{0})dt\\\
+\phi(\gamma_{0},x)\int_{[\tau^{*},\tau_{0}]}S_{u}(t|z;\beta_{0},\Lambda_{0})F_{C}(dt|x,z)\\\
=\phi(\gamma_{0},x)F_{C}([\tau^{*},\infty)|x,z)S_{u}(\tau^{*}|z;\beta_{0},\Lambda_{0})\\\
=\mathbb{P}(T_{0}\wedge C\geq\tau^{*},B=1).$

In the limit case of no cure,
$d(x,z;\tau^{*})=H_{1}([\tau^{*},\infty)|x,z)+H_{0}([\tau^{*},\infty)|x,z)$.
By construction we have $Y^{*}=\min(Y,\tau^{*}),$ and

$\mathbb{P}(Y^{*}=\tau^{*},\Delta^{*}=1|X=x,Z=z)=d(x,z;\tau^{*}).$

Next, let

$H_{k}^{*}([0,t]|x,z)=\mathbb{P}(Y^{*}\leq t,{\Delta^{*}=k}|X=x,Z=z),\quad
k\in\\{0,1\\},\quad t\in[0,\infty),$

and let $H_{k}^{*}(dt|x,z)$ be the associated conditional measures. This means
for any $t\in[0,\tau^{*})$,

$H^{*}_{1}(dt|x,z)=H_{1}(dt|x,z)\quad\text{and}\quad
H^{*}_{0}(dt|x,z)=H_{0}(dt|x,z).$

Moreover,

$H_{1}^{*}(\\{\tau^{*}\\}|x,z)=H_{1}^{*}([\tau^{*},\infty)|x,z)=d(x,z;\tau^{*}),$

and

$H_{0}^{*}(\\{\tau^{*}\\}|x,z)=H_{0}^{*}([\tau^{*},\tau_{0}]|x,z)=\\{1-\phi(\gamma_{0},x)\\}F_{C}([\tau^{*},\infty)|x,z).$

Now, according to the inversion formulae of [6], without any reference to a
model, one can solve the set of equations

$\displaystyle H^{*}_{1}(dt|x,z)$ $\displaystyle=$
$\displaystyle\phi^{*}(x,z)F^{*}_{C}([t,\infty)|x,z)F^{*}_{u}(dt|x,z),$
$\displaystyle H_{0}^{*}(dt|x,z)$ $\displaystyle=$
$\displaystyle\\{1-\phi^{*}(x,z)+\phi^{*}(x,z)S^{*}_{u}(t|x,z)\\}F^{*}_{C}(dt|x,z),$

where $F_{u}^{*}=1-S_{u}^{*}$. Solving (A.1) for $F^{*}_{C}$, $S^{*}_{u}$ and
$\phi^{*}$, the functional $S^{*}_{u}$ is a proper survival function which
puts mass only on sets where $H_{1}^{*}$ does. Note that solving the similar
system with $H_{1},H_{0}$ instead of $H_{1}^{*},H_{0}^{*}$, one gets the true
$F_{C}$, $S_{u}$ and $\phi$. If $\Lambda_{C}^{*}$ denotes the cumulative
hazard function associated to the solution $F^{*}_{C}$, then

$\Lambda_{C}^{*}(dt|x,z)=\frac{H_{0}^{*}(dt|x,z)}{H^{*}_{1}((t,\infty)|x,z)+H_{0}^{*}([t,\infty)|x,z)},\quad
t\geq 0,$

and thus, by construction, we have $F_{C}(dt|x,z)=F^{*}_{C}(dt|x,z)$ on
$[0,\tau^{*})$, for any $x,z$. Then, by (A6) and the second equation in (A.1)
we deduce

$\phi^{*}(x,z)F^{*}_{u}(t|x,z)=\phi(\gamma_{0},x)F_{u}(t|z;\beta_{0},\Lambda_{0}),\quad\forall
t\in[0,\tau^{*}),\forall x,z.$

Next, taking into account that $S^{*}_{u}(t|x,z)=0$, $\forall t\geq\tau^{*}$,
$\forall x,z$, and integrating the second equation (A.1) on
$[\tau^{*},\infty)$, we obtain

$\\{1-\phi^{*}(x,z)\\}F^{*}_{C}([\tau^{*},\infty)|x,z)=H_{0}^{*}(\\{\tau^{*}\\}|x,z)=\\{1-\phi(x,z)\\}F_{C}([\tau^{*},\infty)|x,z).$

Since $F^{*}_{C}([0,\tau^{*})|x,z)=F_{C}([0,\tau^{*})|x,z)$, we deduce that
$\phi^{*}(x,z)=\phi(\gamma_{0},x)$ and thus

$F^{*}_{u}(t|x,z)=F_{u}(t|z;\beta_{0},\Lambda_{0})=F_{u}(t|z;\beta_{0},\Lambda_{0|\tau^{*}}),\quad\forall
t\in[0,\tau^{*}),\forall x,z.$ (A8)

The second equality in the last display is by the construction of the survival
function from the cumulative hazard function: only the values of $\Lambda_{0}$
on $[0,t]$ contribute to obtain $F_{u}(t|z;\beta_{0},\Lambda_{0})$. Since the
inversion formula necessarily yields $F^{*}_{u}([0,\tau^{*}]|x,z)\equiv 1$, we
deduce

$F^{*}_{u}(\\{\tau^{*}\\}|x,z)=S_{u}(\tau^{*}|z;\beta_{0},\Lambda_{0})=S_{u}(\tau^{*}|z;\beta_{0},\Lambda_{0|\tau^{*}}).$
(A9)

Finally, we can write

$E\left[\ell(Y^{*},\Delta^{*},X,Z;\beta,\Lambda_{|\tau^{*}},\gamma_{0})\right]\\\
=\iiint\log\ell(t,1,x,z;\beta,\Lambda_{|\tau^{*}},\gamma_{0})H^{*}_{1}(dt|x,z)G(dx,dz)\\\
+\iiint\log\ell(t,0,x,z;\beta,\Lambda_{|\tau^{*}},\gamma_{0})H_{0}^{*}(dt|x,z)G(dx,dz).$

To obtain the identifiability result it remains to apply Kullback information
inequality. More precisely, it suffices to notice that here, up to additive
terms which do not depend on the parameters, $\ell$ defined in (A2) considered
with $\beta_{0},\Lambda_{0|\tau^{*}}$ corresponds to the log-density of the
conditional law of $(Y^{*},\Delta^{*})$ given $X=x$ and $Z=z$. (Note that the
dominated measure changed as we introduce jumps at $\tau^{*}$.) This follows
from (A8) and (A9). Thus $\beta_{0},\Lambda_{0|\tau^{*}}$ is solution of the
problem (A4). The unicity of the solution is guaranteed by (A3). ∎

###  A.2 Consistency

###### Proof of Theorem 1..

We follow the idea of [7]. Since we are interested in almost sure convergence,
we work with fixed realizations of the data, $\omega$ that will lie in a set
of probability one. Let $\Omega$ be the abstract probability space where the
random vector $(B,T_{0},C,X,Z)$ is defined (for example we can take
$\Omega=\\{0,1\\}\times[0,\tau_{0}]\times[0,\tau]\times\mathcal{X}\times\mathcal{Z}$
and $(B,T_{0},C,X,Z)(\omega)=\omega$. Let $N\subset\Omega$ be a set of
probability one $\mathbb{P}(N)=1$ and fix $\omega\in N$. We will show that
each subsequence $\hat{\gamma}_{n_{k}}$ has a subsequence that converges to
$\gamma_{0}$. As a bounded sequence in $\mathbb{R}^{p}$,
$\hat{\gamma}_{n_{k}}$ has a convergent subsequence
$\hat{\gamma}_{m_{k}}\to\gamma^{*}$. It suffices to show that
$\gamma^{*}=\gamma_{0}$. Since $\hat{\gamma}_{m_{k}}$ maximizes
$\log\hat{L}_{m_{k},1}$, we have

$\begin{split}0&\leq\frac{1}{m_{k}}\log\hat{L}_{m_{k},1}(\hat{\gamma}_{m_{k}})-\frac{1}{m_{k}}\log\hat{L}_{m_{k},1}(\gamma_{0})\\\
&=\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}\left[\left\\{1-\hat{\pi}(X_{i})\right\\}\log\frac{\phi(\hat{\gamma}_{m_{k}},X_{i})}{\phi(\gamma_{0},X_{i})}+\hat{\pi}(X_{i})\log\frac{1-\phi(\hat{\gamma}_{m_{k}},X_{i})}{1-\phi(\gamma_{0},X_{i})}\right]\\\
&=\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}\left[\left\\{1-\pi_{0}(X_{i})\right\\}\log\frac{\phi(\gamma^{*},X_{i})}{\phi(\gamma_{0},X_{i})}+\pi_{0}(X_{i})\log\frac{1-\phi(\gamma^{*},X_{i})}{1-\phi(\gamma_{0},X_{i})}\right]+o(1)\end{split}$
(A10)

if $N\subset\\{\omega:\sup_{x}\left|\hat{\pi}(x)-\pi_{0}(x)\right|\to 0\\}$.
Note that the remainder term $o(1)$ in the previous display depends on
$\omega$ and converges to zero as $\hat{\pi}$ converges to $\pi_{0}$. Next we
will show that, for an appropriate choice of $N$, the first term converges to

$\mathbb{E}\left[\left\\{1-\pi_{0}(X)\right\\}\log\frac{\phi(\gamma^{*},X)}{\phi(\gamma_{0},X)}+\pi_{0}(x)\log\frac{1-\phi(\gamma^{*},X)}{1-\phi(\gamma_{0},X)}\right]$
(A11)

where the expectation is taken with respect to $X$ and
$\gamma^{*}\in\mathbb{R}^{p}$ (for a fixed $\omega$). Since here we are
dealing with a simple parametric model, this convergence follows easily from
the uniform law of large numbers. However, we follow a longer argument to
explain the idea that will be used also in the proof of Theorem 2 (where the
model is semiparametric). It is obvious, by the law of large numbers, that

$\begin{split}&\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}\left[\left\\{1-\pi_{0}(X_{i})\right\\}\log\phi(\gamma_{0},X_{i})+\pi_{0}(X_{i})\log\left(1-\phi(\gamma_{0},X_{i})\right)\right]\\\
&\to\mathbb{E}\left[\left\\{1-\pi_{0}(X)\right\\}\log\phi(\gamma_{0},X)+\pi_{0}(x)\log\left(1-\phi(\gamma_{0},X)\right)\right]\text{
a.s. }\end{split}$

and, at first sight it seems that the same holds when $\gamma_{0}$ is replaced
by $\gamma^{*}$. However, the proof is more delicate because $\gamma^{*}$
depends on $\omega$ and thus also the event of probability one where the
strong law of large numbers holds for this average. To avoid this we consider
a countable dense subset of $G$, $\\{\tilde{\gamma}_{l}\\}_{l\geq 1}$ (for
example the subset for which all components of $\gamma$ are rational numbers).
Now, consider the countable collection of the probability one sets
$\\{N_{l}\\}_{l\geq 1}$ where

$\begin{split}&\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}\left[\left\\{1-\pi_{0}(X_{i})\right\\}\log\phi(\tilde{\gamma}_{l},X_{i})+\pi_{0}(X_{i})\log\left(1-\phi(\tilde{\gamma}_{l},X_{i})\right)\right]\\\
&\to\mathbb{E}\left[\left\\{1-\pi_{0}(X)\right\\}\log\phi(\tilde{\gamma}_{l},X)+\pi_{0}(x)\log\left(1-\phi(\tilde{\gamma}_{l},X)\right)\right]\quad\forall
l\geq 1.\end{split}$

If $N\subseteq\left(\cap_{l\geq 1}N_{l}\right)$, we can write

$\begin{split}&\left|\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}\left[\left\\{1-\pi_{0}(X_{i})\right\\}\log\phi(\gamma^{*},X_{i})+\pi_{0}(X_{i})\log\left(1-\phi(\gamma^{*},X_{i})\right)\right]\right.\\\
&\quad-\mathbb{E}\left[\left\\{1-\pi_{0}(X)\right\\}\log\phi(\gamma^{*},X)+\pi_{0}(x)\log\left(1-\phi(\gamma^{*},X)\right)\right]\Bigg{|}\\\
&\leq\left|\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}\left[\left\\{1-\pi_{0}(X_{i})\right\\}\log\frac{\phi(\gamma^{*},X_{i})}{\phi(\tilde{\gamma}_{l},X_{i})}+\pi_{0}(X_{i})\log\frac{\left(1-\phi(\gamma^{*},X_{i})\right)}{\left(1-\phi(\tilde{\gamma}_{l},X_{i})\right)}\right]\right|\\\
&\quad+\left|\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}\left[\left\\{1-\pi_{0}(X_{i})\right\\}\log\phi(\tilde{\gamma}_{l},X_{i})+\pi_{0}(X_{i})\log\left(1-\phi(\tilde{\gamma}_{l},X_{i})\right)\right]\right.\\\
&\qquad-\mathbb{E}\left[\left\\{1-\pi_{0}(X)\right\\}\log\phi(\tilde{\gamma}_{l},X)+\pi_{0}(x)\log\left(1-\phi(\tilde{\gamma}_{l},X)\right)\right]\Bigg{|}\\\
&\quad+\left|\mathbb{E}\left[\left\\{1-\pi_{0}(X)\right\\}\log\frac{\phi(\gamma^{*},X)}{\phi(\tilde{\gamma}_{l},X)}+\pi_{0}(x)\log\frac{1-\phi(\gamma^{*},X)}{1-\phi(\tilde{\gamma}_{l},X)}\right]\right|.\end{split}$

Since $\tilde{\gamma}_{l}$ can be taken arbitrarily close to $\gamma^{*}$, by
properties of $\phi$ in assumptions (AC3)-(AC4), it can be easily derived
that, for an appropriate choice of $\tilde{\gamma}_{l}$, the first and the
third term on the right hand side in the previous equation converge to zero.
Moreover, the second term also converges to zero in the set of probability one
that we are considering. As a result, we can conclude that

$\begin{split}0&\leq\frac{\log\hat{L}_{m_{k},1}(\hat{\gamma}_{m_{k}})-\log\hat{L}_{m_{k},1}(\gamma_{0})}{m_{k}}\\\
&=\mathbb{E}\left[\left\\{1-\pi_{0}(X)\right\\}\log\frac{\phi(\gamma^{*},X)}{\phi(\gamma_{0},X)}+\pi_{0}(x)\log\frac{1-\phi(\gamma^{*},X)}{1-\phi(\gamma_{0},X)}\right]+o(1)\end{split}$

For each $x\in\mathcal{X}$, consider the function

$g_{x}(z)=\phi(\gamma_{0},x)\log\frac{z}{\phi(\gamma_{0},x)}+\left\\{1-\phi(\gamma_{0},x)\right\\}\log\frac{1-z}{1-\phi(\gamma_{0},x)},\quad
z\in(0,1).$

It is easy to check that $g_{x}(z)\leq 0$ and the equality holds only if
$z=\phi(\gamma_{0},x)$. Hence, the expectation in (A11) is smaller or equal to
zero. Due to the inequality in (A10), it must be equal to zero, which means
that $\phi(\gamma^{*},X)=\phi(\gamma_{0},X)$. By the identifiability
assumption (5), this is possible only if $\gamma^{*}=\gamma_{0}$. ∎

######  Lemma 3.

Assume (AC2),(AC5) hold and $\tau^{*}$ is such that (20) is satisfied. Then
$\sup_{n}\hat{\Lambda}_{n}({\tau^{*}})<\infty$ almost surely.

###### Proof.

By definition

$\hat{\Lambda}_{n}({\tau^{*}})=\frac{1}{n}\sum_{i=1}^{n}\frac{\Delta_{i}\mathds{1}_{\\{Y_{i}<{\tau^{*}}\\}}}{\frac{1}{n}\sum_{j=1}^{n}\mathds{1}_{\\{Y_{i}\leq
Y_{j}\leq\tau_{0}\\}}\exp(\hat{\beta}^{\prime}_{n}Z_{j})\left\\{{\Delta_{j}}+(1-{\Delta_{j}})g_{j}(Y_{j},\hat{\Lambda}_{n},\hat{\beta}_{n},\hat{\gamma}_{n})\right\\}}.$

From assumptions (AC2) and (AC5) we have

$\begin{split}&\frac{1}{n}\sum_{j=1}^{n}\mathds{1}_{\\{\tau^{*}\leq
Y_{j}\leq\tau_{0}\\}}\exp(\hat{\beta}^{\prime}_{n}Z_{j})\left\\{\Delta_{j}+(1-\Delta_{j})g_{j}(Y_{j},\hat{\Lambda}_{n},\hat{\beta}_{n},\hat{\gamma}_{n})\right\\}\\\
&\geq\frac{1}{n}\sum_{j=1}^{n}\Delta_{j}\mathds{1}_{\\{\tau^{*}\leq
Y_{j}\leq\tau_{0}\\}}\exp(\hat{\beta}^{\prime}_{n}Z_{j})\\\ &\geq
c\frac{1}{n}\sum_{j=1}^{n}\Delta_{j}\mathds{1}_{\\{\tau^{*}\leq
Y_{j}\leq\tau_{0}\\}},\end{split}$

for some $c>0$. Since
$\frac{1}{n}\sum_{j=1}^{n}\Delta_{j}\mathds{1}_{\\{\tau^{*}\leq
Y_{j}\leq\tau_{0}\\}}\xrightarrow{a.s.}\mathbb{P}\left(Y\geq\tau^{*},\Delta=1\right)>0$,
it follows that
$\frac{1}{n}\sum_{j=1}^{n}\Delta_{j}\mathds{1}_{\\{\tau^{*}\leq
Y_{j}\leq\tau_{0}\\}}$ is bounded from below away from zero almost everywhere.
As a result

$\begin{split}\sup_{n}\hat{\Lambda}_{n}(\tau^{*})&\leq\sup_{n}\frac{1}{n}\sum_{i=1}^{n}\frac{\Delta_{i}\mathds{1}_{\\{Y_{i}<\tau^{*}\\}}}{\frac{1}{n}\sum_{j=1}^{n}\mathds{1}_{\\{\tau^{*}\leq
Y_{j}\leq\tau_{0}\\}}\exp(\hat{\beta}^{\prime}_{n}Z_{j})\left\\{\Delta_{j}+(1-\Delta_{j})g_{j}(Y_{j},\hat{\Lambda}_{n},\hat{\beta}_{n},\hat{\gamma}_{n})\right\\}}\\\
&\leq\sup_{n}\frac{1}{n}\sum_{i=1}^{n}\frac{\Delta_{i}\mathds{1}_{\\{Y_{i}<\tau_{0}\\}}}{c\frac{1}{n}\sum_{j=1}^{n}\Delta_{j}\mathds{1}_{\\{\tau^{*}\leq
Y_{j}\leq\tau_{0}\\}}}\\\
&\leq\frac{1}{c}\left(\inf_{n}{\frac{1}{n}\sum_{j=1}^{n}\Delta_{j}\mathds{1}_{\\{\tau^{*}\leq
Y_{j}\leq\tau_{0}\\}}}\right)^{-1}\end{split}$

is bounded almost surely. Note that, if (19) is satisfied, then we can take
$\tau^{*}=\tau_{0}$. ∎

###### Proof of Theorem 2..

Let $0<\tau^{*}\leq\tau_{0}$ and

$\hat{l}_{n}^{*}(\beta,\Lambda_{|\tau^{*}},\hat{\gamma}_{n})=\frac{1}{n}\sum_{i=1}^{n}\ell(Y_{i}^{*},\Delta_{i}^{*},X_{i},Z_{i};\beta,\Lambda_{|\tau^{*}},\hat{\gamma}_{n}),$

with $\ell$ defined in (A2). If we consider the Cox PH model for the
conditional law of $T_{0}$, then

$\hat{l}^{*}_{n}(\beta,\Lambda_{|\tau^{*}},\hat{\gamma}_{n})=\frac{1}{n}\sum_{i=1}^{n}\Delta_{i}\left[\mathds{1}_{\\{Y_{i}<\tau^{*}\\}}\left\\{\log\Delta\Lambda(Y_{i})+\beta^{\prime}Z_{i}-\Lambda(Y_{i})e^{\beta^{\prime}Z_{i}}\right\\}\right]\\\
+\frac{1}{n}\sum_{i=1}^{n}(1-\Delta_{i})\mathds{1}_{\\{Y_{i}<\tau^{*}\\}}\log\left\\{1-\phi(\hat{\gamma}_{n},X_{i})+\phi(\hat{\gamma}_{n},X_{i})\exp\left(-\Lambda(Y_{i})e^{\beta^{\prime}Z_{i}}\right)\right\\}\\\
-\frac{\Lambda(\tau^{*}-)}{n}\sum_{i=1}^{n}\mathds{1}_{\\{Y_{i}\geq\tau^{*}\\}}\mathds{1}_{\\{B_{i}=1\\}}e^{\beta^{\prime}Z_{i}}+\frac{1}{n}\sum_{i=1}^{n}\mathds{1}_{\\{Y_{i}\geq\tau^{*}\\}}\mathds{1}_{\\{B_{i}=0\\}}\log\left\\{1-\phi(\hat{\gamma}_{n},X_{i})\right\\},$

and has to be maximized with respect to $\beta$ and $\Lambda$ in the class of
step functions $\Lambda$ with jumps of size $\Delta\Lambda$ at the event times
in $[0,\tau^{*})$. As in [5], it can be shown that the maximizer
$(\hat{\Lambda}^{*}_{n},\hat{\beta}^{*}_{n})$ of $\hat{l}^{*}_{n}$ exists and
it is finite. Moreover, for $t\in[0,\tau_{0}]$,
$\hat{\Lambda}_{n}^{*}=\Lambda^{*}_{n,\hat{\beta}^{*}_{n},\hat{\gamma}_{n}}$
where

$\Lambda^{*}_{n,\beta,\gamma}(t)=\frac{1}{n}\sum_{i=1}^{n}\frac{\Delta_{i}\mathds{1}_{\\{Y_{i}{\leq
t,Y_{i}<\tau^{*}}\\}}}{\frac{1}{n}\sum_{j=1}^{n}\mathds{1}_{\\{{Y_{j}\geq
Y_{i}}\\}}\exp(\beta^{\prime}Z_{j})\left\\{{\Delta_{j}^{*}}+(1-{\Delta_{j}})\mathds{1}_{\\{Y_{j}<\tau^{*}\\}}g_{j}(Y_{j},\Lambda^{*}_{n,\beta,\gamma},\beta,\gamma)\right\\}},$

$\Delta^{*}=\mathds{1}_{\\{T_{0}^{*}\leq
C^{*}\\}}=\Delta\mathds{1}_{\\{Y_{j}<\tau^{*}\\}}+\mathds{1}_{\\{Y_{i}\geq\tau^{*}\\}}\mathds{1}_{\\{B_{i}=1\\}}$
and $g_{j}(t,\Lambda,\beta,\gamma)$ defined in (17).

Let

$\tilde{\Lambda}_{0,n}(t)=\frac{1}{n}\sum_{i=1}^{n}\frac{\Delta_{i}\mathds{1}_{\\{Y_{i}\leq
t{,Y_{i}<\tau_{0}}\\}}}{\frac{1}{n}\sum_{j=1}^{n}\mathds{1}_{\\{Y_{j}\geq
Y_{i}{,Y_{j}\leq\tau_{0}}\\}}\exp(\beta^{\prime}_{0}Z_{j})\left\\{{\Delta_{j}}+(1-{\Delta_{j}})g_{j}(Y_{j},\Lambda_{0},\beta_{0},\gamma_{0}))\right\\}}.$
(A12)

We want to prove that $\hat{\beta}_{n}{\xrightarrow{a.s.}}\beta_{0}$, and
$\sup_{t\in[0,\bar{\tau}]}|\hat{\Lambda}_{n}(t)-\Lambda_{0}(t)|{\xrightarrow{a.s.}}0$
for any $\bar{\tau}<\tau_{0}$. We suppose that the previous statement is
false, i.e $\hat{\beta}_{n}$ does not converge almost surely to $\beta_{0}$ or
there exists $\bar{\tau}$ such that
$\sup_{t\in[0,\bar{\tau}]}|\hat{\Lambda}_{n}(t)-\Lambda_{0}(t)|$ does not
converge to zero almost surely. This means that, there exist $\epsilon>0$ and
$\bar{\tau}<\tau_{0}$ such that

$\mathbb{P}[A_{1}(\bar{\tau},\epsilon)]>0,\;\;\text{ with
}\;\;A_{1}(\bar{\tau},\epsilon)=\left\\{\limsup_{n\rightarrow\infty}\left[\left\|\hat{\beta}_{n}-\beta_{0}\right\|+\sup_{t\in[0,\bar{\tau}]}\left|\hat{\Lambda}_{n}(t)-\Lambda_{0}(t)\right|\right]>\epsilon\right\\}.$

On the other hand, since $(\hat{\Lambda}_{n},\hat{\beta}_{n})$ maximizes
$\hat{\ell}_{n}(\Lambda,\beta,\hat{\gamma}_{n})$, for any realization $\omega$
of the data we have

$\hat{l}_{n}(\hat{\beta}_{n},\hat{\Lambda}_{n},\hat{\gamma}_{n})-\hat{l}_{n}(\beta_{0},\tilde{\Lambda}_{0,n},\hat{\gamma}_{n})\geq
0.$ (A13)

Then the idea for creating the contradiction is to show that the previous
inequality is not satisfied for any $\omega$ in some event of positive
probability. We argue for a fixed realization $\omega$ of the data. As a
bounded sequence in $\mathbb{R}^{q}$, $\hat{\beta}_{n}$ has a convergent
subsequence $\hat{\beta}_{n_{k}}\to\bar{\beta}$. Let $(\tau_{i})_{i\geq 1}$ be
an increasing sequence such that $\lim_{i\to\infty}\tau_{i}=\tau_{0}$. Since
for all $\tau<\tau_{0}$, $\hat{\Lambda}_{n}(\tau)<\infty$ almost surely (see
Lemma 3), by Helly’s selection theorem ([1]), there exists a subsequence
$\hat{\Lambda}_{m_{k}}$ of $\hat{\Lambda}_{n_{k}}$, converging pointwise to a
function $\bar{\Lambda}$ on $[0,\tau_{1}]$. Repeating the same argument, we
can extract a further subsequence converging pointwise to a function
$\bar{\Lambda}$ on $[0,\tau_{2}]$ and so on. Hence, there exist a subsequence
$\hat{\Lambda}_{r_{k}}$ converging pointwise to a function $\bar{\Lambda}$ on
all compacts of $[0,\tau_{0}]$ that do not include $\tau_{0}$. This defines a
monotone function $\bar{\Lambda}$ on $[0,\tau_{0})$, which could be extended
at $\tau_{0}$ by taking the limit. As in Lemma 2 of [5], it can be shown that
$\bar{\Lambda}$ is absolutely continuous and pointwise convergence of monotone
functions to a continuous monotone function implies uniform convergence on
compacts. Note that the chosen subsequence and the limits $\bar{\beta}$ and
$\bar{\Lambda}$ depend on $\omega$. To keep the notation simple, in what
follows we use the index $n$ instead of the chosen subsequence $r_{k}$. For
any $\tau^{*}<\tau_{0}$, we can write

$\displaystyle\\!\\!\\!0\leq\hat{l}_{n}(\hat{\beta}_{n},\hat{\Lambda}_{n},\hat{\gamma}_{n})-\hat{l}_{n}(\beta_{0},\tilde{\Lambda}_{0,n},\hat{\gamma}_{n})$
$\displaystyle=\hat{l}_{n}^{*}(\hat{\beta}_{n},\hat{\Lambda}_{n|\tau^{*}},\hat{\gamma}_{n})+D_{1n}-\hat{l}^{*}_{n}(\beta_{0},\tilde{\Lambda}_{0,n|\tau^{*}},\hat{\gamma}_{n})-D_{2n}$
$\displaystyle=\mathbb{E}[\ell(Y^{*},\Delta^{*},X,Z;\bar{\beta},\bar{\Lambda}_{|\tau^{*}},\gamma_{0})]+D_{1n}+R_{1n}$
$\displaystyle\quad-\mathbb{E}[\ell(Y^{*},\Delta^{*},X,Z;\beta_{0},\Lambda_{0|\tau^{*}},\gamma_{0})]-D_{2n}-R_{2n},$
(A14)

where

$D_{1n}=\hat{l}_{n}(\hat{\beta}_{n},\hat{\Lambda}_{n},\hat{\gamma}_{n})-\hat{l}_{n}^{*}(\hat{\beta}_{n},\hat{\Lambda}_{n|\tau^{*}},\hat{\gamma}_{n}),$
(A15)
$D_{2n}=\hat{l}_{n}(\beta_{0},\tilde{\Lambda}_{0,n},\hat{\gamma}_{n})-\hat{l}^{*}_{n}(\beta_{0},\tilde{\Lambda}_{0,n|\tau^{*}},\hat{\gamma}_{n}),$
(A16)
$R_{1n}=\hat{l}_{n}^{*}(\hat{\beta}_{n},\hat{\Lambda}_{n|\tau^{*}},\hat{\gamma}_{n})-\mathbb{E}[\ell(Y^{*},\Delta^{*},X,Z;\bar{\beta},\bar{\Lambda}_{|\tau^{*}},\gamma_{0})],$
(A17)
$R_{2n}=\hat{l}^{*}_{n}(\beta_{0},\tilde{\Lambda}_{0,n|\tau^{*}},\hat{\gamma}_{n})-\mathbb{E}[\ell(Y^{*},\Delta^{*},X,Z;\beta_{0},\Lambda_{0|\tau^{*}},\gamma_{0})].$
(A18)

Note that the limit of $(\hat{\beta}_{n},\hat{\Lambda}_{n})$ depends on
$\omega$, but here the expectation is taken with respect to
$(Y^{*},\Delta^{*},X,Z)$ for fixed $(\bar{\beta},\bar{\Lambda})$. We now
define the event $A_{3}(\tau^{*})=\left\\{|R_{1n}-R_{2n}|\rightarrow
0\right\\}.$ By Lemma 4, for any $\tau^{*}<\tau_{0}$, we have
$\mathbb{P}[A_{1}(\bar{\tau},\epsilon)\cap
A_{3}(\tau^{*})]=\mathbb{P}[A_{1}(\bar{\tau},\epsilon)]$. Next, for
$\bar{\tau}<\tau_{0}$ and $\epsilon>0$ such that
$\mathbb{P}[A_{1}(\bar{\tau},\epsilon)]>0$, by Lemma 6 there exist $0<c_{1}<1$
and $\delta>0$ such that we have

$\begin{split}c&=\inf\Bigg{\\{}\mathbb{E}[\ell(Y^{*},\Delta^{*},X,Z;\beta_{0},\Lambda_{0|\tau^{*}},\gamma_{0})]-\mathbb{E}[\ell(Y^{*},\Delta^{*},X,Z;\beta,\Lambda_{|\tau^{*}},\gamma_{0})]:\\\
&\left.\qquad\bar{\tau}+\delta\leq\tau^{*}<\tau_{0},\quad\|\beta-\beta_{0}\|\geq
c_{1}\epsilon/2\quad\text{
or}\quad\sup\limits_{t\in[0,\bar{\tau}]}|\Lambda(t)-\Lambda_{0}(t)|\geq(1-c_{1})\epsilon/2\right\\}>0.\end{split}$

Note that if $\omega\in A_{1}(\bar{\tau},\epsilon)$ and $\bar{\beta}$ and
$\bar{\Lambda}$ are the limits for $\hat{\beta}_{n}$ and $\hat{\Lambda}_{n}$,
respectively, then necessarily, either $\|\bar{\beta}-\beta_{0}\|\geq
c_{1}\epsilon/2,$ or
$\sup\limits_{t\in[0,\bar{\tau}]}|\bar{\Lambda}(t)-\Lambda_{0}(t)|\geq(1-c_{1})\epsilon/2$,
and consequently

$\mathbb{E}[\ell(Y^{*},\Delta^{*},X,Z;\bar{\beta},\bar{\Lambda}_{|\tau^{*}},\gamma_{0})]-\mathbb{E}[\ell(Y^{*},\Delta^{*},X,Z;\beta_{0},\Lambda_{0|\tau^{*}},\gamma_{0})]\geq
c,\quad\forall\bar{\tau}{+\delta}\leq\tau^{*}<\tau_{0}.$

Finally, we define

$A_{2}(\tau^{*})=\left\\{\limsup_{n\rightarrow\infty}\left|D_{1n}-D_{2,n}\right|\leq
c/2\right\\},$

with $D_{1n}$ and $D_{2n}$ defined in (A15) and (A16), and choose
$\tau^{*}\in[\bar{\tau}+\delta,\tau_{0})$ such that

$c_{b}\left\\{\mathbb{P}(T_{0}\geq\tau^{*})\log\\{1/\mathbb{P}(T_{0}\geq\tau^{*})\\}+\mathbb{P}(C\in[\tau^{*},\tau_{0}])\right\\}<c/2,$

with $c_{b}$ the constant from Lemma 5. Then we have
$\mathbb{P}[A_{2}(\tau^{*})]=1$. Gathering facts, we deduce that by a suitable
choice of $\tau^{*}\in[\bar{\tau}+\delta,\tau_{0})$, we necessarily have
$\mathbb{P}[A_{1}(\bar{\tau},\epsilon)\cap A_{2}(\tau^{*})\cap
A_{3}(\tau^{*})]>0$. Moreover, with such a suitable $\tau^{*}$, for any
$\omega\in A_{1}(\bar{\tau},\epsilon)\cap A_{2}(\tau^{*})\cap
A_{3}(\tau^{*})$, we have

$\limsup_{n\rightarrow\infty}\left[\hat{l}_{n}(\hat{\beta}_{n},\hat{\Lambda}_{n},\hat{\gamma}_{n})-\hat{l}_{n}(\beta_{0},\tilde{\Lambda}_{0,n},\hat{\gamma}_{n})\right]\leq-c/2<0.$

We deduce that (A13) is violated on an event of positive probability, which by
definition is impossible. Thus $\hat{\beta}_{n}{\xrightarrow{a.s.}}\beta_{0}$,
and
$\sup_{t\in[0,\bar{\tau}]}|\hat{\Lambda}_{n}(t)-\Lambda_{0}(t)|{\xrightarrow{a.s.}}0$
for any $\bar{\tau}<\tau_{0}$.

If condition (19) is satisfied, we want to show in addition that
$|\hat{\Lambda}_{n}(\tau_{0})-\Lambda_{0}(\tau_{0})|{\xrightarrow{a.s.}}0$. In
that case, $\hat{\Lambda}_{n}(\tau_{0})<\infty$ almost surely and as a result,
for any realization $\omega$, there exists a subsequence
$\hat{\Lambda}_{r_{k}}$ converging to some absolutely continuous function
$\bar{\Lambda}$ uniformly on $[0,\tau_{0}]$. Since we already showed that
$|\hat{\Lambda}_{n}(t)-\Lambda_{0}(t)|{\xrightarrow{a.s.}}0$ for any
$t<\tau_{0}$ and
$\Lambda_{0}(\tau_{0})=\lim_{t\uparrow\tau_{0}}\Lambda_{0}(t)$, we necessarily
have $\bar{\Lambda}=\Lambda_{0}$ on the whole interval $[0,\tau_{0}]$. This
concludes the proof of the Theorem. ∎

######  Lemma 4.

Consider a realization of the data $\omega$ and assume that
$\hat{\beta}_{n}(\omega)\to\bar{\beta}$ and
$\hat{\Lambda}_{n}(\omega)(t)\to\bar{\Lambda}(t)$ for any $t\in[0,\tau_{0})$,
for some absolutely continuous function $\bar{\Lambda}$. Let
$0<\tau^{*}<\tau_{0}$ and let $R_{1n}$, $R_{2n}$ be defined as in (A17) and
(A18), respectively. There exists an event $A_{3}(\tau^{*})$ of probability
one such that, for any $\omega\in A_{3}(\tau^{*})$,

$R_{1n}(\omega)-R_{2n}(\omega)\to 0.$

###### Proof.

Let us consider some $0<\tau^{*}<\tau_{0}$. From Theorem 1 and Lemma 2 in [5]
it follows that the event

$A_{3}^{1}{(\tau^{*})}=\left\\{\hat{\gamma}_{n}\to\gamma_{0}\quad\text{ and
}\quad\sup_{t\in[0,\tau^{*}]}|\tilde{\Lambda}_{0,n}(t)-\Lambda_{0}(t)|\to
0\right\\}$

has probability one. Next we argue for the given realization of the data
$\omega{\in A_{3}^{1}(\tau^{*})}$ and will determine the event
$A_{3}{(\tau^{*})}$ appropriately. By the triangular inequality we can write

$\displaystyle|R_{1n}\\!-R_{2n}|$
$\displaystyle\leq\left|\left\\{\hat{l}^{*}_{n}(\hat{\beta}_{n},\hat{\Lambda}_{n|\tau*},\hat{\gamma}_{n})\\!-\hat{l}^{*}_{n}(\beta_{0},\tilde{\Lambda}_{0,n|\tau*},\hat{\gamma}_{n})\right\\}\\!-\\!\left\\{\hat{l}^{*}_{n}(\bar{\beta},\bar{\Lambda}_{|\tau*},\\!\gamma_{0})\\!-\\!\hat{l}^{*}_{n}(\beta_{0},{\Lambda}_{0|\tau*},\\!\gamma_{0})\right\\}\right|$
$\displaystyle\quad+\left|\hat{l}^{*}_{n}(\beta_{0},{\Lambda}_{0|\tau*},\\!\gamma_{0})-\mathbb{E}\left[l(Y^{*},\Delta^{*},X,Z;\beta_{0},\Lambda_{0|\tau^{*}},\\!\gamma_{0})\right]\right|$
$\displaystyle\quad+\left|\hat{l}^{*}_{n}(\bar{\beta},\bar{\Lambda}_{|\tau*},\gamma_{0})-\mathbb{E}\left[l(Y^{*},\Delta^{*},X,Z;\bar{\beta},\bar{\Lambda}_{|\tau^{*}},\gamma_{0})\right]\right|.$
(A19)

Since $\bar{\Lambda}$ is absolutely continuous, it is differentiable almost
everywhere. Let $\bar{\lambda}(t)=\mathrm{d}\bar{\Lambda}(t)/\mathrm{d}t$. By
definition we have

$\begin{split}&\hat{l}^{*}_{n}(\hat{\beta}_{n},\hat{\Lambda}_{n|\tau*},\hat{\gamma}_{n})-\hat{l}^{*}_{n}(\beta_{0},\tilde{\Lambda}_{0,n|\tau*},\hat{\gamma}_{n})\\\
&=\frac{1}{n}\sum_{i=1}^{n}\Delta_{i}\mathds{1}_{\\{Y_{i}<\tau^{*}\\}}\left\\{\log\frac{\Delta\hat{\Lambda}_{n}(Y_{i})}{\Delta\tilde{\Lambda}_{0,n}(Y_{i})}+(\hat{\beta}_{n}-\beta_{0})^{\prime}Z_{i}-\hat{\Lambda}_{n}(Y_{i})e^{\hat{\beta}^{\prime}_{n}Z_{i}}+\tilde{\Lambda}_{0,n}(Y_{i})e^{\beta^{\prime}_{0}Z_{i}}\right\\}\\\
&\qquad+\frac{1}{n}\sum_{i=1}^{n}(1-\Delta_{i})\mathds{1}_{\\{Y_{i}<\tau^{*}\\}}\log\frac{1-\phi(\hat{\gamma}_{n},X_{i})+\phi(\hat{\gamma}_{n},X_{i})\exp\left(-\hat{\Lambda}_{n}(Y_{i})e^{\hat{\beta}^{\prime}_{n}Z_{i}}\right)}{1-\phi(\hat{\gamma}_{n},X_{i})+\phi(\hat{\gamma}_{n},X_{i})\exp\left(-\tilde{\Lambda}_{0,n}(Y_{i})e^{\beta^{\prime}_{0}Z_{i}}\right)}\\\
&\qquad-\frac{1}{n}\sum_{i=1}^{n}\mathds{1}_{\\{Y_{i}\geq\tau^{*}\\}}\mathds{1}_{\\{B_{i}=1\\}}\left\\{\hat{\Lambda}_{n}(\tau^{*})e^{\hat{\beta}^{\prime}_{n}Z_{i}}-\tilde{\Lambda}_{0,n}(\tau^{*})e^{\beta^{\prime}_{0}Z_{i}}\right\\}.\end{split}$

If $\omega\in A^{1}_{3}{(\tau^{*})}$, we obtain

$\begin{split}&\hat{l}^{*}_{n}(\hat{\beta}_{n},\hat{\Lambda}_{n|\tau*},\hat{\gamma}_{n})-\hat{l}^{*}_{n}(\beta_{0},\tilde{\Lambda}_{0,n|\tau*},\hat{\gamma}_{n})\\\
&=\frac{1}{n}\sum_{i=1}^{n}\Delta_{i}\mathds{1}_{\\{Y_{i}<\tau^{*}\\}}\left\\{\log\frac{\bar{\lambda}(Y_{i})}{\lambda_{0}(Y_{i})}+(\bar{\beta}-\beta_{0})^{\prime}Z_{i}-\bar{\Lambda}(Y_{i})e^{\bar{\beta}^{\prime}Z_{i}}+\Lambda_{0}(Y_{i})e^{\beta^{\prime}_{0}Z_{i}}\right\\}\\\
&\qquad+\frac{1}{n}\sum_{i=1}^{n}(1-\Delta_{i})\mathds{1}_{\\{Y_{i}<\tau^{*}\\}}\log\frac{1-\phi(\gamma_{0},X_{i})+\phi(\gamma_{0},X_{i})\exp\left(-\bar{\Lambda}(Y_{i})e^{\bar{\beta}^{\prime}Z_{i}}\right)}{1-\phi(\gamma_{0},X_{i})+\phi(\gamma_{0},X_{i})\exp\left(-\Lambda_{0}(Y_{i})e^{\beta^{\prime}_{0}Z_{i}}\right)}\\\
&\qquad-\frac{1}{n}\sum_{i=1}^{n}\mathds{1}_{\\{Y_{i}\geq\tau^{*}\\}}\mathds{1}_{\\{B_{i}=1\\}}\left\\{\bar{\Lambda}(\tau^{*})e^{\bar{\beta}^{\prime}Z_{i}}-\Lambda_{0}(\tau^{*})e^{\beta^{\prime}_{0}Z_{i}}\right\\}+o(1)\\\
&=\hat{l}^{*}_{n}(\bar{\beta},\bar{\Lambda}_{|\tau*},\gamma_{0})-\hat{l}^{*}_{n}(\beta_{0},{\Lambda}_{0|\tau*},\gamma_{0})+o(1),\end{split}$

where the remainder term depends on $\omega$ and converges to zero. Hence, the
first term on the right hand side of (A.2) converges to zero. Let
$A_{n}^{2}{(\tau^{*})}$ be the event where

$\hat{l}^{*}_{n}(\beta_{0},{\Lambda}_{0|\tau*},\gamma_{0})\to\mathbb{E}\left[l(Y^{*},\Delta^{*},X,Z;\beta_{0},\Lambda_{0|\tau^{*}},\gamma_{0})\right]\quad\text{
as }n\to\infty.$

By the law of large numbers $\mathbb{P}[A^{2}_{n}{(\tau^{*})}]=1$, implying
that also the second term on the right hand side of (A.2) converges to zero if
$\omega\in A^{2}_{n}{(\tau^{*})}$. It remains to deal with the third term.
Note that here $(\bar{\beta},\bar{\Lambda})$ depend on $\omega$ and the
expectation is taken with respect to $(Y,\Delta,X,Z)$ for fixed
$(\bar{\beta},\bar{\Lambda})$. We have the same issue as in the proof of
Theorem 1 when dealing with the terms involving $\bar{\beta}$ and
$\bar{\Lambda}$, so we need to consider approximations by elements of a
countable dense subset of $\mathcal{B}$ and of the space of bounded,
absolutely continuous, increasing functions in $[0,\tau^{*}]$ (is separable,
so such subset exists). The same reasoning is used also in [4, 5, 7]. Hence,
there exists a countable collection of probability one sets
$\\{N_{l}\\}_{l\geq 1}$ where

$\hat{l}^{*}_{n}(\beta_{l},{\Lambda}_{l},\gamma_{0})\to\mathbb{E}\left[l(Y^{*},\Delta^{*},X,Z;\beta_{l},\Lambda_{l},\gamma_{0})\right]\quad\text{
as }n\to\infty$

and $(\beta_{l},\Lambda_{l})$ can be taken arbitrarily close to
$(\bar{\beta},\bar{\Lambda})$. As a result, if $\omega\in
A^{3}_{n}{(\tau^{*})}=\bigcap_{l\geq 1}N_{l}$, then

$\left|\hat{l}^{*}_{n}(\bar{\beta},\bar{\Lambda}_{|\tau*},\gamma_{0})-\mathbb{E}\left[l(Y^{*},\Delta^{*},X,Z;\bar{\beta},\bar{\Lambda}_{|\tau^{*}},\gamma_{0})\right]\right|\to
0.$

To conclude, we define $A_{3}{(\tau^{*})}=A^{1}_{3}{(\tau^{*})}\cap
A^{2}_{3}{(\tau^{*})}\cap A^{3}_{3}{(\tau^{*})}$ and we have
$\mathbb{P}[A_{3}{(\tau^{*})}]=1$. ∎

######  Lemma 5.

Let $D_{1n}$ and $D_{2n}$ be defined as in (A15) and (A16), respectively, for
some $\tau^{*}<\tau_{0}$. Then there exists a constant $c_{b}$ independent of
$\tau^{*}$ such that

$\mathbb{P}\left[\limsup_{n\rightarrow\infty}\left|D_{1n}-D_{2,n}\right|>c_{b}\left\\{\mathbb{P}(T_{0}\geq\tau^{*})\log\\{1/\mathbb{P}(T_{0}\geq\tau^{*})\\}+\mathbb{P}(C\in[\tau^{*},\tau_{0}])\right\\}\right]=0.$

###### Proof.

By definition, for any $\gamma$, $\beta$ and cumulative hazard function