J Dedecker

O Guedj

M L Taupin

Asymptotic confidence interval for R 2 in multiple linear regression

Keywords: Mathematics Subject Classification 2020 : 62H20, 62J05 Multiple correlation coefficient, Asymptotic distribution, Robustness, Heteroscedasticity

Following White's approach of robust multiple linear regression [15], we give asymptotic confidence intervals for the multiple correlation coefficient R 2 under minimal moment conditions. We also give the asymptotic joint distribution of the empirical estimators of the individual R 2 's. Through different sets of simulations, we show that the procedure is indeed robust (contrary to the procedure involving the near exact distribution of the empirical estimator of R 2 is the multivariate Gaussian case) and can be also applied to count linear regression.

Introduction

Let (Y, X (1) , . . . X (p) ) be a random vector with value in R p+1 . The variable Y is the response variable, and (X (1) , . . . , X (p) ) is the vector of explanatory variables. Assume that all the variables are square integrable. The multiple correlation between Y and (X (1) , . . . , X (p) ) can be expressed using the orthogonal projection of Y onto the L 2 -subspace V generated by the variables (X (0) , X (1) , . . . , X (p) ), where we denote by X (0) the constant variable X (0) ≡ 1. Let then Proj V (Y ) be the orthogonal projection of Y onto V . The response variable Y is uncorrelated to the vector (X (1) , . . . , X (p) ) if and only if Proj V (Y ) = E(Y ).

To measure the strength of the correlation between Y and (X (1) , . . . , X (p) ), the usual measure is the coefficient R 2 defined by (assuming that Var(Y ) > 0):

R 2 = Var(Proj V (Y )) Var(Y ) .
The coefficient R 2 is often referred as "the population variance-accounted-for effect size". As we see, it is defined as the proportion of the variance of Y that is explained by the best linear predictor in L 2 based on the variables 1, X (1) , . . . , X (p) .

1

In this paper, we consider the empirical estimator R 2 of the coefficient R 2 based on n independent copies (Y i , X

i , . . . , X (p) i ) 1≤i≤n of the vector (Y, X (1) , . . . , X (p) ). We prove its consistency, asymptotic normality, and we give an asymptotic confidence interval for R 2 . Our approach is a robust one: we will make no additional assumptions on the distribution of the vector (Y, X (1) , . . . , X (p) ), except moment assumptions necessary for the Central Limit Theorem (CLT) to hold, and to estimate consistently the limiting variance.

Most of the existing literature on the distribution of R 2 deals with the case where the vector (Y, X (1) , . . . X (p) ) is normally distributed. In this context, [START_REF] Fisher | The general sampling distribution of the multiple correlation coefficient[END_REF] [START_REF] Fisher | The general sampling distribution of the multiple correlation coefficient[END_REF] obtained the first exact expression of the density function of R 2 , involving a Gauss hypergeometric series. Other exact expressions for the distribution and the cumulative distribution were given by Lee [START_REF] Lee | Some results on the sampling distribution of the multiple correlation coefficient[END_REF]. An exact expression of the cumulative distribution in terms of a series of Gamma distributions is given in Muirhead (1982, Theorem 5.2.5) [START_REF] Muirhead | Aspects of multivariate statistical theory[END_REF]. The first two principal terms of the latter development were given earlier by Lee [START_REF] Lee | Some results on the sampling distribution of the multiple correlation coefficient[END_REF]. Some quantile tables were derived by Lee [START_REF] Lee | Tables of upper percentage points of the multiple correlation coefficient[END_REF] from this second-order approximation of the cumulative distribution function. Lee's results were then implemented by Kelley [6] in his R MBESS package.

In the non-Gaussian case, there are relatively few results. Let us cite Muirhead [START_REF] Muirhead | Aspects of multivariate statistical theory[END_REF] Theorem 5.1.6, who gave the asymptotic normality of R 2 assuming that (Y, X (1) , . . . X (p) ) has an elliptical distribution. The paper closest to ours is that of Ogasawara [START_REF] Ogasawara | Asymptotic expansion and conditional robustness for the sample multiple correlation coefficient under nonnormality[END_REF], who gave an Edgeworth expansion for √ n( R 2 -R 2 ), under some conditions on the distributions of the variables, assuming in particular that all the variables Y, X (1) , . . . X (p) have a moment of order 8. We will compare our result to that of Ogasawara in more details in Remark 3.1 of Section 3. Let us also mention the results of [START_REF] Ogasawara | Asymptotic expansions of the distribution of the estimator for the generalized partial correlation under non-normality[END_REF] [START_REF] Ogasawara | Asymptotic expansions of the distribution of the estimator for the generalized partial correlation under non-normality[END_REF] for partial correlation.

The article is organised as follows: in Section 2, we recall White's results [START_REF] White | A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity[END_REF] concerning the least squares estimator of the coefficients of the orthogonal projection Proj V (Y ) = α 0 + α 1 X (1) + • • • + α p X (p) (consistency, asymptotic normality and estimation of the limiting covariance matrix). This serves two purposes: it recalls the robust procedure described by White [START_REF] White | A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity[END_REF] in the linear model (without the assumption of normality or homoscedasticity), and it gives an initial idea of the proofs that will enable us to obtain an asymptotic confidence interval for R 2 . In Section 3 we prove the consistency and asymptocic normality of R 2 . We also give an estimator of the limiting variance, which together with asymptotic normality provide an asymptotic confidence interval for R 2 . In Section 4, we give the asymptotic distribution of the joint law of the estimators of the individual R 2 's. In Section 5, we evaluate the performance of our confidence interval on different sets of simulations, and we compare this performance with that of the confidence interval that assumes that the vector (Y, X (1) , . . . X (p) ) is Gaussian. Finally, in Section 6, we give the extension of our procedure to the case of vector-valued random variables, and we briefly discuss the connections between our results and sensitivity analysis or screening methods.

Previous known results

Let V i be the L 2 -subspace generated by the variables (X

(0) i , X (1) i , . . . , X (p) i ). Let ε = Y -Proj V (Y ) and ε i = Y i -Proj V i (Y i ). We can then write Y i = Proj V i (Y i ) + ε i .
where the variables ε i are such that E(X (j) i ε i ) = 0 for any j ∈ {0, . . . , p}. Consider the following assumption: Assumption 2.1. Assume that 1, X (1) , . . . , X (p) are linearly independent (as a family of random variables in L 2 ).

Under Assumption 2.1 one can uniquely write

Y i = α 0 + α 1 X (1) i + • • • + α p X (p) i + ε i .
Let α 0:p = (α 0 , α 1 , . . . , α p ) t ∈ R p+1 , and X be the matrix whose jth column is

X (j) = (X (j) 1 , . . . , X (j) n ) t . Let Y = (Y 1 , . . . , Y n ) t and ε = (ε 1 , . . . , ε n ) t . One can write Y = Xα 0:p + ε.
Let M be the matrix such that M j, = E(X (j) X ( ) ). By Assumption (2.1) the matrix M is invertible and one can easily show that

α 0:p = M -1 (E(Y ), E(X (1) Y ) . . . , E(X (p) Y )) t .
2.1 Least square estimator of α 0:p

Let

M = 1 n X t X .
By the strong law of large numbers

lim n→∞ M = M almost surely. (2.1)
Under Assumption 2.1, M is invertible and for large enough n, X t X is also invertible, and one can define

α 0:p = ( α 0 , . . . , α p ) t = (X t X) -1 X t Y.
Let V p+1 be the sub-space of R n generated by the p + 1 columns X (0) , X (1) , . . . , X (p) of the matrix X, and Π V p+1 (Y) be the orthogonal projection of Y on V p+1 (with respect to the euclidean norm). We have

α 0 X (0) + α 1 X (1) + • • • + α p X (p) = Π V p+1 (Y),
and

α 0:p -α 0:p = ( α 0 , . . . , α p ) t -(α 0 , . . . , α p ) t = (X t X) -1 X t ε = M -1 1 n X t ε . (2.2)
2.2 Consistency and asymptotic normality of α 0:p Following White [START_REF] White | A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity[END_REF], we obtain the consistency and asymptotic normality of the least square estimators. The proofs of these results being simple and enlightening, we have chosen to recall them because they may make it easier to understand the proofs of Section 3. By definition of ε, for any j ∈ {0, . . . , p}, E(X (j) ε) = 0. Applying the strong law of large numbers to the jth coordinate of n -1 X t ε, we get

