# A 2-step estimation procedure for semiparametric mixture cure models

Eni<EMAIL_ADDRESS>[ Valentin
<EMAIL_ADDRESS>[ Ingrid Van
<EMAIL_ADDRESS>[ Korteweg de Vries Institute
for Mathematics, University of Amsterdam, Netherlands,  CREST, Ensai, France,
ORSTAT, KU Leuven, Belgium,

(0000)

###### Abstract

Cure models have been developed as an alternative modelling approach to
conventional survival analysis in order to account for the presence of cured
subjects that will never experience the event of interest. Mixture cure
models, which model separately the cure probability and the survival of
uncured subjects depending on a set of covariates, are particularly useful for
distinguishing curative from life-prolonging effects. In practice, it is
common to assume a parametric model for the cure probability and a
semiparametric model for the survival of the susceptibles. Because of the
latent cure status, maximum likelihood estimation is performed by means of the
iterative EM algorithm. Here, we focus on the cure probabilities and propose a
two-step procedure to improve upon the performance of the maximum likelihood
estimator when the sample size is not large. The new method is based on the
idea of presmoothing by first constructing a nonparametric estimator and then
projecting it into the desired parametric class. We investigate the
theoretical properties of the resulting estimator and show through an
extensive simulation study for the logistic-Cox model that it outperforms the
existing method. Practical use of the method is illustrated through two
melanoma datasets.

62N02,

cure model,

logistic model,

presmoothing,

survival analysis,

###### doi:

10.1214/154957804100000000

###### keywords:

[class=MSC]

###### keywords:

††volume: 0

, and

##  1 Introduction

Cure models are used to analyze time until occurrence of an event of interest
when a proportion of the study population is immune to that event (cured).
They are recently becoming increasingly popular in oncology as curative
treatments are now a possibility, meaning that some patients will not
experience cancer relapse/death (see for example [15, 25]). More broadly, cure
models find applications in studies of fertility ([37]), hospitalization of
COVID-19 patients ([28]), equipment failure in engineering ([22]), credit
scoring in economics ([9, 8]), etc.. What makes statistical modeling and
estimation challenging when not all subjects are susceptible to the event of
interest, is the unobserved cure status. As a consequence of a limited follow-
up period, all cured subject are observed as censored, hence mixed with the
uncured ones.

There are two main families of cure models: promotion time models and mixture
cure models (see [2] and [32] for an overview). The latter ones are
particularly attractive in practice because, by separately modeling the uncure
probability (incidence) and the survival of the susceptibles (latency) given
possibly different sets of covariates, they are able to distinguish a curative
from a life-prolonging effect. Early works on mixture cure models were fully
parametric approaches ([11, 41, 14]), while more recently semi-parametric
([35, 31, 16, 43]) and non-parametric ([40, 27, 3, 18, 19]) models have been
proposed. Among them, the semiparametric models are often used in practice as
a reasonable compromise between flexibility and simplicity. These models
assume a parametric form of the incidence and a semiparametric form for the
latency, with the most common choice being the mixture of the logistic with
the Cox proportional hazards model (e.g. [42, 33, 39]).

Estimation in the logistic-Cox or in general semi-parametric mixture cure
models is mostly carried out via the Expectation-Maximization algorithm
because of the latent cure status. Such estimators were proposed in [31] and
[35] for the logistic-Cox model; in [16, 43, 21] for the logistic-accelerated
failure time model. The procedure is implemented in the R-package smcure
([5]). However, for limited sample sizes which are common in practice, such
iterative procedures are characterized by large mean-squared-error (MSE),
convergence problems and instability of the estimators for the incidence
component depending on which variables are included in the latency model (see
for e.g. [24]). This might lead to incorrect conclusions regarding significant
effects.

Here we propose a new second stage estimator based on presmoothing with the
aim of improving upon an initially available estimator, that can for example
be the smcure estimator, for small and moderate sample sizes. The initial
estimator is used to construct a one-dimensional covariate, conditional on
which we compute a nonparametric estimator of the cure probabilities.
Afterwards, the nonparametric estimator is projected on the desired parametric
class (for example logistic). This allows for direct estimation of the
parametric incidence component despite the latent cure status. We focus on the
cure fraction, but once that is estimated, one can also fit a semiparametric
model to the latency component. Compared to the method proposed in [24], this
approach does not restrict us to a one-dimensional covariate and does not
require multidimensional smoothing, which is essential for practical purposes.
Apart from the cure model setting, the idea of constructing a parametric
estimator by nonparametric estimation has been previously proposed in the
context of linear regression, variable selection and functional linear
regression ([7, 1, 12]). The novelty of our method lies in using presmoothing
as a second stage estimator where a preliminary available estimator is used to
reduce the covariate dimension to one. In this way we only need to choose one
bandwidth independently of the number of covariates and still profit from the
advantages of presmoothing: lower MSE and more stable estimators.

The paper is organized as follows. In Sections 2 and 3, we describe the model
and the estimation procedure. In Section 4 we show that the resulting
estimator is consistent and square-root-n convergent with a Gaussian limit
distribution, provided that the initial estimator is consistent. As a
particular case, we focus on the logistic/Cox mixture cure model in Section 5
and illustrate through an extensive simulation study that the proposed
estimator outperforms the smcure estimator by significantly reducing its mean
squared error. In addition, the second step using presmoothing makes the
estimator more stable towards misspecifications in the latency model. Finally,
in Section 6, we apply the method to two medical datasets and show that in
practice it can lead to different conclusions compared to the smcure
estimator.

##  2 The semiparametric mixture cure model

Suppose we are interested in the time $T$ until a certain event happens for a
mixed population of cured ($T=\infty$) and uncured ($T=T_{0}<\infty$)
subjects. Let $B$ be a $0$-$1$ random variable indicating the uncured status:
$B=1$ for susceptible individuals and $B=0$ otherwise. Due to the limited
follow-up period, we cannot actually observe $T$ and $B$. Instead we observe a
finite follow-up time $Y=\min(T,C)$ and a censoring indicator
$\Delta=\mathds{1}_{\\{T\leq C\\}}$, where $C$ denotes the censoring time. As
a result, for all the censored observations, the cure status is unknown. In
the mixture cure model, the survival function of $T$ given two covariate
vectors $X\in\mathbb{R}^{p}$, $Z\in\mathbb{R}^{q}$, is given by

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

where $S_{u}(t|z)=\mathbb{P}(T>t|Z=z,B=1)$ is the survival function of the
susceptibles and $\pi_{0}(x)=\mathbb{P}(B=0|X=x)$ denotes the cure
probability. Using two covariate vectors $X$ and $Z$ for modeling the
incidence and the latency allows the cure probability and the survival of the
uncured to be affected by different variables. However it does not exclude
situations in which the two vectors $X$ and $Z$ are exactly the same or share
some components.

In the context of mixture cure models, the classical survival analysis
assumption of independent censoring, means that $T_{0}\perp(C,X)|Z$ and
$B\perp(C,T_{0},Z)|X$, which imply that

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

(see Lemma 1 in the supplementary material of [24]). As a result, we also have

$\mathbb{P}(T=\infty|X,Z)=\mathbb{P}(T=\infty|X)\quad\text{ and
}\quad\mathbb{P}(T_{0}\leq t|X,Z)=\mathbb{P}(T_{0}\leq t|Z).$ (2)

Among various modeling approaches for the incidence and the latency, the most
common choice is a parametric (logistic) model for the incidence and a
semiparametric (Cox or accelerated failure time) model for the latency ([27,
4, 15, 42, 33, 39]). The popularity of such a choice is due to the simplicity
and ease of interpretation. We focus on this type of model and assume that

$\pi_{0}(x)=1-\phi(\gamma_{0}^{T}x),$

where $\phi:\mathbb{R}\to[0,1]$ is a known function,
$\gamma_{0}\in\mathbb{R}^{p+1}$ and $\gamma^{T}_{0}$ denotes the transpose of
the vector $\gamma_{0}$. Here the first component of $x$ is taken to be equal
to 1 and the first component of $\gamma$ corresponds to the intercept. In
particular, for the logistic model, we have

$\phi(u)=\frac{e^{u}}{1+e^{u}}.$ (3)

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 [3] or use
the test proposed in [23] which is currently developed only for a one-
dimensional covariate.

For the latency, we assume a semiparametric model
$S_{u}(t|z)=S_{u}(t|z;\beta_{0},\Lambda_{0})$ depending on a finite-
dimensional parameter $\beta_{0}\in\mathcal{\mathbb{R}}^{q}$, and an infinite-
dimensional parameter $\Lambda_{0}$. The main examples we keep in mind are the
Cox proportional hazards (PH) model

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

and the accelerated failure time model (AFT)

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

where $\Lambda_{0}$ is the baseline cumulative hazard.

The goal is to estimate the true parameters $\gamma_{0}$, $\beta_{0}$ and
$\Lambda_{0}$ on the basis of $n$ i.i.d. observations
$(Y_{1},\Delta_{1},X_{1},Z_{1}),\dots,(Y_{n},\Delta_{n},X_{n},Z_{n})$. The
general conditions under which the semiparametric mixture cure model is
identifiable, meaning that different parameter values lead to different
distributions of the observed variables $(Y,\Delta,X,Z)$, were derived in [26]
and are the following:

  * (I1)

if $\phi(\gamma,X)=\phi(\tilde{\gamma},X)$ almost surely, then
$\gamma=\tilde{\gamma},$

  * (I2)

the function $S_{u}(\cdot|z)$ has support $[0,\tau(z)]$,

  * (I3)

$\mathbb{P}(C>\tau(Z)|X,Z)>0$ for almost all $X$ and $Z$,

  * (I4)

if, for all $t\geq 0$, we have
$S_{u}(t|Z;\Lambda,\beta)=S_{u}(t|Z;\tilde{\Lambda},\tilde{\beta})$ almost
surely, then $\Lambda=\tilde{\Lambda}$ and $\beta=\tilde{\beta}$,

In the particular case of the logistic-Cox model the conditions become:

  * (I1’)

for all $x$, $0<\phi(\gamma_{0}^{T}x)<1$,

  * (I2’)

the function $S_{u}$ has support $[0;\tau_{0}]$ for some $\tau_{0}<\infty$,

  * (I3’)

$P(C>\tau_{0}|X;Z)>0$ for almost all X and Z,

  * (I4’)

the matrices $Var(X)$ and $Var(Z)$ are positive definite,

(see Proposition 1 and 2 in [26]). Conditions I3 and I3’ are of particular
importance in the context of mixture cure models and essentially tell us that,
in order to correctly identify the cure proportion, we need sufficiently long
follow-up beyond the time when the events occur. In practice, this can be
evaluated based on the plateau of the Kaplan-Meier estimator and the expert
(medical) knowledge.

##  3 The 2-step estimation procedure

Estimation in semiparametric mixture cure models is usually performed via the
expectation maximization algorithm because of the latent cure status. Such
method has been proposed by [31, 35] for the logistic-Cox mixture cure model,
and by [16, 43] for the logistic-AFT model. The procedure is implemented in
the R package smcure ([5]). Despite the simplicity of the method, simultaneous
computation of $\gamma_{0}$, $\beta_{0}$ and $\Lambda_{0}$ through an
iterative procedure leads to several problems for finite, not large sample
sizes which are commonly encountered in practice. This has been previously
reported and illustrated in [24, 4, 13]. The main concerns are the large MSE,
convergence problems and instability of the estimator for the incidence
component depending on which variables are included in the latency model. In
particular, if the latency model is misspecified, even the estimators of the
incidence parameters suffer from induced bias. To alleviate these problems, we
propose the following 2-step estimation procedure that makes use of
presmoothing.

We start with some preliminary estimator $\tilde{\gamma}_{n}$ of $\gamma_{0}$.
This can be any estimator that satisfies the conditions described in Section 4
and in particular for the logistic-Cox or logistic-AFT model we can use the
smcure estimator. We use this preliminary estimator to construct the one-
dimensional index $\tilde{V}=\tilde{\gamma}^{\prime}_{n}X$ estimating
$V=\gamma^{\prime}_{0}X$. Based on this new one-dimensional covariate, we
compute a nonparametric estimator of the cure probability for each subject
defined as follows

$\hat{\pi}_{n}(x)=\prod_{t\in\mathbb{R}}\left(1-\frac{\hat{H}_{1,\tilde{\gamma}_{n}}\left(\mathrm{d}t|\tilde{\gamma}^{T}_{n}x\right)}{\hat{H}_{\tilde{\gamma}_{n}}\left([t,\infty)|\tilde{\gamma}^{T}_{n}x\right)}\right){,}$
(5)

where
$\hat{H}_{\tilde{\gamma}_{n}}([t,\infty)|u)=\hat{H}_{1,\tilde{\gamma}_{n}}([t,\infty)|u)+\hat{H}_{0,\tilde{\gamma}_{n}}([t,\infty)|u)$,
$\hat{H}_{1,\tilde{\gamma}_{n}}(\mathrm{d}t|u)=\hat{H}_{1,\tilde{\gamma}_{n}}((t-\mathrm{d}t,t]|u)$
for small $\mathrm{d}t$ and

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

are estimators of

$H_{l}([t,\infty)|u)=\mathbb{P}\left(Y\geq
t,\Delta=l|{\gamma}^{T}_{0}X=u\right),$ (6)

and $H([t,\infty)|u)=H_{1}([t,\infty)|u)+H_{0}([t,\infty)|u)$. Here $k$ is a
one-dimensional kernel function, $b=b_{n}$ is a bandwidth sequence and
${k}_{{b}}(\cdot)=k(\cdot/b)/b$.

The estimator $\hat{\pi}_{n}(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,\tilde{\gamma}_{n}}\left(\mathrm{d}t|\tilde{\gamma}^{T}_{n}x\right)$
is different from zero only at the observed event times, computation of
$\hat{\pi}_{n}(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^{T}X_{i})^{1-\hat{\pi}_{n}(X_{i})}[1-\phi(\gamma^{T}X_{i})]^{\hat{\pi}_{n}(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}_{n}(X_{i})\right]\log\phi(\gamma^{T}X_{i})+\hat{\pi}_{n}(X_{i})\log\left[1-\phi(\gamma^{T}X_{i})\right]\Big{\\}}{\tau(X_{i})}.$
(7)

We introduce a trimming function $\tau(\cdot)\geq 0$ to avoid regions where
the density function of the index $\gamma^{T}X$, for $\gamma$ in a
neighborhood of $\gamma_{0}$, approaches zero (as done for example in [17]).
We discuss possible choices of $\tau(\cdot)$ in Section 4. Existence and
uniqueness of $\hat{\gamma}_{n}$ hold under the same conditions as for the
maximum likelihood estimator in the binary outcome regression model where
$1-\hat{\pi}_{n}(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. Estimation of the latency component can then be performed
by maximizing the likelihood of the mixture model

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

with respect to $\beta$ and $\Lambda$ for $\gamma=\hat{\gamma}_{n}$. Here
$f_{u}(t|Z;\beta,\Lambda)=-\frac{\mathrm{d}}{\mathrm{d}t}S_{u}(t|Z;\beta,\Lambda)$.
In practice this would mean performing the EM algorithm (as in the smcure
package) only on the latency component, i.e. keeping $\gamma=\hat{\gamma}_{n}$
fixed and updating $\beta$ and $\Lambda$ in each iteration.

We call this a 2-step estimator because it relies on a preliminary estimator
$\tilde{\gamma}_{n}$, which is used to construct the one-dimensional covariate
$\tilde{\gamma}^{T}_{n}X$. In this way, independently of the dimension of $X$,
the kernel estimator requires only one bandwidth parameter. The idea of a
single-index structure is used in several papers to avoid multidimensional
regression (e.g. [34, 17]), but has not previously been exploited in the
context of cure models and presmoothing. A nonparametric estimator for the
cure probability in (5) could be obtained using any nonparametric estimator of
the conditional survival function as done for example in [30, 29] for
estimation of the default probability. Here we use the Beran estimator since
it is easy to compute and exhibits a good behavior.

###  3.1 Rationale behind the new approach

By definition we have

$\pi_{0}(x)=\mathbb{P}(T=\infty\mid
X=x)=\mathbb{P}\left(T=\infty\mid\gamma^{T}_{0}X=\gamma^{T}_{0}x\right).$

Moreover, since from our model it follows that $T\perp
C\mid\gamma^{T}_{0}X,\beta_{0}^{T}Z$, we have

$\displaystyle H_{1}\left(\mathrm{d}t|\gamma^{T}_{0}x,{\beta_{0}^{T}z}\right)$
$\displaystyle=F_{C}\left([t,\infty)|\gamma^{T}_{0}x,{\beta_{0}^{T}z}\right)F_{T}\left(\mathrm{d}t|\gamma^{T}_{0}x,{\beta_{0}^{T}z}\right),$
$\displaystyle H\left([t,\infty)|\gamma^{T}_{0}x,{\beta_{0}^{T}z}\right)$
$\displaystyle=F_{T}\left([t,\infty]|\gamma^{T}_{0}x,{\beta_{0}^{T}z}\right)F_{C}\left([t,\infty)|\gamma^{T}_{0}x,{\beta_{0}^{T}z}\right),$

which yields

$\frac{F_{T}\left(\mathrm{d}t|\gamma^{T}_{0}x,{\beta_{0}^{T}z}\right)}{F_{T}\left([t,\infty]|\gamma^{T}_{0}x,{\beta_{0}^{T}z}\right)}=\frac{H_{1}\left(\mathrm{d}t|\gamma^{T}_{0}x,{\beta_{0}^{T}z}\right)}{H\left([t,\infty)|\gamma^{T}_{0}x,{\beta_{0}^{T}z}\right)}.$

As in [27], we obtain that

$\mathbb{P}\left(T=\infty\mid\gamma^{T}_{0}X=\gamma^{T}_{0}x,{\beta_{0}^{T}Z=\beta_{0}^{T}z}\right)=\prod_{t\in\mathbb{R}}\left\\{1-\frac{H_{1}\left(\mathrm{d}t|\gamma^{T}_{0}x,{\beta_{0}^{T}z}\right)}{H\left([t,\infty)|\gamma^{T}_{0}x,{\beta_{0}^{T}z}\right)}\right\\},$

where $\prod_{t\in\mathbb{R}}$ denotes the product integral. By the first part
of (2), the product integral is also equal to
$P(T=\infty|\gamma_{0}^{T}X=\gamma_{0}^{T}x)$. Similarly, under the slightly
stronger assumption that $T\perp C\mid\gamma^{T}_{0}X$, we obtain

$\mathbb{P}\left(T=\infty\mid\gamma^{T}_{0}X=\gamma^{T}_{0}x\right)=\prod_{t\in\mathbb{R}}\left\\{1-\frac{H_{1}\left(\mathrm{d}t|\gamma^{T}_{0}x\right)}{H\left([t,\infty)|\gamma^{T}_{0}x\right)}\right\\},$

This justifies the definition of our estimator in (5). This assumption is
satisfied if $C\perp(X,Z)\mid\gamma_{0}^{T}X$. We restrict ourselves to this
case for simplicity in order to have conditioning on only one index. However,
the method can be used in general conditioning on both $\gamma_{0}^{T}X$ and
$\beta_{0}^{T}Z$. We illustrate this through one of the simulation settings in
Section 5.

##  4 Asymptotic results

In this section we focus on models where the survival function of the
susceptibles $S_{u}$ has fixed support $[0,\tau_{0}]$ such that

$\inf_{\gamma\in G}\mathbb{P}(C>\tau_{0}\mid\gamma^{T}X)>0\,\text{ almost
surely}.$ (8)

This is the case when the latency follows a Cox regression model and the
identifiability assumption (I3’) is satisfied. Moreover, we assume that
$T\perp C\mid\gamma^{T}_{0}X$ but the results can be generalized as mentioned
in the previous section.

Let us sketch the arguments we will use to obtain asymptotic properties of
$\hat{\gamma}_{n}$. Note that $\hat{\gamma}_{n}$ is the maximizer of the
Bernoulli-type log-likelihood in (7). Hence, the main issue is dealing with
the nonparametric estimator $\hat{\pi}_{n}$, which replaces the latent binary
outcome. By definition we have

$\hat{\pi}_{n}(x)=\hat{g}_{\tilde{\gamma}_{n}}(\tilde{\gamma}^{T}_{n}x),$

where

$\hat{g}_{\gamma}(u)=1-\hat{F}_{n}(\tau_{0}\mid\gamma^{T}X=u)=\prod_{t\in\mathbb{R}}\left(1-\frac{\hat{H}_{1,{\gamma}}(\mathrm{d}t|u)}{\hat{H}_{{\gamma}}([t,\infty)|u)}\right).$
(9)

Note that the product integral actually over $t\in\mathbb{R}$ is the same as
over $t\leq\tau_{0}$ because $\hat{H}_{1,{\gamma}}(\mathrm{d}t|u)=0$ for
$t>Y_{(m)}$, where $Y_{(m)}\leq\tau_{0}$ is the last observed event time. We
can also write

$\pi_{0}(x)=g_{\gamma_{0}}(\gamma^{T}_{0}x),$

with

${g}_{\gamma}(u)=1-{F}_{T}(\tau_{0}\mid\gamma^{T}X=u)=\prod_{t\in\mathbb{R}}\left(1-\frac{{H}_{1}(dt|u)}{H([t,\infty)|u)}\right).$
(10)

Hence

$\displaystyle\hat{\pi}_{n}(x)-\pi_{0}(x)$
$\displaystyle=\hat{g}_{\tilde{\gamma}_{n}}\left(\tilde{\gamma}^{T}_{n}x\right)-g_{\gamma_{0}}\left(\gamma^{T}_{0}x\right)$
(11)
$\displaystyle=\left\\{\hat{g}_{{\gamma}_{0}}\left({\gamma}^{T}_{0}x\right)-g_{\gamma_{0}}\left(\gamma^{T}_{0}x\right)\right\\}+\left\\{\hat{g}_{\tilde{\gamma}_{n}}\left(\tilde{\gamma}^{T}_{n}x\right)-\hat{g}_{\gamma_{0}}\left(\gamma^{T}_{0}x\right)\right\\}$
$\displaystyle=\left\\{\hat{F}_{n}\left(\tau_{0}\mid\gamma^{T}_{0}X=\gamma^{T}_{0}x\right)-{F}_{T}\left(\tau_{0}\mid\gamma^{T}_{0}X=\gamma^{T}_{0}x\right)\right\\}$
$\displaystyle\qquad+\left\\{\hat{g}_{\tilde{\gamma}_{n}}\left(\tilde{\gamma}^{T}_{n}x\right)-\hat{g}_{\gamma_{0}}\left(\gamma^{T}_{0}x\right)\right\\}.$

The first term on the right hand side of the equation, can be dealt with as
usual being the difference between $F_{n}$ and $F$ conditionally on a one
dimensional covariate $\gamma^{T}_{0}X$. The second term results from using
$\tilde{\gamma}_{n}$ instead of $\gamma_{0}$ when constructing the one-
dimensional covariate $\tilde{\gamma}^{T}_{n}X$. The behaviour of this term
depends on the properties of the preliminary estimator $\tilde{\gamma}_{n}$.
We first formulate the results for a general prelimary estimator
$\tilde{\gamma}_{n}$ and a general parametric function $\phi$. We then show
that, for the logistic-Cox model, the maximum likelihood estimator satisfies
the required conditions.

The following assumptions are needed for consistency of $\hat{\gamma}_{n}$.

  * (C1)

The preliminary estimator is consistent, i.e.
$\tilde{\gamma}_{n}-\gamma_{0}=o_{P}(1)$.

  * (C2)

The parameter $\gamma_{0}$ lies in the interior of a compact set
$G\subset\mathbb{R}^{p}$.

  * (C3)

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

$\left|\phi(\gamma_{1}^{T}x)-\phi(\gamma_{2}^{T}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 norm and
$\mathcal{X}\subset\mathbb{R}^{p}$ is the support of $X$.

  * (C4)

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

  * (C5)

For any $\gamma\in G$, the support $\mathcal{X}_{\gamma}$ of
$\gamma^{\prime}\mathcal{X}$ is a bounded convex subset of $\mathbb{R}$. The
density $f_{X_{\gamma}}(\cdot)$ of $X_{\gamma}$ is twice differentiable with a
bounded second derivative.

  * (C6)

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

  * (C7)

The kernel $k$ is a twice continuously differentiable, symmetric probability
density function with compact support.

  * (C8)

(i) The functions $H([0,t]|u)$, $H_{1}([0,t]|u)$ defined in (6) are twice
differentiable with respect to $u$, with uniformly bounded derivatives for all
$t\leq\tau_{0}$, $u\in\mathcal{X}_{\gamma_{0}}$. 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}]$,
$u\in\mathcal{X}_{\gamma_{0}}$,

$\begin{split}\left|H_{c}(t|u)-H_{c}(s|u)\right|&\leq\left|L_{1}(t)-L_{1}(s)\right|,\\\
\left|H_{1c}(t|u)-H_{1c}(s|u)\right|&\leq\left|L_{1}(t)-L_{1}(s)\right|,\\\
\left|\frac{\partial H_{c}(t|u)}{\partial u}-\frac{\partial
H_{c}(s|u)}{\partial u}\right|&\leq\left|L_{2}(t)-L_{2}(s)\right|,\\\
\left|\frac{\partial H_{1c}(t|u)}{\partial u}-\frac{\partial
H_{1c}(s|u)}{\partial
u}\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 number of jump points for the distribution function $F_{C}(t|u)$ of
the censoring times given the index $\gamma_{0}^{T}X=u$, are finite and the
same for all $u$. The partial derivative of $F_{C}(t|u)$ with respect to $u$
exists and is uniformly bounded for all $t\leq\tau_{0}$,
$u\in\mathcal{X}_{\gamma_{0}}$. Moreover, the partial derivative with respect
to $u$ of $F_{T}(t|u)$ (distribution function of the survival times $T$ given
$\gamma_{0}^{T}X=u$) exists and is uniformly bounded for all $t\leq\tau_{0}$,
$u\in\mathcal{X}_{\gamma_{0}}$.

  * (C9)

The function $(x,\gamma)\mapsto g_{\gamma}(\gamma^{T}x)$ is continuously
differentiable with respect to $\gamma$ and the vector
$\nabla_{\gamma}g_{\gamma}(\gamma^{T}x)$ is continuous with respect to
$(x,\gamma)$.

Assumptions (C2)-(C4), (C6)-(C8) are standard assumptions (see for example
[27, 24, 38]. Assumptions (C5) and (C9) are needed because we compute the
nonparametric estimator using the index $\tilde{\gamma}_{n}^{T}X$ instead of
$\gamma_{0}^{T}X$. Such assumptions appear for example in [17].

A possible choice of the trimming function $\tau(\cdot)$ in (7) could be
$\tau(x)=\mathds{1}_{\tilde{\mathcal{X}}}(x)$ if we know a set
$\tilde{\mathcal{X}}$ such as

$\inf_{\gamma\in
G}\inf_{x\in\tilde{\mathcal{X}}}f_{X_{\gamma}}(\gamma^{T}x)=c>0.$ (12)

Otherwise, as shown in [17], one can take
$\tau(x)=\mathds{1}_{\\{f_{X_{\gamma_{0}}}(\gamma_{0}^{T}x)\geq c\\}}$ for
some $c>0$, which is asymptotically equivalent to the previous proposal. In
practice, we can use
$\tau(x)=\mathds{1}_{\\{\hat{f}_{X_{\tilde{\gamma}}}(\tilde{\gamma}_{n}^{T}x)\geq
c\\}}$ based on the preliminary estimator $\tilde{\gamma}_{n}$.

######  Theorem 1.

Assume that conditions (C1)-(C9) are satisfied. Then

$\hat{\gamma}_{n}-\gamma_{0}=o_{P}(1).$

In order to obtain asymptotic normality of $\hat{\gamma}_{n}$ at rate
$\sqrt{n}$, we need the following additional assumptions.

  * (N1)

For each $x\in\mathcal{X}$, the function $\gamma\mapsto\phi(\gamma^{T}x)$ is
twice continuously differentiable with uniformly bounded derivatives in
$G\times\mathcal{X}$.

  * (N2)

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

  * (N3)

The preliminary estimator $\tilde{\gamma}_{n}$ is $\sqrt{n}$ consistent and
such that there exists a function $\zeta$ such that

$\tilde{\gamma}_{n}-\gamma_{0}=\frac{1}{n}\sum_{i=1}^{n}\zeta(T_{i},\Delta_{i},X_{i},Z_{i})+R_{n},$

with $\|R_{n}\|=o_{P}(n^{-1/2})$ and $\mathbb{E}[\zeta(T,\Delta,X,Z)]=0$.

Again (N1)-(N2) are standard assumptions, while (N3) arises from the use of
the index $\tilde{\gamma}_{n}^{T}X$ instead of $\gamma_{0}^{T}X$. As a result
the asymptotic variance of $\hat{\gamma}_{n}$ will also depend on the
asymptotic variance of the preliminary estimator $\tilde{\gamma}_{n}$.

######  Theorem 2.

Assume that conditions (C1)-(C9), (N1)-(N3) are satisfied. Then

$\sqrt{n}(\hat{\gamma}_{n}-\gamma_{0})\xrightarrow{d}N(0,\Sigma_{\gamma}),$

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

Given the complicated form of the covariance matrix $\Sigma_{\gamma}$, we
suggest using a bootstrap procedure for estimating the standard errors as also
done for the maximum likelihood estimator of a semi-parametric mixture cure
model.

If we consider the particular case of a logistic-Cox mixture cure model and
take the maximum likelihood estimator as a preliminary estimator
$\tilde{\gamma}_{n}$, then assumptions (C2)-(C4), (C9), (N1)-(N2) are
obviously satisfied for the logistic model. Morevoer (C1) and (N2) are
satisfied if the cumulative baseline function $\Lambda_{0}$ is strictly
increasing and continuously differentiable under the condition

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

(see Theorem 2 and Theorem 3 in [20]). Then, from Theorems 1 and 2 it follows
that the 2-step estimator is also consistent and $\sqrt{n}$-convergent. If we
continue estimating the latency sub-model using this estimator of
$\gamma_{0}$, then the resulting estimator of $\beta_{0}$ and $\Lambda_{0}$
have the desired asymptotic behavior as in Theorems 2 and 4 in [24]. The proof
remains the same given that they only use consistency and the asymptotic
i.i.d. expression of the estimator as in assumption (N3).

##  5 Simulation study

In this section we investigate the finite-sample behaviour of the 2-step
approach in the logistic/Cox mixture cure model and compare it with the
maximum likelihood estimator implemented in the R package smcure. We use the
smcure estimator as preliminary estimator $\tilde{\gamma}_{n}$ for the new
method.

We make some standard and common choices when computing the nonparametric
estimator in (5). The kernel function $k$ is taken to be the Epanechnikov
kernel $k(u)=(3/4)(1-u^{2})\mathds{1}_{\\{|u|\leq 1\\}}$. Using the
preliminary estimator $\tilde{\gamma}_{n}$, we compute the smoothing bandwidth
by cross-validation as implemented in the R package np for kernel estimators
of conditional distribution functions, in our case for the estimation of
$H=H_{0}+H_{1}$ given $\tilde{\gamma}^{T}X$. In addition, we restrict
ourselves 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 (5) 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 estimating
the conditional distribution $H(t|\tilde{\gamma}^{T}x)$ for $t\leq Y_{(m)}$.

We could use a trimming function

$\tau(x)=\mathds{1}_{\\{\hat{f}_{X_{\tilde{\gamma}}}(\tilde{\gamma}_{n}^{T}x)\geq
c\\}},$

for some small value of $c$ as proposed in Section 4. However, we observe that
in practice this does not affect the results since $c$ can be chosen as
arbitrarily small. Hence, we do not do any trimming so that we do not have to
worry about the choice of the trimming constant. The trimming is mainly
introduced for the asymptotic study in order to avoid the assumption that the
density of the index is bounded from below by a positive constant.

We consider four different models and for each of them, three scenarios,
covering a wide range of settings with different number and choice of
covariates (continuous and discrete), different cure and censoring rate, and
different censoring mechanisms (independent of covariates, depending on the
same index as the incidence model, depending on both indexes of the incidence
and latency). The models are as follows.

Model 1. Both incidence and latency depend on two independent continuous
covariates $X_{1}=Z_{1}\sim N(0,1)$ and $X_{2}=Z_{2}\sim\text{Unif}(-1,1)$. We
generate the cure status $B$ as a Bernoulli random variable with success
probability $\phi(\gamma^{T}X)$ where $\phi$ is the logistic function in (3)
and $\gamma=(\gamma_{0},1.5,1.5)$. 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\left(\beta^{T}z\right)\right),$

and are truncated at $\tau_{0}=15$ for $\rho=0.75$, $\mu=1.5$ and
$\beta=(0.5,0.3)$. 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=17$.

Model 2. Both incidence and latency depend on three independent covariates
$X_{1}=Z_{1}\sim N(0,1)$, $X_{2}=Z_{2}\sim\text{Bernoulli}(0.3)$ and
$X_{3}=Z_{3}\sim\text{Bernoulli}(0.7)$. The cure status and the survival times
for the uncured observations are generated as in Model 1 for
$\gamma=(\gamma_{0},-1,1,-0.3)$, $\beta=(-0.8,1.5,-0.5)$, $\rho=0.75$,
$\mu=1.5$, and $\tau_{0}=7$. The censoring times are generated according to a
Weibull proportional hazards model

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

for various choices of $\lambda_{C}$ and are truncated at $\tau=9$.

Model 3. For the incidence we consider four independent covariates: $X_{1}\sim
N(0,1)$, $X_{2}\sim\text{Unif}(-1,1)$, $X_{3}$ and $X_{4}$ are Bernoulli
random variables with parameters $0.4$ and $0.6$ respectively. The latency
depends on three covariates: $Z_{1}\sim N(0,1)$, $Z_{2}=X_{2}$ and
$Z_{3}=X_{4}$. The cure status and the survival times for the uncured
observations are generated as in Model 1 for
$\gamma=(\gamma_{0},-0.3,0.8,0.5,-1)$, $\rho=0.75$, $\mu=1.5$,
$\beta=(0.1,0.4,-0.2)$ and $\tau_{0}=10$. The censoring times are generated
generated according to a Weibull proportional hazards model

$S_{C}(t|x)=\exp\left(-\lambda_{C}\mu
t^{\rho}\exp\left(0.4\gamma^{T}x+0.5\beta^{T}z\right)\right),$

for various choices of $\lambda_{C}$ and are truncated at $\tau=12$.

Model 4. For the incidence we consider five independent covariates: $X_{1}\sim
N(0,1)$, $X_{2}\sim\text{Unif}(-1,1)$, $X_{3}$ is Binomial with parameters $2$
and $0.5$, $X_{4}$ and $X_{5}$ are Bernoulli random variables with parameters
$0.4$ and $0.6$ respectively. The latency depends on three covariates:
$Z_{1}\sim N(0,1)$, $Z_{2}=X_{3}$ and $Z_{3}=X_{4}$. The cure status and the
survival times for the uncured observations are generated as in Model 1 for
$\gamma=(\gamma_{0},-0.8,0.3,-0.4,0.5,0.6)$, $\rho=0.75$, $\mu=1.5$,
$\beta=(0.2,-0.5,0.3)$ and $\tau_{0}=7$. The censoring times are independent
from $X$, $Z$ and $T$. They are generated from the exponential distribution
with parameter $\lambda_{C}$ and are truncated at $\tau=9$.

For the four models we choose the values of the unspecified parameters
$\gamma_{0}$ and $\lambda_{C}$ in such a way that the cure rate is around
$20\%$, $40\%$ or $60\%$ and the difference between the cure and the censoring
rate is around $5\%$, $10\%$ or $15\%$. The specification of the parameters
and the corresponding censoring and cure rates are given Table 1. 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 (13) 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. Model 3 illustrates the behavior of the method when the
censoring times depends on both the indexes of the incidence and latency
models (which was assumed for simplicity in the theoretical study).

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

Model Scenario $\gamma_{0}$ $\lambda_{C}$ Cens. rate Cure rate $1$ $2$ $0.4$
$36\%$ $20\%$ 1 $2$ $0.6$ $0.4$ $50\%$ $40\%$ $3$ $-0.5$ $0.3$ $63\%$ $58\%$
$1$ $1.6$ $1/35$ $30\%$ $20\%$ 2 $2$ $0.4$ $1/20$ $45\%$ $40\%$ $3$ $-0.6$ $1$
$75\%$ $60\%$ $1$ $2$ $1/9$ $35\%$ $20\%$ 3 $2$ $0.9$ $1/7$ $50\%$ $40\%$ $3$
$-0.1$ $1/7$ $65\%$ $60\%$ $1$ $1.5$ $0.6$ $25\%$ $20\%$ 4 $2$ $0.3$ $0.3$
$55\%$ $40\%$ $3$ $-0.6$ $0.4$ $70\%$ $60\%$

We consider samples of size $n=200$ and $n=400$ since we aim to provide a
method that improves upon the maximum likelihood estimator for small and
moderate sample size. For each configuration, $1000$ datasets were generated
and the estimators of $\beta$ and $\gamma$ were computed through smcure and
the proposed 2-step approach. We report the bias, variance and mean squared
error (MSE) of the estimators, computed over the iterations for which the
smcure procedure converges, in Tables 2-4. In some scenarios, mainly
corresponding to the ones with $15\%$ additional censoring compared to the
cure rate and smaller sample size, the iterative procedure of the EM algorithm
in smcure does not converge. The most problematic setting in this regard is
Model 2 scenario 3, for which $73/1000$ iterations do not converge for $n=200$
and $36/1000$ for $n=400$. The boxplots of the estimators for both methods in
these non-convergent iterations are shown in Figure 1. In the other settings,
only around $2\%$ or less of the iterations do not converge.

Table 2: Bias, variance and MSE of $\hat{\gamma}$ and $\hat{\beta}$ for smcure
(second rows) and the 2-step approach (first rows) in Models 1 and 2.

Scenario 1 Scenario 2 Scenario 3 Mod. n Par. Bias Var. MSE Bias Var. MSE Bias
Var. MSE 1 $200$ $\gamma_{1}$ $-0.007$ $0.246$ $0.246$ $-0.013$ $0.062$
$0.062$ $0.088$ $0.049$ $0.057$ $0.243$ $0.534$ $0.593$ $0.058$ $0.092$
$0.095$ $0.096$ $0.053$ $0.062$ $\gamma_{2}$ $-0.009$ $0.217$ $0.217$ $-0.035$
$0.116$ $0.117$ $-0.043$ $0.089$ $0.091$ $0.201$ $0.302$ $0.343$ $0.096$
$0.146$ $0.155$ $0.049$ $0.099$ $0.102$ $\gamma_{3}$ $-0.050$ $0.406$ $0.408$
$-0.028$ $0.205$ $0.206$ $-0.031$ $0.152$ $0.153$ $0.189$ $0.717$ $0.753$
$0.108$ $0.242$ $0.253$ $0.063$ $0.174$ $0.178$ $\beta_{1}$ $-0.006$ $0.015$
$0.015$ $0.005$ $0.024$ $0.024$ $0.015$ $0.034$ $0.034$ $-0.009$ $0.016$
$0.016$ $-0.002$ $0.025$ $0.025$ $0.009$ $0.034$ $0.034$ $\beta_{2}$ $0.008$
$0.035$ $0.035$ $0.008$ $0.050$ $0.050$ $0.012$ $0.071$ $0.072$ $0.003$
$0.036$ $0.036$ $-0.002$ $0.052$ $0.052$ $0.002$ $0.073$ $0.073$ $400$
$\gamma_{1}$ $-0.004$ $0.109$ $0.109$ $-0.002$ $0.036$ $0.036$ $0.000$ $0.026$
$0.026$ $0.129$ $0.161$ $0.178$ $0.040$ $0.043$ $0.045$ $0.004$ $0.027$
$0.027$ $\gamma_{2}$ $-0.015$ $0.094$ $0.095$ $-0.023$ $0.059$ $0.059$
$-0.012$ $0.047$ $0.048$ $0.101$ $0.108$ $0.118$ $0.059$ $0.059$ $0.063$
$0.057$ $0.052$ $0.055$ $\gamma_{3}$ $-0.008$ $0.201$ $0.201$ $-0.019$ $0.110$
$0.110$ $-0.010$ $0.084$ $0.084$ $0.116$ $0.247$ $0.261$ $0.065$ $0.115$
$0.119$ $0.060$ $0.088$ $0.091$ $\beta_{1}$ $0.005$ $0.008$ $0.008$ $0.013$
$0.012$ $0.012$ $0.011$ $0.017$ $0.017$ $0.003$ $0.008$ $0.008$ $0.009$
$0.012$ $0.012$ $0.006$ $0.017$ $0.017$ $\beta_{2}$ $0.007$ $0.019$ $0.019$
$0.012$ $0.026$ $0.027$ $0.015$ $0.034$ $0.034$ $0.004$ $0.019$ $0.019$
$0.005$ $0.027$ $0.027$ $0.008$ $0.035$ $0.035$ $2$ $200$ $\gamma_{1}$ $0.049$
$0.248$ $0.250$ $0.011$ $0.114$ $0.114$ $0.000$ $0.244$ $0.244$ $0.103$
$0.268$ $0.278$ $0.021$ $0.121$ $0.122$ $0.097$ $0.389$ $0.398$ $\gamma_{2}$
$0.017$ $0.073$ $0.074$ $0.010$ $0.047$ $0.047$ $0.090$ $0.091$ $0.099$
$-0.081$ $0.088$ $0.094$ $-0.040$ $0.051$ $0.052$ $-0.154$ $0.169$ $0.193$
$\gamma_{3}$ $0.025$ $0.947$ $0.948$ $-0.009$ $0.171$ $0.171$ $-0.102$ $0.299$
$0.310$ $0.167$ $1.895$ $1.923$ $0.035$ $0.177$ $0.178$ $0.084$ $0.444$
$0.451$ $\gamma_{4}$ $0.007$ $0.273$ $0.273$ $0.005$ $0.150$ $0.150$ $0.031$
$0.270$ $0.271$ $-0.021$ $0.303$ $0.303$ $-0.011$ $0.162$ $0.162$ $-0.040$
$0.459$ $0.461$ $\beta_{1}$ $-0.003$ $0.014$ $0.014$ $-0.010$ $0.021$ $0.021$
$-0.074$ $0.069$ $0.075$ $0.002$ $0.014$ $0.014$ $-0.007$ $0.021$ $0.021$
$-0.016$ $0.085$ $0.086$ $\beta_{2}$ $0.029$ $0.050$ $0.051$ $0.043$ $0.069$
$0.071$ $0.099$ $0.236$ $0.245$ $0.027$ $0.051$ $0.051$ $0.042$ $0.069$
$0.071$ $0.074$ $0.269$ $0.274$ $\beta_{3}$ $-0.009$ $0.044$ $0.044$ $-0.014$
$0.052$ $0.052$ $-0.034$ $0.220$ $0.221$ $-0.008$ $0.045$ $0.045$ $-0.013$
$0.052$ $0.052$ $-0.028$ $0.252$ $0.253$ $400$ $\gamma_{1}$ $0.029$ $0.111$
$0.112$ $0.010$ $0.055$ $0.055$ $0.017$ $0.107$ $0.107$ $0.055$ $0.111$
$0.114$ $0.013$ $0.056$ $0.056$ $0.058$ $0.129$ $0.133$ $\gamma_{2}$ $0.016$
$0.040$ $0.040$ $0.006$ $0.022$ $0.022$ $0.061$ $0.046$ $0.049$ $-0.043$
$0.039$ $0.041$ $-0.025$ $0.022$ $0.023$ $-0.082$ $0.064$ $0.071$ $\gamma_{3}$
$-0.030$ $0.155$ $0.156$ $-0.017$ $0.079$ $0.080$ $-0.111$ $0.142$ $0.154$
$0.030$ $0.160$ $0.161$ $0.012$ $0.084$ $0.084$ $0.006$ $0.163$ $0.163$
$\gamma_{4}$ $0.010$ $0.114$ $0.114$ $0.015$ $0.066$ $0.067$ $0.030$ $0.128$
$0.129$ $-0.006$ $0.119$ $0.119$ $0.007$ $0.068$ $0.068$ $-0.003$ $0.167$
$0.167$ $\beta_{1}$ $0.004$ $0.007$ $0.007$ $-0.001$ $0.008$ $0.008$ $-0.037$
$0.030$ $0.032$ $0.007$ $0.007$ $0.007$ $0.001$ $0.008$ $0.008$ $0.004$
$0.036$ $0.036$ $\beta_{2}$ $0.010$ $0.024$ $0.024$ $0.018$ $0.029$ $0.029$
$0.063$ $0.104$ $0.108$ $0.009$ $0.024$ $0.024$ $0.017$ $0.029$ $0.029$
$0.043$ $0.112$ $0.114$ $\beta_{3}$ $0.001$ $0.021$ $0.021$ $-0.005$ $0.026$
$0.026$ $-0.015$ $0.100$ $0.100$ $0.001$ $0.022$ $0.022$ $-0.005$ $0.026$
$0.026$ $-0.012$ $0.112$ $0.112$

Table 3: Bias, variance and MSE of $\hat{\gamma}$ and $\hat{\beta}$ for smcure
(second rows) and the 2-step approach (first rows) in Models 3-4 for $n=200$.

Scenario 1 Scenario 2 Scenario 3 Mod. n Par. Bias Var. MSE Bias Var. MSE Bias
Var. MSE 3 $200$ $\gamma_{1}$ $-0.056$ $0.287$ $0.290$ $-0.016$ $0.123$
$0.123$ $-0.010$ $0.097$ $0.097$ $0.431$ $4.843$ $5.028$ $0.068$ $0.144$
$0.149$ $0.011$ $0.098$ $0.099$ $\gamma_{2}$ $0.083$ $0.041$ $0.048$ $0.043$
$0.033$ $0.035$ $0.022$ $0.030$ $0.031$ $-0.029$ $0.071$ $0.072$ $-0.014$
$0.040$ $0.040$ $-0.012$ $0.032$ $0.032$ $\gamma_{3}$ $-0.049$ $0.177$ $0.179$
$-0.021$ $0.105$ $0.106$ $0.004$ $0.093$ $0.093$ $0.105$ $0.268$ $0.279$
$0.049$ $0.120$ $0.123$ $0.043$ $0.097$ $0.099$ $\gamma_{4}$ $-0.080$ $0.216$
$0.222$ $-0.041$ $0.131$ $0.133$ $-0.047$ $0.122$ $0.124$ $0.159$ $1.233$
$1.259$ $0.046$ $0.151$ $0.154$ $0.010$ $0.131$ $0.131$ $\gamma_{5}$ $0.123$
$0.325$ $0.340$ $0.051$ $0.158$ $0.161$ $0.021$ $0.136$ $0.137$ $-0.359$
$4.817$ $4.946$ $-0.062$ $0.183$ $0.187$ $-0.036$ $0.140$ $0.141$ $\beta_{1}$
$0.003$ $0.011$ $0.011$ $0.003$ $0.016$ $0.016$ $0.016$ $0.025$ $0.025$
$0.003$ $0.012$ $0.012$ $0.003$ $0.017$ $0.017$ $0.016$ $0.025$ $0.025$
$\beta_{2}$ $0.022$ $0.034$ $0.034$ $0.017$ $0.045$ $0.046$ $0.028$ $0.075$
$0.076$ $0.010$ $0.036$ $0.036$ $0.008$ $0.048$ $0.048$ $0.023$ $0.077$
$0.077$ $\beta_{3}$ $-0.048$ $0.044$ $0.046$ $-0.042$ $0.060$ $0.061$ $-0.030$
$0.093$ $0.094$ $-0.032$ $0.048$ $0.049$ $-0.027$ $0.063$ $0.064$ $-0.023$
$0.095$ $0.095$ $4$ $200$ $\gamma_{1}$ $0.086$ $0.237$ $0.245$ $0.037$ $0.238$
$0.239$ $0.010$ $0.229$ $0.229$ $0.110$ $0.273$ $0.285$ $0.077$ $0.311$
$0.317$ $0.013$ $0.310$ $0.310$ $\gamma_{2}$ $0.010$ $0.058$ $0.058$ $0.042$
$0.065$ $0.067$ $0.006$ $0.065$ $0.069$ $-0.069$ $0.064$ $0.069$ $-0.091$
$0.077$ $0.085$ $-0.073$ $0.068$ $0.073$ $\gamma_{3}$ $-0.025$ $0.131$ $0.132$
$-0.018$ $0.122$ $0.122$ $-0.021$ $0.122$ $0.122$ $0.000$ $0.153$ $0.153$
$0.033$ $0.162$ $0.164$ $0.027$ $0.155$ $0.156$ $\gamma_{4}$ $0.008$ $0.095$
$0.095$ $0.007$ $0.108$ $0.109$ $-0.007$ $0.110$ $0.110$ $-0.025$ $0.110$
$0.111$ $-0.049$ $0.138$ $0.140$ $-0.075$ $0.156$ $0.161$ $\gamma_{5}$
$-0.026$ $0.212$ $0.212$ $-0.055$ $0.199$ $0.202$ $-0.038$ $0.193$ $0.194$
$0.019$ $0.245$ $0.245$ $0.020$ $0.287$ $0.287$ $0.045$ $0.268$ $0.270$
$\gamma_{6}$ $-0.024$ $0.176$ $0.177$ $-0.026$ $0.172$ $0.172$ $-0.038$
$0.163$ $0.165$ $0.031$ $0.203$ $0.204$ $0.068$ $0.213$ $0.217$ $0.063$
$0.209$ $0.213$ $\beta_{1}$ $0.004$ $0.008$ $0.008$ $0.000$ $0.019$ $0.019$
$-0.001$ $0.035$ $0.035$ $0.004$ $0.008$ $0.008$ $0.003$ $0.019$ $0.019$
$0.001$ $0.037$ $0.037$ $\beta_{2}$ $-0.011$ $0.019$ $0.019$ $-0.001$ $0.049$
$0.049$ $0.010$ $0.095$ $0.096$ $-0.011$ $0.019$ $0.019$ $0.004$ $0.053$
$0.053$ $0.021$ $0.103$ $0.104$ $\beta_{3}$ $0.000$ $0.033$ $0.033$ $0.014$
$0.083$ $0.083$ $0.021$ $0.144$ $0.144$ $-0.001$ $0.033$ $0.033$ $0.006$
$0.088$ $0.088$ $0.007$ $0.155$ $0.155$

Figure 1: Boxplots of estimates of the four $\gamma$ parameters in Model 2,
scenario 3 for the iterations in which smcure does not converge. The
horizontal lines correspond to the true values of the parameters. Left panel:
for $n=200$. Right panel: for $n=400$. Table 4: Bias, variance and MSE of
$\hat{\gamma}$ and $\hat{\beta}$ for smcure (second rows) and the 2-step
approach (first rows) in Models 3-4 for $n=400$.

Scenario 1 Scenario 2 Scenario 3 Mod. n Par. Bias Var. MSE Bias Var. MSE Bias
Var. MSE 3 $400$ $\gamma_{1}$ $-0.045$ $0.149$ $0.151$ $-0.009$ $0.069$
$0.070$ $-0.020$ $0.046$ $0.047$ $0.124$ $0.213$ $0.228$ $0.043$ $0.071$
$0.073$ $-0.003$ $0.045$ $0.045$ $\gamma_{2}$ $0.063$ $0.021$ $0.025$ $0.034$
$0.016$ $0.017$ $0.023$ $0.015$ $0.015$ $-0.018$ $0.030$ $0.030$ $-0.008$
$0.018$ $0.018$ $-0.003$ $0.016$ $0.016$ $\gamma_{3}$ $-0.041$ $0.098$ $0.099$
$-0.025$ $0.058$ $0.059$ $0.001$ $0.050$ $0.050$ $0.052$ $0.116$ $0.119$
$0.019$ $0.061$ $0.062$ $0.026$ $0.052$ $0.052$ $\gamma_{4}$ $-0.088$ $0.111$
$0.119$ $-0.052$ $0.069$ $0.071$ $-0.020$ $0.059$ $0.059$ $0.029$ $0.149$
$0.149$ $0.007$ $0.076$ $0.076$ $0.019$ $0.060$ $0.060$ $\gamma_{5}$ $0.103$
$0.165$ $0.175$ $0.047$ $0.083$ $0.085$ $0.033$ $0.068$ $0.069$ $-0.079$
$0.226$ $0.232$ $-0.030$ $0.085$ $0.086$ $-0.013$ $0.066$ $0.067$ $\beta_{1}$
$-0.001$ $0.005$ $0.005$ $0.001$ $0.007$ $0.007$ $0.003$ $0.010$ $0.010$
$-0.001$ $0.005$ $0.005$ $0.001$ $0.007$ $0.007$ $0.003$ $0.010$ $0.010$
$\beta_{2}$ $0.007$ $0.017$ $0.017$ $0.006$ $0.023$ $0.023$ $0.002$ $0.034$
$0.034$ $0.000$ $0.018$ $0.018$ $0.001$ $0.024$ $0.024$ $-0.001$ $0.035$
$0.035$ $\beta_{3}$ $-0.017$ $0.021$ $0.021$ $-0.012$ $0.028$ $0.028$ $-0.007$
$0.042$ $0.042$ $-0.002$ $0.022$ $0.022$ $-0.003$ $0.029$ $0.029$ $-0.002$
$0.043$ $0.043$ $4$ $400$ $\gamma_{1}$ $0.042$ $0.115$ $0.116$ $0.016$ $0.105$
$0.105$ $0.003$ $0.118$ $0.118$ $0.056$ $0.124$ $0.127$ $0.036$ $0.119$
$0.121$ $0.007$ $0.125$ $0.125$ $\gamma_{2}$ $0.025$ $0.027$ $0.027$ $0.045$
$0.036$ $0.038$ $0.047$ $0.036$ $0.039$ $-0.029$ $0.028$ $0.029$ $-0.040$
$0.032$ $0.034$ $-0.038$ $0.032$ $0.034$ $\gamma_{3}$ $0.002$ $0.064$ $0.064$
$-0.022$ $0.062$ $0.062$ $-0.018$ $0.063$ $0.063$ $0.024$ $0.072$ $0.072$
$0.009$ $0.073$ $0.073$ $0.012$ $0.074$ $0.074$ $\gamma_{4}$ $0.010$ $0.045$
$0.045$ $0.016$ $0.054$ $0.054$ $0.018$ $0.061$ $0.061$ $-0.014$ $0.049$
$0.049$ $-0.019$ $0.058$ $0.059$ $-0.018$ $0.069$ $0.070$ $\gamma_{5}$ $0.000$
$0.101$ $0.101$ $-0.019$ $0.097$ $0.097$ $-0.012$ $0.107$ $0.108$ $0.034$
$0.111$ $0.112$ $0.033$ $0.111$ $0.112$ $0.038$ $0.119$ $0.120$ $\gamma_{6}$
$-0.039$ $0.084$ $0.086$ $-0.045$ $0.081$ $0.084$ $-0.041$ $0.090$ $0.092$
$-0.002$ $0.094$ $0.094$ $0.015$ $0.091$ $0.092$ $0.018$ $0.098$ $0.098$
$\beta_{1}$ $0.000$ $0.004$ $0.004$ $0.003$ $0.009$ $0.009$ $-0.002$ $0.015$
$0.015$ $0.000$ $0.004$ $0.004$ $0.005$ $0.009$ $0.009$ $0.001$ $0.015$
$0.015$ $\beta_{2}$ $-0.003$ $0.009$ $0.009$ $-0.007$ $0.023$ $0.023$ $-0.011$
$0.045$ $0.045$ $-0.002$ $0.009$ $0.009$ $-0.004$ $0.024$ $0.024$ $-0.006$
$0.047$ $0.047$ $\beta_{3}$ $0.005$ $0.016$ $0.017$ $0.012$ $0.039$ $0.039$
$0.015$ $0.072$ $0.072$ $0.004$ $0.017$ $0.017$ $0.006$ $0.040$ $0.040$
$0.006$ $0.074$ $0.074$

Simulations show that the 2-step approach improves considerably upon smcure
for estimation of $\gamma$ when $n=200$ and the censoring rate among the
uncured observations is higher. In almost all scenarios the 2-step approach
has a smaller variance, which is expected due to presmoothing, but it also
often exhibits a lower bias. As the sample size increases or the censoring
rate decreases, we see less difference between the two methods but still, the
2-step approach is usually better. In terms of $\beta$ estimators, both
approaches give very similar results. In addition, the boxplots in Figure 1,
indicate that, even when smcure does not converge, the 2-step approach still
gives more reasonable estimates.

Since the second step of the new method does not depend on the latency model,
we expect it to be more stable than smcure with respect to misspecifications
of the latency model. We investigate this issue by considering two additional
settings: one corresponding to a non-Cox latency model (Model 5 below) for
which we still apply the two methods as if the Cox model was true and one
corresponding to a logistic-Cox model but in which we don’t use the correct
covariates. For the latter, we use Model 4, scenario 2 described above but fit
a latency model with covariates $X_{1},\dots,X_{5}$ instead of
$Z_{1},Z_{2},Z_{3}$. In particular, this means that we are including
covariates $X_{1},X_{2},X_{5}$ that actually do not have any effect and are
excluding $Z_{1}$ which affects the survival of the uncured.

Model 5. Both incidence and latency depend on three independent covariates
$X_{1}=Z_{1}\sim\text{Unif}(-1,1)$, $X_{2}=Z_{2}\sim\text{Bernoulli}(0.4)$ and
$X_{3}=Z_{3}\sim\text{Bernoulli}(0.6)$. The cure status and the survival times
for the uncured observations are generated as in Model 1 for
$\gamma=(1.4,2,1,-1)$, $\rho$ depending on $z$,
$\rho(z)=0.75+\exp(\beta^{T}z)$, $\mu=1.5$, $\beta=(1,0.4,-0.6)$. In
particular this means that the latency model does not satisfy the proportional
hazards assumption. For an observation with covariate $z$, the event time is
truncated at $\tau_{0}(z)$ equal to the $97\%$ quantile of the Weibull
distribution with parameters $\rho(z)$ and $\mu^{-1/\rho(z)}$. The censoring
times are generated according to a Weibull proportional hazards model

$S_{C}(t|x)=\exp\left(-\frac{1}{22}\mu
t^{\tilde{\rho}}\exp\left(\gamma^{T}x\right)\right),$

with $\mu=1.5$ and $\tilde{\rho}=2.5$, truncated at
$\tau=\max_{z}\tau_{0}(z)+2$. This scenario corresponds to a cure rate of
$30\%$ and a censoring rate of $45\%$.

Results for sample size $200$ and $400$, reported in Table 5, show that when
the true latency model is not a Cox proportional hazards model, even $\gamma$
estimates are biased. However, the 2-step approach has lower bias and MSE,
hence suffers less from the misspecification of the latency. On the other
hand, misspecification of the latency covariates when the model is still Cox,
seems to be less critical. It leads to a slight increase in bias and variance
compared to the results in Table 3-4 but again the 2-step approach performs
better.

Table 5: Bias, variance and MSE of $\hat{\gamma}$ and $\hat{\beta}$ for smcure
(second rows) and the 2-step approach (first rows) in Model 5 and Model 4 with
misspecification.

$n=200$ $n=400$ Mod. Par. Bias Var. MSE Bias Var. MSE 5 $\gamma_{1}$ $-0.154$
$0.187$ $0.211$ $-0.087$ $0.100$ $0.107$ $0.514$ $0.393$ $0.657$ $0.437$
$0.152$ $0.343$ $\gamma_{2}$ $-0.304$ $0.277$ $0.370$ $0.213$ $0.144$ $0.189$
$0.623$ $0.464$ $0.852$ $0.507$ $0.184$ $0.441$ $\gamma_{3}$ $-0.232$ $0.203$
$0.257$ $-0.166$ $0.108$ $0.136$ $0.297$ $0.655$ $0.744$ $0.235$ $0.176$
$0.231$ $\gamma_{4}$ $0.207$ $0.188$ $0.231$ $0.132$ $0.112$ $0.130$ $-0.326$
$0.377$ $0.483$ $-0.281$ $0.166$ $0.245$ $\beta_{1}$ $-0.281$ $0.051$ $0.130$
$-0.296$ $0.025$ $0.112$ $-0.354$ $0.056$ $0.181$ $-0.353$ $0.027$ $0.151$
$\beta_{2}$ $-0.137$ $0.043$ $0.062$ $-0.144$ $0.022$ $0.043$ $-0.174$ $0.046$
$0.076$ $-0.169$ $0.023$ $0.052$ $\beta_{3}$ $0.221$ $0.051$ $0.099$ $0.210$
$0.023$ $0.067$ $0.257$ $0.054$ $0.120$ $0.235$ $0.024$ $0.079$ $4$
$\gamma_{1}$ $0.036$ $0.250$ $0.251$ $0.016$ $0.116$ $0.116$ $0.097$ $0.360$
$0.370$ $0.041$ $0.133$ $0.135$ $\gamma_{2}$ $0.059$ $0.066$ $0.070$ $0.050$
$0.036$ $0.039$ $-0.089$ $0.092$ $0.100$ $-0.038$ $0.038$ $0.040$ $\gamma_{3}$
$-0.029$ $0.150$ $0.151$ $-0.023$ $0.075$ $0.075$ $0.025$ $0.222$ $0.223$
$0.008$ $0.091$ $0.091$ $\gamma_{4}$ $0.005$ $0.109$ $0.109$ $0.018$ $0.055$
$0.055$ $-0.059$ $0.149$ $0.152$ $-0.018$ $0.061$ $0.062$ $\gamma_{5}$
$-0.057$ $0.204$ $0.207$ $-0.022$ $0.101$ $0.102$ $0.022$ $0.307$ $0.308$
$0.032$ $0.120$ $0.121$ $\gamma_{6}$ $-0.029$ $0.204$ $0.205$ $-0.047$ $0.103$
$0.105$ $0.075$ $0.287$ $0.293$ $0.014$ $0.122$ $0.122$ $\beta_{1}$ $-0.025$
$0.028$ $0.029$ $-0.017$ $0.011$ $0.012$ $-0.004$ $0.030$ $0.030$ $-0.002$
$0.012$ $0.012$ $\beta_{2}$ $0.017$ $0.068$ $0.069$ $0.007$ $0.032$ $0.032$
$0.010$ $0.074$ $0.074$ $0.002$ $0.033$ $0.033$ $\beta_{3}$ $0.004$ $0.054$
$0.054$ $-0.001$ $0.024$ $0.024$ $0.011$ $0.058$ $0.059$ $0.002$ $0.025$
$0.025$ $\beta_{4}$ $0.013$ $0.094$ $0.095$ $0.012$ $0.043$ $0.043$ $0.004$
$0.103$ $0.103$ $0.005$ $0.044$ $0.044$ $\beta_{5}$ $-0.001$ $0.112$ $0.112$
$0.011$ $0.047$ $0.047$ $-0.019$ $0.119$ $0.119$ $0.000$ $0.049$ $0.049$

##  6 Application

In this section we illustrate the practical use of the method through two
medical datasets for melanoma cancer patients and compare the results with
those provided by the smcure package. Melanoma is a common skin cancer type
for which nowadays it is expected that a considerable fraction of the patients
get cured as a consequence of medical advances in diagnostics and treatment.
Therefore, it is important to account for the presence of cured patients in
the statistical analysis of melanoma survival data and to evaluate new
treatments focusing on cure and not only survival prolongation.

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

The ECOG phase III clinical trial e1684 aimed at evaluating the effect of
treatment (high dose interferon alpha-2b regimen) as the postoperative
adjuvant therapy for melanoma patients. The corresponding dataset, consisting
of $284$ observations (after deleting missing data), is available in the
smcure package [5]. 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). Around $30\%$ of the observations are
censored. The Kaplan-Meier curve is shown in Figure 2.

We fit a logistic-Cox mixture cure model by using the maximum likelihood
principle (smcure package) and the proposed 2-step approach. For our method we
use the smcure estimator as a preliminary estimator. In both cases, standard
errors are computed through $500$ naive bootstrap samples. The resulting
parameter estimates, standard errors and corresponding p-values for the Wald
test are reported in Table 6.

Figure 2: 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 6: Results for the incidence (logistic
component) and the latency (Cox PH component) from the ECOG data.

smcure package 2-step approach Covariates Estimates SE p-value Estimates SE
p-value incidence Intercept $1.3649$ $0.3457$ $8\cdot 10^{-5}$ $1.8351$
$0.4924$ $0.0002$ Age $0.0203$ $0.0159$ $0.2029$ $0.0388$ $0.0191$ $0.0418$
Gender $-0.0869$ $0.3347$ $0.7949$ $-0.0864$ $0.3447$ $0.8126$ Treatment
$-0.5884$ $0.3706$ $0.1123$ $-1.1096$ $0.5283$ $0.0357$ latency Age $-0.0077$
$0.0069$ $0.2663$ $-0.0103$ $0.0068$ $0.1319$ Gender $0.0994$ $0.1932$
$0.6067$ $0.0629$ $0.1831$ $0.7313$ Treatment $-0.1535$ $0.1715$ $0.3707$
$-0.0526$ $0.1904$ $0.7825$

We observe that, despite exhibiting the same effect directions for all
covariates, the two approaches give quite different results in terms of
treatment effect. Age and treatment are both found to have a significant
effect on the cure fraction when using the 2-step method, while smcure does
not detect any significant effect. We also compare the two methods in terms of
prediction accuracy for the incidence in the following way. As in [3], we
split the data into a training and a test set (at a 2:1 ratio), fit the model
in the training set and then compute the prediction error for the test set
according to the formula

$PE=-\sum_{j\in\text{test
set}}\log\left[\phi\left(\hat{\gamma}_{n}^{T}X_{j}\right)^{\hat{w}_{j}}\\{1-\phi\left(\hat{\gamma}_{n}^{T}X_{j}\right)\\}^{1-\hat{w}_{j}}\right],$

where $\hat{\gamma}_{n}$ are the parameter estimates from the training set and
$\hat{w}_{j}$ are the predicted uncure probabilities given the observations,
i.e.

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

We repeat this procedure 1000 times, for random selection of the train and
test set. The boxplot of the difference between the PE of the new method and
the PE of smcure, over these 1000 iteration, is given in Figure 3. We observe
that the 2-step approach leads to lower PE (negative difference) in more that
$50\%$ of the cases and the improvement in PE for the new method is usually
larger compared to the cases in which smcure does better.

In addition, we expect the new approach to be more stable with respect to the
latency model since that does not influence the second step of the estimation.
To illustrate this point, we also fit a cure model with only gender as
covariate for the survival of uncured patients (see Table 7) and see that in
that case, smcure also detects the effect of the treatment to be significant.

Figure 3: Boxplot of the difference between the PE of the new method and the
PE of smcure, over these 1000 iterations for the ECOG data. Table 7: Results
for the incidence (logistic component) and the latency (Cox PH component) from
the ECOG data.

smcure package 2-step approach Covariates Estimates SE p-value Estimates SE
p-value incidence Intercept $1.4000$ $0.2791$ $5\cdot 10^{-7}$ $1.9073$
$0.5225$ $0.0002$ Age $0.0165$ $0.0121$ $0.1709$ $0.0357$ $0.0174$ $0.0399$
Gender $-0.0538$ $0.3101$ $0.8623$ $-0.0429$ $0.3979$ $0.9141$ Treatment
$-0.6765$ $0.3118$ $0.0300$ $-1.2590$ $0.5809$ $0.0302$ latency Gender
$0.0637$ $0.1935$ $0.7421$ $0.0235$ $0.1924$ $0.9028$

###  6.2 Surveillance, Epidemiology and End Results database

Here we consider melanoma data extracted from the SEER database to illustrate
the performance of the method for more than one continuous covariate. The SEER
database collects cancer incidence data from population-based cancer
registries in US. We select the database ‘Incidence - SEER 18 Regs Research
Data’ and, in order to have a reasonable sample size, we extract the melanoma
cancer data for the county of San Francisco in California during the period
$2005-2010$. We consider only patients with known follow-up time and tumor
size (in the range 1-90 mm) and restrict the study to white people because of
the very small number of cases from other races. The event of interest is
death because of melanoma. This cohort consists of $384$ observations out of
which $228$ are male. The age ranges from $23$ to $101$ years old, the follow-
up from $1$ to $143$ months with no events observed after $108$ months.
Because of the high expected cure rate, $89\%$ of the observations are
censored. We consider as covariates in the model: gender (0=male, 1=female),
age and tumor size (continuous). The use of cure models is justified from the
presence of a long plateau containing around $25\%$ of the observations (see
the Kaplan-Meier curve in Figure 4).

Figure 4:  Kaplan-Meier survival curve for the SEER data.

As in the previous section, we compute parameter estimates, standard errors
and corresponding p-values for both methods (see Table 8). We observe that
both methods agree on the directions of the effects and give similar parameter
estimates. Note that the sample size in this case is larger than in the
previous data example. However, smcure only finds age to be significant while
the 2-step approach also detects the tumor size. For the latency, none of the
covariates is found significant with both methods. Also in this case, most of
the time the new method leads to an improvement in terms of prediction errors,
computed according to the procedure described in the previous subsection (see
Figure 5). In addition, we also observe that if we remove the covariate tumor
size from the latency model, the 2-step approach gives similar results while
this time smcure also detects tumor size as significant for the incidence
component (see Table 9). Once more, this behavior reflects the strong
dependence of the incidence estimates on the latency model for smcure.

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

smcure package 2-step approach Covariates Estimates SE p-value Estimates SE
p-value incidence Intercept $-4.4952$ $0.9670$ $3\cdot 10^{-6}$ $-4.9634$
$1.1269$ $10^{-5}$ Age $0.0411$ $0.0142$ $0.0037$ $0.0455$ $0.0167$ $0.0063$
Gender $-0.6805$ $0.4551$ $0.1349$ $-0.7902$ $0.5162$ $0.1259$ Tumor size
$0.0214$ $0.0136$ $0.1147$ $0.0295$ $0.0131$ $0.0239$ latency Age $-0.0161$
$0.0152$ $0.2874$ $-0.0171$ $0.0151$ $0.2572$ Gender $-0.3599$ $0.4217$
$0.3934$ $-0.3125$ $0.4184$ $0.4551$ Tumor size $0.0308$ $0.0224$ $0.1695$
$0.0282$ $0.0217$ $0.1927$

Figure 5: Boxplot of the difference between the PE of the new method and the
PE of smcure, over these 1000 iterations for the SEER data. Table 9:  Results
for the incidence (logistic component) and the latency (Cox PH component) from
the SEER data.

smcure package 2-step approach Covariates Estimates SE p-value Estimates SE
p-value incidence Intercept $-4.7297$ $0.9213$ $2\cdot 10^{-7}$ $-4.5477$
$1.0829$ $2\cdot 10^{-5}$ Age $0.0423$ $0.0137$ $0.0020$ $0.0392$ $0.0159$
$0.0134$ Gender $-0.6393$ $0.4579$ $0.1627$ $-0.7021$ $0.4430$ $0.1130$ Tumor
size $0.0288$ $0.0129$ $0.0251$ $0.0250$ $0.0129$ $0.0529$ latency Age
$-0.0165$ $0.0147$ $0.2626$ $-0.0145$ $0.0146$ $0.3201$ Gender $-0.3208$
$0.4282$ $0.4538$ $-0.2559$ $0.4092$ $0.5316$

## Appendix

##  Appendix A Proofs

###### Proof of Theorem 1.

We start by showing that, when $\tilde{\gamma}_{n}-\gamma_{0}=o_{P}(1)$, we
have

$\sup_{x\in\tilde{\mathcal{X}}}\left|\hat{\pi}(x)-\pi_{0}(x)\right|=o_{P}(1).$
(14)

Note that, independently of the choice of the trimming function,

$\mathbb{E}\left[\left\\{\left(1-\pi_{0}(X)\right)\log{\phi(\gamma^{T}X)}+\pi_{0}(X)\log\left(1-\phi(\gamma^{T}X)\right)\right\\}\tau(X)\right],$

is maximized at $\gamma=\gamma_{0}$ because of condition (I1) and the fact
that, for any $x\in\mathcal{X}$, the function

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

is strictly negative for $z\neq\phi(\gamma_{0}^{T}x)$ and
$g_{x}(\phi(\gamma_{0}^{T}x))=0$.

Next, we show uniform consistency of $\hat{\pi}_{n}$ using the decomposition
(11) and restricting to $\tilde{\mathcal{X}}$ as in (12). We have

$\displaystyle\|\hat{\pi}_{n}-\pi_{0}\|_{\infty}$
$\displaystyle=\sup_{x\in\tilde{\mathcal{X}}}\left|\hat{g}_{\tilde{\gamma}_{n}}\left(\tilde{\gamma}^{T}_{n}x\right)-g_{\gamma_{0}}\left(\gamma^{T}_{0}x\right)\right|$
(15)
$\displaystyle=\sup_{x\in\tilde{\mathcal{X}}}\left|\hat{g}_{{\gamma}_{0}}\left({\gamma}^{T}_{0}x\right)-g_{\gamma_{0}}\left(\gamma^{T}_{0}x\right)\right|+\sup_{x\in\tilde{\mathcal{X}}}\left|\hat{g}_{\tilde{\gamma}_{n}}\left(\tilde{\gamma}^{T}_{n}x\right)-\hat{g}_{\gamma_{0}}\left(\gamma^{T}_{0}x\right)\right|$
$\displaystyle\leq\sup_{x\in\tilde{\mathcal{X}}}\sup_{t\leq\tau_{0}}\left|\hat{F}_{n}\left(t\mid\gamma^{T}_{0}X=\gamma^{T}_{0}x\right)-{F}_{T}\left(t\mid\gamma^{T}_{0}X=\gamma^{T}_{0}x\right)\right|$
$\displaystyle\qquad+\|\tilde{\gamma}^{T}_{n}-\gamma_{0}\|\sup_{x\in\tilde{\mathcal{X}}}\sup_{\gamma\in
G}\left|\nabla_{\gamma}\hat{g}_{{\gamma}}\left({\gamma}^{T}x\right)\right|.$

The first term on the right-hand side of the equation converges to zero by
Theorem 4.1 in [38]. The second term converges to zero because of assumption
(C1) and the fact that

$\sup_{x\in\tilde{\mathcal{X}}}\sup_{\gamma\in
G}\left|\nabla_{\gamma}\hat{g}_{{\gamma}}\left({\gamma}^{T}x\right)\right|=O_{P}(1).$
(16)

Indeed, (16) follows from assumption (C9) and

$\sup_{x\in\tilde{\mathcal{X}}}\sup_{\gamma\in
G}\left|\nabla_{\gamma}\hat{g}_{{\gamma}}\left({\gamma}^{T}x\right)-\nabla_{\gamma}g_{\gamma}\left({\gamma}^{T}x\right)\right|=o_{P}(1),$
(17)

which can be proved as in Lemma A.2 in [17]. Their estimator
$\hat{G}_{\theta}(t|\lambda(\theta,x))$ is the same as our
$\hat{F}_{n}(t|\gamma^{T}X=\gamma^{T}x)$ if we replace $\theta$ by $\gamma$,
consider $\lambda(\theta,x)=\theta^{T}x$ and exchange $T$ with $C$ (they are
interested in the conditional distribution of $C$ given $\lambda(\theta,X)$).
This concludes the proof of (14).

The consistency of $\hat{\gamma}_{n}$ then follows similarly to Theorem 1 in
[24]. There it is required that
$\sup_{x\in\tilde{\mathcal{X}}}\left|\hat{\pi}(x)-\pi_{0}(x)\right|\to 0$
almost surely in order to obtain strong consistency of $\hat{\gamma}_{n}$.
Here we restrict to convergence in probability and conclude via Theorem 1 in
[6] formulated for M-estimators. ∎

###### Proof of Theorem 2.

The result follows from Theorem 2 in [24] once we show that the assumptions
(AN1)-(AN4) of that paper hold. Note that the introduction of the trimming
function $\tau(\cdot)$ would not change anything in the proof. It just allows
us to restrict ourselves to the set $\tilde{\mathcal{X}}$ where the density of
the index is bounded from below, in order to apply standard results from the
literature. The assumptions (AN1) and (AN3) of [24] are the same as
assumptions (C2), (N1) and (N2) here. It remains to verify assumptions (AN2)
and (AN4) which for completeness we state below:

  * (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\tilde{\mathcal{X}}}|\pi(x)|$.

  * (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}^{T}X)}+\frac{1}{1-\phi(\gamma_{0}^{T}X)}\right)\phi^{\prime}(\gamma_{0}^{T}X)X{\tau(X)}\right]\qquad\qquad\\\
