A General Approach for Cure Models in Survival Analysis

TL;DR

This paper introduces a flexible, general approach for estimating mixture cure models in survival analysis with censored data, combining parametric cure probability with nonparametric survival functions.

Contribution

It proposes a novel inversion-based estimation method that allows flexible modeling of the survival function and proves its statistical properties.

Findings

The model is identifiable and statistically consistent.

Kernel estimators improve nonparametric survival function estimation.

The approach outperforms Cox mixture models in simulations.

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.

Tables

Table 1. Table 1: Bias and MSE of β ^ 1 subscript ^ 𝛽 1 \widehat{\beta}_{1} and β ^ 2 subscript ^ 𝛽 2 \widehat{\beta}_{2} for two sample sizes, three values of γ 2 subscript 𝛾 2 \gamma_{2} , three bandwidths of the form h = c ​ n − 2 / 7 ℎ 𝑐 superscript 𝑛 2 7 h=cn^{-2/7} and the bandwidth h C ​ V subscript ℎ 𝐶 𝑉 h_{CV} obtained from cross-validation. Here, P ​ ( c ​ u ​ r ​ e ​ d ) = 0.2 𝑃 𝑐 𝑢 𝑟 𝑒 𝑑 0.2 P(cured)=0.2 and P ​ ( c ​ e ​ n ​ s ​ o ​ r ​ i ​ n ​ g ) = 0.4 𝑃 𝑐 𝑒 𝑛 𝑠 𝑜 𝑟 𝑖 𝑛 𝑔 0.4 P(censoring)=0.4 for γ 2 = 0 subscript 𝛾 2 0 \gamma_{2}=0 . The Cox model is satisfied for γ 2 = 0 subscript 𝛾 2 0 \gamma_{2}=0 .

			$c=2$		$c=3$		$c=4$		$h_{CV}$		Cox
$n$	${\gamma }_2$	Par.	Bias	MSE	Bias	MSE	Bias	MSE	Bias	MSE	Bias	MSE
150	0	${\beta }_1$	.092	.230	.012	.183	-.075	.162	.044	.224	.216	.523
		${\beta }_2$	.002	.540	-.231	.438	-.506	.510	-.158	.536	.291	1.15
	1	${\beta }_1$	.063	.123	-.021	.101	-.094	.094	.017	.120	.147	.191
		${\beta }_2$	-.099	.340	-.334	.331	-.605	.505	-.246	.397	.278	.635
	2	${\beta }_1$	.045	.100	-.029	.086	-.099	.083	-.005	.101	.124	.145
		${\beta }_2$	-.109	.242	-.356	.277	-.632	.497	-.302	.382	.263	.476
300	0	${\beta }_1$	.021	.100	-.029	.088	-.081	.081	.004	.099	.099	.139
		${\beta }_2$	-.081	.266	-.252	.268	-.461	.363	-.148	.288	.135	.365
	1	${\beta }_1$	-.004	.060	-.048	.055	-.097	.055	-.013	.061	.097	.092
		${\beta }_2$	-.107	.181	-.278	.208	-.482	.328	-.150	.201	.215	.302
	2	${\beta }_1$	-.015	.050	-.059	.046	-.107	.049	-.030	.052	.077	.074
		${\beta }_2$	-.124	.157	-.295	.198	-.498	.329	-.197	.217	.181	.247

Table 2. Table 2: Bias and MSE of β ^ 1 subscript ^ 𝛽 1 \widehat{\beta}_{1} and β ^ 2 subscript ^ 𝛽 2 \widehat{\beta}_{2} for two sample sizes, three values of γ 2 subscript 𝛾 2 \gamma_{2} , three bandwidths of the form h = c ​ n − 2 / 7 ℎ 𝑐 superscript 𝑛 2 7 h=cn^{-2/7} and the bandwidth h C ​ V subscript ℎ 𝐶 𝑉 h_{CV} obtained from cross-validation. Here, P ​ ( c ​ u ​ r ​ e ​ d ) = 0.3 𝑃 𝑐 𝑢 𝑟 𝑒 𝑑 0.3 P(cured)=0.3 and P ​ ( c ​ e ​ n ​ s ​ o ​ r ​ i ​ n ​ g ) = 0.5 𝑃 𝑐 𝑒 𝑛 𝑠 𝑜 𝑟 𝑖 𝑛 𝑔 0.5 P(censoring)=0.5 for γ 2 = 0 subscript 𝛾 2 0 \gamma_{2}=0 . The Cox model is satisfied for γ 2 = 0 subscript 𝛾 2 0 \gamma_{2}=0 .

			$c=2$		$c=3$		$c=4$		$h_{CV}$		Cox
$n$	${\gamma }_2$	Par.	Bias	MSE	Bias	MSE	Bias	MSE	Bias	MSE	Bias	MSE
150	0	${\beta }_1$	.082	.139	.003	.116	-.060	.105	.031	.130	.116	.189
		${\beta }_2$	-.017	.418	-.244	.365	-.525	.477	-.217	.419	.215	.618
	1	${\beta }_1$	.057	.076	-.011	.063	-.068	.059	.010	.071	.086	.099
		${\beta }_2$	-.107	.253	-.328	.277	-.608	.478	-.294	.352	.227	.434
	2	${\beta }_1$	.049	.065	-.016	.056	-.071	.054	.006	.066	.065	.078
		${\beta }_2$	-.129	.202	-.361	.260	-.640	.493	-.329	.358	.196	.329
300	0	${\beta }_1$	.059	.074	-.010	.060	-.058	.056	.031	.073	.060	.083
		${\beta }_2$	-.068	.216	-.237	.226	-.453	.326	-.153	.241	.103	.257
	1	${\beta }_1$	-.013	.036	-.043	.034	-.074	.035	-.024	.035	.050	.047
		${\beta }_2$	-.117	.135	-.282	.175	-.486	.306	-.187	.160	.159	.197
	2	${\beta }_1$	-.049	.028	-.080	.030	-.037	.030	-.034	.030	.035	.037
		${\beta }_2$	-.295	.168	-.496	.306	-.244	.194	-.217	.182	.128	.156

Table 3. Table 3: Bias and MSE of the conditional quantiles of order p = 0.25 , 0.50 𝑝 0.25 0.50 p=0.25,0.50 and 0.75 0.75 0.75 at x = 0.25 𝑥 0.25 x=0.25 for two sample sizes, three bandwidths of the form h = c ​ n − 2 / 7 ℎ 𝑐 superscript 𝑛 2 7 h=cn^{-2/7} and the bandwidth h C ​ V subscript ℎ 𝐶 𝑉 h_{CV} obtained from cross-validation. Here, P ​ ( c ​ u ​ r ​ e ​ d ) = 0.2 𝑃 𝑐 𝑢 𝑟 𝑒 𝑑 0.2 P(cured)=0.2 and P ​ ( c ​ e ​ n ​ s ​ o ​ r ​ i ​ n ​ g ) = 0.4 𝑃 𝑐 𝑒 𝑛 𝑠 𝑜 𝑟 𝑖 𝑛 𝑔 0.4 P(censoring)=0.4 for γ 2 = 0 subscript 𝛾 2 0 \gamma_{2}=0 . The Cox model is satisfied for γ 2 = 0 subscript 𝛾 2 0 \gamma_{2}=0 .

			$c=2$		$c=3$		$c=4$		$h_{CV}$		Cox
