Statistique en Grande Dimension

F. Proïa

Master 2 Ingénierie Statistique

Notes de cours

M2 Ingénierie Statistique Statistique en Grande Dimension 1 / 97

Références du cours

Introduction to High-Dimensional Statistics (C. Giraud).

Chapman & Hall/CRC Monographs on Statistics & Applied Probability (2014).
Statistical Learning with Sparsity : The Lasso and Generalizations
(T. Hastie, R. Tibshirani and M. Wainwright).
Chapman & Hall/CRC Monographs on Statistics & Applied Probability (2015).

M2 Ingénierie Statistique Statistique en Grande Dimension 2 / 97

 Concepts généraux

Table des matières

1 Concepts généraux
Grande dimension
Retour sur la régression linéaire
Fléau de la dimension
Critères AIC et BIC
Validation croisée

2 Régressions pénalisées

3 Régressions PCR et PLS

4 Tests multiples

M2 Ingénierie Statistique Statistique en Grande Dimension 3 / 97

 Concepts généraux

Concepts généraux

On commence ce cours par quelques concepts généraux, indispensables en vue de

comprendre le contenu du module.

Qu'entend-on par grande dimension ? Quels outils statistiques prennent le

relai des procédures usuelles en petite dimension ?

Comment intercepter seulement une petite quantité de variables inuentes

parmi une multitude ?

Quelles dicultés inconnues de la petite dimension surviennent en grande

dimension ?

On se focalisera souvent sur le modèle de régression linéaire gaussienne

Y = X β + ε, ε ∼ Nn (0, σ 2 In ).

Cas d'école facile à appréhender pour introduire de nouvelles notions.

L'estimation de β s'attache à séparer le signal Xβ du bruit ε.

On appelle D le dataset : ici D = (Y , X ).

M2 Ingénierie Statistique Statistique en Grande Dimension 4 / 97

 Concepts généraux Grande dimension

Données en grande dimension

L'émergence de la grande dimension (le big data) est due pour l'essentiel à une
véritable révolution technologique.

Progression très rapide des techniques modernes d'acquisition des données.

Comment analyser de telles données ? Les méthodes statistiques usuelles ne

s'appliquent pas toujours.

Quelques exemples :

En médecine : données omiques, séquençage, imagerie médicale.

Dans le numérique : e-commerce, GPS, moteurs de recherche,

vidéosurveillance, chatbot, intelligence articielle.

Dans la nance : technologie blockchain, détection de fraudes.

Dans le monde industriel : consommation électrique (Linky), robotique.

M2 Ingénierie Statistique Statistique en Grande Dimension 5 / 97

 Concepts généraux Grande dimension

Données en grande dimension

Les puces à ADN (microarrays ) peuvent mesurer l'expression de dizaines de

milliers de gènes simultanément.