&\qquad=\frac{1}{n}\sum_{i=1}^{n}\Psi(Y_{i},\Delta_{i},X_{i},Z_{i})+R_{n},\end{split}$

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

Step 1. For (AN2) we use the class of functions

$\Pi=\left\\{f:\mathcal{X}\to[0,1],\,f(x)=g(\gamma^{T}x)\text{ for some
}\gamma\in G\text{ and }g\in\tilde{\Pi}_{\gamma}\right\\},$

where $\tilde{\Pi}_{\gamma}$ is the space of continuously differentiable
functions $g$ from $\mathcal{X}_{\gamma}$ to $[0,1]$ such that
$\sup_{u\in\mathcal{X}_{\gamma}}|g^{\prime}(u)|\leq M$ and

$\sup_{\begin{subarray}{c}u_{1}\neq u_{2}\\\
u_{1},u_{2}\in\mathcal{X}_{\gamma}\end{subarray}}\frac{|g^{\prime}(u_{1})-g^{\prime}(u_{2})|}{|u_{1}-u_{2}|^{\xi}}\leq
M,$

for some $M>0$ and $\xi\in(0,1]$ independent of $\gamma$. The norm that we
consider on $\Pi$ is the sup norm. By assumption (A2) we have $\pi_{0}\in\Pi$.
Next we compute the $\epsilon$-covering number of $\Pi$. Let
$\gamma_{1},\dots,\gamma_{N_{1}}$ be a $\delta$-covering of the compact
$G\in\mathbb{R}^{p}$ with respect to the $l_{1}$ norm. We have
$N_{1}\leq\frac{K_{1}}{\delta^{p}}$ for some constant $K_{1}>0$. Consider the
class $\tilde{\Pi}$ of continuously differentiable functions $g$ from
$\mathcal{U}=\\{\gamma^{T}x\mid\gamma\in{G},x\in\mathcal{X}\\}$ to $[0,1]$
such that $\sup_{u\in\mathcal{U}}|g^{\prime}(u)|\leq M$ and