$n$	${\gamma }_2$	$p$	Bias	MSE	Bias	MSE	Bias	MSE	Bias	MSE	Bias	MSE
150	0	$.25$	.044	.061	.025	.035	.030	.032	.033	.049	.031	.027
		$.50$	.022	.049	.003	.031	.003	.025	.011	.043	.028	.023
		$.75$	.024	.083	.006	.055	.001	.045	.011	.071	.032	.038
	1	$.25$	.083	.053	.102	.039	.155	.053	.096	.049	.267	.104
		$.50$	.060	.051	.092	.041	.144	.049	.073	.048	.251	.093
		$.75$	.072	.089	.114	.075	.154	.077	.091	.085	.254	.117
	2	$.25$	.126	.060	.218	.082	.325	.139	.202	.085	.592	.401
		$.50$	.098	.058	.189	.073	.291	.121	.172	.081	.513	.308
		$.75$	.112	.107	.210	.120	.296	.156	.187	.124	.492	.322
300	0	$.25$	.007	.032	-.007	.023	-.018	.018	.002	.030	.010	.013
		$.50$	-.001	.034	-.013	.023	-.026	.018	-.006	.031	.008	.013
		$.75$	-.003	.053	-.019	.034	-.030	.026	-.009	.046	.010	.023
	1	$.25$	.027	.026	.049	.021	.081	.021	.033	.026	.252	.081
		$.50$	.028	.031	.054	.025	.086	.025	.034	.030	.240	.076
		$.75$	.031	.055	.053	.040	.086	.038	.033	.052	.235	.086
	2	$.25$	.063	.031	.129	.037	.216	.065	.096	.039	.580	.366
		$.50$	.055	.033	.121	.040	.200	.062	.086	.038	.498	.274
		$.75$	.055	.058	.109	.056	.191	.074	.084	.061	.452	.248

Table 4. Table 4: Bias and MSE of the conditional quantiles of order p = 0.25 , 0.50 𝑝 0.25 0.50 p=0.25,0.50 and 0.75 0.75 0.75 at x = 0.25 𝑥 0.25 x=0.25 for two sample sizes, three bandwidths of the form h = c ​ n − 2 / 7 ℎ 𝑐 superscript 𝑛 2 7 h=cn^{-2/7} and the bandwidth h C ​ V subscript ℎ 𝐶 𝑉 h_{CV} obtained from cross-validation. Here, P ​ ( c ​ u ​ r ​ e ​ d ) = 0.3 𝑃 𝑐 𝑢 𝑟 𝑒 𝑑 0.3 P(cured)=0.3 and P ​ ( c ​ e ​ n ​ s ​ o ​ r ​ i ​ n ​ g ) = 0.5 𝑃 𝑐 𝑒 𝑛 𝑠 𝑜 𝑟 𝑖 𝑛 𝑔 0.5 P(censoring)=0.5 for γ 2 = 0 subscript 𝛾 2 0 \gamma_{2}=0 . The Cox model is satisfied for γ 2 = 0 subscript 𝛾 2 0 \gamma_{2}=0 .

			$c=2$		$c=3$		$c=4$		$h_{CV}$		Cox
$n$	${\gamma }_2$	$p$	Bias	MSE	Bias	MSE	Bias	MSE	Bias	MSE	Bias	MSE
150	0	$.25$	.074	.095	.020	.049	.021	.044	.046	.071	.039	.038
		$.50$	.049	.069	.005	.041	-.003	.031	.023	.057	.031	.029
		$.75$	.042	.098	.007	.063	-.005	.050	.024	.084	.034	.046
	1	$.25$	.107	.072	.116	.053	.151	.059	.119	.066	.274	.116
		$.50$	.076	.060	.103	.049	.142	.054	.088	.057	.255	.100
		$.75$	.099	.113	.118	.083	.154	.086	.109	.102	.257	.126
	2	$.25$	.159	.085	.232	.097	.322	.144	.233	.109	.586	.402
		$.50$	.120	.066	.211	.089	.296	.131	.195	.092	.508	.307
		$.75$	.143	.135	.228	.142	.308	.178	.219	.153	.489	.325
300	0	$.25$	.032	.042	-.005	.029	-.029	.023	.016	.038	.018	.017
		$.50$	.018	.039	-.011	.026	-.035	.021	.007	.034	.013	.016
		$.75$	.015	.060	-.018	.037	-.037	.028	-.001	.054	.016	.022
	1	$.25$	.026	.033	.046	.025	.071	.023	.034	.031	.252	.083
		$.50$	.030	.036	.052	.027	.080	.026	.037	.034	.234	.074
		$.75$	.033	.060	.050	.042	.081	.040	.038	.055	.232	.087
	2	$.25$	.070	.040	.132	.042	.212	.064	.103	.046	.570	.357
		$.50$	.060	.038	.125	.044	.201	.064	.089	.043	.492	.269
		$.75$	.062	.066	.115	.061	.193	.078	.090	.070	.451	.251

Equations

$$F_T((t,\infty ]\mid x)$$

$$\{\varphi (\cdot ,\beta ):\beta \in B\},$$

$$C\perp T\mid X\text{ and }\mathbb{P}(C<\infty )=1.$$

$$\mathbb{P}(T=C)=0,$$

$$H_1((-\infty ,t]\mid x)$$

$$H_0((-\infty ,t]\mid x)$$

$$H((-\infty ,\infty )\mid x)=1,\forall x\in \mathcal{X}.$$

$$F_C((-\infty ,t]\mid x)=\mathbb{P}(C\le t\mid X=x)\text{ and }{F}_T((-\infty ,t]\mid x)=\mathbb{P}(T\le t\mid X=x),x\in \mathcal{X}.$$

$$H_1(dt\mid x)=F_C([t,\infty )\mid x)F_T(dt\mid x),H_0(dt\mid x)=F_T([t,\infty ]\mid x)F_C(dt\mid x),$$

$${\Lambda }_T(dt\mid x)=\frac{F_T(dt\mid x)}{F_T([t,\infty ]\mid x)}\text{ and }{\Lambda }_C(dt\mid x)=\frac{F_C(dt\mid x)}{F_C([t,\infty )\mid x)},x\in \mathcal{X}.$$

$${\Lambda }_T(dt\mid x)=\frac{H_1(dt\mid x)}{H([t,\infty )\mid x)}\text{ and }{\Lambda }_C(dt\mid x)=\frac{H_0(dt\mid x)}{H([t,\infty )\mid x)}.$$

$$F_T((t,\infty ]\mid x)$$

$$F_C((t,\infty )\mid x)$$

$$\mathbb{P}(T>{\tau }_{H_1}(x)\mid x)=\prod _{t\in (-\infty ,{\tau }_{H_1}(x)]}\{1-{\Lambda }_T(dt\mid x)\},$$

$$\mathbb{P}(T>{\tau }_{H_1}(x)\mid x)=\mathbb{P}(T=\infty \mid x),$$

$$\mathbb{P}(C=\infty \mid x)=\prod _{t\in \mathbb{R}}\{1-{\Lambda }_C(dt\mid x)\}=0,\forall x\in \mathcal{X}.$$

$$H([t,\infty )\mid x)$$

$$F_T([t,\infty )\mid x)=\frac{H([t,\infty )\mid x)-\mathbb{P}(T=\infty \mid x)F_C([t,\infty )\mid x)}{F_C([t,\infty )\mid x)}.$$

$${\Lambda }_{T,0}(dt\mid x)\stackrel{def}{=}\frac{F_{T,0}(dt\mid x)}{F_{T,0}([t,\infty )\mid x)}=\frac{F_T(dt\mid x)}{F_T([t,\infty )\mid x)}$$

$${\Lambda }_{T,0}(dt\mid x)$$

