Raphaël Coudret
email: raphael.coudret@math.u-bordeaux1.fr

Benoit Liquet
email: benoit.liquet@isped.u-bordeaux2.fr

Jérôme Saracco
email: jerome.saracco@math.u-bordeaux1.fr

Comparison of sliced inverse regression approaches for underdetermined cases

Keywords: dimension reduction, high-dimensional data, semiparametric regression, sparsity Mots-clés : grande dimension, regression semi-paramétrique, réduction de dimension, sparsité AMS 2000 subject classifications: 62H12, 62F07

Among methods to analyze high-dimensional data, the sliced inverse regression (SIR) is of particular interest for non-linear relations between the dependent variable and some indices of the covariate. When the dimension of the covariate is greater than the number of observations, classical versions of SIR cannot be applied. Various upgrades were then proposed to tackle this issue such as RSIR and SR-SIR, to estimate the parameters of the underlying model and to select variables of interest. In this paper, we introduce two new estimation methods respectively based on the QZ algorithm and on the Moore-Penrose pseudo-inverse. We also describe a new selection procedure of the most relevant components of the covariate that relies on a proximity criterion between submodels and the initial one. These approaches are compared with RSIR and SR-SIR in a simulation study. Finally we applied SIR-QZ and the associated selection procedure to a genetic dataset in order to find eQTL.

Résumé : Parmi les méthodes pour analyser des données de grande dimension, la régression inverse par tranches (sliced inverse regression ou SIR en anglais) est particulièrement intéressante si des relations non-linéaires existent entre la variable à expliquer et des combinaisons linéaires des prédicteurs (appelées indices). Lorsque la dimension de ces prédicteurs est plus grande que le nombre d'observations, les versions classiques de SIR ne peuvent plus être utilisées. Des améliorations diverses comme RSIR et SR-SIR (versions régularisées de SIR) ont été proposées dans la litérature pour résoudre ce problème, estimer les paramètres du modèle sous-jacent et enfin réaliser une sélection des prédicteurs les plus pertinents (en un certain sens). Dans cet article, nous introduisons deux nouvelles procédures d'estimation basées respectivement sur l'algorithme QZ et sur l'inverse généralisé de Moore-Penrose. Nous décrivons également une méthode qui repose sur un critère de proximité entre des sous-modèles et le modèle intial pour sélectionner les prédicteurs les plus pertinents. Ces approches sont ensuite comparées avec RSIR et SR-SIR par le biais de simulations. Enfin, nous illustrons, sur un jeu de données génetiques, l'intérêt de l'approche SIR-QZ proposée et de l'algorithme de sélection de prédicteurs associé pour trouver des eQTL.

Introduction

For a univariate response variable y and a multivariate covariate x ∈ R p , the semiparametric regression model

y = f (x β 1 , . . . , x β K , ε) (1)
is an attractive dimension-reduction approach to model the effect of the p-dimensional covariates x on y. Let µ = E(x) and Σ = V(x). The error term ε is assumed to be independent of x. Since the link function f (•) is an unknown smooth function, the parameters β k ∈ R p are not entirely identifiable, only the linear subspace spanned by the β k 's can be identified without additional assumptions. Duan and Li [START_REF] Duan | Slicing regression: a link-free regression method[END_REF] and Li [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF] called this subspace the effective dimension reduction (EDR) subspace. Moreover any direction belonging to this subspace is called an EDR direction.

If the β k 's are assumed linearly independent, the EDR subspace is then a K-dimensional linear subspace of R p . Other authors refer to this subspace as the dimension reduction subspace (DRS) or the central subspace (which is defined as the smallest DRS), see Cook [START_REF] Cook | Principal hessian directions revisited (with discussion)[END_REF] for more details. When the dimension p of x is high and when we have little knowledge about the structure of the relationship between the response and the covariates, this semiparametric regression model is a nice alternative to parametric modeling (since it is really difficult to have knowledge about the structure of the relationship between the response and the covariates) and non-parametric modeling (which suffers from the well-known curse of dimensionality due to the data sparseness in the domain of x). The idea of dimension reduction in model ( 1) is intuitive because it aims at constructing a low dimensional projection of the covariate without losing information to predict the response y. If the dimension K of the EDR subspace is sufficiently small, it facilitates data visualization and explanation and it alleviates the curse of the dimensionality to non-parametrically estimate f with usual approaches such as kernel or splines smoothing (when the error term is additive).

In this semiparametric regression model [START_REF] Amato | Dimension reduction in functional regression with applications[END_REF], an important purpose is to estimate the EDR subspace from a sample {(x i , y i ), i = 1, . . . , n}. Most of the existing approaches are usually based on the eigendecomposition of a specific matrix of interest. The most popular one is the sliced inverse regression (SIR) introduced by Duan and Li [START_REF] Duan | Slicing regression: a link-free regression method[END_REF] and Li [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF] , respectively for single index models (K = 1) and multiple indices models (K ≥ 1). Among alternative methods there are SIR-II, see [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF][START_REF] Yin | Asymptotic distributions for dimension reduction in the sir-ii method[END_REF] for instance, and sliced average variance estimation (SAVE), see [START_REF] Zhu | On kernel method for sliced average variance estimation[END_REF][START_REF] Li | Asymptotics for sliced average variance estimation[END_REF] for example. These approaches require the inverse of Σ. Then, from a practical point of view, it is necessary to inverse an estimate Σ of Σ.

Define xi = (x i -μ) ∈ R p for i = 1, . . . , n, with μ = 1 n ∑ n i=1 x i . A usual (biased) estimate of Σ is

Σ = 1 n n ∑ i=1 (x i -μ)(x i -μ) = 1 n ( x1 , . . . , xn )( x1 , . . . , xn ) . (2)
Clearly, the rank of the p × p matrix Σ is at most equal to n -1 since ∑ n i=1 xi = 0 p where 0 p stands for the null vectors of R p . From this remark on the rank of Σ, this matrix is singular when n < p. Moreover, it is also often ill-conditioned when n ≈ p.

Therefore, SIR, SIR-II or SAVE methods only work well when the sample size n is greater than the dimension p of the covariate x, but naturally fail when n < p. In this underdetermined case, the standard estimate of Σ is not invertible even if the components of x are independent.

In the following, we only focus on the SIR approach. We describe it in Section 2.1 when n > p. The goal of this paper is then twofold. On one hand, we present methods to tackle the issue n < p. On the other hand, we also provide procedures in order to select which components of x have an effect on y.

In Sections 2.2-2.3, we consider two different regularizations added to the SIR method, proposed by Zhong et al [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] and Li and Yin [START_REF] Li | Sliced inverse regression with regularizations[END_REF], to find EDR estimates when n < p. Moreover, the SIR method can be seen as a generalized eigenvalue problem and linear algebra algorithms exist to solve this kind of problem without requiring any matrix inversion. The QZ algorithm (see [START_REF] Moler | An algorithm for generalized matrix eigenvalue problems[END_REF] for instance) is one of them and will be used in the SIR context in Section 2. [START_REF] Bai | A chi-square test for dimensionality for non-gaussian data[END_REF]. In Section 2.5, we also adapt an approach introduced in functional sliced inverse regression (i.e., when x is an explanatory functional variable), based on the Moore-Penrose pseudo-inverse.

Concerning the selection of useful predictors in the indices, Zhong et al [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] use a chi-square test to find which components of x affect y (Section 3.1), while the approach of Li and Yin [START_REF] Li | Sliced inverse regression with regularizations[END_REF] relies on a Lasso penalization (Section 3.2). In Section 3.3, we propose another procedure. We choose randomly some submodels (i.e., using a number p 0 < p of components of x) and we measure how close they are from the initial one with all the p components of x. The latter model is thus taken as a benchmark. Components of x that appear the most in submodels that are the closest to the benchmark are kept. We naturally consider that the other components of x do not affect y.

