arXiv:1701.03769v1 [] 13 Jan 2017

A General Approach for Cure Models
in Survival Analysis
Valentin Patilea

∗

Ingrid Van Keilegom

§

January 16, 2017

Abstract
In survival analysis it often happens that some subjects under study do not experience the event of interest; they are considered to be ‘cured’. The population is thus
a mixture of two subpopulations : the one of cured subjects, and the one of ‘susceptible’ subjects. When covariates are present, a so-called mixture cure model can be
used to model the conditional survival function of the population. It depends on two
components : the probability of being cured and the conditional survival function of
the susceptible subjects.
In this paper we propose a novel approach to estimate a mixture cure model when
the data are subject to random right censoring. We work with a parametric model for
the cure proportion (like e.g. a logistic model), while the conditional survival function
of the uncured subjects is unspecified. The approach is based on an inversion which
allows to write the survival function as a function of the distribution of the observable
random variables. This leads to a very general class of models, which allows a flexible
and rich modeling of the conditional survival function. We show the identifiability of
the proposed model, as well as the weak consistency and the asymptotic normality of
the model parameters. We also consider in more detail the case where kernel estimators
are used for the nonparametric part of the model. The new estimators are compared
with the estimators from a Cox mixture cure model via finite sample simulations.
Finally, we apply the new model and estimation procedure on two medical data sets.

Key Words: Asymptotic normality; bootstrap; kernel smoothing; logistic regression; mixture
cure model; semiparametric model.
MSC2010: 62N01, 62N02, 62E20, 62F12, 62G05.

∗

CREST (Ensai), France. V. Patilea acknowledges support from the research program New Challenges
for New Data of LCL and Genes. Email address: patilea@ensai.fr.
§
ORSTAT, Katholieke Universiteit Leuven, Belgium. I. Van Keilegom acknowledges support from the
European Research Council (2016-2021, Horizon 2020 / ERC grant agreement No. 694409), and from IAP
Research Network P7/06 of the Belgian State. Email address: ingrid.vankeilegom@kuleuven.be.

1

1

Introduction

Driven by emerging applications, over the last two decades there has been an increasing
interest for time-to-event analysis models allowing the situation where a fraction of the right
censored observed lifetimes corresponds to subjects who will never experience the event. In
biostatistics such models including covariates are usually called cure models and they allow
for a positive cure fraction that corresponds to the proportion of patients cured of their
disease. For a review of these models in survival analysis, see for instance Maller & Zhou
(2001) or Peng & Taylor (2014). Economists sometimes call such models split population
models (see Schmidt & Witte 1989), while the reliability engineers refer to them as limitedfailure population life models (Meeker 1987).
At first sight, a cure regression model is nothing but a binary outcome, cured versus
uncured, regression problem. The difficulty comes from the fact that the cured subjects are
unlabeled observations among the censored data. Then one has to use all the observations,
censored and uncensored, to complete the missing information and thus to identify, estimate
and make inference on the cure fraction regression function. We propose a general approach
for this task, a tool that provides a general ground for cure regression models. The idea is
to start from the laws of the observed variables and to express the quantities of interest,
such as the cure rate and the conditional survival of the uncured subjects, as functionals of
these laws. These general expressions, that we call inversion formulae and we derive with
no particular constraint on the space of the covariates, are the vehicles that allow for a wide
modeling choice, parametric, semiparametric and nonparametric, for both the law of the
lifetime of interest and the cure rate. Indeed, the inversion formulae allow to express the
likelihood of the binary outcome model as a function of the laws of the observed variables.
The likelihood estimator of the parameter vector of the cure fraction function is then simply
the maximizer of the likelihood obtained by replacing the laws of the observations by some
estimators. With at hand the estimate of the parameter of the cure fraction, the inversion
formulae will provide an estimate for the conditional survival of the uncured subjects. For
the sake of clarity, we focus on the so-called mixture cure models with a parametric cure
fraction function, the type of model that is most popular among practitioners. Meanwhile,
the lifetime of interest is left unspecified.
The paper is organized as follows. In Section 2 we provide a general description of the
mixture cure models and next we introduce the needed notation and present the inversion
formulae on which our approach is built. We finish Section 2 by a discussion of the identification issue and some new insight on the existing approaches in the literature on cure
models. Section 3 introduces the general maximum likelihood estimator, while in Section 4
we derive the general asymptotic results. A simple bootstrap procedure for making feasible
inference is proposed. Section 4 ends with an illustration of the general approach in the case
where the conditional law of the observations is estimated by kernel smoothing. In Sections
2

5 and 6 we report some empirical results obtained with simulated and two real data sets.
Our estimator performs well in simulations and provides similar or more interpretable results
in applications compared with a competing logistic/proportional hazards mixture approach.
The technical proofs are relegated to the Appendix.

2
2.1

The model
A general class of mixture cure models

Let T denote (a possible monotone transformation of) the lifetime of interest that takes values
in (−∞, ∞]. A cured observation corresponds to the event {T = ∞}, and in the following
this event is allowed to have a positive probability. Let X be a covariate vector with support
X belonging to a general covariate space. The covariate vector could include discrete and
continuous components. The survival function FT ((t, ∞] | x) = P(T > t | X = x) for t ∈ R
and x ∈ X can be written as
FT ((t, ∞] | x) = 1 − φ(x) + φ(x)FT,0 ((t, ∞) | x),
where φ(x) = P(T < ∞ | X = x) and FT,0 ((t, ∞) | x) = P(T > t | X = x, T < ∞).
Depending on which model is used for φ(x) and FT,0 (· | x), one obtains a parametric,
semiparametric or nonparametric model, called a ‘mixture cure model’. In the literature, one
often assumes that φ(x) follows a logistic model, i.e. φ(x) = exp(a + x⊤ b)/[1 + exp(a + x⊤ b)]
for some (a, b⊤ )⊤ ∈ Rd+1 . Recently, semiparametric models (like a single-index model as in
Amico et al. 2017) or nonparametric models (as in Xu & Peng 2014 or López-Cheda et al.
2017) have been proposed. As for the survival function FT,0 (· | x) of the susceptible subjects,
a variety of models have been proposed, including parametric models (see e.g. Boag 1949,
Farewell 1982), semiparametric models based on a proportional hazards assumption (see e.g.
Kuk & Chen 1992, Sy & Taylor 2000, Fang et al. 2005, Lu 2008; see also Othus et al. 2009)
or nonparametric models (see e.g. Taylor 1995, Xu & Peng 2014).
In this paper we propose to model φ(x) parametrically, i.e. we assume that φ(·) belongs
to the family of conditional probability functions
{φ(·, β) : β ∈ B},
where φ(·, β) takes values in the interval (0, 1), β is the parameter vector of the model and B
is the parameter set. This family could be the logistic family or any other parametric family.
For the survival function FT,0 (· | x) we do not impose any assumptions in order to have a
flexible and rich class of models for FT (· | x) to choose from. Later on we will see that for the
estimation of FT,0 (· | x) any estimator that satisfies certain minimal conditions can be used,
and hence we allow for a large variety of parametric, semiparametric and nonparametric
estimation methods.
3

As is often the case with time-to-event data, we assume that the lifetime T is subject to
random right censoring, i.e. instead of observing T , we only observe the pair (Y, δ), where
Y = T ∧ C, δ = 1{T ≤ C} and C is a non-negative random variable, called the censoring
time. Some identification assumptions are required to be able to identify the conditional law
of T from the observed variables Y and δ. Let us assume that
C⊥T |X

P(C < ∞) = 1.

and

(2.1)

The conditional independence between T and C is an usual identification assumption in
survival analysis in the presence of covariates. The zero probability at infinity condition for
C implies that P(C < ∞ | X) = 1 almost surely (a.s.). This latter mild condition is required
if we admit that the observations Y are finite, which is the case in the common applications.
For the sake of simplicity, let us also consider the condition
P(T = C) = 0,

(2.2)

which is commonly used in survival analysis, and which implies that P(T = C | X) = 0 a.s.

2.2

Some notations and preliminaries

We start with some preliminary arguments, which are valid in general without assuming any
model on the functions φ, FT and FC .
The observations are characterized by the conditional sub-probabilities
H1 ((−∞, t] | x) = P(Y ≤ t, δ = 1 | X = x)
H0 ((−∞, t] | x) = P(Y ≤ t, δ = 0 | X = x),

t ∈ R, x ∈ X .

def

Then H((−∞, t] | x) = P(Y ≤ t | X = x) = H0 ((−∞, t] | x) + H1 ((−∞, t] | x). Since we
assume that Y is finite, we have
H((−∞, ∞) | x) = 1,

∀x ∈ X .

(2.3)

For j ∈ {0, 1} and x ∈ X , let τHj (x) = sup{t : Hj ([t, ∞) | x) > 0} denote the right endpoint
of the support of the conditional sub-probability Hj . Let us define τH (x) in a similar way
and note that τH (x) = max{τH0 (x), τH1 (x)}. Note that τH0 (x), τH1 (x) and τH (x) can equal
infinity, even though Y only takes finite values. For −∞ < t ≤ ∞, let us define the
conditional probabilities
FC ((−∞, t] | x) = P(C ≤ t | X = x)

and

FT ((−∞, t] | x) = P(T ≤ t | X = x),

x ∈ X.

Let us show how the probability of being cured could be identified from the observations
without any reference to a model for this probability. Under conditions (2.1)-(2.2) we can
write
H1 (dt | x) = FC ([t, ∞) | x)FT (dt | x),

H0 (dt | x) = FT ([t, ∞] | x)FC (dt | x),
4

and H([t, ∞) | x) = FT ([t, ∞] | x)FC ([t, ∞) | x). These equations could be solved and thus
they allow to express the functions FT (· | x) and FC (· | x) in an unique way as explicit
transformations of the functions H0 (· | x) and H1 (· | x). For this purpose, let us consider the
conditional cumulative hazard measures
ΛT (dt | x) =

FT (dt | x)
FT ([t, ∞] | x)

ΛC (dt | x) =

and

FC (dt | x)
,
FC ([t, ∞) | x)

x ∈ X.

The model equations yield
ΛT (dt | x) =

H1 (dt | x)
H([t, ∞) | x)

and

ΛC (dt | x) =

H0 (dt | x)
.
H([t, ∞) | x)

Then, we can write the following functionals of H0 (· | x) and H1 (· | x) :
Y
FT ((t, ∞] | x) =
{1 − ΛT (ds | x)},

(2.4)

−∞<s≤t