lim n→∞ 1 n X t ε j = lim n→∞ 1 n n i=1 X (j) i ε i = 0 almost surely. (2.5)
To prove the asymptotic normality, we need moment assumptions.

Assumption 2.2. Assume that E(X (j) 2 ε 2 ) < ∞ for all j = 0, . . . , p.

Proposition 2.2. Under Assumptions 2.1 and 2.2, let M(ε) be the matrix defined by

M(ε) j, = E(X (j) X ( ) ε 2 ) for any 0 ≤ j, ≤ p. (2.6)
Then √ n( α 0:p -α 0:p ) L -→ n→∞ N p+1 (0, M -1 M(ε)M -1 ). Remark 2.1. Assume that ε 2 -E(ε 2
) is orthogonal to the space S generated by the variables X (j) X ( ) for all 0 ≤ j ≤ ≤ p (which is true, for instance, if ε is independent of (X (1) , . . . , X (p) )). Then the limiting variance matrix in Proposition 2.2 writes

M -1 M(ε)M -1 = E(ε 2 )M -1 ,
and is the same as in the case where the conditional distribution of Y given (X (1) , . . . , X (p) ) is Gaussian with Var(Y |X (1) , . . . , X (p) ) = σ 2 (homoscedastic case; in that case E(ε 2 ) = σ 2 ). This simple observation enabled White [START_REF] White | A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity[END_REF] to formulate his famous homoscedasticity test, which in fact consists of testing whether

ε 2 -E(ε 2 ) is orthogonal to S (H 0 hypothesis) or not.
Proof of Proposition 2.2. Starting from (2.2) and using (2.3), it suffices to prove that

1 √ n X t ε L -→ n→∞ N p+1 (0, M(ε)). (2.7)
Now (2.7) follows from a direct application of the central limit theorem in R p+1 , since for any 0 ≤ j ≤ p, E(X (j) ε) = 0 and since Assumption 2.2 holds.

Estimation of the limiting covariance matrix

Under Assumption 2.2, by the strong law of large numbers, for any j, ∈ {0, . . . , p},

lim n→∞ 1 n n i=1 X (j) i X ( ) i ε 2 i = M(ε) j, almost surely. (2.8)
Of course, the quantity on right hand is not an estimator of M(ε) j, since the ε i 's are not observed. White [START_REF] White | A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity[END_REF] proposed then to replace the ε i 's by the residuals.

ε i = Y i -(X α 0:p ) i = Y i -α 0 X (0) i -• • • -α p X (p) i . Let then M(ε) j, = 1 n n i=1 X (j) i X ( ) i ε 2 i .
To prove the consistency of M(ε) j, , we shall assume that Assumption 2.2 holds and that:

Assumption 2.3. E(X (j) 4 ) < ∞ for all j ∈ {0, . . . , p}.

Proposition 2.3. Under Assumptions 2.2 and 2.3, for any j, ∈ {0, . . . , p},

lim n→∞ M(ε) j, = M(ε) j, almost surely.
Consequently, under Assumptions 2.1-2.3,

lim n→∞ M -1 M(ε) M -1 = M -1 M(ε)M -1 almost surely.
Remark 2.2. The Assumptions 2.2 and 2.3 are equivalent to: E(X (j) 2 Y 2 ) < ∞ and E(X (j) 4 ) < ∞ for any j ∈ {0, . . . , p}.

Remark 2.3. As usual, Propositions 2.2 and 2.3 enable to obtain confidence regions or test procedures for the coefficients α i . For instance, if one wants to test

H 0 : α i 1 = • • • = α i k = 0 for some 0 ≤ i 1 < i 2 < • • • < i k ≤ p
one can proceed as follows. Let C k be the matrix with k rows and p columns, such that all coordinates of the jth row are 0 except the coordinate i j which is equal to 1. Then, by Proposition 2.2,

√ n(( α i 1 , . . . , α i k ) t -(α i 1 , . . . , α i k ) t ) L -→ n→∞ N k (0, C k M -1 M(ε)M -1 C t k ).
Now, under Assumptions 2.1-2.3, by applying Proposition 2.3

V k = C k M -1 M(ε) M -1 C t k converges almost surely to C k M -1 M(ε)M -1 C t k . If M(ε) is invertible, then the matrix C k M -1 M(ε)M -1 C t k is also invertible and V -1/2 k is well defined (for n large enough). Let then (ξ 1 , . . . , ξ k ) t = √ n V -1/2 k ( α i 1 , . . . , α i k ) t . Under H 0 , we see that the test statistic k i=1 ξ 2 i is such that k i=1 ξ 2 i L -→ n→∞ χ 2 (k) . Proof of Proposition 2.3. Write ε 2 i = ε 2 i + R i , with R i = 2ε i ( ε i -ε i ) + ( ε i -ε i ) 2 .
The following upper bound holds

|R i | ≤ 2 p j=0 |( α j -α j )X (j) i ε i | + (p + 1) p j=0 ( α j -α j ) 2 X (j) i 2 .
Tacking into account (2.8), to prove Proposition 2.3, it suffices to show that, for 1 ≤ j, k ≤ p and 1 ≤ i ≤ n,

lim n→∞ 1 n n i=1 |X (j) i X (k) i R i | = 0 almost surely.
(2.9)

Since | α i -α i | converges almost surely to 0, (2.9) will be satisfied provided, for 1 ≤ j, k, ≤ p,

lim sup n→∞ 1 n n i=1 |X (j) i X (k) i X ( ) i ε i | < ∞ and lim sup n→∞ 1 n n i=1 |X (j) i X (k) i X ( ) i 2 | < ∞ almost surely,
which is true by applying once again the strong law of large numbers under Assumptions 2.2 and 2.3.

3 Confidence interval for R 2

Assuming that Var(Y ) > 0, the mutliple correlation coefficient R 2 is defined by:

R 2 = Var(Proj V (Y )) Var(Y )
.

For ε i = Y i -Proj V i (Y i )
, recall that, under Assumption 2.1,

Y i = Proj V i (Y i ) + ε i = α 0 1 + α 1 X (1) i + • • • + α p X (p) i + ε i .
For any vector z = (z 1 , . . . , z n ) t in R n , the empirical variance of z is given by

Var n (z) = 1 n n i=1 (z i -z) 2 .
Let 1 n = (1, . . . , 1) t and let Π E be the orthogonal projection (with respect to the the euclidean norm) on a subspace E of R n . Let V p+1 be the subspace of R n generated by (1 n , X (1) , . . . ,

X (p) ). The empirical estimator R 2 of R 2 is then R 2 = Var n (Π V p+1 (Y)) Var n (Y) . (3.1)
Our first goal is to give a simple expression of this estimator. Let α 1:p = (α 1 , . . . , α p ) t and α 1:p = ( α 1 , . . . , α p ) t ; for the sake of simplicity, we shall omit the indexes and write α = α 1:p and α = α 1:p . Let also θ = (θ 1 , . . . , θ p ) t and θ = ( θ 1 , . . . , θ p ) t , where

θ k = Cov(Y, X (k) ) Var(Y ) and θ k = n i=1 (Y i -Y)(X (k) i -X (k) ) n i=1 (Y i -Y) 2 . (3.2)
The following lemma gives the expression of R 2 and R 2 in terms of θ, α, θ and α.

Lemma 3.1. Under Assumption 2.1, R 2 = θ t α and R 2 = θ t α.
Proof of Lemma 3.1. We prove the second point, the first point being proved in the same way. Write

Π V p+1 (Y) -Π V p+1 (Y)1 n = α 1 (X (1) -X (1) 1 n ) + • • • + α p (X (p) -X (p) 1 n ) . Now, if S p is the subspace of R n generated by ((X (1) -X (1) 1 n )), . . . , (X (p) -X (p) 1 n ))) it is easy to check that α 1 (X (1) -X (1) 1 n )) + • • • + α p (X (p) -X (p) 1 n )) = Π Sp (Y -Y1 n ).
Since by Assumption 2.1 the variables 1, X (1) , . . . , X (p) are linearly independent in L 2 , the variance matrix M 0 of (X (1) , . . . , X (p) ) t defined by

(M 0 ) j, = Cov(X (j) , X ( ) ) j, ∈ {1, . . . , p} (3.3)
is invertible. Denote by X 0 the matrix defined by