In Section 4, we compare in a simulation study the numerical behavior of the described methods to estimate EDR directions. We also evaluate the different procedures of selection of the useful components of x. In Section 5, we apply the most efficient one on real data from a genetic framework. Finally, some concluding remarks are given in Section 6.

SIR in determined and underdetermined cases

Brief review of usual SIR

Let β be a p × K matrix defined by β = (β 1 , . . . , β K ). The EDR subspace is thus spanned by β .

Inverse regression step.

The basic principle of the SIR method is to reverse the role of y and x, that is, instead of regressing the univariate variable y on the multivariate variable x, the covariable x is regressed on the response variable y. The price we have to pay to succeed in inverting the role of x and y is an additional assumption on the distribution of x, named the linearity condition (described hereafter).

Usual SIR estimate is based on the first moment E(x|y). It has been initially introduced by Duan and Li [START_REF] Duan | Slicing regression: a link-free regression method[END_REF] for single index model and by Li [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF] for multiple indices model. SIR approaches have been extensively studied, see for instance [START_REF] Carroll | Measurement error regression with unknown link: dimension reduction and data visualization[END_REF][START_REF] Chen | Can SIR be as popular as multiple linear regression?[END_REF][START_REF] Zhu | On hybrid methods of inverse regression-based algorithms[END_REF][START_REF] Bercu | A new approach of recursive and non recursive sir methods[END_REF][START_REF] Azais | Optimal quantization applied to sliced inverse regression[END_REF] among others.

Let us now recall the geometric property on which SIR is based. Let us introduce the linearity condition:

(LC) : ∀b ∈ R p , E(x b|x β 1 , . . . , x β K ) is linear in x β 1 , . . . , x β K . (3)
Note that this condition is satisfied when x is elliptically distributed (for instance normally distributed). The reader can find an interesting discussion on this linearity condition in [START_REF] Chen | Can SIR be as popular as multiple linear regression?[END_REF]. Assuming model ( 1) and (LC), Li [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF] showed that the centered inverse regression curve is contained in the linear subspace spanned by the K vectors Σβ 1 , . . . , Σβ K . Let T denote a monotonic transformation of y. He considered the eigendecomposition of the Σ-symmetric matrix Σ -1 M where M = V(E(x|T (y))). Straightforwardly the eigenvectors associated with the largest K eigenvalues of Σ -1 M are some EDR directions. Slicing step. To easily estimate the matrix M, Li [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF] proposed a transformation T , called a slicing, which categorizes the response y into a new response with H > K levels (in order to avoid an artificial reduction of dimension). The support of y is partitioned into H non-overlapping slices s 1 , . . . , s h , . . . , s H . With such transformation T , the matrix of interest M can be now written as

M = ∑ H h=1 p h (m h -µ)(m h -µ)
where p h = P(y ∈ s h ) and m h = E(x|y ∈ s h ).

Estimation process. When a sample {(x i , y i ), i = 1, . . . , n} is available, matrices Σ and M are estimated by substituting empirical versions of the moments for their theoretical counterparts. Let

M = H ∑ h=1 ph ( mh -μ)( mh -μ) , (4)
where ph

= 1 n ∑ n i=1 I[y i ∈ s h ] and mh = 1 n ph ∑ n i=1 x i I[y i ∈ s h ].
Therefore the estimated EDR directions are the eigenvectors associated with the K largest eigenvalues of Σ -1 M. They span the K-dimensional estimated EDR subspace The convergence at rate √ n and the asymptotic normality of estimated EDR directions have been obtained, see [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF][START_REF] Saracco | An asymptotic theory for Sliced Inverse Regression[END_REF] for instance.

The choice of the slicing T is discussed in [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF][START_REF] Kötter | Smoothing and regression. Approaches, computation and application[END_REF][START_REF] Saracco | Pooled slicing methods versus slicing methods[END_REF] but, theoretically, there is no optimal one. In practice, we fix the number of observations per slice to n/H where a stands for the integer part of a. If the sample size n is not proportional to the number H of slices, some slices will then contain [n/H] + 1 observations. Note that, in order to avoid the choice of a slicing, alternative SIR methods have been investigated. For instance, one can mention kernel-based methods of SIR proposed by Zhu and Fang [START_REF] Zhu | Asymptotics for kernel estimate of sliced inverse regression[END_REF] or Aragon and Saracco [START_REF] Aragon | Sliced Inverse Regression (SIR): an appraisal of small sample alternatives to slicing[END_REF]. However, these methods are hard to implement and are computationally slow. Moreover, Bura [START_REF] Bura | L 1 -Statistical Procedures and Related Topics, chapter Dimension reduction via parametric inverse regression[END_REF] and Bura and Cook [START_REF] Bura | Estimating the structural dimension of regressions via parametric inverse regression[END_REF] proposed a parametric version of SIR.

Concerning the determination of the dimension K of the EDR subspace (which is unknown in practice), several works are available in the literature, see for example [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF][START_REF] Schott | Determining the dimensionality in sliced inverse regression[END_REF][START_REF] Ferré | Determining the dimension in sliced inverse regression and related methods[END_REF][START_REF] Bai | A chi-square test for dimensionality for non-gaussian data[END_REF][START_REF] Liquet | Application of the bootstrap approach to the choice of dimension and the α parameter in the SIR α method[END_REF] among others.