$\sup_{\begin{subarray}{c}u_{1}\neq u_{2}\\\
u_{1},u_{2}\in\mathcal{U}\end{subarray}}\frac{|g^{\prime}(u_{1})-g^{\prime}(u_{2})|}{|u_{1}-u_{2}|^{\xi}}\leq
M.$

We have in particular that if $g\in\tilde{\Pi}$, then the restriction of $g$
to $\mathcal{X}_{\gamma}$ belongs to the class $\tilde{\Pi}_{\gamma}$ for any
$\gamma\in\mathcal{G}$. Let $g_{1},\dots,g_{N_{2}}$ be a $\rho$-covering of
$\tilde{\Pi}$ with respect to the sup norm. Since $\mathcal{U}$ is bounded and
convex ($G$ can be chosen to be a small neighborhood of $\gamma_{0}$), from
Theorem 2.7.1 in [36], we have

$\log N_{2}\leq K_{2}\left(\frac{1}{\rho}\right)^{1/(1+\xi)},$

for some constant $K_{2}>0$. For $1\leq i\leq N_{1}$ and $1\leq j\leq N_{2}$,
define $f_{i,j}(x)=g_{j}(\gamma^{T}_{i}x)$. We show that $\\{f_{i,j}\\}$,
$1\leq i\leq N_{1}$ and $1\leq j\leq N_{2}$ is an $\epsilon$-covering of
$\Pi$. Let $f\in\Pi$, $f(x)=g(\gamma^{T}x)$. From Whitney’s theorem it follows
that $g\in\tilde{\Pi}_{\gamma}$ can be extended to a function
$\bar{g}\in\tilde{\Pi}$, i.e. $\bar{g}(u)=g(u)$ for
$u\in\mathcal{X}_{\gamma}$. Let $i$ and $j$ be such that
$\sup_{u\in\mathcal{U}}|\bar{g}(u)-g_{j}(u)|\leq\rho$ and
$\|\gamma-\gamma_{i}\|\leq\delta.$ We then have

