Statistique en grande dimension

Frédéric Lavancier

Université de Nantes
M2 Ingénierie Statistique

2020/2021
Statistique en grande dimension F. Lavancier

Références

• ESL : "The elements of statistical learning", T. Hastie, R. Tibshirani,

J. Friedman.

• "An introduction to statistical learning with applications in R", G. James,

D. Witten, T. Hastie, R. Tibshirani.

• "Régression avec R", P-A. Cornillon, E. Matzner-Løber

• "Introduction to high-dimensional statistics", C. Giraud.

La plupart des figures de ce document sont tirées du premier ou-

vrage ESL ci-dessus, avec la permission des auteurs.

1
Contents

1 Introduction 4

2 Rappels sur la régression linéaire 7

2.1 La régression linéaire est une projection . . . . . . . . . . . . . 7
2.2 Modélisation statistique . . . . . . . . . . . . . . . . . . . . . 10
2.3 Qualité de l’estimation et validation du modèle . . . . . . . . 11
2.4 Aspects en grande dimension . . . . . . . . . . . . . . . . . . . 14

3 Choix de modèles 17
3.1 Erreurs de prévision . . . . . . . . . . . . . . . . . . . . . . . 18
3.2 Critères usuels . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2.1 Cp de Mallows . . . . . . . . . . . . . . . . . . . . . . . 22
3.2.2 Critère AIC . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2.3 Critère BIC . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2.4 Comparaison des critères . . . . . . . . . . . . . . . . . 26
3.3 Validation croisée . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.1 Hold out . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.2 Leave-one-out (LOO) . . . . . . . . . . . . . . . . . . . 28
3.3.3 Leave-k-out . . . . . . . . . . . . . . . . . . . . . . . . 30
3.3.4 K-fold . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.4 Sélection automatique d’un sous-modèle . . . . . . . . . . . . 32
3.5 Estimation du modèle retenu . . . . . . . . . . . . . . . . . . . 34

4 Réduction de dimension et régression sous contraintes 35

4.1 Réduction de dimension : PCR et PLS . . . . . . . . . . . . . 36
4.1.1 Régression sur composantes principales (PCR) . . . . . 36
4.1.2 Régression des moindres carrés partiels (PLS) . . . . . 39
4.2 Régression ridge . . . . . . . . . . . . . . . . . . . . . . . . . . 41

2
Statistique en grande dimension F. Lavancier

4.3 Régression Lasso . . . . . . . . . . . . . . . . . . . . . . . . . 43

4.3.1 Principe . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.3.2 Résolution . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.3.3 Aparté R : comparaison entre glmnet et lars . . . . . 52
4.3.4 Aspects théoriques . . . . . . . . . . . . . . . . . . . . 55
4.4 Quelques généralisations des régressions Ridge et Lasso . . . . 60
4.4.1 Elastic net . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.4.2 Gauss-Lasso . . . . . . . . . . . . . . . . . . . . . . . . 63
4.4.3 Adaptive Lasso . . . . . . . . . . . . . . . . . . . . . . 63
4.4.4 Logistic Lasso . . . . . . . . . . . . . . . . . . . . . . . 64
4.4.5 Group-Lasso . . . . . . . . . . . . . . . . . . . . . . . . 65

5 Tests multiples 67
5.1 Présentation du problème et notations . . . . . . . . . . . . . 68
5.2 Principe des tests multiples : contrôler le F W ER ou le F DR 71
5.2.1 Procédure de Bonferroni . . . . . . . . . . . . . . . . . 72
5.2.2 Procédure de Benjamini-Hochberg . . . . . . . . . . . . 73
5.2.3 p-values ajustées . . . . . . . . . . . . . . . . . . . . . 78

3
Chapter 1

Introduction

Lors d’une étude statistique, on dispose généralement de n individus pour

lesquels on a observé/relevé p variables (quantitatives ou qualitatives).
Qu’est-ce que la grande dimension? Plusieurs situations sont possibles
selon que n et/ou p sont grands.

n grand, p de taille raisonnable

Exemple 1: Les clients d’une banque. Leur nombre n est très grand, mais
leurs caractéristiques individuelles sont en nombre limité.
Exemple 2: "Data Lake" Renault. Chaque concessionnaire du monde
entier saisit ses factures du jour en ligne. Le "Data Lake" de Renault contient
ainsi les factures de chaque véhicule, pour chaque concessionnaire et chaque
jour de l’année. Si on s’intéresse par exemple aux factures, leur nombre
n est immense, mais les variables les caractérisant (nom du véhicule, type
d’intervention, montant, etc) sont en nombre p limité (quelques dizaines au
maximum).
C’est en général une situation favorable d’un point de vue théorique, car
plus il y a d’individus et meilleures sont les estimations. De plus, si n est très
grand, on peut utiliser sans risque les résultats asymptotiques disponibles,
principalement le théorème limite central, pour quantifier la qualité des esti-
mations ou des prévisions via des intervalles de confiance.
En fait, le principal problème lorsque n est très grand est d’ordre infor-
matique. Le stockage peut ne plus être possible sur un simple ordinateur,
ce qui complique l’accès aux données, et les calculs sont infaisables sur toute
la population en un temps raisonnable. Des solutions informatiques de type

4
Statistique en grande dimension F. Lavancier

"Big Data" sont alors nécessaires (clusters de serveurs, données distribuées,

technologie Hadoop) mais ce n’est pas le sujet de ce cours.

n de taille raisonnable, p grand

Exemple 1: Text Mining. Pour l’analyse des textes (par exemple issus de la
satisfaction clients ou de réponses à un questionnaire), on liste tous les "mots"
utilisés dans l’ensemble des textes, et on forme un tableau récapitulant pour
chaque texte le nombre d’apparitions de chaque mot. Chaque colonne du
tableau correspond à un mot et chaque ligne à un texte analysé. Dans ce
contexte, n correspond au nombre de textes analysés et p au nombre total
de mots, qui est généralement très grand.
Exemple 2: En génomique. Pour chaque individu, on dispose de "l’expression"
de certains de ses gènes. Il s’agit d’une valeur numérique par gène. Le nombre
p de gènes "exprimés" est de l’ordre de 105 ou 106 . Etant donné n individus
observés, l’objectif peut-être d’essayer de former des groupes homogènes en
terme d’expression des gènes (classification non supervisée), ou en supposant
qu’une partie des individus est malade et l’autre non, de tenter de trouver
les gènes s’exprimant différemment pour les individus malades.
Exemple 3: En analyse de signaux. Lorsqu’on observe le résultat d’une
analyse par spectrométrie, la réponse est une valeur de spectre par longueur
d’ondes. Il y a en général plusieurs milliers de longueurs d’onde mesurées
par individu.
Dans ces situations, on dispose d’un très grand nombre de variables p
par rapport à n, par exemple p > n/2, voire même p >> n. Cette sit-
uation pose problème pour les analyses statistiques classiques, notamment
pour la régression linéaire. En effet, supposons que l’on souhaite déterminer
le lien entre une variable d’intérêt et les autres variables à disposition, on ne
peut raisonnablement pas espérer estimer un modèle ayant plus de variables
que d’observations disponibles. Deux éléments clés permettent néanmoins
l’analyse statistique:
1. On suppose que parmi les nombreuses variables explicatives à dispo-
sition, seule une petite partie d’entre elles est en réalité utile pour le
modèle (sans que l’on sache lesquelles a priori). Il s’agit de l’hypothèse
fondamentale de parcimonie.
2. Des méthodes adaptées (choix de modèles, modèles contraints) perme-
ttent d’exploiter la parcimonie et d’estimer convenablement le modèle.

5
Statistique en grande dimension F. Lavancier

Ce contexte forme ce que l’on appelle généralement la statistique en grande

dimension, et il s’agit de la situation dont traite principalement ce cours.

n et p grand
Il s’agit bien évidemment de la situation la plus compliquée : chaque individu
est associé à un très grand nombre de variables, et il y a un très grand nombre
d’individus observés. On peut penser aux exemples précédents (text mining,
génomique, spectres) en imaginant dans chaque cas n très grand.
Dans cette situation, on utilise à la fois des outils informatiques adaptés
au Big Data pour le stockage et l’accès aux données, et des outils statistiques
de grande dimension pour intégrer le grand nombre de variables dans les
modèles.

6
Chapter 2

Rappels sur la régression linéaire

L’objectif de ce chapitre n’est pas d’introduire la régression linéaire, que l’on

suppose avoir déjà été étudiée, mais de rappeler ses principes de base, avec
un regard particulier pour la situation où le nombre de variables explicatives
est important, voire supérieur au nombre d’observations.

2.1 La régression linéaire est une projection

On dispose de n observations y1 , . . . , yn de la variable à expliquer et de n
observations x1 , . . . , xn des variables explicatives. Pour tout i, yi ∈ R est une
valeur numérique réelle et xi ∈ Rp est un vecteur de taille p.
On note
     
y1 x01 x11 . . . xp1
Y =  ...  et X =  ...  =  ... .. ..  .
    
. . 
0
yn xn x1n . . . xpn

Les p variables explicatives, observées sur les n individus, sont donc

   
x11 xp1
 ..   .. 
X1 =  .  , · · · , Xp =  .  .
x1n xpn

On suppose que X, de dimension (n, p), est de plein rang, c’est à dire
rang(X) = min(n, p). Lorsque p ≤ n, cela signifie qu’il n’existe aucune
combinaison linéaire entre les variables explicatives. Lorsque p > n, cela

7
Statistique en grande dimension F. Lavancier

signifie que la famille des p variables explicatives génère Rn , autrement dit il

existe au moins n variables parmi les p qui sont non linéairement liées entre
elles.

On note [X] = {Xα, α ∈ Rp } = {v ∈ Rn , ∃α ∈ Rp , v = Xα} l’espace

vectoriel engendré par les colonnes de X. Il est de dimension rang(X).

• Si p < n, [X] est un sous espace vectoriel de Rn de dimension p.

• Si p ≥ n, [X] = Rn .

On note [X]⊥ l’espace vectoriel orthogonal à [X] dans Rn , c’est à dire

[X] = {v ∈ Rn , X 0 v = 0}.
⊥

• Si p < n, il est de dimension n − p.

• Si p ≥ n, [X]⊥ = {0}.

Régresser Y ∈ Rn sur les colonnes de X consiste à projeter Y sur [X]. On

note P[X] la matrice de cette projection. On note également P[X]⊥ la matrice
de projection sur [X]⊥ . On a la relation P[X] + P[X]⊥ = In , où In désigne la
matrice identité de taille n

• Si p < n, P[X] = X(X 0 X)−1 X 0 . Dans ce cas, X 0 X, matrice carrée de

taille p, est de rang p et est donc bien inversible. De plus, P[X]⊥ =
In − P[X] = In − X(X 0 X)−1 X 0 .

• Si p ≥ n, P[X] = In et P[X]⊥ = In − P[X] = 0n (la matrice nulle de taille

n).

Le projeté est noté Ŷ , i.e. Ŷ = P[X] Y . Il est obtenu par la méthode des
moindres carrés ordinaires (MCO) au sens où

Ŷ = argminZ∈[X] kY − Zk2

ou de façon équivalent Ŷ = X β̂ avec

β̂ = argminβ∈Rp kY − Xβk2 .

8
Statistique en grande dimension F. Lavancier

Figure 2.1: Figure extraite de l’ouvrage ESL montrant le nuage des points
(xi , yi ) en rouge lorsque xi est de dimension p = 2, autrement dit X =
(X1 , X2 ). Le plan représenté est celui des moindres carrés, auquel appartient
les ŷi .

• Si p < n, Ŷ = X β̂ où β̂ = (X 0 X)−1 X 0 Y . Cette solution est unique.

• Si p ≥ n, Ŷ = Y = X β̂ où β̂ représente simplement les coordonnées

de Y dans la famille génératrice de Rn formée des colonnes de X. Si
p = n, β̂ est unique car cette famille est alors une base de Rn . Si p > n,
il y a une infinité de façon de choisir β̂ (par exemple en annulant p − n
coordonnées).

Les figures 2.1 et 2.2 illustrent le principe de la régression linéaire dans le

cas p = 2.

9
Statistique en grande dimension F. Lavancier

Figure 2.2: Figure extraite de l’ouvrage ESL. Le plan en jaune représente

[X] lorsque p = 2. Le vecteur Y est projeté sur [X] pour donner Ŷ .

2.2 Modélisation statistique

Dans la partie précédente, il n’est pas question de loi de probabilité ni de
statistique : étant donné un vecteur Y , on se contente de le projeter sur
l’espace vectoriel [X].
Le cadre statistique est le suivant : on suppose que chaque observation yi
est la réalisation d’une variable aléatoire ayant pour espérance x0i β pour un
certain β ∈ Rp et ayant même variance σ 2 . De plus les yi sont supposés non
corrélées entre eux. Ces hypothèses reviennent à affirmer que

Y = Xβ +  (2.1)

où  = Y − Xβ est un vecteur aléatoire d’espérance nulle et de covariance

σ 2 In . Dans ce formalisme les variables explicatives X sont supposées non
aléatoires.
Dit autrement, on suppose que E(Y ) ∈ [X] mais Y n’appartient pas à [X]
(sauf si [X] = Rn , le cas p ≥ n). Supposer que Y ∈ [X] serait une hypothèse
extrêmement forte et peu réaliste. Dans ce cadre, régresser Y sur [X] revient
donc à estimer les coordonnées de E(Y ) dans [X].

10
Statistique en grande dimension F. Lavancier

L’intérêt d’une telle modélisation peut être de deux types :

• Descriptif : comprendre le lien entre E(Y ) et les variables explicatives

(quelles sont les variables qui, en moyenne, influencent le plus Y et
comment?). Cela nécessite une estimation précise de β dans le sens où
Ekβ̂ − βk2 doit être la plus petite possible.

• Prédictif : étant donné un nouvel individu o dont on connait les valeurs

xo des variables explicatives, on prédit yo via l’estimation de E(yo ) c’est
à dire ŷo = x0o β̂. Une bonne prédiction doit minimiser E(yo − ŷo )2 . Il est
à noter qu’une prévision
1
Pn naïve de E(yo ) aurait pu être simplement la
moyenne empirique n i=1 yi . Cette prévision est en général fortement
biaisée car les yi n’ont pas même espérance.

Remarque 2.2.1. Si la relation (2.1) n’est pas vérifiée, la régression de Y

sur [X] peut toujours être effectuée pour obtenir Ŷ = P[X] Y . Mais dans ce
cas Ŷ est généralement un mauvais estimateur de E(Y ). Trouver la bonne
matrice X (c’est à dire choisir les bonnes variables explicatives) telle que
(2.1) est vérifiée est donc primordial. C’est l’objet de la sélection de modèles
abordée en section 2.3 et discutée plus en détail dans le chapitre suivant.

2.3 Qualité de l’estimation et validation du mod-

èle
On suppose dans cette partie que la relation (2.1) est vérifiée avec  centré
de variance σ 2 In . On suppose de plus que p < n et X est de plein rang, i.e.
rang(X) = p. Ces hypothèses sont le cadre d’étude classique (favorable) de
la régression linéaire. On résume les principales propriétés qui en découlent
ci-dessous. Le cas p ≥ n est discuté en section 2.4.

Dans le cadre précédent, on a les propriétés suivantes:

• β̂ = (X 0 X)−1 X 0 Y est sans biais, i.e. E(β̂) = β, et de variance V ar(β̂) =

σ 2 (X 0 X)−1 .

• β̂ est le meilleur estimateur linéaire sans biais de β au sens du coût

quadratique (Théorème de Gauss-Markov).

11
Statistique en grande dimension F. Lavancier

• Si l’on suppose de plus que  suit une loi normale, alors β̂ correspond
à l’estimateur du maximum de vraisemblance de β.

Les résidus de la régression sont ˆ = Y − Ŷ . Leur norme correspond à

la somme des carrés des résidus : SCR = kˆk2 . Elle permet d’estimer sans
biais la variance σ 2 de  :
n
2 SCR 1 X
σ̂ = = (yi − ŷi )2 .
n−p n − p i=1

Les résidus sont à la base de toutes les procédures de test et d’analyse de la

qualité de la modélisation, résumées ci-dessous.

Les principaux tests :

• Le test de Student de significativité du paramètre βj consiste à tester

H0 : βj = 0 contre H1 : βj 6= 0. Il utilise la statistique

β̂j
T =
σ̂β̂j
q
où σ̂β̂j désigne l’écart-type de β̂j qui vaut σ̂β̂j = σ̂ (X 0 X)−1
jj , c’est à

dire la racine carrée du j-ème terme de la diagonale de V ar(β̂).

La région critique du test au niveau α est RCα = {|T | > tn−p (1−α/2)}
où tn−p (1 − α/2) est le quantile d’ordre 1 − α/2 d’une loi de Student de
degré de liberté n − p. Le niveau du test est exactement α si la loi de 
est une loi normale, et vaut asymptotiquement α (i.e. lorsque n → ∞)
sinon (et sous certaines hypothèses de régularité).

• Le test de Fisher de contraintes linéaires sur les paramètres consiste

à tester H0 : Rβ = 0 contre H1 : Rβ 6= 0 où R est une matrice de
q contraintes de taille (q, p). Un cas particulier courant est le test
des modèles emboités où les contraintes sont la nullité simultanée de q
paramètres dans le modèle général. La statistique de test est
n − p SCRc − SCR
F = ,
q SCR

12
Statistique en grande dimension F. Lavancier

où SCR désigne la SCR du modèle sans contrainte et SCRc celle

du modèle contraint. Si  suit une loi normale, une région critique au
niveau α est RCα = {F > fq,n−p (1 − α)} où fq,n−p (1 − α) est le quantile
d’ordre 1−α d’une loi de Fisher de degré de liberté (q, n−p). De même
le niveau vaut asymptotiquement α si  ne suit pas une loi normale (et
sous certaines hypothèses de régularité).
• Sans entrer dans les détails, nous mentionnons le test d’homoscédasticité
de Breusch-Pagan (la matrice de variance de  est-elle bien constante
sur la diagonale?), et les tests d’auto-corrélation de Durbin-Watson et
de Breusch-Godfrey (la matrice de variance de  est-elle bien diago-
nale?).

Les mesures usuelles de la qualité de modélisation :