Standardized version for SIR. Another way to obtain a basis of the EDR subspace is to consider the eigendecomposition of Σ -1/2 MΣ -1/2 , that is the eigendecomposition of M * = V(E(z|T (y)) where z = Σ -1/2 (xµ) is the standardized version of the covariate x. For the multiple indices model (1), we then focus on the first K eigenvectors η 1 , . . . , η K associated with the largest K eigenvalues of the I p -symmetric matrix M * . Transforming back to the original scale, the vectors Σ -1/2 η k , k = 1, . . . , K are in the EDR subspace. Their estimation procedure is a straightforward replication of the previous estimation process using

M * = Σ -1/2 M Σ -1/2 .

RSIR: A modified estimated variance matrix

As previously mentioned, the rank of Σ implies that this matrix is singular when n < p and ill-conditioned when n ≈ p. To tackle these cases Zhong et al [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] introduce an upgrade of the SIR method, called RSIR, that relies on a modification of Σ such that the result can be inverted. This leads to the following estimate of Σ: where s is a positive real parameter and I p is the p × p identity matrix. For a given matrix A, let A 2 = Trace(A A). To find a suitable s, Zhong et al [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] propose to minimize the mean squared error

Σ(s) = Σ + sI p ,
L(s) = K ∑ k=1 Trace(V( βk (s))) + K ∑ k=1 E( βk (s)) -β k 2 ,
where β (s) = ( β1 (s), . . . , βK (s)) is the matrix of the K first generalized eigenvector of M and Σ(s), which is built such that the constraint βk (s) Σ(s) βk (s) = I[k = k] is verified for all (k, k) ∈ {1, . . . , K} 2 . More details about generalized eigenvectors can be found in Section 2.4. Because β k is unknown, the authors replaced it with E( βk (s 0 )) in the expression of L(s), to obtain an approximation L(s). Note that the parameter s 0 has to be sufficiently small in order for E( βk (s 0 )) to be close to β k . In practice, s 0 is chosen equal to 0. Variances and expectations in L(s) are then estimated with bootstrap samples, which leads to an estimate L(s) of L(s). Remark that estimating E( βk (s 0 )) for s 0 = 0, implies using SIR with Σ. To do so, Zhong et al [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] apply the QZ algorithm (see Section 2.4 for details). The optimal regularization parameter is then given by

s opt = arg min s L(s).
The corresponding matrix of estimated EDR directions is finally defined by βRSIR = β (s opt ).

SR-SIR: A ridge sliced inverse regression.

We describe here the SR-SIR method from Li and Yin [START_REF] Li | Sliced inverse regression with regularizations[END_REF]. When Σ is invertible, let ẑi = Σ -1/2 (x i -μ) for i = 1, . . . , n and let introduce the p × n matrix Z = (ẑ 1 , . . . , ẑn ). Let W = diag( p1 , . . . , pH ) and let S be a n × H matrix made of elements S i,h defined as

S i,h = 1 n I[y i ∈ s h ] ph .
We can thus write M * = Z S W S Z . Let η1 , . . . , ηH be the first H eigenvectors of M * . They also appear in the following singular value decomposition:

Z S W 1/2 = ( η1 , . . . , ηH )Λ 1/2 V ,
where Λ 1/2 and V are respectively an r × r matrix and an H × r matrix, and r is the rank of Z S W 1/2 . If r > K, let V be the K first columns of Λ 1/2 V and write η = ( η1 , . . . , ηK ). From Eckart-Young theorem we have η, V = arg min

u,v Z S W 1/2 -uv 2 , (5)
where the minimum is taken over the respective sets of p × K matrices and K × H matrices. If we write v = (v 1 , . . . , v H ) and z S = (z 1 , . . . , zH ), then we have It is thus easy to check that ( 5) is equivalent to

Z S W
η, V W -1/2 = arg min u,v H ∑ h=1 ph zh -uv h 2 , (6)
as described in [START_REF] Cook | Testing predictor contributions in sufficient dimension reduction[END_REF]. Li and Yin [START_REF] Li | Sliced inverse regression with regularizations[END_REF] showed that ( 6) is equivalent to

β , V W -1/2 = arg min u,v H ∑ h=1 ph ( mh -μ) -Σuv h 2 , (7)
where β = Σ -1/2 η. Note that this equation is also defined when Σ is not invertible. From [START_REF] Bottolo | ESS ++ : a C ++ objected-oriented algorithm for bayesian stochastic search model exploration[END_REF], the authors proposed thus a ridge version of the estimator β for a given regularization parameter s:

β (s), V (s) W -1/2 = arg min u,v G s (u, v). (8)
where

G s (u, v) = H ∑ h=1 ph ( mh -μ) -Σuv h 2 + s u 2
In practice, β (s) can be then obtained from [START_REF] Bottolo | Evolutionary stochastic search for bayesian model exploration[END_REF] with an alternating least-squares algorithm even when n < p. The SR-SIR method rely on a generalized crossvalidation criterion to find the optimal regularization parameter s opt (see [START_REF] Li | Sliced inverse regression with regularizations[END_REF] for details). Finally the matrix of estimated EDR directions is defined by βSR-SIR = β (s opt ).

Remark 1. The exisence of a solution for [START_REF] Bottolo | ESS ++ : a C ++ objected-oriented algorithm for bayesian stochastic search model exploration[END_REF] is not proved as explained by Bernard-Michel et al [START_REF] Bernard-Michel | A note on sliced inverse regression with regularizations[END_REF]. Indeed, assume that β (s), V (s) W -1/2 is such a solution and that β (s) is not the null vector, we then have

G s 1 2 β (s), 2 V (s) W -1/2 < G s β (s), V (s) W -1/2 ,
which contradicts the fact that β (s), V (s) W -1/2 verifies (7). This encourages Bernard-Michel et al [START_REF] Bernard-Michel | A note on sliced inverse regression with regularizations[END_REF] to replace [START_REF] Bottolo | ESS ++ : a C ++ objected-oriented algorithm for bayesian stochastic search model exploration[END_REF] with the following optimization problem:

β (s), V (s) W -1/2 = arg min u,v H ∑ h=1 ph ( mh -μ) -Σuv h 2 + s uv W 1/2 2 .
The value of β (s) in this problem is actually the estimate of the RSIR method, for a regularization parameter s.

SIR-QZ: Solving the generalized eigenvalues problem in SIR

When Σ is regular, usual SIR estimates of the EDR directions are eigenvectors of Σ -1 M. This eigendecomposition is actually a special case of a generalized eigenvalues problem which consists in finding real numbers λ and non-null vectors v such that: When Σ is singular the generalized eigenvalue problem can still be solved if the function λ → Mλ Σ behave properly. We call this function a matrix pencil. In this section, we present the QZ algorithm which allows us to find couples (λ , v) that verify (9) for a wide range of matrix pencils including some with singular matrices Σ. The QZ algorithm can be viewed as an extension of the QR algorithm and was proposed by Moler and Stewart [START_REF] Moler | An algorithm for generalized matrix eigenvalue problems[END_REF]. The reader can refer to chapter 7 of [START_REF] Golub | Matrix computations[END_REF] for details. A brief description of this algorithm is provided in the following.

Mv = λ Σv. (9)
Notice that if we have two invertible matrices Q and Z, then finding λ and v in ( 9) is equivalent to find λ and w in

Q MZw = λ Q ΣZw, (10)
and to set v = Zw. Similarly to the QR algorithm that is designed to find the Schur decomposition of a matrix in order to compute its eigenvalues, the QZ algorithm aims at finding unitary matrices Q and Z such that Q MZ and Q ΣZ are upper triangular, for square matrices M and Σ. Such a transformation is called a generalized Schur decomposition. When working with complex matrices, Q and Z always exist (see [START_REF] Golub | Matrix computations[END_REF], Theorem 7.7-1). Possible values of λ that verify ( 9) are such that det( Mλ Σ) = 0, and such that det(Q( Mλ Σ)Z) = 0. The latter determinant is the product of the diagonal elements of Q( Mλ Σ)Z since it is an upper triangular matrix. Hence, the generalized eigenvalues of ( 9) are the ratios of the diagonal elements of Q MZ to the ones of Q ΣZ, provided that the diagonal elements of Q ΣZ are not equal to zero. More specifically, this can be seen with the following formula ([19], Theorem 7.7-1)

det( M -λ Σ) = det(Q Z ) p ∏ j=1 (t j -λ u j ), (11)
where t 1 , . . . ,t p and u 1 , . . . , u p are the respective diagonal elements of Q MZ and Q ΣZ. Notice that the generalized Schur decomposition only produces complex upper triangular matrices Q MZ and Q ΣZ. However, there is a similar available decomposition for real matrices M and Σ (see Appendix B.1). Equation [START_REF] Carroll | Measurement error regression with unknown link: dimension reduction and data visualization[END_REF] implies that if it exists j ∈ {1, . . . , p} such that t j = u j = 0, then det( Mλ Σ) = 0 for all λ ∈ C, and trying to choose eigenvectors corresponding to the greatest eigenvalues to estimate the EDR directions does not make sense. Numerically, for any j ∈ {1, . . . , p}, due to rounding errors, t j and u j are almost always different from 0, but if both their absolute value are too small, det( Mλ Σ) is sufficiently unstable to call its value in question. As a consequence, every computed λ j can be wrong. For similar reasons, if |u j | is too small for a given j ∈ {1, . . . , p}, t j /u j should not be considered as an eigenvalue. These remarks and the regularization procedure in [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] lead to the algorithm of Appendix B.2 to find estimate βQZ of the EDR directions.

Let X = (x 1 , . . . , x n ). When n is sufficiently smaller than p and H > K, a generalized eigenvector v of ( 9) is such that the n indices X v only takes H distinct values, as explained in Appendix B.3. In practice, the regularization parameter s of the algorithm of Appendix B.2 is small and we can distinguish easily H clusters in the values of X βQZ in Figure 8 in Appendix B.3. This is a drawback of this approach when n < p since the values of X β are a priori distinct. To circumvent this shortcoming, we compute several βQZ with different number of slices H 1 , . . . , H N H . Let denote the corresponding estimates βQZ,1 , . . . , βQZ,N H . Thus, we choose γ as an estimate of X β with γ, δ = arg min γ,δ X ( βQZ,1 , . . . , βQZ,N H )γδ where the minimum is taken over the respective sets of n × K matrices and K × KN H matrices. The solution of ( 12) is then given by a principal component analysis. We call the whole approach SIR-QZ.

Remark 2. When n is smaller enough than p, finding a satisfying estimate of β may not be possible. For example, if K = 1, n < p and if the columns of X are not linearly dependent, there are infinitely many solutions u of the system X u = X β for a given estimate β of β . Recalling the underlying model [START_REF] Amato | Dimension reduction in functional regression with applications[END_REF], there is no reason why β should be a better estimate of β than any of these solutions. That is why, when n < p, we focus on estimates of X β rather than on β itself, in [START_REF] Chen | Can SIR be as popular as multiple linear regression?[END_REF].

SIR-MP: a generalization of the inverse for singular matrices

We describe in this section a method which mimics the SIR approach developed for a functional covariate.

In the functional SIR context, x is an explanatory functional variable (assumed square integrable in order to have its covariance operator well-defined) while y is still a real response variable. In this context, while the covariance operator of x is invertible, it has unbounded inverse so that its estimator is ill-conditioned. Then several methods have been proposed when the covariance operator does not need to be inverted.

One of them consists in using the eigendecomposition of M + Σ instead of Σ -1 M, where M + is the Moore-Penrose generalized inverse of M, also called Moore-Penrose pseudoinverse of M. In the particular context of functional sliced inverse regression, the reader can find a discussion on the fact that the eigenvectors of Σ -1 M are eigenvectors of M + Σ, in Ferré and Yao [START_REF] Ferré | Reply to the paper : "a note on smoothed functional inverse regression[END_REF] and references cited therein.

Let us now focus on an alternative approach introduced by Amato et al. [START_REF] Amato | Dimension reduction in functional regression with applications[END_REF]. They used the fact that Σ -1/2 MΣ -1/2 is a finite rank operator, where Σ (resp. M) stands here for the covariance operator of x (resp. E(x|T (y))) in this functional context. The eigenvectors of this operator are eigenvectors of Σ 1/2 M + Σ 1/2 . The authors claimed that the reason of their approach is that a smooth estimate of M produces more stable estimates of the eigenvalue decomposition of M than that of the empirical estimate of Σ. Thus the eigenfunctions η 1 , . . . , η K associated with the smallest K eigenvalues α 1 , . . . , α K of Σ 1/2 M + Σ 1/2 are also the eigenfunctions associated with the largest eigenvalues of Σ -1/2 MΣ -1/2 equal to 1/α k for k = 1, . . . , K. In order to transform back to the original scale, we can not use the transformation Σ -1/2 η k . A basis of the (functional) EDR space is instead given by

b k = M + Σ 1/2 η k for k = 1, . . . , K. (13)
We provide in Appendix A a brief proof of this result.

In the context of our paper (that is, n < p < ∞), we will evaluate how the functional SIR procedure behaves in the multivariate framework. To do this, we simply substitute the operators of covariance M and Σ by the estimates M and Σ previously defined in (2) and (4). The resulting estimated directions are: where the η k 's are the eigenvectors of Σ 1/2 M + Σ 1/2 associated with the smallest eigenvalues (among those not structurally equal to zero, see Remark 3 for details).

b k = M + Σ 1/2 η k for k = 1, . . . , K, (14)
Remark 3. For a p-dimensional covariate, the p × p matrix Σ 1/2 M + Σ 1/2 is symmetric positive semidefinite and its rank r is at most equal to H -1 when H < n < p. Therefore, the eigenvalues of Σ 1/2 M + Σ 1/2 are such that α1 ≥ • • • ≥ αr > 0, and the geometric multiplicity of the eigenvalue zero is equal to pr by construction. Thus we are interested in the eigenvectors η k associated with the K eigenvalues αr , . . . , αr-K+1 .

3. Selecting relevant components of x which are linked with y Let β j,k denote the jth element of the EDR direction β k , for k = 1, . . . , K and j = 1, . . . , p. If β j,. = (β j,1 , . . . , β j,K ) is the null vector then the jth component of x does not have any effect on y. Finding such components is an important concern when n < p because it allows x ∈ R p to be reduced to x ∈ R p , where p < p, without any loss of information. If in addition p is less enough than n, the EDR directions can then be accurately estimated with a classical SIR procedure applied on y and x . In Section 3.1, we describe the method from [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] to determine which β j,. are null. In Section 3.2, we explain the Lasso procedure detailed in [START_REF] Li | Sliced inverse regression with regularizations[END_REF] to obtain such conclusions. In Section 3.3, we introduce another method to solve this problem based on proximity measures between models with only a few components of x and the initial model (in which every component of x is taken into account).

RSIR: Bootstrap estimates and a chi-squared test

Let βRSIR, j,k be the elements of the matrix βRSIR . Zhong et al [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] claim that for j = 1, . . . , p, the vector βRSIR, j,. = ( βRSIR, j,1 , . . . , βRSIR, j,K ) follows asymptotically a multivariate normal distribution with mean β j,. (s opt ) and covariance matrix Γ j . Provided that Γ j can be inverted, if

β j,.
(s opt ) is the null vector then β RSIR, j,. Γ -1 j βRSIR, j,. follows asymptotically a chi-squared distribution with K degrees of freedom. This encouraged Zhong et al [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] to use a chi-squared test on β RSIR, j,. Γ -1 j βRSIR, j,. to select which components of x have effect on y, where Γ j is an estimate of Γ j computed from bootstrap estimates of β j,. (s opt ).

Difficulties using this procedure could arise when n < p because the distribution under the null hypothesis of this test is asymptotic. In addition, it leads to inferences about the vectors β j,. (s opt ) for j = 1, . . . , p which are not necessary the same than for β j,. .

SR-SIR: A Lasso method

From β (s opt ) and ( ṽ1 , . . . , ṽH ) = V (s opt ) P -1/2 , given in equation ( 8), Li and Yin [START_REF] Li | Sliced inverse regression with regularizations[END_REF] propose to minimize the following expression under a constraint on the L 1 -norm of the vector φ :

G(φ ) = H ∑ h=1 ph ( mh -μ) -Σdiag(φ ) β (s opt ) ṽh 2 .
Soumis au Journal de la Société Française de Statistique File: "SIRn<p-2012 _v2".tex, compiled with jsfds, version : 2009/12/09 date: November 30, 2012 This leads to the following optimization problem for a parameter τ > 0:

φτ = arg min φ {G(φ )} , s.t. |φ | ≤ τ
where the minimum is taken over the set of vectors φ of length p. The Lasso procedure [START_REF] Tibshirani | Regression Shrinkage and Selection via the Lasso[END_REF] can be used to find φτ . For most positive values of τ 1 and τ 2 such that τ 1 > τ 2 , the number of zeros in φτ 1 is less than the number of zeros in φτ 2 . An example for which it is not true can be found in [START_REF] Efron | Least angle regression[END_REF], Section 3.1.

Li and Yin [START_REF] Li | Sliced inverse regression with regularizations[END_REF] consider that a component of x that corresponds to a zero in φτ does not have any effect on y. Let p τ be the number of non-null elements of φτ . In practice, choosing τ implies choosing the amount of selection provided by φτ . To do so, Li and Yin [START_REF] Li | Sliced inverse regression with regularizations[END_REF] propose to use classical model selection criteria. More specifically, this involves minimizing one of the following expression over a set of tested values of τ:

AIC = pH log G( φτ ) pH + 2p τ , BIC = pH log G( φτ ) pH + log(pH)p τ , RIC = (pH -p τ ) log G( φτ ) pH -p τ + p τ (log(pH) -1) + 4 pH -p τ -2 .

CSS based on SIR: Closest submodel selection for SIR methods

The idea of the procedure described here is to select submodels of (1) with only a given number p 0 of components of x which are the closest to the initial one. The components of x that appear the most in these submodels are asserted to have an effect on y.

Let Y = (y 1 , . . . , y n ) . To do this, we propose the following algorithm.

Initialize p 0 ∈]1, p[, N 0 ∈ N * and ζ ∈]0, 1[ or ρ ∈]0, 1[.
Step 1. Compute the estimated indices γ ∈ R n on Y and the whole covariate matrix X using SIR-QZ. Let a = 1.

Step 2. Select randomly p 0 components of x and build the corresponding matrix X For example, in our simulation study, we set N 0 = 10 4 and ζ = 10% to determine the closest N 1 = ζ N 0 submodels. In our real data application, we use N 0 = 9 × 10 5 and ρ = 0.9 to select the top N 1 submodels.

Note that choosing p 0 < n allows us to use classical SIR instead of SIR-QZ in Step 3 which significantly improves the computational time. In addition, any SIR approach that provides estimates of the indices (when n ≤ p and n > p) could be used in the whole algorithm instead of SIR-QZ.

A simulation study

In Sections 2-3, we presented 4 methods to estimate EDR directions (or indices) and 3 procedures to select which components of x have effects on y. In Section 4.1, we illustrate them on a single simulated data set. To compare their numerical performances, we then study them on several replications in Section 4.2.

Analysis of a single data set

Simulated model

We consider the following single index model

y = (x β ) 3 + ε (15)
where x and β are a p-dimensional vectors defined hereafter and ε ∼ N (0, σ 2 ) with σ = 10 -3 . Let p = 200 and p = 20. We choose β = 1 10 (I(1 ≤ p ), . . . , I(p ≤ p )) , so that p is the number of non-null components of β , that is the number of components of x that affect y. We construct x = (x 1 , . . . , x p ) as follows: for j = 1, . . . , p , generate σ 2 j from the law U ([0.05, 0.1]) and x j from the law N (0, σ 2 j ). For j = p + 1, . . . , p, set

σ 2 j = 12 -( j -1)/p ( j -1)/p 2 σ 2 ( j-1) mod p +1 ,
when mod denotes the modulo operation. Generate then x j from the law N (0, σ 2 j ) and set x j = x ( j-1) mod p +1 + x j . This ensures that cor(x j , x ( j-1) mod p +1 ) = ( j -1)/p /12.

Estimation of EDR indices

We simulate an independent and identically distributed sample (X,Y ) of size n = 100 from model [START_REF] Duan | Slicing regression: a link-free regression method[END_REF]. We plot Y versus the true indices X β in Figure 1. We analyze X and Y with the various methods presented in Section 2.2-2.5.

-For the RSIR method we evaluate L(s) for s ∈ {0, 10 -10 , 10 -9 , . . . , 10 5 } with 50 bootstrap samples and H = 10. In Figure 2(a) we plot the values of Y against the indices provided by the RSIR method. The structure of Figure 1 can not be discerned in Figure 2(a). The regularization parameter that RSIR provides is equal to 10 5 and thus the RSIR procedure is equivalent to an eigendecomposition of M. -Concerning the SR-SIR method, the chosen regularization parameter s opt is chosen in {10 -10 , 10 -9 , . . . , 10 5 }. The number of iterations of the alternating least square algorithm of SR-SIR is set to 50 and we take H = 10. For this example, we find s opt = 10 3 . In Figure 2(b), we draw the values of Y against the estimated indices X βSR-SIR . The points of this graphic do not form the same shape as the points in Figure 1. -We run SIR-QZ for {H 1 , . . . , H N H } = {5, . . . , 15}. In Figure 2(c), we plot y against the corresponding estimated indices. This graphic exhibits a structure which is similar to the one in Figure 1. -We finally apply SIR-MP with H = 10. We observe in Figure 2(d), which shows how Y and the indices produced by the SIR-MP method are related, that this method also fails to recover the shape of Figure 1. Thus, Figure 2 shows that for this data set, SIR-QZ provides better estimations of the indices than RSIR, SR-SIR and SIR-MP.

To quantify such conclusions, we can use a criterion that measures how X β and X β are close form each other, for a given estimate β . Let X = ( x1 , . . . , xn ) and let P be the projector on the subspace of R n spanned by X β . More precisely, we have

P = X β (β X X β ) -1 β X. (16)
Similarly, define P RSIR , P SR-SIR , and P SIR-MP by respectively replacing β by βRSIR , βSR-SIR , and βSIR-MP in [START_REF] Efron | Least angle regression[END_REF]. Let us also define P SIR-QZ by replacing X β by Īn γ in [START_REF] Efron | Least angle regression[END_REF], where Īn is a matrix that centers γ (see Appendix B.3). For a given method m, we use the squared trace correlation between the subspaces spanned by X β and by X βm as a measure of the closeness between these subspaces. It is defined by

R(m) = 1 K Trace(PP m ). (17)
Notice that if R(m) is the squared cosine of the angle between the vectors X β and X βm . This quality measure belongs to [0, 1] and the higher its value is, the better the indices are estimated. In Table 1, we present values of R(m) for the four considered methods. For SIR-QZ, R(m) is clearly higher than for RSIR, SR-SIR and SIR-MP (0.74 versus less than 0.1). Notice also that the computational time is a lot greater for RSIR and SR-SIR than for SIR-MP and SIR-QZ.

Selection of components of x

We rely on the true positive rate (TPR) and on the false positive rate (FPR) to evaluate procedures that find which elements of β are equal to 0. The TPR is the number of selected components of x that actually affect y divided by the total number of components of x that affect y. The FPR is the number of selected components of x that do not affect y divided by the total number of components of x that do not affect y.

For the RSIR selection method the returned p-values are ordered and the components that correspond to the first p-values are selected. For the CSS procedure, the components related to the greatest number of occurrence are selected. We evaluate the results of both methods by selecting the best 10 (resp. 20, 50, 100, 150) components. For the SR-SIR method, we use the different criteria (AIC, BIC and RIC) proposed in [START_REF] Li | Sliced inverse regression with regularizations[END_REF] to determine the appropriate Lasso parameter.

The number of bootstrap sample generated for the RSIR method is set to 10 3 . It is great enough so that increasing it does not improve significantly the quality criteria. This method relies on a regularization parameter which can also be provided by the Algorithm 1 to find βQZ (see Appendix B.2). The estimate βQZ can thus be plugged in the RSIR selection method. For the CSS method, we choose N 0 = 10 4 , ζ = 10% and p 0 = 50. While increasing N 0 could lead to better quality criteria, the computational time is sufficiently large not to choose it greater. The tested values of the Lasso parameter are the SR-SIR method are in {1, 2, . . . , 100}. This algorithm needs an estimate of β in input. We use βSR-SIR but, because of the poor results of Table 1 for this estimate, we also take βQZ with H = 10, and the true EDR direction β .

Results of the corresponding TPR and FPR are displayed in Table 2. The SR-SIR method performs very well if an accurate estimate of β is provided, while the results are really bad otherwise because no selected component of x has any effect of y. Concerning RSIR and the CSS method, their FPR are similar, but the TPR are greater for the latter. The results for RSIR with β QZ are slightly better than with the full RSIR procedure. The CSS method also seems to need good estimates of the indices since working with the true ones produces better TPR than using γ.

To get more insights about the numerical performances of the various procedures tested in this example, we run them in several replications in the following section.

General behaviors of the estimates over several replications

We generate 100 samples of size n = 100 from the model [START_REF] Duan | Slicing regression: a link-free regression method[END_REF]. For each of them, we launch the RSIR, SR-SIR, SIR-QZ and SIR-MP procedures with the same parameters as in Section 4.1.2, and compute the quality criterion given in [START_REF] Ferré | Determining the dimension in sliced inverse regression and related methods[END_REF]. We display the values of the criterion in Figure 3. The trend that is exhibited in Table 1 is confirmed in this graphic since the values of R(SIR-QZ) are clearly greater than the others. Various methods to select important components of x are then run in the 100 samples drawn from model [START_REF] Duan | Slicing regression: a link-free regression method[END_REF]:

-RSIR with βRSIR , 1000 bootstrap samples and the following significance levels of the corresponding test: 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, -SR-SIR with βSR-SIR and τ ∈ {1, 2, . . . , 100}, -The CSS procedure with γ from SIR-QZ and parameters N 0 = 10 4 , ζ = 10% and p 0 = 50 for various number of selected components: 10, 20, 30, 40, 60, 80, 100, 120, 140, 160. The mean ROC curves over the 100 replications are displayed in Figure 4. The CSS method outperforms RSIR while SR-SIR provides poor results. Notice that for RSIR, the values of the FPR are close to the chosen levels of test, in spite of the fact that this test is asymptotic.

Real data application

Description of the Dataset

We illustrate our developed approach on a genetic dataset which contains transcriptomic data and genomic data. In this study, we aim at finding genetic causes of variation in the expression of genes, that is eQTL (expression Quantitative Trait Loci). In this context, the gene expression data are treated as a quantitative phenotype and the genotype data (SNPs) are considered as predictors. In this illustration, we study the Hopx gene analyzed in [START_REF] Petretto | New insights into the genetic control of gene expression using a bayesian multi-tissue approach[END_REF]. We investigate the ability of SIR-QZ combined with the CSS selection procedure to find a parsimonious index that explains the variability of Hopx gene expression in the heart tissue using p = 770 SNPs from n = 29 inbred line rats.

SIR-QZ and CSS results

We first run SIR-QZ for {H 1 , . . . , H N H } = {2, . . . , 6}. In Figure 5(a), we plot the dependent variable Hopx versus the index based on the whole set of SNPs (p = 770). This graphic clearly exhibits a link between the phenotype and the index estimated by a smooth kernel method. In this illustration, this link is almost linear. In contrast, the plot (not given here) of Hopx gene versus the second EDR index does not show any structure. From this graphical diagnostic, it appears that only one EDR direction provides relevant information to explain the variability of the gene expression.

Then, in order to find a parsimonious index, we run our CSS selection procedure with N 0 = 900.000. The examined values of p 0 are in {10, 20} while ρ takes value in {0.75, 0.80, 0.85, 0.90}. In Table 3, we present the number of selected SNPs for each combination of this two parameters. The threshold used for the selection will be detailed below. We can observe that, not surprisingly, the numbers of selected SNPs increases with p 0 and decreases with ρ. Moreover for a given value of p 0 , we specify, in this table, the number of selected SNPs in common with those selected with ρ = 0.9 (corresponding to the parsimonious model). We also indicate, for a given ρ, the number of SNPs in common with those selected when p 0 = 10 and when p 0 = 20. This table highlights an overlap of 10 SNPs among all the sets of the selected SNPs for the various couple (p 0 ,ρ). Note that, the smallest set contains 11 SNPs when p 0 = 10 and ρ = 0.9 which comforts us about the stability of the CSS procedure.

In eQTL study, it is known that only a few number of SNPs can explain the variation of the gene expression. Thus, from the expertise of the biologists, we decide to select the sparsest model, that is with (p 0 , ρ) = (10, 0.9). Figure 6 exhibits the selected 11 SNPs for this choice of p 0 and ρ. The threshold (horizontal red line in the figure) is defined as follows: N 1

p 0 p + u 1-α/2 p N 1 p 0 p (1 -p 0 p ) where u 1-α/2 p is the quantile of order (1 -α/2
p ) of the standard normal distribution. It corresponds to the upper bound of the prediction interval of the occurrence of a SNP in the selected model under the hypothesis that none of the SNPs are associated with the gene expression. The level of this interval is fixed at 1α = 0.95 and is corrected by a Bonferroni approach.

On Figure 5(b), we plot the dependent variable Hopx versus the index based on the 10 SNPs selected according to our previous comments on Table 3. The linear correlation between this index and the one estimated on all the SNPs which is equal to 0.956, highlights the good behaviour of our CCS strategy to select the relevant SNPs. Thus, not surprisingly, we observe the same relation between Hopx gene expression and the estimated indices.

Comparison methods

We compare our approach with three popular multivariate methods for analyzing high-dimensionnal datasets: A Lasso approach, the sparse Partial Least Squares (sPLS) and a Bayesian variable selection regression (ESS++). The Lasso method [START_REF] Tibshirani | Regression Shrinkage and Selection via the Lasso[END_REF], which is a penalized regression method proposed to improve OLS (ordinary Least square), often performs poorly in both prediction and interpretation especially when n is small and p is large. This technique tends to shrink the regression coefficients towards zero in order to select a sparse subset of covariates and provide a better prediction performance. sPLS [START_REF] Cao | Sparse canonical methods for biological data integration: application to a cross-platform study[END_REF] seeks for the best linear combination of SNPs to predict the outcome. To ensure sparsity, sPLS includes a penalty function on some loading coefficients which is equivalent to a restriction on the number of loading vectors and on the number of SNPs, in each vector, that have a non-null coefficient. Both Lasso and sPLS approaches require a preliminary calibration of the tuning parameters which directly affects the number of selected variables, the estimate of the model parameters and therefore the statistical performances of the models. Calibration procedures usually involve the minimization of the mean square error of prediction through V-fold cross validation. In this illustration, we used the leave-one-out crossvalidation method to choose the tuning parameter for both methods. We finally compare our results with ESS++ a Bayesian variable selection approach for linear regression that can analyze single and multiple responses ( [START_REF] Bottolo | Evolutionary stochastic search for bayesian model exploration[END_REF][START_REF] Bottolo | ESS ++ : a C ++ objected-oriented algorithm for bayesian stochastic search model exploration[END_REF]). ESS++ exploits recent developments in MCMC search algorithms to explore the 2 p -dimensional model space. The performances of this method have been, among others, illustrated on eQTL studies ( [START_REF] Petretto | New insights into the genetic control of gene expression using a bayesian multi-tissue approach[END_REF]). Figure 7 presents the Venn diagram of the sets of SNPs selected by the different approaches. Two SNPs (D14Mit3 and D2Cebr204s17) are selected by the four methods. D14Mit3 (chromosome 14) is clearly the first SNP in the list of the SNPs selected by the CSS procedure (see Figure 6) and D2Cebr204s17 (chromosome 2) is at the third position. Moreover, our proposed approach reveals 4 SNPs (D18Ucsf1, Fgg, D2Utr9, D16Rat75) not selected by the other methods. Although three of them are very close to our proposed threshold (red line in Figure 6), while the SNP D18Ucsf1 (chromosome 18) is clearly selected by our procedure.

The main advantage of our approach is the opportunity to reveal a non-linear link between the gene expression and a parsimonious index while the other compared approaches are based on a linear model. However, theses methods could treat multiple correlated phenotypes (multiple continuous responses). For example, ESS++ has been used to study the joint variability of gene expression in seven tissues (adrenal, gland, aorta, fat, heart, kidney, liver, skeletal) from inbred line rats [START_REF] Bottolo | ESS ++ : a C ++ objected-oriented algorithm for bayesian stochastic search model exploration[END_REF]. An extension of SIR-QZ for multivariate response is under investigation.

Concluding remarks

Although regularizing the estimated covariance matrix or constraining the optimization problem are natural ways to extend SIR to underdetermined cases (n < p), it may not be clear which one should be chosen and how to set the related parameter for each procedure. For the RSIR method, we illustrated that in such a context, the corresponding parameter should rather be determined with respect to the stability of the linear algebra algorithm (as in SIR-QZ) than with a statistical criterion. Moreover, the SIR-QZ approach introduced in this paper produces better results in simulation than the SR-SIR method that constraint the underlying optimization problem. In addition, the poor performances of SIR-MP suggest that adapting properties of the pseudo-inverse from the functional SIR to our high dimensional context is not well-adapted. We assumed that the dimension K of the EDR subspace was known in our simulation study. While, in the application, an empirical argument was given to determine this dimension, its estimation remains to be done with care. We also proposed the CSS method that searches which submodels are the most informative to select relevant components of x. This procedure relies on a quality measure of the estimated indices, and it outperforms, in simulation, RSIR and SR-SIR selection procedures which are both based on estimates of β . We thus explain these results by pointing out that, when n < p, only the indices can be properly estimated, in general cases. Note that the space of the submodels may not be browsed optimally by our proposed CSS algorithm. Improvements should be made using optimization techniques such as genetic algorithms.

All the developed methods have been implemented in R language and are available upon request from the corresponding author. An R package is currently under development.

The illustration of genetic data highlights the opportunity of SIR-QZ combined by the CSS procedure to reveal a few number of SNPs which can explain the variability of the expression of the Hopx gene. Then, a non linear link between the gene expression and the parsimonious index could be estimated. The choice of the number of SNPs to keep is still a topic of concern since alternatives to our threshold could be considered.

In genetic datasets, the response variable is often multivariate. For instance, it could represent several phenotypes as in eQTL studies. Some approaches already handle such datasets. Since univariate results for SIR-QZ and the CSS selection procedure are consistent with the other methods presented in our application, it thus appears interesting to extend them to the multivariate case.

be the element of Q ΣZ at the same location than t j in Q MZ and define similarly ũ j 1 , j 2 for each t j 1 , j 2 . Thus, we have

det( M -λ Σ) = det(Q Z ) ∏ j∈J ( t j -λ ũ j ) ∏ j∈J c
det t j, jλ ũ j, j t j, j+1λ ũ j, j+1 t j+1, j t j+1, j+1λ ũ j+1, j+1 .

Hence, for j ∈ J, if ũ j = 0, then λ j = t j / ũ j is a real generalized eigenvalue. In addition, for j ∈ J c , Moler and Stewart [START_REF] Moler | An algorithm for generalized matrix eigenvalue problems[END_REF] succeeded in finding ( t j , t j+1 ) ∈ C 2 , and ( ũ j , ũ j+1 ) ∈ R\{0} such that λ j = t j / ũ j and λ j+1 = t j+1 / ũ j+1 are generalized eigenvalues of M and Σ. This leads to vectors t = ( t1 , . . . , tp ) and ũ = ( ũ1 , . . . , ũp ) that are sent back by the QZ algorithm in order to provide generalized eigenvalues.

B.2. Algorithm

As explained in Section 2.4, the QZ algorithm has to be controlled when dealing with singular pencils. The following pseudocode in Scilab language allows the user to do so in the context of sliced inverse regression, for underdetermined cases. Typical values for s min , c and ε chosen in the simulation study of Section 4 are respectively 10 -16 , 10 and 10 -10 . The QZ algorithm is implemented in Matlab through the eig function. With Scilab, one should use spec, which is based on the LAPACK library. The R software is able to call functions from this library. Thus the QZ algorithm can be easily tested with this software. Notice that Algorithm 1 is designed to handle real values of ũ j and t j but, as mentioned in Appendix B.1, they can be complex. In that case, knowing if the blocks made of t j, j , t j+1, j , t j, j+1 and t j+1, j+1 and Soumis au Journal de la Société Française de Statistique File: "SIRn<p-2012 _v2".tex, compiled with jsfds, version : 2009/12/09 date: November 30, 2012 of ũ j, j , ũ j+1, j , ũ j, j+1 and ũ j+1, j+1 produce unstable eigenvalues is more difficult. As explained in Section 5 of [START_REF] Moler | An algorithm for generalized matrix eigenvalue problems[END_REF], the QZ algorithm aims at finding stable λ j and λ j+1 corresponding to these 2 × 2 blocks. Because we do not control this procedure, we simply report if the QZ algorithm send back complex values in t. We never encounter this case in the simulation study of Section 4.

B.3. The sliced indices issue

Hereafter, we describe why the Algorithm 1 produces clustered indices as in Figure 8. Recall that Σ = 1 n X X . Define Īn = I n -1 n 1 n,n , where every element of the n × n matrix 1 n,n is equal to 1. Notice that X = X Īn and then Σ = 1 n X Īn X . We can also write M = X Īn S W S Īn X , where S and W are defined in Section 2.3. Because of the structure of S, for any H × α matrix A, SA has at most H distinct rows, so has Īn SA. Let w be the first generalized eigenvector of Īn S W S Īn and Īn associated with the eigenvalue λ . This means that Īn S W S Īn w = λ Īn w which implies that w = Īn w has at most H distinct values and then that w has also at most H distinct values.

Assume that X has full column rank, which is likely to happen when p > n. Then, X + w is a generalized eigenvector of M and Σ and X X + w = w has at most H distinct values. The eigenvalue that is related to X + w is equal to nλ .

If it exists β1 = X + w, a generalized eigenvector of M and Σ such that the generalized eigenvalue which is related to β1 is greater than the one corresponding to X + w, then we should have

2 ( 12 )

 212 Soumis au Journal de la Société Française de Statistique File: "SIRn<p-2012 _v2".tex, compiled with jsfds, version : 2009/12/09 date:November 30, 2012

(a) . Step 3 .

 3 Compute the SIR-QZ indices γ(a) ∈ R n based on Y and X (a) . Step 4. Calculate the linear correlation between the indices γ and γ(a) . Let us denote by ĉ(a) the square of this correlation. Let a = a + 1. Repeat N 0 times steps 2-4. Step 5. Consider the submodels corresponding to the N 1 largest correlations ĉ(a) .Either the user set ζ ∈]0, 1[ and then gets N 1 = ζ N 0 , or the user chose a value for ρ and then N 1 is the number of submodels such that ĉ(a) > ρ.Step 6. Count the number of occurrences of each component of x in these N 1 submodels. The components that affect y are the ones that have the greater number of occurrences.

FIGURE 1 .

 1 FIGURE 1. Plot of Y versus X β generated from model (15), with n = 100 and p = 200. The horizontal scale was standardized.

FIGURE 2 .

 2 FIGURE 2. Plot of Y versus estimates indices obtained with (a) the RSIR method, (b) the SR-SIR method, (c) the SIR-QZ method and (d) the SIR-MP method. The horizontal scale was standardized.

  Soumis au Journal de la Société Française de Statistique File: "SIRn<p-2012 _v2".tex, compiled with jsfds, version : 2009/12/09 date:November 30, 2012

FIGURE 3 .

 3 FIGURE 3. Boxplots of 100 values of R(m) for various estimation procedures m and samples of size n = 100 generated from model (15) with p = 200.

FIGURE 4 .

 4 FIGURE 4. Mean ROC curves for various procedures to select components of x that affect y in the model of Section 4.1, over 100 replications, with n = 100 and p = 200. Solid circles: RSIR, crosses: SR-SIR, empty circles: SIR-CSS.

FIGURE 5 .

 5 FIGURE 5. Plots of the dependent variable Hopx versus the index based (a) on the whole SNPs, (b) on the 10 selected SNPs. The linear correlation between these two indices (evaluated on n = 29 rats) is 0.956. The dotted (red) line is the estimated linear model, the solid (blue) line is the kernel estimate of the link function with a bandwidth chosen by cross-validation.

FIGURE 7 .

 7 FIGURE 7. Venn diagram of the sets of SNPs selected by Lasso, sPLS, ESS++ and SIR-QZ (combined with CSS) approaches

β 1 X 1 X

 11 Īn S W S Īn X β1 β Īn X β1 > w Īn S W S Īn w w Īn w .But, because w is the first generalized eigenvector of Īn S W S Īn and Īn , it maximizes u Īn S W S Īn u u Īn u over the vectors u of length n, which contradicts the latter equation. Hence, such a β1 does not exist, and the first generalized eigenvector of M and Σ is X + w.In this paragraph, we show that it exists H -1 orthogonal vectors w 1 , . . . , w H-1 such that, for k = 1, . . . , H -1,w k Īn S W S Īn w k w k Īn w k = w Īn S W S Īn w w Īn w ,which means that the first K generalized eigenvectors of M and Σ are the vectors X + w k for k = 1, . . . , K ≤ H -1. We assume n > H. Let sort Y increasingly and reorder the columns of X such that each column corresponds to the appropriate element of Y . This transformation implies that S W S is block diagonal with H blocks. For h = 1, . . . , H, the size of the block h is equal to n ph and each element it contains is equal to 1n 2 ph. Hence, the rank of each block is equal to 1 and each block provides a positive eigenvalue for S W S , which is equal to n phn 2 ph = 1 n . The corresponding eigenvector is made of elements I[y i ∈s h ] √ n ph for i = 1, . . . , n.We have now H orthonormal eigenvectors of S W S for the eigenvalue 1 n , and then we can find H -1 orthonormal centered eigenvectors w1 , . . . , wH-1 of S W S for this eigenvalue. In addition, for k = 1, . . . , H -1, wk is also a generalized eigenvector of Īn S W S Īn and Īn because Īn wk = wk , and wk maximizes u Īn S W S Īn u u Īn u , over the vectors u of length n. Finally, we have that the first K generalized eigenvectors Soumis au Journal de la Société Française de Statistique File: "SIRn<p-2012 _v2".tex, compiled with jsfds, version : 2009/12/09 date: November 30, 2012

FIGURE 8 .

 8 FIGURE 8. Plot of Y versus X βQZ with H = 10. The horizontal scale was standardized.

TABLE 1 .

 1 Quality measure and computational time of various estimates of the indices X β for the simulated sample of size n = 100, with p = 200.

	m	RSIR	SR-SIR	SIR-QZ SIR-MP
	R(m)	0.051	0.088	0.741	0.000
	Computational time (s)	101.88 1,234.70	7.85	0.25

Soumis au Journal de la Société Française de Statistique File: "SIRn<p-2012 _v2".tex, compiled with jsfds, version : 2009/12/09 date:

November 30, 2012

TABLE 2 .

 2 True positive rate (TPR), false positive rate (FPR) and computational time of various methods run on the simulated sample of size n = 100, with p = 200, to determine which components of β of model (15) are not null.

	Soumis au Journal de la Société Française de Statistique
	File: "SIRn<p-2012 _v2".tex, compiled with jsfds, version : 2009/12/09
	date: November 30, 2012

TABLE 3 .

 3 Plot of the occurrence of the SNPs by the CSS procedure (p 0 = 10, ρ = 0.9). The horizontal solid red line represent the threshold to select the most relevant SNPs. Results on selected SNPs for various values of p 0 and ρ

	Chromosome

  Soumis au Journal de la Société Française de Statistique File: "SIRn<p-2012 _v2".tex, compiled with jsfds, version : 2009/12/09 date:November 30, 2012

  Soumis au Journal de la Société Française de Statistique File: "SIRn<p-2012 _v2".tex, compiled with jsfds, version : 2009/12/09 date:November 30, 2012

  Because it is based on generalized real Schur decompositions, the notations involved are related to Appendix B.1 rather than Section 2.4. a s t r u n o f t h e QZ a l g o r i t h m t h a t c o r r e s p o n d s t o t h e K / / g r e a t e s t v a l u e s o f t/ ũ . ALGORITHM 1. A procedure to estimate EDR directions with the QZ algorithm

	e n
	k e e p G o i n g = %F;
	e l s e
	s = s * c ;
	end
	end
	/ / The e s t i m a t e d EDR d i r e c t i o n s a r e t h e e i g e n v e c t o r s s e n t by t h e
	/ / l

/ / I n i t i a l i z e s min , c and ε . s = s min ; k e e p G o i n g = %T; w h i l e k e e p G o i n g / / Use t h e QZ a l g o r i t h m on M and Σ(s) t o f i n d / / v e c t o r s ũ and t . i f ( sum ( abs ( ũ ) < ε ) + sum ( abs ( t ) < ε ) == 0 ) & ( l e n g t h ( ũ ) -sum ( abs ( ũ ) < ε ) >= K ) t h

Acknowledgements

The authors wish to gratefully thank Marie Chavent for the first discussions that led to this work.

Appendix A: Proof of (14)

We have, for k = 1, . . . , K, Σ 1/2 M + Σ 1/2 η k = α k η k . By pre-multiplying by the matrix (M + Σ 1/2 ), we obtain:

Using the definition of b k given in [START_REF] Cook | Testing predictor contributions in sufficient dimension reduction[END_REF], we get:

From the comments on the functional SIR context, provided in Section 2.5, the proof is complete.

Appendix B: Details about SIR-QZ

B.1. Generalized real Schur decomposition

We work here with real matrices M and Σ. Similarly to the generalized Schur decomposition introduced in Section 2.4, that produces complex matrices Q and Z, the generalized real Schur decomposition (see Theoresm 7.7-2 of [START_REF] Golub | Matrix computations[END_REF]) ensures that we can find Q and Z such that Q MZ is an upper quasi-triangular real matrix and Q ΣZ is an upper triangular real one. An upper quasitriangular matrix can be defined as the sum of an upper triangular matrix and of a block diagonal matrix where the sizes of the block are 1 × 1 or 2 × 2. For a 1 × 1 diagonal block of the matrix Q MZ, its unique element is called t j if it is located at the jth row and at the jth column of Q MZ. We write t j 1 , j 2 the element of a 2 × 2 diagonal block, in the j 1 th row and in the j 2 th column of Q MZ. An example of such an upper quasi-triangular matrix is given below:

, where * denotes some real values. Let J be made of the elements j ∈ {1, . . . , p} such that t j exists and J c be the set made of j ∈ {1, . . . , p -1} such that t j, j and t j+1, j+1 exist. For each j ∈ J, let ũ j