$\displaystyle\|f-f_{i,j}\|_{\infty}$
$\displaystyle\leq\sup_{x\in\mathcal{X}}\left|g\left(\gamma^{T}x\right)-g_{j}\left(\gamma^{T}x\right)\right|+\sup_{x\in\mathcal{X}}\left|g_{j}\left(\gamma^{T}x\right)-g_{j}\left(\gamma^{T}_{i}x\right)\right|$
$\displaystyle\leq\sup_{u\in\mathcal{U}}|\bar{g}(u)-g_{j}(u)|+K_{3}\|\gamma-\gamma_{i}\|_{1}$
$\displaystyle\leq\rho+MK_{3}\delta.$

Hence, if we take $\rho=\epsilon/2$ and $\delta=\epsilon/MK_{3}$ then we get
an $\epsilon$-covering of $\Pi$. It follows that

$\log N(\epsilon,\Pi,\|\cdot\|_{\infty})\leq
c_{1}-p\log\epsilon+c_{2}\left(\frac{1}{\epsilon}\right)^{1/(1+\xi)}\qquad
c_{1},c_{2}>0,$

and as a result (AN2) is satisfied.

Step 2. For assumption (AN4))(i), note that
$\hat{\pi}_{n}(x)=\hat{g}_{\tilde{\gamma}_{n}}(\tilde{\gamma}^{T}_{n}x)$ with
$\hat{g}_{\gamma}$ as in (9). From consistency of $\tilde{\gamma}_{n}$ it
follows that $\tilde{\gamma}_{n}\in G$. Let
${\hat{g}}^{\prime}_{\tilde{\gamma}_{n}}$ denote the derivative of the
function $u\mapsto\hat{g}_{\tilde{\gamma}_{n}}(u)$. Let $g_{\gamma}$ be as in
(10). Next we show that