• Le coefficient de corrélation multiple (ou coefficient de détermination),
dans le cas où X contient la colonne constante égale à 1 (autrement dit
si le modèle contient une constante) :
SCE SCR
R2 = =1−
SCT SCT
où SCE = kŶ − Ȳ k2 est la somme des carrés expliqués et SCT =
kY − Ȳ k2 est la somme des carrés totaux. Lorsque le modèle ne contient
pas de constante, R2 = 1 − SCR/kY k2 . Plus le R2 est proche de 1,
plus le modèle s’ajuste bien aux données.
• Le R2 corrigé (ajusté) :
n − 1 SCR
Ra2 = 1 −
n − p SCT
dans le cas où le modèle contient une constante. Le Ra2 s’interprète
comme le R2 sans souffrir de son principal défaut, qui est que le R2
croît nécessairement lorsqu’on ajoute une variable au modèle même si
cette dernière est non significative.

Les critères usuels de choix de modèles : si on hésite entre plusieurs mod-

èles, le but étant de trouver le modèle le plus parcimonieux qui vérifie (2.1),
on peut comparer les critères suivants (voir le chapitre suivant pour plus de
détails)

13
Statistique en grande dimension F. Lavancier

• Ra2 : on privilégie le modèle ayant le Ra2 le plus élevé.

• Le Cp de Mallows
SCR
Cp = − n + 2p
σ̂ 2
où σ̂ 2 est un estimateur sans biais de σ 2 (en général celui calculé sur le
plus gros modèle) et on privilégie le modèle ayant le Cp le plus faible.

• le critère AIC
SCR
AIC = n log + 2(p + 1)
n
et on privilégie le modèle ayant le AIC le plus faible.

• le critère BIC
SCR
BIC = n log + (p + 1) log n
n
et on privilégie le modèle ayant le BIC le plus faible.

Enfin, on peut procéder à un choix automatique basé sur l’un des critères
précédents à l’aide d’une procédure pas à pas ascendante (forward stepwise
selection) de la façon suivante :

1. On part du modèle le plus petit (contenant uniquement la constante)

2. On ajoute la variable conduisant à la meilleure amélioration possible

selon le critère choisi (par exemple la plus grande diminution du BIC)

3. On recommence jusqu’à ce qu’aucun nouvel ajout de variable n’améliore

plus le critère.

De la même manière, on peut effectuer une procédure pas à pas descendante

(backward stepwise selection) en partant du plus gros modèle et en enlevant
successivement la "pire" variable.

2.4 Aspects en grande dimension

Il est intéressant d’observer ce qui se passe en grande dimension : si p > n,
le modèle linéaire
Y = Xβ + 

14
Statistique en grande dimension F. Lavancier

est mal défini puisqu’il existe une infinité de coefficients β conduisant au

même vecteur v = Xβ. Par exemple v peut s’exprimer en fonction des n
premières colonnes de X, ou des n dernières, etc.
Chercher à estimer β n’a donc pas vraiment de sens. D’ailleurs, β̂ obtenu
par les MCO n’est pas non plus unique, voir la section 2. Aborder l’analyse
d’un tel modèle en toute généralité parait donc vain.
La situation devient plus raisonnable si l’on suppose que dans la relation
linéaire précédente, la plupart des coefficient βj sont en fait nuls. Pour fixer
les idées, supposons que seuls p∗ coefficients βj sont non nuls, avec p∗ << n.
Il s’agit d’une hypothèse de parcimonie (sparsity) cruciale permettant
l’analyse statistique en grande dimension. Elle revient à supposer que Xβ
appartient en réalité à un sous-espace vectoriel [X ∗ ] de [X] de dimension
p∗ << n, où X ∗ désigne la sous matrice de X ne contenant que les p∗ colonnes
associées aux coefficients βj non-nuls.
Si on connaissait X ∗ , on serait dans la situation favorable de la section 2.3.
La difficulté réside dans le fait que l’on ne connait pas p∗ et encore moins
X ∗ . Il convient donc d’utiliser des outils pour trouver X ∗ à partir de X et
estimer le lien linéaire associé.
Concernant les critères standards de sélection de modèles présentés dans
la partie précédente, ils s’appuient tous sur les résidus. Mais si l’on part d’un
modèle ayant p ≥ n variables, alors Ŷ = Y et donc ˆ = 0, et ce même en
absence de lien réel entre Y et X. En conséquence, pour un tel modèle :
• R2 = 1,

• Ra2 n’a plus de sens (car n − p < 0),

• Cp n’a plus de sens (car à la fois le SCR du modèle estimé et σ̂ 2 calculé

sur le modèle le plus gros sont nuls),

• AIC = −∞,

• BIC = −∞.
Toutes les méthodes de sélection de modèles présentées précédemment ne
sont donc plus utilisables.
Il y a néanmoins une exception : la sélection forward. En effet cette
dernière part du plus petit modèle, pour lequel les critères de sélection sont
calculables (à part le Cp de Mallows qui s’appuie sur l’estimation de σ̂ 2 à
partir du plus gros modèle). Ces critères restent calculables et comparables

15
Statistique en grande dimension F. Lavancier

pour les premières étapes de l’algorithme pour lesquelles un nombre limité

de variables a été ajouté. Par l’hypothèse de parcimonie, on peut espérer
que cet algorithme s’arrête avant d’avoir inclus un nombre trop important
de variables dans le modèle.

Ainsi en présence d’un grand nombre de variables explicatives (p > n

mais également p ≤ n avec p grand), le modèle étant supposé par ailleurs
parcimonieux, il est nécessaire d’introduire des méthodes alternatives aux
MCO afin d’estimer au mieux le modèle linéaire. C’est ce que nous verrons au
chapitre 4. Mais quelle que soit la méthode utilisée, il conviendra d’adopter
une stratégie de sélection de modèles. Cette étape est commune à toute
modélisation statistique, qu’elle concerne un modèle linéaire, logistique, ou
autre. Le chapitre suivant clarifie les méthodes usuelles pour y parvenir.

16
Chapter 3

Choix de modèles

Pour fixer les idées on se place dans le cadre d’un modèle de régression
linéaire. Comme on l’a vu dans le premier chapitre, la présence de nom-
breuses variables explicatives (p grand par rapport à n) rend l’estimateur
par moindres carrés peu fiable car sa variance explose (en effet, si X de di-
mension (n, p) est de rang n < p, alors X 0 X est une matrice (p, p) de rang n et
est donc non inversible). En particulier les prévisions basées sur l’estimateur
par MCO sont mauvaises. Ce phénomène n’est pas propre à la régression
linéaire mais s’observe dans la plupart des modèles statistiques en grande di-
mension. Il est assez naturel dans ce contexte de supposer de la parcimonie
: un faible nombre de variables sont en réalité pertinentes dans le modèle,
mais on ne sait pas a priori lesquelles.
La façon la plus naturelle de procéder dans ce contexte est de choisir au
mieux un sous-modèle. Cela nécessite d’estimer pour chaque modèle l’erreur
associée, ce qui est détaillée dans la partie 3.1. Différentes approches exis-
tent pour estimer cette erreur : des critères, comme AIC ou BIC, dont la
construction est rappelée dans la partie 3.2, et la validation croisée discutée
dans la partie 3.2.4. Les méthodes de sélection consistent à parcourir un
grand nombre de sous-modèles possibles (idéalement tous) et de retenir celui
dont l’erreur estimée est la plus faible. La partie 3.4 résume les approches
classiques.
De façon alternative à la sélection d’un sous-modèle, d’autres techniques
permettent d’aborder l’estimation de modèles en grande dimension en ex-
ploitant la parcimonie. Le chapitre suivant traitera de la réduction de di-
mension et des estimations sous contraintes.

17
Statistique en grande dimension F. Lavancier

3.1 Erreurs de prévision

Pour fixer les idées on se place dans le cadre d’un modèle de régression
linéaire, mais les notions développées sont valables pour la plupart des mod-
èles considérés en statistique.
On suppose que le vrai modèle (inconnu) expliquant Y est

Y = X ∗β ∗ + 

où X ∗ est une matrice non aléatoire de taille (n, p∗ ) contenant p∗ variables

∗
explicatives, β ∗ ∈ Rp est le paramètre inconnu et  ∈ Rn est l’erreur de
modélisation, aléatoire, supposée centrée et de variance σ 2 In .
On ne connait pas a priori les "vraies" variables explicatives et on effectue
donc la régression linéaire de Y sur p variables explicatives regroupées dans
la matrice X de taille (n, p). En toute généralité, p n’est pas nécessairement
égale à p∗ , et même lorsque c’est le cas les variables explicatives dans X et
X ∗ ne sont pas nécessairement les mêmes. On obtient donc

Ŷ = P[X] Y = X β̂.

Il est à noter que β̂ n’a pas nécessairement la même taille que β ∗ . Pour évaluer
la qualité d’estimation on s’intéresse donc à X β̂ − X ∗ β ∗ dont les deux termes
ont la même dimension. Mais dans une optique de prévision on s’intéresse
plutôt aux résidus Y − Ŷ . Plusieurs erreurs quadratiques moyennes peuvent
être considérées.

1. L’erreur sur l’échantillon d’apprentissage (la MSE usuelle) est

1 1
Rtrain = E(kY − Ŷ k2 ) = E(kY − X β̂k2 ).
n n
Cette erreur nous informe de la qualité de l’ajustement sur les données
utilisées, mais ne nous dit rien sur ses qualités de prévision sur de
nouvelles données (échantillon test). Or c’est ce dernier point qui est
généralement intéressant.

2. L’erreur test "in" est l’erreur de prévision sur un nouvel échantillon

construit ainsi : on considère les mêmes valeurs des variables explica-
tives que dans l’échantillon d’apprentissage, mais un nouveau vecteur
Ỹ , indépendant de Y et de même loi, est associé. Autrement dit,

18
Statistique en grande dimension F. Lavancier

Ỹ = X ∗ β ∗ + ˜ où ˜ est indépendant et de même loi que . Comme

précédemment, on ne connait pas X ∗ et la prévision se base sur les
mêmes variables explicatives X que dans l’échantillon d’apprentissage.
La prévision de Ỹ est donc la même que celle de Y , i.e. X β̂, la différence
importante étant que β̂, construit à partir de Y , est indépendant de Ỹ .
L’erreur "in" est donc
1
Rin = E(kỸ − X β̂k2 ),
n
où Ỹ est indépendant et de même loi que Y .
3. Enfin l’erreur de prévision la plus générale, appelée erreur test, s’intéresse
à la qualité de prévision en présence à la fois de nouvelles valeurs X̃ ∗
et X̃ pour les variables explicatives, et d’une réalisation Ỹ indépen-
dante de Y . Les matrices X̃ ∗ et X̃, de taille respective (n, p∗ ) et (n, p),
contiennent les mêmes variables explicatives que X ∗ et X mais pour n
valeurs différentes. Autrement dit Ỹ = X̃ ∗ β ∗ + ˜ où ˜ est indépendant
et de même loi que , et la prévision se base sur X̃. La prévision de Ỹ
est alors X̃ β̂ et l’erreur associée
1
Rtest = E(kỸ − X̃ β̂k2 ).
n
Cette erreur dépend de X̃ et X̃ ∗ , i.e. Rtest = Rtest (X̃, X̃ ∗ ), et on
remarque que Rin = Rtest (X, X ∗ ).
Proposition 3.1.1. On a les relations suivantes.
1 1
Rtest = σ 2 + kX̃ ∗ β ∗ − E(X̃ β̂))k2 + tr(V ar(X̃ β̂))
n n
n
1 X
= (σ 2 + Bi2 + Vi )
n i=1
0
où Bi = x̃∗i β ∗ − E(x̃0i β̂) et Vi = V ar(x̃0i β̂) sont respectivement un terme de
biais et un terme de variance pour la prévision de ỹi .
Si p ≤ n,
1 n−p 2 1 p
Rtrain = kP[X]⊥ X ∗ β ∗ k2 + σ , Rin = σ 2 + kP[X]⊥ X ∗ β ∗ k2 + σ 2
n n n n

19
Statistique en grande dimension F. Lavancier

Si p ≥ n,
1 1
Rtrain = kP ⊥ X ∗ β ∗ k2 , Rin = 2σ 2 + kP ⊥ X ∗ β ∗ k2 .
n [X] n [X]
Ainsi
min(n, p) 2
Rin = Rtrain + 2 σ . (3.1)
n
Ces résultats conduisent aux observations suivantes.

• Le risque le plus général Rtest est la somme de trois termes : une

erreur incompressible σ 2 , due à la présence du bruit ˜ dans les nouvelles
données, qui est imprévisible; un terme de biais dû au fait qu’on a
effectué la régression sur les mauvaises variables X au lieu de X̃ (ce
terme disparaît si le bon choix a été effectué); et un terme de variance dû
à l’estimateur β̂. Le meilleur modèle satisfait donc un bon compromis
biais-variance.

• Le risque Rin , qui est un cas particulier de Rtest lorsque X̃ = X et

X̃ ∗ = X ∗ , fait intervenir de la même manière un terme de biais et de
variance, qui sont rendus plus explicites dans ce contexte.

• En analysant Rin , on s’aperçoit qu’un gros modèle, dans le sens où

X ∗ β ∗ ∈ [X], a un biais nul, car alors P[X]⊥ X ∗ β ∗ = 0. En revanche, il
aura une forte variance car p est alors grand. Au contraire, un petit
modèle (p petit) aura une faible variance mais un biais potentiellement
élevé.

• Pour choisir un bon modèle, il conviendrait d’estimer Rin ou mieux

Rtest . Mais la quantité naturelle à laquelle on a accès est le risque
empirique SCR = kY − Ŷ k2 . Ce dernier est un estimateur sans biais
de nRtrain , i.e. E(SCR/n) = Rtrain par définition, et non des risques de
prévision plus réalistes Rin et Rtest .

• Le risque Rtrain n’a pas la même dynamique que Rin et Rtest : son terme
de biais disparaît bien lorsque le modèle est gros (X ∗ β ∗ ∈ [X]), mais
son terme de variance aussi (il décroit avec p). Ainsi, baser le choix du
modèle sur Rtrain conduirait à choisir le plus gros modèle possible, celui
qui interpole les données et implique SCR = 0.

20
Statistique en grande dimension F. Lavancier

• Grâce à la relation (3.1), on peut proposer un estimateur sans biais de

Rin : SCR/n + 2 min(n, p) σ̂ 2 /n, pourvu que σ̂ 2 soit un estimateur sans
biais de σ 2 . C’est la motivation à la base du Cp de Mallows, voir la
section suivante. Le risque Rtest est par contre trop général pour être
estimé convenablement à partir de SCR.
Remarque 3.1.2. L’égalité (3.1) est vraie pour toute matrice X fixée déter-
ministe et motive le Cp de Mallows. Néanmoins, une fois le modèle choisi
grâce au Cp de Mallows (ou par tout autre critère basé sur les données),
l’égalité (3.1) devient fausse, autrement dit l’estimation de Rin a posteriori
(c’est à dire après avoir sélectionné le modèle) en utilisant SCR/n + 2pσ̂ 2 /n
devient biaisé, même si σ̂ 2 est un estimateur sans biais. En effet, choisir
un modèle, c’est à dire X, par un critère de sélection, rend le choix de X
aléatoire et dépendant de β̂ (puisque les critères en dépendent). Ainsi dans
le calcul de Rin conduisant à (3.1), la propriété E(P[X]⊥ Y ) = P[X]⊥ E(Y ) de-
vient fausse si X a été choisie a posteriori. En fait, l’estimation de Rin a
posteriori en utilisant le Cp de Mallows sous estime Rin . Ce phénomène est
appelé biais de sélection.
Proof. On a

nRtrain = E(kP[X]⊥ Y k2 )
= E(kP[X]⊥ X ∗ β ∗ + P[X]⊥ k2 )
= kP[X]⊥ X ∗ β ∗ k2 + E(kP[X]⊥ k2 )
= kP[X]⊥ X ∗ β ∗ k2 + E((P[X]⊥ )0 (P[X]⊥ ))
= kP[X]⊥ X ∗ β ∗ k2 + E(tr((P[X]⊥ )(P[X]⊥ )0 ))
= kP[X]⊥ X ∗ β ∗ k2 + tr(P[X]⊥ E(0 )P[X]⊥ )
= kP[X]⊥ X ∗ β ∗ k2 + σ 2 tr(P[X]⊥ )
(
kP[X]⊥ X ∗ β ∗ k2 + (n − p)σ 2 si p ≤ n,
=
kP[X]⊥ X ∗ β ∗ k2 si p ≥ n.

nRtest = E(k(Ỹ − E(Ỹ )) + (E(Ỹ ) − E(X̃ β̂)) + (E(X̃ β̂) − X̃ β̂)k2 )

En développant le carré, on observe que les termes croisés s’annulent car par
indépendance de Ỹ et β̂

E[(Ỹ − E(Ỹ ))0 (E(X̃ β̂) − X̃ β̂)] = [E(Ỹ − E(Ỹ ))]0 [E(E(X̃ β̂) − X̃ β̂)] = 0

21
Statistique en grande dimension F. Lavancier

et les autres termes croisés sont le produit d’un terme déterministe par
l’espérance d’une variable centrée. Ainsi
nRtest = E(kỸ − E(Ỹ )k2 ) + kE(Ỹ ) − E(X̃ β̂))k2 + E(kE(X̃ β̂) − X̃ β̂k2 )
= nσ 2 + kX̃ ∗ β ∗ − E(X̃ β̂))k2 + tr(V ar(X̃ β̂)).

En choisissant X̃ = X et X̃ ∗ = X ∗ dans l’expression de Rtest , on obtient

