Sylvain Takerkart

Guillaume Auzias

Bertrand Thirion

Daniele Schön

Liva Ralaivola

Graph-based inter-subject classification of local fMRI patterns

Keywords: fMRI, classification, graphs, kernels, inter-subject, variability

Classification of medical images in multi-subjects settings is a difficult challenge due to the variability that exists between individuals. Here we introduce a new graph-based framework designed to deal with inter-subject functional variability present in fMRI data. A graphical model is constructed to encode the functional, geometric and structural properties of local activation patterns. We then design a specific graph kernel, allowing to conduct SVM classification in graph space. Experiments conducted in an inter-subject classification task of patterns recorded in the auditory cortex show that it is the only approach to perform above chance level, among a wide range of tested methods.

Introduction

Since brain functions are considered to arise from activities distributed across networks, multivariate pattern recognition methods are adapted to study information processing in the brain. Multi-voxel pattern analysis (MVPA) of functional MRI allows to study the organization of distributed representations in the brain by using a classification framework where one attempts to predict the category of the input stimuli from the data [START_REF] Haxby | Distributed and overlapping representations of faces and objects in ventral temporal cortex[END_REF]. MVPA has been largely applied within individual subjects [START_REF] Norman | Beyond mind-reading: multi-voxel pattern analysis of fMRI data[END_REF].

However, because of the challenge posed by the large inter-individual variability, very few studies describe inter-subject decoding, i.e successfull prediction on data from a subject that was not part of the training set. Among these, most rely on global (i.e full-brain) analysis, using large-scale features ( [START_REF] Michel | A supervised clustering approach for fMRI-based inference of brain states[END_REF], [START_REF] Poldrack | Decoding the Large-Scale structure of brain function by classifying mental states across individuals[END_REF]). The few studying fine-scale local patterns achieve very low inter-subject generalization performances [START_REF] Clithero | Within-and crossparticipant classifiers reveal different neural coding of information[END_REF], or use implicit abstract models [START_REF] Haxby | A common, High-Dimensional model of the representational space in human ventral temporal cortex[END_REF]. None attempt to characterize the inter-subject functional variability (i.e the fact that the correlation between cortical folding and the underlying functional organization vary between subjects [START_REF] Van Essen | Surface-Based and probabilistic atlases of primate cerebral cortex[END_REF]) and use such characterization in the classification process.

In this paper, we introduce a new framework specifically aimed at tackling the challenges offered by inter-subject functional variability by modelling its spatial properties. Our approach integrates the fact that the geometric properties of local functional features, as well as their levels of activation, can vary across subjects, under the assumption that the underlying spatial structure of the local activation pattern is consistent. We therefore design a graphical model to represent such patterns with their properties of interest: the nodes of the graphs represent small activation patches; their attributes carry the relevant features, such as their position and activation level; the edges of the graph (given by spatial adjacency) encode the spatial structure of the pattern. The classification is then performed directly in graph-space with Support Vector Machines (SVM), using a graph kernel specifically designed to take into account all properties of our graphical model. Since the topographic organization of primary sensory areas ensures the structural consistency of activation patterns across subjects, we validate our framework on fMRI data recorded in the primary auditory cortex during a tonotopy experiment [START_REF] Humphries | Tonotopic organization of human auditory cortex[END_REF].

This paper is organized as follows: we detail the construction of our graphical model in section 2 and the kernel design in section 3; we then describe our experiments in section 4 and present our results in section 5.

Graphical Model of Activation Patterns

We here describe the design of a graphical model that summarizes the information relevant to characterize a pattern a pattern of activation measured with fMRI within a contiguous region of interest (ROI).

Graph nodes. Assuming that the ROI admits an underlying subdivision into a set of smaller and functionally relevant sub-regions, the first step to construct our graphical model estimates a parcellation of the ROI, i.e. a partition into a set of sub-regions or parcels [START_REF] Flandin | Improved detection sensitivity of functional MRI data using a brain parcellation technique[END_REF]. Specifically, for a contiguous ROI R, we compute the parcellation V = {V i } q i=1 of R, so that the q parcels verify: ∪ q i=1 V i = R and V i ∩ V j = ∅ whenever i = j. We use V as the set of vertices of our graphical model, each parcel corresponding to a node of the graph.

Nodes attributes: functional features. Let f be the real-valued function describing the Bold activation. In a parcel V i , the activations values {f (v)} v∈Vi are summarized to form an n-dimensional feature vector F (V i ); we note

F the n × q matrix F = [F (V 1 ) • • • F (V q )]
of functional features, where each column i of F is F (V i ). A simple example of such F , which we make use of later on, is to compute the mean activation value within a parcel.

Nodes attributes: geometric features. Let x be a coordinate system defined in R. We summarize the geometric information of parcel V i through a mdimensional feature vector X(V i ) , computed from the locations {x(v)} v∈Vi . We note

X the m × q matrix X = [X(V 1 ) • • • X(V q )]
of geometric features. X may contain information on the shape or location of the parcels.

Graph edges: structural information. The set of edges is represented by a binary adjacency matrix A = (a ij ) ∈ R q×q , where a ij = 1 if parcels V i and V j are spatially adjacent (i.e if ∃v i ∈ V i , ∃v j ∈ V j so that v i and v j are neighbors, for a given neighborhood definition), and a ij = 0 otherwise. This adjacency matrix encodes the spatial structure of the activation pattern.

Full graphical model. Using these definitions, we have defined an attributed region adjacency graph [START_REF] Pavlidis | Structural pattern recognition[END_REF] (hereafter noted a-RAG) G = (V, A, F, X ), which represents the fMRI activation pattern within the ROI R by encoding functional, geometric and structural information.

Kernel design

The fully generic family of convolution kernels [START_REF] Haussler | Convolution kernels on discrete structures[END_REF] is defined as:

K(G, H) = g⊆G,h⊆H t k t (g, h), (1)
where t ∈ N * is the, usually small, number of base kernels k t , which act on subgraphs g and h. To design a kernel for our a-RAGs, we need to choose the type of subgraphs, the value of t and to instantiate each base kernel k t . Subgraphs used to design such kernels include paths and random walks [START_REF] Gärtner | Exponential and geometric kernels for graphs[END_REF], tree motives [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF] etc. For simplicity reasons, we will here use paths of length two; note that the definitions below are directly extendable to other choices. Given the characteristics of the a-RAGs, we need to define t = 3 elementary kernels k f , k g and k s , respectively acting on functional, geometric and structural features.

Let G = (V G , A G , F G , X G ) and H = (V H , A H , F H , X H
) be two a-RAGs: we note g ij = {i, j} and h kl = {k, l} two pairs of nodes in G and H, respectively; let q G and q H be the number of nodes in G and H, respectively -note that q G and q H may be different.

Functional kernel. Kernel k f aims at measuring the similarity of the activation in parcels of g ij and h kl . We therefore propose the following product (and thus positive definite) kernel

k f (g ij , h kl ) = e -F G i -F H k 2 /2σ 2 f • e -F G j -F H l 2 /2σ 2 f , (2)
where