FC ((t, ∞) | x) =

Y

{1 − ΛC (ds | x)},

t ∈ R,

(2.5)

−∞<s≤t

Q
where s∈A stands for the product-integral over the set A (see Gill and Johansen 1990).
Moreover, if τH1 (x) < ∞, then
Y
P(T > τH1 (x) | x) =
{1 − ΛT (dt | x)},
t∈(−∞,τH1 (x)]

but there is no way to identify the conditional law of T beyond τH1 (x). Therefore, we will
impose
P(T > τH1 (x) | x) = P(T = ∞ | x),
(2.6)
Q
Q
i.e. t∈R {1 − ΛT (dt | x)} = −∞<t≤τH (x) {1 − ΛT (dt | x)}. Note that if τH1 (x) = ∞,
1
condition (2.6) is no longer an identification restriction, but just a simple consequence of the
definition of ΛT (· | x). Finally, the condition that P(C < ∞) = 1 in (2.1) can be re-expressed
by saying that we assume that H0 (· | x) and H1 (· | x) are such that
Y
P(C = ∞ | x) =
{1 − ΛC (dt | x)} = 0,
∀x ∈ X .
(2.7)
t∈R

Let us point out that this condition is satisfied only if τH1 (x) ≤ τH0 (x). Indeed, if τH1 (x) >
τH0 (x) then necessarily τH0 (x) < τH (x) and so H([τH0 (x), ∞) | x) > 0. Hence, ΛC (R | x) =
ΛC ((−∞, τH0 (x)] | x) < ∞, and thus P(C = ∞ | x) > 0, which contradicts (2.7).
It is important to understand that any two conditional sub-probabilities H0 (· | x) and
H1 (· | x) satisfying conditions (2.1) and (2.2) define uniquely FT (· | x) and FC (· | x). Indeed,
FT (· | x) is precisely the probability distribution of T given X = x with all the mass beyond
τH1 (x) concentrated at infinity. In general, FT (· | x) and FC (· | x) are only functionals of
H0 (· | x) and H1 (· | x).
We will assume conditions (2.1), (2.2) and (2.6) throughout the paper.
5

2.3

A key point for the new approach: the inversion formulae

Write
H([t, ∞) | x) = FT ([t, ∞] | x)FC ([t, ∞) | x)
= FT ([t, ∞) | x)FC ([t, ∞) | x) + P(T = ∞ | x)FC ([t, ∞) | x),
and thus
FT ([t, ∞) | x) =

H([t, ∞) | x) − P(T = ∞ | x)FC ([t, ∞) | x)
.
FC ([t, ∞) | x)

(2.8)

Consider the conditional cumulative hazard measure for the finite values of the lifetime of
interest:
FT (dt | x)
FT,0 (dt | x)
def
=
ΛT,0 (dt | x) =
FT,0 ([t, ∞) | x)
FT ([t, ∞) | x)
for t ∈ R. Since H1 (dt | x) = FC ([t, ∞) | x)FT (dt | x), using relationship (2.8) we obtain
ΛT,0(dt | x) =

H1 (dt | x)
.
H([t, ∞) | x) − P(T = ∞ | x)FC ([t, ∞) | x)

Next, using the product-integral we can write
Y
FT,0 ((t, ∞) | x) =
{1 − ΛT,0 (ds | x)},

t ∈ R, x ∈ X .

(2.9)

(2.10)

−∞<s≤t

Let us recall that FC (· | x) can be written as a transformation of H0 (· | x) and H1 (· | x),
see equations (2.4) and (2.5). This representation is not surprising since we can consider C
as a lifetime of interest and hence T plays the role of a censoring variable. Hence, estimating
the conditional distribution function FC (· | x) should not be more complicated than in a
classical conditional Kaplan-Meier setup, since the fact that T could be equal to infinity
with positive conditional probability is irrelevant when estimating FC (· | x).
Finally, the representation of FC (· | x) given in equation (2.5), plugged into equation
(2.9), allows to express ΛT,0 (· | x), and thus FT,0 (· | x), as maps of P(T = ∞ | x) and the
measures H0 (· | x) and H1 (· | x). This will be the key element for providing more insight in
the existing approaches and the starting point of our new approach.

2.4

Model identification issues

Let us now investigate the identification issue. Recall that our model involves the functions
FT (· | x), FC (· | x) and φ(·, β), and the assumptions (2.1), (2.2) and (2.6). For a fixed value
of the parameter β, and for t ∈ R and x ∈ X , let
ΛβT,0 (dt | x) =

H1 (dt | x)
,
H([t, ∞) | x) − [1 − φ(x, β)]FC ([t, ∞) | x)
6

(2.11)

and
β
FT,0
((t, ∞) | x) =

Y n

−∞<s≤t

o
1 − ΛβT,0 (ds | x) .

(2.12)

Let FY,δ (·, · | x) denote the conditional law of (Y, δ) given X = x. Moreover, let
β
β
FY,δ
(dt, 1 | x) = φ(x, β)FC ((t, ∞) | x)FT,0
(dt | x)

and
β
β
FY,δ
(dt, 0 | x) = [FT,0
((t, ∞) | x)φ(x, β) + 1 − φ(x, β)]FC (dt | x).

These equations define a conditional law for the observations (Y, δ) based on the model.
More precisely, for a choice of FT (· | x), FC (· | x) and β, the model yields a conditional law
for (Y, δ) given X = x. If the model is correctly specified, there exists a value β0 such that
β0
FY,δ (·, · | x) = FY,δ
(·, · | x),

∀x ∈ X .

(2.13)

The remaining question is whether the true value of the parameter is identifiable. In other
words, one should check if, given the conditional subdistributions H0 (· | x) and H1 (· | x),
x ∈ X , there exists a unique β0 satisfying condition (2.13). For this purpose we impose the
following mild condition:

and we show that

e almost surely
φ(X, β) = φ(X, β),

e
⇒ β = β,

(2.14)

βe
β0
e ∀x ∈ X .
(·, · | x), ∀x ∈ X ⇒ φ(x, β0 ) = φ(x, β),
FY,δ
(·, · | x) = FY,δ
e

β0
β
Indeed, if FY,δ
(·, · | x) = FY,δ
(·, · | x), then for any x,

β0
e C ((t, ∞) | x)F βe (dt | x),
φ(x, β0 )FC ((t, ∞) | x)FT,0
(dt | x) = φ(x, β)F
T,0

for all t ∈ (−∞, τH (x)] ∩ R. Our condition (2.7) guarantees that τH (x) = τH0 (x) so that
FC ((t, ∞) | x) should be necessarily positive for t ∈ (−∞, τH (x)), and thus could be simplified in the last display. Deduce that
β0
e βe (dt | x),
φ(x, β0 )FT,0
(dt | x) = φ(x, β)F
T,0

for all t ∈ (−∞, τH (x)). Finally, recall that by construction, in the model we consider,
β
β
FT,0
((−∞, ∞) | x) = 1 for any β such that FY,δ
(·, · | x) coincides with the conditional law of
(Y, δ) given X = x, for any x. Thus taking integrals on (−∞, ∞) on both sides of the last
e Let us gather these facts in the following statement.
display we obtain φ(·, β0 ) = φ(·, β).
Theorem 2.1. Under conditions (2.1), (2.2), (2.6) and (2.14) the model is identifiable.
7

2.5

Interpreting the previous modeling approaches

We suppose here that the function φ(x) follows a logistic model, and comment on several
models for FT,0 that have been considered in the literature.
2.5.1

Parametric and proportional hazards mixture model

In a parametric modeling, one usually supposes that τH0 (x) = τH1 (x) = ∞ and that ΛT,0 (· |
x) belongs to a parametric family of cumulative hazard functions, like for instance the Weibull
model; see Farewell (1982).
Several contributions proposed a more flexible semiparametric proportional hazards (PH)
approach; see Fang et al. (2005), Lu (2008) and the references therein. In such a model one
imposes a PH structure for the ΛT,0 (· | x) measure. More precisely, it is supposed that
ΛT,0(dt | x) =

H1 (dt | x)
= exp(x⊤ γ)Λ0 (dt),
H([t, ∞) | x) − [1 − φ(x, β)]FC ([t, ∞) | x)

where γ is some parameter to be estimated and Λ0 (·) is an unknown baseline cumulative
hazard function. Our inversion formulae reveal that in this approach the parameters γ and
Λ0 depend on the observed conditional measures H0 (· | x) and H1 (· | x), but also on the
parameter β. The same is true for the parametric models.
2.5.2

Kaplan-Meier mixture cure model

Taylor (1995) suggested to estimate FT,0 using a Kaplan-Meier type estimator. With such an
approach one implicitly assumes that the law of T given X and given that T < ∞ does not
depend on X. This is equivalent to supposing that ΛT,0(· | x) = ΛT,0(·). Next, to estimate
ΛT,0 (·) one has to modify the unconditional version of the usual inversion formulae (2.4)
to take into account the conditional probability of the event {T = ∞}. Following Taylor’s
approach we rewrite (2.9) as
ΛT,0(dt) =

H1 (dt | x)
o
n
.
1−φ(x,β)
H1 ([t, ∞) | x)+ [t,∞) 1 − φ(x,β)FT,0 ([s,∞))+1−φ(x,β) H0 (ds | x)
R

Next, assume that the last equality remains true if H0 (dt | x) and H1 (dt | x) are replaced
by their unconditional versions, that is assume that
ΛT,0(dt) =

n
H1 ([t, ∞)) + [t,∞) 1 −
R

H1 (dt)
1−φ(x,β)
φ(x,β)FT,0 ([s,∞))+1−φ(x,β)

o

.

(2.15)

H0 (ds)

See equations (2) and (3) in Taylor (1995). The equation above could be solved iteratively
(m)
(m+1)
by a EM-type procedure: for a given β and an iteration FT,0 (·), build ΛT,0 (dt) and the
(m+1)
updated estimate FT,0 (·); see Taylor (1995) for the details. Let us point out that even if
8

(T, C) is independent of X and thus H1 (· | x) does not depend on x, the subdistribution
H0 (· | x) still depends on x, since
H0 (dt | x) = FT ((t, ∞] | x)FC (dt | x)
= [FT ((t, ∞) | x) + P(T = ∞ | x)]FC (dt | x).
Hence, a more natural form of equation (2.15) is
ΛT,0 (dt) =

n
H1 ([t, ∞)) + [t,∞) 1 −
R

H1 (dt)
1−φ(x,β)
φ(x,β)FT,0 ([s,∞))+(1−φ(x,β))