$$F_{T,0}((t,\infty )\mid x)=\prod _{-\infty <s\le t}\{1-{\Lambda }_{T,0}(ds\mid x)\},t\in \mathbb{R},x\in \mathcal{X}.$$

$${\Lambda }_{T,0}^{\beta }(dt\mid x)=\frac{H_1(dt\mid x)}{H([t,\infty )\mid x)-[1-\varphi (x,\beta )]F_C([t,\infty )\mid x)},$$

$$F_{T,0}^{\beta }((t,\infty )\mid x)=\prod _{-\infty <s\le t}\left\{1-{\Lambda }_{T,0}^{\beta }(ds\mid x)\right\}.$$

$$F_{Y,\delta }^{\beta }(dt,1\mid x)=\varphi (x,\beta )F_C((t,\infty )\mid x)F_{T,0}^{\beta }(dt\mid x)$$

$$F_{Y,\delta }^{\beta }(dt,0\mid x)=[F_{T,0}^{\beta }((t,\infty )\mid x)\varphi (x,\beta )+1-\varphi (x,\beta )]F_C(dt\mid x).$$

$$F_{Y,\delta }(\cdot ,\cdot \mid x)=F_{Y,\delta }^{{\beta }_0}(\cdot ,\cdot \mid x),\forall x\in \mathcal{X}.$$

$$\varphi (X,\beta )=\varphi (X,\tilde{\beta }),\text{almost surely}\Rightarrow \beta =\tilde{\beta },$$

$$F_{Y,\delta }^{{\beta }_0}(\cdot ,\cdot \mid x)=F_{Y,\delta }^{\tilde{\beta }}(\cdot ,\cdot \mid x),\forall x\in \mathcal{X}\Rightarrow \varphi (x,{\beta }_0)=\varphi (x,\tilde{\beta }),\forall x\in \mathcal{X}.$$

$$\varphi (x,{\beta }_0)F_C((t,\infty )\mid x)F_{T,0}^{{\beta }_0}(dt\mid x)=\varphi (x,\tilde{\beta })F_C((t,\infty )\mid x)F_{T,0}^{\tilde{\beta }}(dt\mid x),$$

$$\varphi (x,{\beta }_0)F_{T,0}^{{\beta }_0}(dt\mid x)=\varphi (x,\tilde{\beta })F_{T,0}^{\tilde{\beta }}(dt\mid x),$$

$${\Lambda }_{T,0}(dt\mid x)=\frac{H_1(dt\mid x)}{H([t,\infty )\mid x)-[1-\varphi (x,\beta )]F_C([t,\infty )\mid x)}=\exp (x^{\top }\gamma ){\Lambda }_0(dt),$$

Full text

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.

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 limited-failure 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 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 The model

2.1 A general class of mixture cure models

Let $T$ denote (a possible monotone transformation of) the lifetime of interest that takes values in $(-\infty,\infty]$. A cured observation corresponds to the event $\{T=\infty\}$, and in the following this event is allowed to have a positive probability. Let $X$ be a covariate vector with support $\mathcal{X}$ belonging to a general covariate space. The covariate vector could include discrete and continuous components. The survival function $F_{T}((t,\infty]\mid x)=\mathbb{P}(T>t\mid X=x)$ for $t\in\mathbb{R}$ and $x\in\mathcal{X}$ can be written as

[TABLE]