σ f ∈ R * + and F G p (resp. F H p ) is the pth column of F G (resp. F H ).
Using such a product of Gaussian kernels allows dealing with the inter-subject functional variability.

Geometric kernel. The second base kernel k g acts on the geometric features. To allow for inter-subject variability, we follow the same principle as for the functional kernel, which gives:

k g (g ij , h kl ) = e -X G i -X H k 2 /2σ 2 g • e -X G j -X H l 2 /2σ 2 g , (3)
where σ g ∈ R * + , and X G p (resp. X H p ) is the pth column of X G (resp. X H ). Structural kernel. The base kernel k s aims at valuing the structural similarity of G and H. Since our main hypothesis is that the structure is consistent across patterns recorded in different subjects, we adopt a decision function (by opposition to the smooth functional and geometric kernels), by using the linear kernel on binary entries a G ij and a H kl of the adjacency matrices A G and A H , which encodes the fact that pairs g ij and h kl are both edges:

k s (g ij , h kl ) = a G ij .a H kl ( 4
)
Resulting kernel. With the definitions of k f , k g and k s , we may define the resulting kernel (with parameters σ g and σ f ): 4 Real experiment and data analysis fMRI tonotopy experiment. In order to test this framework, we used a dataset that was acquired to study the tonotopic property of the human auditory cortex. This property states that neighboring neurons in the cortex respond to auditory stimuli of neighboring frequencies. This results in a spatially organized mapping of the auditory cortex. Typically higher frequencies are represented in the medial part of the primary auditory cortex (A1) while lower frequencies are represented more laterally, although recent high resolution studies suggest a more complex patterns with mirror-symmetric frequency gradients [START_REF] Humphries | Tonotopic organization of human auditory cortex[END_REF].

K(G, H) = q G i,j=1 q H k,l=1 k g (g ij , h kl ) • k f (g ij , h kl ) • k s (g ij , h kl ), (5)
For each of the nine subjects, a T1 image was acquired (1mm isotropic voxels). Each stimulus consisted of a 8s sequence of 60 isochronous tones covering a narrow bandwidth around a central frequency f . There were five types of sequences (i.e. conditions), each one centered around a different frequency f (f ∈ {300Hz, 500Hz, 1100Hz, 2200Hz, 4000Hz}), with no overlap between the bandwidths covered by any two types of stimuli. Five functional sessions were acquired, each containing six sequences per condition presented in a pseudorandom order. Echo-planar images (EPI) were acquired with slices parallel to the sylvian fissure (repetition time=2.4s, voxel size=2x2x3mm, matrix size 128x128). Fieldmaps were recorded to allow EPI geometric distortion correction.

Data processing and classification. The preprocessing of the functional data, carried on in SPM8 (www.fil.ion.ucl.ac.uk/spm), consisted in realignment, slice timing correction and "fieldmap" unwarping. Then, a generalized linear model was performed (in nipy: nipy.sourceforge.net) with one regressor of interest per stimulus. The weight (beta) maps of these regressors served as estimates of the response size for each stimulus. For the spatial proximity to have an anatomical meaning, one has to work on the 2D cortical surface. We therefore used freesurfer (surfer.nmr.mgh.harvard.edu) to extract the cortical surface from the T1 image and automatically define the primary auditory cortex (Fig. 1A) as part of Heschl's gyrus (thus defining our cortical ROI R in each subject). The beta maps are then projected onto the cortex, and slightly smoothed (equivalent fwhm of 3mm) in cortical space, which defines the function {f (v)} v∈R (Fig. 1B). A 2D local coordinates system is defined through a conformal mapping of R onto a rectangle [START_REF] Auzias | Model-driven harmonic parameterization of the cortical surface[END_REF], defining {x(v)} v∈R (Fig. 1C). One then need a parcellation technique that produces homogeneous parcels. We use Ward's hierarchical clustering algorithm on anatomo-functional features {f (v), x(v)}, with an added spatial constraint: the merging criterion of two adjacent parcels consists in minimzing the variance across all parcels [START_REF] Michel | A supervised clustering approach for fMRI-based inference of brain states[END_REF] (Fig. 1D), and the neighboring criterion is given by the neighboring of the vertices of the cortical mesh. Finally, we define one functional feature F , the mean activation in a parcel (normalized between 0 and 1) and a 2D geometric feature vector X as the coordinates of the barycenter of each parcel.

Using the kernel trick, one can then directly perform Support Vector Classification in graph space (G-SVC) to guess the class of the input stimulus (i.e the tone frequency) from a given activation pattern. In practice, the Gram matrix was computed with custom-designed code, and given as input to the SVC function of the scikit-learn python module (scikit-learn.org), which uses Lib-SVM (www.csie.ntu.edu.tw/ ~cjlin/libsvm/) to implement SVM. A leaveone-subject-out cross-validation scheme was used to measure the classification accuracy. For each fold of this cross-validation, we estimated the kernel parameters σ f and σ g as the median euclidean distance between the functional and geometric (respectively) feature vectors of the parcels between all nodes of all example-graphs in the training set, which corresponds to a standard heuristic [START_REF] Caputo | Appearance-based object recognition using svms: which kernel should i use?[END_REF]. Two types of analyses were conducted. First, we performed simulations with a fixed number of nodes q for all graphs and all subjects; several simulations were conducted, with q ∈ {5, 6, 7, 8, 9, 10, 15, 20, 25, 30}. Second, we imposed different number of nodes for each subject, randomly taken in three intervals I q = [5, 13], [START_REF] Pavlidis | Structural pattern recognition[END_REF]22] and [START_REF] Flandin | Improved detection sensitivity of functional MRI data using a brain parcellation technique[END_REF]21]; for each of these intervals, we ran 18 simulations with different random sets of node numbers (each time with a different q for each subject).

In order to compare the performances of our framework to vector-based methods, one need to define a vertex to vertex correspondance across subjects (which is not needed with our framework). We therefore aligned the anatomy of all subjects to a standard spherical cortical space using freesurfer and projected the beta maps to this standard space. R is now common to all subjects and we can use {f (v)} v∈R as a feature set. Several classification algorithms were then used, each time with several values of their respective hyper-parameters: 1) linear SVC; 2) non linear SVC, with gaussian (with σ ∈ {10 -n } n∈{-2,-1,0,1} ) and polynomial (of order n ∈ {2, 3, 4}) kernels; 3) k-nearest neigbors (with k ∈ {3, 5, 7, 11, 15}); 4) logistic regression with l 1 and l 2 regularization (with weight λ ∈ {10 n } n∈{2,3,4} ).

Results and discussion