Sur la photo ci-dessus : plus de 41000 transcriptions de gènes humains (quantité

d'ARNm associée à un gène dans une cellule) !

1
Et encore, même ces puces sont devenues has been : RNA-Seq (séquençage
haut-débit) et Single-Cell ont déjà pris le relai à peine 20 ans plus tard.

1. © P. Lacroix
M2 Ingénierie Statistique Statistique en Grande Dimension 6 / 97
 Concepts généraux Grande dimension

Données en grande dimension

Le jeu de données Hopx contient l'expression de 770 marqueurs génétiques

répartis sur 20 chromosomes. Il contient aussi une série de mesures faites sur 4
tissus spéciques, le tout sur 29 rats de laboratoire.

Heatmap ci-dessus : seulement 3 chromosomes parmi les 20 à disposition. Une

multitude de données pour expliquer un phénomène, à l'aide de 29 mesures.

M2 Ingénierie Statistique Statistique en Grande Dimension 7 / 97

 Concepts généraux Grande dimension

Données en grande dimension

On cherche à analyser D = (Y , X ) de manière statistique : les covariables X

expliquent-elles la réponse Y?
De manière volontairement péjorative, on pourrait diérencier deux classes de
statistiques.

ème
La statistique du 20 siècle :

Il y a plus d'observations que de covariables : n>p et p petit.

On établit des résultats à n et p xés ou, à défaut, lorsque n → +∞
(consistance, normalité asymptotique, ...)
Régressions usuelles, tests de signicativité des coecients β, etc.

ème
La statistique du 21 siècle :

Des milliers de données pour chaque individu : p très grand et p ≫ n.

Nombreuses complexités mathématiques et computationnelles.
Très dicile de travailler à n et p xés. Si l'on considère n → +∞,
alors on aura également p → +∞ (a fortiori plus rapidement que n).

M2 Ingénierie Statistique Statistique en Grande Dimension 8 / 97

 Concepts généraux Retour sur la régression linéaire

Modèle de régression linéaire

L'essentiel de ce module est dédié à la régression linéaire en grande dimension, il

convient donc de faire quelques rappels.

Supposons p xé et petit (p ⩽ n), et posons

Y = Xβ + ε

avec Y ∈ Rn (la réponse), X ∈ Rn×p (les covariables), β ∈ Rp (les coecients)

2
et ε ∈ Rn un bruit de variance σ > 0.

Rem : on peut ajouter une colonne de 1 dans X pour tenir compte d'un intercept
β0 , auquel cas p ↔ p + 1.
Tout estimateur OLS de β satisfait

βb ∈ arg minp ∥Y − X β∥2 .

β∈R

Si rg(X ) = p, alors βb = (X T X )−1 X T Y est l'unique solution.

M2 Ingénierie Statistique Statistique en Grande Dimension 9 / 97

 Concepts généraux Retour sur la régression linéaire

Modèle de régression linéaire

Le prédicteur de Y est Yb = X βb = P[X ] Y , où P[X ] = X (X T X )−1 X T est la

matrice de projection orthogonale de Y sur [X ] = Vect{X1 , . . . , Xp }.

Il peut exister plusieurs βb solutions du problème de minimisation (si rg(X ) < p)

mais Yb est unique.

M2 Ingénierie Statistique Statistique en Grande Dimension 10 / 97

 Concepts généraux Retour sur la régression linéaire

Modèle de régression linéaire

L'estimation est sans biais,

E[β]
b =β et E[Yb ] = X β.

La variance de βb est minimale dans la classe des estimateurs linéaires (Thm. de

Gauss-Markov), elle vaut
b = σ 2 (X T X )−1 .
V(β)

L'erreur quadratique moyenne d'estimation qui en découle est

E ∥βb − β∥2 = σ 2 tr((X T X )−1 ).

 

Si le design est orthonormal (X

T
X = Ip ), elle vaut pσ 2 . Donc elle grandit
linéairement avec p!

M2 Ingénierie Statistique Statistique en Grande Dimension 11 / 97

 Concepts généraux Retour sur la régression linéaire

Modèle de régression linéaire

Si de plus le bruit est gaussien, ε ∼ Nn (0, σ 2 In ), alors

βb ∼ Np (β, σ 2 (X T X )−1 ).

Le test de signicativité du k -ème coecient, H0 : “βk = 0” vs H1 : “βk ̸= 0”

(test de Student), peut être conduit grâce au fait que

βbk − βk
√ ∼ t(n − p), pour k = 1, . . . , p ,
σb 2 ck

où ck est le k -ème élément diagonal de (X T X )−1 et b2

σ l'estimateur sans biais
usuel de σ2 .
Rem : si le bruit n'est pas gaussien, on retrouve des résultats similaires mais
asymptotiquement (lorsque n → +∞).

M2 Ingénierie Statistique Statistique en Grande Dimension 12 / 97

 Concepts généraux Retour sur la régression linéaire

Modèle de régression linéaire

Dans ce contexte, on peut sélectionner un modèle (c'est-à-dire détecter les

covariables inuentes) selon des protocoles bien connus :

On conserve un sous-ensemble de covariables marginalement signicatives

qui optimisent un certain critère (ex : algorithmes backward/forward ).

Test de Fisher des modèles emboîtés : les covariables sont considérées dans
leur ensemble et pas seulement marginalement.

Attention au coecient R2 bien connu des statisticiens : il n'est pas adapté

pour comparer des modèles emboîtés (nested models ). [Pourquoi ?]

On fera des rappels sur l'AIC et le BIC dans les slides suivants, également
sur la validation croisée.

M2 Ingénierie Statistique Statistique en Grande Dimension 13 / 97

 Concepts généraux Retour sur la régression linéaire

Modèle de régression linéaire

Exemple sur le dataset eucalyptus (n = 1429, p = 2).

M2 Ingénierie Statistique Statistique en Grande Dimension 14 / 97

 Concepts généraux Retour sur la régression linéaire

Modèle de régression linéaire

Exemple sur le dataset eucalyptus (n = 1429, p = 3).

M2 Ingénierie Statistique Statistique en Grande Dimension 15 / 97

 Concepts généraux Retour sur la régression linéaire

Modèle de régression linéaire

Quelles sont les nouveaux enjeux ? De manière non-exhaustive :

Si p > n, rg(X ) ⩽n∧p =n n'est pas de rang p : l'OLS n'est plus unique.

L'analyse asymptotique lorsque n → +∞ doit être repensée pour intégrer le

fait que p → +∞, plus rapidement que n.
Entre autres eets indésirables, la matrice des covariances empiriques
X T X /n n'est plus able pour estimer les liens entre covariables (centrées).

L'aspect informatique doit être entièrement revu : algorithmes

d'optimisation en grande dimension.

M2 Ingénierie Statistique Statistique en Grande Dimension 16 / 97

 Concepts généraux Retour sur la régression linéaire

Modèle de régression linéaire

L'hypothèse fondamentale des modèles en grande dimension est la suivante : il

peut y avoir beaucoup de covariables mais beaucoup sont peu informatives.

Peu de covariables actives et beaucoup de covariables inactives : c'est la notion de

sparsité (voir slide suivant).

On va s'attarder sur deux approches :

On peut forcer la procédure d'estimation à ne pas utiliser de nombreuses

variables (ex : régressions pénalisées qui peuvent permettre, dans certaines
conditions, d'annuler beaucoup de coecients simultanément).

On peut réduire la dimension en travaillant sur des axes factoriels (ex :

régressions PCR et PLS).

M2 Ingénierie Statistique Statistique en Grande Dimension 17 / 97

 Concepts généraux Retour sur la régression linéaire

Sparsité (modèles parcimonieux)

| {z } | {z } | {z }
not-sparse sparse group-sparse

| {z } | {z } | {z }
sparse-group-sparse fused-sparse pattern-sparse

M2 Ingénierie Statistique Statistique en Grande Dimension 18 / 97

 Concepts généraux Fléau de la dimension

Fléau de la dimension (curse of dimensionality)

Lorsque p grandit...

Le concept de voisinage n'a plus beaucoup de sens car les points sont de
plus en plus isolés : voir slides suivants.

L'accumulation de petites uctuations peut avoir un eet conséquent. De

même, l'accumulation d'évènements rares peut ne pas être rare.

Les temps de calcul dans les procédures d'optimisation peuvent exploser :

l'algorithmie doit être méticuleusement pensée.

La géométrie des espaces de grande dimension est contre-intuitive.

M2 Ingénierie Statistique Statistique en Grande Dimension 19 / 97

 Concepts généraux Fléau de la dimension

Une géométrie peu intuitive

Expérience simple : on génère n = 100 points U1 , . . . , Un uniformément dans

[0, 1]p . Les distances entre les points sont représentées pour diérents p.

On voit que :

La distance au plus proche voisin augmente : c'est attendu.

Les distances s'équilibrent : le voisin le plus proche et le voisin le plus éloigné

sont à des distances similaires : c'est contre-intuitif !

La notion de voisinage a-t-elle encore une quelconque utilité ?

M2 Ingénierie Statistique Statistique en Grande Dimension 20 / 97

 Concepts généraux Fléau de la dimension

Une géométrie peu intuitive

On voudrait trouver n de sorte que pour tout x ∈ [0, 1]p , il existe une observation
Ui telle que d(x, Ui ) ⩽ ϵ. En d'autres termes, combien d'observations faudrait-il
pour garantir que l'hypercube est globalement visité ?

Le volume d'une boule de dimension p et de rayon r >0 vaut

π p/2
Vp (r ) = r p.
Γ(p/2 + 1)
√ √
Rem : avec Γ(3/2) = π/2, Γ(2) = 1 et Γ(5/2) = 3 π/4, que trouve-t-on pour
p = 1, 2, 3 ?
On a aussi, lorsque p → +∞,

2
2p
cp

1 2πer p+1
Vp (r ) ∼ √ ≍ = e− 2 ln p+o(p ln p) .
πp p p (1+p)/2

Vp (r ) → 0 très rapidement avec p : plus qu'exponentiellement vite.

M2 Ingénierie Statistique Statistique en Grande Dimension 21 / 97

 Concepts généraux Fléau de la dimension

Une géométrie peu intuitive

Il doit alors exister des boules Bp de rayon ϵ recouvrant l'hypercube,

n
[
[0, 1]p ⊂ Bp (Ui , ϵ)
i=1

d'où l'on déduit que nVp (ϵ) ⩾ Vol([0, 1]p ) = 1. Il faudrait donc

 2
− 2p
1 √ 2πeϵ p+1
n⩾ = πp ≍ c −p p (1+p)/2 = e 2 ln p+o(p ln p) .
Vp (ϵ) p

Rem : si p = 100 (`petite' grande dimension), n ∼ 1040 !

M2 Ingénierie Statistique Statistique en Grande Dimension 22 / 97

 Concepts généraux Fléau de la dimension

Une géométrie peu intuitive

Soit Cp (x, r ) = Bp (x, r )\Bp (x, fr ) avec 0 ⩽ f ⩽ 1, c'est-à-dire la fraction f de la

sphère située à sa surface. Alors,

Vol(Cp (x, r )) Vp (r ) − Vp (fr )

= = 1 − f p.
Vol(Bp (x, r )) Vp (r )

Rem : fp tend très vite vers 0, le volume de l'hypersphère est essentiellement

concentré à sa surface ! (f = 0.99 ci-dessous).

Morale de ces exemples : attention avec les espaces en grande dimension, les
notions qui nous semblent intuitives ne le sont plus.

M2 Ingénierie Statistique Statistique en Grande Dimension 23 / 97

 Concepts généraux Critères AIC et BIC

Sélection de modèle

Reprenons le modèle de régression linéaire Y = Xβ + ε avec ε ∼ Nn (0, σ 2 In ). On

notera X[m] = Vect{Xk , k ∈ m} la matrice formée des covariables d'un modèle

m ⊂ {1, . . . , p} de coecients β[m] = (βk )k ∈ m .

Le `vrai' modèle m∗ , correspondant à la réalité, s'écrit donc

X
Y = X[m∗ ] β[m∗ ] + ε = βk Xk + ε.
k ∈ m∗

Le prédicteur Yb = X[m∗ ] βb[m∗ ] est obtenu par projection orthogonale de Y sur

l'espace [Xm∗ ] des covariables de m∗ .
On choisit le risque quadratique comme critère de qualité,

h i
2
rm∗ = E X[m∗ ] βb[m∗ ] − X β .

C'est la variance associée à l'estimation du signal Xβ par m∗ .

M2 Ingénierie Statistique Statistique en Grande Dimension 24 / 97

 Concepts généraux Critères AIC et BIC

Sélection de modèle

Le modèle m∗ étant inconnu, il faut le chercher parmi une collection de modèles

candidats M. On peut :
Déterminer rm par projection de Y sur [Xm ] pour tout m ∈ M.
Retenir nalement
m0 ∈ arg min {rm }.
m∈M

Il s'agit d'un oracle : il n'a aucune portée pratique car dépendant de β, il joue un
rôle théorique de référence.

On peut chercher à estimer le risque rm et considérer plutôt

b ∈ arg min {b
m rm }
m∈M

an de comparer les performances de m

b avec celles, idéales, données par m0 .

M2 Ingénierie Statistique Statistique en Grande Dimension 25 / 97

 Concepts généraux Critères AIC et BIC

Sélection de modèle

Lemme

Le risque quadratique associé au modèle m vaut

2
rm = (I − P[Xm ] )X β + dm σ 2
| {z } | {z }
biais variance

où [Xm ] = Vect{Xk , k ∈ m} et dm est sa dimension.

Preuve

L'oracle m0 correspond donc au meilleur compromis biais-variance.

M2 Ingénierie Statistique Statistique en Grande Dimension 26 / 97

 Concepts généraux Critères AIC et BIC

Sélection de modèle

Posons
2
rbm = Yb[m] − Y + (2dm − n)σ 2 .

où Yb[m] = X[m] βb[m] est la projection de Y sur [Xm ].

Lemme

L'estimateur rbm est sans biais pour estimer le risque rm .

Preuve

M2 Ingénierie Statistique Statistique en Grande Dimension 27 / 97

 Concepts généraux Critères AIC et BIC

Sélection de modèle

Dans ce contexte, le choix de modèle b ∈ arg minm ∈ M {b

m rm } se ramène en fait
au critère d'information d'Akaike (AIC),

2
+ 2dm σ 2 .

b AIC ∈ arg min
m Yb[m] − Y
m∈M

Défaut de l'AIC : les gros modèles seront globalement favorisés, car ∥Yb[m] − Y ∥2
peut être très faible en cas d'overtting et la pénalisation n'est pas toujours
susante pour contrebalancer.

M2 Ingénierie Statistique Statistique en Grande Dimension 28 / 97

 Concepts généraux Critères AIC et BIC

Sélection de modèle

Rem : si X[m] est de rang plein (donc en petite dimension), on a dm = |m| ce qui
fait qu'on pénalise par le nombre de covariables incluses dans le modèle.

En poursuivant cette logique, on peut se demander par quelle pénalité pen(m)

remplacer 2dm pour obtenir des performances optimales. Soit

p p 2
pen(m) ∝ dm + −2 ln πm

où (πm )m ∈ M est une probabilité sur M.

On montre que le modèle choisi par cette pénalité atteint des performances
comparables à l'oracle m0 à un terme ∝ − ln πm près (exemple en TD).

M2 Ingénierie Statistique Statistique en Grande Dimension 29 / 97

 Concepts généraux Critères AIC et BIC

Sélection de modèle

Si maintenant l'on adopte le point de vue bayésien


 m ∼ π(m) = πm (sur M)
f = X[m] β[m] | m ∼ π(f | m) (sur [Xm ])
Y |f ,m ∼ Nn (f , σ 2 In ),


alors on obtient pour probabilité a posteriori de chaque modèle m,

∥Y −f ∥2
Z
π(m | Y ) ∝ πm e− 2σ 2 dπ(f | m).
f ∈ [Xm ]

Lemme

Sous réserve de quelques hypothèses de régularité sur la probabilité π(· | m), pour
tous m, m′ ∈ M, lorsque n → +∞,

π(m | Y ) b[m′ ] − Y ∥2 − ∥Y
∥Y b[m] − Y ∥2 + σ 2 (dm′ − dm ) ln n πm
ln = + ln + O(1).
π(m′ | Y ) 2σ 2 πm ′

M2 Ingénierie Statistique Statistique en Grande Dimension 30 / 97

 Concepts généraux Critères AIC et BIC

Sélection de modèle

Cela suggère de choisir

∥Yb[m] − Y ∥2 + σ 2 dm ln n − 2σ 2 ln πm .

b ∈ arg min
m
m∈M

Avec πm = 1/|M|, un a priori uniforme sur M, on retrouve le critère

d'information bayésien (BIC),

∥Yb[m] − Y ∥2 + σ 2 dm ln n .

b BIC ∈ arg min
m
m∈M

Les gros modèles sont plus pénalisés par le BIC que par l'AIC.

La notion de pénalisation est très importante pour la suite du cours. Lorsqu'on

pénalise fortement les gros modèles, on ne retient que peu de covariables parmi la
multitude de la grande dimension.

M2 Ingénierie Statistique Statistique en Grande Dimension 31 / 97

 Concepts généraux Validation croisée

Validation croisée (cross-validation)

Pour limiter l'overtting (trop s'attacher aux données d'apprentissage et perdre

ses capacités prédictives), une stratégie indispensable est la validation croisée.

Le principe en est :

Splitter le dataset en un ensemble d'apprentissage (training set ) et un

ensemble de validation (validation set ).

Les données de validation ne sont pas utilisées dans l'apprentissage.

M2 Ingénierie Statistique Statistique en Grande Dimension 32 / 97

 Concepts généraux Validation croisée

Validation croisée (cross-validation)

Beaucoup de procédures possibles, les deux plus connues sont :

Hold-out CV : le dataset est préalablement séparé en apprentissage et

validation (très sensible au découpage, mais parfois on ne peut pas faire
autrement, cf. séries chronologiques).

V -fold CV : le dataset est séparé en V sous-ensembles, puis chaque

sous-ensemble est successivement retiré pour servir de validation.

Rem : si V = n, c'est un leave-one-out (LOO).

M2 Ingénierie Statistique Statistique en Grande Dimension 33 / 97

 Concepts généraux Validation croisée

Validation croisée (cross-validation)

Considérons une partition {G1 , . . . , GV } des individus et D−k = (Y−k , X−k ) les
données privées du k -ème sous-ensemble Gk .
La procédure V -fold consiste, pour k = 1, . . . , V , à :

Prédire Ybi pour tout i ∈ Gk à l'aide des données D−k .

Calculer l'erreur quadratique moyenne des prédictions de cette étape,

1 X
2
Rbk = Yi − Ybi .
|Gk |
i ∈ Gk

À l'issue des V répétitions, le risque estimé de la procédure est

V
1 X
RbCV = Rbk .
V
k=1

M2 Ingénierie Statistique Statistique en Grande Dimension 34 / 97

 Concepts généraux Validation croisée

Validation croisée (cross-validation)

Si plusieurs modèles sont en concurrence ou si un hyperparamètre doit être

calibré, chaque étape contiendra autant d'évaluations de Rbk qu'il y a de
comparaisons à eectuer.

Le modèle retenu sur cette base sera celui qui minimisera RbCV . Ex : calibration
(tuning ) du Lasso (cf. section suivante).

Défaut principal de la V -fold : dicile à réaliser si V et/ou le nombre de

comparaisons K sont trop élevés :

Il faut envisager VK estimations et nVK prédictions.

Tout dépend du temps nécessaire à l'estimation : on ne traitera pas de la

même manière une régression linéaire et un réseau de neurones.

La procédure est accélérée si V est petit, mais on risque de perdre l'eet

stabilisateur de l'échantillonnage répété.

Compromis à trouver entre le temps de calcul et la stabilité de la procédure. En

pratique : V entre 5 et 10 (si n est susant).

M2 Ingénierie Statistique Statistique en Grande Dimension 35 / 97

 Régressions pénalisées

Table des matières

1 Concepts généraux

2 Régressions pénalisées
Best-Subset
Lasso
Ridge et Elastic-Net
Group-Lasso et autres généralisations du Lasso

3 Régressions PCR et PLS

4 Tests multiples

M2 Ingénierie Statistique Statistique en Grande Dimension 36 / 97

 Régressions pénalisées

Préambule technique

Pour tout a ⩾ 1, on dénit

p
!1a
X
a
|β|a = |βk | .
k=1

Il s'agit d'une norme (généralement appelée ℓa ) si a ⩾ 1, alors que si 0 <a<1 ce

n'est pas une norme. [Pourquoi ?]

La norme innie ℓ∞ est dénie par |β|∞ = max{|βk |, k = 1, . . . , p}, elle vérie

|β|∞ = lim |β|a .

a → +∞

On note |β|0 = |supp(β)|, soit le nombre de coordonnées non-nulles de β. Ce

n'est pas une norme [Pourquoi ?] mais par analogie on la nommera ℓ0 .

M2 Ingénierie Statistique Statistique en Grande Dimension 37 / 97

 Régressions pénalisées Best-Subset

Régressions pénalisées

On se place à nouveau dans le contexte de la régression linéaire Y = Xβ + ε avec

un vecteur de coecients β ∈ Rp qui peut être de grande dimension. On suppose
que les covariables sont normalisées (∥Xk ∥ = 1, k = 1 à p ).
La première approche étudiée consiste à pénaliser de sorte que l'estimation de β
soit sparse (peu de coecients actifs).

Considérons une sélection de modèle basée sur le critère

 2
m
bλ ∈ arg min Yb[m] − Y + λ|m|
m∈M
 2
= arg min min Y − X β + λ|β|0
m ∈ M β : Sβ = m

où Sβ = supp(β) est l'ensemble des coordonnées non-nulles de β (et donc en

particulier |β|0 = |m| si Sβ = m), et λ ⩾ 0.
Similaire à l'AIC et au BIC : plus le modèle est gros, plus il est pénalisé (en
proportion de λ).

M2 Ingénierie Statistique Statistique en Grande Dimension 38 / 97

 Régressions pénalisées Best-Subset

Régressions pénalisées

Dénition

Pour λ ⩾ 0, toute solution du problème

 2
βbλ ∈ arg minp Y − Xβ + λ|β|0
β∈R

est un estimateur Best-Subset.

La question semble résolue... sauf que β 7→ |β|0 n'étant pas convexe, le problème
d'optimisation n'a pas de solution facilement atteignable (il est NP-hard ).

En d'autres termes : soit |m| est petit et on explore toutes les possibilités, soit |m|
est trop grand et il faut trouver une heuristique et/ou une autre approche.

Idée : on relaxe le problème en remplaçant la pénalisation par la norme convexe la

`plus proche' de |β|0 , c'est-à-dire |β|1 .

M2 Ingénierie Statistique Statistique en Grande Dimension 39 / 97

 Régressions pénalisées Lasso

Lasso

Dénition

Pour λ ⩾ 0, toute solution du problème convexe

 2
βbλ ∈ arg minp Y − Xβ + λ|β|1
β∈R

est un estimateur Lasso (Least Absolute Shrinkage and Selection Operator).

La solution existe mais elle n'est pas nécessairement unique. Le prédicteur

X βbλ est bien unique (preuve en TD).

Il n'y a pas d'écriture explicite de βbλ , on va étudier un algorithme

d'approximation (la descente de coordonnées).

Rien ne permet d'armer que le modèle induit par βbλ correspond au modèle
optimal pour la pénalité ℓ0 , mais il existe des garanties théoriques dont on
parlera succinctement.

M2 Ingénierie Statistique Statistique en Grande Dimension 40 / 97

 Régressions pénalisées Lasso

Lasso

Remarque importante : selon les références (livres, cours en ligne, etc.), le Lasso
peut aussi être déni par

2
( )
Y − Xβ
βbλ ∈ arg minp + λ|β|1
β∈R 2

ou même avec n ou 2n à la place de 2 au dénominateur, si les colonnes de X ne

sont pas normalisées.

Il est évident que les procédures restent inchangées, il sura d'interpréter λ

correctement (λ ↔ 2λ ↔ nλ ↔ 2nλ.)

M2 Ingénierie Statistique Statistique en Grande Dimension 41 / 97

 Régressions pénalisées Lasso

Lasso

Même si l'on n'utilisera pas tout ce bagage théorique dans le module, il est
important de discuter quelque peu des propriétés des fonctions convexes.

Tout ceci est la base des solutions apportées à ces nouvelles problématiques !

On rappelle qu'une application F : Rp → R est convexe si pour tous x, y ∈ Rp et

a ∈ [0, 1],
F (ax + (1 − a)y ) ⩽ aF (x) + (1 − a)F (y ).
On dénit le sous-gradient ∂F de F au point x par

∂F (x) = w ∈ Rp : ∀ y ∈ Rp , F (y ) ⩾ F (x) + ⟨w , y − x⟩ .


En particulier si x∗ minimise F, alors 0 ∈ ∂F (x ∗ ). [Pourquoi ?]

Si F est convexe et dérivable en x, alors ∂F (x) = {∇F (x)} et pour tout y ∈ Rp ,

F (y ) ⩾ F (x) + ⟨∇F (x), y − x⟩.

De ce qui précède, on voit que si x∗ minimise F et que F est dérivable en x∗ alors

∇F (x ∗ ) = 0. Et si F n'est pas dérivable en x∗ ?
M2 Ingénierie Statistique Statistique en Grande Dimension 42 / 97
 Régressions pénalisées Lasso

Lasso

Un exemple convexe de grande importance : Fa (β) = |β|a pour a ⩾ 1.

Lemme

Le sous-gradient de Fa au point β ∈ Rp est donné par

∂|β|a = w ∈ Rp : ⟨w , β⟩ = |β|a

et |w |b ⩽ 1

où b est l'exposant conjugué de a (1/a + 1/b = 1).

Preuve

M2 Ingénierie Statistique Statistique en Grande Dimension 43 / 97

 Régressions pénalisées Lasso

Lasso

Pour a = 1, c'est la pénalisation du Lasso.

Corollaire

Le sous-gradient de F1 au point β ∈ Rp est donné par

∂|β|1 = w ∈ Rp : wk = sgn(βk )

si βk ̸= 0 et wk ∈ [−1, 1] si βk = 0 .

M2 Ingénierie Statistique Statistique en Grande Dimension 44 / 97

 Régressions pénalisées Lasso

Lasso

Ainsi, dans le cas du Lasso, le sous-gradient de la fonction convexe à minimiser

(notons-la ϕ pour simplier) est donné par

∂ϕ(β) = − 2X T (Y − X β) + λw : w ∈ ∂|β|1 .


Comme βbλ minimise ϕ, 0 ∈ ∂ϕ(βbλ ) d'où l'on déduit l'existence d'un b ∈ Rp

w qui
satisfait

λ
b = 2X T (Y − X βbλ )
λw soit encore X T X βbλ = X T Y − w
b.
2

De plus, w
b est tel que, pour k = 1, . . . , p ,

w
bk = sgn(βbλ, k ) si βbλ, k ̸= 0
bk ∈ [−1, 1]
w sinon.

M2 Ingénierie Statistique Statistique en Grande Dimension 45 / 97

 Régressions pénalisées Lasso

Lasso

De la relation X T X βbλ = X T Y − λw
b /2, on tire

D λ E
0 ⩽ ∥X βbλ ∥2 = βbλ , X T Y − w b
2
X  λ 
= βbλ, k XkT Y − sgn(βbλ, k )
2
k ∈m
bλ

où m
b λ = supp(βbλ ) est le modèle engendré par βbλ .
Il s'ensuit que tout choix de λ ⩾ 2|X T Y |∞ conduit à βbλ = 0.
Dans le cas contraire, il n'y a pas de formule explicite, sauf si le design est
orthogonal (donc en petite dimension). On traitera ce cas en TD, en particulier
les opérateurs de hard-thresholding et de soft-thresholding.

Mais... certaines coordonnées sont exactement nulles ! C'est donc adapté à la

sélection de variables.

M2 Ingénierie Statistique Statistique en Grande Dimension 46 / 97

 Régressions pénalisées Lasso

Lasso

En eet, par dualité de Lagrange, on montre que pour tout λ > 0, il existe un
rayon rλ > 0 tel que

 2 2
βbλ ∈ arg minp Y − Xβ + λ|β|1 } = arg min Y − Xβ
β∈R β ∈ B1 (rλ )

où B1 (rλ ) = {β ∈ Rp : |β|1 ⩽ rλ } est la boule ℓ1 de rayon rλ .

Évidemment, λ et rλ évoluent en sens contraire : rλ ↑ quand λ↓ et
réciproquement.

La géométrie des boules ℓ1 fait que, pour λ susamment grand, certaines

coordonnées de βbλ seront exactement nulles.

M2 Ingénierie Statistique Statistique en Grande Dimension 47 / 97

 Régressions pénalisées Lasso

Lasso

Problème algorithmique : comment calculer βbλ ?

Par convexité, βk 7→ ϕ(β) est minimisée en tout βk∗ annulant la k -ème dérivée
partielle ∂k ϕ(β), si un tel point existe, sinon en βk∗ = 0.

Lemme

Si les covariables sont normalisées ( ∥Xk ∥ = 1), on a

   
λ X
βk∗ = Rk 1 − avec Rk = XkT Y − βj Xj .
2|Rk | + j ̸= k

Preuve

M2 Ingénierie Statistique Statistique en Grande Dimension 48 / 97

 Régressions pénalisées Lasso

Lasso

Cela donne l'algorithme de la descente de coordonnées :

Fixer β∗ arbitrairement.

Répéter jusqu'à convergence (critère à xer) :

Pour k = 1, . . . , p ,
   
λ X
βk∗ = Rk 1 − avec Rk = XkT Y − βj∗ Xj .
2|Rk | + j ̸= k

Renvoyer βbλ = β ∗ .
Convergence vers arg minβ ∈ Rp ϕ(β) assurée par la convexité et la séparabilité.

Les coordonnées peuvent être exactement nulles, d'où la sparsité de l'estimation

et la sélection de variables.

Rem : d'autres algorithmes existent, comme FISTA ou LARS (voir les packages
glmnet ou lars pour des mises en pratique).

M2 Ingénierie Statistique Statistique en Grande Dimension 49 / 97

 Régressions pénalisées Lasso

Lasso

On sait donc résoudre le problème à λ xé. Mais comment choisir λ?

On considère un ensemble de λ ∈ Λ, grille de valeurs possibles.

Pour tout λ ∈ Λ, on détermine βbλ et donc le modèle m

bλ qui lui est associé
(contenant les covariables dont le coecient est non-nul).

On choisit λ
b par validation croisée sur cet ensemble de modèles.

Finalement, on peut :

Conserver βbλb comme estimateur de β si seule la sélection de variables nous

intéresse : il sera fortement biaisé en raison de la pénalisation (eet de
shrinkage ) mais son support nous sut.

Réestimer β sur le sous-ensemble de variables, de petite dimension,

sélectionnées par m
b λb pour corriger le biais (Gauss-Lasso).

M2 Ingénierie Statistique Statistique en Grande Dimension 50 / 97

 Régressions pénalisées Lasso

Lasso

Exemple sur le dataset diabetes du package lars (n = 442, p = 10).

M2 Ingénierie Statistique Statistique en Grande Dimension 51 / 97

 Régressions pénalisées Lasso

Lasso

Exemple sur le dataset diabetes du package lars (n = 442, p = 10).

M2 Ingénierie Statistique Statistique en Grande Dimension 52 / 97

 Régressions pénalisées Lasso

Lasso

Une coupure verticale donne les variables sélectionnées pour chaque λ (ex : on
veut les 4 variables les plus inuentes, on retient {3, 9, 4, 7}).
Choix de λ par validation croisée :

M2 Ingénierie Statistique Statistique en Grande Dimension 53 / 97

 Régressions pénalisées Lasso

Lasso

Les garanties théoriques sont très diciles à obtenir dans de tels contextes. Elles
reposent sur de nombreux raisonnements techniques.

De façon très simpliée, sous certaines hypothèses et à condition que λ soit bien
calibré, avec une probabilité proche de 1,

r
|β|0 ln p
βbλ − β ≲ .
| {z } n
estimation

Quelques informations essentielles à en tirer :

Pour que l'estimation soit able, p = eo(n) est la borne supérieure que peut
atteindre le nombre de covariables.

La sparsité dans β fait diminuer les bornes d'erreur (|β|0 petit).

Les constantes cachées dans ≲ peuvent être grandes et nalement remettre

en question la portée réelle de la garantie théorique.

On peut obtenir des bornes similaires pour les erreurs de prédiction et de sélection.

M2 Ingénierie Statistique Statistique en Grande Dimension 54 / 97

 Régressions pénalisées Lasso

Lasso

Une lacune du Lasso : il résiste mal à la présence de multicolinéarité dans X.

Supposons que deux covariables sont identiques (X1 = X2 ) et considérons les
modèles
Y = β1 X1 + ε et Y = α1 X1 + α2 X2 + ε.

Pour un λ donné, soit βbλ un estimateur Lasso dans le premier modèle et, pour
tout u ∈ [0, 1], soit

α
bλ, 1 = u βbλ, 1 et bλ, 2 = (1 − u)βbλ, 1 .
α

Alors,

2 2
Y − βbλ, 1 X1 = Y −α
bλ, 1 X1 − α
bλ, 2 X2 et βbλ 1
= α
bλ 1 .

Il y a une innité de solutions au Lasso dans le modèle 2 ! [Pourquoi ?] Attention

donc à l'interprétation des résultats de l'algorithme d'estimation en présence de
covariables sélectionnées très corrélées.

M2 Ingénierie Statistique Statistique en Grande Dimension 55 / 97

 Régressions pénalisées Ridge et Elastic-Net

Ridge

Il existe de nombreuses autres pénalisations possibles que celle du Lasso.

Dénition

Pour λ ⩾ 0, toute solution du problème convexe

2
+ λ|β|22

βbλ ∈ arg minp Y − Xβ
β∈R

est un estimateur Ridge.

Par rapport au Lasso, le Ridge a principalement un avantage et un inconvénient

majeurs :

βbλ est unique (si λ > 0) et il admet une écriture explicite, le problème est
très facile à résoudre.

L'estimation n'est pas sparse (aucun coecient estimé est nul), il n'est pas
exploitable pour la sélection de variables.

M2 Ingénierie Statistique Statistique en Grande Dimension 56 / 97

 Régressions pénalisées Ridge et Elastic-Net

Ridge

Lemme

Si λ > 0, l'estimateur Ridge est donné de manière unique par

−1
βbλ = X T X + λIp XTY.

Preuve

M2 Ingénierie Statistique Statistique en Grande Dimension 57 / 97

 Régressions pénalisées Ridge et Elastic-Net

Ridge

Exemple sur le dataset diabetes du package lars (n = 442, p = 10).

M2 Ingénierie Statistique Statistique en Grande Dimension 58 / 97

 Régressions pénalisées Ridge et Elastic-Net

Ridge

On remarque bien que, comme attendu, toutes les variables sont incluses dans le
modèle quel que soit λ : inadapté à la sélection de variables.

Choix de λ par validation croisée :

M2 Ingénierie Statistique Statistique en Grande Dimension 59 / 97

 Régressions pénalisées Ridge et Elastic-Net

Elastic-Net

Il existe de nombreuses autres pénalisations possibles que celle du Lasso.

Dénition

Pour λ⩾0 et µ ⩾ 0, toute solution du problème convexe

2
+ λ|β|22 + µ|β|1

βbλ, µ ∈ arg minp Y − Xβ
β∈R

est un estimateur Elastic-Net.

Son utilité principale est de corriger le Lasso lorsqu'il existe des covariables
trop corrélées.
Idée : au lieu de sélectionner une variable parmi les variables corrélées (Lasso), il aura
tendance à sélectionner soit tout le groupe de variables corrélées, soit aucune (grâce au
lissage du Ridge).

La présence de la pénalité ℓ1 rend à nouveau possible la sparsité dans

l'estimation et donc la sélection de variables.

M2 Ingénierie Statistique Statistique en Grande Dimension 60 / 97

 Régressions pénalisées Ridge et Elastic-Net

Elastic-Net

Il est facile de mettre à jour la descente de coordonnées du Lasso pour l'adapter à

l'Elastic-Net. Si les covariables sont normalisées (∥Xk ∥ = 1) :

Fixer β∗ arbitrairement.

Répéter jusqu'à convergence (critère à xer) :

Pour k = 1, . . . , p ,
   
Rk µ X
βk∗ = 1− avec Rk = XkT Y − βj∗ Xj .
1+λ 2|Rk | +
j ̸= k

Renvoyer βbλ, µ = β ∗ .

Les coordonnées peuvent être exactement nulles, d'où la sparsité de l'estimation

et la sélection de variables.

M2 Ingénierie Statistique Statistique en Grande Dimension 61 / 97

 Régressions pénalisées Ridge et Elastic-Net

Elastic-Net

On retrouve le Lasso pour λ=0 et le Ridge pour µ = 0. Mais on peut réécrire le

problème de minimisation sous la forme

 1 − α 
2
βbλ ∈ arg minp Y − Xβ +λ |β|22 + α|β|1 .
β∈R 2

C'est l'approche suivie par glmnet du package du même nom : combinaison du

Lasso (α = 1) et du Ridge (α = 0 et λ ↔ λ/2).

M2 Ingénierie Statistique Statistique en Grande Dimension 62 / 97

 Régressions pénalisées Ridge et Elastic-Net

Elastic-Net

Exemple sur le dataset diabetes du package lars (n = 442, p = 10), α = 0.8.

M2 Ingénierie Statistique Statistique en Grande Dimension 63 / 97

 Régressions pénalisées Ridge et Elastic-Net

Elastic-Net

Une coupure verticale donne les variables sélectionnées pour chaque λ (ex : on
veut les 4 variables les plus inuentes, on retient {3, 9, 4, 7}).
Choix de λ par validation croisée, pour α = 0.8 :

M2 Ingénierie Statistique Statistique en Grande Dimension 64 / 97

 Régressions pénalisées Ridge et Elastic-Net

Sparsité vs Shrinkage

M2 Ingénierie Statistique Statistique en Grande Dimension 65 / 97

 Régressions pénalisées Group-Lasso et autres généralisations du Lasso

Group-Lasso

On s'intéresse désormais à un contexte très répandu où les covariables ont une

structure de groupes.

Lorsque p est très grand, il peut être pertinent de travailler à l'échelle supérieure.
Exemples classiques :

Les gènes peuvent être regroupés en chromosomes.

Les villes peuvent être regroupées en pays, les pays en continents.

Les minutes peuvent être regroupées en heures, les heures en jours.

Suite à un clustering, on peut travailler sur les clusters de variables plutôt

que sur les variables.

M2 Ingénierie Statistique Statistique en Grande Dimension 66 / 97

 Régressions pénalisées Group-Lasso et autres généralisations du Lasso

Group-Lasso

Soit {G1 , . . . , Gg } une partition des variables {1, . . . , p} en g groupes. On note

βG k les coecients de β associés au groupe Gk (k = 1 à g ).
L'hypothèse fondamentale est désormais : il peut y avoir beaucoup de groupes
mais beaucoup sont peu informatifs.

Quelle pénalisation utiliser pour que, telle l'action du Lasso sur les variables, les
groupes dans leur ensemble soient pénalisés ?

Dénition

Pour λ = (λ1 , . . . , λg ) ⩾ 0, toute solution du problème convexe

 g 
2
X
βbλ ∈ arg minp Y − Xβ + λk |βGk |2
β∈R
k=1

est un estimateur Group-Lasso.

M2 Ingénierie Statistique Statistique en Grande Dimension 67 / 97

 Régressions pénalisées Group-Lasso et autres généralisations du Lasso

Group-Lasso

Quelques remarques issues de cette dénition :

En pratique, et surtout pour la procédure de validation croisée, on choisira

λk = λ. On ne peut pas se permettre d'optimiser g hyperparamètres !

Chaque groupe recevant la même pondération, il peut être pertinent qu'ils

soient globalement homogènes (en tailles et ordres de grandeur).

Si g =p et Gk = {k}, on retrouve la pénalisation du Lasso.

Attention : une pénalisation en |βGk |1 au lieu de |βGk |2 ne conduit pas à une

sélection à l'échelle des groupes ! [Pourquoi ?]

Au contraire, avec |βGk |2 , la géométrie sous-jacente engendre des groupes

soit entièrement nuls, soit entièrement non-nuls. C'est donc adapté à la
sélection de groupes.

La procédure d'estimation peut être conduite comme pour la descente de

coordonnées du Lasso : la descente se fait cette fois par blocs.

M2 Ingénierie Statistique Statistique en Grande Dimension 68 / 97

 Régressions pénalisées Group-Lasso et autres généralisations du Lasso

Group-Lasso

Exemple sur le dataset colon du package gglasso (n = 62, p = 100, g = 20).

M2 Ingénierie Statistique Statistique en Grande Dimension 69 / 97

 Régressions pénalisées Group-Lasso et autres généralisations du Lasso

Group-Lasso

Exemple sur le dataset colon du package gglasso (n = 62, p = 100, g = 20).

M2 Ingénierie Statistique Statistique en Grande Dimension 70 / 97

 Régressions pénalisées Group-Lasso et autres généralisations du Lasso

Group-Lasso

On observe que les coecients d'un même groupe s'activent ensemble.

Choix de λ par validation croisée :

M2 Ingénierie Statistique Statistique en Grande Dimension 71 / 97

 Régressions pénalisées Group-Lasso et autres généralisations du Lasso

Autres généralisations du Lasso

Le Sparse-Group-Lasso :

Structure en groupes, chaque groupe est sparse.

La sélection est multi-échelle : des groupes sont sélectionnés, puis les

variables dans ces groupes.

Pour λ = (λ1 , . . . , λg ) ⩾ 0 et µ ⩾ 0,
 g 
2
X
βbλ, µ ∈ arg minp Y − Xβ + λk |βGk |2 + µ|β|1 .
β∈R
k=1

Adapté pour les datasets de très grande dimension qu'on peut clusteriser,
mais en pratique dicile à calibrer.

M2 Ingénierie Statistique Statistique en Grande Dimension 72 / 97

 Régressions pénalisées Group-Lasso et autres généralisations du Lasso

Autres généralisations du Lasso

Le Fused-Lasso :

Les covariables sont ordonnées (dans le temps, dans l'espace, ...), de sorte
qu'on aimerait que βk ≈ βk+1 .
Pour λ ⩾ 0,
 p−1 
2
X
βbλ ∈ arg minp Y − Xβ +λ |βk+1 − βk | .
β∈R
k=1

Il s'agit simplement d'un Lasso dans le modèle reparamétrisé par γ = ∆β .

M2 Ingénierie Statistique Statistique en Grande Dimension 73 / 97

 Régressions pénalisées Group-Lasso et autres généralisations du Lasso

Autres généralisations du Lasso

L'Adaptive-Lasso :

Pour réduire l'eet de shrinkage du Lasso, on peut adapter la pénalisation

en donnant plus de poids aux variables considérées comme inuentes.

Supposons qu'un estimateur βbλ soit déjà disponible (ex : Lasso, Ridge,
Gauss-Lasso, OLS, etc.)

Pour µ ⩾ 0,
 p 
2
X |βk |
βbλ, µ ∈ arg minp Y − Xβ +µ .
β∈R |βbλ, k |
k=1

[Que se passe-t-il si des composantes de βbλ sont nulles ?]

Idée : si βbλ est de qualité raisonnable, on se rapproche du problème initial

pénalisé par ℓ0 . [Pourquoi ?]

M2 Ingénierie Statistique Statistique en Grande Dimension 74 / 97

 Régressions PCR et PLS

Table des matières

1 Concepts généraux

2 Régressions pénalisées

3 Régressions PCR et PLS

Régression PCR
Régression PLS

4 Tests multiples

M2 Ingénierie Statistique Statistique en Grande Dimension 75 / 97

 Régressions PCR et PLS

Changement de paradigme

On ne cherche plus à pénaliser l'estimation pour la contraindre à atteindre la

sparsité, désormais on tente de réduire la dimension par une transformation
préalable des données.

La régression est conduite sur des axes factoriels. Deux approches étudiées :

On fait une ACP sur X à l'issue de laquelle on ne retient que r ≪p

composantes principales Z1 , . . . , Zr , puis on considère Y = Z α + ε.
On construit s≪p variables latentes W1 , . . . , Ws sur la base d'une
corrélation maximale avec Y, puis on considère Y = W α + ε.

M2 Ingénierie Statistique Statistique en Grande Dimension 76 / 97

 Régressions PCR et PLS Régression PCR

Rappels sur l'ACP

Supposons que la matrice X ∈ Rn×p est normalisée (X̄k =0 et ∥Xk ∥ = 1). On

i
note X sa i -ème ligne.

Proposition (ACP normée  1/2)

Soit la décomposition en valeurs singulières (SVD) de X,

r0
X
X = σk uk vkT
k=1

où r0 = rg(X ) et σ1 ⩾ . . . ⩾ σr0 > 0. Pour tout r ⩽ r0 , l'espace vectoriel qui

minimise
n r0
2
X X
S : dim(S) ⩽ r 7→ X i − PS X i = σk2
i=1 k=r +1

est donné par Sr = Vect{v1 , . . . , vr }.

M2 Ingénierie Statistique Statistique en Grande Dimension 77 / 97

 Régressions PCR et PLS Régression PCR

Rappels sur l'ACP

Proposition (ACP normée  2/2)

Pour k = 1, . . . , r , les vecteurs

Zk = Xvk
sont les composantes principales. Ce sont sont des combinaisons linéaires des

variables initiales, et elles sont non-corrélées : ⟨Zk , Zℓ ⟩ = 0 pour k ̸= ℓ.

M2 Ingénierie Statistique Statistique en Grande Dimension 78 / 97

 Régressions PCR et PLS Régression PCR

Régression PCR

Le contexte est donc tout trouvé pour traiter par régression PCR (Principal
Component Regression ) le modèle Y = Xβ + ε :

Par une ACP, réduire les p colonnes de X aux r composantes principales de

Z (intérêt pratique lorsque r ≪ p ).
Considérer le modèle de régression Y = Zα + ε et estimer

b = (Z T Z )−1 Z T Y .
α

Remonter à un estimateur de β grâce à la relation

βb = Vr α
b

où la matrice Vr ∈ Rp×r est formée des v1 , . . . , vr issus de l'ACP (les

loadings ). En eet, comme Z α = XVr α, ce choix est naturel (Z α ≈ X β ).
On peut sélectionner r soit par les méthodes usuelles de l'ACP (screeplot de
l'inertie), soit par validation croisée.

M2 Ingénierie Statistique Statistique en Grande Dimension 79 / 97

 Régressions PCR et PLS Régression PCR

Régression PCR

Exemple sur le dataset Hitters du package ISLR (n = 263, p = 16).

M2 Ingénierie Statistique Statistique en Grande Dimension 80 / 97

 Régressions PCR et PLS Régression PCR

Régression PCR

La méthode paraît très simple et s'exonère des éventuels problèmes de colinéarités

auxquels le Lasso, par exemple, est très sensible.

Néanmoins, un détail d'importance : on n'utilise pas Y dans la réduction de

dimension ! En calculant Vr , seules les covariables de X sont exploitées.
Le fait de négliger les p−r composantes principales restantes repose sur un critère
(l'inertie projetée) intégralement issu de X, et non sur les éventuels eets sur Y.
Est-ce totalement pertinent pour expliquer les liens entre X et Y ?

M2 Ingénierie Statistique Statistique en Grande Dimension 81 / 97

 Régressions PCR et PLS Régression PLS

Régression PLS

Pour corriger cette lacune, considérons la régression PLS (Partial Least Squares
Regression ) an de traiter le modèle Y = Xβ + ε :

On va construire W1 , . . . , Ws , combinaisons linéaires orthogonales des

covariables (intérêt pratique lorsque s ≪ p ).
Contrairement à une simple ACP sur X, on va exploiter Y !

Les loadings de Vs sont déterminés par la relation W = XVs (remarquer

l'analogie avec Z = XVr dans la régression PCR) :
La première composante v1 maximise l'application v 7→ ⟨Y , Xv ⟩2 sous
la contrainte ∥v ∥ = 1. [Analogie avec l'ACP ?]
Pour 2 ⩽ k ⩽ s , les autres composantes maximisent l'application
v 7→ ⟨Y , Xv ⟩2 sous les contraintes ∥v ∥ = 1 et ⟨Xv , Xvℓ ⟩ = 0, ℓ < k .
Les s étapes de cet algorithme sont très facilement menées. La matrice
W ∈ Rn×s ainsi formée contient des composantes orthogonales.

M2 Ingénierie Statistique Statistique en Grande Dimension 82 / 97

 Régressions PCR et PLS Régression PLS

Régression PLS

Le contexte pour la régression PLS est donc :

Construire W selon les contraintes énoncées. Contrairement à la régression

PCR, la réponse Y est exploitée.

Considérer le modèle de régression Y = Wα + ε et estimer

b = (W T W )−1 W T Y .
α

Remonter à un estimateur de β grâce à la relation

βb = Vs α
b

où Vs ∈ Rp×s contient les loadings. En eet, comme W α = XVs α, ce choix

est naturel (W α ≈ X β ).

On peut sélectionner s par validation croisée.

M2 Ingénierie Statistique Statistique en Grande Dimension 83 / 97

 Régressions PCR et PLS Régression PLS

Régression PLS

Exemple sur le dataset Hitters du package ISLR (n = 263, p = 16).

M2 Ingénierie Statistique Statistique en Grande Dimension 84 / 97

 Tests multiples

Table des matières

1 Concepts généraux

2 Régressions pénalisées

3 Régressions PCR et PLS

4 Tests multiples
Contexte statistique
Contrôler le FDR

M2 Ingénierie Statistique Statistique en Grande Dimension 85 / 97

 Tests multiples Contexte statistique

Tests multiples

On termine ce module par un aperçu des tests multiples, c'est-à-dire que l'on sort
du contexte de la régression.

On s'interroge sur la abilité d'une conclusion issue de beaucoup d'expériences

statistiques, chacune munie d'un risque d'erreur.

Comme exemple introductif, supposons que l'expression d'un gène est mesurée
dans une population saine A et dans une population atteinte d'une certaine
maladie B à laquelle on soupçonne ce gène d'être relié.

On peut alors tester H0 : “µA = µB ” vs H1 : “µA ̸= µB ”, où µA et µB sont

l'expression moyenne du gène dans chaque population. Dans ce cas, un rejet de
H0 conrmerait notre intuition.

Sous les hypothèses classiques (normalité et indépendance), un test de Student

pourrait être conduit, calibré pour que P0 (rejeter H0 ) = α.

M2 Ingénierie Statistique Statistique en Grande Dimension 86 / 97

 Tests multiples Contexte statistique

Tests multiples

Mais il est bien trop simpliste de supposer qu'un seul gène puisse être responsable
d'une maladie, au contraire les réseaux de gènes sont tellement complexes qu'il
faudrait évaluer l'expression d'une grande quantité d'entre eux.

Avec les technologies modernes (ou déjà has been, cf. slide 6), on peut mesurer
simultanément l'expression de très nombreux gènes (jusqu'à ∼ 105 voire plus).

Si P0 (rejeter H0, g ) = α pour chaque gène g, à combien de faux-positifs (FP)

doit-on s'attendre ? Si l'on appelle G0 l'ensemble des gènes inactifs,

X
E[FP] = P0 (rejeter H0, g ) = α|G0 |.
g ∈ G0

Sur 10000 gènes inactifs, on pourrait en détecter ∼ 500 à tort ! Procédure

désastreuse. Comment contrôler le risque dans ce contexte de tests multiples ?

M2 Ingénierie Statistique Statistique en Grande Dimension 87 / 97

 Tests multiples Contexte statistique

Uniformité des p-values

Dans un souci de généralisation, considérons des tests paramétriques

H0, i : “θ ∈ Θ0, i ” vs H1, i : “θ ∈ Θ1, i ”

pour i = 1, . . . , m. Chaque test (unilatéral ici) prend la forme 1{Si ⩾ si } , où Si est

une statistique de test et si un seuil de rejet.

La p-value du i -ème test est donnée par

pi = sup Tθ (Si )
θ ∈ Θ0, i

où Tθ (t) = Pθ Z ⩾ t avec Z de même loi que Si . Il s'agit de la probabilité
maximale d'observer une valeur plus extrême que Si lorsque H0, i est vraie.

M2 Ingénierie Statistique Statistique en Grande Dimension 88 / 97

 Tests multiples Contexte statistique

Uniformité des p-values

Lemme

Sous H0, i , la p-value pi est stochastiquement dominée par une variable aléatoire

uniforme, c'est-à-dire que

∀ u ∈ [0, 1], sup Pθ pi ⩽ u ⩽ u.
θ ∈ Θ0, i

Preuve

M2 Ingénierie Statistique Statistique en Grande Dimension 89 / 97

 Tests multiples Contexte statistique

Tests multiples

On considère m tests pour lesquels on dispose des réalisations pb1 , . . . , pbm .

Dénition

Soient Rb ⊂ I = {1, . . . , m} l'ensemble des tests rejetés et I0 ⊂ I l'ensemble des

tests i ∈ I pour lesquels H0, i est vraie. Les ensembles

Rb ∩ I0 et Rb \I0

contiennent respectivement les faux-positifs et les vrais-positifs.

On notera FP = |Rb ∩ I0 | et TP = |Rb \I0 | la taille de ces ensembles.

Idéalement, on devrait baser la construction de Rb sur des arguments garantissant

que FP soit petit et TP grand.

M2 Ingénierie Statistique Statistique en Grande Dimension 90 / 97

 Tests multiples Contexte statistique

Correction de Bonferroni

Cela suggère une première approche très brutale, adaptée seulement lorsque m
n'est pas très grand. Dans le cas contraire, on limite FP mais TP restera petit.

La correction de Bonferroni consiste à choisir

n αo
Rb = i ∈ I : pbi ⩽ .
m

Lemme

Avec la correction de Bonferroni, P(FP > 0) ⩽ α.

Preuve

On appelle FWER cette probabilité (Family Wise Error Rate ).

M2 Ingénierie Statistique Statistique en Grande Dimension 91 / 97

 Tests multiples Contrôler le FDR

Taux FDR

Dénition

Le taux FDR (False Discovery Rate) correspond à la valeur moyenne du taux de

faux-positifs parmi les tests positifs,

 
FP
FDR =E (convention 0/0 = 0).
FP + TP

Une procédure qui contrôlerait le FDR permettrait de maîtriser à la fois FP et TP.

En particulier, on voudrait un numérateur petit et un dénominateur grand.

Basons une telle procédure sur la règle de décision

Rb = {i ∈ I : pbi ⩽ s}.

Comment choisir s pour que FDR ⩽α (numérateur petit) tout en maximisant |R|
b
(dénominateur grand) ?

M2 Ingénierie Statistique Statistique en Grande Dimension 92 / 97

 Tests multiples Contrôler le FDR

Contrôler le FDR

On voit que

X
E[FP] = Pθ (b
pi ⩽ s)
i ∈ I0
X
⩽ pi ⩽ s) ⩽ s|I0 | ⩽ sm.
sup Pθ (b
θ ∈ Θ0, i
i ∈ I0

Si l'on pose s = pb(k) , la k -ème des p-values triées par ordre croissant, on a

|R|
b =k et E[FP] ⩽ mb p(k) , et donc

mb
p(k)
FDR ⩽ .
k
Il s'ensuit qu'un choix de seuil pertinent devrait satisfaire pb(k) ⩽ αk/m pour que,
à son tour, FDR ⩽ α. En conséquence, posons

n αkb o n αk o
Rb = i ∈ I : pbi ⩽ où kb = max k ∈ I : pb(k) ⩽ .
m m

M2 Ingénierie Statistique Statistique en Grande Dimension 93 / 97

 Tests multiples Contrôler le FDR

Contrôler le FDR

Un exemple avec m = 100, la droite est d'équation k 7→ αk/m avec α = 5%.

Seules les p-values sous la droite rejettent leurs tests respectifs.

Cette règle de décision est appelée procédure de Benjamini-Hochberg (BH).

[Pourquoi αk/m et non αk/|I0 | qui est une meilleure majoration ?]

M2 Ingénierie Statistique Statistique en Grande Dimension 94 / 97

 Tests multiples Contrôler le FDR

Contrôler le FDR

Plus généralement, avec la règle de décision

n αβ(k)
b o n αβ(k) o
Rb = i ∈ I : pbi ⩽ où kb = max k ∈ I : pb(k) ⩽
m m
et β : I → R+ une fonction croissante, on a le contrôle suivant.

Proposition

Pour la procédure décrite ci-dessus, on a

α|I0 | X β(j ∧ m)
FDR ⩽ .
m j(j + 1)
j ⩾1

Il semble donc pertinent de choisir β de sorte que

X β(j ∧ m)
⩽ 1.
j(j + 1)
j ⩾1

M2 Ingénierie Statistique Statistique en Grande Dimension 95 / 97

 Tests multiples Contrôler le FDR

Contrôler le FDR

Cela n'implique pas le contrôle FDR ⩽α pour BH où β(k) = k [Pourquoi ?], mais
on montre qu'il est valide sous certaines hypothèses supplémentaires.

En revanche, il est assuré pour Bonferroni où β(k) = 1.

La procédure de Benjamini-Yekutieli (BY) consiste à choisir

m
k X 1
β(k) = où Hm = .
Hm k
k=1

Là encore, le contrôle est assuré.

M2 Ingénierie Statistique Statistique en Grande Dimension 96 / 97

 Tests multiples Contrôler le FDR

Illustration de Benjamini-Hochberg

Exemple sur le dataset golub du package multtest (m = 3051 gènes mesurés

sur nA = 27 patients atteints de leucémie de type ALL et nB = 11 de type AML).

Procédure BH en trait plein, α = 5% et correction de Bonferroni en pointillés

(inadaptées ici). Rejets selon BH en rouge et non-rejets en vert.

M2 Ingénierie Statistique Statistique en Grande Dimension 97 / 97