X 0 = ((X (j) i -X (j) )) 1≤i≤n,1≤j≤p . (3.4)
By the strong law of large numbers

lim n→∞ 1 n X t 0 X 0 = M 0 almost surely. Since M 0 defined in (3.3) is invertible, it follows that X t 0 X 0 is also invertible for n large enough, in such a way that α = (X t 0 X 0 ) -1 X t 0 (Y -Y1 n ), (3.5)
and

Π Sp ((Y -Y1 n )) = X 0 (X t 0 X 0 ) -1 X t 0 (Y -Y1 n )
. Let • eu be the euclidean norm on R n . The following equality hold:

nVar n (Π V p+1 (Y)) = Π Sp (Y -Y1 n ) 2 eu = (Y -Y1 n ) t X 0 (X t 0 X 0 ) -1 X t 0 X 0 (X t 0 X 0 ) -1 X t 0 (Y -Y1 n ) = (Y -Y1 n ) t X 0 (X t 0 X 0 ) -1 X t 0 (Y -Y1 n ) = (Y -Y1 n ) t X 0 α . (3.6)
The result follows from (3.1), (3.6) and (3.7) since

(Y -Y1 n ) t X 0 nVar n (Y) = θ t . (3.7)
Now, to build a confidence interval for R 2 , we need to describe the asymptotic distribution of

√ n [ α : θ] -[α : θ] := √ n ( α 1 , . . . , α p , θ 1 , . . . , θ p ) t -(α 1 , . . . , α p , θ 1 , . . . , θ p ) t .
Let us then define the two matrices involved in this asymptotic distribution. For k ∈ {1, . . . , p}, let

e (k) = X (k) -Proj W (X (k) ) and e (k) i = X (k) i -Proj W i (X (k) i ), (3.8)
where W is the subspace of L 2 generated by 1 and Y , and W i is the subspace of L 2 generated by 1 and Y i . Let also e (k) = (e

(k) 1 , . . . , e (k)
n ) t . Assumption 3.1. E((X (j) -E(X (j) )) 2 ε 2 ) < ∞ and E((Y -E(Y )) 2 e (j)
2 ) < ∞ for all j ∈ {1, . . . , p}.

Let then A be the 2p × 2p symmetric matrix defined as follows:

• If (j, k) ∈ {1, . . . , p} 2 then A j,k = E((X (j) -E(X (j) ))(X (k) -E(X (k) ))ε 2 ); • If (j, k) ∈ {p + 1, . . . , 2p} 2 then A j,k = E((Y -E(Y )) 2 e (j-p) e (k-p) ); • (j, k) ∈ {1, . . . , p} × {p + 1, . . . , 2p} then A j,k = E((X (j) -E(X (j) ))(Y -E(Y ))εe (k-p) ).
Let δ j,k = 0 if j = k and δ j,j = 1 and let B be the 2p × 2p symetric matrix defined as follows:

• If (j, k) ∈ {1, . . . , p} 2 then B j,k = (M -1 0 ) j,k (see (3.3) for the definition of M 0 );

• If (j, k) ∈ {p + 1, . . . , 2p} 2 then B j,k = (Var(Y )) -1 δ j,k ;
• If (j, k) ∈ {1, . . . , p} × {p + 1, . . . , 2p} then B j,k = 0.

Proposition 3.1. Under Assumptions 2.1 and 3.1

√ n [ α : θ] -[α : θ] L -→ n→∞ N 2p (0, BAB).
Proof of Proposition 3.1. Let B be the 2p × 2p symetric matrix defined as follows:

• If (j, k) ∈ {1, . . . , p} 2 then B j,k = n(X t 0 X 0 ) -1 j,k ; • If (j, k) ∈ {p + 1, . . . , 2p} 2 then B j,k = (Var n (Y)) -1 δ j,k ; • If (j, k) ∈ {1, . . . , p} × {p + 1, . . . , 2p} then B j,k = 0.
Starting from (3.2) and (3.5), and noting that

Y i -Y = α 1 (X (1) i -X (1) ) + • • • + α p (X (p) i -X (p) ) + (ε i -ε) , (X (k) i -X (k) ) = θ k (Y i -Y) + (e (k)
i -e (k) ) for k ∈ {1, . . . , p} , we see that

√ n [ α : θ] -[α : θ] = B 1 √ n (ε t X 0 , (Y -Y1 n ) t E) t , (3.9)
where E is the n × p matrix such that, for (i, j) ∈ {1, . . . , n} × {1, . . . p}, (E i,j ) = e (j)

i . Denote by X 0 the matrix defined by X0 = (X (j)

i -E(X (j) )) 1≤i≤n,1≤j≤p .

(3.10)

One can easily check that

1 √ n (ε t X 0 , (Y 1 -Y, . . . , Y n -Y)E) t -(ε t X0 , (Y 1 -E(Y ), . . . , Y n -E(Y ))E) t P -→ n→∞ 0 . (3.11)
Now, by the multivariate central limit theorem,

1 √ n (ε t X0 , (Y 1 -E(Y ), . . . , Y n -E(Y ))E) t L -→ n→∞ N 2p (0, A). (3.12)
The result follows from (3.9), (3.11) and (3.12), since by the strong law of large numbers,

lim n→∞ B = B almost surely. (3.13)
As a consequence, we get Proposition 3.2. Under Assumptions 2.1 and 3.1

√ n R 2 -R 2 L -→ n→∞ N (0, V ).
where V = (θ 1 , . . . , θ p , α 1 , . . . , α p )BAB(θ 1 , . . . , θ p , α 1 , . . . , α p ) t .

(3.14)

Remark 3.1. Ogasawara [START_REF] Ogasawara | Asymptotic expansion and conditional robustness for the sample multiple correlation coefficient under nonnormality[END_REF] has given an Edgeworth expansion of

√ n( R 2 -R 2 )
, which is a more precise result than Proposition 3.2, but requires some conditions on the distribution of (Y, X (1) , • • • , X (p) ). In particular, he required that the all the variables have a moment of order 8. Note that, if we do not look for an Edgeworth expansion, the method of Ogasawara consists in expressing R 2 as a differentiable function of S = ((Cov(X (i) , X (j) )) 1≤i≤j≤p , (Cov(Y, X (i) )) 1≤i≤p , Var(Y )), proving the CLT for √ n( S -S) ( S being the empirical estimator of S), and applying the delta method. The differences with our approach are the following :

• To prove the CLT for √ n( S -S), one needs moments of order 4 for all variables, which is a more restrictive condition than the moment conditions of our Proposition 3.2.

• Applying the delta method to a function of S implies that the limiting variance V will be expressed as a function of a q × q matrix, where q = (p + 1)(p + 2)/2 (the matrix Ω in [START_REF] Ogasawara | Asymptotic expansion and conditional robustness for the sample multiple correlation coefficient under nonnormality[END_REF]), while we can express V as a function of a 2p × 2p matrix (the matrix BAB, see (3.14)). Note that, with our expression of V we are able to give simple sufficient conditions ensuring that V > 0 (see Lemma 3.2 below).

• Finally our approach can be easily extended to the the case where the variables Y, X (1) , • • • , X (p) are vector-valued (see Section 6.1).

Proof of Proposition 3.2. Starting from Lemma 3.1 and Proposition 3.1, it suffices to apply the delta-method to the function φ : R p × R p → R defined by φ(x, y) = x t y .

The proof Proposition 3.2 is complete by evaluating the differential Dφ x,y of φ at point (x, y):

Dφ x,y (h 1 , h 2 ) = y t h 1 + x t h 2 = (y, x) t (h 1 , h 2 ) .
To build a confidence interval for R 2 , it remains to find a consistent estimator of V . We shall simply replace each element in the definition of V by its empirical counterpart.

Let then A be the 2p × 2p symmetric matrix defined as follows:

• If (j, k) ∈ {1, . . . , p} 2 then A j,k = 1 n n i=1 (X (j) i -X (j) )(X (k) i -X (k) ) ε 2 i ; • If (j, k) ∈ {p + 1, . . . , 2p} 2 then A j,k = 1 n n i=1 (Y i -Y) 2 e (j-p) i e (k-p) i
, where e

(k-p) i = (X (k-p) i -X (k-p) ) -θ k-p (Y i -Y);
• If (j, k) ∈ {1, . . . , p} × {p + 1, . . . , 2p} then Hence, it remains to prove that A converges to A almost surely.