Table 1 contains the performances of our G-SVC framework vs. the different benchmark vector-based methods. The scores are the mean classification accuracies across all folds of the cross-validation; for the benchmark methods, the reported score is the highest one across simulations ran with all values of their hyper-parameters; for G-SVC, it is the mean performances across all values of q. All vector-based methods performed similarly, but none of them performed significantely above chance level (equal to 0.2 since there were five types of auditory stimuli). Our G-SVC framework was the only method to perform significantely above chance level (one sample t-test, p < 0.05, indicated by a ⋆), in both the left and right Heschl's gyri (HG). It also performed significantely better than the best benchmark method (linear SVC for the right HG, logistic regression for the left HG; paired t-test with matched left-out subject, p < 0.05; indicated by a ⋄). Table 2 describes the accuracy of our G-SVC depending on the number of graph nodes q. Our framework performed above chance level in all cases (one sample t-test, p < 0.05; ⋆), and significantely better than the best benchmark method (paired t-test, p < 0.05; ⋄) in all but two cases (left HG, q = 10, 15). Results of the simulations using a randomized number of nodes q for each subject are given in Table 3. In all 108 simulations, our G-SVC framework performed above chance level (one sample t-test, p < 0.05, ⋆). This shows that the proposed kernel is efficient for graphs having different number of nodes. Moreover, in 80% of the simulations (86/108), it outperformed the best benchmark method (paired t-test, p < 0.05, ⋄). Overall, the performances are very similar to the ones obtained when the number of nodes is fixed for all subjects.

Interval Iq [START_REF] Flandin | Improved detection sensitivity of functional MRI data using a brain parcellation technique[END_REF][START_REF] Norman | Beyond mind-reading: multi-voxel pattern analysis of fMRI data[END_REF] [14,22] [5,21] right HG 0.43 (18 ⋆ , 18 ⋄ ) 0.45 (18 ⋆ , 18 ⋄ ) 0.45 (18 ⋆ , 18 ⋄ ) left HG 0.35 (18 ⋆ , 12 ⋄ ) 0.36 (18 ⋆ , 9 ⋄ ) 0.36 (18 ⋆ , 11 ⋄ ) Table 3. Mean G-SVC performances with variable q over 18 simulations (number that resulted in accuracy above chance level ⋆, above the best benchmark method ⋄) Discussion. In this paper, we designed a graphical model and a graph-kernel that allows to compare and discriminate high dimensional patterns of fMRI activation across subjects. To our knowledge, this is the first time that a graph-based pattern recognition approach is used to study local fMRI activation patterns. Specifically, we have demonstrated the power of this approach in a classification task that aims at predicting experimental variables from fMRI patterns when the classifier is trained on data from a set of subjects and its generalization performance is tested on data from a different subject.

Our G-SVC framework performed above chance level in all cases, whereas none of the benchmark methods did. This validates that this framework is appropriate to deal with inter-subject variability that is not accounted for by a state of the art group alignment, such as the one offered in freesurfer. However, it is important to use the anatomical information at our disposal, as was attempted in this study by i) working on the cortical surface, and ii) defining a local coordinates system within our region of interest for each subject, thus using the most of the local individual anatomy. Indeed, we also tested G-SVC using the raw three-dimensional voxel coordinates for x, and the performances (not reported here) were systematically lower.

A major asset of using a graph kernel such as the one defined here is that the graphs can have different number of nodes, without having to strictly solve a graph matching problem, i.e to assign nodes to one another across instances. This is emphasized by the results of our simulations where the number of nodes was randomly chosen and forced to be different for each subject. In those, our framework still outperformed vector-based methods, demonstrating the ability to deal with graphs with different node numbers. We also noted a few simulations where the classification accuracy was higher than with a fixed q for all subjects. This suggests that there might exist an "optimal" number of nodes for each subject and that performances could increase by adequately choosing q for each subject. We will address this in future work, for instance by using a model selection criterion (BIC, AIC etc.) in the parcellation. Another source of potential improvement lies in selecting the hyper-parameter σ f and σ g of the kernel in a nested cross-validation, rather than using a heuristic estimator.

Conclusion

We have designed a graphical model carrying functional, geometric, and structural features to represent spatial patterns. With a custom-designed graph kernel providing a natural metric, this framework seems particularly attractive to perform classification tasks in multi-subject settings, as evidenced by the conclusive results obtained here when dealing with inter-subject functional variability in fMRI data. This modelling strategy therefore opens the possibility to characterize and discriminate populations based on their functional patterns of activation.

Acknowledgment

Thanks to the CNRS Neuro-IC program for funding, and to the Centre IRMf de Marseille and its staff for data acquisition.

Fig. 1 .

 1 Fig. 1. Graph construction process. A. Inflated cortical surface of the right hemisphere for one subject; Heschl's gyrus is highlighted. B. A normalized activation pattern C. The two coordinates fields computed within the ROI. D. Example of parcellation with 10 parcels (arbitrary colors for illustration only). E. Example of an attributed graph: the nodes are located at the barycenter of the parcels and the color encodes the average level of activation in the parcel.

Table 2 .

 2 ⋆⋄ 0.45 ⋆⋄ 0.45 ⋆⋄ 0.43 ⋆⋄ 0.47 ⋆⋄ 0.43 ⋆⋄ 0.44 ⋆⋄ 0.48 ⋆⋄ 0.47 ⋆⋄ 0.44 ⋆⋄ left HG 0.33 ⋆⋄ 0.37 ⋆⋄ 0.35 ⋆⋄ 0.38 ⋆⋄ 0.35 ⋆⋄ 0.35 ⋆ 0.34 ⋆ 0.37 ⋆⋄ 0.39 ⋆⋄ 0.38 ⋆⋄ Influence of the number of graph nodes q in G-SVC

		G-SVC		lin. SVC		non lin. SVC	k-NN	log. reg.
	right HG	0.45 ⋆⋄		0.25		0.23		0.27	0.27
	left HG	0.36 ⋆⋄		0.26		0.24		0.25	0.24
		Table 1. G-SVC vs. benchmark methods
	G-SVC q nodes	5	6	7	8	9	10	15	20	25	30
	right HG	0.48


Elodie Brunel

ESTIMATION/IMPUTATION STRATEGIES FOR MISSING DATA IN SURVIVAL ANALYSIS

Keywords:

Introduction

We consider the problem of estimation from right-censored data, when the censoring indicator is possibly missing. Let T be a random variable representing the time to death from the cause of interest. Let C denote a right-censoring random time. Under usual random censorship, the observation is Y = T ∧ C and δ = 1 I(T ≤ C) and it is assumed that T and C admit densities respectively denoted by f T and g. In addition, C is assumed to be independent of T .

When the cause of death is not recorded, the censoring indicator is missing: this is the missing censoring indicator (MCI) model, see [START_REF] Subramanian | Survival analysis for the missing censoring indicator model using kernel density estimation techniques[END_REF], which is defined as follows. Let ξ be the missingness indicator, that is ξ = 1 if δ is observed and ξ = 0 otherwise. The observed data are then given for individual i ∈ {1, . . . , n}:

(Y i , δ i , ξ i = 1) or (Y i , ξ i = 0).
We consider the so-called MAR context, which means that we assume that the indicator is Missing At Random i.e. ξ and δ are independent conditionally to Y (contrary to the MCAR case where it is assumed that the indicator are Missing Completely At Random, i.e. ξ is independent of T and C). In this paper, we mainly concentrate on different strategies for dealing with missing indicators.