o

.
H0 (ds | x)

The investigation of a EM-type procedure based on the latter equation will be considered
elsewhere.

3

Maximum likelihood estimation

Let (Yi , δi , Xi ) (i = 1, . . . , n) be a sample of n i.i.d. copies of the vector (Y, δ, X).
We use a likelihood approach based on formulae (2.9) and (2.5) to build an estimator of
β
b k (· | x) of the subdistriφ(·, β) and FT,0
(· | x). To build the likelihood we use estimates H
butions Hk (· | x), k ∈ {0, 1}. These estimates are constructed with the sample of (Y, δ, X),
without reference to any model for the conditional probability P(T < ∞ | x). At this stage
b k (· | x). To derive the asymptotic results
it is not necessary to impose a particular form for H
β
we will only impose that these estimators satisfy some mild conditions. Let FbT,0
(· | x) be
b 0 (· | x) and H
b 1 (· | x) instead of H0 (· | x)
defined as in equations (2.11) and (2.12) with H
and H1 (· | x), that is
)
(
Y
b
H
(ds
|
x)
1
β
,
FbT,0
((t, ∞) | x) =
1−
b
bC ([s, ∞) | x)
H([s,
∞)
|
x)
−
[1
−
φ(x,
β)]
F
−∞<s≤t
b 0 (· | x)
where FbC (· | x) is the estimator obtained from equations (2.4) and (2.5) but with H
b 1 (· | x), i.e.
and H
(
)
Y
b 0 (ds | x)
H
FbC ((t, ∞) | x) =
1−
.
b
H([s,
∞)
|
x)
−∞<s≤t

Let fX denote the density of the covariate vector with respect to some dominating measure. The contribution of the observation (Yi , δi , Xi ) to the likelihood when δi = 1 is then
β
φ(Xi , β)FbT,0
({Yi } | Xi )fX (Xi )FbC ((Yi , ∞) | Xi ),

while the contribution when δi = 0 is

β
FbC ({Yi } | Xi )fX (Xi )[φ(Xi , β)FbT,0
((Yi , ∞) | Xi ) + 1 − φ(Xi , β)].

9

Since the laws of the censoring variable and of the covariate vector do not carry information of
the parameter β, we can drop the factors FbC ((Yi , ∞) | Xi )fX (Xi ) and FbC ({Yi } | Xi )fX (Xi ).
bn (β) where
Hence the criterion to be maximized with respect to β ∈ B is L
bn (β) =
L

n n
oδi n
o1−δi
Y
β
β
φ(Xi , β)FbT,0
({Yi } | Xi )
φ(Xi , β)FbT,0
((Yi , ∞) | Xi ) + 1 − φ(Xi , β)
.
i=1

The estimator we propose is

bn (β).
βb = arg max log L
β∈B

(3.1)

Let us review the identification issue in the context of the likelihood estimation approach.
If conditions (2.1), (2.2) and (2.6) hold true and the parametric model for the conditional
probability of the event {T = ∞} is correct and identifiable in the sense of condition (2.14),
the true parameter β0 is the value identified by condition (2.13). This is the conclusion of
Theorem 2.1 above. It remains to check that the proposed likelihood approach allows to
consistently estimate β0 .
Let
log p(t, δ, x; β) = δ log p1 (t, x; β) + (1 − δ) log p0 (t, x; β),
(3.2)
with
p1 (t, x; β) =

β
φ(x, β)FT,0
(dt | x)FC ((t, ∞) | x)
,
H1 (dt | x)

β
FC (dt | x)[φ(x, β)FT,0
((t, ∞) | x) + 1 − φ(x, β)]
p0 (t, x; β) =
.
H0 (dt | x)

(3.3)
(3.4)

Following a common notational convention, see for instance Gill (1994), here we treat dt not
just as the length of a small time interval [t, t + dt) but also as the name of the interval
itself. Moreover, we use the convention 0/0 = 1. Let us notice that, up to additive terms
not containing β, the function β 7→ E[log p(Y, δ, X; β)] is expected to be the limit of the
bn (·). Hence, a minimal condition for guaranteeing the consistency
random function log L
of the likelihood estimation approach is that β0 is the maximizer of the limit likelihood
criterion E[log p(Y, δ, X; ·)]. This is proved in the following proposition using a Bernoulli
sample likelihood inequality. The proof is given in the Appendix.
Proposition 3.1. Suppose that conditions (2.1), (2.2), (2.6) and (2.14) hold true. If β0 is
the value of the parameter defined by equation (2.13), then for any β 6= β0 ,
E [log p(Y, δ, X; β)] < E [log p(Y, δ, X; β0 )] .

4

General asymptotic results

Little assumptions were needed for our analysis so far. To proceed further with the asymptotic results we need to be more specific with respect to several aspects. In order to prove
10

bn (β) along sequences of values
consistency, we have to control the asymptotic behavior of L
of the parameter β. Such a control requires a control of denominators like
b
H([t,
∞) | x) − (1 − φ(x, β))FbC ([t, ∞) | x)

on the support of H1 (· | x), uniformly with respect to x. A usual way to deal with this
technical difficulty is to consider a finite threshold τ beyond which no uncensored lifetime is
observed, i.e.
inf H1 ((−∞, τ ] | x) = 1.
(4.1)
x

Moreover, to be able to keep denominators away from zero, we require the condition
inf H1 ({τ } | x)H0 ((τ, ∞) | x) > 0.

x∈X

(4.2)

In particular, this condition implies
τH0 (x) > τH1 (x) = τ,

x ∈ X.

Moreover, given condition (2.2), necessarily H0 ({τ } | x) = 0, ∀x. This means that FC ({τ } |
x) = 0, ∀x. This constraint on H0 ({τ } | x) could be relaxed at the expense of suitable
adjustments of the inversion formulae. For simplicity, we keep condition (2.2). Let us also
β
notice that condition (4.2) implies inf x FC ((τ, ∞) | x) > 0, and inf x FT,0
({τ } | x) > 0, ∀β.
Conditions like in equations (4.1)-(4.2) are more or less explicitly used in the literature
of cure models. Sometimes τ is justified as representing a total follow-up of the study. For
instance, Lu (2008) supposes that Y = min{T, min(C, τ )} and δ = 1{T ≤ min(C, τ )}, where
T = ηT ∗ +(1−η)∞, with T ∗ < ∞ and η ∈ {0, 1}. The conditional probability of being cured
is precisely the conditional probability of the event {η = 0}. Next, Lu (2008) supposes that
inf x P(τ ≤ T ≤ C | x) > 0, and Λ0 (τ ) < ∞, where Λ0 (·) is the cumulative hazard function
of T ∗ . All these conditions together clearly imply our conditions (4.1)-(4.2).
Fang et al. (2005) implicitly restrict the uncensored lifetimes to some compact interval
[0, τ ] and suppose E(δ1{Y ≥ τ }) > 0. This could be possible only if H1 ({τ } | x) > 0 for a set
of values x with positive probability. In a proportional hazards context with the covariates
taking values in a bounded set, as assumed by Fang et al. (2005), this is equivalent to
H1 ({τ } | x) ≥ c > 0 for almost all x, for some constant c.
The fact that technical conditions similar to our conditions (4.1)-(4.2) could be traced in
the cure models literature is not unexpected in view of our Section 2.5. Indeed, the existing
approaches could be interpreted through our inversion formulae and thus the technical problems we face in the asymptotic investigation are expected to be also present in the alternative
approaches.

11

4.1

Consistency

Let us sketch the arguments we use in the proof of Theorem 4.1 below for deriving the
b On one hand, if the conditional subdistributions Hk (· | x) are given, one
consistency of β.
can build the purely parametric likelihood
Ln (β) =

n n
Y

β
φ(Xi , β)FT,0
({Yi } | Xi )

i=1

oδin
o1−δi
β
φ(Xi , β)FT,0
((Yi , ∞) | Xi )+1−φ(Xi , β)
. (4.3)

bn (β)
By construction, Ln (β) is a functional of H0 (· | x) and H1 (· | x), x ∈ X , while L
is a functional of the estimated versions of H0 (· | x) and H1 (· | x). Hence, a prerequisite
condition for deriving the consistency of our semiparametric estimator βb is the consistency
of the infeasible maximum likelihood estimator
βe = arg max log Ln (β).
β∈B

A necessary condition for the consistency of βb is

bn (β) − log Ln (β)| = oP (1).
sup | log L

(4.4)

β∈B

We then have that

b ≥ log Ln (β)
e − oP (1).
log Ln (β)

From this we will derive the consistency of βb using Section 5.2 in van der Vaart (1998). See
the proof in the Appendix for details.
b k − Hk as
To prove condition (4.4), we have to guarantee the uniform convergence of H
stated in Assumption (AC1) below. Indeed, this uniform convergence will imply
sup sup |FbC ([t, ∞) | x) − FC ([t, ∞) | x)| = oP (1),

(4.5)

x∈X t∈(−∞,τ ]

and

sup sup

β
β
sup |FbT,0
([t, ∞) | x) − FT,0
([t, ∞) | x)| = oP (1).

(4.6)

β∈B x∈X t∈(−∞,τ ]

See Lemma 7.1 in the Appendix. The uniform convergence in equation (4.4) then follows.
b we need the following assumptions :
To prove the consistency of β,

(AC1) For τ appearing in conditions (4.1)-(4.2),

b k ([t, ∞) | x) − Hk ([t, ∞) | x)| = oP (1),
sup sup |H
x∈X t∈(−∞,τ ]

(AC2) The parameter set B ⊂ Rp is compact.

12

k ∈ {0, 1}.

(AC3) There exist some constants a > 0 and c1 > 0 such that
|φ(x, β) − φ(x, β ′ )| ≤ c1 kβ − β ′ ka ,

∀β, β ′ ∈ B, ∀x, x′ ∈ X .

(AC4) inf β∈B inf x∈X φ(x, β) > 0.
Now we can state our consistency result.
Theorem 4.1. Assume that (AC1)-(AC4) and (2.1), (2.2), (2.6), (2.14), (4.1) and (4.2)
hold true. Moreover, assume that there exists a unique value β0 in the parameter space B
such that (2.13) is true. Then, βb − β0 = oP (1).

Let us point out that the consistency result is stated in terms of the subdistributions of
the observations and the conditional probability model {φ(x, ·) : β ∈ B}. If the identification
b x)
assumptions used in Proposition 3.1 hold true and the model is correctly specified, φ(β,
consistently estimates the cure probability P(T = ∞ | x) for all x in the support of X. Let
us also notice that condition (AC3) guarantees the Glivenko-Cantelli property for certain
classes of functions. It could be significantly weakened, but in the applications our condition
(AC3) will cover the common modeling situations. Condition (AC4) is a weak condition on
the model φ(x, β) and is e.g. satisfied for the logistic model if X and B are compact.

4.2

Asymptotic normality

For the asymptotic normality we will use the approach in Chen et al. (2003). For this
bn (β) with respect to β.
purpose we use the derivative of log L
bn (β) with respect
First note that the vector of partial derivatives of the log-likelihood log L
to the components of β equals
1X n
b 1 ({Yi } | Xi )]
b
∇β δi log φ(Xi , β) + δi log[H
∇β log Ln (β) =
n i=1
h
i
b
b
−δi log H([Yi , ∞) | Xi ) − (1 − φ(Xi , β))FC ([Yi , ∞) | Xi )
n

(1)
+δi log FbT,β ([Yi , ∞) | Xi )
h
io
(1)
+(1 − δi ) log φ(Xi , β)Fb ([Yi , ∞) | Xi ) + 1 − φ(Xi , β) .
T,β

=

1
n

n
X
i=1

b0, H
b 1 ),
∇β qi (β, H

where qi is defined in the proof of Theorem 4.1.
b we embed the nuisance functions
To develop the asymptotic normality of our estimator β,
Hk ([·, ∞) | ·) (k = 0, 1) in a functional space H, which is equipped with a pseudo-norm k·kH .
Both the space H and its pseudo-norm k · kH will be chosen depending on the estimators
13

b k ([·, ∞) | ·), k = 0, 1, and have to satisfy certain conditions, which we give below. The
H
true vector of nuisance functions is
η0 (t, x) = (η01 (t, x), η02 (t, x)) = (H0 ([t, ∞) | x), H1 ([t, ∞) | x)).
For each x ∈ X and for each η ∈ H, let η(dt, x) = (η1 (dt, x), η2 (dt, x)) be the measures
associated to the non-increasing functions η(·, x) = (η1 (·, x), η2 (·, x)), and define
n

1X
m(Yi , δi , Xi ; β, η)
Mn (β, η) =
n i=1

(4.7)

M(β, η) = E[m(Y, δ, X; β, η)],

(4.8)

and
where
δ∇β φ(x, β) δ∇β φ(x, β)T2 (η)(t, x)
−
+ δT4 (β, η)(t, x)
φ(x, β)
T1 (β, η)(t, x)
∇β φ(x, β)T3(β, η)(t, x) + φ(x, β)(T3 T4 )(β, η)(t, x) − ∇β φ(x, β)
,
+ (1 − δ)
φ(x, β)T3 (β, η)(t, x) + 1 − φ(x, β)

m(t, δ, x; β, η) =

where
T1 (β, η)(t, x) = [η1 (t, x) + η2 (t, x)] − (1 − φ(x, β))T2 (β, η)(t, x)

Y 
η1 (ds, x)
T2 (β, η)(t, x) ≡ T2 (η)(t, x) =
1−
(η1 + η2 )(s, x)
−∞<s<t


Y
η2 (ds, x)
T3 (β, η)(t, x) =
1−
T1 (β, η)(s, x)
−∞<s<t
T4 (β, η)(t, x) = ∇β log T3 (β, η)(t, x)
Z
= −∇β φ(x, β)

(−∞,t)

T2 (η)(s, x)η2 (ds, x)
.
T1 (β, η)(s, x)[T1(β, η)(s, x) − η2 ({s}, x)]

Note that for η = η0 , we have
T1 (β, η0 )(t, x) = H([t, ∞) | x) − (1 − φ(x, β))FC ([t, ∞) | x),
β
T2 (η0 )(t, x) = FC ([t, ∞) | x) and T3 (β, η0)(t, x) = FT,0
([t, ∞) | x). Hence, we have that

M(β0 , η0 ) = 0

and

b ηb)k = infβ∈B kMn (β, ηb)k,
kMn (β,

b 0 ([t, ∞) | x), H
b 1 ([t, ∞) | x)).
where ηb(t, x) = (b
η1 (t, x), ηb2 (t, x)) = (H
If in addition β = β0 ,
β0
T1 (β0 , η0 )(t, x) = FT,0
([t, ∞) | x)FC ([t, ∞) | x),

14

t ∈ (−∞, τ ],

(4.9)

and thus for any x, the map t 7→ T1 (β0 , η0 )(t, x) is decreasing on (−∞, τ ]. Moreover, by
condition (4.2),
inf T1 (β0 , η0 )(τ, x) > 0.
x∈X

We need this lower bound to be valid on a neighborhood around β0 . Hence, let us consider
B0 a neighborhood of β0 such that
inf inf T1 (β, η0 )(τ, x) > 0.

β∈B0 x∈X

(4.10)

The existence of B0 is guaranteed by condition (4.2) and the regularity of the function
φ(·, ·); see assumption (AN3) below that strengthens assumption (AC3). Finally, let us
note that by construction for any t ∈ (−∞, τ ), H({t} | x) = H0 ({t} | x) + H1 ({t} | x) ≥
(1 − φ(x, β))FC ({t} | x) + H1 ({t} | x) and thus
H([t, ∞) | x) − (1 − φ(x, β))FC ([t, ∞) | x) − H1 ({t} | x)
≥ H((t, ∞) | x) − (1 − φ(x, β))FC ((t, ∞) | x).
Then, by the arguments guaranteeing the existence of a set B0 as in equation (4.10),
inf inf

inf

[T1 (β, η0 )(t, x) − H1 ({t} | x)] > 0.

β∈B0 x∈X t∈(−∞,τ )

(4.11)

Further, define the Gâteaux derivative of M(β, η0 ) in the direction [η − η0 ] by
h
i
−1
∇η M(β, η0 )[η − η0 ] = lim τ
M(β, η0 + τ (η − η0 )) − M(β, η0 ) ,
τ →0

and in a similar way the Gâteaux derivatives ∇η Tj (β, η0)[η − η0 ] are defined.
We need the following assumptions :

(AN1) The matrix ∇β M(β, η) exists for β in a neighborhood B0 of β0 and is continuous in β
for β = β0 . Moreover, ∇β M(β0 , η0 ) is non-singular.
(AN2) Hk ([·, ∞) | ·) ∈ H for k = 0, 1.
(AN3) The function β → φ(x, β) is continuously differentiable for all x ∈ X , and the derivative
is bounded uniformly in x ∈ X and β ∈ B0 . Moreover, B0 is compact and β0 belongs
to the interior of B0 .
b k ([·, ∞) | ·) satisfies the following :
(AN4) For k = 0, 1, the estimator H
b k ([·, ∞) | ·) ∈ H) → 1
(i) P(H
b k − Hk )([·, ∞) | ·)kH = oP (n−1/4 )
(ii) k(H
15

(iii) There exist functions Ψ1 and Ψ2 , such that
1
X
k=0


Z
E ψ1k (Y, X)
∗

−∞<u<Y

=

1
n

n
X



b k − Hk )([u, ∞) | X)
ψ2k (u, X)d (H

Ψ1 (Yi , δi , Xi ) + R1n ,

i=1

and
1
X

k,ℓ=0

∗

E



ψ3k (Y, X)

Z

ψ4k (u, X)ψ5k

−∞<u<Y
n

=





b
(Hk − Hk )([u, ∞) | X), X dHℓ (u | X)

1X
Ψ2 (Yi , δi , Xi ) + R2n ,
n i=1

where E∗ denotes conditional expectation given the data (Yi , δi , Xi ), 1 ≤ i ≤ n,
the functions ψjk are defined in (7.13) in the Appendix, and where
E[Ψℓ (Y, δ, X)] = 0,

Rℓn = oP (n−1/2 )

(j = 1, . . . , 5; k = 0, 1; ℓ = 1, 2). Note that the above expectations are conditionally on the sample and are taken with respect to the generic variables Y, δ, X
which have the same law as the sample.
R∞p
log N(ε, H, k · kH )dε < ∞, where N(ε, H, k · kH ) is the
(AN5) The class H satisfies 0
ε-covering number of the space H with respect to the norm k · kH , i.e. the smallest
number of balls of k · kH -radius ε needed to cover the space H.
Theorem 4.2. Assume that βb − β0 = oP (1) and that (AN1)-(AN5) and (2.1), (2.2), (2.6),
(2.14), (4.1) and (4.2) hold true. Then,

n1/2 βb − β0 ⇒ N (0, Ω),

where


Ω = ∇β M(β0 , η0 )


V ∇β M(β0 , η0 )

and V = Var m(Y, δ, X; β0 , η0 ) + Ψ1 (Y, δ, X) + Ψ2 (Y, δ, X) .

4.3

−1

−1

Bootstrap consistency

Although in principle one can use Theorem 4.2 above for making inference, the asymptotic
variance Ω has a complicated structure, and the estimation of Ω would not only be cumbersome, but its precision for small samples could moreover be rather poor. We continue
16

this section by showing that a bootstrap procedure can be used to estimate the asymptotic
b to approximate the whole distribution of βb or to construct confidence intervals
variance of β,
or test hypotheses regarding β0 .
Here, we propose to use a naive bootstrap procedure, consisting in drawing triplets
∗ ∗
(Yi , δi , Xi∗ ), 1 ≤ i ≤ n, randomly with replacement from the data (Yi , δi , Xi ), 1 ≤ i ≤ n.
b ∗ be the same estimator as H
b k (k = 0, 1) but based on the bootstrap data, and for
Let H
k
P
each (β, η) let Mn∗ (β, η) = n−1 ni=1 m(Yi∗ , δi∗ , Xi∗ ; β, η). Define the bootstrap estimator βb∗
to be any sequence that satisfies
b ηb)k = infβ∈B kM ∗ (β, ηb∗) − Mn (β,
b ηb)k,
kMn∗ (βb∗ , ηb∗ ) − Mn (β,
n

b ∗ ([t, ∞) | x), H
b ∗ ([t, ∞) | x)).
where ηb∗ (t, x) = (b
η1∗ (t, x), ηb2∗ (t, x)) = (H
0
1
The following result shows that the bootstrap works, in the sense that it allows to recover
correctly the distribution of n1/2 (βb − β0 ).

Theorem 4.3. Assume that βb − β0 = oP (1) and that (AN1)-(AN5) hold true. Moreover,
assume that ∇β M(β, η) is continuous in η (with respect to k · kH ) at (β, η) = (β0 , η0 ), and
b k − Hk replaced by H
b∗ − H
b k (k = 0, 1) in P∗ -probability. Then,
that (AN4) holds true with H
k


b ≤ u − P n1/2 (βb − β0 ) ≤ u = oP (1),
sup P∗ n1/2 (βb∗ − β)

u∈Rp

where P∗ denotes probability conditionally on the data (Yi , δi , Xi ), i = 1, . . . , n, and where
the inequality sign means the component-wise inequality for vectors.

4.4

Verification of the assumptions for kernel estimators

We finish this section with an illustration of the verification of the assumptions of our
asymptotic results when the conditional subdistributions Hk are estimated by means of
kernel smoothing.
Consider the case where X is composed of continuous and discrete components, that is
X = (Xc , Xd ) ∈ Xc × Xd ⊂ Rdc × Rdd , with dc + dd = d ≥ 1. For simplicity, assume that the
support Xd of the discrete subvector Xd is finite. We also assume that the life time T has
not been transformed by a logarithmic or other transformation, so that its support is [0, ∞].
The subdistributions Hk ([t, ∞) | x) could be estimated by means of a kernel estimator :
b k ([t, ∞) | x) =
H

where for any (xc , xd ) ∈ Xc × Xd ,

n
X
i=1

e hn (Xi − x)
K
I(Yi ≥ t, δi = k),
Pn e
j=1 Khn (Xj − x)

e hn (Xi − x) = Khn (Xc,i − xc )I(Xd,i = xd ),
K
17

hn is a bandwidth sequence, Kh (·) = K(·/h)/hdc , K(u) = k(u1 ) · . . . · k(udc ) and k is a
probability density function.
Nonparametric smoothing of continuous covariates is possible for dimensions dc larger
than 1. However, the technical arguments necessary to verify the assumptions used for the
asymptotic results are tedious. Therefore, in the following we consider dc = 1. The discrete
covariates do not contribute to the curse of dimensionality, and therefore dd could be larger
than 1. However, for simplicity, below we do not consider discrete covariates.
To satisfy assumption (AN4), we need to impose the following conditions :
−1
(C1) The sequence hn satisfies nh4n → 0 and nh3+ζ
→ ∞ for some ζ > 0.
n (log n)

(C2) The support X of X is a compact subset of R.
(C3) The probability density function K has compact support,
twice continuously differentiable.

R

uK(u)du = 0 and K is

Further, let F1 be the space of functions from [0, τ ] to [0, 1] with variation bounded by
M, and let F2 be the space of continuously differentiable functions f from X to [−M, M]
that satisfy supx∈X |f ′(x)| ≤ M and supx1 ,x2 ∈X |f ′(x1 ) − f ′ (x2 )|/|x1 − x2 |ǫ ≤ M for some
M < ∞ and 0 < ǫ < 1. Let
n
∂
η(·, x) ∈ F1 for all x ∈ X ,
H = (t, x) → η(t, x) : η(·, x) ∈ F1 ,
∂x
o
and η(t, ·) ∈ F2 for all 0 ≤ t ≤ τ .
We define the following norm associated with the space H : for η ∈ H, let
kηkH = sup sup |η(t, x)|.
0≤t≤τ x∈X

Then, it follows from Propositions 1 and 2 in Akritas and Van Keilegom (2001) that
−1
b
P(Hk ∈ H) → 1 provided nh3+ζ
→ ∞, with ζ > 0 as in condition (C1). Moreover,
n (log n)
b
supx,t |Hk ([t, ∞) | x)−Hk ([t, ∞) | x)| = OP ((nhn )−1/2 (log n)1/2 ) = oP (n−1/4 ) (see Proposition
1 in Akritas and Van Keilegom 2001). The class H satisfies assumption (AN5) thanks to
Lemma 6.1 in Lopez (2011). It remains to show the validity of assumption (AN4)(iii). We
will show the first statement, the second one can be shown in a similar way. Note that the

18

left hand side equals
Z
1
h
X
E ψ1k (Y, X)


X
1
ψ2k (u, X) fX−1 (X)
Kh (Xi − X) d I(Yi ≥ u, δi = k)
nh
0<u<Y
i=1

i
−1/2
−Hk ([u, ∞) | X) Yi , δi , Xi + oP (n
)

k=0

n

n
1
n
h
X
X
1 −1
Kh (Xi − X) − ψ2k (Yi , X)I(Yi ≤ Y, δi = k)
=
E ψ1k (Y, X) fX (X)
nh
i=1
k=0
Z
o
i
−
ψ2k (u, X) dHk ([u, ∞) | X) Yi , δi , Xi + oP (n−1/2 )

=−

1
X
k=0

0<u<Y
n
X
−1

n

−

Z

i=1

h
n
E ψ1k (Y, Xi ) ψ2k (Yi , Xi )I(Yi ≤ Y, δi = k)
ψ2k (u, Xi ) dHk ((−∞, u] | Xi )

0<u<Y

o

i
Yi , δi , Xi , X = Xi + oP (n−1/2 + h2 ),

which is of the required form.

5

Simulations

In this section we will investigate the small sample performance of our estimation method.
We consider the following model. The covariate X is generated from a uniform distribution
on [−1, 1], and the conditional probability φ(x, β) of not being cured follows a logistic model :
φ(x, β) =

exp(β1 + β2 x)
,
1 + exp(β1 + β2 x)

for any −1 ≤ x ≤ 1. We will work with β0 = (β01 , β02 ) = (1.75, 2) and (1.1, 2), corresponding
to an average cure rate of 20% respectively 30%. The conditional distribution function
FT,0 (·|x) of the uncured individuals is constructed as follows. For a given X, we draw T
from an exponential distribution with mean equal to exp[−(γ0 + γ1 x + γ2 /(1 + 2x2 ))], where
γ0 = γ1 = 0.5 and γ2 ∈ {0, 1, 2}. Next, in order to respect condition (4.2), we truncate
this distribution at τ , which is the quantile of order 0.97 of an exponential distribution with
mean E{exp[−(γ0 + γ1 X + γ2 /(1 + 2X 2 ))]}, i.e.
FT,0 ([0, t]|x) = 1 − exp[− exp(γ0 + γ1 x + γ2 /(1 + 2x2 ))t]I(0 ≤ t ≤ τ ).
Note that this is the distribution function corresponding to a Cox model with baseline hazard
equal to I(0 ≤ t ≤ τ ), and exponential factor equal to exp(γ0 + γ1 x + γ2 /(1 + 2x2 )).
Next, we generate the censoring variable C independently of X from an exponential
distribution with mean equal to 1.65 when β0 = (1.75, 2), and with mean 1.45 when β0 =
(1.1, 2). In this way we have respectively 40% and 50% of censoring when γ2 = 0.
19

In what follows we will compare our estimator of β with the estimator proposed by Lu
(2008) which assumes a Cox model for the uncured individuals. The exponential factor in
the Cox model is assumed to be linear in the covariate X, and hence the Cox model will only
be verified when γ2 = 0. The estimated β coefficients under the Cox model are obtained
using the R package smcure.
For our estimation procedure we used the kernel estimators given in Section 4.4, and we
programmed βb using the optimization procedure optim in R. As starting values we used the
estimator obtained from a logistic model based on the censoring indicator (as a surrogate for
the unobserved cure indicator). However, due to the non-concavity of our likelihood function
and due to the inconsistency of this vector of starting values, the procedure optim often ends
up in a local maximum instead of the global maximum. To circumvent this problem, we
added the following intermediate step to the estimation procedure. Based on the initial
starting values, we estimate β from a logistic model based on the nonparametric estimator
P
FbT ([0, ∞)|x), so we maximize the log-likelihood ni=1 {(1 − FbT ([0, ∞)|Xi )) log(φ(Xi , β)) +
FbT ([0, ∞)|Xi ) log(1 − φ(Xi , β))}. Since this log-likelihood is concave it has a unique local
and global maximum, expected to be close to the maximizer of our likelihood. We now use
this intermediate estimate as starting value for our likelihood maximization.
The results of this two-step maximization procedure are given in Table 1 for the case
where β0 = (1.75, 2), and in Table 2 for the case where β0 = (1.1, 2). A total of 500
samples of size n = 150 and n = 300 are generated, and the tables show the bias and
mean squared error (MSE) of the estimators βb1 and βb2 obtained under the Cox model
and from our procedure. The kernel function K is taken equal to the Epanechnikov kernel :
b k (·|Xi ) (k = 0, 1)
K(u) = (3/4)(1−u2)I(|u| ≤ 1). The bandwidth h of the kernel estimators H
is taken proportional to n−2/7 so as to verify regularity condition (C1), i.e. h = cn−2/7 for
several values of c, namely c = 2, 3 and 4. In addition, we also used the cross-validation
(CV) procedure proposed by Li, Lin and Racine (2013) for kernel estimators of conditional
distribution functions. The CV procedure is implemented in the package np in R. For each
b 0 and H
b 1 and used the average
sample in our simulation, we calculated these bandwidths for H
of these two bandwidths.
The tables show that our estimator outperforms the one that is based on the Cox model,
even when the Cox model is correct. They also show that our estimator is only mildly
sensitive to the bandwidth, which could be explained by the fact that we average out the
effect of the bandwidth. We also see that the CV selection of the bandwidth is working rather
well, in the sense that the MSE is close to the smallest value among the MSE’s corresponding
to the three fixed bandwidths.
Next, we look at the estimation of the quartiles of the distribution FT,0 (·|x) when x =
βb
0.25. We estimate these quartiles by means of our nonparametric estimator FbT,0
(·|x) and
by means of the Cox model studied in Lu (2008). The results given in Tables 3 and 4 show
that, as could be expected, when the Cox model is not satisfied (i.e. when γ2 = 1 or 2), the
20

n γ2
150 0
1
2
300

0
1
2

Par.
β1
β2
β1
β2
β1
β2
β1
β2
β1
β2
β1
β2

c=2
Bias MSE
.092 .230
.002 .540
.063 .123
-.099 .340
.045 .100
-.109 .242
.021 .100
-.081 .266
-.004 .060
-.107 .181
-.015 .050
-.124 .157

c=3
Bias MSE
.012 .183
-.231 .438
-.021 .101
-.334 .331
-.029 .086
-.356 .277
-.029 .088
-.252 .268
-.048 .055
-.278 .208
-.059 .046
-.295 .198

c=4
Bias MSE
-.075 .162
-.506 .510
-.094 .094
-.605 .505
-.099 .083
-.632 .497
-.081 .081
-.461 .363
-.097 .055
-.482 .328
-.107 .049
-.498 .329

hCV
Bias MSE
.044 .224
-.158 .536
.017 .120
-.246 .397
-.005 .101
-.302 .382
.004 .099
-.148 .288
-.013 .061
-.150 .201
-.030 .052
-.197 .217

Cox
Bias MSE
.216 .523
.291 1.15
.147 .191
.278 .635
.124 .145
.263 .476
.099 .139
.135 .365
.097 .092
.215 .302
.077 .074
.181 .247

Table 1: Bias and MSE of βb1 and βb2 for two sample sizes, three values of γ2 , three bandwidths
of the form h = cn−2/7 and the bandwidth hCV obtained from cross-validation. Here,
P (cured) = 0.2 and P (censoring) = 0.4 for γ2 = 0. The Cox model is satisfied for γ2 = 0.

n γ2
150 0
1
2
300

0
1
2

Par.
β1
β2
β1
β2
β1
β2
β1
β2
β1
β2
β1
β2

c=2
Bias MSE
.082 .139
-.017 .418
.057 .076
-.107 .253
.049 .065
-.129 .202
.059 .074
-.068 .216
-.013 .036
-.117 .135
-.049 .028
-.295 .168

c=3
Bias MSE
.003 .116
-.244 .365
-.011 .063
-.328 .277
-.016 .056
-.361 .260
-.010 .060
-.237 .226
-.043 .034
-.282 .175
-.080 .030
-.496 .306

c=4
Bias MSE
-.060 .105
-.525 .477
-.068 .059
-.608 .478
-.071 .054
-.640 .493
-.058 .056
-.453 .326
-.074 .035
-.486 .306
-.037 .030
-.244 .194

hCV
Bias MSE
.031 .130
-.217 .419
.010 .071
-.294 .352
.006 .066
-.329 .358
.031 .073
-.153 .241
-.024 .035
-.187 .160
-.034 .030
-.217 .182

Cox
Bias MSE
.116 .189
.215 .618
.086 .099
.227 .434
.065 .078
.196 .329
.060 .083
.103 .257
.050 .047
.159 .197
.035 .037
.128 .156

Table 2: Bias and MSE of βb1 and βb2 for two sample sizes, three values of γ2 , three bandwidths
of the form h = cn−2/7 and the bandwidth hCV obtained from cross-validation. Here,
P (cured) = 0.3 and P (censoring) = 0.5 for γ2 = 0. The Cox model is satisfied for γ2 = 0.

MSE of the quartiles obtained under the Cox model is much higher than the corresponding
MSE obtained from our procedure. This shows the importance of having a model that does
21

n γ2
150 0

1

2

300

0

1

2

p
.25
.50
.75
.25
.50
.75
.25
.50
.75
.25
.50
.75
.25
.50
.75
.25
.50
.75

c=2
Bias MSE
.044 .061
.022 .049
.024 .083
.083 .053
.060 .051
.072 .089
.126 .060
.098 .058
.112 .107
.007 .032
-.001 .034
-.003 .053
.027 .026
.028 .031
.031 .055
.063 .031
.055 .033
.055 .058

c=3
Bias MSE
.025 .035
.003 .031
.006 .055
.102 .039
.092 .041
.114 .075
.218 .082
.189 .073
.210 .120
-.007 .023
-.013 .023
-.019 .034
.049 .021
.054 .025
.053 .040
.129 .037
.121 .040
.109 .056

c=4
Bias MSE
.030 .032
.003 .025
.001 .045
.155 .053
.144 .049
.154 .077
.325 .139
.291 .121
.296 .156
-.018 .018
-.026 .018
-.030 .026
.081 .021
.086 .025
.086 .038
.216 .065
.200 .062
.191 .074

hCV
Bias MSE
.033 .049
.011 .043
.011 .071
.096 .049
.073 .048
.091 .085
.202 .085
.172 .081
.187 .124
.002 .030
-.006 .031
-.009 .046
.033 .026
.034 .030
.033 .052
.096 .039
.086 .038
.084 .061

Cox
Bias MSE
.031 .027
.028 .023
.032 .038
.267 .104
.251 .093
.254 .117
.592 .401
.513 .308
.492 .322
.010 .013
.008 .013
.010 .023
.252 .081
.240 .076
.235 .086
.580 .366
.498 .274
.452 .248

Table 3: Bias and MSE of the conditional quantiles of order p = 0.25, 0.50 and 0.75 at
x = 0.25 for two sample sizes, three bandwidths of the form h = cn−2/7 and the bandwidth
hCV obtained from cross-validation. Here, P (cured) = 0.2 and P (censoring) = 0.4 for
γ2 = 0. The Cox model is satisfied for γ2 = 0.

not impose any assumptions on the distribution of the uncured individuals and which still
provides very accurate estimators for the logistic part of the model.
We also verify how close the distributions of βb1 and βb2 are to a normal distribution. We
know thanks to Theorem 4.2 that the estimators converge to a normal limit when n tends to
infinity. Figure 1 shows that for n = 150 the distribution is rather close to a normal limit,
especially for βb1 . The figure is based on 1000 samples generated from the above model with
P (cured) = 0.2 and P (censoring) = 0.4. The results for n = 300 (not shown here for space
constraints) are close to a straight line, showing that the results improve when n increases.
Finally, we verify the accuracy of the naive bootstrap proposed in Section 4.3. We
consider the above model, but restrict attention to n = 150 and to the case where P (cured) =
0.2 and P (censoring) = 0.4. Figure 2 shows boxplots of the variance of βb1 and βb2 obtained
from 250 bootstrap resamples for each of 500 samples. The bandwidth is h = 3n−2/7 . The
empirical variance of the 500 estimators of β1 and β2 is also added, and shows that the
bootstrap variance is well centered around the corresponding empirical variance.
22

n γ2
150 0

1

2

300

0

1

2

p
.25
.50
.75
.25
.50
.75
.25
.50
.75
.25
.50
.75
.25
.50
.75
.25
.50
.75

c=2
Bias MSE
.074 .095
.049 .069
.042 .098
.107 .072
.076 .060
.099 .113
.159 .085
.120 .066
.143 .135
.032 .042
.018 .039
.015 .060
.026 .033
.030 .036
.033 .060
.070 .040
.060 .038
.062 .066

c=3
Bias MSE
.020 .049
.005 .041
.007 .063
.116 .053
.103 .049
.118 .083
.232 .097
.211 .089
.228 .142
-.005 .029
-.011 .026
-.018 .037
.046 .025
.052 .027
.050 .042
.132 .042
.125 .044
.115 .061

c=4
Bias MSE
.021 .044
-.003 .031
-.005 .050
.151 .059
.142 .054
.154 .086
.322 .144
.296 .131
.308 .178
-.029 .023
-.035 .021
-.037 .028
.071 .023
.080 .026
.081 .040
.212 .064
.201 .064
.193 .078

hCV
Bias MSE
.046 .071
.023 .057
.024 .084
.119 .066
.088 .057
.109 .102
.233 .109
.195 .092
.219 .153
.016 .038
.007 .034
-.001 .054
.034 .031
.037 .034
.038 .055
.103 .046
.089 .043
.090 .070

Cox
Bias MSE
.039 .038
.031 .029
.034 .046
.274 .116
.255 .100
.257 .126
.586 .402
.508 .307
.489 .325
.018 .017
.013 .016
.016 .022
.252 .083
.234 .074
.232 .087
.570 .357
.492 .269
.451 .251

Table 4: Bias and MSE of the conditional quantiles of order p = 0.25, 0.50 and 0.75 at
x = 0.25 for two sample sizes, three bandwidths of the form h = cn−2/7 and the bandwidth
hCV obtained from cross-validation. Here, P (cured) = 0.3 and P (censoring) = 0.5 for
γ2 = 0. The Cox model is satisfied for γ2 = 0.

6

Data analysis

Let us now apply our estimation procedure on two medical data sets. The first one is
about 286 breast cancer patients with lymph-node-negative breast cancer treated between
1980 and 1995 (Wang et al. (2005)). The event of interest is distant-metastasis, and the
associated survival time is the distant metastasis-free survival time (defined as the time to
first distant progression or death, whichever comes first). 107 of the 286 patients experience
a relapse from breast cancer. The plot of the Kaplan-Meier estimator of the data is given in
Figure 3(a) and shows a large plateau at about 0.60. Furthermore, a large proportion of the
censored observations is in the plateau, which suggests that a cure model is appropriate for
these data. As a covariate we use the age of the patients, which ranges from 26 to 83 years
and the average age is about 54 years.
We estimate β using our estimator and using the estimator based on the Cox model. The
bandwidth h is selected using cross-validation, as in the simulation section. The estimated
23

0

1

2

3

1

2

2.5

3

−3 −2 −1

Theoretical Quantiles

0

1

2

3

Theoretical Quantiles

−3 −2 −1

0

1

2

Theoretical Quantiles

3

−3 −2 −1

0

1

2

Theoretical Quantiles

3

4
3
2
1

Sample Quantiles

5
4
3
2
1

Sample Quantiles

8
6
4
2
0

Sample Quantiles

5

Theoretical Quantiles

0

1.5

Sample Quantiles
−3 −2 −1

0.5

2.5
1.5
0.5

Sample Quantiles

4
3
2
1

Sample Quantiles

0
−3 −2 −1

−3 −2 −1

0

1

2

3

Theoretical Quantiles

Figure 1: QQ-plots of βb1 (first row) and βb2 (second row) for 1000 samples of size n = 150.
The first column corresponds to γ2 = 0, the second to γ2 = 1 and the third to γ2 = 2. The
bandwidth is h = 3n−2/7 .
intercept is -0.224 (with standard deviation equal to 0.447 obtained using a naive bootstrap
procedure), and the estimated slope parameter is -0.005 (with standard deviation equal
to 0.008). Under the Cox model the estimated intercept and slope are respectively 0.063
and -0.010. A 95% confidence interval is given by (−1.100, 0.653) for the intercept and
(−0.021, 0.011) for the slope, where the variance is again based on the naive bootstrap
procedure. The graph of the two estimators of the function φ(x) is given in Figure 3(b).
The estimated coefficients and curves are quite close to each other, suggesting that the Cox
model might be valid. This is also confirmed by Figure 3(c)-(d), which shows the estimation
of the survival function 1 − FT,0 (·|x) of the uncured patients for x = 48 and x = 60 based
on our estimation procedure and the procedure based on the Cox model. The figure shows
that the two estimators are close for both values of x.
Next, we analyse data provided by the Medical Birth Registry of Norway (see
http://folk.uio.no/borgan/abg-2008/data/data.html). The data set contains information on
births in Norway since 1967, related to a total of 53,558 women. We are interested in the
time between the birth of the first and the second child, for those mothers whose first child
died within the first year (n = 262). The covariate of interest is age (X), which is the age
of the mother at the birth of the first child. The age ranges from 16.8 to 29.8 years, with
an average of 23.2 years. The cure rate is the fraction of women who gave birth only once.
24

0.12

0.18

0.08

0.14
0.10
0.40

0.04

0.06

0.4
0.3
0.2
0.1

0.04

0.10

0.08

0.20

0.30

0.12

0.7
0.5
0.3
0.1

Figure 2: Boxplots of the variance of βb1 (first row) and βb2 (second row) obtained from 250
bootstrap resamples for each of 500 samples of size n = 150. The first column corresponds
to γ2 = 0, the second to γ2 = 1 and the third to γ2 = 2. The bandwidth is h = 3n−2/7 . The
empirical variance of the 500 estimators of β1 and β2 is also added (dashed line).

Figure 4(a) shows the Kaplan-Meier estimator, and suggests that a cure fraction is present.
As we did for the first data set, we analyse these data using the approach proposed in
this paper, and also using the Cox mixture cure model. The estimated intercept equals 1.952
using our model and 0.034 using the Cox model. The bootstrap confidence interval for the
intercept is (−1.577, 5.481) (the estimated standard deviation equals 1.801). The estimated
slope equals -0.041 respectively 0.052 using the two models. For our estimation procedure
the confidence interval is given by (−0.193, 0.111) (with estimated standard deviation equal
to 0.078). Figure 4(b) shows that the two estimators of the function φ(x) are quite different,
and have opposite slopes. Moreover, the survival function 1−FT,0 (·|x) of the uncured patients
is given in Figure 4(c)-(d) for x = 21 and x = 25. We see that the estimator based on the
Cox model is quite different from ours, suggesting that the Cox model might not be valid for
these data, although a formal test would need to confirm this. This is however beyond the
scope of this paper. Also note that the estimator of the cure proportion 1−φ(x) is increasing
under our model and decreasing under the Cox model. It seems however natural to believe
that the probability of having no second child (so the cure proportion) is increasing with
age, which is again an indication that the Cox model is not valid for these data.

25

0.0

0.25

0.35

1−P(cured|X)

0.4

0.45

(b)

0.8

(a)

1000

2000

3000

4000

5000

30

40

50

60

t

X

(c)

(d)

70

80

4000

5000

0.0

0.0

0.4

0.4

0.8

0.8

0

0

1000

2000

3000

4000

5000

0

t

1000

2000

3000
t

Figure 3: Analysis of the breast cancer data : (a) Kaplan-Meier estimator; (b) Graph of the
proposed estimator of φ(x) (solid curve) and of the estimator based on the Cox model (dashed
curve); (c) Estimation of 1 − FT,0 (·|x) using the proposed estimator (solid curve) and using
the estimator based on the Cox model (dashed curve) when x = 48; (d) Idem when x = 60.

7

Appendix : Proofs

Proof of Proposition 3.1. By the properties of the likelihood of a Bernoulli random variable
given that Y ∈ dt and X = x, we have
"
β
φ(x, β)FT,0
(dt | x)FC ((t, ∞) | x)
E δ log
H1 (dt | x)
h
i

β
FC (dt | x) φ(x, β)FT,0 ((t, ∞) | x) + 1 − φ(x, β)
+(1 − δ) log
Y ∈ dt, X = x ≤ 0.
H0 (dt | x)
Integrate with respect to Y and X and deduce that

E [log p(Y, δ, X; β)] ≤ E [log p(Y, δ, X; β0 )] .

26

If there exists some β 6= β0 such that the last inequality becomes an equality, then necessarily
p1 (Y, X; β) = 1 almost surely. Then
β
β0
FC ((t, ∞) | x)φ(x, β)FT,0
(dt | x) = FC ((t, ∞) | x)φ(x, β0 )FT,0
(dt | x), ∀ − ∞ < t ≤ τH (x),

for almost all x ∈ X . By Theorem 2.1, we deduce that necessarily β = β0 .
Lemma 7.1. Let conditions (4.2), (AC1) and (AC4) hold true. Then,
sup

sup |FbC ([t, ∞) | x) − FC ([t, ∞) | x)| = oP (1),

x∈X t∈(−∞,τ ]

sup sup

b0, H
b 1 )(t, x) − T1 (β, H0 , H1 )(t, x)| = oP (1),
sup |T1 (β, H

β∈B x∈X t∈(−∞,τ ]

where

T1 (β, H0, H1 )(t, x) = H([t, ∞) | x) − (1 − φ(x, β))FC ([t, ∞) | x),

(b)

0.0

0.75
0.65

0.4

1−P(cured|X)

0.8

0.85

(a)

1000

2000

3000

4000

18

20

22

24

t

X

(c)

(d)

26

28

30

0.0

0.0

0.4

0.4

0.8

0.8

0

0

1000

2000

3000

4000

0

t

1000

2000

3000

4000

t

Figure 4: Analysis of the second birth data : (a) Kaplan-Meier estimator; (b) Graph of the
proposed estimator of φ(x) (solid curve) and of the estimator based on the Cox model (dashed
curve); (c) Estimation of 1 − FT,0 (·|x) using the proposed estimator (solid curve) and using
the estimator based on the Cox model (dashed curve) when x = 21; (d) Idem when x = 25.

27

b 0, H
b 1 ) is defined similarly, but with H0 , H1 and FC replaced by H
b0, H
b 1 and FbC ,
and T1 (β, H
respectively. Moreover,
β
β
sup |FbT,0
([t, ∞) | x) − FT,0
([t, ∞) | x)| = oP (1).

sup sup

β∈B x∈X t∈(−∞,τ ]

Proof of Lemma 7.1. Let us first investigate the uniform convergence of the estimated cub C (· | x). For any t ∈ (−∞, τ ] let us write
mulative hazard measure Λ
Z
Z
b 0 (ds | x)
H0 (ds | x)
H
b
−
ΛC ((−∞, t] | x) − ΛC ((−∞, t] | x) =
b
(−∞,t] H([s, ∞) | x)
(−∞,t] H([s, ∞) | x)
#
"
Z
Z
b 0 (ds | x) − H0 (ds | x)
1
H
1
b 0 (ds | x) +
−
H
.
=
b
H([s, ∞) | x)
∞) | x) H([s, ∞) | x)
(−∞,t]
(−∞,t] H([s,

b 0 (ds | x) are well defined, since, with probability tending to
The integrals with respect to H
b 0 ((−∞, s] | x), s ∈ (−∞, τ ], is a function of bounded
1, for each x ∈ X , the map s 7→ H
variation. The uniform convergence Assumption (AC1) implies that
sup sup
x∈X t∈(−∞,τ ]

b C ((−∞, t] | x) − ΛC ((−∞, t] | x)
Λ

≤ c sup sup
x∈X t∈(−∞,τ ]

n

b 0 ([t, ∞) | x)−H0 ([t, ∞) | x) + H([t,
b
H
∞) | x)−H([t, ∞) | x)
H([τ, ∞) | x)2

o

,

for some constant c > 0. Next, by Duhamel’s identity (see Gill and Johansen 1990),
FbC ((t, ∞) | x) − FC ((t, ∞) | x) = −FC ((t, ∞) | x)
Z

FbC ([s, ∞) | x)  b
ΛC (ds | x) − ΛC (ds | x) .
×
(−∞,t] FC ((s, ∞) | x)

b C (· | x)
Then, the uniform convergence of FbC (· | x) follows from the uniform convergence of Λ
b0, H
b 1), and hence we omit
and condition (4.2). The same type of arguments apply for T1 (β, H
the details.
Next, since by conditions (4.2) and (AC4) we have
inf inf

inf

β∈B x∈X t∈(−∞,τ ]

[H([t, ∞) | x) − (1 − φ(x, β))FC ([t, ∞) | x)] > 0,

there exists some constant c > 0 with the property that


h
i
b
P inf inf
inf
H([t,
∞) | x) − (1 − φ(x, β))FbC ([t, ∞) | x) ≥ c > 0 → 1.
β∈B x∈X t∈(−∞,τ ]

β
Hence, the uniform convergence of FbT,0
(· | x) follows.

28

Proof of Theorem 4.1. Let us write
β
FT,0
({Yi }

| Xi ) =

ΛβT,0 ({Yi }

|

β
Xi )FT,0
([Yi , ∞)

β
FT,0
([Yi , ∞) | Xi )
| Xi ) = H1 ({Yi } | Xi )
,
T1 (β, H0 , H1 )(Yi , Xi )

where
T1 (β, H0, H1 )(t, x) = H([t, ∞) | x) − (1 − φ(x, β))FC ([t, ∞) | x).
Moreover, let
qi (β, H0 , H1 ) = q(β, H0, H1 )(Yi , δi , Xi ),
where, for t ∈ R, d ∈ {0, 1}, x ∈ X ,
β
q(β, H0, H1 )(t, d, x) = d{log φ(x, β) + log FT,0
([t, ∞) | x) − log T1 (β, H0 , H1 )(t, x)}
β
+ (1 − d) log{φ(x, β)FT,0
([t, ∞) | x) + 1 − φ(x, β)}.

Let

n

1X
Qn (β, H0 , H1) =
qi (β, H0 , H1 ).
n i=1

b0, H
b 1 ) that is defined as Qn (β, H0, H1 ), but with H0 , H1 ,
Similarly, let us consider Qn (β, H
β
b0, H
b 1 , FbC and Fbβ , respectively. Then the estimator βb in equation
FC and FT,0
replaced by H
T,0
(3.1) becomes
b 0, H
b 1 ).
βb = arg max Qn (β, H
β∈B

The first step is to check that

b0, H
b 1 ) − Qn (β, H0, H1 ) = oP (1).
sup Qn (β, H

(7.12)

β∈B

This follows directly from Lemma 7.1. Next, given our assumptions, it is easy to check that
for any t ∈ (−∞, τ ], d ∈ {0, 1}, x ∈ X ,
|q(β, H0 , H1 )(t, d, x) − q(β ′ , H0 , H1 )(t, d, x)| ≤ Ckβ − β ′ ka ,

∀β, β ′ ∈ B,

with a > 0 from Assumption (AC3) and some constant C depending only on c1 from Assumption (AC3) and the positive values inf x∈X H1 ({τ } | x) and inf x∈X H0 ((τ, ∞) | x). It
follows that the class {(t, d, x) → q(β, H0, H1 )(t, d, x) : β ∈ B} is Glivenko-Cantelli. Hence,
sup |Qn (β, H0 , H1 ) − Q(β, H0, H1 )| = oP (1),
β∈B

where Q = E(Qn ). Finally, Proposition 3.1 guarantees that
β0 = arg max Q(β, H0 , H1 ).
β∈B

Gathering the facts, we deduce that βb − β0 = oP (1).
29

Proof of Theorem 4.2. We show the asymptotic normality of our estimator by verifying the
high-level conditions in Theorem 2 in Chen et al. (2003). First of all, for the consistency we
refer to Section 4.1, whereas conditions (2.1) and (2.2) in Chen et al. (2003) are satisfied by
construction and thanks to assumption (AN1), respectively. Concerning (2.3), first note that
the expression inside the expected value in ∇η M(β, η0 )[η − η0 ] is linear in ∇η Tj (β, η0 )[η − η0 ]
(j = 1, 2, 3, 4). Hence, we will focus attention on the latter Gâteaux derivatives. First,
∇η T1 (β, η0)[η − η0 ](t, x) = (η1 − η01 + η2 − η02 )(t, x) − (1 − φ(x, β))∇η T2 (η0 )[η − η0 ](t, x).
Using Duhamel’s formula (see Gill and Johansen 1990), we can write
Z
1
∇η T2 (η0 )[η − η0 ](t, x) = −T2 (η0 )(t, x)
−∞<u<t (η01 + η02 )(u, x) − η01 ({u}, x)
n
η01 (du, x)(η1 + η2 − η01 − η02 )(u, x) o
.
× (η1 − η01 )(du, x) −
(η01 + η02 )(u, x)
In a similar way, we find that
Z

1
−∞<u<t T1 (β, η0 )(u, x) − η02 ({u}, x)


n
η02 (du, x) T1 (β, η) − T1 (β, η0 ) (u, x) o
.
× (η2 − η02 )(du, x) −
T1 (β, η0 )(u, x)

∇η T3 (β, η0 )[η − η0 ](t, x) = −T3 (β, η0 )(t, x)

Finally,
∇η T4 (β, η0 )[η − η0 ](t, x)
(
Z
∇η T2 (η0 )[η − η0 ](s, x)η02 (ds, x) + T2 (η0 )(s, x)(η2 − η02 )(ds, x)


= −∇β φ(x, β)
T1 (β, η0 )(s, x) T1 (β, η0 )(s, x) − η02 ({s}, x)
(−∞,t)
T2 (η0 )(s, x)η02 (ds, x)∇η T1 (β, η0 )[η − η0 ](s, x)
−
2 

T1 (β, η0 )(s, x) T1 (β, η0)(s, x) − η02 ({s}, x)

)
T2 (η0 )(s, x)η02 (ds, x) ∇η T1 (β, η0 )[η − η0 ](s, x) − (η2 − η02 )({s}, x)
.
−

2
T1 (β, η0 )(s, x) T1 (β, η0 )(s, x) − η02 ({s}, x)

Note that all denominators in ∇η Tj (β, η0 )[η − η0 ](t, x) are bounded away from zero, thanks
to (4.10) and (4.11). By tedious but rather elementary arguments, it follows from these
formulae that
k∇η Tj (β, η0 )[η − η0 ] − ∇η Tj (β0 , η0 )[η − η0 ]kH ≤ Ckβ − β0 kkη − η0 kH ,
for some constant C. Hence, it can be easily seen that ∇η Tj (β, η0)[η − η0 ] satisfies the
second property in assumption (2.3) in Chen et al. (2003), and hence the same holds true
for ∇η M(β, η0 )[η − η0 ]. Similarly, by decomposing Tj (β, η) − Tj (β, η0 ) − ∇η Tj (β, η0 )[η − η0 ]
30

using Taylor-type arguments (in η), the first property in assumption (2.3) is easily seen to
hold true.
Next, conditions (2.4) and (2.6) are satisfied thanks to Assumption (AN4) and because
it follows from the above calculations of ∇η Tj (β, η0 )[η − η0 ] (j = 1, 2, 3, 4) that


Z
X
∇η M(β0 , η0 )[η − η0 ] =
E ψ1k (Y, X)
ψ2k (u, X)d ((ηk − η0k )(u, X))
+

X



E ψ3k (Y, X)

k,ℓ∈{0,1}

−∞<u<Y

k∈{0,1}

Z

−∞<u<Y




ψ4k (u, X)ψ5k (ηk − η0k )(u, X) dHℓ (u | X) (7.13)

for certain measurable functions ψjk (j = 1, . . . , 5; k = 0, 1).
It remains to verify condition (2.5). Note that
|m(t, δ, x; β2 , η2 ) − m(t, δ, x; β1 , η1 )| ≤ C1 (t, δ, x)kβ2 − β1 k + C2 (t, δ, x)kη2 − η1 kH
for some functions Cj satisfying E[Cj2 (Y, δ, X)] < ∞ (j = 1, 2), and hence (2.5) follows from
assumption (AN5) and Theorem 3 in Chen et al. (2003). This finishes the proof.
Proof of Theorem 4.3. To prove this theorem we will check the conditions of Theorem B in
Chen et al. (2003), which gives high level conditions under which the naive bootstrap
is consistent. The only difference between their setting and our setting is that we are
proving bootstrap consistency in P-probability, whereas their result holds true a.s. [P]. As
a consequence, in their high level conditions we can replace all a.s. [P] statements by the
corresponding statements in P-probability.
First of all, it follows from assumption (AN1) that condition (2.2) in Chen et al. (2003)
holds with η0 replaced by any η in a neighborhood of η0 , and from the proof of Theorem 4.2
it follows that the same holds true for condition (2.3). Next, conditions (2.4B) and (2.6B)
in Chen et al. follow from the fact that we assume that assumption (AN4) continues to hold
b k − Hk by H
b∗ − H
b k (k = 0, 1). It remains to verify condition (2.5’B) in
true if we replace H
k
Chen et al. This follows from Theorem 3 in Chen et al., whose conditions have been verified
already for our Theorem 4.2.

References
[1] Akritas, M.G. & Van Keilegom, I. (2001). Nonparametric estimation of the residual
distribution. Scand. J. Statist. 28, 549–568.
[2] Amico, M., Legrand, C. & Van Keilegom, I. (2017). The single-index/Cox mixture
cure model (submitted).
[3] Boag, J.W. (1949). Maximum likelihood estimates of the proportion of patients cured
by cancer therapy. J. Roy. Statist. Soc. - Series B 11, 15–53.
31

[4] Chen, X., Linton, O. & Van Keilegom, I. (2003). Estimation of semiparametric
models when the criterion function is not smooth. Econometrica 71, 1591–1608.
[5] Fang, H.B., Li, G. & Sun, J. (2005). Maximum likelihood estimation in a semiparametric logistic/proportional-hazards mixture model. Scand. J. Statist. 32, 59–75.
[6] Farewell, V.T. (1982). The use of mixture models for the analysis of survival data
with long-term survivors. Biometrics 38, 1041–1046.
[7] Gill, R.D. (1994). Lectures on survival analysis. Lectures on probability theory: Ecole
d’été de probabilités de Saint-Flour XXII. Lecture notes in mathematics 1581. Springer.
[8] Gill, R.D. & Johansen, S. (1990). A survey of product-integration with a view toward
application in survival analysis. Ann. Statist. 18, 1501–1555.
[9] Kuk, A.Y.C. & Chen, C.-H. (1992). A mixture model combining logistic regression
with proportional hazards regression. Biometrika 79, 531–541.
[10] Li, Q., Lin, J. & Racine, J.S. (2013). Optimal bandwidth selection for nonparametric conditional distribution and quantile functions. J. Buss. Econ. Statist. 31, 57–65.
[11] López-Cheda, A., Cao, R., Jácome M.A. & Van Keilegom, I. (2017). Nonparametric incidence estimation and bootstrap bandwidth selection in mixture cure models.
Comput. Statist. Data Anal. 105, 144–165.
[12] Lopez, O. (2011). Nonparametric estimation of the multivariate distribution function
in a censored regression model with applications. Commun. Stat. - Theory Meth. 40,
2639–2660.
[13] Lu, W. (2008). Maximum likelihood estimation in the proportional hazards cure model.
Ann. Inst. Stat. Math. 60, 545–574.
[14] Maller, R.A. & Zhou, S. (1996). Survival Analysis with Long Term Survivors.
Wiley, New York.
[15] Meeker, W.Q. (1987). Limited failure population life tests: Application to integrated
circuit reliability. Technometrics 29, 51–65.
[16] Othus, M., Li, Y. & Tiwari, R.C. (2009). A class of semiparametric mixture cure
survival models with dependent censoring. J. Amer. Statist. Assoc. 104, 1241–1250.
[17] Peng, Y. & Taylor, J.M.G. (2014). Cure models. In: Klein, J., van Houwelingen, H.,
Ibrahim, J. G., and Scheike, T. H., editors, Handbook of Survival Analysis, Handbooks
of Modern Statistical Methods series, chapter 6, pages 113-134. Chapman & Hall, Boca
Raton, FL, USA.
32

[18] Schmidt, P. & Witte, A.D. (1989). Predicting criminal recidivism using split population survival time models. J. Econometrics 40, 141–159.
[19] Sy, J.P. & Taylor, J.M.G. (2000). Estimation in a Cox proportional hazards cure
model. Biometrics 56, 227–236.
[20] Taylor, J.M.G. (1995). Semi-parametric estimation in failure time mixture models.
Biometrics 51, 899–907.
[21] van der Vaart, A.D. (1998). Asymptotic Statistics. Cambridge University Press.
[22] Wang, Y., Klijn, J.G.M., Sieuwerts, A.M., Look, M.P., Yang, F., Talantov, D., Timmermans, M., Meijet-van Gelder, M.E.M., Yu, J., Jatkoe, T.,
Berns, E.M.J.J., Atkins, D. & Foekens, J.A. (2005). Gene-expression profiles to
predict distant metastasis of lymph-node-negative primary breast cancer. The Lancet
365, 671–679.
[23] Xu, J. & Peng, Y. (2014). Nonparametric cure rate estimation with covariates. Canad.
J. Statist. 42, 1–17.

33