A j,k = 1 n n i=1 (X (j) i -X (j) )(Y i -Y) ε i e (k-p) i . Proposition 3.3. Let V n = ( θ 1 , . . . ,
To prove this point, we first introduce the matrix A defined as follows:

• If (j, k) ∈ {1, . . . , p} 2 then A j,k = 1 n n i=1 (X (j) i -X (j) )(X (k) i -X (k) )ε 2 i ; • If (j, k) ∈ {p + 1, . . . , 2p} 2 then A j,k = 1 n n i=1 (Y i -Y) 2 e (j-p) i e (k-p) i ; • If (j, k) ∈ {1, . . . , p} × {p + 1, . . . , 2p} then A j,k = 1 n n i=1 (X (j) i -X (j) )(Y i -Y)ε i e (k-p) i
.

Since E(Y 4 ) < ∞ and E(X (j) 4 ) < ∞, we deduce from the strong law of large numbers that

lim n→∞ A = A almost surely.
Let us now prove that, for (j, k) ∈ {1, . . . , p} 2 ,

lim n→∞ | A j,k -A j,k | = 0 almost surely.
We write

ε 2 i = ε 2 i + R i , with R i = 2ε i ( ε i -ε i ) + ( ε i -ε i ) 2 .
One can easily see that

|R i | ≤ 2 p =0 |( α -α )X ( ) i ε i | + (p + 1) p =0 ( α -α ) 2 X ( ) i 2 . (3.16) Moreover | A j,k -A j,k | ≤ 1 n n i=1 |(X (j) i -X (j) )(X (k) i -X (k) )R i | . (3.17)
Combining (3.16) and (3.17), and using that | α j -α j | converges almost surely to 0, we infer that A j,k -A j,k converges almost surely to 0 as soon as lim sup

n→∞ 1 n n i=1 |X (j) i X (k) i X ( ) i ε i | < ∞ and lim sup n→∞ 1 n n i=1 |X (j) i X (k) i X ( ) i 2 | < ∞ almost surely,
which is true (applying once again the strong law of large numbers) since E(Y 4 ) < ∞ and E(X (j) 4 ) < ∞ for any j ∈ {1, . . . , p}.

Let us now prove that, for (j, k) ∈ {p+1, . . . , 2p} 2 , A j,k -A j,k converges almost surely to 0 as n → ∞. We write e

(j-p) i e (k-p) i = e (j-p) i e (k-p) i + T i , with |T i | ≤ |e (j-p) i ( e (k-p) i -e (k-p) i )| + |e (k-p) i ( e (j-p) i -e (j-p) i )| + 1 2 ( e (k-p) i -e (k-p) i ) 2 + 1 2 ( e (j-p) i -e (j-p) i ) 2 . (3.18) For ∈ {1, . . . , p}, let β = E(X ) -θ E(Y ) and β = X ( ) -θ Y. With these notations X ( ) k = β + θ Y k + e ( ) k
and X ( )

k = β + θ Y k + e ( ) k . ( 3
|T i | ≤ |e (j-p) i ( β k-p -β k-p )| + |e (j-p) i Y i ( θ k-p -θ k-p )| + |e (k-p) i ( β j-p -β j-p )| + |e (k-p) i Y i ( θ j-p -θ j-p )| + ( β k-p -β k-p ) 2 + ( β j-p -β j-p ) 2 + ( θ k-p -θ k-p ) 2 Y 2 i + ( θ j-p -θ j-p ) 2 Y 2 i . (3.20) Moreover | A j,k -A j,k | ≤ 1 n n i=1 |(Y i -Y) 2 T i | . ( 3
|Y 2 i e (j-p) i | < ∞, lim sup n→∞ 1 n n i=1 |Y 3 i e (j-p) i | < ∞, and lim sup n→∞ 1 n n i=1 Y 4 i < ∞,
which is true (applying once again the strong law of large numbers) since E(Y 4 ) < ∞ and E(X (j) 4 ) < ∞ for any j ∈ {1, . . . , p}.

The fact that, for (j, k) ∈ {1, . . . , p} × {p + 1, . . . , 2p}, A j,k -A j,k converges almost surely to 0 as n → ∞ may be proved exactly in the same way.

As an immediate Corollary, we get Corollary 3.1. Assume that Assumptions 2.1 and 2.3 hold, and that E(Y 4 ) < ∞. Assume also that V > 0, where V is defined by (3.14). For δ ∈ (0, 1), let c 1-(δ/2) be the quantile of order 1 -δ of the N (0, 1)-distribution, and let V n be defined by (3.15).

Then   R 2 - c 1-(δ/2) V n √ n , R 2 + c 1-(δ/2) V n √ n   is an asymptotic confidence interval for R 2 of level 1 -δ. Remark 3.2.
As is sensitivity analysis, one can also define

R 2 T i = 1 - Var(Proj V (-i) (Y ) ) Var(Y ) ,
where V (-i) is the linear space of L 2 generated by the variables 1, X (1) , . . . , X (i-1) , X (i+1) , . . . , X (p) . Then R 2 T i is the proportion of the variance of Y that is not explained by the best linear predictor based on the variables 1, X (1) , . . . , X (i-1) , X (i+1) , . . . , X (p) . Denoting by R 2 (-i) the proportion of the variance of Y that is explained by the best linear predictor based on 1, X (1) , . . . , X (i-1) , X (i+1) , . . . , X (p) , we immediately see that R 2

T i = 1 -R 2 (-i) .
Consequently, a confidence interval for R 2 T i is directly obtained from a confidence interval for R 2 (-i) .

To be complete, we discuss the positivity of the variance term V .

Lemma 3.2. Assume that M 0 is invertible, that E((X (j) -E(X (j) ) 2 ε 2 ) < ∞ and Assumption 3.1 holds.

1. If R 2 = 1 then R 2 = 1 almost surely, and V = 0.

2. If R 2 = 0 then V = 0, and n R 2 converges in distribution to (G 1 , . . . , G p )(G p+1 , . . . , G 2p ) t , where (G 1 , . . . , G 2p
) is the Gaussian random vector with covariance matrix BAB.

3. If R 2 > 0 and if the family

((Y -E(Y )) 2 , ((Y -E(Y ))(X (j) -E(X (j) ))) 1≤j≤p , ((X (j) -E(X (j) ))(X (j) -E(X (j) ))) 1≤i≤j≤p ) (3.22)
is linearly independent, then V > 0.

Remark 3.3. The fact that the family (3.22) is linearly independent implies that (Y, X (1) , . . . , X (p) ) is also linearly independent, so that R 2 < 1.

Remark 3.4. Note that one can easily check that the variable (G 1 , . . . , G p )(G p+1 , . . . , G 2p ) t of Item 2 is non negative. Indeed, in that particular case, (G 1 , . . . , G p ) t = M -1 0 (Z 1 , . . . , Z p ) t and (G p+1 , . . . , G 2p ) t = (Var(Y )) -1 (Z 1 , . . . , Z p ) t , where (Z 1 , . . . , Z p ) t is a Gaussian vector with covariance matrix (A i,j ) 1≤i,j≤p . The asymptotic distribution of Item 2 may be used to test H 0 : R 2 = 0 against H 1 : R 2 = 0, but it is simpler to test the equivalent hypothesis H 0 : α 1 = • • • = α p = 0 against H 1 : α i = 0 for some i ∈ {1, . . . , p} (see Remark 2.3).

Proof of Lemma 3.2. Item 1 is clear

: R 2 = 1 if and only if Y = α 0 + α 1 X (1) + • • • + α p X (p)
almost surely. In that case R 2 = 1 almost surely. Consequently R 2 -R 2 = 0 almost surely, and V = 0. We now prove Item 2. R 2 = 0 if and only if Y -E(Y ) is orthogonal to the space generated by (X (1) -E(X (1) ), . . . , (X (p) -E(X (p) )). This is equivalent to α 1 = • • • = α p = 0, which is also equivalent to θ 1 = • • • = θ p = 0. This implies that V = 0. Applying Proposition 3.1, we see that

√ n( α 1 , . . . , α p , θ 1 , . . . , θ p ) t L -→ n→∞ N 2p (0, B A B).
The last assertion of Item 2 follows by noting that

n R 2 = √ n( θ 2 , . . . , θ p ) √ n( α 2 , . . . , α p ) t .
We now prove Item 3. Since R 2 > 0, it follows that the two vectors (α 1 , . . . , α p ) t and (θ 1 , . . . , θ p ) t are not equal to (0, . . . , 0) t . Let (a 1 , . . . , a p , b 1 , . . . , b p ) t = B(θ 1 , . . . , θ p , α 1 , . . . , α p ) t .

By definition of the matrix B, and since M 0 is invertible, we infer that the two vectors (a 1 , . . . , a p ) t and (b 1 , . . . , b p ) t are not equal to (0, . . . , 0) t . Now, by definition, V = (a 1 , . . . , a p , b 1 , . . . , b p ) A(a 1 , . . . , a p , b 1 , . . . , b p ) t . Since A is the covariance matrix of the vector (X (1) -E(X (1) ))ε, . . . , (X (p) -E(X (p) ))ε, (Y -E(Y ))e (1) , . . . , (Y -E(Y ))e (p) t ,

we can write

V = Var a 1 (X (1) -E(X (1) ))ε + • • • + a p (X (p) -E(X (p) ))ε + b 1 (Y -E(Y ))e (1) + • • • + b p (Y -E(Y ))e (p)
.

Hence V = 0 if and only if

a 1 (X (1) -E(X (1) ))ε+• • •+a p (X (p) -E(X (p) ))ε = -(b 1 (Y -E(Y ))e (1) +• • •+b p (Y -E(Y ))e (p) ) almost surely. (3.23) Recall now that ε = (Y -E(Y )) -α 1 (X (1) -E(X (1) )) -• • • -α p (X (p) -E(X (p) )) e (k) = (X (k) -E(X (k) )) -θ k (Y -E(Y )) for k ∈ {1, . . . , p}.
Without loss of generality, assume that α 1 = 0. Then, on the left side of (3.23) the terms

-a 1 α 1 (X (1) -E(X (1) )) 2 , -a 2 α 1 (X (1) -E(X (1) ))(X (2) -E(X (2) )),
. . . , -a p α 1 (X (1) -E(X (1) ))(X (p) -E(X (p) ))

appear, but they do not appear in the right side of (3.23). Since we assumed that the family (3.22) is linearly independent, we infer that V = 0 implies a 1 = • • • = a p = 0. Since we know that (a 1 , . . . , a p ) t is not equal to (0, . . . , 0) t , we conclude that V > 0.

Asymptotic joint distribution of individual R 2 's

For k ∈ {1, . . . , p}, let U (k) be the L 2 -subspace generated by 1 and X (k) , and let (assuming that Var(Y ) > 0)

R 2 (k) = Var(Proj U (k) (Y )) Var(Y ) . From Lemma 3.1, we have R 2 (k) = τ k θ k
, where θ k is defined by (3.2) and

τ k = Cov(Y, X (k) ) Var(X (k) ) . (4.1)
The empirical estimator of R 2 (k) is then

R 2 (k) = τ k θ k , where τk = n i=1 (Y i -Y)(X (k) i -X (k) ) n i=1 (X (k) i -X (k) ) 2 . (4.2)
The asymptotic distribution of

√ n( R 2 (k) -R 2 (k)
) is given in Proposition 3.2. In this section, we give the asymptotic joint distribution of

√ n ( R 2 (1) . . . , R 2 (p) ) t -(R 2 (1) , . . . , R 2 (p) ) t .
Let τ = (τ 1 , . . . , τ p ) t and τ = ( τ 1 , . . . , τ p ) t . The first step is to identify the limit distribution of

√ n [ τ : θ] -[τ : θ] .
Let us then define the two matrices involved in this asymptotic distribution. For k ∈ {1, . . . , p}, let

ε (k) = Y -Proj U (k) (Y ) and ε (k) i = Y i -Proj U (k) i (Y i ),
where

U (k) i
is the subspace of L 2 generated by 1 and

X (k) i . Let also ε (k) = (ε (k) 1 , . . . , ε (k) n ) t .
Recall that e (k) and e (k) i have been defined in (3.8).

Assumption 4.1. E((X (j) -E(X (j) )) 2 ε (j) 2 ) < ∞ and E((Y -E(Y )) 2 e (j)
2 ) < ∞ for all j ∈ {1, . . . , p}.

Let then C be the 2p × 2p symmetric matrix defined as follows:

• If (j, k) ∈ {1, . . . , p} 2 then C j,k = E((X (j) -E(X (j) ))(X (k) -E(X (k) ))ε (j) ε (k) ); (4.3) • If (j, k) ∈ {p + 1, . . . , 2p} 2 then C j,k = E((Y -E(Y )
) 2 e (j-p) e (k-p) );

• (j, k) ∈ {1, . . . , p} × {p + 1, . . . , 2p} then

C j,k = E((X (j) -E(X (j) ))(Y -E(Y ))ε (j) e (k-p )).
Let D be the 2p × 2p diagonal matrix defined as follows :

• If j ∈ {1, . . . , p} then D j,j = (Var(X (j) )) -1 .

• If j ∈ {p + 1, . . . , 2p} then D j,j = (Var(Y )) -1 .
Proposition 4.1. Under Assumption 4.1 and assuming that Var(X (j) ) > 0 for all j ∈ {1, . . . , p}, we have

√ n [ τ : θ] -[τ : θ] L -→ n→∞ N 2p (0, DCD).
Proof of Proposition 4.1. Let D be the 2p × 2p diagonal matrix defined as follows :

• If j ∈ {1, . . . , p} then D j,j = (Var n (X (j) )) -1 .

• If j ∈ {p + 1, . . . , 2p} then D j,j = (Var n (Y)) -1 .

Starting from (3.2) and (3.5), and noting that, for k ∈ {1, . . . , p},

Y i -Y = τ k (X (k) i -X (k) ) + (ε (k) i -ε (k) ) , (X (k) i -X (k) ) = θ k (Y i -Y) + (e (k)
i -e (k) ) ,

we see that (recall that X 0 and X0 have been defined in (3.4) and (3.10))

√ n [ τ : θ] -[τ : θ] = D 1 √ n (ε (1) t X (1) 0 , . . . , ε (p) t X (p) 0 , (Y -Y1 n ) t e (1) , . . . , (Y -Y1 n ) t e (p) ) t . (4.4)
One can easily check that

1 √ n (ε (1) t X (1) 0 , . . . , ε (p) t X (p) 0 , (Y -Y1 n ) t e (1) , . . . , (Y -Y1 n ) t e (p) ) t -(ε (1) t X(1) 0 , . . . , ε (p) t X(p) 0 , (Y -E(Y )1 n ) t e (1) , . . . , (Y -E(Y )1 n ) t e (p) ) t P -→ n→∞ 0. (4.5)
Now, by the multivariate central limit theorem, The result follows from (4.4), (4.5), (4.6) and (4.7).

1 √ n (ε (1) t X(1) 0 , . . . , ε (p) t X(p) 0 , (Y -E(Y )1 n ) t e (1) , . . . , (Y -E(Y )1 n ) t e (p) ) t L -→ n→∞ N 2p (0, C). ( 4
Let H be the matrix with p rows and 2p column such that

• If (i, j) ∈ {1, . . . , p} 2 then H i,j = θ i δ i,i .

• If i ∈ {1, . . . , p} and j ∈ {p + 1, . . . , 2p} then H i,j = τ i δ i,p+i .

As a consequence of Proposition 4.1, we have Proposition 4.2. Under Assumption 4.1 and assuming that Var(X (j) ) > 0 for all j ∈ {1, . . . , p}, we have

√ n ( R 2 (1) . . . , R 2 (p) ) t -(R 2 (1) , . . . , R 2 (p) ) t L -→ n→∞ N (0, HDCDH t ).
Remark 4.1. Other quantities of interest are

S 2 (k) = Var(Proj U (k) (Y )) Var(Proj V (Y )) for k ∈ {1, . . . , p},
which are the proportions of the variance of Proj V (Y ) that are explained by the best linear predictors in L 2 based on the variables 1, X (k) . Note that

S 2 (k) = R 2 (k) R 2 = τ k θ k θ t α .
Then, the asymptotic distribution of √ n ( S 2 (1) . . . , S

) t = φ([ τ : θ]) -φ([τ : θ]) ,
it suffices to apply the delta method to the function φ. The results follows from Proposition 4.1 and the fact that the differential Dφ τ ,θ of φ at point (τ , θ) is given by Dφ τ ,θ (h) = Hh .

To conclude this section, let us a give a consistent estimator of the matrix HDCDH t . We shall simply replace each matrix by its empirical counterpart.

Let then C be the 2p × 2p symmetric matrix defined as follows:

• If (i, j) ∈ {1, . . . , p} 2 then C j,k = 1 n n i=1 (X (j) i -X (j) )(X (k) i -X (k) ) ε (j) i ε (k) i , where ε (k)
i = (Y i -Y) -τ k (X (k) i -X (k) ); • If (j, k) ∈ {p + 1, . . . , 2p} 2 then C j,k = 1 n n i=1 (Y i -Y) 2 e (j-p) i e (k-p) i ;
• If (j, k) ∈ {1, . . . , p} × {p + 1, . . . , 2p} then

C j,k = 1 n n i=1 (X (j) i -X (j) )(Y i -Y) ε (j) i e (k-p) i . Proposition 4.3. Assume that Var(X (j) ) > 0, that E(Y 4 ) < ∞ and that E(X (j) 4 ) < ∞ for any j ∈ {1, . . . , p}. Then lim n→∞ H D C D H t = HDCDH t almost surely.
Proof of Proposition 4.3. By the strong law of large numbers, H and D converge almost surely to H and D respectively. Hence, it remains to prove that C converges to C almost surely, which can be done as in the proof of Proposition 3.3.

Simulations

We shall consider different models. For each model, we shall estimate the coverage level of two confidence intervals of level 95% for the R 2 , for n = 200 to n = 1000, via a basic Monte-Carlo procedure (with N = 3000 repetitions).

n 200 300 400 500 600 700 800 900 1000 CI1 0.942 0.945 0.952 0.949 0.951 0.947 0.947 0.952 0.95 CI2 0.949 0.952 0.954 0.95 0.953 0.951 0.949 0.951 0.951 Table 1: Estimated coverage levels of CI1 and CI2 at level 95% for model (5.1) with ε ∼ N (0, 1) and Proj V (Y ) = 0.5 + 0.5X (1) + X (2) . n 200 300 400 500 600 700 800 900 1000 CI1 0.934 0.937 0.941 0.942 0.947 0.949 0.95 0.945 0.95 CI2 0.936 0.932 0.933 0.936 0.937 0.936 0.939 0.935 0.934 Table 2: Estimated coverage levels of CI1 and CI2 at level 95% for model (5.1) with ε ∼ St [START_REF] Muirhead | Aspects of multivariate statistical theory[END_REF] and Proj V (Y ) = 0.5 + 0.5X (1) + X (2) .

The first confidence interval (CI1) is the one described in Corollary 3.1 (with the small change that we take the quantile t n,0.975 of the Student distribution St(n) instead of the quantile c 0.975 ; this has no theoretical justification, but can improve a bit the coverage level for small n). The second confidence interval (CI2) is the non-asymptotic confidence interval ci.R2 of the R package MBESS, as described by Kelley [START_REF] Kelley | Confidence Intervals for Standardized Effect Sizes: Theory, Application, and Implementation[END_REF], and based on a precise approximation of the distribution of R 2 given by Lee [START_REF] Lee | Some results on the sampling distribution of the multiple correlation coefficient[END_REF], [START_REF] Lee | Tables of upper percentage points of the multiple correlation coefficient[END_REF], when the vector (Y, X (1) , . . . , X (p) ) is Gaussian. It is clearly indicated in Kelley [START_REF] Kelley | Sample Size Planning for the Squared Multiple Correlation Coefficient: Accuracy in Parameter Estimation via Narrow Confidence Intervals[END_REF] (Discussion, page 552-553) that CI2 is not robust to the non-normality of (Y, X (1) , . . . , X (p) ), which will be confirmed by the simulations.

5.1

An example where (Y, X (1) , X (2) ) is Gaussian

We consider here the model Y = 0.5 + 0.5X (1)

+ X (2) + ε , (5.1)
where X 1 , X 2 , ε are i.i.d. with N (0, 1) distribution. Let V be the sub-space of L 2 generated by 1, X (1) , X (2) , so that Proj V (Y ) = 0.5 + 0.5X (1) + X (2) . One can then easily check that R 2 = 5/9. The estimated coverage levels for CI1 and CI2 are given in Table 1 below. We see that the estimated coverage level of CI2 is always close to 0.95, which is not a surprise since CI2 is based on a precise approximation of the distribution of R 2 in the case where (Y, X (1) , X (2) ) is Gaussian. We see that CI1 also gives very good results, with estimated coverage levels between 0.945 and 0.952 as soon as n ≥ 300.

An example where the error term has a Student distribution

We consider here a slight modification of the model (5.1), where ε ∼ St [START_REF] Muirhead | Aspects of multivariate statistical theory[END_REF]. Again, V is the sub-space of L 2 generated by 1, X (1) , X (2) , so that Proj V (Y ) = 0.5 + 0.5X (1) + X (2) . One can then easily check that R 2 = 0.5.

The estimated coverage levels for CI1 and CI2 are given in Table 2 below. We see that the estimated coverage level of CI2 is always between 0.932 and 0.94, and is not getting closer to 0.95 as n increases. This confirms that CI2 is not robust to non normality, even when the distribution of the error is symmetric. We see that CI1 is always better than CI2 as soon as n ≥ 300, with a coverage level between 0.94 and 0.95 when n ≥ 400. n 200 300 400 500 600 700 800 900 1000 CI1 0.91 0.926 0.928 0.935 0.937 0.939 0.944 0.944 0.946 CI2 0.878 0.878 0.878 0.877 0.877 0.887 0.884 0.885 0.887 Table 3: Estimated coverage levels of CI1 and CI2 at level 95% for the heteroscedastic model (5.2) with Proj V (Y ) = 0.5 + 0.5X (1) + X (2) .

n

An heteroscedastic example

We continue with a modification of the model (5.1). We consider the model Y = 0.5 + 0.5X (1)

+ X (2) + 0.2 + 0.8X (1) 2 e , (5.2)
where X (1) , X (2) , e are i.i.d. with N (0, 1) distribution. Let V be the sub-space of L 2 generated by 1, X (1) , X (2) , so that Proj V (Y ) = 0.5 + 0.5X (1) + X (2) . Again, one can easily check that R 2 = 5/9. The estimated coverage levels for CI1 and CI2 are given in Table 3 below. We see that the estimated coverage level of CI2 is always around 0.88. Again, this confirms that CI2 is not robust to non normality. We see that CI1 is always better than CI2, with a coverage level greater than 0.935 for n ≥ 500, and greater than 0.94 for n ≥ 800.

An example where the model is misspecified

We consider the model

Y = X 2 + e , (5.3)
where X ∼ N (0, 1), e ∼ N (0, 1) and e is independent of X. Let V be the sub-space of L 2 generated by 1, X, |X|. Elementary computations show that Proj V (Y ) = α 0 + α 1 X + α 2 |X| with α 0 -0.752, α 1 = 0, α 2 2.196, and R 2 0.584. The model is misspecified in the sense that X 2 does not belong to V . It follows that the error term

ε = Y -Proj V (Y ) is such that E(ε|X) = X 2 -α 0 -α 2 |X| = 0 .
The estimated coverage levels for CI1 and CI2 are given in Table 4 below. We see that the estimated coverage level of CI2 is always around 0.87. Again, this confirms that CI2 is not robust to non normality. We see that CI1 is always better than CI2, with a coverage level greater than 0.94 for n ≥ 300.

We assume in the following that all the variables Y, X (1) , . . . X (p) belong to L 2 (H). Let V be the subspace of L 2 (H) generated by the variables (X (1) , . . . , X (p) ), and let Proj V (Y ) be the orthogonal projection of Y onto V . To measure the quality of approximation of Y by Proj V (Y ), one can consider the coefficient R 2 defined by (assuming that E( Y 2 H ) > 0):

R 2 = E Proj V (Y ) 2 H E Y 2 H . Let (Y i , X (1)
i , . . . X (p) i ) 1≤i≤n be n independent copies of the vector (Y, X (1) , . . . X (p) ), and let V i be the L 2 (H)-subspace generated by the variables (X

(1) i , . . . , X (p) i ). Let ε = Y -Proj V (Y ) and ε i = Y i - Proj V i (Y i ). We can then write Y i = Proj V i (Y i ) + ε i .
where the variables ε i are such that E( X

i , ε i ) = 0 for any j ∈ {0, . . . , p}. Assume that X (1) , . . . , X (p) are linearly independent (as a family of random variables in L 2 (H)). Then one can uniquely write

Y i = α 1 X (1) i + • • • + α p X (p) i + ε i . Let α = (α 1 , . . . , α p ) t ∈ R p ,
and M be the matrix such that M j, = E( X (j) , X ( ) ) . (

Since X (1) , . . . , X (p) are linearly independent, M is invertible and one can easily show that α = M -1 (E( X (1) , Y ) . . . , E( X (p) , Y )) t .

Let now M be the matrix such that M j, = n -1 n i=1 X Since M is invertible, for n large enough, M is invertible. Consequently, one can define

α = ( α 1 , . . . , α p ) t = M -1 1 n n i=1 X (1)
i , Y i , . . . ,

1 n n i=1 X (p) i , Y i t .
It follows that

α -α = ( α 1 , . . . , α p ) t -(α 1 , . . . , α p ) t = M -1 1 n n i=1 X (1)
i , ε i , . . . ,

1 n n i=1 X (p) i , ε i t . (6.3)
Then, proceeding as in Section 2.2, we see that α converges almost surely to α, and that: if for any j ∈ {0, . . . , p}, one has E( X (j) , ε 2 ) < ∞, then

√ n( α -α) L -→ n→∞ N p (0, M -1 M(ε)M -1 ) ,
where M(ε) is the matrix defined by M(ε) j, = E( X (j) , ε X ( ) , ε ) for any 0 ≤ j, ≤ p.

One can then estimate the matrix M(ε) as in Section 2.3. Define the residuals ε i by:

ε i = Y i -α 1 X (1) i -• • • -α p X (p) i .
Let then

M(ε) j, = 1 n n i=1 X (j) i , ε i X ( ) i , ε i .
As in Proposition 2.3, the following consistency result holds: assume that E(

X (j) 2 H ε 2 H ) < ∞ and E( X (j) 4
H ) < ∞ for any j ∈ {0, . . . , p}. Then, for any j, ∈ {0, . . . , p},

lim n→∞ M(ε) j, = M(ε) j, almost surely.
Consequently, since M is invertible,

lim n→∞ M -1 M(ε) M -1 = M -1 M(ε)M -1 almost surely. The empirical estimator of R 2 is R 2 = n i=1 Proj V i (Y i ) 2 H n i=1 Y i 2 H . Let θ = (θ 1 , . . . , θ p ) t , with θ = ( θ 1 , . . . , θ p ) t θ j = E( X (j) , Y ) E( Y 2 H ) and θ j = n i=1 X (j) i , Y i n i=1 Y i 2 H
for j ∈ {1, . . . , p}.

Then, as in Lemma 3.1 one can prove that R 2 = θ t α and R 2 = θ t α. Let e

(k) = X (k) -Proj W (X (k) ) and e (k) i = X (k) i -Proj W i (X (k) i ), (6.4)
where W is the subspace of L 2 (H) generated by Y , and W i is the subspace of L 2 (H) generated by Y i . Let also e (k) = (e

(k) 1 , . . . , e (k)
n ) t . As in Proposition 3.2, one can prove the following result: if E( X (j) , ε 2 ) < ∞ and E( Y, e (j) 2 ) < ∞ for all j ∈ {1, . . . , p}. Then

√ n R 2 -R 2 L -→ n→∞ N (0, V ).
where V = (θ 1 , . . . , θ p , α 1 , . . . , α p )BAB(θ 1 , . . . , θ p , α 1 , . . . , α p ) t .

Here, the matrices A and B are as follows:

• If (j, k) ∈ {1, . . . , p} 2 then A j,k = E( X (j) , ε X (k) , ε ); • If (j, k) ∈ {p + 1, . . . , 2p} 2 then
A j,k = E( Y, e (j-p) Y, e (k-p) );

• (j, k) ∈ {1, . . . , p} × {p + 1, . . . , 2p} then

A j,k = E( X (j) , ε Y, e (k-p) ).
Let δ j,k = 0 if j = k and δ j,j = 1 and let B be the 2p × 2p matrix defined as follows:

• If (j, k) ∈ {1, . . . , p} 2 then B j,k = (M -1 ) j,k (see (6.1) for the definition of M);

• If (j, k) ∈ {p + 1, . . . , 2p} 2 then B j,k = (E( Y 2 H )) -1 δ j,k ; • If (j, k) ∈ {1, . . . , p} × {p + 1, . . . , 2p} then B j,k = 0.
Now the matrices A and B can be consistently estimated by their empirical counterparts based on the residuals, provided E( Y 4 H ) < ∞ and E( X (j) 4 H ) < ∞ for any j ∈ {1, . . . , p} (as done in the proof Proposition 3.3). As in Proposition 3.3 this provides a consistent estimator of the limiting variance V .

An alternative to the usual sensitivity analysis

In sensitivity analysis, one usually assumes that Y = m(X (1) , . . . , X (p) ) (

(or that Y = m(X (1) , . . . , X (p) ) + ε as in [START_REF] Huet | Metamodel construction for sensitivity analysis ESAIM[END_REF]) where (X (1) , . . . , X (p) ) is a p random vector with a known distribution = P 1 ⊗ • • • ⊗ P p on a compact subset X of R p . The number of variables p may be large, and the function m is an unknown function from R p to R, which may present strong non-linearities and high order interaction effects between its coordinates. On the basis of a n-sample (Y i , X

i , . . . , X

i ), for i = 1, . . . , n, we want to measure the impact of the main effects or interactions, i.e. to determine the influence of each variable or group of variables on the output variable Y . The usual sensitivity analysis is mainly based on Hoeffding's decomposition, strongly related to the independency between the inputs variables X (1) , . . . , X (p) . In that latter case, if m is square integrable, one may consider the classical Hoeffding-Sobol decomposition [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF] that leads to write m according to its ANOVA functionnal expansion: m(X (1) , . . . , X

(p) ) = m 0 + i 1 < i 2 < • • • < i k , i i , . . . , i k ∈ {1, . . . , p} m i 1 ,...,i k (X (i 1 ) , . . . , X (i k ) ) (6.6)
where the functions m i 1 ,...,i k are centered and orthogonal in L 2 involving conditional expectation that is

m 0 = E(m(X (1) , . . . , X (p) )), m j (X (j) ) = E(m(X (1) , . . . , X (p) )|X (j) ) -m 0 m j,k (X (j) , X (k) ) = E(m(X (1) , . . . , X (p) )|X (j) , X (k) ) -m j (X (j) ) -m k (X (k) ) -m 0 , . . .
Under the assumption of independence on the X (j) 's this decomposition is unique, with E(m i 1 ,...,i k (X (i 1 ) , . . . , X (i k ) )) = 0 and for all (i 1 , . . . , i k ) = (i 1 , . . . , i ), E m i 1 ,...,i k (X (i 1 ) , . . . , X (i k ) )m i 1 ,...,i (X (i 1 ) , . . . , X (i ) ) = 0 .

Then, given a threshold γ n we keep the covariates such that:

|cor(Y, X (j) )| ≥ γ n .

To assess the performance of a screening method, whether it is based on computing p correlations as we mentioned or not, one has to ask two questions:

1. Is the method able to retrieve all the covariates that are associated with the outcome Y ? This question is known in the literature as the Sure Screening Property (see Fan and Lv [1]).

2. Is the method specific enough, i.e. can we ensure that not too many covariates are wrongly selected? This question directly relates to the control of the False Discovery Rate (FDR).

The questions of the sure screening property and the control of the FDR in a model-free environment is a growing field of research (see for instance Fan and Lv [START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF], Zhao and Li [START_REF] Zhao | Principled sure independence screening for Cox models with ultrahigh-dimensional covariates[END_REF], Fan and Song [START_REF] Fan | Sure independence screening in generalized linear models with NPdimensionality[END_REF], W. Pan, X. Wang, W. Xiao, H. Zhu [START_REF] Pan | A Generic Sure Independence Screening Procedure[END_REF]) To answer both questions, the choice of the threshold γ n is crucial. One paper which is often cited is the paper by Fan and Lv [START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF]. Their context is the following: they suppose a true joint linear model:

Y = p j=1 X (j) β j + ε
where ε is a standard Gaussian error. Moreover, they suppose that for each j = 1, . . . , p, the X

i are i.i.d (respectively to i = 1, . . . , n). They also suppose that X (j) are standardized with mean 0 and standard deviation 1. They consider the screening set:

M = j : |cor(Y, X (j) )| ≥ n log(n)
.

Ensuring the Sure Screening Property in that specific case consists in showing that the support of β is at least included in M with a probability close to one. However, the threshold chosen γ n = n/ log(n) does not limit the number of false positives, i.e., the covariates that are not in the support of β but are selected marginally. This threshold, although not having any theoretical support for controlling the FDR, is widely used in the screening literature. The two questions 1. and 2. above are addressed for instance in Zhao and Li [START_REF] Zhao | Principled sure independence screening for Cox models with ultrahigh-dimensional covariates[END_REF] in the context of the Cox model. In this paper the authors propose to choose M = j :

√ n|α j | vj ≥ γ ,
where αj corresponds to the estimator in the marginal (and possibly misspecified) model, and vj is a consistent estimator of the asymptotic variance of √ n(α j -α j ). According to Zhao and Li [START_REF] Zhao | Principled sure independence screening for Cox models with ultrahigh-dimensional covariates[END_REF], denoting by Φ the distribution function of a standard Gaussian and by q the expected false positive rate, the idea is to choose

γ = Φ -1 1 - q 2 .
It is common in papers dealing with screening to consider that independence beetwen X (j) being truly associated with Y in the joint model is equivalent to being associated with Y marginally (see Condition Our idea. First, we highlight the fact that even in a model-free context, the question of screening the covariates is still relevant as a lowering dimension step before trying to fit a joint association between the remaining covariates and the outcome.

Recall from Section 4 the expression of τ j and τj given in Equation (4.1) and (4.2)

τ j = Cov(Y, X (j) )
Var(X (j) ) and

τ j = n i=1 (Y i -Y)(X (j) i -X (j) ) n i=1 (X (j)
i -X (j) ) 2 .

Then, Proposition 4.1 gives √ n( τ j -τ j )

L -→ n→∞ N (0, v j )
where v j = E((X (j) -E(X (j) )) 2 ε (j) 2 ) (Var(X (j) )) 2 .

Moreover, it has been shown in Section 4 that if we let

v j = 1 n n i=1 (X (j)
i -X (j) ) 2

(Var n (X (j) )) 2 then, for all j = 1, . . . , p, lim n→∞ v j = v j almost surely.

Similarly to Zhao and Li [START_REF] Zhao | Principled sure independence screening for Cox models with ultrahigh-dimensional covariates[END_REF], a reasonnable screening set rule should be M = j :

√ n| τ j | v j ≥ γ with γ = Φ -1 1 - q 2 .
The next step is to show that this rule does indeed allow us to obtain the sure screening property, as well as satisfactory control of the FDR. We believe that moment assumptions on the variables should be sufficient to show these two properties, and we plan to study these issues in future work.

  .21) Combining (3.20) and (3.21), and using that | β -β | and | θ -θ | converge almost surely to 0, we infer that A j,k -Ãj,k converges almost surely to 0 as soon as, almost surely

. 6 )

 6 Moreover, by the strong law of large numbers, lim n→∞ D = D almost surely (4.7)

  Proof of Proposition 4.2. Let φ : R 2p → R p be such that φ(x 1 , . . . , x p , y 1 , . . . , y p ) = (x 1 y 1 , . . . , x p y p ) t .

	Since
	( R 2 (1) . . . , R 2 (p) ) t -(R 2 (1) , . . . , R 2 (p)

2 (p) ) t -(S 2 (1) , . . . , S 2 (p) ) t may be derived from the asymptotic distribution of √ n [ τ : θ : α] -[τ : θ : α] via the delta method.

Table 4 :

 4 Estimated coverage levels of CI1 and CI2 at level 95% for model(5.3) with Proj V (Y ) -0.752 + 2.196|X|.

	200	300	400	500	600	700	800	900	1000
	CI1 0.935 0.941 0.943 0.943 0.943 0.948 0.947 0.945 0.947
	CI2 0.875 0.871 0.885 0.868 0.861 0.87 0.866 0.874 0.878

in[START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF] and Assumption 8 in Zhao and Li[START_REF] Zhao | Principled sure independence screening for Cox models with ultrahigh-dimensional covariates[END_REF])). However, if we are working in a model-free context as in Section

[START_REF] Huet | Metamodel construction for sensitivity analysis ESAIM[END_REF], these assumptions are unnecessary and it seems possible to define a threshold γ using the distribution function Φ.

An example of polynomial regression

We consider again the model (5.3), but now V is the sub-space of L 2 generated by 1, X, X 2 . It is then

The estimated coverage levels for CI1 and CI2 are given in Table 5 below. We see that the estimated coverage level of CI2 is always around 0.75. Again, this confirms that CI2 is not robust to non normality. We see that CI1 is always much better than CI2, with a coverage level greater than 0.93 for n ≥ 500 and greater than 0.94 for n ≥ 800.

An example of Poisson regression

In this example, V is the subspace of L 2 generated by 1, X (1) , X (2) , where X (1) , X (2) are independent, X (1) is uniformly distributed over [0, 1], and X (2) is exponentially distributed with parameter 1. The response variable Y is a count variable, whose conditional distribution given X (1) , X (2) is a Poisson distribution with parameter 0.5 + X (1) + X (2) . In that case, we have Proj V (Y ) = E(Y |X (1) , X (2) (2) with α 0 = 0.5, α 1 = 2, α 2 = 1, and R 2 = 13/37. Note that this is again an heteroscedastic model, since Var(ε|X (1) , X (2) ) = Var(Y -Proj V (Y )|X (1) , X (2) ) = 0.5 + X (1) + X (2) .

The estimated coverage levels for CI1 and CI2 are given in Table 6 below. We see that the estimated coverage level of CI2 is always around 0.86. Again, this confirms that CI2 is not robust to non normality. We see that CI1 is always much better than CI2, with a coverage level greater than 0.93 for n ≥ 300 and greater than 0.94 for n ≥ 600.

6 Extension to higher dimension and discussions

Extension to vector-valued random variables

In this section, we assume that (Y, X (1) , . . . X (p) ) is a random vector with value in H p+1 , where H is a real separable Hilbert space with inner product •, • and norm • H . Let L 2 (H) be the space of H-valued random variables Z such that E( Z 2 H ) < ∞, and recall that L 2 (H) is also a separable Hilbert space with inner product E( •, • ). This leads to the following variance decomposition:

Var m i 1 ,...,i k (X (i 1 ) , . . . , X (i k ) ) .

The Sobol sensitivity indices introduced by Sobol [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF] are defined for any group (X (i 1 ) , . . . , X (i k ) ) by

They quantify the contribution of a subset of variables (X (1) , . . . , X (p) ) to the output Y . Several approaches are available for estimating these sensitivity indices, see for example Iooss and Lemaître [START_REF] Iooss | Uncertainty management in Simulation-Optimization of Complex Systems: Algorithms and Applications[END_REF] for a recent review. An alternative approach, which is much less computationally expensive, is as follows. Recall that the X (i) 's are bounded and let X(i) = X (i) -E(X (i) ) (this centering may be done, since the distribution of (X (1) , . . . , X (p) ) is known). Since the variables X (j) 's are independent, the family

is an orthogonal family of L 2 . We can then compute the individual R 2 's for this family, that is

where V (i 1 ,...,i k ) is the subspace of L 2 generated by the product X(i 1 )

The main difference between the two approaches is that the sum of all the indices R i 1 ,...,i k is not equal to 1, unless Y belongs to the subspace V F of L 2 generated by the family F. More precisely, we have that

The advantage of using these R 2 's is that they are very easy to estimate (based on n independent observations (Y i , X

i , . . . , X

i ) 1≤i≤n ). In addition, as we have seen, we can give confidence intervals for these quantities.

A first step toward robust screening

In this section, we outline a strategy for robust screening based on the results of Section 4.

Quick context

The basic idea of screening is to find, among a set of covariates X (1) , . . . , X (p) , the ones having an association with an outcome Y . One usual way which is often cited consists in computing the p correlations: cor(Y, X (j) ) for 1 ≤ j ≤ p.