This model has been considered by several authors in the last decade. For a complete overview on missing data, we refer to the monography by [START_REF] Little | Statistical analysis with missing data[END_REF]. Most papers are interested in survival function and cumulative hazard rate estimation. In particular, van der Laan and McKeague van der Laan and McKeague (1998) build a sieved nonparametric maximum likelihood estimator of the survival function in the MAR case and prove its efficiency. Their estimator is a generalization of the Kaplan-Meier estimator to this context and is the first proposal reaching the efficiency bound. [START_REF] Subramanian | Asymptotically efficient estimation of a survival function in the missing censoring indicator model[END_REF] also proposes an estimator of the survival function in the MAR case; he proves his estimate to be efficient as well.

Kernel methods have also been used to build different estimators in the MAR context. Subramanian [START_REF] Subramanian | Survival analysis for the missing censoring indicator model using kernel density estimation techniques[END_REF] estimates the cumulative hazard rate with a ratio of kernel estimators. He provides an almost sure representation, and a Central Limit Theorem (CLT). He deduces results of the same type for the survival function. A study in a similar context is also (1) : I3M, UMR 5149 CNRS, Université Montpellier 2, 34095 Montpellier cedex 5, FRANCE, tel +33 4 67 14 39 57, fax +33 4 67 14 35 58, email: ebrunel@math.univ-montp2.fr (2) : MAP5, UMR 8145 CNRS, Université Paris Descartes, FRANCE, email: fabienne.comte@parisdescartes.fr (3) : LSTA, Université Paris 6, FRANCE, email: agathe.guilloux@upmc.fr (4) : Centre de Recherche Saint-Antoine (UMR S 938) (5) :This work is supported by French Agence Nationale de la Recherche (ANR) ANR Grant "Prognostic" ANR-09-JCJC-0101-01. 1 provided by [START_REF] Wang | Asymptotically efficient product-limit estimators with censoring indicators missing at random[END_REF]. Recently, [START_REF] Wang | Probability density estimation for survival data with censoring indicators missing at random[END_REF] proposed density estimator based on kernels and Kaplan Meier-type corrections of censoring. They prove a CLT and suggest a bandwidth selection strategy. Extensions of these works to conditional functions (both cumulative hazard and survival functions) in the presence of covariates is developed in Wang and Shen [START_REF] Wang | Estimation and confidence bands of a conditional survival function with censoring indicators missing at random[END_REF].

Our aim here is to study nonparametric hazard rate estimation and more specifically, regression and model selection estimation of the function defined, for all t > 0 by:

(1)

λ(t) = f T (t) 1 -F T (t)
,

where f T and F T are respectively the probability density function (p.d.f.) and the cumulative distribution function (c.d.f.) of the survival time T . Two contrasts and estimators are proposed, and both involve the estimation of a specific nuisance functional parameter, π or ζ:

π(y) = E(ξ|Y = y) = P(ξ = 1|Y = y). (2) ζ(y) = E(δ|Y = y) = P(δ = 1|Y = y), (3)
We compare different estimation-imputation strategies for recovering π and ζ. More precisely, we propose either a parametric strategy following a logistic modelstandard for π, but involving imputation steps for ζ -, or a pure nonparametric regression strategy. We provide theoretical properties and numerical comparisons for these procedures.

Model and strategies

2.1. Model assumptions. The unknown function λ to be estimated is hazard rate of the random variable T , as defined by (1). We shall denote by G C the c.d.f. of C. We define the conditional expectations of ξ and δ by (2) and (3).

The global model assumption is denoted (A0) and has several parts, specified hereafter.

(A0-1) The random vectors (Y i , C i ) are independent copies, for i = 1, . . . , n, of (Y, C). (A0-2) For i = 1, . . . , n, we observe

Y i = T i ∧ C i , ξ i ∈ {0, 1}, and δ i = 1 I(T i ≤ C i ) if ξ i = 1, otherwise ξ i = 0. (A0-3) C is independent of T .
(A0-4) ξ and δ are independent given Y .

Our estimation strategies are based on penalized contrasts. Namely, we obtain different collections of estimator by minimization of empirical criteria on functional spaces spanned by orthonormal bases defined on an interval [0, τ ]. We consider in the following two specific and classical examples of bases and model collections:

(1) Trigonometric bases. They are defined by ϕ 0