Rin et on en déduit
nRin = nσ 2 + kX ∗ β ∗ − E(X β̂))k2 + tr(V ar(X β̂))
= nσ 2 + kE(Y ) − E(P[X] Y )k2 + tr(V ar(P[X] Y ))
0
= nσ 2 + kE(P[X]⊥ Y )k2 + tr(P[X] V ar(Y )P[X] )
= nσ 2 + kP[X]⊥ X ∗ β ∗ k2 + σ 2 tr(P[X] )
(
nσ 2 + kP[X]⊥ X ∗ β ∗ k2 + pσ 2 si p ≤ n,
=
nσ 2 + kP[X]⊥ X ∗ β ∗ k2 + nσ 2 si p ≥ n.

3.2 Critères usuels

3.2.1 Cp de Mallows
Pour choisir un bon modèle de prévision, l’idéal est de minimiser l’erreur
test Rtest . Cette erreur très générale est difficile à estimer. De façon plus
raisonnable, on peut estimer Rin qui correspond au cas particulier de Rtest
lorsqu’on considère les mêmes valeurs des variables explicatives X. D’après
(3.1), on utilise l’estimateur naturel
SCR min(n, p) 2
R̂in = +2 σ̂
n n
qui est un estimateur sans biais dès que σ̂ 2 est sans biais. Pour cette raison,
σ̂ 2 n’est pas calculé à partir du modèle testé (d’où est issue la SCR) mais
provient généralement du plus gros modèle disponible Xmax . Le point clé est
que ce dernier vérifie X ∗ β ∗ ∈ [Xmax ]. En supposant que ce modèle contient
pmax variables avec pmax < n, et en notant SCRmax la SCR correspondante,
cela donne
SCRmax
σ̂ 2 = .
n − pmax

22
Statistique en grande dimension F. Lavancier

Cet estimateur est bien sans biais car E(σ̂ 2 ) = nRtrain /(n − pmax ) et d’après
la proposition 3.1.1, nRtrain = kP[Xmax ]⊥ X ∗ β ∗ k2 + (n − pmax )σ 2 = (n − pmax )σ 2
car X ∗ β ∗ ∈ [Xmax ].
Dans le cas pmax < n (et donc p < n), minimiser R̂in revient à minimiser

SCR
Cp = − n + 2p
σ̂ 2
qui est l’expression usuelle du Cp de Mallows implémenté dans les logiciels.
Si pmax ≥ n, l’estimation de σ 2 à partir du plus gros modèle ne peut pas se
faire comme précédemment, car alors SCRmax = 0. Le Cp n’est pas défini
dans ce cas.
Le défaut principal du Cp de Mallows est la nécessité de trouver un es-
timateur sans biais de σ 2 . L’approche précédente requiert pmax < n, ce qui
empêche le calcul du critère en grande dimension pmax ≥ n, même si le mod-
èle testé ne contient qu’un nombre p faible de variables. De plus, si pmax < n
mais pmax est grand, l’estimateur σ̂ 2 sera certes sans biais, mais sa variance
sera élevée ce qui impliquera une mauvaise estimation de Rin et donc un
mauvais choix de modèle. Le Cp de Mallows est donc à utiliser uniquement
en présence d’un nombre total modéré de variables explicatives.

3.2.2 Critère AIC

L’AIC (Akaike Infomation Criterion) est un critère construit à partir de la
vraisemblance. Supposons que le modèle soit Gaussien ( suit une loi Gaussi-
enne). Pour un modèle donné faisant intervenir comme variables explicatives
la matrice X, l’estimateur du maximum de vraisemblance maximise la log-
vraisemblance log VX (Y ; θ) où VX désigne la vraisemblance du modèle con-
sidéré et θ = (β, σ 2 ) est le vecteur des paramètres inconnus du modèle. On
obtient comme solution dans le cas Gaussien l’estimateur des MCO β̂ et
σ̂ 2 = SCR/n. Ainsi − log VX (Y ; θ̂) est minimal pour ce modèle.
Pour mesurer la qualité du modèle, l’idée est de considérer un échantillon
test indépendant Ỹ et d’observer ce que vaut en moyenne − log VX (Ỹ ; θ̂), où
θ̂ est l’estimateur calculé sur l’échantillon d’apprentissage Y . C’est une façon
de mesurer la qualité de l’estimation sur un échantillon test. On retiendra
finalement le modèle (c’est à dire la matrice X) pour lequel E(− log VX (Ỹ ; θ̂))
est minimal. Il s’agira en pratique d’estimer E(− log VX (Ỹ ; θ̂)), ce qui est
l’objectif du critère AIC.

23
Statistique en grande dimension F. Lavancier

Cette approche est dans le même esprit que la minimisation de Rin ayant
conduit au Cp de Mallows, mais ici c’est E(− log VX (Ỹ ; θ̂)) qui joue le rôle
de la fonction de perte et non Rin . En fait lorsque σ 2 est supposé connu
(hypothèse peu réaliste), le critère AIC et le Cp de Mallows coincident.

Proposition 3.2.1. Dans le modèle de régression Gaussien, si σ 2 est connu,

alors le critère AIC est équivalent au Cp de Mallows.

Proof. Puisque σ 2 est connu, le seul √ paramètre inconnu est β. Dans le cas
Gaussien − log VX (Ỹ ; β̂) = n log(σ 2π) + kỸ − X β̂k /(2σ 2 ) de telle sorte
2

que le modèle qui minimise E(− log VX (Ỹ ; θ̂)) minimise de façon équivalente
Rin . Le critère AIC et le critère Cp coincident : il suffit de minimiser un
estimateur sans biais de Rin , c’est à dire SCR/n + 2 min(n, p)σ 2 /n, où ici σ 2
n’a pas à être estimé puisqu’il est supposé connu.
Dans le cas général (σ 2 inconnu), H. Akaike a montré dans un article de
1973 que sous certaines conditions, lorsque n est grand,
2 p+1
E(− log VX (Ỹ ; θ̂)) ≈ − E(− log VX (Y ; θ̂)) + 2 , (3.2)
n n
où p + 1 correspond à la dimension du vecteur θ = (β, σ 2 ). Le critère AIC
consiste donc à minimiser −2 log VX (Y ; θ̂) + 2(p + 1) qui est (à un facteur n
près qui ne change rien à la minimisation) un estimateur asymptotiquement
sans biais du risque. Dans le cas Gaussien,

n n kY − X β̂k2
− log VX (Y ; θ̂) =log(2π) + log(σ̂ 2 ) +
2 2 2σ̂ 2
n n SCR
= log(2π) + log(SCR/n) +
2 2 2SCR/n

de telle sorte que le critère AIC revient à minimiser

SCR
AIC = n log + 2(p + 1).
n
Comme cela a été remarqué dans le chapitre précédent, lorsque p ≥ n
SCR = 0 donc AIC = −∞. Le critère AIC n’est donc pas adapté à la
grande dimension.

24
Statistique en grande dimension F. Lavancier

3.2.3 Critère BIC

Le critère BIC (Bayesian Information Criterion) est motivé différemment
que par la minimisation de l’erreur test. Comme son nom l’indique, il s’agit
d’une approche bayésienne. Dans ce contexte on associe à chaque modèle une
probabilité a priori et on choisit le modèle ayant la probabilité a posteriori
(sachant les observations) maximale.
Formellement, on note

• MX le modèle linéaire associée à la matrice X

• P(MX ) : la probabilité a priori que ce modèle soit valide, c’est à dire

que X contiennent les bonnes variables explicatives, sans considération
de la valeur du paramètre associé. On supposera par la suite que cette
probabilité est identique pour chaque modèle, de telle sorte qu’aucun
modèle n’est privilégié a priori.

• P(θ|MX ) : la loi a priori du paramètre θ dans le modèle MX .

• P(Y |MX , θ) : la loi de Y sachant que le modèle est MX associé au

paramètre θ. Vu comme une fonction de θ, il s’agit simplement de la
vraisemblance, notée VX (Y ; θ) dans la partie précédente.

La probabilité a posteriori de MX sachant les observations Y vaut d’après

la formule de Bayes

P(Y |MX )P(MX )

P(MX |Y ) = ,
P(Y )

où P(Y ) est la loi de Y intégrée sur tous les modèles possibles. Puisqu’on a
supposé P(MX ) constant, maximiser P(MX |Y ) en MX revient à maximiser
P(Y |MX ). Cette quantité vaut d’après la formule des probabilités totales
Z
P(Y |MX ) = P(Y |MX , θ)P(θ|MX )dθ.

Lorsque n est grand, et sous certaines conditions, G. Schwarz (1978) a montré

que l’on peut approcher le logarithme de l’intégrale précédente de la façon
suivante
p+1
log P(Y |MX ) ≈ log P(Y |MX , θ̂) − log n,
2

25
Statistique en grande dimension F. Lavancier

où θ̂ est l’estimateur du maximum de vraisemblance de θ dans le modèle

MX et (p + 1) est sa dimension. Un résultat remarquable est que cette
approximation ne dépend pas du choix de la loi a priori P(θ|MX ) sur les
paramètres.
Le critère BIC est défini comme valant −2 fois la quantité précédente, qu’il
s’agit donc de minimiser parmi les modèles testés. On remarque que la forme
est similaire au critère AIC puisque c’est la somme de la log-vraisemblance
calculée en θ̂ plus un terme de pénalité, qui fait intervenir ici log n au lieu de
2 dans AIC.
Si l’on suppose enfin que le modèle est Gaussien, la log-vraisemblance
dans cette formule se simplifie comme pour le AIC et on obtient
SCR
BIC = n log + (p + 1) log n.
n
Pour les mêmes raisons que pour AIC, le critère BIC n’est pas adapté à la
grande dimension.

3.2.4 Comparaison des critères

Tous les critères précédents consistent à minimiser une expression de la forme

f (SCR) + c(n)p

où f est une fonction croissante de SCR et c(n)p est un terme pénalisant

les modèles ayant beaucoup de variables. La fonction c(n) est appelée la
pénalité. Elle vaut 2 pour Cp et AIC et elle vaut log(n) pour BIC. D’autres
pénalités ont été proposées, privilégiant de façon plus ou moins forte les
modèles les plus parcimonieux (p faible).
On observe en particulier que dès que log(n) > 2, le critère BIC pénalise
davantage la complexité du modèle que AIC. Le Cp de Mallows est quant à
lui très similaire à AIC, comme l’atteste la proposition 3.2.1.
Lors de la sélection de variables dans un modèle de régression linéaire,
les critères usuelles s’ordonnent de la manière suivante en fonction de leur
propension à sélectionner le modèle le plus parcimonieux :

BIC < F test < Cp ≈ AIC < Ra2 < R2

26
Statistique en grande dimension F. Lavancier

3.3 Validation croisée

Le principe de la validation croisée est d’estimer le risque test Rtest en
confrontant le modèle à un échantillon test qui n’a pas été utilisé pour
l’ajustement du modèle. Il y a différentes manières de construire un échan-
tillon test, décrites ci-après.
Les avantages de la validation croisée par rapport à l’utilisation des critères
de sélection est double :

• La validation croisée estime Rtest (X̃) pour de nouvelles valeurs X̃ et

pas seulement Rin = Rtest (X). (L’estimation est néanmoins limitée
aux valeurs X̃ présentes dans l’échantillon test et il faut donc veiller à
ne pas trop généraliser son interprétation).

• La procédure est plus générale puisque l’estimation de Rtest ne s’appuie

pas sur l’hypothèse que les données suivent un modèle particulier (con-
trairement par exemple à AIC ou BIC qui s’appuient pour leur con-
struction sur la vraisemblance du modèle).

Le principal inconvénient de la validation croisée est son temps de calcul,

qui peut parfois être prohibitif.
L’objectif est d’estimer Rtest . Comme toute estimation, les procédures
suivantes souffrent d’un biais et d’une variance, que l’on souhaite les plus
petits possibles. Un biais trop élevé signifie que l’on sous-estime ou sur-estime
le risque associé. Par ailleurs, une procédure d’estimation trop variable rend
peu fiable la comparaison des risques entre deux modèles : le fait qu’un risque
estimé est inférieur à un autre peut être dû aux fluctuations des estimateurs.
Pour choisir la bonne méthode de validation croisée, on veillera à ce que le
biais et la variance soient les moins élevés possibles, tout en respectant un
temps de calcul raisonnable.

3.3.1 Hold out

Le principe est simple : on sépare l’échantillon en deux parties, une pour
estimer le modèle, l’autre pour estimer l’erreur de prévision. Concernant
la taille de chaque échantillon, on rencontre deux choix courants : soit la
même taille n/2 pour les deux échantillons, soit 3n/4 pour l’échantillon
d’apprentissage et n/4 pour l’échantillon test.
Avantages :

27
Statistique en grande dimension F. Lavancier

• Très simple à mettre en oeuvre

• Un seul ajustement à effectuer par modèle.
Inconvénients :
• La taille de l’échantillon d’apprentissage (n/2 ou 3n/4) rend l’estimation
moins précise que si on utilisait les n observations. Cela conduit à une
sur-estimation de l’erreur test par rapport à celle issue d’une estimation
sur n individus.
• L’évaluation sur un seul échantillon test rend l’estimation de l’erreur
très variable (la variance est élevée).
Pour minimiser ce dernier inconvénient, on peut recommencer la procé-
dure Hold-out plusieurs fois (une dizaine de fois), en considérant plusieurs
séparations aléatoires de l’échantillon, toujours de même taille. Le risque
final, obtenu par la moyenne des risques issus de chaque découpage, est alors
moins variable. Cette généralisation fait néanmoins perdre les avantages ini-
tiaux, en particulier le coût est supérieur puisqu’il faut estimer et évaluer
chaque modèle autant de fois qu’il y a de découpages différents.
L’idée d’itérer les découpages est en fait à la base des méthodes suivantes.

3.3.2 Leave-one-out (LOO)

La validation croisée Leave-one-out consiste à mettre de côté une observation,
d’estimer le modèle sur les n−1 autres observations, puis de calculer le risque
de prévision sur l’observation mise de côté. Ce processus est répété pour
chaque observation de l’échantillon, de telle sorte que n erreurs de prévision
sont finalement calculées. Leur moyenne est une estimation de l’erreur test.
(−i)
De façon plus formelle, soit ŷi la prévision de yi à partir du modèle
estimé sur toutes les observations sauf i. Le risque estimé par LOO est
n
1X (−i)
CVLOO = (yi − ŷi )2 .
n i=1

Avantages :
• Chaque échantillon d’apprentissage a une taille n − 1, très proche de
la vraie taille n de l’échantillon total. Ainsi le risque estimé sera peu
biaisé.

28
Statistique en grande dimension F. Lavancier

Inconvénients :

• Peut-être très couteux à mettre en oeuvre car il nécessite en toute

généralité n ajustements de chaque modèle testé. Il y a une exception
notable : dans le cas de la régression linéaire, CVLOO s’exprime directe-
ment grâce à la seule estimation du modèle sur les n observations, cf
la proposition 3.3.1.

• Tous les échantillons d’apprentissage sont très semblables (ils différent

que par une observation) et donc le risque final par LOO ne moyenne
pas suffisamment de situations différentes. En conséquences l’estimation
du risque a une forte variance.

Proposition 3.3.1. L’estimation de l’erreur test par LOO d’un modèle de

régression linéaire Y = Xβ +  est simplement
n  2
1X yi − ŷi
CVLOO = ,
n i=1 1 − hi

où ŷi est la prévision de yi à partir du modèle estimé sur toutes les obser-
vations et hi est l’effet levier, défini par le ième élément de la diagonale de
P[X] , c’est à dire hi = x0i (X 0 X)−1 xi .

Proof. On note Y(−i) , respectivement (−i) , le vecteur Y , resp. , privé de son

i-ème élément yi , resp. i , et X(−i) la matrice X privée de sa i-ème ligne x0i .
Le modèle sans l’observation i s’écrit donc Y(−i) = X(−i) β + (−i) . On note
enfin β̂(−i) l’estimation de β par les MCO dans ce modèle.
(−i)
Il s’agit de montrer que yi − ŷi = ˆi /(1 − hi ). Ceci est une conséquence
du lemme suivant.
Lemme 3.3.2.
1
β̂(−i) = β̂ − (X 0 X)−1 xi ˆi .
1 − hi
En effet, on a alors

(−i) x0i (X 0 X)−1 xi hi ˆi

yi − ŷi = yi − x0i β̂(−i) = yi − x0i β̂ + ˆi = ˆi + ˆi = .
1 − hi 1 − hi 1 − hi

29
Statistique en grande dimension F. Lavancier

Preuve du lemme. On sait que

0
β̂(−i) = (X(−i) X(−i) )−1 X(−i)
0
Y(−i) .
0
On peut facilement vérifier que X(−i) X(−i) = X 0 X − xi x0i et que X(−i)
0
Y(−i) =
0
X Y − xi yi . Par ailleurs la formule de Sherman-Morrison donne l’inverse de
M + uv 0 pour toute matrice inversible M de taille p et tous vecteurs u et v
dans Rp :
M −1 uv 0 M −1
(M + uv 0 )−1 = M −1 − .
1 + u0 M −1 v
On en déduit que

0 1
(X(−i) X(−i) )−1 = (X 0 X)−1 + (X 0 X)−1 xi x0i (X 0 X)−1 .
1 − hi
Ainsi
 
1
β̂(−i) = (X 0 X)−1 + 0 −1 0 0
(X X) xi xi (X X)−1
(X 0 Y − xi yi )
1 − hi
(X 0 X)−1 xi x0i β̂ (X 0 X)−1 xi hi yi
= (X 0 X)−1 X 0 Y − (X 0 X)−1 xi yi + −
1 − hi 1 − hi
(X 0 X)−1 xi
= β̂ − ((1 − hi )yi − x0i β̂ + hi yi )
1 − hi
1
= β̂ − (X 0 X)−1 xi ˆi .
1 − hi

3.3.3 Leave-k-out
L’idée est la même que le Leave-one-out, sauf qu’au lieu de mettre de côté une
seule observation, on en met k et on calcule le risque de prévision sur les k
observations mises de  côté.
 On réitère le processus pour tous les découpages
n
possibles, il y en a = n!/(k!(n − k)!), et on moyenne tous les risques
k
pour obtenir l’estimation finale.
Si k = 1, on retrouve évidemment le LOO. Pour k plus grand, l’estimation
du risque est un peu plus biaisé qu’avec le LOO (car on ajuste le modèle sur
moins d’observations), mais la variance est plus faible car les échantillons

30
Statistique en grande dimension F. Lavancier

d’apprentissage sont plus diverses. Ainsi il existe un choix de k garantissant

un bon compromis biais-variance (ce choix dépend de n et le consensus est
de le prendre de l’ordre de n/10 ou n/5). Cependant, le défaut majeur de
cette démarche est son coût : elle requiert d’estimer le modèle sur un nombre
prohibitif d’échantillons d’apprentissage. A titre d’exemple, pour n = 100
et k = 5, ce qui est loin de constituer un cas de grande dimension, cela
représente plus de 75 millions d’échantillons.
Avantage :

• Bonne estimation du risque test pour un choix de k approprié.

Inconvénient :

• Généralement impossible à mettre en oeuvre en pratique.

3.3.4 K-fold
On coupe l’échantillon en K parties égales (ou à peu près si n n’est pas
multiple de K) de taille respective n/K. On apprend le modèle en utilisant
K − 1 parties de l’échantillon et on le teste sur la partie restante pour fournir
une estimation du risque. En répétant cette opération pour tous les cas
possibles, cela représente K risques à calculer, que l’on moyenne pour fournir
l’estimation du risque final.
Dans cette procédure, les K échantillons d’apprentissage ont donc une
taille de n − n/K et le risque final est la moyenne des K risques tests, chacun
calculé sur un échantillon test de n/K individus.
Lorsque K = n, il s’agit de la procédure LOO.
Lorsque K = 2, cela ressemble à une procédure Hold out dans laquelle
l’échantillon d’apprentissage et l’échantillon test ont même taille. Ce n’est
cependant pas tout à fait équivalent car ici on inverse le rôle joué par chaque
partie pour fournir deux estimations du risque que l’on moyenne, contraire-
ment à Hold out où une seule estimation du risque est effectuée.
Si K est choisi trop grand, il y aura un faible biais car les échantil-
lons d’apprentissage auront une taille très proche de n, mais par contre
l’estimation du risque souffrira d’une grande variance comme pour le LOO.
Le désavantage majeure dans ce cas est similaire au LOO : le calcul peut-être
très long car il faut estimer K modèles. A l’inverse, si K est petit, le calcul
est rapide, la variance d’estimation est moindre car on moyenne des situa-
tions où les échantillons sont très différents, mais l’estimation risque d’être

31
Statistique en grande dimension F. Lavancier

biaisée car la taille de l’échantillon d’apprentissage est beaucoup plus petite

que la vraie taille n.
En pratique, les valeurs K = 5 à K = 10 sont recommandées.
Avantage :

• Bonne estimation du risque test pour un choix de K approprié (K = 5

à K = 10)

• Relativement rapide à mettre en oeuvre pour les choix de K précédents.

Inconvénient :

• L’estimation du risque n’est pas aussi bonne qu’avec une procédure

leave-k-out optimale.

3.4 Sélection automatique d’un sous-modèle

L’idée est de parcourir un ensemble de sous-modèles, d’estimer pour chacun
d’entre eux le risque de prévision associé à l’aide d’une des méthodes dis-
cutées précédemment, puis de retenir le sous-modèle dont l’erreur estimée
est minimale. Il existe quatre façons classiques de parcourir un ensemble de
sous-modèles :

• La recherche exhaustive. On parcourt tous les 2p sous-modèles possi-

bles. Cette approche est impossible en grande dimension. Par exemple,
pour p = 20 (ce qui reste raisonnable), cela représente plus d’un million
de modèles à tester.

• La sélection stepwise backward. On part du plus gros modèle possible

à p variables et on élimine la variable la moins importante, dans le sens
où le modèle sans cette dernière a le risque le plus faible parmi les mod-
èles à p − 1 variables. On continue jusqu’à ce que l’élimination d’une
variable n’améliore plus le risque. Il y a au maximum 1 + p(p + 1)/2
modèles testés. Par exemple, pour p = 20, cela représente 211 mod-
èles, ce qui est beaucoup plus raisonnable que la recherche exhaustive.
Par contre cette démarche ne garantit pas de trouver le sous-modèle
optimal. En effet une variable éliminée en début d’algorithme n’a plus
aucune chance d’être réintégrée par la suite, alors qu’il est possible
qu’elle appartienne en réalité au sous-modèle optimal. Cet algorithme

32
Statistique en grande dimension F. Lavancier

est en général inutilisable en grande dimension car les plus gros modèles
sont trop coûteux à estimer, et les critères de sélection lorsque p > n
deviennent dégénérés (R2 = 1, AIC = −∞, etc).

• La sélection stepwise forward. On part du plus petit modèle ne con-

tenant que la constante et on ajoute la variable la plus importante. On
continue jusqu’à ce que l’ajout d’une variable n’améliore plus le modèle.
Comme pour la sélection backward, il y a au maximum 1 + p(p + 1)/2
modèles testés, sans garantie de trouver le sous-modèle optimal (une
variable introduite en début d’algorithme y reste jusqu’à la fin). Con-
trairement à la sélection backward, la sélection forward est générale-
ment possible en grande dimension car on part des plus petits modèles
dont l’estimation est faisable.

• La sélection stepwise hybride (backward ou forward). Pour commencer

on part du plus petit modèle (pour la hybride forward) et on ajoute
la variable la plus importante. Dans les étapes suivantes, on estime
l’erreur associée à l’ajout de chaque variable (comme en forward clas-
sique) mais on estime également l’erreur associée au retrait de chaque
variable déjà présente dans le modèle. Ainsi le modèle suivant retenu
peut contenir une variable de plus ou une variable de moins. On con-
tinue ainsi jusqu’à ce que l’ajout ou le retrait d’une variable n’améliore
plus la modélisation. La hybride backward est similaire sauf que l’on
part du plus gros modèle. Le coût est un peu plus élevé que les
sélections backward et forward car à chaque étape les p variables sont
testées. Mais cette méthode parcourt davantage de modèles.

En présence d’un grand nombre de variables, la méthode hybride forward est

la plus recommandée, car elle part des plus petits modèles dont l’estimation
est faisable et elle a un coût raisonnable tout en parcourant une grande variété
de sous-modèles possibles.

Mise en oeuvre sous R (exemple avec lm et les critères AIC et BIC) :

En supposant que le tableau de données se nomme tab, on peut utiliser
la fonction regsubsets de la librairie leaps pour effectuer une sélection
exhaustive ou backward ou forward dans un modèle de régression linéaire,
de la façon suivante.
res=regsubsets(y~.,data=tab)
res=regsubsets(y~.,data=tab,method="backward")

33
Statistique en grande dimension F. Lavancier

res=regsubsets(y~.,data=tab,method="forward")
La valeur des SCR, BIC, Cp pour chaque modèle testé est disponible dans
summary(res), voir aussi plot(res,scale="Cp") (par exemple).
De façon alternative, on peut utiliser la fonction step qui s’applique à
des modèles plus généraux (issus de la fonction glm notamment). Pour cela
on commence par estimer le plus gros modèle et le plus petit modèle.
fit_full=lm(y~.,data=tab)
fit0=lm(y~1,data=tab)
Pour la sélection backward, en supposant qu’il y a n observations :
step(fit_full,direction="backward") #critère AIC
step(fit_full,direction="backward",k=log(n)) #critère BIC
Pour la sélection backward hybride (par défaut sous R) :
step(fit_full) #critère AIC
step(fit_full,k=log(n)) #critère BIC
Pour la sélection forward :
step(fit0,scope=formula(fit_full),direction="forward") #critère AIC
step(fit0,scope=formula(fit_full),direction="forward",k=log(n)) #critère BIC
Pour la sélection forward hybride :
step(fit0,scope=formula(fit_full)) #critère AIC
step(fit0,scope=formula(fit_full),k=log(n)) #critère BIC

3.5 Estimation du modèle retenu

Une fois le modèle sélectionné par l’une des approches précédentes, il convient
de l’estimer sur les n individus de l’échantillon total avant de l’utiliser pour
des prévisions.
Cette précision concerne surtout un modèle retenu par validation croisée
: son estimation s’est faite sur un sous-échantillon d’apprentissage afin de le
confronter à un sous-échantillon test. Maintenant que ce modèle est retenu, il
faut l’estimer sur tout l’échantillon à disposition, ce qui ne peut qu’améliorer
la qualité d’estimation et donc les prévisions associées.

34
Chapter 4

Réduction de dimension et
régression sous contraintes

En présence de nombreuses variables explicatives, on suppose généralement

que peu d’entre elles sont pertinentes pour modéliser Y . On peut relever
trois familles de méthodes qui permettent d’exploiter cette parcimonie:

• Les méthodes de sélection, abordées dans le chapitre précédent. L’idée

est de choisir le sous-modèle dont l’estimation du risque de prévision
est minimale.

• La réduction de dimension. L’idée est de projeter les p variables, évolu-

ant en toute généralité dans un espace vectoriel de dimension p, dans un
sous-espace vectoriel de dimension beaucoup plus petit. On peut alors
régresser Y sur ce sous-espace afin d’éviter les écueils de la grande di-
mension. C’est l’objet de la section 4.1 qui présente les méthodes PCR
et PLS.

• Les estimations contraintes. L’idée est d’utiliser une méthode d’estimation

qui contraint les paramètres à ne pas exploser (contrairement aux MCO
en grande dimension). Ainsi l’estimation est moins variable et les prévi-
sions plus fiables. Parmi ces méthodes, les techniques de type "Lasso"
conduisent à estimer certains coefficients par 0, auquel cas une sélection
de variables s’opère dans le même temps. Ces méthodes sont abordées
dans les sections 4.2-4.4.

35
Statistique en grande dimension F. Lavancier

4.1 Réduction de dimension : PCR et PLS

4.1.1 Régression sur composantes principales (PCR)
L’analyse en composantes principales (ACP ou PCA en anglais) consiste
à trouver une base orthogonale de l’espace vectoriel [X] dont les vecteurs
Z1 , . . . , Zp , appelés composantes principales ou axes principaux, sont con-
struits de telle sorte à garder le plus "d’information" possible contenue dans
les vecteurs initiaux X1 , . . . , Xp . L’information est quantifiée à l’aide de la
variance des coordonnées des n individus sur chaque nouvel axe : plus la
variance est élevée et plus le nuage de points se répartit bien sur l’axe, gar-
dant ainsi le maximum de la diversité initiale du nuage, contrairement au
cas extrême inverse où tous les points seraient projetés au même endroit
(conduisant à une variance nulle).
Un exemple de construction des 2 axes principaux d’un nuage de points
en 2D est donné dans la figure 4.1. L’intérêt de considérer des projections
sur les premiers axes principaux pour résumer un nuage de points en illustré
dans la figure 4.2. Les données sont l’image 3D d’un chameau. Le premier
axe principal est celui qui traverse le chameau sur sa longueur, le second sur
sa hauteur et le troisième sur sa largeur. La projection sur le plan formé
par les deux premiers axes (à droite de la figure) suffit à bien représenter les
données. Cette projection est plus informative que celle sur les deux derniers
axes (à gauche de la figure), à partir de laquelle il est plus difficile de deviner
la nature des données initiales. Travailler à partir du plan factoriel de droite
permet donc une réduction de dimension (de 3 à 2 dans cet exemple) sans
grande perte d’information.
Cette approche est sensible à l’échelle : si une variable Xj est d’un ordre
de grandeur 1000 fois supérieur aux autres variables, l’essentiel de la variance
du nuage de points reposera sur Xj et le premier axe sera donc très proche
de Xj . Pour éviter ce phénomène, les variables initiales sont centrées et
réduites. On suppose donc dans la suite que chaque variable Xj est centrée
et de variance 1.
Concrètement, étant donnée une matrice X de variables centrées ré-
duites, le premier axe Z1 est choisi comme étant la combinaison linéaire
de X1 , . . . , Xp de variance maximale : Z1 = Xα1 avec kα1 k = 1 et Var(Xα1 )
maximale parmi tous les vecteurs de la forme Xα. Dans ce formalisme,
α1 ∈ Rp représente la direction du premier axe principal et Xα1 ∈ Rn est
l’ensemble des coordonnées du nuage de points sur cet axe. Le second axe

36
Statistique en grande dimension F. Lavancier

Figure 4.1: Construction des deux axes principaux d’un nuage de points en
2D.

37
Statistique en grande dimension F. Lavancier

Figure 4.2: Chameau (en 3D) observé via sa projection sur le plan formé
des axes principaux 2 et 3 de l’ACP (à gauche) et sur le plan formé des
deux premiers axes de l’ACP (à droite). La projection 2D de droite est plus
informative et suffit à comprendre l’objet 3D.

est choisi de la même manière, avec la contrainte supplémentaire d’être or-

thogonal à Z1 . Cette dernière contrainte s’écrit à l’aide du produit scalaire
(Xα2 )0 Xα1 = 0, soit α20 X 0 Xα1 = 0. De façon générale, pour j = 1, . . . , p, le
j-ème axe est Zj = Xαj où

αj = argmax Var(Xα) (4.1)

α∈Rp

sous les contraintes kαk = 1 et α0 X 0 Xαl = 0 pour tout l = 1, . . . , j − 1.

Dans (4.1), puisque les variables sont supposées centrées, on a Var(Xα) =
1 0 0
n
α X Xα.
Par construction tous les axes sont orthogonaux et ils sont ordonnés du
plus "informatif" Z1 , au moins informatif Zp , au sens où la variance des
coordonnées des individus sur les axes décroit. L’algorithme pour résoudre
le problème d’optimisation précédent est simple : les directions αj corre-
spondent aux vecteurs propres (normalisés) de la matrice X 0 X, ordonnés par
ordre décroissant de leur valeur propre associée.

Le principe de la régression sur composantes principales (PCR) est le

suivant : plutôt que de régresser Y sur X1 , . . . , Xp , on régresse Y sur les

38
Statistique en grande dimension F. Lavancier

premiers axes principaux Z1 , . . . , ZM (en incluant de plus une constante),

où le nombre de composantes M retenues est à choisir. Cette approche
permet de réduire la dimension du problème de p à M , tout en garantissant
l’orthogonalité des variables explicatives retenues. Le modèle s’écrit ainsi
M
X
Y = γ0 + γj Zj + 
j=1

et l’on peut montrer que les estimateurs par moindres carrés associés sont
Y 0 Zj
γ̂0 = Ȳ , γ̂j = . (4.2)
Zj0 Zj

L’interprétation des M composantes principales retenues et des coeffi-

cients estimés γ̂j de la PCR n’est pas toujours aisée. Néanmoins, il est
possible de revenir aux variables initiales puisque Zj = Xαj . On obtient

Ŷ = Ȳ + X β̂ (4.3)
PM
avec β̂ = j=1 γ̂j αj .

La mise en oeuvre d’une PCR requiert le choix du nombre M de com-

posantes retenues. Ceci se fait par validation croisée (voir les différentes
versions dans le chapitre précédent) : on estime l’erreur de prévision de la
PCR lorsque M varie de 1 à p (ou une valeur maximale inférieure à p), et on
retient la valeur de M associée à l’erreur minimale.

Mise en oeuvre sous R : fonction pcr de la libraire pls. La sélection par

validation croisée est proposée en option de la fonction.

4.1.2 Régression des moindres carrés partiels (PLS)

Les composantes d’une ACP sont construites pour représenter au mieux
l’information contenue dans les variables X1 , . . . , Xp . Cette construction ne
tient pas compte du lien avec Y , qui est pourtant l’objectif premier d’un
modèle de régression. La régression PLS remédie à ce défaut en constru-
isant les composantes de telle sorte qu’elles représentent au mieux les vari-
ables X1 , . . . , Xp et qu’elles soient dans le même temps le plus possible cor-
rélées avec Y . Ainsi, contrairement à (4.1), les composantes d’une PLS sont

39
Statistique en grande dimension F. Lavancier

Z1 , . . . , Zp où pour tout j, Zj = Xαj et

αj = argmax Cov(Xα, Y ) (4.4)
α∈Rp

sous les contraintes kαk = 1 et αX 0 Xαl = 0 pour tout l = 1, . . . , j − 1. Dans

cette expression Cov(Xα, Y ) = α0 X 0 Y /n car les variables sont supposées
centrées. Puisque Cov(Xα, Y )2 = Corr(Xα, Y )2 Var(Xα)Var(Y ), où Corr
désigne la corrélation, (4.4) peut aussi se récrire
αj = argmax Corr(Xα, Y )2 Var(Xα),
α∈Rp

ce qui montre bien que chaque composante est choisie comme le compromis
qui d’une part maximise la variance de la projection du nuage de points formé
par les variables explicatives X1 , . . . , Xp (comme pour l’ACP) et d’autre part
maximise la corrélation entre la composante et Y .
L’algorithme pour construire les axes de la PLS est le suivant : partant
de la matrice X(1) = X centrée réduite et de j = 1, tant que j ≤ p et αj 6= 0
0 0
1. αj = X(j) Y /kX(j) Y k.
2. Zj = X(j) αj
3. X(j+1) ← X(j) − Zj Zj0 X(j) /(Zj0 Zj ) et retour à 1 avec j ← j + 1.
Si αj = 0 pour un certain j, Zk = 0 pour tout k ≥ j. La première étape est
la solution de (4.4) sachant que la contrainte d’orthogonalité (αX 0 Xαl = 0
pour tout l = 1, . . . , j − 1) est nécessairement vérifiée grâce à l’étape 3.
En effet cette opération effectue la projection P[Zj ]⊥ X(j) , autrement dit elle
construit de nouvelles variables (regroupées dans la nouvelle matrice X(j+1) )
orthogonales à l’axe Zj . Ainsi tout vecteur X(j+1) α considéré dans l’itération
suivante sera nécessairement orthogonal aux axes déjà construits.

Une fois les axes obtenus, la régression PLS est similaire à la PCR : on
effectue la régression sur les M premiers axes, en choisissant M par valida-
tion croisée. De part l’orthogonalité des axes construits, les coefficients de
régression s’obtiennent de la même manière que pour la PCR, voir (4.2). De
même, puisque chaque axe est combinaison linéaire des variables X1 , . . . , Xp ,
on peut exprimer Ŷ en fonction de X, comme dans (4.3) pour la PCR.

Mise en oeuvre sous R : fonction plsr de la libraire pls. La sélection par

validation croisée est proposée en option de la fonction.

40
Statistique en grande dimension F. Lavancier

4.2 Régression ridge

En présence d’un grand nombre de variables, la matrice X 0 X n’est plus
forcément de rang p et n’est donc plus inversible, rendant l’estimateur des
moindres carrés incalculable. De même, en présence de variables fortement
corrélées entre elles (problème de multicolinéarité), la matrice peut être in-
versible mais son inverse, et par suite l’estimateur des MCO, sont très in-
stables, dans le sens où une légère modification des données (par exemple
la suppression d’un individu, c’est à dire d’une ligne de X) peut conduire
à une matrice inverse radicalement différente, et de même pour la valeur
de l’estimateur. Ce phénomène est clairement reflété dans la variance de
l’estimateur par MCO qui vaut σ 2 (X 0 X)−1 .
Lorsque X 0 X est non inversible ou d’inverse instable, cela se lit dans ses
valeurs propres : certaines sont nulles ou quasiment nulles. Une manière de
régulariser l’inversion consiste à ajouter une petite valeur λ > 0 aux valeurs
propres. On remplace donc X 0 X par X 0 X + λIp où Ip est la matrice identité
de taille p. Il s’agit de la régularisation de Tikhonov, connue en statistique
sous le nom de régression ridge.
L’estimateur ridge de β dans un modèle de régression est ainsi
β̂ridge = (X 0 X + λIp )−1 X 0 Y.
La proposition suivante montre que l’introduction du paramètre λ biaise
l’estimation, mais en contre-partie diminue sa variance (c’est l’effet de la
stabilisation), de telle sorte que pour des valeurs petites de λ, l’estimateur
ridge est préférable à l’estimateur par MCO au sens du coût quadratique.
Proposition 4.2.1. Si Y = Xβ +  avec E() = 0 et Var() = σ 2 In :
• l’estimateur ridge est biaisé, précisément E(β̂ridge ) = (X 0 X+λIp )−1 X 0 Xβ

• sa variance vaut Var(β̂ridge ) = σ 2 (X 0 X + λIp )−1 X 0 X(X 0 X + λIp )−1

• sa matrice EQM est inférieure à celle de l’estimateur par MCO (dans
le sens où la différence est définie positive) dès que 0 < λ < 2σ 2 /kβk2 .
Proof. Le calcul de l’espérance et de la variance est immédiat. On en déduit
après quelques calculs
EQM (β̂ridge ) = (E(β̂ridge ) − β)(E(β̂ridge ) − β)0 + Var(β̂ridge )
= (X 0 X + λIp )−1 (σ 2 X 0 X + λ2 ββ 0 )(X 0 X + λIp )−1 .

41
Statistique en grande dimension F. Lavancier

Si X 0 X n’est pas inversible, l’estimateur par MCO n’est pas défini de façon
unique et sa variance est infinie. L’EQM de β̂ridge est donc toujours inférieure
dans ce cas. On suppose dans la suite que X 0 X est inversible. En notant β̂
l’estimateur des MCO de β, on a
EQM (β̂) = Var(β̂) = σ 2 (X 0 X)−1 ,
qui peut se récrire
EQM (β̂) = σ 2 (X 0 X + λIp )−1 (X 0 X + 2λIp + λ2 (X 0 X)−1 )(X 0 X + λIp )−1 .
En notant A = (X 0 X + λIp )−1 , qui est symétrique, on en déduit que
EQM (β̂) − EQM (β̂ridge ) = A(2λσ 2 Ip + λ2 σ 2 (X 0 X)−1 − λ2 ββ 0 )A0 .
Cette matrice est définie positive si et seulement si (2λσ 2 Ip + λ2 σ 2 (X 0 X)−1 −
λ2 ββ 0 ) l’est. Ceci peut se vérifier en utilisant le fait que A est inversible
et qu’une matrice M est définie positive ssi il existe N inversible telle que
M = N 0 N (décomposition de Choleski). De plus la somme de deux matrices
définies positives est définie positive. Puisque λ2 σ 2 (X 0 X)−1 est définie pos-
itive, une condition suffisante est donc que la matrice (2λσ 2 Ip − λ2 ββ 0 ) soit
définie positive.
Or la matrice ββ 0 est de rang 1, son image étant engendrée par β, et
son unique valeur propre non nulle vaut kβk2 . En effet tout vecteur v de
l’image s’écrit v = ββ 0 u = (β 0 u)β pour un certain u ∈ Rp et on vérifie que
(ββ 0 )β = βkβk2 . Cela implique que les valeurs propres de (2λσ 2 Ip − λ2 ββ 0 )
valent 2λσ 2 (pour p − 1 d’entre elles) et 2λσ 2 − λ2 kβk2 . Ces valeurs propres
sont strictement positives si et seulement si 0 < λ < 2σ 2 /kβk2 , qui est
donc une condition suffisante pour que EQM (β̂) − EQM (β̂ridge ) soit définie
positive.
L’estimateur ridge peut-être vu d’une façon alternative. Il est solution
du problème de minimisation sous contrainte suivant
p
X
2
β̂ridge = argminβ∈Rp kY − Xβk sous la contrainte βj2 ≤ κ
j=1

pour un certain κ > 0. Autrement dit, l’estimateur ridge est construit comme
l’estimateur par MCO mais sous la contrainte que sa norme (donc ses coef-
ficients) n’explose pas. Le Lagrangien associé à ce problème s’écrit
kY − Xβk2 + λ(kβk2 − κ)

42
Statistique en grande dimension F. Lavancier

où λ est le multiplicateur de Lagrange. En annulant le gradient par rapport

à β on retrouve la solution

β̂ridge = (X 0 X + λIp )−1 X 0 Y.

Le paramètre λ est lié à κ par l’identité kβ̂ridge k2 = κ (en annulant la dérivée

du Lagrangien par rapport à λ), mais cette relation ne conduit pas à une écri-
ture simple de β̂ridge en fonction de κ. L’utilisation en pratique de l’estimateur
β̂ridge nécessite de choisir le paramètre κ ou de façon équivalente λ. Etant
donné la forme explicite en fonction de λ, c’est cette paramétrisation qui
est retenue, et le problème en pratique consiste donc à trouver la meilleure
valeur de λ possible.
La proposition 4.2.1 assure que pour λ suffisamment petit, β̂ridge est
préférable à l’estimateur des MCO. Néanmoins si on choisit λ trop petit,
l’intérêt est nul puisque pour λ = 0, β̂ridge coincide avec l’estimateur des
MCO. Par ailleurs la borne supérieure 2σ 2 /kβk2 est inutilisable en pratique
puisqu’on ne connait ni σ 2 ni β. Le choix de λ se fait donc par validation
croisée (voir le chapitre précédent pour les différentes stratégies possibles)
: on choisit le paramètre λ qui conduit à une estimation de l’erreur test
minimale.

Mise en oeuvre sous R :

Fonction glmnet de la librairie du même nom en choisissant α = 0 dans
les options de la fonction. La validation croisée pour choisir λ peut se faire
avec la fonction cv.glmnet.
Ces fonctions commencent par centrer et réduire les variables X avant
d’estimer β sous contraintes. Cette démarche est naturelle
Pp car si les vari-
2
ables ne sont pas à la même échelle, la contrainte j=1 βj ≤ κ impactera
différemment chaque variable. Les estimations retournées par les fonctions
sont néanmoins dans l’échelle initiale.

4.3 Régression Lasso

4.3.1 Principe
La régression Lasso (pour Least Absolute Shrinkage and Selection Operator)
suit la même idée que la régression ridge, mais au lieu de contraindre la

43
Statistique en grande dimension F. Lavancier

norme `2 de β à ne pas exploser, elle contraint sa norme `1 :

p
X
β̂Lasso = argminβ∈Rp kY − Xβk2 sous la contrainte |βj | ≤ κ.
j=1

Cette modification semble minime mais elle a des conséquences importantes :

alors que l’estimateur Ridge a tendance à réduire la valeur des coefficients (ef-
fet shrinkage) pour satisfaire la contrainte, l’estimateur Lasso en annule car-
rément certains tout en réduisant les autres. Ainsi, contrairement au Ridge
qui garde toutes les variables, le Lasso effectue indirectement une sélection
des variables. Cet effet, dû à la forme de la contrainte Lasso, est illustré dans
la figure 4.3 et est rendu explicite par la résolution du problème exposée dans
la partie suivante.
Le Lagrangien associé au problème d’optimisation s’écrit
p
X
2
kY − Xβk + 2λ( |βj | − κ) (4.5)
j=1

où 2λ est le multiplicateur de Lagrange lié à κ par la contrainte pj=1 |βj | = κ

P
(obtenu en annulant la dérivée par rapport à λ). Le choix 2λ plutôt que λ
n’a aucune incidence et est effectué pour simplifier quelques calculs dans
la résolution du problème (voir partie suivante). Comme pour l’estimateur
Ridge, on a le choix de paramétrer le problème par κ ou par λ et c’est ce
dernier choix qui est préféré pour des raisons pratiques.
Le point négatif du Lasso est qu’il n’y a pas de formule explicite pour
la solution β̂Lasso à λ > 0 donné. En particulier, on ne peut pas dériver le
Lagrangien par rapport à β puisque la présence des valeurs absolues le rend
non dérivable sur Rp . La partie suivante présente néanmoins des algorithmes
très efficaces pour la résolution du problème d’optimisation.
D’un point de vue pratique, comme pour les méthodes précédentes (PCR,
PLS, Ridge), on commence par centrer et réduire les variables pour éviter les
effets d’échelle, puis on choisit le paramètre de régularisation λ par validation
croisée.

Mise en oeuvre sous R :

Fonction glmnet de la librairie du même nom en choisissant α = 1 dans
les options de la fonction (choix par défaut). La validation croisée pour
choisir λ peut se faire avec la fonction cv.glmnet. La résolution du Lasso se

44
Statistique en grande dimension F. Lavancier

Figure 4.3: Figure extraite de l’ouvrage ESL. Parmi l’ensemble des valeurs
possibles (β1 , β2 ) lorsque p = 2, on observe le point β̂ représentant
l’estimateur des MCO qui minimise SCR(β) = kŶ − Xβk2 . Les ellipses
rouges sont les lignes de niveau de SCR(β), dont la valeur augmente au fur
et à mesure que l’on s’éloigne de β̂. La solution Lasso (à gauche) est le point
de contact de ces ellipses avec l’ensemble vérifiant la contrainte souhaitée
|β̂1 | + |β̂2 | ≤ κ. De même à droite le point de contact est la solution ridge
vérifiant β̂12 + β̂22 ≤ κ. La forme de la contrainte Lasso facilite des solutions
sur les "coins" de l’ensemble, pour lesquels une coordonnée est nulle.

45
Statistique en grande dimension F. Lavancier

fait dans glmnet par l’algorithme de descente par coordonnée (voir la partie
suivante).
Fonction lars de la librairie du même nom. La validation croisée se fait
avec cv.lars. La résolution dans lars se fait par l’algorithme LARS-Lasso
(voir la partie suivante).
Les fonctions précédentes normalisent les variables automatiquement par
défaut. La partie 4.3.3 donne plus de précisions concernant ces deux fonctions
et leur comparaison.

4.3.2 Résolution
D’après l’expression (4.5) du Lagrangien, on en déduit que pour λ > 0 donné,
l’estimateur Lasso β̂ est le vecteur qui minimise
p
1 X
L(β) = kY − Xβk2 + λ |βj |. (4.6)
2 j=1

En d’autres termes β̂ = argminβ∈Rp L(β).

Proposition 4.3.1. Pour tout λ > 0, il existe toujours une solution Lasso
β̂ au problème d’optimisation précédent. Cette dernière n’est pas nécessaire-
ment unique mais la prévision X β̂ est unique.

Proof. L’existence d’une solution est garantie par la convexité de L(β). Mais
comme cette fonction n’est pas différentiable, la solution n’est pas nécessaire-
ment unique. Soit deux solutions différentes β̂1 et β̂2 et posons β̄ = (β̂1 +β̂2 )/2
leur moyenne. Par convexité de l’application x 7→ kxk2 , on sait que pour tout
x, y, k(x + y)/2k2 ≤ (kxk2 + kyk2 )/2 et l’inégalité est stricte ssi x 6= y. Si
X β̂1 6= X β̂2 , en utilisant le résultat précédent et l’inégalité triangulaire pour
la norme `1 on obtient

L(β̄) < (L(β̂1 ) + L(β̂2 ))/2.

Or par définition de la solution Lasso, pour tout β ∈ Rp , L(β̂1 ) ≤ L(β) et de

même L(β̂2 ) ≤ L(β). En particulier pour β = β̄, on aboutit à la contradiction
L(β̄) < L(β̄), ce qui montre que X β̂1 = X β̂2 .
La proposition suivante établit deux propriétés importantes de l’estimateur
Lasso qui sont à la base des algorithmes de résolution présentés ensuite.

46
Statistique en grande dimension F. Lavancier

Proposition 4.3.2. On a les deux propriétés suivantes.

1. β̂ minimise (4.6) si et seulement si, pour j = 1, . . . , p,
Xj0 (Y − X β̂) = λ sign(β̂j ) si β̂j 6= 0 (4.7)
et
|Xj0 (Y − X β̂)| ≤ λ si β̂j = 0. (4.8)
Ceci implique que pour les variables actives (associées à un β̂j 6= 0),

|Xj0 (Y − X β̂)| = λ (4.9)

tandis que pour les variables non actives (associées à un β̂j = 0), seule
l’inégalité (4.8) est vérifiée.
2. En supposant que chaque variable Xj est centrée et réduite (standard-
isation classique lors d’une procédure Lasso), la j-ème coordonnée β̂j
de la solution Lasso s’exprime en fonction des autres coordonnées de la
façon suivante

0
 si |Rj | ≤ λ,
β̂j = (Rj − λ)/n si Rj > λ, (4.10)

(Rj + λ)/n si Rj < −λ,


où Rj = Xj0 (Y − k6=j Xk β̂k ).

P

Proof. On rappelle les notions d’optimisation suivantes. Soit f une fonction

convexe de Rd dans R. Un sous-gradient s ∈ Rd de f au point x ∈ Rd vérifie
f (x + h) ≥ f (x) + s0 h, ∀h ∈ Rd .
Si la fonction f est différentiable en x, alors le sous-gradient de f en x est
unique et correspond à son gradient en x. Dans le cas contraire, on note
∂f (x) l’ensemble (éventuellement vide) des sous-gradients possibles de f en
x. Selon la règle de Fermat (ou condition d’optimalité du premier ordre), le
point x∗ est le minimum d’une fonction convexe f ssi 0 ∈ ∂f (x∗ ).
La solution Lasso minimise la fonction convexe L(β). On cherche donc
son sous-gradient. Il vaut
p
!
X
∂L(β) = −X 0 Y + X 0 Xβ + λ ∂ |βj | .
j=1

47
Statistique en grande dimension F. Lavancier

Or le sous-gradient de x 7→ |x|, pour x ∈ R, vaut sign(x) si x 6= 0 et corre-

spond à l’intervalle [−1, 1] si x = 0. On en déduit que la j-ème composante
de ∂L(β) vaut
(
X 0 (Xβ − Y ) + λsign(βj ) si βj 6= 0,
∂L(β)j =  j 0 0

Xj (Xβ − Y ) − λ; Xj (Xβ − Y ) + λ si βj = 0.

Donc β̂ est une solution Lasso ssi 0 ∈ ∂L(β̂) ssi les conditions (4.7) et
(4.8) sont réalisées pour j = 1, . . . , p. La condition nécessaire (mais non
suffisante) (4.9) s’en déduit immédiatement.
Concernant le second point de la proposition, on remarque que
X
Xj0 (Y − X β̂) = Xj0 Y − Xj0 Xj β̂j − Xj0 Xk β̂k = Rj − Xj0 Xj β̂j = Rj − nβ̂j ,
k6=j

où la dernière égalité provient de la normalisation Xj0 Xj = n. Ainsi lorsque

β̂j = 0, |Rj | = |Xj0 (Y − X β̂)| ≤ λ d’après (4.8), ce qui montre la première
égalité dans (4.10). Lorsque β̂j 6= 0, nβ̂j = Rj − λsign(β̂j ) en utilisant
(4.7). Si β̂j > 0, on déduit de la première relation que Rj = nβ̂j + λ > 0
et de même si β̂j < 0, Rj < 0, donc β̂j et Rj ont même signe. Ainsi si
Rj ≥ λ > 0, nβ̂j = Rj − λsign(β̂j ) = Rj − λsign(Rj ) = Rj − λ et si Rj < 0,
nβ̂j = Rj + λ.

Algorithme de descente par coordonnée

Cet algorithme utilise la seconde propriété de la proposition 4.3.2 : si on
connait toutes les coordonnées de β̂ sauf la j-ème β̂j , cette dernière s’obtient
explicitement grâce à (4.10). On initialise donc β̂ à une certaine valeur β̂ (0)
puis on met à jour chaque coordonnée itérativement jusqu’à convergence
(c’est à dire jusqu’à ce que la fonction à minimiser L(β) ne décroît plus à
une précision  près) :
(0)
1. Pour j = 1, . . . , p, on calcule Rj = Xj0 (Y − k6=j Xk β̂k ) et on met à
P

jour β̂j par (4.10).

2. Si |L(β̂) − L(β̂ (0) )| < , la convergence est atteinte et l’algorithme

s’arrête, sinon retour à 1 avec β̂ (0) ← β̂.

48
Statistique en grande dimension F. Lavancier

Puisque la fonction L(β) est convexe, la descente par coordonnée converge

nécessairement vers le minimum global de L(β) (ou l’un d’entre eux s’il n’est
pas unique). Il s’agit de la méthode utilisée dans glmnet pour la résolution
du Lasso.

Algorithme LARS-Lasso
Cet algorithme fournit toutes les solutions β̂(λ) pour tout λ ≥ 0, ce que l’on
appelle le "chemin Lasso", pour un coût algorithmique faible. Il exploite les
propriétés suivantes :

• pour λ = ∞, β̂(∞) = 0,

• pour λ = 0, β̂(0) = β̂ M CO où β̂ M CO est l’estimateur par MCO,

• entre ces deux extrêmes, (4.7) montre que β̂(λ) évolue de façon linéaire
par morceaux en fonction de λ. Précisément l’évolution est linéaire
tant que les variables actives restent les mêmes et un changement de dy-
namique s’opère uniquement lorsqu’une nouvelle variable devient active
(ou une active devient inactive dans le cas de la version LARS-Lasso).

L’algorithme LARS (Least angle regression) construit une solution β̂(λ), pour
tout λ ≥ 0, vérifiant (4.9) et (4.8). Puisque seul (4.9) est assuré et non
(4.7), il ne garantit pas d’obtenir la solution Lasso (mais une solution très
proche). L’algorithme LARS-Lasso, présenté dans un second temps, est une
modification de l’algorithme LARS permettant d’obtenir la solution Lasso.

Algorithme LARS. On construit β̂(λ), pour tout λ ≥ 0, en faisant décroitre

λ de +∞ à 0.
Soit λ1 = maxj |Xj0 Y |.
Pour tout λ ≥ λ1 , β̂(λ) = 0.
Pour λ < λ1 , β̂(λ) est obtenu par l’algorithme suivant. On l’initialise en
posant k = 1 et en définissant la première variable active comme étant celle
réalisant |Xj0 Y | = λ1 .

1. Partant de λ = λk , tant que pour toutes les variables non-actives Xj

la condition |Xj0 (Y − X β̂(λ))| < λ reste vérifiée, on fait décroitre λ en
définissant
λ λ
β̂(λ) = β̂(λk ) + (1 − )β̂kM CO , (4.11)
λk λk

49
Statistique en grande dimension F. Lavancier

où β̂kM CO est l’estimateur des MCO de Y sur les k variables actives de

l’étape courante, complété avec des 0 pour les variables non-actives.

2. On nomme λk+1 la dernière valeur de λ issue de l’étape précédente.

La variable non-active Xj ayant réalisée la condition de sortie |Xj0 (Y −
X β̂(λk+1 ))| = λk+1 devient la (k + 1)-ème variable active.

3. Retour à l’étape 1 avec k ← k + 1.

Proposition 4.3.3. L’algorithme LARS fournit une solution β̂(λ) vérifiant,

pour tout λ ≥ 0, (4.9) et (4.8).

Proof. Lorsque λ ≥ λ1 , β̂(λ) = 0 et dans ce cas la condition (4.8) s’écrit

|Xj0 Y | ≤ λ pour tout j. Cette condition est bien vérifiée par définition de
λ1 . De plus, tous les coefficients β̂j (λ) étant nuls, (4.9) n’a pas besoin d’être
vérifiée. On remarque par ailleurs que si λ > λ1 , |Xj0 Y | < λ.
On montre que les propriétés (4.9) et (4.8) sont vraies lorsque λ ≤ λ1
en les montrant à chaque étape de l’algorithme par récurrence sur k. On
vérifiera également à chaque étape que l’inégalité dans (4.8) est stricte pour
les variables non-actives.
(Dans la suite de la preuve, nous supposerons pour simplifier qu’une seule
variable à la fois peut devenir active. La preuve s’adapte sans problème
majeur au cas où plusieurs variables peuvent être actives en même temps.)
A l’étape k = 1, il n’y a qu’une seule variable active. Supposons sans
perte de généralité qu’il s’agit de Xp . Par définition λ1 = |Xp0 Y | et pour
toutes les variables non-actives (c’est à dire autres que Xp ), |Xj0 Y | < λ1 .
En λ = λ1 , on a donc bien l’inégalité stricte dans (4.8) qui est vérifée pour
les variables non actives. Elle le restera tant que l’on reste dans l’étape
courante, par construction de l’algorithme LARS (dès que l’inégalité devient
une égalité, on passe à l’étape suivante). On a par ailleurs β̂(λ1 ) = 0 et
X β̂1M CO = Xp (Xp0 Xp )−1 Xp0 Y . Ainsi
 
λ λ
Xp (Y − X β̂(λ)) = Xp Y − (1 − )Xp (Xp Xp ) Xp Y = Xp0 Y
0 0 0 −1 0
λ1 λ1

dont la valeur absolue vaut λ, ce qui montre (4.9).

Supposons à présent que (4.9) et (4.8) sont vraies jusqu’à l’étape k − 1
et montrons-les pour l’étape k. Toutes les variables inactives satisfont (4.8)
avec une inégalité stricte tout au long de cette étape par construction de

50
Statistique en grande dimension F. Lavancier

l’algorithme LARS. Concernant les k variables actives à l’étape k, notons X̃k

la sous matrice de X qui les contient. Par définition β̂kM CO est obtenu par
projection de Y sur les k variables actives ce qui implique X β̂kM CO = P[X̃k ] Y .
Par définition de cette projection X̃k0 X β̂kM CO = X̃k0 Y et donc
  
0 0 λ λ M CO
X̃k (Y − X β̂(λ)) = X̃k Y − X β̂(λk ) + (1 − )β̂k
λk λk
λ λ
= X̃k0 Y − X̃k0 X β̂(λk ) − (1 − )X̃k0 X β̂kM CO
λk λk
λ λ
= X̃k0 Y − X̃k0 X β̂(λk ) − (1 − )X̃k0 Y
λk λk
λ
= X̃k0 (Y − X β̂(λk )).
λk
Cela signifie que pour toute les variables actives Xj de l’étape k,

λ 0
Xj0 (Y − X β̂(λ)) = X (Y − X β̂(λk )).
λk j

Par hypothèse de récurrence, |Xj0 (Y − X β̂(λk ))| = λk pour les (k − 1) vari-

ables actives de l’étape précédente. Cette égalité est également vraie pour
la nouvelle variable active de l’étape k puisque il s’agit précisément de la
condition de son activation. Ainsi (4.9) est vérifiée.

L’algorithme LARS ne répond que partiellement au problème LASSO car

il garantit (4.9) et (4.8) mais pas (4.7). La différence a lieu lorsqu’un des
coefficients actifs β̂j change de signe, autrement dit lorsque β̂j (λ) = 0 pour
un certain λ, alors que la variable associée est bien active. Pour prendre
en compte ce cas et répondre au problème Lasso, il suffit d’ajouter le test
suivant à la première étape de l’algorithme LARS
Algorithme LARS-Lasso:

1’. Si pour une variable active Xj de l’étape k, son coefficient β̂j (λ) devient
nul pour un certain λ0 < λk , lorsque β̂(λ) évolue comme décrit en
(4.11), alors on met à jour la dynamique de la façon suivante :

– cette variable Xj devient inactive,

51
Statistique en grande dimension F. Lavancier

– on fait décroitre λ en définissant β̂(λ) comme en (4.11), mais en

s’appuyant uniquement sur les variables actives restantes :
λ λ
β̂(λ) = β̂(λ0 ) + (1 − )β̂ M CO ,
λ0 λ0

où β̂ M CO est l’estimateur des MCO de Y sur les variables actives

restantes de l’étape courante, complété avec des 0 pour les vari-
ables non-actives.

On admet que cette modification permet effectivement de garantir (4.7).

La différence entre LARS et Lasso est illustrée dans la figure 4.4. Les
solutions sont identiques jusqu’à ce qu’un des coefficients actifs revienne à 0,
la modification étant effectuée par l’étape 1’ décrite ci-dessus.

4.3.3 Aparté R : comparaison entre glmnet et lars

Les fonctions glmnet et lars calculent le chemin Lasso complet pour tout
λ. Plus précisément :

• La fonction lars avec l’option par défaut type=’lasso’ calcule le

chemin Lasso en utilisant l’algorithme LARS-Lasso. La solution est
exacte pour tout λ.

• La fonction glmnet avec l’option par défaut α = 1 approche le chemin

Lasso pour une grille de valeurs de λ (choisie automatiquement par
défaut), en lesquelles elle approche la solution Lasso par l’algorithme
de descente par coordonnées. La précision de l’approximation en ces λ
est gérée par l’option thresh qui correspond à la précision  du critère
d’arrêt de l’algorithme. Par défaut thresh vaut 1e − 7. L’approche de
glmnet n’est donc pas aussi exacte qu’avec la fonction lars mais elle
est dans certains cas plus rapide, et surtout glmnet autorise d’autres
pénalités que le Lasso en jouant sur le paramètre α ("elastic net", voir
la partie suivante) et à d’autres modèles que la régression linéaire, dont
en particulier la régression logistique, ce que ne permet pas lars.

Dans une utilisation classique des ces fonctions, on détermine le meilleur λ par
validation croisée (fonction cv.lars pour lars et fonction cv.glmnet pour
glmnet). On accède alors aux coefficients Lasso et aux prévisions associés à

52
Statistique en grande dimension F. Lavancier

Figure 4.4: Figure extraite de l’ouvrage ESL illustrant l’estimation par LARS
(à gauche) et Lasso (à droite) dans un exemple faisant intervenir p = 6
variables explicatives. Chaque figure montre l’évolution des 6 coefficients
β̂j (λ) lorsque λ varie (de gauche à droite) de λ = λ1 , le cas où β̂ = 0, à
λ = 0, le cas où β̂ = β̂ M CO . La légende indique "L1 Arc Length" à la
place de λ : il s’agit d’une quantité valant 0 lorsque λ = λ1 pour croître
linéairement (jusqu’à environ 20 dans cet exemple) lorsque λ décroit jusque
0. La différence entre LARS et Lasso survient lorsqu’un des coefficients actifs
(celui en bleu foncé) revient en 0.

53
Statistique en grande dimension F. Lavancier

ce λ grâce aux fonctions coef et predict en spécifiant l’option s=λ dans ces
fonctions. Une normalisation des données est effectuée par défaut dans lars
et glmnet, mais la transformation inverse est automatiquement effectuée au
moment de retourner les coefficients et les prévisions, qui sont donc à la
bonne échelle des données initiales.
Mais glmnet et lars n’utilisent pas la même définition par défaut du
paramètre de régularisation λ et ne normalisent pas de la même manière
les données. Dans l’utilisation classique qu’on a de ces fonctions, cela n’a
pas d’importance car le choix de λ renvoyé par cv.lars (respectivement
par cv.glmnet) est interprété correctement par lars (respectivement par
glmnet), et les normalisations sont bien sûr gérées de façon cohérente par
chaque fonction.
Cependant il est difficile d’estimer exactement le même modèle (i.e. avec
la même pénalité) avec ces deux fonctions, ce qui pourrait être utile afin de
comparer leurs approches. En fait, la fonction lars minimise
p
1 2
X
kY − Xβk + λ |βj |, (4.12)
2 j=1

où Y et X sont les versions normalisées de lars tandis que la fonction glmnet

minimise p
1 2
X
kY − Xβk + λ |βj |,
2n j=1

où Y et X sont les versions normalisées de glmnet (différentes de lars).

Si on souhaite estimer exactement le même modèle, pour la même pénal-
ité, avec ces deux fonctions, il faut :

• Travailler sans standardisation (option normalize=F pour lars et standardize=F pour

glmnet). Il est toujours possible de normaliser comme on le souhaite
les données au préalable.

• Préciser mode="lambda" pour lars lorsqu’on spécifie la valeur λ souhaitée

dans coef ou predict, pour que la pénalisation ait le même sens
qu’avec glmnet.

• Si on a choisi la pénalité λ pour lars, il faut choisir la pénalité λ/n

pour glmnet (d’après la forme des fonctions à minimiser ci-dessus).

54
Statistique en grande dimension F. Lavancier

Exemple : pour estimer un modèle associé à la pénalité λ = 4 dans (4.12)

fit.lars=lars(x, y, normalize = F)
beta.lars=coef(fit.lars, s=4, mode="lambda")

fit.glmnet=glmnet(x, y, standardize = F)
beta.glmnet=coef(fit.glmnet, s = 4/n)

Pour affiner l’approximation de glmnet (si cette dernière est trop différente de
la solution lars), on peut imposer un critère d’arrêt plus strict dans glmnet
avec l’option thresh, par exemple thresh = 1e-20. Cela nécessitera peut-
être d’augmenter le nombre maximal d’itérations autorisées pour atteindre la
convergence avec l’option maxit. De même l’estimation des coefficients pour
le λ choisi (ici 4/n) correspond à une interpolation linéaire des coefficients
estimés aux λ les plus proches dans la grille choisie dans glmnet. En ajoutant
l’option exact=T à coef, une nouvelle grille de λ est calculée au voisinage
du λ choisi pour affiner l’estimation des coefficients. Dans ce cas, il faut
également rappeler en option les données x=x,y=y. Ce raffinement donnerait
donc

fit.glmnet=glmnet(x, y, standardize = F, thresh = 1e-20)

beta.glmnet=coef(fit.glmnet, s = 4/n, exact = T, x=x, y=y)

4.3.4 Aspects théoriques

L’étude théorique de l’estimateur Lasso est délicate. Par exemple son biais
et sa variance sont non explicites en général.
En supposant que le vrai modèle ne contient qu’un nombre limité de
variables explicatives (hypothèse de parcimonie), deux types de questions
motivent la théorie :

• la qualité d’estimation du Lasso : les coefficients estimés sont-ils de

bonnes approximations des vrais coefficients ?

• le pouvoir de sélection du Lasso : à quel point l’estimateur Lasso est-il

capable de retrouver les bonnes variables ?

Ces deux objectifs ne se réalisent pas nécessairement en même temps. Par

exemple, dans un modèle de régression Gaussien, l’estimateur par MCO est
consistant dès que (X 0 X)−1 → 0 mais il ne sélectionne aucune variable car

55
Statistique en grande dimension F. Lavancier

aucune de ses coordonnées n’est nulle presque surement (sa loi est Gaussi-
enne).
Pour illustrer les principales propriétés du Lasso, on se restreint dans la
proposition suivante au cas particulier où les variables explicatives sont cen-
trées, réduites et non-corrélées (i.e. orthogonales) entre elles. Cela revient
à supposer que X 0 X = nIp . (Cette situation est largement artificielle car
en pratique il y a toujours une corrélation empirique non nulle, même min-
ime, entre les variables explicatives. De plus cette hypothèse implique que
nécessairement p ≤ n car il n’est pas possible d’avoir plus de n vecteurs or-
thogonaux 2 à 2 dans Rn ). Dans ce cas, la solution Lasso à λ > 0 fixé est
explicite (voir TD) et vaut, pour j = 1, . . . , p,

0
 si |β̂jM CO | ≤ λ/n,
β̂(λ) = β̂jM CO − λ/n si β̂jM CO > λ/n,

 M CO
β̂j + λ/n si β̂jM CO < −λ/n,

où β̂jM CO désigne l’estimateur MCO de la régression. Les deux caractéris-

tiques principales de l’approche Lasso apparaissent : elle sélectionne les
variables les plus pertinentes (celles pour lesquelles |β̂jM CO | > λ/n dans le
cas précédent) et elle réduit les coefficients MCO des variables sélectionnées
en rapprochant leur valeur de 0, ce qui est appelé la propriété de "seuillage
doux" du Lasso.
La proposition suivante montre que si le paramètre de régularisation λ
tend√vers l’infini avec n à la bonne vitesse (moins vite que n mais plus vite
que n), l’estimateur Lasso est consistant et sélectionne les bonnes variables
asymptotiquement. Elle souligne également un défaut bien connu du Lasso
: à cause de sa propriété de seuillage doux, les coefficients des variables
"importantes" ont tendance à être sous-estimés (en valeur absolue).
Dans le cas général (au delà de X 0 X = nIp ), le même type de propriétés
a été établi dans les années 2000-2010, mais sous une condition importante
assez restrictive appelée "condition d’irrépresentabilité". Cette dernière sup-
pose que les variables non-pertinentes ne doivent pas être trop corrélées avec
les variables pertinentes (condition évidemment vérifiée si X 0 X = nIp ). Cette
condition témoigne d’un autre défaut du Lasso : cette méthode est sensible
au problème de multicolinéarité, comme les MCO (voir TD à ce sujet).
Proposition 4.3.4. Dans un modèle de régression Gaussien de paramètre
β pour lequel X 0 X = nIp , l’estimateur Lasso β̂(λn ) associé au paramètre de
régularisation λn vérifie les propriétés suivantes.

56
Statistique en grande dimension F. Lavancier

i) Si λn /n → 0, il est consistant : β̂(λn ) → β en probabilité lorsque

n → ∞.
√
ii) Si λn /n → 0 et λn / n → ∞, il sélectionne asymptotiquement les
bonnes variables dans le sens où, en notant J l’ensemble des indices j
pour lesquels βj = 0,
 
P ∀j ∈ J , β̂j (λn ) = 0 et ∀j ∈
/ J , β̂j (λn ) 6= 0 → 1.

iii) L’estimation des coefficients non nuls est biaisée vers 0, dans le sens
où |E(β̂j (λn ))| < βj dès que |βj | > λn /n.

Proof. i) Pour la consistance on écrit, pour tout  > 0, pour tout j = 1, . . . , p,

 
 
M CO λn
P |β̂j (λn ) − βj | >  = P |β̂j (λn ) − βj | > , β̂j >
n
 
M CO λn
+ P |β̂j (λn ) − βj | > , β̂j <−
n
 
M CO λn
+ P |β̂j (λn ) − βj | > , |β̂j |< . (4.13)
n

On nomme T1 + T2 + T3 la décomposition précédente et on va montrer que

chaque terme de cette somme tend vers 0. Pour T1 ,
 
M CO λn
T1 = P |β̂j (λn ) − βj | > , β̂j >
n
 
M CO λn M CO λn
= P |β̂j − − βj | > , β̂j >
n n
 
λn
≤ P |β̂jM CO − − βj | > 
n
   
M CO λn M CO λn
= P β̂j − βj >  + + P β̂j − βj < − +
n n
   
≤ P β̂jM CO − βj >  + P β̂jM CO − βj < −/2

où la dernière inégalité est vraie dès que λn /n < /2 ce qui est garanti à partir
d’un certain n par l’hypothèse λn /n → 0. Puisque (X 0 X)−1 = n−1 Ip → 0,
β̂jM CO est consistant et les deux probabilités dans la majoration ci-dessus

57
Statistique en grande dimension F. Lavancier

tendent vers 0. Cela montre que T1 → 0. Un raisonnement similaire conduit

de même à T2 → 0. Pour le troisième terme :
   
M CO λn M CO λn
T3 = P |β̂j (λn ) − βj | > , |β̂j |< = P |βj | > , |β̂j |< .
n n

Si βj = 0, |βj | >  est impossible et T3 = 0. Si βj 6= 0, en supposant que n

est suffisamment grand pour garantir λn /n < |βj |/2 (ce qui est possible car
λn /n → 0 par hypothèse), on a
 
M CO λn  
T3 ≤ P |β̂j |< < P |β̂jM CO | < |βj |/2 .
n

Or |β̂jM CO | < |βj |/2 ssi −|βj |/2 − βj < β̂jM CO − βj < |βj |/2 − βj , ce qui
implique |β̂jM CO − βj | > |βj |/2 (considérer les cas βj > 0 et βj < 0). Ainsi
 
T3 ≤ P |β̂jM CO − βj | > |βj |/2

qui tend vers 0 par consistance de β̂jM CO .

La consistance de l’estimateur Lasso est ainsi montrée.
ii) Pour la propriété de sélection de variables du Lasso,
 
P ∀j ∈ J , β̂j (λn ) = 0 et ∀j ∈/ J , β̂j (λn ) 6= 0
 
M CO λn M CO λn
= P ∀j ∈ J , |β̂j |< et ∀j ∈
/ J , |β̂j |> .
n n

Or le modèle étant Gaussien avec X 0 X = nIp , β̂ M CO ∼ N (β, σ 2 Ip /n). Cela

montre que les coordonnées β̂jM CO sont indépendantes entre elles. Ainsi
 
P ∀j ∈ J , β̂j (λn ) = 0 et ∀j ∈ / J , β̂j (λn ) 6= 0
   
Y
M CO λn Y M CO λn
= P |β̂j |< P |β̂j |> . (4.14)
j∈J
n n
j ∈J
/

√
Pour tout j ∈ J , n β̂jM CO suit la loi Z ∼ N (0, σ 2√ ). Donc toutes les
probabilités dans le premier produit valent P(|Z| < λn / n), qui tend vers 1

58
Statistique en grande dimension F. Lavancier

√ √
/ J , n (β̂jM CO −
car par hypothèse λn / n → ∞. Par ailleurs, pour tout j ∈
βj ) ∼ Z et dans ce cas on obtient
   
M CO λn M CO λn
P |β̂j |> = 1 − P |β̂j |≤
n n
√ √
 
λn λn
= 1 − P − √ − nβj ≤ Z ≤ √ − nβj
n n
√ |βj |
 
≥ 1 − P |Z| ≥ n
2

dès que λn /n < |βj |/2, ce qui arrive nécessairement à partir d’un certain n
par l’hypothèse λn /n → 0. Cette
 dernière probabilité tend vers 0, ce qui
montre que P |β̂jM CO | > λn /n tend vers 1 pour tout j ∈ / J . Ainsi (4.14)
tend vers 1 lorsque n → ∞ ce qui montre le résultat annoncé.
iii) Enfin pour le calcul du biais, on utilise la décomposition

E(β̂j (λn ))
     
= E β̂j (λn )1|β̂ M CO |≤ λn + E β̂j (λn )1β̂ M CO > λn + E β̂j (λn )1β̂ M CO <− λn
j n j j
   n  n

M CO λ n M CO λ n
=0+E β̂j − 1β̂ M CO > λn + E β̂j + 1β̂ M CO <− λn .
n j n n j n

On va montrer que si βj > λn /n, alors 0 ≤ E(β̂j (λn )) < βj . Le cas contraire
βj < −λn /n se traite de la même manière pour montrer que βj < E(β̂j (λn )) ≤
0. On suppose donc βj > λn /n et on note f la densité de β̂jM CO , c’est à dire
d’une loi N (βj , σ 2 /n), et f ∗ la densité d’une loi N (0, σ 2 /n). On a d’une part
Z ∞ Z − λnn
λn λn
E(β̂j (λn )) = (x − )f (x)dx + (x + )f (x)dx
λn
n
n −∞ n
Z ∞ Z − λn −βj
λn ∗ n λn
= (x + βj − )f (x)dx + (x + βj + )f ∗ (x)dx
λn
n
−βj n −∞ n
Z 0 Z ∞ Z − λn −βj
n
∗ ∗
≥ xf (x)dx + xf (x)dx + xf ∗ (x)dx
λn
n
−βj 0 −∞

en utilisant le fait que βj − λn /n ≥ 0 et βj + λn /n ≥ 0 et en décomposant le

premier intervalle d’intégration en deux. Le changement de variables x → −x

59
Statistique en grande dimension F. Lavancier

pour les intégrales de gauche et de droite conduit finalement à

Z βj + λnn
E(β̂j (λn )) ≥ xf ∗ (x)dx ≥ 0.
βj − λnn

D’autre part, en reprenant la décomposition initiale,

E(β̂j (λn ))
   λ   λn
 
λn

M CO n M CO M CO
= E β̂j 1 − 1|β̂ M CO |< λn − P β̂j > − P β̂j <−
j n n n n
Z λn /n     
λn λn λn
= βj − xf (x)dx − P β̂jM CO > − P β̂jM CO < − .
−λn /n n n n
(4.15)

Puisque βj > λn /n,

   
M CO λn 
M CO
 
M CO

M CO λn
P β̂j <− < P β̂j < βj = P β̂j > βj < P β̂j >
n n

ce qui montre que le dernier terme dans (4.15) est négatif. Concernant le
terme intégrale, en utilisant le fait que f est croissante sur ] − ∞; λn /n]
lorsque βj > λn /n, on en déduit que pour tout x < 0, xf (x) ≥ xf (0) et pour
x ∈ [0; λn /n], xf (x) ≥ xf (0), d’où
Z λn /n Z 0 Z λn /n
xf (x)dx = xf (x)dx + xf (x)dx
−λn /n −λn /n 0
Z 0 Z λn /n
≥ f (0) xdx + f (0) xdx = 0.
−λn /n 0

On connait ainsi le signe de chaque terme dans (4.15), d’où l’on déduit que
E(β̂j (λn )) < βj .

4.4 Quelques généralisations des régressions Ridge

et Lasso
Les qualités et défauts principaux de la régression Ridge sont :

60
Statistique en grande dimension F. Lavancier

• Qualités : robuste à la multicolinéarité, robuste à la grande dimension

• Défauts : ne sélectionne aucune variable mais rétrécit les coefficients

Et pour la régression Lasso :
• Qualités : robuste à la grande dimension, sélectionne les variables

• Défauts : sensible à la multicolinéarité, biais vers 0 des coefficients

importants.

4.4.1 Elastic net

Cette méthode est un compromis entre Ridge et Lasso. L’idée est de tirer
partie des qualités de sélection de l’estimateur Lasso, tout en garantissant une
meilleure robustesse en cas de multicolinéarité, propriété propre à l’estimateur
Ridge. L’estimateur elastic-net est solution du problème :
p
X
2
(1 − α)βj2 + α|βj | ≤ κ,


β̂en = argminβ∈Rp kY −Xβk sous la contrainte
j=1

ou en l’écrivant à l’aide du Lagrangien

p
X
2
(1 − α)βj2 + α|βj | ,


β̂en = argminβ∈Rp kY − Xβk + λ
j=1

où λ est le paramètre de régularisation (lié à κ) qu’il convient de choisir en

pratique.
On remarque que la contrainte est un compromis entre la contrainte Lasso
(le cas α = 1) et la contrainte Ridge (le cas α = 0). Une représentation du
domaine des contraintes en dimension 2 pour différentes valeurs de α est
donnée dans la figure 4.5.
D’un point de vue pratique, pour α donné, on choisit le paramètre λ par
validation croisée puis β̂en peut être obtenu par un algorithme de descente
par coordonnée. Le paramètre α peut également être choisi par validation
croisée.

Mise en oeuvre sous R :

Fonctions glmnet et cv.glmnet de la librairie glmnet, en choisissant en
option la valeur de α.

61
Statistique en grande dimension F. Lavancier

α=1 α = 0.8 α = 0.6

α = 0.4 α = 0.2 α=0

Figure 4.5: Domaine des contraintes elastic-net en dimension 2, i.e. (1 −

α)(β12 + β22 ) + α(|β1 | + |β2 |) ≤ κ, pour différentes valeurs de α. Le cas α = 1
correspond au Lasso et α = 0 au Ridge.

62
Statistique en grande dimension F. Lavancier

4.4.2 Gauss-Lasso
Pour réduire le biais dû au seuillage doux du Lasso, l’idée de la procédure
Gauss-Lasso est d’effectuer l’estimation en deux étapes :
1. On effectue une régression Lasso pour connaître le support, c’est à dire
les variables associées à des coefficients non nuls.
2. On effectue une régression standard par MCO sur les variables sélec-
tionnées à la première étape, sous réserve que leur nombre soit petit
devant n.
On fait ainsi confiance à la propriété de sélection du Lasso, mais on s’appuie
sur les MCO pour l’estimation des coefficients importants. Cette procédure
en deux étapes est assez répandue en pratique.
Attention, si un grand nombre de variables a été sélectionné par Lasso,
cette démarche n’est pas raisonnable car on retombe alors dans les écueils
de l’estimation par MCO en présence d’un trop grand nombre de variables
: la variance d’estimation sera élevée. De plus, bien que la démarche puisse
sembler naturelle, il n’y a (à ma connaissance) aucune garantie théorique
permettant d’appuyer cette procédure et les études par simulation (voir TP)
conduisent à des résultats peu convaincants.
Mise en oeuvre sous R :
Une procédure Lasso avec glmnet ou lars, suivie d’une estimation par
MCO standard avec lm.

4.4.3 Adaptive Lasso

Cet estimateur est construit comme le Lasso, mais en permettant une pé-
nalité différente pour chaque coefficient, cette dernière étant d’autant plus
importante que le coefficient semble proche de 0. Cet a priori sur la valeur
des coefficients est donné par un estimateur préliminaire β̂. L’estimateur
Lasso adaptatif est précisément défini comme solution de
p
2
X |βj |
β̂AL = argminβ∈Rp kY − Xβk + λ .
j=1 |β̂j |

La pénalité propre au coefficient βj est donc λ/|β̂j |. Elle est d’autant plus
forte que |β̂j | est faible, ce qui conduira la procédure Lasso à annuler plus
facilement ce coefficient.

63
Statistique en grande dimension F. Lavancier

Pour l’estimateur préliminaire β̂, les choix usuels sont l’estimateur par
MCO s’il est calculable, l’estimateur Ridge ou l’estimateur Lasso classique.
L’intérêt de cette méthode est sa capacité à mieux sélectionner les vari-
ables pertinentes que les méthodes précédentes.
La solution peut s’obtenir à l’aide d’un algorithme de type LARS ou d’une
descente par coordonnées.

Mise en oeuvre sous R :

Une fois le vecteur des estimations préliminaires β̂ obtenu, il suffit de pré-
ciser penalty.factor=1/abs(β̂) en option des fonctions glmnet et cv.glmnet.

4.4.4 Logistic Lasso

Pour rappel, la régression logistique consiste à modéliser une variable binaire
Y en fonction des variables explicatives X de la façon suivante. En supposant
que Y prend les valeurs 0 et 1, le modèle s’écrit pour l’individu i
0
exi β
P(yi = 1| X•i = x0i ) = 0 .
1 + exi β
En notant π(xi ) la probabilité précédente, on en déduit que les Yi sont in-
dépendants et suivent une loi de Bernoulli B(π(xi )), d’où la vraisemblance
n
Y
V (β; Y ) = π(xi )yi (1 − π(xi ))1−yi .
i=1

En utilisant la forme de π(xi ), on obtient la log-vraisemblance

n
0
X
L(β; Y ) = yi x0i β − ln(1 + exi β ). (4.16)
i=1

L’estimateur du maximum de vraisemblance β̂ maximise la log-vraisemblance,

ou de façon équivalente minimise son opposé

β̂ = argminβ∈Rp − L(β; Y ).

La solution n’est pas explicite et s’obtient par un algorithme d’optimisation

(dans la fonction glm sous R, il s’agit par défaut de l’algorithme IRLS, pour
"Iteratively Reweighted Least Squares").

64
Statistique en grande dimension F. Lavancier

La régression logistique Lasso suit le même principe que la régression

Lasso : dans le problème d’optimisation définissant β̂, on ajoute une con-
trainte sur la norme `1 des coefficients.
p
X
β̂Logit−Lasso = argminβ∈Rp − L(β; Y ) sous la contrainte |βj | ≤ κ,
j=1

où L est donnée par (4.16). De façon équivalente

p
X
β̂Logit−Lasso = argminβ∈Rp − L(β; Y ) + λ |βj |,
j=1

où λ est le paramètre de régularisation (lié à κ) à choisir par validation

croisée.
La résolution du problème d’optimisation précédent peut se faire par un
algorithme de descente par coordonnées.

Mise en oeuvre sous R :

Fonction glmnet avec l’option family="binomial".

4.4.5 Group-Lasso
Dans certains modèles, il peut être préférable de sélectionner les variables par
groupe. Par exemple, une variable qualitative à k modalités est encodée dans
un modèle de régression linéaire (ou logistique) par k−1 variables indicatrices.
Si l’on souhaite supprimer la variable qualitative du modèle, cela revient à
supprimer le groupe des k − 1 variables indicatrices correspondantes.
En supposant qu’on a formé G groupes de variables, chacun contenant
kg variables, la régression Group-Lasso pénalise les groupes de la manière
suivante
G
X
2
β̂Group−Lasso = argminβ∈Rp kY − Xβk + λ kg kβg k
g=1

où βg désigne le sous-vecteur de β associé aux kg coefficients du groupe g

et kβg k est sa norme euclidienne. La forme de la pénalité conduit ainsi à
une solution qui annule certaines normes kβg k et donc tous les coefficients
du groupe g.

65
Statistique en grande dimension F. Lavancier

p sont de taille kg = 1, on
Dans le cas particulier où tous les groupes
retrouve le Lasso classique car alors kβg k = βg2 = |βg |.
La même idée se généralise à la régression logistique Lasso.

Mise en oeuvre sous R :

Fonction gglasso de la librairie du même nom. Elle permet une régres-
sion linéaire Group-Lasso par défaut, mais également une régression logis-
tique Group-Lasso avec l’option loss = "logit". Les groupes de variables
sont spécifiés par l’option group sous la forme d’un vecteur d’indices. Par ex-
emple group=c(1,1,2,3,2) signifie que les deux premières variables forment
un groupe, que la seconde et la dernière variable forment un autre groupe,
et que la quatrième variable est seule.

66
Chapter 5

Tests multiples

Lors d’une analyse statistique, il n’est pas rare d’effectuer de nombreux tests
statistiques.
Exemple 1: En génomique. On souhaite tester si l’expression d’un gêne
est différent entre deux conditions expérimentales (en faisant typiquement
un test d’égalité de moyennes). Cette procédure est effectuée pour tous les
gênes mesurés, ce qui peut représenter des dizaines de milliers de tests. His-
toriquement la génomique est le domaine dans lequel s’est le plus développé
la théorie des tests multiples.
Exemple 2: En régression linéaire lorsque p est grand. On est amené dans
ce contexte à tester la significativité de chaque variable, ce qui représente p
tests.
Exemple 3: Corrélation fortuite. En présence de p variables, on peut
tester la significativité des corrélations linéaires entre chaque variable. Cela
représente p(p − 1)/2 corrélations à tester. Certaines risquent à tort d’être
considérées comme étant significatives, ce sont des corrélations fortuites.
Lorsqu’on effectue de nombreux tests statistiques, chacun étant associé
à un risque de première espèce α, on peut s’attendre à détecter de nom-
breux faux positifs. Ce chapitre expose en détail ce problème et les solutions
courantes pour contrôler ce risque.

67
Statistique en grande dimension F. Lavancier

5.1 Présentation du problème et notations

On suppose qu’on a m tests à effectuer, chacun étant associé à une hypothèse
nulle H0,i et une hypothèse alternative H1,i pour i = 1, . . . , m. On note I0
l’ensemble des indices i pour lesquels H0,i est vrai et on note m0 son cardinal.
Evidemment I0 et m0 ne sont pas connus en pratique.

Dans la plupart des cas, les hypothèses sont de la forme

H0,i : µi = 0 contre H1,i : µi > 0,

ou H0,i : µi = 0 contre H1,i : µi 6= 0,

pour certains paramètres inconnus µ1 , . . . , µm . Par exemple µi représente

la différence d’expression d’un gêne entre deux conditions expérimentales,
ou µi représente un paramètre dans une régression linéaire, ou µi représente
une corrélation linéaire entre deux variables. Cela motive les terminologies
suivantes :

• Les positifs P sont les indices i pour lesquels on a rejeté l’hypothèse

nulle H0,i

• Les faux-positifs F P sont les indices i positifs à tort (autrement dit

H0,i est rejeté alors que i ∈ I0 ).

• Les vrais-positifs V P sont les indices i positifs à raison (autrement dit

H0,i est rejeté et i ∈
/ I0 ).

On a évidemment P = F P + V P .

Lors de la mise en oeuvre des m tests, on obtient m p-values p̂i . On fera

l’hypothèse non restrictive suivante :
Hypothèse : Pour tout i ∈ I0 , p̂i ∼ U([0, 1]).
Exemples : si la statistique de test admet une loi continue de fonction de
répartition F et la région critique est de la forme

• RCα = {T > F −1 (1 − α)}, alors p̂ = 1 − F (T ) ∼ U([0, 1]).

• RCα = {T < F −1 (α/2)}∪{T > F −1 (1−α/2)}, alors p̂ = 2 min(F (T ), 1−

F (T )) ∼ U([0, 1]).

68
Statistique en grande dimension F. Lavancier

Preuve des exemples. Pour le premier cas, en notant F −1 l’inverse de F :

P(p̂ ≤ p) = P(1 − F (T ) ≤ p) = P(T ≥ F −1 (1 − p)) = 1 − F (F −1 (1 − p)) = p.

Pour le second cas,

P(p̂ ≤ p) = P(p̂ ≤ p, F (T ) ≤ 1 − F (T )) + P(p̂ ≤ p, F (T ) > 1 − F (T ))

   
−1 p −1 1 −1 p −1 1
= P T ≤ F ( ), T ≤ F ( ) + P T ≥ F (1 − ), T > F ( ) .
2 2 2 2

Comme p ≤ 1, on sait que p/2 ≤ 1/2 et 1 − p/2 ≥ 1/2. Puisque F −1 est

croissante, cela implique F −1 ( p2 ) ≤ F −1 ( 21 ) et F −1 (1 − p2 ) ≥ F −1 ( 12 ). D’où
 p   p  p p
P(p̂ ≤ p) = P T ≤ F −1 ( ) + P T ≥ F −1 (1 − ) = + (1 − (1 − )) = p.
2 2 2 2

On souhaite se donner une règle de décision sur les m p-values p̂i . On

peut la construire de deux manières
• A partir de quel seuil τ (pour les p̂i ) rejette-t-on les hypothèses nulles?
Cela aboutit à la règle : si p̂i < τ , on rejette H0,i .

• De façon alternative, si on ordonne les p̂i de p̂(1) à p̂(m) : à partir de quel

rang k̂ considère-t-on que les p-values p̂(1) à p̂(k̂) conduisent à rejeter
les hypothèses nulles correspondantes, tandis que pour p̂(k̂+1) à p̂(m) on
ne les rejette pas?
Ces deux points de vue diffèrent dans la mesure où τ n’est pas aléatoire alors
que k̂ oui. Le lien entre ces deux points de vue est illustré dans la figure 5.1.

L’objectif d’une bonne procédure de tests multiples est de trouver une

règle de décision comme ci-dessus qui maximise le nombre de V P tout en
minimisant les F P .

Si on se fixe un seuil τ indépendant de m, par exemple τ = 5%, on aura

en moyenne
m
!
X X
E(F P ) = E 1{i∈I0 , p̂i <τ } = P(p̂i < τ ) = m0 τ. (5.1)
i=1 i∈I0

69
Statistique en grande dimension F. Lavancier

●
0.30

●
0.15

●
● ●

τ
● ●
● ● ● ●
●
● ●
0.00

● ●
● ● ●

^
5 k 10 15 20

Figure 5.1: Représentation de 20 p-values ordonnées de la plus petite p̂(1) à

la plus grande p̂(20) . Si on se fixe un seuil τ > 0, on rejette tous les tests
ayant une p-value inférieure à τ . De façon équivalente, on rejette les tests
associés aux k̂ plus petites p-values. Ici k̂ = 8.

70
Statistique en grande dimension F. Lavancier

Par exemple, si on fait m = 1000 tests qui sont en théorie tous négatifs
(m = m0 ) en utilisant le seuil τ = 5%, on obtient en moyenne 50 F P . Pour
contrôler l’apparition des F P , il faut donc affiner la règle de décision en
tenant compte de m.

5.2 Principe des tests multiples : contrôler le

F W ER ou le F DR
Rappel : lorsqu’on effectue un seul test (H0 contre H1 ), on adopte générale-
ment le principe de Neyman-Pearson. On se donne la probabilité de décider
H1 à tort (faux positif), c’est à dire le risque de première espèce α, et on
privilégie, si on a le choix, le test qui, à α donné, maximise la puissance, c’est
à dire la probabilité de décider H1 à raison (vrai positif).

On trouve deux façons principales d’étendre ce principe au cas des tests

multiples :
1. par le F W ER (Family-Wise Error Rate, taux d’erreur par famille) :
on décide de fixer
F W ER = P(F P > 0)
à α, ou autrement dit P(F P = 0) à 1−α, et on privilégie les procédures
qui, à α fixé, maximisent le nombre de V P .
2. par le F DR (False Discovery Rate, taux de faux positifs) : on note
F DP la proportion de F P parmi P , avec la convention F DP = 0 si
P = 0, c’est à dire F DP = (F P/P ) 1P ≥1 . On décide alors de fixer
F DR = E(F DP ) = E((F P/P ) 1P ≥1 )
à α et on privilégie les procédures qui, à α fixé, maximisent le nombre
de V P .
La proposition suivante montre que lorsque m = 1, les deux critères coïnci-
dent et correspondent au risque de première espèce usuel, tandis que dans le
cas général F DR ≤ F W ER. Donc F DR est un critère moins exigeant que
F W ER dans le sens où une procédure assurant F W ER ≤ α assure automa-
tiquement F DR ≤ α, mais le contraire est faux. En particulier on fournit
généralement plus de V P avec une procédure assurant F DR ≤ α qu’avec la
même procédure assurant F W ER ≤ α.

71
Statistique en grande dimension F. Lavancier

Proposition 5.2.1. On a les propriétés suivantes :

(i) En présence d’un seul test (m = 1), F W ER = F DR correspond au

risque usuel de première espèce.

(ii) Dans le cas général (m ≥ 1), on a F DR ≤ F W ER.

Proof. (i) Si m = 1, F W ER = P(F P = 1) correspond exactement à la

probabilité de rejeter H0 à tort, c’est à dire au risque de première espèce
usuel. Par ailleurs, toujours si m = 1, F DP vaut 1 lorsqu’il y a un faux
positif (situation où F P = P = 1) et 0 sinon, donc son espérance F DR vaut
également P(F P = 1).
(ii) Dans le cas général F P 1P ≥1 ≤ P 1F P ≥1 (ce qui se vérifie en con-
sidérant les cas possibles P = 0 ou P ≥ 1, et F P = 0 ou F P ≥ 1), d’où
(F P/P )1P ≥1 ≤ 1F P ≥1 ce qui implique le résultat en passant à l’espérance.

5.2.1 Procédure de Bonferroni

La procédure de Bonferroni est très simple : pour assurer F W ER ≤ α, il
suffit de fixer le seuil τ = α/m.
La force de cette procédure est qu’elle garantit un contrôle du F W ER
quelle que soit la dépendance entre les p-values des m tests. Son gros défaut
est qu’elle est beaucoup trop conservative si m est grand, dans le sens où le
seuil τ devient tellement petit qu’il y a très peu de positifs détectés, donc
forcément très peu de F P (ce qui est un bon point) mais aussi très peu de
V P.

Proposition 5.2.2. La procédure de Bonferroni assure F W ER ≤ α.

Proof. On observe un F P si pour un test i, p̂i ≤ α/m alors que H0,i est vrai,
ce qui signifie i ∈ I0 . Donc
X
F W ER = P(F P > 0) = P(∃i ∈ I0 , p̂i ≤ α/m) ≤ P(p̂i ≤ α/m)
i∈I0

car la probabilité d’une union est inférieure à la somme des probabilités.

Puisque pour i ∈ I0 , p̂i ∼ U([0, 1]), ces dernières probabilités valent α/m et
on aboutit à F W ER ≤ m0 α/m ≤ α.

72
Statistique en grande dimension F. Lavancier

5.2.2 Procédure de Benjamini-Hochberg

La procédure de BH contrôle le FDR. Il s’agit de la procédure de tests mul-
tiples la plus utilisée en pratique.
On rappelle que le choix d’un seuil τ pour les p-values conduit à P = k̂
positifs, qui sont associés aux k̂ plus petites p-values avec p̂(k̂) ≤ τ < p̂(k̂+1) ,
voir la figure 5.1.
L’heuristique de la procédure de BH est la suivante : on sait que pour un
seuil τ donné E(F P ) = m0 τ ≤ mτ , voir (5.1). On peut donc s’attendre à
ce qu’en moyenne F P ≤ mτ . Tant que τ reste dans l’intervalle [p̂(k̂) , p̂(k̂+1) [,
les conclusions des tests restent identiques et conduisent à k̂ positifs. Donc
pour τ = p̂(k̂) , on peut s’attendre à ce qu’en moyenne F P ≤ mp̂(k̂) , soit
F P/P ≤ mp̂(k̂) /k̂, autrement dit F DR ≤ mp̂(k̂) /k̂. Evidemment cette iné-
galité n’est pas juste car l’aléa et l’espérance n’ont pas été manipulés conven-
ablement. Néanmoins, si on se base sur cette inégalité, on aura F DR ≤ α
dès que mp̂(k̂) /k̂ ≤ α. Il y a plusieurs choix possibles de k̂ conduisant à
cette condition. Parmi ces possibilités, on souhaite maximiser le nombre de
positifs k̂, donc on choisit pour k̂ le plus grand k vérifiant p̂(k) ≤ αk/m.
Suivant cette heuristique, la région de rejet de BH correspond donc à
toutes les p-values p̂i inférieures à p̂(k̂) où k̂ vérifie

k̂ = max{k, p̂(k) ≤ αk/m}. (5.2)

Comme l’illustre la figure 5.2, l’interprétation graphique est simple : p̂(k̂) est
la dernière p-value sous la droite de pente α/m, dans le nuage de points des
p-values ordonnées. En comparaison, la procédure de Bonferroni (pour le
même seuil α) rejette nécessairement moins d’hypothèses.
La proposition suivante montre que le F DR associé à la procédure BH
est bien inférieure à α lorsque les p-values sont indépendantes entre elles. Si
les p-values ne sont pas indépendantes entre elles, le résultat n’est pas vrai en
général mais il le reste sous l’hypothèse que l’inégalité (5.4) présentée dans
le lemme ci-dessous reste vraie. Cela est en particulier vérifié en présence de
certaines dépendances positives entre les p-values.

Proposition 5.2.3. Si les p-values sont indépendantes entre elles, la procé-

dure de Benjamini-Hochberg assure F DR ≤ α.

Proof. Par définition il y a k̂ positifs et l’hypothèse H0,i est rejetée si p̂i est

73
Statistique en grande dimension F. Lavancier

●
0.30

●
0.15

●
● ●

● ● BH
● ● ● ●
●
● ●
Bonferroni
0.00

● ●
● ● ●

5 10 15 20

Figure 5.2: Représentation de 20 p-values ordonnées de la plus petite p̂(1)

à la plus grande p̂(20) . La région de rejet de Bonferroni correspond aux p-
values inférieures à α/m (ligne rouge), tandis que la région de rejet de BH
correspond aux k̂ plus petites p-values, où p̂(k̂) est la dernière p-value sous la
droite de pente α/m (ligne bleu). Ici pour m = 20 et α = 0.05, on a k̂ = 5
tandis que Bonferroni ne rejette qu’un seul test.

74
Statistique en grande dimension F. Lavancier

inférieure à p̂(k̂) . Il s’agit d’un faux positif si de plus i ∈ I0 . Ainsi

P !
i∈I0 1p̂i ≤p̂(k̂) X  1p̂i ≤αk̂/m
  
FP
F DR = E 1k̂≥1 = E 1k̂≥1 ≤ E 1k̂≥1
k̂ k̂ i∈I0 k̂

où l’on a utilisé pour la dernière inégalité le fait que p̂(k̂) ≤ αk̂/m pour la
procédure de BH. En décomposant sur toutes les valeurs possibles de k̂, on
obtient
m
!
X 1p̂i ≤αk̂/m X
F DR ≤ E 1k̂≥1 1k̂=k
i∈I0 k̂ k=0
m  
XX 1p̂i ≤αk/m
= E 1k̂=k
i∈I k=1
k
0
m
XX 1  
= P p̂i ≤ αk/m, k̂ = k
i∈I k=1
k
0
m
XX 1  
= P k̂ = k|p̂i ≤ αk/m P (p̂i ≤ αk/m)
i∈I k=1
k
0
m
α XX  
= P k̂ = k|p̂i ≤ αk/m (5.3)
m i∈I k=1
0

en utilisant le fait que p̂i ∼ U([0, 1]) lorsque i ∈ I0 . Le lemme suivant,

démontré plus bas, est la clé de la preuve.
Lemme 5.2.4. Si les p-values p̂1 , . . . , p̂m sont indépendantes alors pour tout
i = 1, . . . , m,
   
P k̂ ≤ k|p̂i ≤ αk/m ≤ P k̂ ≤ k|p̂i ≤ α(k + 1)/m . (5.4)
 
En notant ak = P k̂ ≤ k|p̂i ≤ α(k + 1)/m , on a pour tout k ≥ 1,
     
P k̂ = k|p̂i ≤ αk/m = P k̂ ≤ k|p̂i ≤ αk/m − P k̂ = k − 1|p̂i ≤ αk/m
 
= P k̂ ≤ k|p̂i ≤ αk/m − ak−1
≤ ak − ak−1

75
Statistique en grande dimension F. Lavancier

en utilisant (5.4). En injectant cette inégalité dans (5.3), on obtient

m
α XX αX
F DR ≤ (ak − ak−1 ) = (am − a0 ).
m i∈I k=1 m i∈I
0 0

On remarque que a0 = 0 car si p̂i ≤ α/m, cela implique que p̂(1) ≤ α/m
 et
donc d’après (5.2) que k̂ ≥ 1, autrement dit a0 = P k̂ ≤ 0|p̂i ≤ α/m = 0.
Ainsi, en utilisant de plus le fait que am ≤ 1 et que m0 ≤ m, on conclut que
αX α
F DR ≤ am ≤ m0 ≤ α.
m i∈I m
0

Preuve du Lemme 5.2.4. Sans perte de généralité on suppose que i = 1 pour

la preuve. On rappelle que par définition k̂ dépend de toutes les p-values et
vaut k̂ = k(p̂1 , . . . , p̂m ) = max{j, p̂(j) ≤ αj/m}. On a
 
  P k(p̂1 , . . . , p̂m ) ≤ k, p̂1 ≤ αk/m
P k(p̂1 , . . . , p̂m ) ≤ k | p̂1 ≤ αk/m = .
P(p̂1 ≤ αk/m)
En notant fi la densité de p̂i , le numérateur vaut par indépendance des p̂i
Z
1{k(x1 ,...,xm )≤k,x1 ≤αk/m} f1 (x1 ) . . . fm (xm )dx1 . . . dxm
Z Z 
= 1{k(x1 ,...,xm )≤k,x1 ≤αk/m} f1 (x1 )dx1 f2 (x2 ) . . . fm (xm )dx2 . . . dxm
Z  
= P k(p̂1 , x2 , . . . , xm ) ≤ k, p̂1 ≤ αk/m f2 (x2 ) . . . fm (xm )dx2 . . . dxm

Ainsi en divisant par P(p̂1 ≤ αk/m) on obtient

 
P k(p̂1 , . . . , p̂m ) ≤ k | p̂1 ≤ αk/m
Z  
= P k(p̂1 , x2 , . . . , xm ) ≤ k | p̂1 ≤ αk/m f2 (x2 ) . . . fm (xm )dx2 . . . dxm .

Pour montrer l’inégalité du lemme, il suffit donc de prouver que pour tout
x2 , . . . , x m ,
 
P k(p̂1 , x2 , . . . , xm ) ≤ k | p̂1 ≤ αk/m
 
≤ P k(p̂1 , x2 , . . . , xm ) ≤ k | p̂1 ≤ α(k + 1)/m .

76
Statistique en grande dimension F. Lavancier

En notant g la fonction p 7→ k(p, x2 , . . . , xm ) définie sur [0, 1], on doit

donc prouver que
   
P g(p̂1 ) ≤ k | p̂1 ≤ αk/m ≤ P g(p̂1 ) ≤ k | p̂1 ≤ α(k + 1)/m .

D’après la définition de k̂, la fonction g ne peut valoir que deux valeurs,

soit g(0) soit g(1), avec un éventuel saut en une valeur p0 de [0, 1] (p0 dépen-
dant de x2 , . . . , xm ) avec la convention que p0 = 0 si g(0) = g(1). En effet,
augmenter une p-value ne peut que diminuer l’indice k̂ retenu pour le rejet des
hypothèses. Ce changement s’opère au plus un fois auquel cas g(1) = g(0)−1,
sinon g(0) = g(1). La fonction est par ailleurs continue à gauche.
On note g −1 (k) = sup{p, g(p) ≥ k} avec la convention que g −1 (k) = −∞
si k > g(0). Il s’agit de l’inverse à droite de g dans le sens où g(g −1 (k)) = k
pour tout k ∈ Im(g). Il est facile de vérifier que {p, g(p) < k} = {p >
g −1 (k)}. En effet

• Si k > g(0) alors g −1 (k) = −∞ et {g(p) < k} = {p ≥ 0} = {p >

g −1 (k)}.

• Si k = g(0) alors g −1 (k) = p0 et {g(p) < k} = {p > p0 } = {p >

g −1 (k)}.

• Si k ≤ g(1), alors g −1 (k) = 1 et {g(p) < k} = ∅ = {p > g −1 (k)}.

Ainsi, quel que soit k,

   
P g(p̂1 ) ≤ k | p̂1 ≤ αk/m = P g(p̂1 ) < k + 1 | p̂1 ≤ αk/m
 
= P p̂1 > g −1 (k + 1) | p̂1 ≤ αk/m
 
−1
= 1 − P p̂1 ≤ g (k + 1) | p̂1 ≤ αk/m
−1 (k+1))
(
1 ≤g
1 − P(p̂P(p̂ 1 ≤αk/m)
si g −1 (k + 1) ≤ αk/m
=
0 sinon.

77
Statistique en grande dimension F. Lavancier

Puisque P(p̂1 ≤ αk/m) ≤ P(p̂1 ≤ α(k + 1)/m)

P(p̂1 ≤g −1 (k+1))
(
  1 − P(p̂ 1 ≤α(k+1)/m)
si g −1 (k + 1) ≤ αk/m
P g(p̂1 ) ≤ k | p̂1 ≤ αk/m ≤
0 sinon,
P(p̂1 ≤g −1 (k+1))

−1
1 − P(p̂1 ≤α(k+1)/m) si g (k + 1) ≤ αk/m

= 0 si αk/m < g −1 (k + 1) ≤ α(k + 1)/m

0 sinon,

P(p̂1 ≤g −1 (k+1))
(
1 − P(p̂ 1 ≤α(k+1)/m)
si g −1 (k + 1) ≤ α(k + 1)/m
≤
0 sinon,
 
= P g(p̂1 ) ≤ k | p̂1 ≤ α(k + 1)/m .

5.2.3 p-values ajustées

Lors de la réalisation d’un test statistique (m = 1), on calcule généralement
la p-value p̂ car cette démarche nous évite de choisir le niveau du test α au
préalable, contrairement à la construction d’une région critique. Le choix de
α peut se faire a posteriori grâce à la règle : si p̂ < α, on rejette l’hypothèse
nulle au niveau α. On peut ainsi observer directement pour quels choix de α
la décision du test se trouve modifiée.
Cet avantage est perdu avec les procédures de Bonferroni et de BH présen-
tées ci-dessus, car même si elles s’appuient sur les m p-values initiales p̂1 , . . . , p̂m ,
elles nécessitent de se fixer un seuil α au préalable (correspondant respective-
ment au contrôle du F W ER ou du F DR). On peut néanmoins modifier les
p-values initiales en p̃1 , . . . , p̃m en s’appuyant sur la correction de Bonfer-
roni (respectivement de BH) pour aboutir à la simple règle de décision : si
p̃i ≤ α, alors on rejette le test i pour la procédure de Bonferroni associée à
F W ER ≤ α (respectivement on rejette le test i pour la procédure de BH
associée à F DR ≤ α).

Proposition 5.2.5.

• Les p-values ajustées au sens de Bonferroni valent, pour i = 1, . . . , m,

p̃i = mp̂i .

78
Statistique en grande dimension F. Lavancier

Le test i est rejeté par la procédure de Bonferroni associée à F W ER ≤

α ssi p̃i ≤ α.

• Les p-values ajustées au sens de BH valent, pour i = 1, . . . , m,

 
m
p̃(i) = min p̂(j) .
i≤j≤m j

La présence des parenthèses en indice indique qu’il s’agit de l’ajustement

de la i-ème plus petite p-value p̂(i) . Pour obtenir p̃i , il suffit de remettre
les p-values dans l’ordre initial.
Le test i est rejeté par la procédure de BH associée à F DR ≤ α ssi
p̃i ≤ α.

Proof. La procédure de Bonferroni associée à F W ER ≤ α rejette le test i

ssi p̂i ≤ α/m ⇔ mp̂i ≤ α ⇔ p̃i ≤ α.
Pour BH, on montre l’équivalence (sans perte de généralité) sur les p-
values ordonnées p̂(1) , . . . , p̂(m) . La procédure de BH associée à F DR ≤ α
rejette le test associé à p̂(i) ssi i ≤ k̂ où
   
j m
k̂ = argmax p̂(j) ≤ α = argmax p̂(j) ≤ α .
1≤j≤m m 1≤j≤m j

Par définition m p̂ ≤ α donc p̃(k̂) ≤ m

k̂ (k̂)
p̂ ≤ α. Ainsi, si la procédure de
k̂ (k̂)
BH rejette le test associé à p̂(i) , alors p̃(i) ≤ p̃(k̂) ≤ α.
Réciproquement, si p̃(i) ≤ α, alors il existe j ≥ i tel que mj p̂(j) ≤ α, ce qui
montre que k̂ ≥ i et donc que le test associé à p̂(i) est rejeté par la procédure
de BH.

Mise en oeuvre sous R : Etant donné le vecteur des p-values initiales

p̂1 , . . . , p̂m , la fonction p.adjust retourne les p-values ajustées p̃1 , . . . , p̃m au
sens de BH avec l’option method="BH".
D’autres corrections de p-values sont proposées dans cette fonction. Soit
elles correspondent à des méthodes contrôlant la FEWR comme la méthode
de Bonferroni (method="bonferroni", method="holm", method="hochberg",
method="hommel"), soit elles correspondent à des méthodes contrôlant la
FDR comme la méthode de BH (method="BH", method="BY").

79