$\sup_{x\in{\tilde{\mathcal{X}}}}|\hat{g}^{\prime}_{\tilde{\gamma}_{n}}(\tilde{\gamma}^{T}_{n}x)-g^{\prime}_{\gamma_{0}}(\gamma^{T}_{0}x)|=o_{P}(1),$
(18)

and

$\sup_{\begin{subarray}{c}x_{1},x_{2}\in{\tilde{\mathcal{X}}}\\\
\tilde{\gamma}^{T}_{n}x_{1}\neq\tilde{\gamma}^{T}_{n}x_{2}\end{subarray}}\frac{\left|\hat{g}^{\prime}_{\tilde{\gamma}_{n}}\left(\tilde{\gamma}^{T}_{n}x_{1}\right)-g^{\prime}_{\gamma_{0}}\left(\gamma^{T}_{0}x_{1}\right)-\hat{g}^{\prime}_{\tilde{\gamma}_{n}}\left(\tilde{\gamma}^{T}_{n}x_{2}\right)+g^{\prime}_{\gamma_{0}}\left(\gamma^{T}_{0}x_{2}\right)\right|}{|\tilde{\gamma}^{T}_{n}x_{1}-\tilde{\gamma}^{T}_{n}x_{2}|^{\xi}}=o_{P}(1),$
(19)