(x) = (1/ √ τ )1 I [0,τ ] (x), ϕ 2j+1 (x) = 2/τ sin(2πjx/τ )1 I [0,τ ] (x), ϕ 2j (x) = 2/τ cos(2πjx/τ )1 I [0,τ ] (x).
Considering (ϕ j ) 0≤j≤m-1 yields spaces S m spanned by them, with odd dimensions m. We denote by S n the nesting space of the collection, i.e. the space corresponding to the maximal dimension for S m and we set D n := dim(S n ) ≤ n/ log(n). (2) Histogram bases. They are defined by ϕ j (x) = 2 m /τ 1 I [(j-1)τ /2 m ,jτ /2 m [ (x), for j = 1, . . . , 2 m so that S m = span(ϕ j , j = 1, . . . , 2 m ) and D m = 2 m . We shall take m ≤ [log 2 (n)/2] where [z] denotes the integer part of z and log 2 (x) = log(x)/ log(2). We denote by S n the nesting space of the collection and we set D n := dim(S n ) ≤ n/ log(n).

These bases are representative examples of localized bases for the second one (as piecewise polynomials, wavelets) or bounded non localized bases for the first one.

2.2. Strategy involving knowledge of ζ. The crucial property for the construction of an estimation procedure is the following: for any integrable function h, we have

E(ζ(Y )h(Y )) = E[E(δ|Y )h(Y )] = E(δh(Y )) = E[1 I(T ≤ C)h(T )] = h(t)(1 -G C )(t)f T (t)dt.
This yields the equality

(4) E(ζ(Y )h(Y )) = E(δh(Y )) = h(y)λ(y)(1 -L(y))dy with 1 -L(y) = L(y) := P(Y ≥ y) = (1 -F T (y))(1 -G C (y)).
Next, we note that, with assumption (A0-4) and the definition of ζ, we have

E(δ i ξ i + (1 -ξ i )ζ(Y i )|Y i ) = E(δ i |Y i )E(ξ i |Y i ) + E[(1 -ξ i )E(δ i |Y i )|Y i ] = E(E(δ i |Y i )(ξ i + (1 -ξ i ))|Y i ) = E(δ i |Y i ). (5)
If ζ was known, we would consider the contrast:

Γ th n,1 (h) = 1 n n i=1 1 0 h 2 (y)1 I(Y i ≥ y)dy - 2 n n i=1 (ξ i δ i + (1 -ξ i )ζ(Y i )) h(Y i ),
which is a particular case of the contrast introduced in Brunel et al. [START_REF] Brunel | Nonparametric estimation for survival data with censoring indicators missing at random[END_REF]. Indeed, if we compute the expectation of this theoretical contrast, we obtain, under the MAR assumption and using (4) and ( 5),

E Γ th n,1 (h) = h 2 µ 1 -2 h(y)λ(y)dµ 1 (y) = h -λ 2 µ 1 -λ 2 µ 1 ,
with dµ 1 (y) = (1 -L(y))dy. Clearly, the above quantity is small if h is near of λ, and the measure denoted by µ 1 plays the role of a reference weighting norm. This explains why minimizing Γ th n,1 over an appropriate set of functions would be a relevant strategy to estimate λ.

As ζ is unknown, we must substitute another quantity in the contrast. The first strategy we could think about, is to use an estimator ζ of ζ. Consequently, we consider

(6) Γ est n,1 (h) = 1 n n i=1 1 0 h 2 (y)1 I(Y i ≥ y)dy - 2 n n i=1 ξ i δ i + (1 -ξ i ) ζ(Y i ) h(Y i ).
Estimators of ζ(y) are constructed below. This strategy of estimation of the unknown hazard rate λ, via an estimation of ζ, is also considered in Wang et al. [START_REF] Wang | Probability density estimation for survival data with censoring indicators missing at random[END_REF]. Such a regression strategy was also considered in [START_REF] Brunel | Penalized contrast estimation of density and hazard rate with censored data[END_REF] without missing indicator.

Strategy involving knowledge of π.

A similar reasoning can be lead in case π is known. Analogously, we consider the contrast

Γ th n,2 (h) = 1 n n i=1 1 0 h 2 (y)π(y)1 I(Y i ≥ y)dy - 2 n n i=1 δ i ξ i h(Y i ).
Indeed, we can compute the expectation of Γ th n,2 (h). First, we have

E(δ i ξ i h(Y i )) = E(E(δ i ξ i |Y i )h(Y i ))
and with (A0-4), we get

E(δ i ξ i h(Y i )) = E(ζ(Y i )π(Y i )h(Y i )). Now, applying (4) to hπ (instead of h) yields E(δ i ξ i h(Y i )) = h(y)π(y)λ(y)(1 -L(y))dy.
Consequently,

E Γ th n,2 (h) = h 2 (y)dµ 2 (y) -2 h(y)λ(y)dµ 2 (y) = h -λ 2 µ 2 -λ 2 µ 2 , with dµ 2 (y) = π(y)(1 -L Y (y))dy.
As previously, this explains why minimizing Γ th n,1 should allow us to estimate λ. As π is unknown, we replace it by an estimator π in the contrast. Thus, we consider

(7) Γ est n,2 (h) = 1 n n i=1 1 0 h 2 (y)π(y)1 I(Y i ≥ y)dy - 2 n n i=1 δ i ξ i h(Y i ).
2.4. Estimation of ζ or π: logit or nonparametric regression. In both ( 6) and ( 7), a step has to be completed: we must make proposals for ζ and π. Two strategies can be naturally considered. First, coherently with the nonparametric setting of the first step, we can use a nonparametric regression method.

Nonparametric regression for estimation of ζ(y) or π(y). The definition of ζ and π as conditional expectations given Y (see ( 3) and ( 2)), allows us to estimate these functions by nonparametric mean square contrast minimization. More precisely, we obtain a collection of estimators of ζ or π on an interval A by minimizing for T varying over spaces S m , the contrasts

γn,1 (T ) = 1 n n i=1 [ξ i T 2 (Y i ) -2ξ i δ i T (Y i )] or γn,2 (T ) = 1 n n i=1 [T 2 (Y i ) -2ξ i T (Y i )]
respectively.

On the other hand, a natural strategy to estimate conditional indicators is the parametric logistic regression. Logistic regression. If we add a logistic assumption

ζ(y) = exp(α 0 + α 1 y) 1 + exp(α 0 + α 1 y) , π(y) = exp(β 0 + β 1 y) 1 + exp(β 0 + β 1 y) ,
then we may use a maximum likelihood (M.L.) procedure. Note that, nevertheless, ζ and π are different from this point of view. Indeed, concerning the estimation of π, ξ and Y are observed and therefore, β 0 and β 1 correspond to standard logit-estimation. But this is not true for ζ: we cannot directly compute the M.L. estimators for (α 0 , α 1 ) since there are censoring indicators δ i missing. Thus, we follow a procedure proposed by [START_REF] Rubin | Multiple imputation for nonresponse in surveys[END_REF]. For a binary classification variable, a new logistic regression model is simulated from the posterior predictive distribution of the parameters and is used to impute the missing values for each variable .

2.5. Computing the hazard estimators. We consider that we estimate the hazard rate on a compact set A = [0, τ ]. We set standard assumptions about boundedness from above and below.

(A1.1) ∀y ∈ A, λ(y) ≤ λ A,∞ < +∞. (A1.2) ∀y ∈ A, 1 -L Y (y) ≥ L0 = min y∈A 1 -L Y (y) > a > 0 for fixed positive constant a. (A1.3) ∀y ∈ A, π(y) ≥ π 0 = inf y∈A π(y) > b > 0 for fixed positive constant b. (A1.4) ∀y ∈ A, 0 < f * 0 ≤ f Y (y) ≤ f * 1 < +∞. Assumptions (A1.1
) and (A1.2) are common natural assumptions, and assumption (A1.3) is specific to the strategy involving π in the reference measure. Assumption (A1.4) is a common technical assumption.

First, we define two types of estimators λ(i) m for i = 1, 2 on the space S m by: λ

(i) m = arg min h∈Sm Γ est n,i (h) where S m = span{ϕ j , j = 1, . . . , D m }
The ϕ j 's constitute an L 2 -orthonormal basis, and the function h is of the form h = j a j ϕ j . The estimators are obtained by minimization of the contrasts ( 6)-( 7) given in Sections 2.2 and 2.3. For the histogram basis, we can give their explicit expression in the MAR setting :

λ(i),[H] m (y) = D m j=1 â(i),[H] j ϕ j (y) with the coefficient â(i),[H] j
having the following form, provided the denominators are non zero:

â(1),[H] j = n i=1 [δ i ξ i + (1 -ξ i ) ζ(Y i )]ϕ j (Y i ) n i=1 1 0 ϕ 2 j (y)1 I (Y i ≥y) dy , â(2),[H] j = n i=1 δ i ξ i ϕ j (Y i ) n i=1 1 0 ϕ 2 j (y)π(y)1 I(Y i ≥ y)dy .
More generally, let us define the matrices and vectors:

-→ â (i) = â(i) j 1≤j≤Dm , Θ (1) m := 1 n n i=1 ϕ j (y)ϕ j ′ (y)1 I {Y i ≥y} dy 1≤j,j ′ ≤Dm , ∆ (1) m = vec 1 n n i=1 δ i ξ i + (1 -ξ i ) ζ(Y i ) ϕ j (Y i )) 1≤j≤Dm , Θ (2) m := 1 n n i=1 ϕ j (y)ϕ j ′ (y)π(y)1 I {Y i ≥y} dy 1≤j,j ′ ≤Dm ∆ (2) m = 1 n n i=1 δ i ξ i ϕ j (Y i ) 1≤j≤Dm .
Then the coefficients of the estimators must fulfill the matrix constraint:

Θ (i) m -→ â (i) m = ∆ (i) m . It follows that the estimator is well defined if Θ (i)
m is invertible. We define ρ(M ) as the spectral radius of a matrix M , i.e. the largest eigenvalue in modulus of M . We set

(8) -→ â (i) m = (Θ (1) m ) -1 ∆ (1) m if ρ(Θ (i) m ) ≥ max(ĉ i , n -1/2 ) and
-→ â (i) m = 0 otherwise, with ĉ1 = L0 /3 and ĉ2 = 4π 0 L0 /9.

The quantities L0 and π0 are estimators of L0 = min y∈A 1 -L Y (y) > 0 under (A1.2) and π 0 = min y∈A π(y) > 0 under (A1.3). An estimator L0 is defined in Comte et al. (2009), and proved to have the required properties, namely, to satisfy: (A2.1) For any integer k ≥ 1, there exists a constant C

(L 0 ) k > 0 such that P(| L0 -L0 | > L0 /2) ≤ C (L 0 ) k /n k .
For π0 , we can study min π(y) but we state results in term of the conditions π must fulfill:

(A2.2) P π -π π ∞,A > 1/4 := P(Ω c π ) ≤ C n k for k ≥ 4 and C a given constant.
After this stage, the coefficients â(i) j , for i = 1, 2 provide the development of an estimator

(9) λ(i) m (y) = j â(i) j ϕ j (y).
The model selection device is now based on the following criterion:

for i = 1, 2, (10) m(i) = arg min m∈Mn (Γ est n,i ( λm ) + pen (i) (m))
where

M n = {m ∈ N, dim(S m ) ≤ N n }. Moreover (11) pen (i) (m) = κ i λ ∞,A dim(S m ) n ,
where Φ 0 is such that ∀t ∈ S m , t ∞ ≤ Φ 0 t and λ = λm 0 is an estimator in the collection, on a space S m 0 with dimension D m 0 such that log(n) ≤ D m 0 ≤ n 1/4 . Moreover, the dimension D m of the spaces S m has to be larger than log(n).

The computation of ζ and π is similar and the vector of their coefficients in the basis (ϕ j ) 1≤j≤dm is defined by the inverse of matrices

Υ (1) m = 1 n n i=1 ξ i ϕ j (Y i )ϕ j ′ (Y i ) 1≤j,j ′ ≤Dm , Υ (2) m = 1 n n i=1 ϕ j (Y i )ϕ j ′ (Y i ) 1≤j,j ′ ≤Dm
, multiplied respectively by vectors

Ξ (1) m = ∆ (2) m , Ξ (2) m = 1 n n i=1 ξ i ϕ j (Y i ) 1≤j≤Dm
with adequate thresholds.

2.6. Theoretical Results. We consider that we estimate the hazard rate on a compact set

A = [0, τ ].
The empirical reference norms are defined by

h 2 n,1 = 1 n n i=1 1 0 h 2 (y)1 I(Y i ≥ y)dy, h 2 n,2 = 1 n n i=1 1 0 h 2 (y)π(y)1 I(Y i ≥ y)dy
associated with the natural scalar product h 1 , h 2 n,i obtained by polarization identity.

Let us also define another integral norm with respect to d̺(y) = f Y (y)dy where f Y is the pdf of Y 1 , that is

(12) ψ 2 ̺ = ψ 2 (y)d̺(y) = ψ 2 (y)f Y (y)dy.
The result we can obtain is the following, for estimators ζ or π computed on an independent sample of observations.

Theorem 1. Let λ(i) m(i) be the estimator defined by ( 8)-( 9)-( 10)-( 11) for i = 1, 2, under Assumptions (A0). Assume that the collection of models is nested and the estimators ζ and π are computed on independent samples.

(1) If (A1.1)-(A1.2)-(A1.4) and (A2.1) hold, there exists a choice of κ 1 such that, for n large enough

(13) E( λ1 I A - λ(1) m(1) 2 n,1 ) ≤ C 1 inf m∈Mn ( λ1 I A -λ m 2 + λ A,∞ D m n ) + C ′ 1 E( ζ -ζ 2 ̺ ) + C ′′ 1 n ,
where C 1 is a numerical constant and C ′ 2 , C ′′ 2 are constants depending on the constants of the problem.

(2) If (A1.1)-(A1.4) and (A2.1)-(A2.2) hold, there exists a choice of κ 2 such that, for n large enough

(14) E( λ1 I A - λ(2) m(2) 2 n,2 ) ≤ C 2 inf m∈Mn ( λ1 I A -λ m 2 + λ A,∞ D m n ) + C ′ 2 E( π -π 2 ̺ ) + C ′′ 2 n ,
where C 2 is a numerical constant and C ′ 2 , C ′′ 2 are constants depending on the constants of the problem.

Inequality (13) is proved in [START_REF] Brunel | Nonparametric estimation for survival data with censoring indicators missing at random[END_REF] in a more general setting, and a sketch of the proof of ( 14) is given in appendix.

Note that the larger κ i , the larger pen (i) (m) in the upper bound.

The results stated in ( 13)-( 14) involve three terms:

• the first terms are common, inf m∈Mn ( λ1 I A -λ m 2 + λ A,∞ D m /n and correspond to the squared-bias ( λ1 On the other hand, with the logit strategy, we have ζ(y) = e α0 +y α1 1 + e α0 +y α1 , π(y) = e β0 +y β1

I A -λ m 2 )/variance ( λ A,∞ D m /n) compromise, • the second term is E( ζ -ζ 2 ̺ ) or E( π -π 2 ̺ ),
1 + e β0 +y β1 , and we get

E( ζ -ζ 2 ̺ ) ≤ 2 E((α 0 -α 0 ) 2 ) + E(Y 2 1 )E((α 1 -α 1 ) 2 ) , E( π -π 2 ̺ ) ≤ 2 E(( β0 -β 0 ) 2 ) + E(Y 2
1 )E(( β1 -β 1 ) 2 ) . All these terms are of order 1/n as soon as all parameters are estimated with parametric rate.

The questions that we have to study empirically are: For K = 500 replications over different paths, we compute the (empirical) average MISE of the penalized estimators λ over a grid of size 100:

with θ =
MISE = 1 K K k=1 τ k 100 100 i=1 λ(t i ) -λ(k) (t i ) 2 ,
where τ k is the inter-quantile interval length associated with the 10% and 90% empirical quantiles of the Y i 's . The value of the constant κ appearing in the penalty has been calibrated over the three models and fixed to 2 (resp. 3.5) for the contrast Γ est n,1 (resp. Γ est n,2 ). We also give in parenthesis the standard deviation value of the MISE evaluated over the 500 samples. The results are summarized in Tab. 1-4. For each model, we give the MISEs of our penalized estimators λ(1) m(1) for the estimation strategy involving ζ in column 1-2, λ(2) m(2) for the estimation strategy involving π column 3-4 and λ(3) m(3) for the imputation strategy in column 5-6. For each strategy and each model, we consider various sample sizes n = 200, 1000 and 5000 and two levels of missing indicators namely 15% and 45%. The abbeviations "NP" and "logit" stand for the way of estimating the auxiliary functions ζ or π according to the strategy under study (estimation of ζ with Γ est n,1 and Γ imp n or estimation of π with Γ est n,2 ). First, whatever the model we considered, we observe that the strategy involving π gives very poor performances compared to the other ones. The difference between the strategies involving π or ζ is of great importance for a numerical point of view. In fact, as we can see by examining the sequence of coefficients â(1),[H] j and â(2),[H] j , the estimator π appear in the denominator (while the estimator ζ does not) and it results that a small estimation error for π can strongly increase the final error on the hazard estimator. We can see on Figure 1 that despite a quite good estimation of π, the MISE in Table 1-3 is always the worst for the estimation strategy involving π. Consequently, we focus on the other strategies. Both imputation or estimation strategy involving the function ζ have very similar behavour and perform as well. We also notice that the logit estimator of ζ versus the nonparametric one is nearly impossible to beat at least for sample sizes n = 200, 1000. Of course, when the logit assumption is violated, for large sample (n = 5000) and important missing rate 45%, the results are deteriorating, see 3.

(h) = ν n,1 (h) -ν n,2 (h) with ν n,1 (h) = 1 n n i=1 (δ i ξ i h(Y i ) -h, λ µ 2 ), (16) ν n,2 (h) = 1 n n i=1 h(y)λ(y)π(y)1 I {Y i ≥y} dy -h, λ µ 2 . ( 17
)
We define a set Ω = ∆ ∩ ΩL 0 ∩ Ω π where Ω π is defined in (A2.2) and Ĝ

= m∈Mn Ĝm , Ĝm = {min Sp(Θ (2) m ) ≥ max(4π 0 L0 /9, n -1/2 )}, ∆ = ∀h ∈ S n , h 2 n,2 h µ 2 -1 ≤ 1 2 and ΩL 0 = | L0 L0 -1| ≤ 1 2 .
It is easy to see that (see [START_REF] Brunel | Nonparametric estimation for survival data with censoring indicators missing at random[END_REF]):

Lemma 1. Under the assumptions of Theorem 1, for n large enough, the following embedding holds:

∆ ∩ ΩL 0 ∩ Ω π ⊂ Ĝ ∩ ΩL 0 ∩ Ω π .
As a consequence, for all m ∈ M n , the matrices Θ

(2) m are invertible on Ω.

Let us denote by λ A = λ1 I A and π A = π1 I A . We start with the decomposition

Γ est n,2 (h) -Γ est n,2 (ℓ) = h -λ A 2 n,2 -ℓ -λ A 2 n,2 -2ν n (h -ℓ) + 1 n n i=1 (h -λ A ) 2 (y)1 I {Y i ≥y} (π -π A )(y)dy - 1 n n i=1 (ℓ -λ A ) 2 (y)1 I {Y i ≥y} (π -π A )(y)dy + 2 n n i=1 (h -ℓ)(y)λ(y)1 I {Y i ≥y} (π -π A )(y)dy.
Then, let us denote by λ m the orthogonal projection of λ A on S m and use the definition of λ m on Ω (see Lemma 1). We assume first that λ A,∞ is known and denote by pen(m) = κ λ A,∞ D m /n, and write Γ est n,2 ( λ m) + pen( m) ≤ Γ est n,2 (λ m ) + pen(m). Then we get, on Ω,

λ m -λ A 2 n,2 ≤ λ m -λ A 2 n,2 + pen(m) + 2ν n ( λ m -λ m ) -pen( m) - 1 n n i=1 ( λ m -λ A ) 2 (y)1 I {Y i ≥y} (π -π A )(y)dy + 1 n n i=1 (λ m -λ A ) 2 (y)1 I {Y i ≥y} (π -π A )(y)dy - 2 n n i=1 ( λ m -λ m )(y)λ(y)1 {Y i ≥y} (π -π A )(y)dy. (18) First, we write 2|ν n ( λ m -λ m )| ≤ 1 8 λ m -λ m 2 µ 2 + 8 sup h∈B m, m(0,1) ν 2 n (h), where B m,m ′ (0, 1) = {h ∈ S m + S m ′ , h µ 2 ≤ 1}. It follows that 2|ν n ( λ m -λ m )|1 I Ω ≤ 1 4 λ m -λ m 2 n,2 1 I Ω + 16 sup h∈B m, m(0,1) ν 2 n,1 (h) -p(m, m) + +16p(m, m) + 16 sup h∈B m, m(0,1) ν 2 n,2 (h) ≤ 1 2 λ m -λ 2 n,2 1 I Ω + 1 2 λ m -λ 2 n,2 + 16 sup h∈B m, m(0,1) ν 2 n (h) -p(m, m) + +16p(m, m) + 16 sup h∈B m, m(0,1) ν 2 n,2 (h), (19)
where p(m, m ′ ) is defined below.

Next, using the definition of Ω π and assumption (A2.2), we have

1 n n i=1 (λ m -λ A ) 2 (y)1 I {Y i ≥y} (π -π A )(y)dy 1 I Ωπ ≤ 1 4n n i=1 (λ m -λ A ) 2 (y)1 I {Y i ≥y} π A (y)dy ≤ 1 4 λ m -λ A 2 n,2 . (20)
In the same way,

(21) 1 n n i=1 ( λ m -λ) 2 (y)1 I {Y i ≥y} (π -π)(y)dy 1 I Ωπ ≤ 1 4 λ m -λ 2 n,2 .
Moreover, using first the Schwarz Inequality and then 2xy ≤ x 2 /a + ay 2 for x, y and a > 0 (here a = 16), we get We can use Talagrand Inequality (see [START_REF] Talagrand | New concentration inequalities in product spaces[END_REF]) to prove the proposition below, and we refer to [START_REF] Brunel | Nonparametric estimation for survival data with censoring indicators missing at random[END_REF] for details. For the method allowing to replace λ A,∞ by an estimator, we refer to Comte et al. (2008). 2 6. Bibliography

2 n n i=1 ( λ m -λ m )(y)λ(y)1 I {Y i ≥y} (π -π A )(y)dy ≤ 1 16n n i=1 ( λ m -λ m ) 2 (y)1 I {Y i ≥y} π(y)dy + 16 n n i=1 [(π -π A ) 2 (y)/π(y)]λ 2 (y)1 I A (y)dy ≤ 1 16 λ m -λ m 2 n,2 + 16 λ ∞,A nπ 0 n i=1 (π -π A ) 2 (

Figure 1 .

 1 Figure 1. Strategy involving estimation of π: Left-side: True hazard rate (black line) for Model 2 and its estimators λ(2) m(2) with both nonparametric (magenta dotted line) and logit (cyan plain line) estimators of π. Right-Side: True auxiliary function π (black line) and its nonparametric (magenta) and logit (cyan) estimators. n = 1000 with ≃ 30% censoring and 45% missing indicators.

Figure 2 .

 2 Figure 2. Strategies involving estimation of ζ: Left-side: True hazard rate (black line) for Model 3 and its estimators λ(1) m(1) with nonparametric (magenta dotted line) and logit (cyan plain line) estimators of ζ and λ(3) m(3) with nonparametric estimator of ζ (red dotted line). Right-Side: True auxiliary function ζ (black line) and its nonparametric (magenta) and logit (cyan) estimators. n = 1000 with ≃ 30% censoring and 45% missing indicators.

  y)λ(y)1 I A (y)dy Therefore, as π is computed on an independent sample of observations, -λ m )(y)λ(y)1I {Y i ≥y} (π -π)(y)dy 1 I Ω h∈B m, m(0,1) ν 2 n,1 (h) -p(m, m) + + 16E(p(m, m)pen( m)). (23) Now write ν n,2 (h) = A h(y)λ(y)( L(y) -E( L(y))π(y)dy with L(y) = (1/n) n i=1 1 I Y i ≥y . y)( L(y) -E( L(y)) 2 dµ 2 (y).Therefore, since E[( L(y) -E( L(y)) 2 ] = L(y)(1 -L(y))/n ≤ 1/(4n) and π(y) ≤ 1

Proposition 1 .

 1 Under the Assumptions of Theorem 1, there exists a numerical constant κ such forp(m, m ′ ) = (κ/16) λ 2 A,∞ D m + D m ′ n , h) -p(m, m)) + ≤ C nand this, together with (23) and (24) yields, as 16p(m,m ′ ) ≤ pen(m) + pen(m ′ )study E( λ m -λ 2 n,2 1 I Ω c ).We have first (see Brunel et al. Brunel et al. (2013)): Lemma 2. Under the assumptions of Theorem 1, we have λ m 2 n,2 ≤ Cn 3 for a constant C depending on the basis. Moreover (A2.2) ensures that P(Ω c π ) ≤ C/n k , (A2.1) implies that P(Ω c L0 ) ≤ C (L 0 ) k /n k and Proposition 4 in Comte et al. (2008) can be adapted here to get P(∆ c ) ≤ C (∆) k /n k , under the condition that D n 2 ≤ n/ log 2 (n) if basis (1) is used and D n ≤ n/ log 2 (n) if basis (2) is used. Gathering these elements yields E( λ m -λ 2 n 1 I Ω c ) ≤ C/n by choosing k = 4.

  that is the mean-square risk of the estimator of ζ or π on A. These quantities can have different orders, depending on the second step strategy, • the last terms C ′′ 1 /n or C ′′ 2 /n are negligible. We can obtain adaptive nonparametric rate with the same model selection principle for ζ and π, if we apply the nonparametric strategy (see Theorem 2 for ζ in Brunel et al. (2013) Brunel et al. (2013)).

  1 (resp. 0.32) to obtain 85 % (resp. 55 %) non missing indicators, and α 1 = -2, α 2 = 5.2. Note that, the hazard rate h(y) coincides with the auxiliary function ζ(y) in this model.

	Model 3 (Non-monotone). Consider the same schema as Model 2, with
		ζ(y) = 0.9| sin(π/2 -1.5y)| + 0.1.
	Weibull/	Estimation strategy	Estimation strategy	Imputation strategy
	Logistic	with Γ est n,1		with Γ est n,2		Γ imp n
	Missing rate 15%	45 %	15%	45 %	15 %	45 %
	n = 200
	NP	0.055	0.061	0.120	0.269	0.055	0.067
		(0.032)	(0.043)	(0.081)	(0.162)	(0.032)	(0.051)
	logit	0.055	0.053	0.151	0.286	0.053	0.055
		(0.029)	(0.032)	(0.123)	(0.160)	(0.031)	(0.032)
	n = 1000
	NP	0.017	0.021	0.035	0.086	0.018	0.024
		(0.010)	(0.014)	(0.018)	(0.037)	(0.010)	(0.015)
	logit	0.016	0.016	0.062	0.092	0.016	0.019
		(0.008)	(0.009)	(0.036)	(0.026)	(0.008)	(0.012)
	n = 5000
	NP	0.008	0.009	0.011	0.025	0.008	0.009
		(0.007)	(0.008)	(0.010)	(0.023)	(0.006)	(0.009)
	logit	0.007	0.008	0.023	0.024	0.007	0.008
		(0.005)	(0.006)	(0.022)	(0.023)	(0.006)	(0.006)
	Table 1. Average and standard deviation (in parenthesis) of the MISE over 500
	replicated samples for hazard rate estimators of Model 1 (Weibull-Logistic), for
	estimation strategies with Γ est n,1 and Γ est n,2 and imputation strategy Γ imp n . Censoring
	rate ≃ 30%.

Table 2 .

 2 Average and standard deviation (in parenthesis) of the MISE over 500 replicated samples for hazard rate estimators of the Model 2, for estimation strategies with Γ est n,1 and Γ est n,2 and imputation strategy Γ imp

n . Censoring rate ≃ 30%.

Table 3 .

 3 Average

	and standard deviation (in parenthesis) of the MISE over 500
	replicated samples for hazard rate estimators of Model 3, for estimation strategies
	with Γ est n,1 and Γ est n,2 and imputation strategy Γ imp n . Censoring rate ≃ 30%.
				5. Proofs
	5.1. Proof of Theorem 1. The first inequality is obtained as a particular case of the bound
	proved in Brunel et al. Brunel et al. (2013), see also Comte et al. ?. Below, we give steps for
	obtaining the second inequality.
	The reference centered empirical process is
	ν n (h) =	1 n	n i=1	δ

i ξ i h(Y i ) -h(y)λ(y)π(y)1 I {Y i ≥y} dy .

(1) If the logit assumptions are fulfilled, is the nonparametric strategy much worse than the parametric one? (2) If the logit model is not satisfied, can the logit approximation still correctly estimate π or ζ? (3) How much does the quality of the estimator of these quantities influences that of λ?

Imputation based strategy

The second approach, we aim to explore, is inspired from stochastic regression imputation. The idea is no longer to substitute an estimator of the unknown quantity ζ but rather to substitute a random number γ i to ζ(Y i ) which has the same expectation.

Let γ i generating from a Bernoulli distribution with parameter ζ(Y i ) given ξ i = 0. Here, the quantity ζ(Y i ) has also to be estimated as in Section 2.4, and we consider again both nonparametric and logit strategies. Then, we are in position to define the following contrast:

and the imputation strategy gives the following estimator:

We do not provide the theoretical study of the estimator defined by this contrast since

We propose a numerical comparison of this strategy with the other strategies developped in Section 2.

Numerical Comparison

To illustrate the behaviors of the different estimators, we consider three models to simulate data from the MCI, i.e. n independent replications of

In each model, the parameters are set in order to obtain around 30 % of censoring and 85 % (resp. 55 %) of non missing indicators.

Model 1 (Weibull-Logistic). Each T i is drawn from a Weibull distribution with scale a = 10 and shape b = 4, each C i from an exponential distribution with mean µ = 25. For each individual i, the missing indicator is simulated from a Bernoulli distribution with parameter π(y i ) with

we set β 1 = 0.1 and β 2 = 0.25 (resp. 0.005) to obtain 85 % (resp. 55 %) non missing indicators. In this model:

Model 2 (Monotone). It is based on an example from [START_REF] Subramanian | The multiple imputations based Kaplan-Meier estimator[END_REF]. For each i, Y i is drawn from an exponential distribution with parameter 1. The censoring and missing indicators are simulated from Bernoulli distributions with parameters resp. ζ(y i ) and π(y i ) with: π(y) = 1exp(-θ exp(y)) and ζ(y) = 1 1 + exp(-α 1 -α 2 y) ,