where $\phi(x)=\mathbb{P}(T<\infty\mid X=x)$ and $F_{T,0}((t,\infty)\mid x)=\mathbb{P}(T>t\mid X=x,T<\infty)$. Depending on which model is used for $\phi(x)$ and $F_{T,0}(\cdot\mid x)$, one obtains a parametric, semiparametric or nonparametric model, called a ‘mixture cure model’. In the literature, one often assumes that $\phi(x)$ follows a logistic model, i.e. $\phi(x)=\exp(a+x^{\top}b)/[1+\exp(a+x^{\top}b)]$ for some $(a,b^{\top})^{\top}\in\mathbb{R}^{d+1}$. Recently, semiparametric models (like a single-index model as in Amico et al. 2017) or nonparametric models (as in Xu & Peng 2014 or López-Cheda et al. 2017) have been proposed. As for the survival function $F_{T,0}(\cdot\mid 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 $\phi(x)$ parametrically, i.e. we assume that $\phi(\cdot)$ belongs to the family of conditional probability functions

[TABLE]

where $\phi(\cdot,\beta)$ takes values in the interval $(0,1)$, $\beta$ 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 $F_{T,0}(\cdot\mid x)$ we do not impose any assumptions in order to have a flexible and rich class of models for $F_{T}(\cdot\mid x)$ to choose from. Later on we will see that for the estimation of $F_{T,0}(\cdot\mid 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.

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,\delta)$, where $Y=T\wedge C$, $\delta=\textbf{1}\{T\leq 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 $\delta$. Let us assume that

[TABLE]

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 $\mathbb{P}(C<\infty\mid 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

[TABLE]

which is commonly used in survival analysis, and which implies that $\mathbb{P}(T=C\mid 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 $\phi$, $F_{T}$ and $F_{C}$.

The observations are characterized by the conditional sub-probabilities

[TABLE]

Then $H((-\infty,t]\mid x)\stackrel{{\scriptstyle def}}{{=}}\mathbb{P}(Y\leq t\mid X=x)=H_{0}((-\infty,t]\mid x)+H_{1}((-\infty,t]\mid x).$ Since we assume that $Y$ is finite, we have

[TABLE]

For $j\in\{0,1\}$ and $x\in\mathcal{X},$ let $\tau_{H_{j}}(x)=\sup\{t:H_{j}([t,\infty)\mid x)>0\}$ denote the right endpoint of the support of the conditional sub-probability $H_{j}$. Let us define $\tau_{H}(x)$ in a similar way and note that $\tau_{H}(x)=\max\{\tau_{H_{0}}(x),\tau_{H_{1}}(x)\}$. Note that $\tau_{H_{0}}(x)$, $\tau_{H_{1}}(x)$ and $\tau_{H}(x)$ can equal infinity, even though $Y$ only takes finite values. For $-\infty<t\leq\infty,$ let us define the conditional probabilities

[TABLE]

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

[TABLE]

and $H([t,\infty)\mid x)=F_{T}([t,\infty]\mid x)F_{C}([t,\infty)\mid x).$ These equations could be solved and thus they allow to express the functions $F_{T}(\cdot\mid x)$ and $F_{C}(\cdot\mid x)$ in an unique way as explicit transformations of the functions $H_{0}(\cdot\mid x)$ and $H_{1}(\cdot\mid x).$ For this purpose, let us consider the conditional cumulative hazard measures

[TABLE]

The model equations yield

[TABLE]

Then, we can write the following functionals of $H_{0}(\cdot\mid x)$ and $H_{1}(\cdot\mid x)$ :

[TABLE]

where $\prod_{s\in A}$ stands for the product-integral over the set $A$ (see Gill and Johansen 1990).

Moreover, if $\tau_{H_{1}}(x)<\infty$, then

[TABLE]

but there is no way to identify the conditional law of $T$ beyond $\tau_{H_{1}}(x).$ Therefore, we will impose

[TABLE]

i.e. $\prod_{t\in\mathbb{R}}\{1-\Lambda_{T}(dt\mid x)\}=\prod_{-\infty<t\leq\tau_{H_{1}}(x)}\{1-\Lambda_{T}(dt\mid x)\}$. Note that if $\tau_{H_{1}}(x)=\infty,$ condition (2.6) is no longer an identification restriction, but just a simple consequence of the definition of $\Lambda_{T}(\cdot\mid x)$. Finally, the condition that $\mathbb{P}(C<\infty)=1$ in (2.1) can be re-expressed by saying that we assume that $H_{0}(\cdot\mid x)$ and $H_{1}(\cdot\mid x)$ are such that

[TABLE]

Let us point out that this condition is satisfied only if $\tau_{H_{1}}(x)\leq\tau_{H_{0}}(x)$. Indeed, if $\tau_{H_{1}}(x)>\tau_{H_{0}}(x)$ then necessarily $\tau_{H_{0}}(x)<\tau_{H}(x)$ and so $H([\tau_{H_{0}}(x),\infty)\mid x)>0.$ Hence, $\Lambda_{C}(\mathbb{R}\mid x)=\Lambda_{C}((-\infty,\tau_{H_{0}}(x)]\mid x)<\infty$, and thus $\mathbb{P}(C=\infty\mid x)>0$, which contradicts (2.7).

It is important to understand that any two conditional sub-probabilities $H_{0}(\cdot\mid x)$ and $H_{1}(\cdot\mid x)$ satisfying conditions (2.1) and (2.2) define uniquely $F_{T}(\cdot\mid x)$ and $F_{C}(\cdot\mid x)$. Indeed, $F_{T}(\cdot\mid x)$ is precisely the probability distribution of $T$ given $X=x$ with all the mass beyond $\tau_{H_{1}}(x)$ concentrated at infinity. In general, $F_{T}(\cdot\mid x)$ and $F_{C}(\cdot\mid x)$ are only functionals of $H_{0}(\cdot\mid x)$ and $H_{1}(\cdot\mid x)$.

We will assume conditions (2.1), (2.2) and (2.6) throughout the paper.

2.3 A key point for the new approach: the inversion formulae

Write

[TABLE]

and thus

[TABLE]

Consider the conditional cumulative hazard measure for the finite values of the lifetime of interest:

[TABLE]

for $t\in\mathbb{R}$. Since $H_{1}(dt\mid x)=F_{C}([t,\infty)\mid x)F_{T}(dt\mid x)$, using relationship (2.8) we obtain

[TABLE]

Next, using the product-integral we can write

[TABLE]

Let us recall that $F_{C}(\cdot\mid x)$ can be written as a transformation of $H_{0}(\cdot\mid x)$ and $H_{1}(\cdot\mid x)$, see equations (2.4) and (2.2). 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 $F_{C}(\cdot\mid 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 $F_{C}(\cdot\mid x)$.

Finally, the representation of $F_{C}(\cdot\mid x)$ given in equation (2.2), plugged into equation (2.9), allows to express $\Lambda_{T,0}(\cdot\mid x)$, and thus $F_{T,0}(\cdot\mid x)$, as maps of $\mathbb{P}(T=\infty\mid x)$ and the measures $H_{0}(\cdot\mid x)$ and $H_{1}(\cdot\mid 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 $F_{T}(\cdot\mid x)$, $F_{C}(\cdot\mid x)$ and $\phi(\cdot,\beta)$, and the assumptions (2.1), (2.2) and (2.6). For a fixed value of the parameter $\beta,$ and for $t\in\mathbb{R}$ and $x\in\mathcal{X}$, let

[TABLE]

and

[TABLE]

Let $F_{Y,\delta}(\cdot,\cdot\mid x)$ denote the conditional law of $(Y,\delta)$ given $X=x$. Moreover, let

[TABLE]

and

[TABLE]

These equations define a conditional law for the observations $(Y,\delta)$ based on the model. More precisely, for a choice of $F_{T}(\cdot\mid x)$, $F_{C}(\cdot\mid x)$ and $\beta$, the model yields a conditional law for $(Y,\delta)$ given $X=x$. If the model is correctly specified, there exists a value $\beta_{0}$ such that

[TABLE]

The remaining question is whether the true value of the parameter is identifiable. In other words, one should check if, given the conditional subdistributions $H_{0}(\cdot\mid x)$ and $H_{1}(\cdot\mid x),$ $x\in\mathcal{X},$ there exists a unique $\beta_{0}$ satisfying condition (2.13). For this purpose we impose the following mild condition:

[TABLE]

and we show that

[TABLE]

Indeed, if $F_{Y,\delta}^{\beta_{0}}(\cdot,\cdot\mid x)=F_{Y,\delta}^{\widetilde{\beta}}(\cdot,\cdot\mid x),$ then for any $x$,

[TABLE]

for all $t\in(-\infty,\tau_{H}(x)]\cap\mathbb{R}.$ Our condition (2.7) guarantees that $\tau_{H}(x)=\tau_{H_{0}}(x)$ so that $F_{C}((t,\infty)\mid x)$ should be necessarily positive for $t\in(-\infty,\tau_{H}(x)),$ and thus could be simplified in the last display. Deduce that

[TABLE]

for all $t\in(-\infty,\tau_{H}(x))$. Finally, recall that by construction, in the model we consider, $F_{T,0}^{\beta}((-\infty,\infty)\mid x)=1$ for any $\beta$ such that $F_{Y,\delta}^{\beta}(\cdot,\cdot\mid x)$ coincides with the conditional law of $(Y,\delta)$ given $X=x$, for any $x$. Thus taking integrals on $(-\infty,\infty)$ on both sides of the last display we obtain $\phi(\cdot,\beta_{0})=\phi(\cdot,\widetilde{\beta})$. Let us gather these facts in the following statement.

Theorem 2.1**.**

Under conditions (2.1), (2.2), (2.6) and (2.14) the model is identifiable.

2.5 Interpreting the previous modeling approaches

We suppose here that the function $\phi(x)$ follows a logistic model, and comment on several models for $F_{T,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 $\tau_{H_{0}}(x)=\tau_{H_{1}}(x)=\infty$ and that $\Lambda_{T,0}(\cdot\mid 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 $\Lambda_{T,0}(\cdot\mid x)$ measure. More precisely, it is supposed that

[TABLE]

where $\gamma$ is some parameter to be estimated and $\Lambda_{0}(\cdot)$ is an unknown baseline cumulative hazard function. Our inversion formulae reveal that in this approach the parameters $\gamma$ and $\Lambda_{0}$ depend on the observed conditional measures $H_{0}(\cdot\mid x)$ and $H_{1}(\cdot\mid x),$ but also on the parameter $\beta$. The same is true for the parametric models.

2.5.2 Kaplan-Meier mixture cure model

Taylor (1995) suggested to estimate $F_{T,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<\infty$ does not depend on $X$. This is equivalent to supposing that $\Lambda_{T,0}(\cdot\mid x)=\Lambda_{T,0}(\cdot)$. Next, to estimate $\Lambda_{T,0}(\cdot)$ 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=\infty\}$. Following Taylor’s approach we rewrite (2.9) as

[TABLE]

Next, assume that the last equality remains true if $H_{0}(dt\mid x)$ and $H_{1}(dt\mid x)$ are replaced by their unconditional versions, that is assume that

[TABLE]

See equations (2) and (3) in Taylor (1995). The equation above could be solved iteratively by a EM-type procedure: for a given $\beta$ and an iteration $F_{T,0}^{(m)}(\cdot),$ build $\Lambda_{T,0}^{(m+1)}(dt)$ and the updated estimate $F_{T,0}^{(m+1)}(\cdot)$; see Taylor (1995) for the details. Let us point out that even if $(T,C)$ is independent of $X$ and thus $H_{1}(\cdot\mid x)$ does not depend on $x$, the subdistribution $H_{0}(\cdot\mid x)$ still depends on $x$, since

[TABLE]

Hence, a more natural form of equation (2.15) is

[TABLE]

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

3 Maximum likelihood estimation

Let $(Y_{i},\delta_{i},X_{i})$ ($i=1,\ldots,n$) be a sample of $n$ i.i.d. copies of the vector $(Y,\delta,X)$.

We use a likelihood approach based on formulae (2.9) and (2.2) to build an estimator of $\phi(\cdot,\beta)$ and $F_{T,0}^{\beta}(\cdot\mid x).$ To build the likelihood we use estimates $\widehat{H}_{k}(\cdot\mid x)$ of the subdistributions $H_{k}(\cdot\mid x),$ $k\in\{0,1\}.$ These estimates are constructed with the sample of $(Y,\delta,X),$ without reference to any model for the conditional probability $\mathbb{P}(T<\infty\mid x).$ At this stage it is not necessary to impose a particular form for $\widehat{H}_{k}(\cdot\mid x)$. To derive the asymptotic results we will only impose that these estimators satisfy some mild conditions. Let $\widehat{F}_{T,0}^{\beta}(\cdot\mid x)$ be defined as in equations (2.11) and (2.12) with $\widehat{H}_{0}(\cdot\mid x)$ and $\widehat{H}_{1}(\cdot\mid x)$ instead of $H_{0}(\cdot\mid x)$ and $H_{1}(\cdot\mid x)$, that is

[TABLE]

where $\widehat{F}_{C}(\cdot\mid x)$ is the estimator obtained from equations (2.4) and (2.2) but with $\widehat{H}_{0}(\cdot\mid x)$ and $\widehat{H}_{1}(\cdot\mid x)$, i.e.

[TABLE]

Let $f_{X}$ denote the density of the covariate vector with respect to some dominating measure. The contribution of the observation $(Y_{i},\delta_{i},X_{i})$ to the likelihood when $\delta_{i}=1$ is then

[TABLE]

while the contribution when $\delta_{i}=0$ is

[TABLE]

Since the laws of the censoring variable and of the covariate vector do not carry information of the parameter $\beta$, we can drop the factors $\widehat{F}_{C}((Y_{i},\infty)\mid X_{i})f_{X}(X_{i})$ and $\widehat{F}_{C}(\{Y_{i}\}\mid X_{i})f_{X}(X_{i}).$ Hence the criterion to be maximized with respect to $\beta\in B$ is $\widehat{L}_{n}(\beta)$ where

[TABLE]

The estimator we propose is

[TABLE]

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=\infty\}$ is correct and identifiable in the sense of condition (2.14), the true parameter $\beta_{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 $\beta_{0}.$

Let

[TABLE]

with

[TABLE]

[TABLE]

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 $\beta,$ the function $\beta\mapsto\mathbb{E}[\log p(Y,\delta,X;\beta)]$ is expected to be the limit of the random function $\log\widehat{L}_{n}(\cdot).$ Hence, a minimal condition for guaranteeing the consistency of the likelihood estimation approach is that $\beta_{0}$ is the maximizer of the limit likelihood criterion $\mathbb{E}[\log p(Y,\delta,X;\cdot)]$. 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 $\beta_{0}$ is the value of the parameter defined by equation (2.13), then for any $\beta\neq\beta_{0}$,

[TABLE]

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 consistency, we have to control the asymptotic behavior of $\widehat{L}_{n}(\beta)$ along sequences of values of the parameter $\beta.$ Such a control requires a control of denominators like

[TABLE]

on the support of $H_{1}(\cdot\mid x),$ uniformly with respect to $x.$ A usual way to deal with this technical difficulty is to consider a finite threshold $\tau$ beyond which no uncensored lifetime is observed, i.e.

[TABLE]

Moreover, to be able to keep denominators away from zero, we require the condition

[TABLE]

In particular, this condition implies

[TABLE]

Moreover, given condition (2.2), necessarily $H_{0}(\{\tau\}\mid x)=0$, $\forall x.$ This means that $F_{C}(\{\tau\}\mid x)=0$, $\forall x.$ This constraint on $H_{0}(\{\tau\}\mid 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}F_{C}((\tau,\infty)\mid x)>0,$ and $\inf_{x}F_{T,0}^{\beta}(\{\tau\}\mid x)>0$, $\forall\beta.$

Conditions like in equations (4.1)-(4.2) are more or less explicitly used in the literature of cure models. Sometimes $\tau$ is justified as representing a total follow-up of the study. For instance, Lu (2008) supposes that $Y=\min\{T,\min(C,\tau)\}$ and $\delta=\mathbf{1}\{T\leq\min(C,\tau)\},$ where $T=\eta T^{*}+(1-\eta)\infty,$ with $T^{*}<\infty$ and $\eta\in\{0,1\}.$ The conditional probability of being cured is precisely the conditional probability of the event $\{\eta=0\}.$ Next, Lu (2008) supposes that $\inf_{x}\mathbb{P}(\tau\leq T\leq C\mid x)>0,$ and $\Lambda_{0}(\tau)<\infty,$ where $\Lambda_{0}(\cdot)$ 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,\tau]$ and suppose $\mathbb{E}(\delta\mathbf{1}\{Y\geq\tau\})>0.$ This could be possible only if $H_{1}(\{\tau\}\mid 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 $H_{1}(\{\tau\}\mid x)\geq 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.

4.1 Consistency

Let us sketch the arguments we use in the proof of Theorem 4.1 below for deriving the consistency of $\widehat{\beta}$. On one hand, if the conditional subdistributions $H_{k}(\cdot\mid x)$ are given, one can build the purely parametric likelihood

[TABLE]

By construction, $\mathcal{L}_{n}(\beta)$ is a functional of $H_{0}(\cdot\mid x)$ and $H_{1}(\cdot\mid x),$ $x\in\mathcal{X},$ while $\widehat{L}_{n}(\beta)$ is a functional of the estimated versions of $H_{0}(\cdot\mid x)$ and $H_{1}(\cdot\mid x).$ Hence, a prerequisite condition for deriving the consistency of our semiparametric estimator $\widehat{\beta}$ is the consistency of the infeasible maximum likelihood estimator

[TABLE]

A necessary condition for the consistency of $\widehat{\beta}$ is

[TABLE]

We then have that

[TABLE]

From this we will derive the consistency of $\widehat{\beta}$ using Section 5.2 in van der Vaart (1998). See the proof in the Appendix for details.

To prove condition (4.4), we have to guarantee the uniform convergence of $\widehat{H}_{k}-H_{k}$ as stated in Assumption (AC1) below. Indeed, this uniform convergence will imply

[TABLE]

and

[TABLE]

See Lemma 7.1 in the Appendix. The uniform convergence in equation (4.4) then follows.

To prove the consistency of $\widehat{\beta}$, we need the following assumptions :

• (AC1)

For $\tau$ appearing in conditions (4.1)-(4.2),

[TABLE]

• (AC2)

The parameter set $B\subset\mathbb{R}^{p}$ is compact.

• (AC3)

There exist some constants $a>0$ and $c_{1}>0$ such that

[TABLE]

• (AC4)

$\inf_{\beta\in B}\inf_{x\in\mathcal{X}}\phi(x,\beta)>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 $\beta_{0}$ in the parameter space $B$ such that (2.13) is true. Then, $\widehat{\beta}-\beta_{0}=o_{\mathbb{P}}(1).$

Let us point out that the consistency result is stated in terms of the subdistributions of the observations and the conditional probability model $\{\phi(x,\cdot):\beta\in B\}.$ If the identification assumptions used in Proposition 3.1 hold true and the model is correctly specified, $\phi(\widehat{\beta},x)$ consistently estimates the cure probability $\mathbb{P}(T=\infty\mid 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 $\phi(x,\beta)$ and is e.g. satisfied for the logistic model if $\mathcal{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 purpose we use the derivative of $\log\widehat{L}_{n}(\beta)$ with respect to $\beta$.

First note that the vector of partial derivatives of the log-likelihood $\log\widehat{L}_{n}(\beta)$ with respect to the components of $\beta$ equals

[TABLE]

where $q_{i}$ is defined in the proof of Theorem 4.1.

To develop the asymptotic normality of our estimator $\widehat{\beta}$, we embed the nuisance functions $H_{k}([\cdot,\infty)\mid\cdot)$ ($k=0,1$) in a functional space ${\cal H}$, which is equipped with a pseudo-norm $\|\cdot\|_{\cal H}$. Both the space ${\cal H}$ and its pseudo-norm $\|\cdot\|_{\cal H}$ will be chosen depending on the estimators $\widehat{H}_{k}([\cdot,\infty)\mid\cdot)$, $k=0,1$, and have to satisfy certain conditions, which we give below. The true vector of nuisance functions is

[TABLE]

For each $x\in{\cal X}$ and for each $\eta\in{\cal H}$, let $\eta(dt,x)=(\eta_{1}(dt,x),\eta_{2}(dt,x))$ be the measures associated to the non-increasing functions $\eta(\cdot,x)=(\eta_{1}(\cdot,x),\eta_{2}(\cdot,x))$, and define

[TABLE]

and

[TABLE]

where

[TABLE]

where

[TABLE]

Note that for $\eta=\eta_{0}$, we have

[TABLE]

$T_{2}(\eta_{0})(t,x)=F_{C}([t,\infty)\mid x)$ and $T_{3}(\beta,\eta_{0})(t,x)=F_{T,0}^{\beta}([t,\infty)\mid x).$ Hence, we have that

[TABLE]

where $\widehat{\eta}(t,x)=(\widehat{\eta}_{1}(t,x),\widehat{\eta}_{2}(t,x))=(\widehat{H}_{0}([t,\infty)\mid x),\widehat{H}_{1}([t,\infty)\mid x))$.

If in addition $\beta=\beta_{0}$,

[TABLE]

and thus for any $x$, the map $t\mapsto T_{1}(\beta_{0},\eta_{0})(t,x)$ is decreasing on $(-\infty,\tau]$. Moreover, by condition (4.2),

[TABLE]

We need this lower bound to be valid on a neighborhood around $\beta_{0}$. Hence, let us consider $B_{0}$ a neighborhood of $\beta_{0}$ such that

[TABLE]

The existence of $B_{0}$ is guaranteed by condition (4.2) and the regularity of the function $\phi(\cdot,\cdot)$; see assumption (AN3) below that strengthens assumption (AC3). Finally, let us note that by construction for any $t\in(-\infty,\tau)$, $H(\{t\}\mid x)=H_{0}(\{t\}\mid x)+H_{1}(\{t\}\mid x)\geq(1-\phi(x,\beta))F_{C}(\{t\}\mid x)+H_{1}(\{t\}\mid x)$ and thus

[TABLE]

Then, by the arguments guaranteeing the existence of a set $B_{0}$ as in equation (4.10),

[TABLE]

Further, define the Gâteaux derivative of $M(\beta,\eta_{0})$ in the direction $[\eta-\eta_{0}]$ by

[TABLE]

and in a similar way the Gâteaux derivatives $\nabla_{\eta}T_{j}(\beta,\eta_{0})[\eta-\eta_{0}]$ are defined.

We need the following assumptions :

• (AN1)

The matrix $\nabla_{\beta}M(\beta,\eta)$ exists for $\beta$ in a neighborhood $B_{0}$ of $\beta_{0}$ and is continuous in $\beta$ for $\beta=\beta_{0}$. Moreover, $\nabla_{\beta}M(\beta_{0},\eta_{0})$ is non-singular.

• (AN2)

$H_{k}([\cdot,\infty)\mid\cdot)\in{\cal H}$ for $k=0,1$.

• (AN3)

The function $\beta\rightarrow\phi(x,\beta)$ is continuously differentiable for all $x\in{\cal X}$, and the derivative is bounded uniformly in $x\in{\cal X}$ and $\beta\in B_{0}$. Moreover, $B_{0}$ is compact and $\beta_{0}$ belongs to the interior of $B_{0}$.

• (AN4)

For $k=0,1$, the estimator $\widehat{H}_{k}([\cdot,\infty)\mid\cdot)$ satisfies the following :

• $(i)$

$\mathbb{P}(\widehat{H}_{k}([\cdot,\infty)\mid\cdot)\in{\cal H})\rightarrow 1$

• $(ii)$

$\|(\widehat{H}_{k}-H_{k})([\cdot,\infty)\mid\cdot)\|_{\cal H}=o_{\mathbb{P}}(n^{-1/4})$

• $(iii)$

There exist functions $\Psi_{1}$ and $\Psi_{2}$, such that

[TABLE]

and

[TABLE]

where $\mathbb{E}^{*}$ denotes conditional expectation given the data $(Y_{i},\delta_{i},X_{i}),$ $1\leq i\leq n$, the functions $\psi_{jk}$ are defined in (7.13) in the Appendix, and where

[TABLE]

($j=1,\ldots,5;k=0,1;\ell=1,2$). Note that the above expectations are conditionally on the sample and are taken with respect to the generic variables $Y,\delta,X$ which have the same law as the sample.

• (AN5)

The class ${\cal H}$ satisfies $\int_{0}^{\infty}\sqrt{\log N(\varepsilon,{\cal H},\|\cdot\|_{\cal H})}d\varepsilon<\infty$, where $N(\varepsilon,{\cal H},\|\cdot\|_{\cal H})$ is the $\varepsilon$-covering number of the space ${\cal H}$ with respect to the norm $\|\cdot\|_{\cal H}$, i.e. the smallest number of balls of $\|\cdot\|_{\cal H}$-radius $\varepsilon$ needed to cover the space ${\cal H}$.

Theorem 4.2**.**

Assume that $\widehat{\beta}-\beta_{0}=o_{\mathbb{P}}(1)$ and that (AN1)-(AN5) and (2.1), (2.2), (2.6), (2.14), (4.1) and (4.2) hold true. Then,

[TABLE]

where

[TABLE]

and V=\mbox{Var}\big{(}m(Y,\delta,X;\beta_{0},\eta_{0})+\Psi_{1}(Y,\delta,X)+\Psi_{2}(Y,\delta,X)\big{)}.

4.3 Bootstrap consistency

Although in principle one can use Theorem 4.2 above for making inference, the asymptotic variance $\Omega$ has a complicated structure, and the estimation of $\Omega$ would not only be cumbersome, but its precision for small samples could moreover be rather poor. We continue this section by showing that a bootstrap procedure can be used to estimate the asymptotic variance of $\widehat{\beta}$, to approximate the whole distribution of $\widehat{\beta}$ or to construct confidence intervals or test hypotheses regarding $\beta_{0}$.

Here, we propose to use a naive bootstrap procedure, consisting in drawing triplets $(Y_{i}^{*},\delta_{i}^{*},X_{i}^{*}),$ $1\leq i\leq n,$ randomly with replacement from the data $(Y_{i},\delta_{i},X_{i}),$ $1\leq i\leq n$. Let $\widehat{H}_{k}^{*}$ be the same estimator as $\widehat{H}_{k}$ ($k=0,1$) but based on the bootstrap data, and for each $(\beta,\eta)$ let $M_{n}^{*}(\beta,\eta)=n^{-1}\sum_{i=1}^{n}m(Y_{i}^{*},\delta_{i}^{*},X_{i}^{*};\beta,\eta)$. Define the bootstrap estimator $\widehat{\beta}^{*}$ to be any sequence that satisfies

[TABLE]

where $\widehat{\eta}^{*}(t,x)=(\widehat{\eta}_{1}^{*}(t,x),\widehat{\eta}_{2}^{*}(t,x))=(\widehat{H}_{0}^{*}([t,\infty)\mid x),\widehat{H}_{1}^{*}([t,\infty)\mid x))$.

The following result shows that the bootstrap works, in the sense that it allows to recover correctly the distribution of $n^{1/2}(\widehat{\beta}-\beta_{0})$.

Theorem 4.3**.**

Assume that $\widehat{\beta}-\beta_{0}=o_{\mathbb{P}}(1)$ and that (AN1)-(AN5) hold true. Moreover, assume that $\nabla_{\beta}M(\beta,\eta)$ is continuous in $\eta$ (with respect to $\|\cdot\|_{\cal H}$) at $(\beta,\eta)=(\beta_{0},\eta_{0})$, and that (AN4) holds true with $\widehat{H}_{k}-H_{k}$ replaced by $\widehat{H}_{k}^{*}-\widehat{H}_{k}$ ($k=0,1$) in ${\mathbb{P}}^{*}$-probability. Then,

[TABLE]

where $\mathbb{P}^{*}$ denotes probability conditionally on the data $(Y_{i},\delta_{i},X_{i})$, $i=1,\ldots,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 $H_{k}$ are estimated by means of kernel smoothing.

Consider the case where $X$ is composed of continuous and discrete components, that is $X=(X_{c},X_{d})\in\mathcal{X}_{c}\times\mathcal{X}_{d}\subset\mathbb{R}^{d_{c}}\times\mathbb{R}^{d_{d}},$ with $d_{c}+d_{d}=d\geq 1.$ For simplicity, assume that the support $\mathcal{X}_{d}$ of the discrete subvector $X_{d}$ 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,\infty]$. The subdistributions $H_{k}([t,\infty)\mid x)$ could be estimated by means of a kernel estimator :

[TABLE]

where for any $(x_{c},x_{d})\in\mathcal{X}_{c}\times\mathcal{X}_{d},$

[TABLE]

$h_{n}$ is a bandwidth sequence, $K_{h}(\cdot)=K(\cdot/h)/h^{d_{c}}$, $K(u)=k(u_{1})\cdot\ldots\cdot k(u_{d_{c}})$ and $k$ is a probability density function.

Nonparametric smoothing of continuous covariates is possible for dimensions $d_{c}$ 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 $d_{c}=1$. The discrete covariates do not contribute to the curse of dimensionality, and therefore $d_{d}$ 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 :

• (C1)

The sequence $h_{n}$ satisfies $nh_{n}^{4}\rightarrow 0$ and $nh_{n}^{3+\zeta}(\log n)^{-1}\rightarrow\infty$ for some $\zeta>0$.

• (C2)

The support ${\cal X}$ of $X$ is a compact subset of $\mathbb{R}$.

• (C3)

The probability density function $K$ has compact support, $\int uK(u)du=0$ and $K$ is twice continuously differentiable.

Further, let ${\cal F}_{1}$ be the space of functions from $[0,\tau]$ to $[0,1]$ with variation bounded by $M$, and let ${\cal F}_{2}$ be the space of continuously differentiable functions $f$ from ${\cal X}$ to $[-M,M]$ that satisfy $\sup_{x\in{\cal X}}|f^{\prime}(x)|\leq M$ and $\sup_{x_{1},x_{2}\in{\cal X}}|f^{\prime}(x_{1})-f^{\prime}(x_{2})|/|x_{1}-x_{2}|^{\epsilon}\leq M$ for some $M<\infty$ and $0<\epsilon<1$. Let

[TABLE]

We define the following norm associated with the space ${\cal H}$ : for $\eta\in{\cal H}$, let

[TABLE]

Then, it follows from Propositions 1 and 2 in Akritas and Van Keilegom (2001) that ${\mathbb{P}}(\widehat{H}_{k}\in{\cal H})\rightarrow 1$ provided $nh_{n}^{3+\zeta}(\log n)^{-1}\rightarrow\infty$, with $\zeta>0$ as in condition (C1). Moreover, $\sup_{x,t}|\widehat{H}_{k}([t,\infty)\mid x)-H_{k}([t,\infty)\mid x)|=O_{\mathbb{P}}((nh_{n})^{-1/2}(\log n)^{1/2})=o_{\mathbb{P}}(n^{-1/4})$ (see Proposition 1 in Akritas and Van Keilegom 2001). The class ${\cal 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 left hand side equals

[TABLE]

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 $\phi(x,\beta)$ of not being cured follows a logistic model :

[TABLE]

for any $-1\leq x\leq 1$. We will work with $\beta_{0}=(\beta_{01},\beta_{02})=(1.75,2)$ and $(1.1,2)$, corresponding to an average cure rate of 20$\%$ respectively 30$\%$. The conditional distribution function $F_{T,0}(\cdot|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[-(\gamma_{0}+\gamma_{1}x+\gamma_{2}/(1+2x^{2}))]$, where $\gamma_{0}=\gamma_{1}=0.5$ and $\gamma_{2}\in\{0,1,2\}$. Next, in order to respect condition (4.2), we truncate this distribution at $\tau$, which is the quantile of order 0.97 of an exponential distribution with mean $\mathbb{E}\{\exp[-(\gamma_{0}+\gamma_{1}X+\gamma_{2}/(1+2X^{2}))]\}$, i.e.

[TABLE]

Note that this is the distribution function corresponding to a Cox model with baseline hazard equal to $I(0\leq t\leq\tau)$, and exponential factor equal to $\exp(\gamma_{0}+\gamma_{1}x+\gamma_{2}/(1+2x^{2}))$.

Next, we generate the censoring variable $C$ independently of $X$ from an exponential distribution with mean equal to 1.65 when $\beta_{0}=(1.75,2)$, and with mean 1.45 when $\beta_{0}=(1.1,2)$. In this way we have respectively 40$\%$ and 50$\%$ of censoring when $\gamma_{2}=0$.

In what follows we will compare our estimator of $\beta$ 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 $\gamma_{2}=0$. The estimated $\beta$ 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 $\widehat{\beta}$ 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 $\beta$ from a logistic model based on the nonparametric estimator $\widehat{F}_{T}([0,\infty)|x)$, so we maximize the log-likelihood $\sum_{i=1}^{n}\{(1-\widehat{F}_{T}([0,\infty)|X_{i}))\log(\phi(X_{i},\beta))+\widehat{F}_{T}([0,\infty)|X_{i})\log(1-\phi(X_{i},\beta))\}$. 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 $\beta_{0}=(1.75,2)$, and in Table 2 for the case where $\beta_{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 $\widehat{\beta}_{1}$ and $\widehat{\beta}_{2}$ obtained under the Cox model and from our procedure. The kernel function $K$ is taken equal to the Epanechnikov kernel : $K(u)=(3/4)(1-u^{2})I(|u|\leq 1)$. The bandwidth $h$ of the kernel estimators $\widehat{H}_{k}(\cdot|X_{i})$ ($k=0,1$) 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 sample in our simulation, we calculated these bandwidths for $\widehat{H}_{0}$ and $\widehat{H}_{1}$ and used the average 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 $F_{T,0}(\cdot|x)$ when $x=0.25$. We estimate these quartiles by means of our nonparametric estimator $\widehat{F}_{T,0}^{\widehat{\beta}}(\cdot|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 $\gamma_{2}=1$ or 2), the 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 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 $\widehat{\beta}_{1}$ and $\widehat{\beta}_{2}$ 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 $\widehat{\beta}_{1}$. 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 $\widehat{\beta}_{1}$ and $\widehat{\beta}_{2}$ 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 $\beta_{1}$ and $\beta_{2}$ is also added, and shows that the bootstrap variance is well centered around the corresponding empirical variance.

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 $\beta$ 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 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 $\phi(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-F_{T,0}(\cdot|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. 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 $\phi(x)$ are quite different, and have opposite slopes. Moreover, the survival function $1-F_{T,0}(\cdot|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-\phi(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.

7 Appendix : Proofs

Proof of Proposition 3.1..

By the properties of the likelihood of a Bernoulli random variable given that $Y\in dt$ and $X=x$, we have

[TABLE]

Integrate with respect to $Y$ and $X$ and deduce that

[TABLE]

If there exists some $\beta\neq\beta_{0}$ such that the last inequality becomes an equality, then necessarily $p_{1}(Y,X;\beta)=1$ almost surely. Then

[TABLE]

for almost all $x\in\mathcal{X}.$ By Theorem 2.1, we deduce that necessarily $\beta=\beta_{0}$. ∎

Lemma 7.1**.**

Let conditions (4.2), (AC1) and (AC4) hold true. Then,

[TABLE]

[TABLE]

where

[TABLE]

and $T_{1}(\beta,\widehat{H}_{0},\widehat{H}_{1})$ is defined similarly, but with $H_{0},$ $H_{1}$ and $F_{C}$ replaced by $\widehat{H}_{0},$ $\widehat{H}_{1}$ and $\widehat{F}_{C},$ respectively. Moreover,

[TABLE]

Proof of Lemma 7.1.

Let us first investigate the uniform convergence of the estimated cumulative hazard measure $\widehat{\Lambda}_{C}(\cdot\mid x)$. For any $t\in(-\infty,\tau]$ let us write

[TABLE]

The integrals with respect to $\widehat{H}_{0}(ds\mid x)$ are well defined, since, with probability tending to 1, for each $x\in\mathcal{X}$, the map $s\mapsto\widehat{H}_{0}((-\infty,s]\mid x),$ $s\in(-\infty,\tau],$ is a function of bounded variation. The uniform convergence Assumption (AC1) implies that

[TABLE]

for some constant $c>0.$ Next, by Duhamel’s identity (see Gill and Johansen 1990),

[TABLE]

Then, the uniform convergence of $\widehat{F}_{C}(\cdot\mid x)$ follows from the uniform convergence of $\widehat{\Lambda}_{C}(\cdot\mid x)$ and condition (4.2). The same type of arguments apply for $T_{1}(\beta,\widehat{H}_{0},\widehat{H}_{1})$, and hence we omit the details.

Next, since by conditions (4.2) and (AC4) we have

[TABLE]

there exists some constant $c>0$ with the property that

[TABLE]

Hence, the uniform convergence of $\widehat{F}_{T,0}^{\beta}(\cdot\mid x)$ follows. ∎

Proof of Theorem 4.1.

Let us write

[TABLE]

where

[TABLE]

Moreover, let

[TABLE]

where, for $t\in\mathbb{R},\;d\in\{0,1\},\;x\in\mathcal{X},$

[TABLE]

Let

[TABLE]

Similarly, let us consider $Q_{n}(\beta,\widehat{H}_{0},\widehat{H}_{1})$ that is defined as $Q_{n}(\beta,H_{0},H_{1}),$ but with $H_{0},$ $H_{1},$ $F_{C}$ and $F_{T,0}^{\beta}$ replaced by $\widehat{H}_{0},$ $\widehat{H}_{1},$ $\widehat{F}_{C}$ and $\widehat{F}_{T,0}^{\beta}$, respectively. Then the estimator $\widehat{\beta}$ in equation (3.1) becomes

[TABLE]

The first step is to check that

[TABLE]

This follows directly from Lemma 7.1. Next, given our assumptions, it is easy to check that for any $t\in(-\infty,\tau],$ $d\in\{0,1\},$ $x\in\mathcal{X},$

[TABLE]

with $a>0$ from Assumption (AC3) and some constant $C$ depending only on $c_{1}$ from Assumption (AC3) and the positive values $\inf_{x\in\mathcal{X}}H_{1}(\{\tau\}\mid x)$ and $\inf_{x\in\mathcal{X}}H_{0}((\tau,\infty)\mid x).$ It follows that the class $\{(t,d,x)\rightarrow q(\beta,H_{0},H_{1})(t,d,x):\beta\in B\}$ is Glivenko-Cantelli. Hence,

[TABLE]

where $Q=\mathbb{E}(Q_{n})$. Finally, Proposition 3.1 guarantees that

[TABLE]

Gathering the facts, we deduce that $\widehat{\beta}-\beta_{0}=o_{\mathbb{P}}(1).$ ∎

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 $\nabla_{\eta}M(\beta,\eta_{0})[\eta-\eta_{0}]$ is linear in $\nabla_{\eta}T_{j}(\beta,\eta_{0})[\eta-\eta_{0}]$ ($j=1,2,3,4$). Hence, we will focus attention on the latter Gâteaux derivatives. First,

[TABLE]

Using Duhamel’s formula (see Gill and Johansen 1990), we can write

[TABLE]

In a similar way, we find that

[TABLE]

Finally,

[TABLE]

Note that all denominators in $\nabla_{\eta}T_{j}(\beta,\eta_{0})[\eta-\eta_{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

[TABLE]

for some constant $C$. Hence, it can be easily seen that $\nabla_{\eta}T_{j}(\beta,\eta_{0})[\eta-\eta_{0}]$ satisfies the second property in assumption (2.3) in Chen et al. (2003), and hence the same holds true for $\nabla_{\eta}M(\beta,\eta_{0})[\eta-\eta_{0}]$. Similarly, by decomposing $T_{j}(\beta,\eta)-T_{j}(\beta,\eta_{0})-\nabla_{\eta}T_{j}(\beta,\eta_{0})[\eta-\eta_{0}]$ using Taylor-type arguments (in $\eta$), 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 $\nabla_{\eta}T_{j}(\beta,\eta_{0})[\eta-\eta_{0}]$ ($j=1,2,3,4$) that

[TABLE]

for certain measurable functions $\psi_{jk}$ ($j=1,\ldots,5;k=0,1$).

It remains to verify condition (2.5). Note that

[TABLE]

for some functions $C_{j}$ satisfying ${\mathbb{E}}[C_{j}^{2}(Y,\delta,X)]<\infty$ ($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 $\mathbb{P}$-probability, whereas their result holds true a.s. $[\mathbb{P}]$. As a consequence, in their high level conditions we can replace all a.s. $[\mathbb{P}]$ statements by the corresponding statements in $\mathbb{P}$-probability.

First of all, it follows from assumption (AN1) that condition (2.2) in Chen et al. (2003) holds with $\eta_{0}$ replaced by any $\eta$ in a neighborhood of $\eta_{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 true if we replace $\widehat{H}_{k}-H_{k}$ by $\widehat{H}_{k}^{*}-\widehat{H}_{k}$ ($k=0,1$). It remains to verify condition (2.5’B) in Chen et al. This follows from Theorem 3 in Chen et al., whose conditions have been verified already for our Theorem 4.2. ∎