from which we can derive that (AN4)(i) is satisfied since
$\sup_{x\in{\tilde{\mathcal{X}}}}|g^{\prime}_{\gamma_{0}}(\gamma^{T}_{0}x)|<\infty$
and

$\sup_{\begin{subarray}{c}x_{1},x_{2}\in{\tilde{\mathcal{X}}}\\\
{\gamma}^{T}_{0}x_{1}\neq{\gamma}^{T}_{0}x_{2}\end{subarray}}\frac{\left|g^{\prime}_{\gamma_{0}}\left(\gamma^{T}_{0}x_{2}\right)-g^{\prime}_{\gamma_{0}}\left(\gamma^{T}_{0}x_{1}\right)\right|}{|{\gamma}^{T}_{0}x_{1}-{\gamma}^{T}_{0}x_{2}|^{\xi}}<\infty.$

From Theorems 4.1 and 4.2 in [38] we have

$\sup_{x\in{\tilde{\mathcal{X}}}}|\hat{g}^{\prime}_{{\gamma}_{0}}({\gamma}^{T}_{0}x)-g^{\prime}_{\gamma_{0}}(\gamma^{T}_{0}x)|=o_{P}(1),$
(20)

and

$\sup_{\begin{subarray}{c}x_{1},x_{2}\in{\tilde{\mathcal{X}}}\\\
{\gamma}^{T}_{0}x_{1}\neq{\gamma}^{T}_{0}x_{2}\end{subarray}}\frac{\left|\hat{g}^{\prime}_{{\gamma}_{0}}\left({\gamma}^{T}_{0}x_{1}\right)-g^{\prime}_{\gamma_{0}}\left(\gamma^{T}_{0}x_{1}\right)-\hat{g}^{\prime}_{{\gamma}_{0}}\left({\gamma}^{T}_{0}x_{2}\right)+g^{\prime}_{\gamma_{0}}\left(\gamma^{T}_{0}x_{2}\right)\right|}{|{\gamma}^{T}_{0}x_{1}-{\gamma}^{T}_{0}x_{2}|^{\xi}}=o_{P}(1),$
(21)

by conditioning on the variable $\gamma_{0}^{T}X$. Next, we need to deal with
the fact that we are using $\tilde{\gamma}_{n}^{T}X$ instead of
$\gamma_{0}^{T}X$.

We can write

$\log\hat{g}_{\tilde{\gamma}_{n}}(u)=\sum_{i=1}^{n}\Delta_{i}\log\left(1-W_{i}(u;\tilde{\gamma}_{n})\right),$

where

