1 Introduction

\[V(\widehat{\beta})=\sigma^2(X'X)^{-1}\] Donc la précision des estimateurs va diminuer dès que le déterminant de \((X'X)\) s’approche de zéro.

2 Régression Ridge

2.1 Equivalence avec une contrainte sur la norme des coefficients

  • Nous allons que l’estimateur ridge est une solution du problème \[\widetilde{\beta}=\mbox{argmin}_{\beta\in \mathbb{R}^p,\;\|\beta\|\leq \delta} \|Y-X\beta\|^2\]

  • Supposons que l’estimateur \(\widehat{\beta}\) des MC ne satisfait pas la contrainte sinon bien évidemment \(\widetilde{\beta}=\widehat{\beta}\)

  • Si \(\delta\) est « petit » et tel que \(\|\widehat{\beta}\|^2>\delta\), il faut calculer l’estimateur obtenu. Pour cela, introduisons le lagrangien du problème : \[\mathcal{L}(\beta,\tau)= \|Y-X\beta\|^2+\tau (\|\beta^2\|-\delta) \]

  • Le minimum recherché satisfait forcément la condition \[\displaystyle\frac{\partial \mathcal{L}}{\partial \beta}=0\]

  • Donc une solution de notre problème va satisfaire le système suivant \[\left\{\begin{array}{rcl} -2X'(Y-X\widetilde{\beta}_{\mbox{ridge}})-2\widetilde{\tau}\widehat{\beta}_{\mbox{ridge}}&=&0\\ \|\widetilde{\beta}_{\mbox{ridge}}\|-\delta&=&0 \end{array} \right. \]

  • On peut vérifier que le couple \(\widehat{\beta}_{\mbox{ridge}}=(X'X+\kappa I_p)^{-1}X'Y\) et \[\widetilde{\tau}=(\widehat{\beta}_{\mbox{ridge}}X'Y-\widehat{\beta}_{\mbox{ridge}}'X'X\widehat{\beta}_{\mbox{ridge}})/\|\widehat{\beta}_{\mbox{ridge}}\|^2 \] est bien un minimum du problème.

2.2 Propriétés statistiques de l’estimateur ridge

Proposition Si \(\widehat{\beta}_{\mbox{ridge}}=(X'X+\kappa I_p)^{-1}X'Y\), alors, \[ \mathbb{E}\left(\widehat{\beta}_{\mbox{ridge}}\right)=\beta − \kappa(X'X + \kappa I_p)^{−1}\beta.\] \[V\left(\widehat{\beta}_{\mbox{ridge}}\right)=\sigma^2 (X'X + \kappa I_p)^{−1} (X'X)^{-1}(X'X + \kappa I_p)^{−1}.\]

Remarques

  • L’estimateur ridge est biaisé,

  • Comme la variance fait intervenir \((X'X + \kappa I_p)^{−1}\), donc la variance de l’estimateur Ridge est plus faible.

Proposition (Comparaison des EQM) Soit \(r\) le nombre de valeurs propres non-nulles de \(X'X\), alors \[ \mbox{tr}[\mbox{EQM}(\widehat{\beta})]=\sigma^2\displaystyle\sum_{j=1}^p \displaystyle\frac{1}{\lambda_j}\] \[ \mbox{tr}[\mbox{EQM}(\widehat{\beta}_{\mbox{ridge}})]=\displaystyle\sum_{j=1}^r \displaystyle\frac{\sigma^2\lambda_j+\kappa[P'\beta]_j^2}{(\lambda_k\kappa)^2}\] Et si \(\kappa\leq \displaystyle\frac{\sigma^2}{\beta'\beta}\) alors, \[ \mbox{tr}[\mbox{EQM}(\widehat{\beta}_{\mbox{ridge}})]\leq \mbox{tr}[\mbox{EQM}(\widehat{\beta})]\]

Remarque : la régression ridge peut être plus précise (dans l’estimation des paramètres) que la régression ordinaire, au sens de la trace de l’EQM. Cependant, la condition \(\kappa\leq \displaystyle\frac{\sigma^2}{\beta'\beta}\) dépend de paramètres inconnus \(\beta\) et \(\sigma^2\) et elle n’est donc pas utilisable pour choisir une valeur de \(\kappa\).

2.3 La régression ridge en pratique

2.3.1 Centrage et réduction

  • En général, la valeur de \(\beta\) dépend de l’échelle de mesure de la variable explicative associée : \(\beta\) sera différent si la variable est mesurée en gramme ou en kilo.

  • Rappelons que la régression ridge contraint la norme au carré de \(\beta\) (\(\|\beta\|^2\)) à être inférieure à une valeur \(\delta\). Lors du calcul de la norme, afin de ne pas pénaliser ou favoriser un coefficient, il est souhaitable que chaque coefficient soit affecté de manière « semblable ». Une manière de réaliser cet équilibre consiste à centrer et réduire toutes les variables.

  • A la différence de la régression classique, où les variables sont en général conservées telles que mesurées, il est d’usage de centrer et réduire les variables explicatives.

  • Une variable centrée-réduite, \(\widetilde{X}_j\) issue de la variable \(X_j\) s’écrit \[\widetilde{X}_j=(X_j-\overline{x}_j\mathbb{1})/\widehat{\sigma}_j\]

  • La matrice \(\widetilde{X}\) contiendra donc des variables centrées réduites.

  • Il est aussi d’usage de centrer la variable \(Y\) \[Y-\overline{y}\mathbb{1}\] Donc le terme constant (intercept) sera complètement inutile

  • Une fois l’estimation Ridge des paramètres a été effectuée on retrouve les valeur ajustées de \(Y\) de la façon suivante \[\widehat{Y}_{\mbox{ridge}}(\kappa)=\widehat{\sigma}\left(\widetilde{X}\widehat{\beta}_{\mbox{ridge}}\right)+\overline{y}\mathbb{1}\]

2.3.2 Choix de \(\widetilde{\kappa}\)

  • Il faut choisir la valeur « optimum » de \(\kappa\), valeur notée \(\widetilde{\kappa}\) (ou la valeur de \(\delta\)).

  • En général cette étape est pratiquement impossible à réaliser a priori. La valeur \(\widetilde{\kappa}\) sera choisie grâce aux données, elle sera donc stochastique

2.3.2.1 Méthode graphique

  • Une première méthode consiste à tracer un diagramme d’évolution des coefficients \(\widehat{\beta}_{\mbox{ridge}}\) en fonction de \(\kappa\).

  • Un diagramme similaire existe, utilisant non plus en abscisses \(\kappa\), mais le nombre effectif de paramètres \[\mbox{tr}\left(\widetilde{X}(\widetilde{X}'\widetilde{X}+\kappa I_p)^{-1}\widetilde{X}'\right)=\displaystyle\sum_{j=1}^p \displaystyle\frac{d_j^2}{d_j^2+\kappa}.\] où \(d_j^2\) représente la je valeur propre de \(\widetilde{X}'\widetilde{X}.\)

  • La valeur de \(\widetilde{\kappa}\) est alors choisie comme la valeur la plus petite avant laquelle tous les coefficients « plongent » vers 0.

2.3.2.2 Critères analytiques

Il y en a qui estiment \(\kappa\) de la façon suivante

\[ \widetilde{\kappa}=\displaystyle\frac{p\widehat{\sigma}^\ast}{(\widehat{\beta}^\ast)'\widehat{\beta}^\ast}\]

où \(\widehat{\beta}^\ast\) et \(\widehat{\sigma}^\ast\) sont les estimateurs classiques des MC en considérant le jeux des données \((\widetilde{Y},\widetilde{X})\)

2.3.2.3 Apprentissage-validation

  • Elle consiste à séparer de manière aléatoire les données en deux parties distinctes \((X_a,\,Y_a)\) (apprentissage) et \((X_v,\,Y_v)\) (validation).

  • Le jeu d’apprentissage est centré-réduit,

  • Une régression ridge est conduite avec le jeu d’apprentissage \((X_a, Y_a)\) pour toutes les valeurs de \(\kappa\) possibles : on choisit une grille de valeurs pour κ, comprises entre 0 et une valeur maximale.

  • On utilise ensuite tous ces modèles et les variables explicatives \(X_v\), les valeurs de la variable à expliquer sont prédites \(\widetilde{Y}_{v}^{\mbox{ridge}}(\kappa)\) pour tous les \(\kappa\).

  • La qualité du modèle est ensuite obtenue en mesurant la distance entre les observations prévues et les vraies observations par un critère. Le plus courant est le PRESS

\[\mbox{PRESS}=\| \widetilde{Y}_{v}^{\mbox{ridge}}(\kappa)- Y_v\|^2.\]

  • Le coefficient optimal \(\widetilde{\kappa}\) choisi est celui qui conduit à la minimisation du critère choisi.

  • Cette procédure semble la plus indiquée mais elle nécessite beaucoup de données puisqu’il en faut suffisamment pour estimer le modèle, mais il faut aussi beaucoup d’observations dans le jeu de validation \((X_v, Y_v)\) pour bien évaluer la capacité de prévision dans de nombreux cas de figure.

2.3.3 Mise en oeuvre sous R

  • Nous sommes en présence de biscuits non cuits pour lesquels on souhaite connaître rapidement et à moindre coût, la composition en quatre ingrédients : les lipides, les sucres, la farine et l’eau.

  • Des méthodes classiques de chimie analytique permettent de mesurer la composition des biscuits mais elles sont assez longues et coûteuses et ne peuvent pas être mises en ligne sur une chaîne de production

  • Il serait souhaitable de pouvoir les remplacer par la mesure d’un spectre d’absorbance dans le domaine proche infrarouge (ou spectre proche infrarouge). Pour savoir si cela est possible, nous allons devoir essayer d’expliquer la composition par le spectre.

  • Nous avons \(n_a = 40\) biscuits non cuits sur lesquels sont mesurés les spectres proches infrarouges : on mesure l’absorbance à une longueur d’onde donnée, pour toutes les longueurs d’ondes entre 1100 et 2498 nanomètres et régulièrement espacées de 2 nanomètres.

> length(seq(1100,2498,by=2))
[1] 700

Nous avons donc 700 variables potentiellement explicatives.

  • Ensuite, pour chaque biscuit, on mesure sa composition par les méthodes traditionnelles. Ici nous allons nous intéresser uniquement au pourcentage de sucres.

  • Nous avons donc \(p = 700\) variables pour \(n_a = 40\) individus. Nous sommes bien dans le cas où l’estimateur des moindres carrés classiques \((X'X)^{−1}X'Y\) n’est pas défini, puisque le rang de \(X'X\) vaut ici 40 et non pas \(p = 700\).

  • Nous souhaitons savoir si l’on peut vraiment expliquer le taux de sucres par le spectre proche infrarouge, nous disposons d’un échantillon de validation pour comparer les méthodes :

    • Cet échantillon de validation comporte \(n_v = 32\) individus et ne sera jamais utilisé pour estimer les coefficients d’un modèle quel qu’il soit. Il sert uniquement à comparer une méthode avec une autre et à connaître, pour une méthode, sa capacité de prévision
    • Cette séparation en deux échantillons de tailles 40 et 32 fait partie du jeu de données et nous ne nous poserons donc pas la question de cette répartition.

  • Les données peuvent être chargées dans le package fds.

> library(fds)
Loading required package: rainbow
Loading required package: MASS
Loading required package: pcaPP
Loading required package: RCurl
Loading required package: bitops
> data(labp)
> data(labc)
> data(nirp)
> data(nirc)
> 
> Xbrut.app <- t(nirc$y) ## La matrice des variables explicatives, échantillon d'apprentissage. 
> colnames(Xbrut.app)=paste("spect_",nirc$x,sep="")
> Ybrut.app <- t(labc) ## La matrices des variables réponses: les variables mesurant le sucre dans les biscruits.
> colnames(Ybrut.app)=rownames(labc)
> Xbrut.val <- t(nirp$y)
> colnames(Xbrut.val)=paste("spect_",nirp$x,sep="")
> Ybrut.val <- t(labp)
> colnames(Ybrut.val)=rownames(labc)
> 
> 
> Yselec <- 2 ## Variable à expliquer
> cookie.app <- cbind.data.frame(Ybrut.app[,Yselec],Xbrut.app)
> colnames(cookie.app)[1] <- "sucres"
> 
> cookie.val <- cbind.data.frame(Ybrut.val[,Yselec],Xbrut.val)
> colnames(cookie.val)[1] <- "sucres"

  • Comme cela est l’usage, la régression ridge sous R centre et réduit toutes les variables explicatives. Elle centre aussi la variable à expliquer mais ne la réduit pas.

  • Nous allons estimer \(\kappa\) et \(\beta\) par \(\widetilde{\kappa}\) et \(\widehat{\beta}_{\mbox{ridge}}\) sur le jeu de données d’apprentissage regroupant \(n_a = 40\) individus et \(p = 700\) variables explicatives.

  • Dans cet exemple nous allons utiliser la validation croisée et diviser les 40 observations en 4 parties de 10 individus, de manière aléatoire.

  • Cette séparation sera toujours la même quelles que soient les méthodes et elle est effectuée en utilisant une fonction du package pls.

> library(pls)

Attaching package: 'pls'
The following object is masked from 'package:stats':

    loadings
> set.seed(87) ## Pour que les valeurs tirées d'une facçon aléatoire soient toujours les mêmes.
> cvseg <- cvsegments(nrow(cookie.app),k=4,type="random")
> cvseg
$V1
 [1]  1 38 31  6 27 28 22 29 17 23

$V2
 [1] 12  3 13 21 16 40  8 18 33 20

$V3
 [1]  2  7 14 34 26 10 15 37 36 19

$V4
 [1] 32 11 25 39  5 24 30  4  9 35

attr(,"incomplete")
[1] 0
attr(,"type")
[1] "random"

  • Nous choisissons un ensemble de valeurs possibles régulièrement espacées pour \(\kappa\) entre \(0\) et \(\kappa_{\mbox{max}}\)

  • Pour chaque valeur de \(\kappa\) nous avons donc un estimateur \(\widehat{\beta}_{\mbox{ridge}}\) calculé sur toutes les observations sauf celle de la \(i^{ième}\) partie.

  • Ensuite nous calculons le PRESS sur les observations de la \(i^{\ième}\) partie. Ces PRESS sont ensuite additionnés pour obtenir le PRESS de validation croisée et nous déduisons la valeur \(\widetilde{\kappa}\) qui minimise le PRESS.

  • Pour ce faire on créé la fonction suivante

> library(MASS)
>  choix.kappa <- function(kappamax,cvseg,nbe=1000) {
+    ## kappamax: le kappa maximal à considérer
+    ## cvseg: sont les échantillons de la cross-validation
+    ## nbe: le nombre de kappa à considérer 
+    press <- rep(0,nbe)
+    for (i in 1:length(cvseg)){
+      valid <- cookie.app[unlist(cvseg[i]),] ## Echantillon de validation 
+      modele <- lm.ridge(sucres~.,data = cookie.app[unlist(cvseg[-i]),],lambda=seq(0,kappamax,length=nbe)) ## Estimation du modeèle de Ridge en faisant varier le kappa de 0 à kappamax et en considérant nbe valeurs.
+  coeff <- coef(modele) ## On récupère les coefficients d'estimation ridge
+  prediction <- matrix(coeff[,1],nrow(coeff),nrow(valid))+ coeff[,-1]%*%t(data.matrix(valid[,-1])) ## 
+  press <- press+rowSums((matrix(valid[,1],nrow(coeff),nrow(valid),byrow=T)-prediction)^2)
+  }
+  kappaet <- seq(0,kappamax,length=nbe)[which.min(press)] ## On considère le plus petit PRESS parmi les quatre échantillons de la cross validation.
+  return(list(kappaet=kappaet,press=press))
+  }
  • Nous pouvons donc regarder l’évolution du PRESS en fonction de \(\kappa\) et choisir la valeur \(\widetilde{\kappa}\) par validation croisée.
> res <- choix.kappa(1,cvseg) ## kappamax est choisie égale  à 1
> kappaet <- res$kappaet
> dt=data.frame(seq(0,1,length=1000),res$press)
> colnames(dt)=c("kappa","PRESS")
> 
> library(ggplot2)
> p<-ggplot(dt,aes(x=kappa,y=PRESS))+geom_line()+xlab(expression(kappa))
> p<-p+geom_vline(xintercept = res$kappaet)
> p+annotate(geom="text",x=c(0.02,res$kappaet,0.95),y=c(res$press[2],res$press[which.min(res$press)],res$press[995]),label=c("MCO","Ridge","Constante"),col="red")+theme_bw()

  • Donc la valeur optimale pour le \(\kappa\) est 0.1081081.

  • Nous allons maintenant comparer l’Erreur quadratrique moyenne de prévision entre le modèle de Ridge et le modèle MCO classique.

> coeff <- coef(lm.ridge(sucres~.,data = cookie.app,
+ lambda=kappaet))
> prediction.ridge <- rep(coeff[1],nrow(cookie.val))+as.vector(coeff[-1]%*% t(data.matrix(cookie.val[,-1])))
> eqmp.ridge<-mean((cookie.val[,1]-prediction.ridge)^2) ## EQMP pour le modèle de Ridge
> eqmp.ridge
[1] 4.951633
> modele.lm=lm(sucres~.,data = cookie.app)
> prediction.mco<-predict(modele.lm,newdata=cookie.val)
> eqmp.lm<-mean((cookie.val[,1]-prediction.mco)^2) ##
> eqmp.lm
[1] 134.5219
  • Ainsi l’EQMP pour le modèle MCO est 134.521903 et le modèle Ridge 4.9516331 ce qui montre que la régression ridge apport une amélioration considérable à la régression linéaire ordinaire. La prévision par proche infrarouge du taux de sucres semble assez satisfaisante, à condition de bien choisir la méthode de régression.

3 Régression sur composantes : PCR et PLS

3.1 Régression sur composantes principales (PCR)

3.1.1 Calcul de l’estimateur PCR.

  • L’écriture classique du modèle de régression \[ Y = X_1\beta_1 + \ldots + X_p\beta_p + \epsilon.\]

  • Nous souhaitons introduire des composantes, c’est-à-dire changer de variable et trouver un modèle équivalent : \[ Y = X_1^\ast\beta_1 + \ldots + X_k^\ast\beta_p + \epsilon.\]

  • Comme \(X'X\) est une matrice symétrique alors, elle peut être diagonalisable de la façon suivante \[X'X=P\Lambda P'\] où
    • \(P\) est la matrice des vecteurs \(P\) normalisée, donc \(P\) est matrice orthogonale, \((P')^{-1}=P\).
    • \(\Lambda\) est la matrice diagonale des valeurs propres \(\lambda_1,\,\lambda_2,\ldots,\lambda_p\).
    • Comme \(X'X\) est semi-définie positive, donc \(\lambda_1\geq \lambda_2\geq \ldots\geq \lambda_p\geq 0\).

  • Si l’on effectue l’analyse en composantes principales (ACP) du tableau \(X\) (ou du triplet \((X, I_p, I_n/n)\)), la matrice \(P\) est la matrice des axes principaux normés à l’unité, mais les valeurs propres de l’ACP sont les \(\{λ_j\}\) avec \(j\) variant de \(1\) à \(p\) divisés par \(n.\)

  • En remplaçant \(X\) par \(XPP'\) nous avons \[Y = XPP'\beta + \epsilon\]

  • Donc \[Y = X^\ast\beta^\ast + \epsilon\]

où \(X^\ast=XP\) et \(\beta^\ast=P'\beta\)

  • Les colonnes de \(X^\ast\) sont les composantes principales. On a aussi \[(X^\ast)'X^\ast=(XP)'XP=P'X'XP=P'P\Lambda P'P=\Lambda\] Cela signifie que les nouvelles variables \(X^\ast_j = XP_j\) constituant les colonnes de \(X^\ast\), sont orthogonales entre elles et de norme \(\lambda_j\).

  • Le but de la régression sur composantes principales consiste à ne conserver qu’une partie des composantes principales, à l’image de ce qui est fait en ACP. Les \(k\) composantes principales conservées seront la part conservée de l’information contenue dans les variables explicatives, alors que les \((p − k)\) éliminées seront la part d’information contenue dans les variables explicatives qui sera éliminée, car considérée comme négligeable.

  • L’information est mesurée en terme d’inertie ou de dispersion et est égale à la valeur propre : plus la valeur propre \(\lambda_j\) est élevée, plus la part d’information apportée par la composante \(j\) est importante. On doit donc conserver que les composantes dont la part d’information associée est grande, à savoir conserver les composantes associées aux \(k\) premières valeurs propres. Les estimateurs des coefficients des \(k\) premières composantes principales retenues seront les moins variables.

  • Si nous souhaitons conserver seulement une composante, la première composante principale sera sélectionnée. La matrice \(\X^ast\) ne sera composée que de la variable \(X^\ast_1\) et nous obtenons l’estimateur suivant \[\widehat{\beta}^\ast_1=\left((X^\ast_1)'X^\ast_1\right) (X^\ast_1)'Y\]. Et la variance de \(\widehat{\beta}^\ast_1\) sera \[ \mbox{Var}(\widehat{\beta}^\ast_1)=\sigma^2\,\left((X^\ast_1)'X^\ast_1\right)^{-1}=\displaystyle\frac{\sigma^2}{\lambda_1}\]

  • Si nous souhaitons conserver les \(k\) premières composantes \[\widehat{\beta}^\ast(k)=\left((X^\ast_{[1:k]})'X^\ast_{[1:k]}\right) (X^\ast_{[1:k]})'Y\].

  • Les variables \(X^\ast_j\) étant orthogonales, nous obtenons

\[\mbox{Cov}\left(\widehat{\beta}^\ast_j(k),\widehat{\beta}^\ast_{j'}(k)\right)=\left\{ \begin{array}{lll} \displaystyle\frac{\sigma^2}{\lambda_1} & \mbox{ si } & j=j'\\ 0 & \mbox{ sinon } & \end{array} \right. \]

  • Si \(k=p\), \(\widehat{\beta}^\ast=\widehat{\beta}^\ast(p)\) alors \[\mbox{V}(\widehat{\beta}^\ast)=\sigma^2\Lambda^{-1}\]

  • L’estimateur \(\widehat{\beta}^\ast\) minimise les moindres carrés puisque les moindres carrés du modèle étoile et du modèle initial sont identiques par construction \[\|Y-X\beta\|^2=\|Y-XPP'\beta\|^2=\|Y-X^\ast\beta^\ast\|^2\]

3.1.2 Choix du nombre de composantes

3.1.2.1 Méthode graphique.

Pour déterminer \(k\), il est possible, à l’image de ce qui est fait en ACP, de tracer le diagramme en tuyaux d’orgue des valeurs propres et de choisir le numéro \(k\) de la valeur propre après laquelle les valeurs propres sont nettement plus petites. En général, cette procédure est adaptée à l’interprétation (c’est-à-dire à l’ACP), mais sélectionne trop peu de composantes pour un modèle utilisé à des fins de prévision.

3.1.2.2 Apprentissage-validation.

  • Elle consiste à séparer de manière aléatoire les données en deux parties distinctes \((X_a,\,Y_a)\) (apprentissage) et \((X_v,\,Y_v)\) (validation).

  • Le jeu d’apprentissage est centré-réduit,

  • Une régression ridge est conduite avec le jeu d’apprentissage \((X_a, Y_a)\) pour toutes les valeurs de \(\kappa\) possibles : on choisit une grille de valeurs pour \(\kappa\), comprises entre 0 et une valeur maximale.

  • On utilise ensuite tous ces modèles et les variables explicatives \(X_v\), les valeurs de la variable à expliquer sont prédites \(\widetilde{Y}^{\mbox{PCR}}_{v}(k)\) pour tous les \(k\).

  • La qualité du modèle est ensuite obtenue en mesurant la distance entre les observations prévues et les vraies observations par un critère. Le plus courant est le PRESS

\[\mbox{PRESS}=\| \widetilde{Y}^{\mbox{PCR}}_{v}(k)- Y_v\|^2.\]

  • Le coefficient optimal \(\widetilde{k}\) choisi est celui qui conduit à la minimisation du critère choisi.

  • Cette procédure semble la plus indiquée mais elle nécessite beaucoup de données puisqu’il en faut suffisamment pour estimer le modèle, mais il faut aussi beaucoup d’observations dans le jeu de validation \((X_v, Y_v)\) pour bien évaluer la capacité de prévision dans de nombreux cas de figure.

3.1.2.3 Validation croisée

  • Pour toutes les valeurs de \(k\) possibles (\(k\) variant de \(1\) à \(K\) fixé, avec \(K \leq \mbox{rang}(X)\)), on supprime une observation (ou un groupe de \(b\) observations) puis on estime le modèle sans cette (ou ces) observation(s).

  • On peut alors prévoir cette (ou ces) observation(s) grâce à ce modèle estimé.

  • Dans le cas d’une seule observation enlevée, la \(i^{\mbox{ième}}\), pour un nombre de composantes \(k\), la prévision est notée \(\widetilde{Y}^{\mbox{PCR}}_{v}(k)\)

  • On peut enfin à l’aide d’un critère, par exemple le PRESS, connaître la capacité de prévision d’un modèle à \(k\) composantes par \[\mbox{PRESS}=\| \widetilde{Y}^{\mbox{PCR}}_v(\kappa)- Y_v\|^2.\]

  • Le nombre optimal \(k\) de composantes est celui qui réalise le minimum du PRESS.

3.1.3 Mise en oeuvre sous R.

  • Nous reprenons l’exemple de la prévision du taux de sucres par un spectre proche infrarouge (700 variables explicatives).

  • Afin d’utiliser la régression sur composantes principales, nous devons déterminer le nombre de composantes à retenir. Ce nombre \(k\) sera toujours déterminé par validation croisée sur 4 groupes de 10 observations. Rappelons la méthode proposée par le package pls.

> library(pls)
> set.seed(87) ## Pour que les valeurs tirées d'une facçon aléatoire soient toujours les mêmes.
> cvseg <- cvsegments(nrow(cookie.app),k=4,type="random")
> cvseg
$V1
 [1]  1 38 31  6 27 28 22 29 17 23

$V2
 [1] 12  3 13 21 16 40  8 18 33 20

$V3
 [1]  2  7 14 34 26 10 15 37 36 19

$V4
 [1] 32 11 25 39  5 24 30  4  9 35

attr(,"incomplete")
[1] 0
attr(,"type")
[1] "random"

  • La régression sur composantes principales est conduite simplement grâce à la fonction pcr

  • Nous pouvons avoir au maximum 40 composantes principales (\(\min(n_a, p) = n_a = 40\))

> n.app <- nrow(cookie.app) ## Taille de l'échantillon d'apprentissage
> stdX.app <- sqrt(apply(cookie.app[,-1],2,var)*(n.app-1)/n.app) ## On calcul les eécart-types des variables.
  • La modélisation PCR
> modele.pcr <- pcr(sucres~.,ncomp=30,data = cookie.app,scale=stdX.app,validation = "CV",segments=cvseg)
> msepcv.pcr <- MSEP(modele.pcr,estimate=c("train","CV")) ## Erreur Quadratique de Prediction.
> x=c(msepcv.pcr$val[1,,],msepcv.pcr$val[2,,])
> y=c(rep("train",length(msepcv.pcr$val[1,,])),rep("cv",length(msepcv.pcr$val[2,,])))
> z=c(0:29,0:29)
> dt=data.frame(x,y,z)
> colnames(dt)=c("msep","sample","comps")
> p<-ggplot(dt,aes(x=comps,y=msep,col=sample))+geom_line()
> p+theme_bw()

  • Le nombre de composantes \(k\) est trouvé numériquement par
> ncomp.pcr <- which.min(msepcv.pcr$val["CV",,])-1
> ncomp.pcr
6 comps 
      6
  • La prévision par le modèle des observations du jeu de validation est obtenue par
> modele.pcr.fin <- pcr(sucres~.,ncomp=ncomp.pcr,data = cookie.app,
+ scale=stdX.app)
> ychap <- predict(modele.pcr.fin,newdata=cookie.val)[,1,ncomp.pcr]
> res.pcr <-cookie.val[,"sucres"]-ychap
> mean(res.pcr^2)
[1] 1.027865

4 Régression aux moindres carrés partiels (PLS)

4.1 Description de l’algorithme

  • A l’image de la régression sur composantes principales, nous sommes intéressés par de nouvelles variables explicatives \(t^{(1)},\, t^{(2)}, \ldots, t^{(k)}\), combinaisons linéaires des variables de départ, \(t^{(j)} = Xc_j\), qui soient orthogonales entre elles et classées par ordre d’importance.

  • Le choix de ces composantes \(t^{(j)}\) doit être dicté, non pas par la part de variabilité qu’elles représentent parmi les variables explicatives originales (comme en régression sur composantes principales), mais par leur lien avec la variable à expliquer.

  • La régression PLS se fait alors d’une façon itérative:

  • 1ère étape, le tableau \(X\) est noté \(X^{(1)}\) et \(Y\) noté \(Y^{(1)}\). La première composante PLS \(t^{(1)}\in \mathbb{R}^n\) est choisie telle que \[t^{(1)}= \mbox{argmax}_{t=X^{(1)}w,\,w\in \mathbb{R}^p,\,\|w\|=1} <t,Y^{(1)}>\] Ensuite nous effectuons la régression univariée de \(Y^{(1)}\) sur \(t^{(1)}\). \[Y^{(1)}= r_1t^{(1)} +\epsilon_1 \]

où \(r_1\in \mathbb{R}^n\) est le coefficient de la régression estimé par MC et \(\epsilon_1=(I- P_{t^{(1)}})Y^{(1)}\) sont les résidus de la régression simple sans constante.

  • 2-ième étape, Soit \(Y^{(2)}=\epsilon_1=(I- P_{t^{(1)}})Y^{(1)}\) la partie non encore expliquée de \(Y\). Soit \(X^{(2)}=(I- P_{t^{(1)}})X^{(1)}\)la partie de \(X^{(1)}\) n’ayant pas encore servi à expliquer. La seconde composante PLS est choisie telle que

\[t^{(2)}= \mbox{argmax}_{t=X^{(2)}w,\,w\in \mathbb{R}^p,\,\|w\|=1} <t,Y^{(2)}>\] Ensuite nous effectuons la régression univariée de \(Y^{(2)}\) sur \(t^{(2)}\). \[Y^{(2)}= r_2t^{(2)} +\epsilon_2 \]

où \(r_2\in \mathbb{R}^n\) est le coefficient de la régression estimé par MC et \(\epsilon_2=(I- P_{t^{(2)}})Y^{(2)}\) sont les résidus de la régression simple sans constante.

Et on continue ainsi jusqu’à la k-ième étape.

  • La régression PLS cherche donc une suite de composantes PLS qui soient orthogonales entre elles et cela par construction. Puisque \(t^{(j)}\) est une combinaison linéaire des colonnes de \(X^{(j)}\), qui est par construction dans l’orthogonal de \(\mathcal{I}(t^{(1)}, \ldots , t^{(j−1)})\), alors \(t^{(j)}\) sera bien orthogonale à \(t^{(1)}, \ldots , t^{(j−1)}\).

  • Ces composantes sont choisies comme maximisant la covariance (empirique) entre \(Y\) et une composante \(t\) quand \(X\) et \(Y\) sont centrées au préalable.

4.2 Mise en oeuvre sous R.

  • Nous reprenons l’exemple de la prévision du taux de sucres par un spectre proche infrarouge

  • Pour la régression PLS, et à l’image de la régression sur composantes principales, nous devons déterminer le nombre \(k\) de composantes PLS.

  • Il sera déterminé par validation croisée sur 4 groupes de 10 observations. En utilisant le package pls nous avons la modélisation PLS jusqu’au nombre maximal de \(K = 28\) composantes.

> modele.pls <- plsr(sucres~.,ncomp=30,data = cookie.app,scale= stdX.app,validation = "CV",segments=cvseg)
  • Le choix du nombre de composantes est réalisé graphiquement par
> msepcv.pls <- MSEP(modele.pls,estimate=c("train","CV")) ## Erreur Quadratique de Prediction.
> x=c(msepcv.pls$val[1,,],msepcv.pls$val[2,,])
> y=c(rep("train",length(msepcv.pls$val[1,,])),rep("cv",length(msepcv.pls$val[2,,])))
> z=c(0:29,0:29)
> dt=data.frame(x,y,z)
> colnames(dt)=c("msep","sample","comps")
> p<-ggplot(dt,aes(x=comps,y=msep,col=sample))+geom_line()
> p+theme_bw()

  • Ce minimum est obtenu simplement par
> ncomp.pls <- which.min(msepcv.pls$val["CV",,])-1
> ncomp.pls 
10 comps 
      10
  • Nous pouvons ensuite évaluer la capacité de prévision par le MSEP sur notre jeu de validation
> modele.pls.fin <- plsr(sucres~.,ncomp=ncomp.pls,data=cookie.app,scale=stdX.app)
> ychap <- predict(modele.pls.fin,newdata=cookie.val)[,1,ncomp.pls]
> res.pls <- cookie.val[,"sucres"]-ychap
> mean(res.pls^2)
[1] 3.997247

5 Plus de lecture et de liens