$W_{i}(u;\tilde{\gamma}_{n})=k\left(\frac{\tilde{\gamma}^{T}_{n}X_{i}-u}{b}\right)\bigg{/}\left(\sum_{j=1}^{n}\mathds{1}_{\\{T_{j}\geq
T_{i}\\}}k\left(\frac{\tilde{\gamma}^{T}_{n}X_{j}-u}{b}\right)\right).$

Hence

$\hat{g}^{\prime}_{\tilde{\gamma}_{n}}(u)=\hat{g}_{\tilde{\gamma}_{n}}(u)\sum_{i=1}^{n}\Delta_{i}\frac{W^{\prime}_{i}(u;\tilde{\gamma}_{n})}{1-W_{i}(u;\tilde{\gamma}_{n})}.$
(22)

Using the triangular inequality we obtain

$\displaystyle\sup_{x\in{\tilde{\mathcal{X}}}}\left|\hat{g}^{\prime}_{\tilde{\gamma}_{n}}(\tilde{\gamma}_{n}^{T}x)-\hat{g}^{\prime}_{{\gamma}_{0}}({\gamma}_{0}^{T}x)\right|$
(23)
$\displaystyle\leq\sup_{x\in{\tilde{\mathcal{X}}}}\left|\hat{g}_{\tilde{\gamma}_{n}}(\tilde{\gamma}_{n}^{T}x)-\hat{g}_{{\gamma}_{0}}({\gamma}_{0}^{T}x)\right|\sup_{x\in{\tilde{\mathcal{X}}}}\left|\sum_{i=1}^{n}\Delta_{i}\frac{W^{\prime}_{i}(\tilde{\gamma}_{n}^{T}x;\tilde{\gamma}_{n})}{1-W_{i}(\tilde{\gamma}_{n}^{T}x;\tilde{\gamma}_{n})}\right|$
$\displaystyle+\sup_{x\in{\tilde{\mathcal{X}}}}|\hat{g}^{\prime}_{{\gamma}_{0}}({\gamma}_{0}^{T}x)|\sup_{x\in\mathcal{X}}\left|\sum_{i=1}^{n}\Delta_{i}\left\\{\frac{W^{\prime}_{i}(\tilde{\gamma}_{n}^{T}x;\tilde{\gamma}_{n})}{1-W_{i}(\tilde{\gamma}_{n}^{T}x;\tilde{\gamma}_{n})}-\frac{W^{\prime}_{i}({\gamma}_{0}^{T}x;{\gamma}_{0})}{1-W_{i}({\gamma}_{0}^{T}x;{\gamma}_{0})}\right\\}\right|.$

As is (15), we have

$\sup_{x\in{\tilde{\mathcal{X}}}}\left|\hat{g}_{\tilde{\gamma}_{n}}(\tilde{\gamma}_{n}^{T}x)-\hat{g}_{{\gamma}_{0}}({\gamma}_{0}^{T}x)\right|=O_{P}(n^{-1/2}).$

Moreover,
$\sup_{x\in{\tilde{\mathcal{X}}}}|\hat{g}^{\prime}_{{\gamma}_{0}}({\gamma}_{0}^{T}x)|=O_{P}(1)$
and, as in the proof of Proposition 4.1 in [38] it can be seen that

$\sup_{x\in{\tilde{\mathcal{X}}}}\left|\sum_{i=1}^{n}\Delta_{i}\frac{W^{\prime}_{i}({\gamma}_{0}^{T}x;{\gamma}_{0})}{1-W_{i}({\gamma}_{0}^{T}x;{\gamma}_{0})}\right|=O_{P}(1).$

and

$\sup_{x\in{\tilde{\mathcal{X}}}}\left|\sum_{i=1}^{n}\Delta_{i}\nabla_{\gamma}\left(\frac{W^{\prime}_{i}({\gamma}^{T}x;{\gamma})}{1-W_{i}({\gamma}^{T}x;{\gamma})}\right)\bigg{|}_{\gamma_{0}}\right|=O_{P}(1/b).$

From a Taylor expansion and $\tilde{\gamma}_{n}-\gamma_{0}=O_{P}(n^{-1/2})$,
we conclude that the left hand side of (23) converges to zero. Together with
(20) we obtain (18).

Next we show that

$\sup_{\begin{subarray}{c}x_{1},x_{2}\in{\tilde{\mathcal{X}}}\\\
\tilde{\gamma}^{T}_{n}x_{1}\neq\tilde{\gamma}^{T}_{n}x_{2}\end{subarray}}\frac{\left|\hat{g}^{\prime}_{\tilde{\gamma}_{n}}\left(\tilde{\gamma}^{T}_{n}x_{1}\right)-\hat{g}^{\prime}_{\gamma_{0}}\left(\gamma^{T}_{0}x_{1}\right)-\hat{g}^{\prime}_{\tilde{\gamma}_{n}}\left(\tilde{\gamma}^{T}_{n}x_{2}\right)+\hat{g}^{\prime}_{\gamma_{0}}\left(\gamma^{T}_{0}x_{2}\right)\right|}{|\tilde{\gamma}^{T}_{n}x_{1}-\tilde{\gamma}^{T}_{n}x_{2}|^{\xi}}=o_{P}(1).$

Note that by (22) and the fact that
$\sup_{x\in\mathcal{X}}|\hat{g}^{\prime}_{{\gamma}}({\gamma}^{T}x)|=O_{P}(1)$
for both $\gamma=\gamma_{0}$ and $\gamma=\tilde{\gamma}_{n}$, it is sufficient
to consider

$\sup_{\begin{subarray}{c}x_{1},x_{2}\in{\tilde{\mathcal{X}}}\\\
\tilde{\gamma}^{T}_{n}x_{1}\neq\tilde{\gamma}^{T}_{n}x_{2}\end{subarray}}\frac{\left|V_{\tilde{\gamma}_{n}}\left(\tilde{\gamma}^{T}_{n}x_{1}\right)-V_{\gamma_{0}}\left(\gamma^{T}_{0}x_{1}\right)-V_{\tilde{\gamma}_{n}}\left(\tilde{\gamma}^{T}_{n}x_{2}\right)+V_{\gamma_{0}}\left(\gamma^{T}_{0}x_{2}\right)\right|}{|\tilde{\gamma}^{T}_{n}x_{1}-\tilde{\gamma}^{T}_{n}x_{2}|^{\xi}},$
(24)

where

$V_{\gamma}(u)=\sum_{i=1}^{n}\Delta_{i}\frac{W^{\prime}_{i}(u;{\gamma})}{1-W_{i}(u;{\gamma})}.$

In addition we can also restrict the supremum over $x_{1},x_{2}$ such that
$\tilde{\gamma}^{T}_{n}x_{1}\neq\tilde{\gamma}^{T}_{n}x_{2}$ and
${\gamma}^{T}_{0}x_{1}\neq{\gamma}^{T}_{0}x_{2}$ because otherwise we would
have $|\tilde{\gamma}^{T}_{n}x_{1}-\tilde{\gamma}^{T}_{n}x_{2}|\leq cn^{-1/2}$
and (24) would obviously hold for $\xi$ sufficiently small by using a Taylor
expansion and $\sup_{u}V^{\prime}_{\gamma}(u)=O_{P}(1/b)$. Note also that

$\sup_{\begin{subarray}{c}x_{1},x_{2}\in{\tilde{\mathcal{X}}}\\\
\tilde{\gamma}^{T}_{n}x_{1}\neq\tilde{\gamma}^{T}_{n}x_{2},\,{\gamma}^{T}_{0}x_{1}\neq{\gamma}^{T}_{0}x_{2}\end{subarray}}\frac{|{\gamma}^{T}_{0}x_{1}-{\gamma}^{T}_{0}x_{2}|}{|\tilde{\gamma}^{T}_{n}x_{1}-\tilde{\gamma}^{T}_{n}x_{2}|}=O_{P}(1),$

hence the denominator can be replaced by
$|{\gamma}^{T}_{0}x_{1}-{\gamma}^{T}_{0}x_{2}|^{\xi}$. By Taylor expansion we
obtain that the largest order term in (24) is

$(\tilde{\gamma}_{n}-\gamma_{0})\sup_{\begin{subarray}{c}x_{1},x_{2}\in{\tilde{\mathcal{X}}}\\\
{\gamma}^{T}_{0}x_{1}\neq{\gamma}^{T}_{0}x_{2}\end{subarray}}\sum_{i=1}^{n}\Delta_{i}\frac{\nabla_{\gamma}\left(\frac{W^{\prime}_{i}({\gamma}^{T}x_{1};{\gamma})}{1-W_{i}({\gamma}^{T}x_{1};{\gamma})}\right)\bigg{|}_{\gamma_{0}}-\nabla_{\gamma}\left(\frac{W^{\prime}_{i}({\gamma}^{T}x_{2};{\gamma})}{1-W_{i}({\gamma}^{T}x_{2};{\gamma})}\right)\bigg{|}_{\gamma_{0}}}{|{\gamma}^{T}_{0}x_{1}-{\gamma}^{T}_{0}x_{2}|^{\xi}}.$

Using Lemma A.1. in [38] and that for a smooth function $l$ (in our case $l$
is the kernel function or its derivatives) we have

$|l(\frac{x-z}{b})-l(\frac{x-z}{b})|/|x-y|^{\xi}=O(b^{-\xi}),$

uniformly over $x,y,z$, we obtain that the expression in the previous equation
is of the order $O_{P}(n^{-1/2}b^{-1-\xi})=o_{P}(1)$. Together with (21) this
yields (19).

Step 3. Assumption (AN4)(ii) can be checked using (15). The first term on the
right-hand side of that equation is of order $o_{P}(n^{-1/4})$ because of
Theorem 4.1 in [38] and the assumption (C6). The second term on the right-hand
side of (15) is also of order $o_{P}(n^{-1/4})$ because of assumptions (N3),
(C9) and (17).

Step 4. For (AN4)(iii), using the decomposition (11) we have

$\displaystyle\mathbb{E}^{*}\left[\left(\hat{\pi}(X)-\pi_{0}(X)\right)\left(\frac{1}{\phi(\gamma_{0}^{T}X)}+\frac{1}{1-\phi(\gamma_{0}^{T}X)}\right)\phi^{\prime}(\gamma_{0}^{T}X)X{\tau(X)}\right]$
(25)
$\displaystyle=\mathbb{E}^{*}\left[\frac{\hat{F}_{n}(\tau_{0}\mid\gamma^{T}_{0}X)-{F}_{T}(\tau_{0}\mid\gamma^{T}_{0}X)}{\phi(\gamma_{0}^{T}X)\left\\{1-\phi(\gamma_{0}^{T}X)\right\\}}\phi^{\prime}(\gamma_{0}^{T}X)X{\tau(X)}\right]$
$\displaystyle\quad+\mathbb{E}^{*}\left[\frac{\hat{g}_{\tilde{\gamma}_{n}}(\tilde{\gamma}^{T}_{n}X)-\hat{g}_{\gamma_{0}}(\gamma^{T}_{0}X)}{\phi(\gamma_{0}^{T}X)\left\\{1-\phi(\gamma_{0}^{T}X)\right\\}}\phi^{\prime}(\gamma_{0}^{T}X)X{\tau(X)}\right].$

Since $\sup_{x\in\mathcal{X}}H([t,\infty)|\gamma^{T}_{0}x)<1$ by condition
(8), from Theorem 3.2 in [10], it follows as in [24] (see proof of Theorem 5)
that the first term on the right-hand side of (LABEL:eqn:expectation*) is
equal to $\frac{1}{n}\sum_{i=1}^{n}\psi(Y_{i},\Delta_{i},X_{i})+R_{n}$ with
$\|R_{n}\|=o_{P}(n^{-1/2})$ and

$\begin{split}\psi(Y,\Delta,X)=-\left\\{\frac{\Delta\mathds{1}_{\\{Y\leq\tau_{0}\\}}}{H\left([Y,\infty)|\gamma^{T}_{0}X\right)}-\int_{0}^{Y\wedge\tau_{0}}\frac{H_{1}\left(ds|\gamma^{T}_{0}X\right)}{H^{2}\left([s,\infty)|\gamma^{T}_{0}X\right)}\right\\}\frac{\phi^{\prime}\left(\gamma_{0}^{T}X\right)}{\phi\left(\gamma_{0}^{T}X\right)}X{\tau(X)}.\end{split}$
(26)

For the second term in (LABEL:eqn:expectation*), by the mean value theorem, we
write

$\displaystyle\hat{g}_{\tilde{\gamma}_{n}}\left(\tilde{\gamma}^{T}_{n}x\right)-\hat{g}_{\gamma_{0}}\left(\gamma^{T}_{0}x\right)$
$\displaystyle=(\tilde{\gamma}_{n}-\gamma_{0})^{T}\nabla_{\gamma}\hat{g}_{\gamma}\left(\gamma^{T}x\right)\big{|}_{\gamma^{*}}$
$\displaystyle=(\tilde{\gamma}_{n}-\gamma_{0})^{T}\nabla_{\gamma}\hat{g}_{\gamma}\left(\gamma^{T}x\right)\big{|}_{\gamma^{*}},$

for some $\gamma^{*}$ such that
$\|\gamma^{*}-\gamma_{0}\|\leq\|\tilde{\gamma}_{n}-\gamma_{0}\|$. Using (17)
and assumption (C9), we obtain

$\hat{g}_{\tilde{\gamma}_{n}}\left(\tilde{\gamma}^{T}_{n}x\right)-\hat{g}_{\gamma_{0}}\left(\gamma^{T}_{0}x\right)=(\tilde{\gamma}_{n}-\gamma_{0})^{T}\nabla_{\gamma}{g}_{\gamma}\left(\gamma^{T}x\right)\big{|}_{\gamma_{0}}+o_{P}(n^{-1/2}),$

where the $o_{P}(n^{-1/2})$ term is uniform with respect to
$x\in{\tilde{\mathcal{X}}}$. It follows that

$\displaystyle\mathbb{E}^{*}\left[\frac{\hat{g}_{\tilde{\gamma}_{n}}(\tilde{\gamma}^{T}_{n}X)-\hat{g}_{\gamma_{0}}(\gamma^{T}_{0}X)}{\phi(\gamma_{0}^{T}X)\left\\{1-\phi(\gamma_{0}^{T}X)\right\\}}\phi^{\prime}(\gamma_{0}^{T}X)X{\tau(X)}\right]$
$\displaystyle=\mathbb{E}\left[\frac{\nabla_{\gamma}{g}_{\gamma}(\gamma^{T}X)^{T}\big{|}_{\gamma_{0}}}{\phi(\gamma_{0}^{T}X)\left\\{1-\phi(\gamma_{0}^{T}X)\right\\}}\phi^{\prime}(\gamma_{0}^{T}X)X{\tau(X)}\right](\tilde{\gamma}_{n}-\gamma_{0})+o_{P}(n^{-1/2}),$

because $E^{*}$ is just the expectation with respect to the variable $X$. By
definition of $g_{\gamma}(\gamma^{T}x)$ and a Taylor expansion we have

$\displaystyle g_{\gamma}\left(\gamma^{T}x\right)$
$\displaystyle=\mathbb{E}\left[1-\phi\left(\gamma_{0}^{T}X\right)|\gamma^{T}X=\gamma^{T}x\right]$
$\displaystyle=1-\phi\left(\gamma^{T}x\right)+(\gamma-\gamma_{0})^{T}\mathbb{E}\left[X|\gamma^{T}X=\gamma^{T}x\right]\phi^{\prime}\left(\gamma^{T}x\right)$
$\displaystyle\quad+\frac{1}{2}(\gamma-\gamma_{0})^{T}\mathbb{E}\left[\phi^{\prime\prime}\left(\gamma^{T}_{*}x\right)XX^{T}|\gamma^{T}X=\gamma^{T}x\right](\gamma-\gamma_{0}),$

for some $\|\gamma_{*}-\gamma_{0}\|\leq\|\gamma-\gamma_{0}\|$. Hence

$\nabla_{\gamma}{g}_{\gamma}\left(\gamma^{T}x\right)\big{|}_{\gamma_{0}}=-\phi^{\prime}\left(\gamma^{T}_{0}x\right)x+\mathbb{E}\left[X|\gamma^{T}_{0}X=\gamma^{T}_{0}x\right]\phi^{\prime}\left(\gamma^{T}_{0}x\right).$

and consequently

$\displaystyle\mathbb{E}\left[\frac{\nabla_{\gamma}{g}_{\gamma}(\gamma^{T}X)^{T}\big{|}_{\gamma_{0}}}{\phi(\gamma_{0}^{T}X)\left\\{1-\phi(\gamma_{0}^{T}X)\right\\}}\phi^{\prime}(\gamma_{0}^{T}X)X{\tau(X)}\right]$
$\displaystyle=-\mathbb{E}\left[\frac{\phi^{\prime}(\gamma^{T}_{0}X)X^{T}}{\phi(\gamma_{0}^{T}X)\left\\{1-\phi(\gamma_{0}^{T}X)\right\\}}\phi^{\prime}(\gamma_{0}^{T}X)X{\tau(X)}\right]$
$\displaystyle\quad+\mathbb{E}\left[\frac{\mathbb{E}[X^{T}|\gamma^{T}_{0}X]\phi^{\prime}(\gamma^{T}_{0}X)}{\phi(\gamma_{0}^{T}X)\left\\{1-\phi(\gamma_{0}^{T}X)\right\\}}\phi^{\prime}(\gamma_{0}^{T}X)X{\tau(X)}\right]$
$\displaystyle=\mathbb{E}\left[\frac{\phi^{\prime}(\gamma^{T}_{0}X)^{2}\left\\{\mathbb{E}[X^{T}|\gamma_{0}^{T}X]\mathbb{E}[X\tau(X)|\gamma_{0}^{T}X]-\mathbb{E}[X^{T}X{\tau(X)}|\gamma_{0}^{T}X]\right\\}}{\phi(\gamma_{0}^{T}X)\left\\{1-\phi(\gamma_{0}^{T}X)\right\\}}\right],$

We denote this expression by $Q$. This yields

$\displaystyle\mathbb{E}^{*}\left[\frac{\hat{g}_{\tilde{\gamma}_{n}}(\tilde{\gamma}^{T}_{n}X)-\hat{g}_{\gamma_{0}}(\gamma^{T}_{0}X)}{\phi(\gamma_{0}^{T}X)\left\\{1-\phi(\gamma_{0}^{T}X)\right\\}}\phi^{\prime}(\gamma_{0}^{T}X)X{\tau(X)}\right]=Q(x)(\tilde{\gamma}_{n}-\gamma_{0})+o_{P}(n^{-1/2}),$

and by assumption (N3)

$\displaystyle\mathbb{E}^{*}\left[\left(\hat{\pi}(X)-\pi_{0}(X)\right)\left(\frac{1}{\phi(\gamma_{0}^{T}X)}+\frac{1}{1-\phi(\gamma_{0}^{T}X)}\right)\phi^{\prime}(\gamma_{0}^{T}X)X{\tau(X)}\right]$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\Psi(Y_{i},\Delta_{i},X_{i})+Q\sum_{i=1}^{n}\zeta(Y_{i},\Delta_{i},X_{i},Z_{i})+R_{n},$

with $\|R_{n}\|=o_{P}(n^{-1/2})$. This means that assumption (AN4)(iii) of
[24] is satisfied with
$\varPsi(Y,\Delta,X,Z)=\psi(Y,\Delta,X,Z)+\zeta(Y,\Delta,X,Z)$.

This concludes the verification of the assumptions of Theorem 3 in [24]. It
also follows that the covariance matrix is given by

$\Sigma_{\gamma}=(\Gamma^{\prime}_{1}\Gamma_{1})^{-1}\Gamma^{\prime}_{1}V\Gamma_{1}(\Gamma^{\prime}_{1}\Gamma_{1})^{-1}=\Gamma_{1}^{-1}V\Gamma_{1}^{-1},$
(27)

where $V=Var(\Psi(Y,\Delta,X,Z))$ and

$\Gamma_{1}=-\mathbb{E}\left[\left(\frac{1}{\phi(\gamma_{0}^{T}X)}+\frac{1}{1-\phi(\gamma_{0}^{T}X)}\right)\phi^{\prime}(\gamma_{0}^{T}X)^{2}XX^{T}{\tau(X)}\right].$

∎

[Acknowledgments] I. Van Keilegom acknowledges financial support from the
European Research Council (2016-2021, Horizon 2020 and grant agreement
694409). V. Patilea gratefully acknowledges support from the Joint Research
Initiative ‘Models and mathematical processing of very large data’ under the
aegis of Risk Foundation, in partnership with MEDIAMETRIE and GENES, France.
For the simulations we used the Lisa cluster of the Dutch national
Supercomputer. The data that were analysed in Section 5 have been provided by
the SEER programme (www.seer.cancer.gov) research data (1973–2014), National
Cancer Institute, Division of Cancer Control and Population Sciences,
Surveillance Research Program.

## References

  * [1] Aerts, M., Hens, N., and Simonoff, J. S. (2010).  Model selection in regression based on pre-smoothing.  _J. Appl. Stat._ 37, 9, 1455–1472.
  * [2] Amico, M. and Van Keilegom, I. (2018).  Cure models in survival analysis.  _Annu. Rev. Stat. Appl._  _5_ , 311–342.
  * [3] Amico, M., Van Keilegom, I., and Legrand, C. (2019).  The single-index/cox mixture cure model.  _Biometrics_ 75, 2, 452–462.
  * [4] Burke, K. and Patilea, V. (2020).  A likelihood-based approach for cure regression models.  _TEST_ , 1–20.
  * [5] Cai, C., Zou, Y., Peng, Y., and Zhang, J. (2012).  smcure: An r-package for estimating semiparametric mixture cure models.  _Comput. Meth. Prog. Bio._ 108, 3, 1255–1260.
  * [6] Chen, X., Linton, O., and Van Keilegom, I. (2003).  Estimation of semiparametric models when the criterion function is not smooth.  _Econometrica_ 71, 5, 1591–1608.
  * [7] Cristobal, J. C., Roca, P. F., and González-Manteiga, W. (1987).  A class of linear regression parameter estimators constructed by nonparametric estimation.  _Ann. Stat._ , 603–609.
  * [8] Dirick, L., Bellotti, T., Claeskens, G., and Baesens, B. (2019).  Macro-economic factors in credit risk calculations: including time-varying covariates in mixture cure models.  _J. Bus. Econ. Stat._ 37, 1, 40–53.
  * [9] Dirick, L., Claeskens, G., and Baesens, B. (2017).  Time to default in credit scoring using survival analysis: a benchmark study.  _J. Oper. Res. Soc._ 68, 6, 652–665.
  * [10] Du, Y. and Akritas, M. (2002).  Uniform strong representation of the conditional kaplan-meier process.  _Math. Methods Stat._ 11, 2, 152–182.
  * [11] Farewell, V. T. (1982).  The use of mixture models for the analysis of survival data with long-term survivors.  _Biometrics_ , 1041–1046.
  * [12] Ferraty, F., González-Manteiga, W., Martínez-Calvo, A., and Vieu, P. (2012).  Presmoothing in functional linear regression.  _Stat. Sin._ , 69–94.
  * [13] Han, X. (2017).  Statistical methods for analysis of genetic and survival data with latent heterogeneity.  Ph.D. thesis, New York University.
  * [14] Kuk, A. Y. and Chen, C.-H. (1992).  A mixture model combining logistic regression with proportional hazards regression.  _Biometrika_ 79, 3, 531–541.
  * [15] Legrand, C. and Bertrand, A. (2019).  Cure models in cancer clinical trials.  In _Textbook of Clinical Trials in Oncology_. Chapman and Hall/CRC, 465–492.
  * [16] Li, C.-S. and Taylor, J. M. (2002).  A semi-parametric accelerated failure time cure model.  _Stat. Med._ 21, 21, 3235–3247.
  * [17] Lopez, O., Patilea, V., and Van Keilegom, I. (2013).  Single index regression models in the presence of censoring depending on the covariates.  _Bernoulli_ 19, 3, 721–747.
  * [18] López-Cheda, A., Cao, R., Jácome, M. A., and Van Keilegom, I. (2017).  Nonparametric incidence estimation and bootstrap bandwidth selection in mixture cure models.  _Computational Statistics & Data Analysis_  _105_ , 144–165.
  * [19] López-Cheda, A., Jácome, M. A., and Cao, R. (2017).  Nonparametric latency estimation for mixture cure models.  _Test_ 26, 2, 353–376.
  * [20] Lu, W. (2008).  Maximum likelihood estimation in the proportional hazards cure model.  _Ann. I. Stat. Math._ 60, 3, 545–574.
  * [21] Lu, W. (2010).  Efficient estimation for an accelerated failure time model with a cure fraction.  _Stat. Sin._  _20_ , 661.
  * [22] Meeker, W. Q. (1987).  Limited failure population life tests: application to integrated circuit reliability.  _Technometrics_ 29, 1, 51–65.
  * [23] Müller, U. U. and Van Keilegom, I. (2019).  Goodness-of-fit tests for the cure rate in a mixture cure model.  _Biometrika_ 106, 1, 211–227.
  * [24] Musta, E., Patilea, V., and Van Keilegom, I. (2022).  A presmoothing approach for estimation in semiparametric cox mixture cure models.  _Bernoulli (to appear)_.
  * [25] Othus, M., Barlogie, B., LeBlanc, M. L., and Crowley, J. J. (2012).  Cure models as a useful statistical tool for analyzing survival.  _Clin. Cancer Res._ 18, 14, 3731–3736.
  * [26] Parsa, M. and Van Keilegom, I. (2022).  Accelerated failure time vs cox proportional hazards mixture cure models: David vs Goliath?  _Stat. Papers (under revision)_.
  * [27] Patilea, V. and Van Keilegom, I. (2020).  A general approach for cure models in survival analysis.  _Ann. Stat._ 48, 4, 2323–2346.
  * [28] Pedrosa-Laza, M., López-Cheda, A., and Cao, R. (2022).  Cure models to estimate time until hospitalization due to covid-19.  _Appl. Intell._ 52, 1, 794–807.
  * [29] Peláez, R., Cao, R., and Vilar, J. M. (2021).  Nonparametric estimation of the probability of default with double smoothing.  _SORT-Statistics and Operations Research Transactions_ , 93–120.
  * [30] Peláez Suárez, R., Cao Abad, R., and Vilar Fernández, J. M. (2021).  Probability of default estimation in credit risk using a nonparametric approach.  _Test_ 30, 2, 383–405.
  * [31] Peng, Y. and Dear, K. B. (2000).  A nonparametric mixture model for cure rate estimation.  _Biometrics_ 56, 1, 237–243.
  * [32] Peng, Y. and Yu, B. (2021).  _Cure Models: Methods, Applications, and Implementation_.  Chapman and Hall/CRC.
  * [33] Stringer, S., Denys, D., Kahn, R. S., and Derks, E. M. (2016).  What cure models can teach us about genome-wide survival analysis.  _Behav. Genet._ 46, 2, 269–280.
  * [34] Strzalkowska-Kominiak, E. and Cao, R. (2014).  Beran-based approach for single-index models under censoring.  _Comp. Stat._ 29, 5, 1243–1261.
  * [35] Sy, J. P. and Taylor, J. M. (2000).  Estimation in a cox proportional hazards cure model.  _Biometrics_ 56, 1, 227–236.
  * [36] van der Vaart, A. W. and Wellner, J. A. (1996).  _Weak convergence and empirical processes_.  Springer Series in Statistics. Springer-Verlag, New York.  With applications to statistics.
  * [37] Van Geloven, N., Van der Veen, F., Bossuyt, P., Hompes, P., Zwinderman, A., and Mol, B. (2013).  Can we distinguish between infertility and subfertility when predicting natural conception in couples with an unfulfilled child wish?  _Hum. Reprod._ 28, 3, 658–665.
  * [38] Van Keilegom, I. and Akritas, M. G. (1999).  Transfer of tail information in censored regression models.  _Ann. Stat._  _27_ , 1745–1784.
  * [39] Wycinka, E. and Jurkiewicz, T. (2017).  Mixture cure models in prediction of time to default: comparison with logit and cox models.  In _Contemporary Trends and Challenges in Finance_. Springer, 221–231.
  * [40] Xu, J. and Peng, Y. (2014).  Nonparametric cure rate estimation with covariates.  _Can. J. Stat._ 42, 1, 1–17.
  * [41] Yamaguchi, K. (1992).  Accelerated failure-time regression models with a regression model of surviving fraction: an application to the analysis of “permanent employment” in Japan.  _J. Am. Stat. Assoc._ 87, 418, 284–292.
  * [42] Yilmaz, Y. E., Lawless, J. F., Andrulis, I. L., and Bull, S. B. (2013).  Insights from mixture cure modeling of molecular markers for prognosis in breast cancer.  _J. Clin. Oncol._ 31, 16, 2047–2054.
  * [43] Zhang, J. and Peng, Y. (2007).  A new estimation method for the semiparametric accelerated failure time mixture cure model.  _Stat. Med._ 26, 16, 3157–3171.