Chapitre 4 Sélection de variables et de modèles

4.1 Introduction

Ce chapitre présente des principes, outils et méthodes très généraux pour choisir un « bon » modèle. Nous allons principalement utiliser la régression linéaire pour illustrer les méthodes en supposant que tout le monde connaît ce modèle de base. Les méthodes présentées sont en revanche très générales et peuvent être appliquées avec n’importe quel autre modèle (régression logistique, arbres de classification et régression, réseaux de neurones, analyse de survie, etc.)

L’expression « sélection de variables » fait référence à la situation où l’on cherche à sélectionner un sous-ensemble de variables à inclure dans notre modèle à partir d’un ensemble de variables \(X_1, \ldots, X_p\). Le terme variable ici inclut autant des variables distinctes que des transformations d’une ou plusieurs variables.

Par exemple, supposons que les variables age, sexe et revenu sont trois variables explicatives disponibles. Nous pourrions alors considérer choisir entre ces trois variables. Mais aussi, nous pourrions considérer inclure \(\mathrm{age}^2\), \(\mathrm{age}^3\), \(\log(\mathrm{age})\), etc. Nous pourrions aussi considérer des termes d’interactions entre les variables, comme \(\mathrm{age} \cdot \mathrm{revenu}\) ou \(\mathrm{age}\cdot\mathrm{revenu}\cdot\mathrm{sexe}\). Le problème est alors de trouver un bon sous-ensemble de variables parmi toutes celles considérées.

L’expression « sélection de modèle » est un peu plus générale. D’une part, elle inclut la sélection de variables car, pour une famille de modèles spécifiques (régression linéaire par exemple), choisir un sous-ensemble de variables revient à choisir un modèle. D’autre part, elle est plus générale car elle peut aussi faire référence à la situation où l’on cherche à trouver le meilleur modèle parmi des modèles de natures différentes. Par exemple, on pourrait choisir entre une régression linéaire, un arbre de régression, une forêt aléatoire, un réseau de neurones, etc.

4.2 Sélection de variables et de modèles selon les buts de l’étude

Nous disposons d’une variable réponse \(Y\) et d’un ensemble de variables explicatives \(X_1, \ldots, X_p\). L’attitude à adopter dépend des buts de l’étude.

1e situation : On veut développer un modèle pour faire des prédictions sans qu’il soit important de tester formellement les effets des paramètres individuels.

Dans ce cas, on désire seulement que notre modèle soit performant pour prédire des valeurs futures de \(Y\). On peut alors baser notre choix de variable (et de modèle) en utilisant des outils qui nous guiderons quant aux performances prédictives futures du modèle (voir \(\mathsf{AIC}\), \(\mathsf{BIC}\) et validation croisée plus loin). On pourra enlever ou rajouter des variables et des transformations de variables au besoin afin d’améliorer les performances prédictives. Les méthodes que nous allons voir concernent essentiellement ce contexte.

2e situation : On veut développer un modèle pour estimer les effets de certaines variables sur notre \(Y\) et tester des hypothèses de recherche spécifiques concernant certaines variables.

Dans ce cas, il est préférable de spécifier le modèle dès le départ selon des considérations scientifiques et de s’en tenir à lui. Faire une sélection de variables dans ce cas est dangereux car on ne peut pas utiliser directement les valeurs-p des tests d’hypothèses (ou les intervalles de confiance sur les paramètres) concernant les paramètres du modèle final car elles ne tiennent pas compte de la variabilité due au processus de sélection de variables.

Une bonne planification de l’étude est alors cruciale afin de collecter les bonnes variables, de spécifier le ou les bons modèles, et de s’assurer d’avoir suffisamment d’observations pour ajuster le ou les modèles désirés.

Si procéder à une sélection de variables est quand même nécessaire dans ce contexte, il est quand même possible de le faire en divisant l’échantillon en deux. La sélection de variables pourrait être alors effectuée avec le premier échantillon. Une fois qu’un modèle est retenu, on pourrait alors réajuster ce modèle avec le deuxième échantillon (sans faire de sélection de variables cette fois-ci). L’inférence sur les paramètres (valeurs-p, etc.) sera alors valide. Le désavantage ici qu’il faut avoir une très grande taille d’échantillon au départ afin d’être en mesure de le diviser en deux.

4.3 Mieux vaut plus que moins

Il est préférable d’avoir un modèle un peu trop complexe qu’un modèle trop simple. Plaçons-nous dans le contexte de la régression linéaire et supposons que le vrai modèle est inclus dans le modèle qui a été ajusté. Il y a donc des variables en trop dans le modèle qui a été ajusté: ce dernier est dit surspécifié.

Par exemple, supposons que le vrai modèle est \(Y=\beta_0+\beta_1X_1+\varepsilon\) mais que c’est le modèle \(Y=\beta_0+\beta_1X_1+\beta_2X_2+\varepsilon\) qui a été ajusté. Dans ce cas, règle générale, les estimateurs des paramètres et les prédictions provenant du modèle sont sans biais. Mais leurs variances estimées seront un peu plus élevées car on estime des paramètres pour des variables superflues.

Supposons à l’inverse qu’il manque des variables dans le modèle ajusté et que le modèle ajusté est sous-spécifié. Par exemple, supposons que le vrai modèle est \(Y=\beta_0+\beta_1X_1+\beta_2X_2+\varepsilon\), mais que c’est le modèle \(Y=\beta_0+\beta_1X_1+\varepsilon\) qui est ajusté. Dans ce cas, généralement, les estimateurs des paramètres et les prédictions sont biaisés.

Ainsi, il est généralement préférable d’avoir un modèle légèrement surspécifié qu’un modèle sous-spécifié. Plus généralement, il est préférable d’avoir un peu trop de variables dans le modèle que de prendre le risque d’omettre une ou plusieurs variables importantes. Il faut faire attention et ne pas tomber dans l’excès et avoir un modèle trop complexe (avec trop de variables inutiles) car il pourrait souffrir de surajustement (over-fitting). Les exemples qui suivent illustreront ce fait.

4.4 Trop beau pour être vrai

Cette section traite de l’optimisme de l’évaluation d’un modèle lorsqu’on utilise les mêmes données qui ont servies à l’ajuster pour évaluer sa performance. Un principe fondamental lorsque vient le temps d’évaluer la performance prédictive d’un modèle est le suivant : si on utilise les mêmes observations pour évaluer la performance d’un modèle que celles qui ont servi à l’ajuster (à estimer le modèle et ses paramètres), on va surestimer sa performance. Autrement dit, notre estimation de l’erreur que fera le modèle pour prédire des observations futures sera biaisée à la baisse. Ainsi, il aura l’air meilleur que ce qu’il est en réalité. C’est comme si on demandait à un cinéaste d’évaluer son dernier film. Comme c’est son film, il n’aura généralement pas un regard objectif. C’est pourquoi on aura tendance à se fier à l’opinion d’un critique.

On cherchera donc à utiliser des outils et méthodes qui nous donneront l’heure juste (une évaluation objective) quant à la performance prédictive d’un modèle.

4.5 Principes généraux

Les idées présentées ici seront illustrées à l’aide de la régression linéaire. Par contre, elles sont valides dans à peu près n’importe quel contexte de modélisation.

Plaçons-nous d’abord dans un contexte plus général que celui de la régression linéaire. Supposons que l’on dispose de \(n\) observations indépendantes sur (\(Y, X_1, \ldots, X_p\)) et que l’on a ajusté un modèle \(\widehat{f}(X_1, \ldots, X_p)\), avec ces données, pour prédire une variable continue \(Y\).

Ce modèle peut être un modèle de régression linéaire, \[\begin{align*} \widehat{f}(X_1, \ldots, X_p) = \widehat{\beta}_0 + \widehat{\beta}_1X_1 + \cdots + \widehat{\beta}_pX_p \end{align*}\] mais il pourrait aussi avoir été construit selon d’autres méthodes (réseau de neurones, arbre de régression, forêt aléatoire, etc.) Une manière de quantifier la performance prédictive du modèle est l’erreur quadratique moyenne de généralisation (generalization mean squared error), \[\begin{align*} \mathsf{EQMG}=\mathsf{E}\left[\left\{(Y-\widehat{f}(X_1, \ldots, X_p)\right\}^2\right] \end{align*}\] lorsque (\(Y, X_1, \ldots, X_p\)) est choisi au hasard dans la population. Cette quantité mesure l’erreur théorique (la différence au carré entre la vraie valeur de \(Y\) et la valeur prédite par le modèle) que fait le modèle en moyenne pour l’ensemble de la population. Plus cette quantité est petite, meilleur est le modèle. Le problème est que l’on ne peut pas la calculer car on n’a pas accès à toute la population. Tout au plus peut-on essayer de l’estimer ou bien d’estimer une fonction qui, sans l’estimer directement, classifiera les modèles dans le même ordre qu’elle.

Une première idée est d’estimer l’erreur quadratique moyenne de l’échantillon d’apprentissage (training mean squared error), \[\begin{align*} \widehat{\mathsf{EQM}}_a= \frac{1}{n}\sum_{i=1}^n \left\{Y_i-\widehat{f}(X_{i1}, \ldots, X_{ip})\right\}^2. \end{align*}\]

Malheureusement, selon le principe fondamental de la section précédente, cette quantité n’est pas un bon estimateur de l’\(\mathsf{EQMG}\). En effet, comme on utilise les mêmes observations que celles qui ont estimé le modèle, l’\(\widehat{\mathsf{EQM}}_a\) aura tendance à toujours diminuer lorsqu’on augmente la complexité du modèle (par exemple, lorsqu’on augmente le nombre de paramètres). L’\(\widehat{\mathsf{EQM}}_a\) tend à surestimer la qualité du modèle en sous-estimant l’\(\mathsf{EQMG}\) et le modèle a l’air meilleur qu’il ne l’est en réalité.

4.5.1 Choix d’un modèle polynomial en régression linéaire

Cet exemple simple servira à illustrer le fait qu’on ne peut utiliser directement les mêmes données qui ont servi à ajuster un modèle pour évaluer sa performance.

Nous disposons de 100 observations sur une variable cible \(Y\) et d’une seule variable explicative \(X\). Le fichier selection1_train contient les données. Nous voulons considérer des modèles polynomiaux (en \(X\)) afin d’en trouver un bon pour prédire \(Y\). Un modèle polynomial est un modèle de la forme \(Y=\beta_0 + \beta_1X+\cdots+\beta_kX^k+\varepsilon\). Le cas \(k=1\) correspond à un modèle linéaire simple, \(k=2\) à un modèle cubique, \(k=3\) à un modèle cubique, etc. Notre but est de déterminer l’ordre (\(k\)) du polynôme qui nous donnera un bon modèle. Voici d’abord le graphe de ces 100 observations de l’échantillon d’apprentissage.

Ces données ont été obtenues par simulation et le vrai modèle sous-jacent (celui qui a généré les données) est le modèle cubique, c’est-à-dire le modèle d’ordre \(k=3\).

J’ai ajusté tour à tour à tour les modèles polynomiaux jusqu’à l’ordre 10, avec l’échantillon d’apprentissage de taille 100. C’est-à-dire, le modèle linéaire avec un polynôme d’ordre \(k=1\) (linéaire), \(k=2\) (quadratique), etc., jusqu’à \(k=10\). J’ai ensuite obtenu la valeur de l’erreur quadratique moyenne d’apprentissage pour chacun de ces modèles. En pratique, on ne pourrait pas calculer l’erreur quadratique moyenne de généralization, mais j’ai approximé cette dernière en simulant 100 000 observations du vrai modèle (selection1_test), en obtenant la prédiction pour chacune de ces 100 000 observations en utilisant le modèle d’ordre \(k\) ajusté sur les données d’apprentissage et en calculant l’erreur moyenne quadratique par la suite. En pratique, on ne peut réaliser cette opération car on ne connaît pas le vrai modèle.

Voici le graphe de l’erreur moyenne quadratique d’apprentissage (\(\widehat{\mathsf{EQM}}_a\)) et de l’erreur moyenne quadratique théorique (\(\mathsf{EQMG}\)) en fonction de l’ordre (\(k\)) du modèle utilisé.

On voit clairement que l’\(\widehat{\mathsf{EQM}}_a\) diminue en fonction de l’ordre sur l’échantillon d’apprentissage: plus le modèle est complexe, plus l’erreur observée sur l’échantillon d’apprentissage est petite. La courbe \(\mathsf{EQM}\) donne l’heure juste, car il s’agit d’une estimation de la performance réelle des modèles sur de nouvelles données. On voit que le meilleur modèle est donc le modèle cubique (\(k=3\)), ce qui n’est pas surprenant puisqu’il s’agit du modèle que utilisé pour générer les données. On peut aussi remarquer d’autres éléments intéressants. Premièrement, on obtient un bon gain en performance (\(\mathsf{EQM}\)) en passant de l’ordre \(2\) à l’ordre \(3\). Ensuite, la perte de performance en passant de l’ordre \(3\) à \(4\), et ensuite à des ordres supérieurs n’est pas si sévère, même si elle est présente. Cela illustre empiriquement qu’il est préférable d’avoir un modèle un peu trop complexe que d’avoir un modèle trop simple. Il serait beaucoup plus grave pour la performance de choisir le modèle avec \(k=2\) que celui avec \(k=4\).

En pratique par contre, on n’a pas accès à la population : les 100 000 observations qui ont servi à estimer l’\(\mathsf{EQM}\) théorique ne seront pas disponible. Si on a seulement l’échantillon d’apprentissage, soit 100 observations dans notre exemple, comment faire alors pour choisir le bon modèle? C’est ce que nous verrons à partir de la section suivante.

Mais avant cela, nous allons discuter un peu plus en détail au sujet de la régression linéaire et d’une mesure très connue, le coefficient de détermination (\(R^2\)). Supposons que l’on a ajusté un modèle de régression linéaire \[\begin{align*} \widehat{f}(X_1, \ldots, X_p) = \widehat{Y}=\widehat{\beta}_0 + \widehat{\beta}_1X_1+ \cdots + \widehat{\beta}_p X_p. \end{align*}\] La somme du carré des erreurs (\(\mathsf{SCE}\)) pour notre échantillon est \[\begin{align*} \mathsf{SCE}=\sum_{i=1}^n (Y_i - \widehat{\beta}_0 - \widehat{\beta}_1X_1 - \cdots - \widehat{\beta}_p X_p)^2 = \sum_{i=1}^n (Y_i-\widehat{Y}_i)^2. \end{align*}\] On peut démontrer que si on ajoute une variable quelconque au modèle, la valeur de la somme du carré des erreurs va nécessairement baisser. Il est facile de se convaincre de cela. En régression linéaire, les estimations sont obtenues par la méthode des moindres carrés qui consiste justement à minimiser la \(\mathsf{SCE}\). Ainsi, en ajoutant une variable \(X_{p+1}\) au modèle, la \(\mathsf{SCE}\) ne peut que baisser car, dans le pire des cas, le paramètre de la nouvelle variable sera \(\widehat{\beta}_{p+1}=0\) et on retombera sur le modèle sans cette variable. C’est pourquoi, la quantité \(\widehat{\mathsf{EQM}}_a=\mathsf{SCE}/n\) ne peut être utilisée comme outil de sélection de modèles en régression linéaire.

Nous venons d’ailleurs d’illustrer cela avec notre exemple sur les modèles polynomiaux. En effet, augmenter l’ordre du polynôme de \(1\) revient à ajouter une variable. Le coefficient de détermination (\(R^2\)) est souvent utilisé, à tort, comme mesure de qualité du modèle. Il peut s’interpréter comme étant la proportion de la variance de \(Y\) qui est expliquée par le modèle.

Le coefficient de détermination est \[\begin{align*} R^2=\{\mathsf{cor}(\boldsymbol{y}, \widehat{\boldsymbol{y}})\}^2 = 1-\frac{\mathsf{SCE}}{\mathsf{SCT}}, \end{align*}\] où \(\mathsf{SCT}=\sum_{i=1}^n (Y_i-\overline{Y})^2\) est la somme des carrés totale calculée en centrant les observations. La somme des carrés totale, \(\mathsf{SCT}\), ne varie pas en fonction du modèle. Ainsi, on voit que le \(R^2\) va méchaniquement augmenter lorsqu’on ajoute une variable au modèle (car la \(\mathsf{SCE}\) diminue). C’est pourquoi on ne peut pas l’utiliser comme outil de sélection de variables.

Le problème principal que nous avons identifié jusqu’à présent afin d’être en mesure de bien estimer la performance d’un modèle est le suivant : si on utilise les mêmes observations pour évaluer la performance d’un modèle que celles qui ont servi à l’ajuster, on va surestimer sa performance.

Il existe deux grandes approches pour contourner ce problème lorsque le but est de faire de la sélection de variables ou de modèle :

  • utiliser les données de l’échantillon d’apprentissage (en échantillon) et pénaliser la mesure d’ajustement (ici \(\widehat{\mathsf{EQM}}_a\)) pour tenir compte de la complexité du modèle (par exemple, à l’aide de critères d’informations).
  • tenter d’estimer l’\(\mathsf{EQM}\) directement sur d’autres données (hors échantillon) en utilisant des méthodes de rééchantillonnage, notamment la validation croisée ou la validation externe (division de l’échantillon).

4.6 Pénalisation et critères d’information

Plaçons-nous dans le contexte de la régression linéaire pour l’instant. Nous avons déjà utilisé les critères \(\mathsf{AIC}\) et \(\mathsf{BIC}\) en analyse factorielle. Il s’agit de mesures qui découlent d’une méthode d’estimation des paramètres, la méthode du maximum de vraisemblance (maximum likelihood).

Il s’avère que les estimateurs des paramètres obtenus par la méthode des moindres carrés en régression linéaire sont équivalents à ceux provenant de la méthode du maximum de vraisemblance si on suppose la normalité des termes d’erreurs du modèle. Ainsi, dans ce cas, nous avons accès aux \(\mathsf{AIC}\) et \(\mathsf{BIC}\), deux critères d’information définis pour les modèles dont la fonction objective est la vraisemblance (qui mesure la probabilité des observations sous le modèle postulé suivant une loi choisie par l’utilisateur). La fonction de vraisemblance \(\mathcal{L}\) et la log-vraisemblance \(\ell\) mesurent l’adéquation du modèle.

Supposons que nous avons ajusté un modèle avec \(p\) paramètres en tout (incluant l’ordonnée à l’origine). En régression linéaire, le critère d’information d’Akaike, \(\mathsf{AIC}\), est \[\begin{align*} \mathsf{AIC} &=-2 \ell(\widehat{\boldsymbol{\beta}}, \widehat{\sigma}^2) +2p=n \ln (\mathsf{SCE}) - n\ln(n) + 2p, \end{align*}\] tandis que le critère d’information bayésien de Schwartz, \(\mathsf{BIC}\), est défini par \[\begin{align*} \mathsf{BIC} &=-2 \ell(\widehat{\boldsymbol{\beta}}, \widehat{\sigma}^2) + p\ln(n)=n \ln (\mathsf{SCE}) - n\ln(n) + p\ln(n) \end{align*}\] Plus la valeur du \(\mathsf{AIC}\) (ou du \(\mathsf{BIC}\)) est petite, meilleur est l’adéquation. Que se passe-t-il lorsqu’on ajoute un paramètre à un modèle? D’une part, la somme du carré des erreurs va méchaniquement diminuer, et donc la quantité \(n \ln (\mathsf{SCE}/n)\) va diminuer. D’autre part, la valeur de \(p\) augmente de \(1\). Ainsi, le \(\mathsf{AIC}\) peut soit augmenter, soit diminuer, lorsqu’on ajoute un paramètre; idem pour le \(\mathsf{BIC}\). Par exemple, le \(\mathsf{AIC}\) va diminuer seulement si la baisse de la somme du carré des erreurs est suffisante pour compenser le fait que le terme \(2p\) augmente à \(2 (p+1)\).

Ces critères pénalisent l’ajout de variables afin de se prémunir contre le surajustement. De plus, le \(\mathsf{BIC}\) pénalise plus que le \(\mathsf{AIC}\). Par conséquent, le critère \(\mathsf{BIC}\) va choisir des modèles contenant soit le même nombre, soit moins de paramètres que le \(\mathsf{AIC}\).

Les critères \(\mathsf{AIC}\) et \(\mathsf{BIC}\) peuvent être utilisés comme outils de sélection de variables en régression linéaire mais aussi beaucoup plus généralement avec d’autres méthodes basées sur la vraisemblance (analyse factorielle, régression logistique, etc.) En fait, n’importe quel modèle dont les estimateurs proviennent de la méthode du maximum de vraisemblance produira ces quantités. Nous donnerons des formules générales pour le \(\mathsf{AIC}\) et le \(\mathsf{BIC}\) dans le chapitre sur la régression logistique.

Le critère \(\mathsf{BIC}\) est le seul de ces critères qui est convergent. Cela veut dire que si l’ensemble des modèles que l’on considère contient le vrai modèle, alors la probabilité que le critère \(\mathsf{BIC}\) choisissent le bon modèle tend vers 1 lorsque \(n\) tend vers l’infini. Il faut mettre cela en perspective : il est peu vraisemblable que \(Y\) ait été généré exactement selon un modèle de régression linéaire, car le modèle de régression n’est qu’une approximation de la réalité. Certains auteurs trouvent que le \(\mathsf{BIC}\) est quelquefois trop sévère (il choisit des modèles trop simples) pour les tailles d’échantillons finies. Dans certaines applications, cette parcimonie est utile, mais il n’est pas possible de savoir d’avance lequel de ces deux critères (\(\mathsf{AIC}\) et \(\mathsf{BIC}\)) sera préférable pour un problème donné.

Avant de revenir à l’exemple, voici la description d’une modification du coefficient de détermination, le \(R^2\) ajusté, qui permet (contrairement au \(R^2\)) de faire de la sélection de variables. En régression linéaire, le \(R^2\) ajusté est \[\begin{align*} R^2_a=1-\frac{\mathsf{SCE}/(n-p)}{\mathsf{SCT}/(n-1)}. \end{align*}\] Lorsqu’on ajoute une variable, la somme du carré des erreurs (\(\mathsf{SCE}\)) diminue mais c’est aussi le cas de la quantité \((n-p)\). Ainsi, le \(R^2\) ajusté peut soit augmenter, soit diminuer lorsqu’on ajoute une variable. On peut donc l’utiliser pour choisir le modèle. Plus \(R^2_a\) est élevé, mieux c’est. Ce critère est moins sévère que le \(\mathsf{AIC}\), Ainsi, en général, il va choisir un modèle avec le même nombre ou bien avec plus de paramètres que le \(\mathsf{AIC}\). Pour résumer, on aura la situation suivante : \[ \#(\mathsf{BIC}) \leq \#(\mathsf{AIC}) \leq \#(R^2_a),\] où \(\#\) représente le nombre de paramètres du modèle linéaire.

Il est facile d’obtenir les quantités \(R^2_a\), \(\mathsf{AIC}\) et \(\mathsf{BIC}\) avec la procédure glmselect dans SAS. Le fichier selection1_intro.sas contient les programmes. La sortie qui suit provient des commandes :

proc glmselect data=multi.selection1_train;
model y=x x*x x*x*x /selection=none;
run;

Il s’agit du modèle cubique (d’ordre 3) en \(x\).

Le tableau qui suit résume ces quantités pour tous les modèles de l’ordre 1 à l’ordre 10.

Tableau 4.1: Mesures d’adéquation du modèle linéaire et estimés de l’erreur
\(\mathsf{EQM}\) \(\widehat{\mathsf{EQM}}_a\) \(R^2\) \(R^2_a\) \(\mathsf{AIC}\) \(\mathsf{BIC}\)
1 3191 3674 0.65 0.65 1111 1119
2 3133 2879 0.73 0.72 1088 1099
3 2697 2620 0.75 0.74 1081 1094
4 2767 2582 0.75 0.74 1081 1097
5 2771 2581 0.75 0.74 1083 1102
6 2780 2578 0.75 0.74 1085 1106
7 2780 2577 0.75 0.74 1087 1111
8 2797 2531 0.76 0.74 1087 1113
9 2811 2528 0.76 0.73 1089 1118
10 2849 2519 0.76 0.73 1091 1122

Les colonnes \(\mathsf{EQM}\) et \(\widehat{\mathsf{EQM}}_a\) ont déjà été expliquées à la section précédente et ont été représentées graphiquement. On voit que l’erreur moyenne quadratique des données d’apprentissage, \(\widehat{\mathsf{EQM}}_a\), diminue toujours à mesure qu’on ajoute des variables (c’est-à-dire, qu’on augmente l’ordre du polynôme). Les critères \(\mathsf{AIC}\) et \(\mathsf{BIC}\) suggèrent le modèle cubique (\(k=3\)), c’est-à-dire le vrai modèle. Le \(R^2\) ajusté quant à lui choisit le modèle d’ordre \(4\) (qui est le deuxième meilleur selon le \(\mathsf{EQM}\)). N’oubliez pas que ces trois critères sont calculés avec l’échantillon d’apprentissage (\(n=100\)), mais en pénalisant l’ajout de variables. On est ainsi en mesure de contrecarrer le problème provenant du fait qu’on ne peut pas utiliser directement le \(\widehat{\mathsf{EQM}}_a\).

Le \(\mathsf{AIC}\) et le \(\mathsf{BIC}\) sont des critères très utilisés et très généraux. Ils sont disponibles dès qu’on utilise la méthode du maximum de vraisemblance est utilisée comme méthode d’estimation. Le \(R^2\) ajusté a une portée plus limitée car il est spécialisé à la régression linéaire.

4.7 Division de l’échantillon et validation croisée

La deuxième grande approche après celle consistant à pénaliser le \(\widehat{\mathsf{EQM}}_a\) consiste à tenter d’estimer le \(\mathsf{EQM}\) directement sans utiliser deux fois les mêmes données. Nous allons voir deux telles méthodes ici, la validation externe (division de l’échantillon) et la validation croisée (cross-validation).

Ces deux méthodes s’attaquent directement au problème qu’on ne peut utiliser (sans ajustement) les mêmes données qui ont servi à estimer les paramètres d’un modèle pour estimer sa performance. Pour ce faire, l’échantillon de départ est divisé en deux, ou plusieurs parties, qui vont jouer des rôles différents.

4.7.1 Validation externe et division de l’échantillon

Cette idée est très simple. Nous avons un échantillon de taille \(n\). Nous pouvons le diviser au hasard en deux parties de tailles respectives \(n_1\) et \(n_2\) (\(n_1+n_2=n\)),

  • un échantillon d’apprentissage (training) de taille \(n_1\) et
  • un échantillon de validation (test) de taille \(n_2\).

L’échantillon d’apprentissage servira à estimer les paramètres du modèle. L’échantillon de validation servira à estimer la performance prédictive (par exemple estimer l’\(\mathsf{EQM}\)) du modèle. Comme cet échantillon n’a pas servi à estimer le modèle lui-même, il est formé de « nouvelles » observations qui permettent d’évaluer d’une manière réaliste la performance du modèle. Comme il s’agit de nouvelles observations, on n’a pas à pénaliser la complexité du modèle et on peut directement utiliser le critère de performance choisi, par exemple, l’erreur quadratique moyenne, c’est-à-dire, la moyenne des erreurs au carré pour l’échantillon de validation. Cette quantité est une estimation valable de l’\(\mathsf{EQM}\) de ce modèle. On peut faire la même chose pour tous les modèles en compétition et choisir celui qui a la meilleure performance sur l’échantillon de validation.

Cette approche possède plusieurs avantages. Elle est facile à implanter. Elle est encore plus générale que les critères \(\mathsf{AIC}\) et \(\mathsf{BIC}\). En effet, ces critères découlent de la méthode d’estimation du maximum de vraisemblance. Plusieurs autres types de modèles ne sont pas estimés par la méthode du maximum de vraisemblance (par exemple, les arbres, les forêts aléatoires, les réseaux de neurones, etc.) La performance de ces modèles peut toujours être estimée en divisant l’échantillon. Cette méthode peut donc servir à comparer des modèles de familles différentes. Par exemple, choisit-on un modèle de régression linéaire, une forêt aléatoire ou bien un réseau de neurones?

Cette approche possède tout de même un désavantage. Elle nécessite une grande taille d’échantillon au départ. En effet, comme on divise l’échantillon, on doit en avoir assez pour bien estimer les paramètres du modèle (l’échantillon d’apprentissage) et assez pour bien estimer sa performance (l’échantillon de validation).

La méthode consistant à diviser l’échantillon en deux (apprentissage et validation) afin de sélectionner un modèle est valide. Par contre, si on veut une estimation sans biais de la performance du modèle choisi (celui qui est le meilleur sur l’échantillon de validation), on ne peut pas utiliser directement la valeur observée de l’erreur de ce modèle sur l’échantillon de validation car elle risque de sous-évaluer l’erreur. En effet, supposons qu’on a 10 échantillons et qu’on ajuste 10 fois le même modèle séparément sur les 10 échantillons. Nous aurons alors 10 estimations différentes de l’erreur du modèle. Il est alors évident que de choisir la plus petite d’entre elles sous-estimerait la vraie erreur du modèle. C’est un peu ce qui se passe lorsqu’on choisit le modèle qui minimise l’erreur sur l’échantillon de validation. Le modèle lui-même est un bon choix, mais l’estimation de son erreur risque d’être sous-évaluée.

Une manière d’avoir une estimation de l’erreur du modèle retenu consiste à diviser l’échantillon de départ en trois (plutôt que deux). Aux échantillons d’apprentissage et de validation, s’ajoute un échantillon « test ». Cet échantillon est laissé de côté durant tout le processus de sélection du modèle qui est effectué avec les deux premiers échantillons tel qu’expliqué plus haut. Une fois un modèle retenu (par exemple celui qui minimise l’erreur sur l’échantillon de validation), on peut alors évaluer sa performance sur l’échantillon test qui n’a pas encore été utilisé jusque là. L’estimation de l’erreur du modèle retenu sera ainsi valide. Il est évident que pour procéder ainsi, on doit avoir une très grande taille d’échantillon au départ.

4.7.2 Validation croisée

Si la taille d’échantillon n’est pas suffisante pour diviser l’échantillon en deux et procéder comme nous venons de l’expliquer, la validation croisée est une bonne alternative. Cette méthode permet d’imiter le processus de division de l’échantillon.

Voici les étapes à suivre pour faire une validation croisée à \(K\) groupes (\(K\)-fold cross-validation) :

  1. Diviser l’échantillon au hasard en \(K\) parties \(P_1, P_2, \ldots, P_K\) contenant toutes à peu près le même nombre d’observations.
  2. Pour \(j = 1\) à \(K\),
    1. Enlever la partie \(j\).
    2. Estimer les paramètres du modèle en utilisant les observations des \(K-1\) autres parties combinées.
    3. Calculer la mesure de performance (par exemple la somme du carré des erreurs) de ce modèle pour le groupe \(P_j\).
  3. Faire la somme des \(K\) estimations de performance pour obtenir une mesure de performance finale et repondérer au besoin.
Illustration de la validation croisée: on scinde l'échantillon d'apprentissage en cinq groupes (abcisse) et à chaque étape, une portion différente des données est mise de côté et ne sert que pour la validation.

Figure 4.1: Illustration de la validation croisée: on scinde l’échantillon d’apprentissage en cinq groupes (abcisse) et à chaque étape, une portion différente des données est mise de côté et ne sert que pour la validation.

On recommande habituellement de prendre entre \(K=5\) et \(10\) groupes (le choix de 10 groupes est celui qui revient le plus souvent en pratique). Si on prend \(K=10\) groupes, alors chaque modèle est estimé avec 90% des données et on prédit ensuite le 10% restant. Comme on passe en boucle les 10 parties, chaque observation est prédite une et une seule fois à la fin. Il est important de souligner que les groupes sont formés de façon aléatoire et donc que l’estimé que l’on obtient peut être très variable, surtout si la taille de l’échantillon d’apprentissage est petite. Il arrive également que le modèle ajusté sur un groupe ne puisse pas être utilisé pour prédire les observations mises de côté, notamment si des variables catégorielles sont présentes mais qu’une modalité n’est présente que dans un des groupes; ce problème se présente en pratique si certaines classes ont peu d’observations. Un échantillonnage stratifié permet de pallier à cette lacune et de s’assurer d’une répartition plus homogène des variables catégorielles.

Le cas particulier \(K=n\) (en anglais leave-one-out cross validation, ou \(\mathsf{LOOCV}\)) consiste à enlever une seule observation, à estimer le modèle avec les \(n-1\) autres et à valider à l’aide de l’observation laissée de côté: on répète cette procédure pour chaque observation. Pour les modèles linéaires, il existe des formules explicites qui nous permettent d’éviter d’ajuster \(n\) régressions par moindre carrés. Cette forme de validation croisée tend à être trop optimiste.

Revenons à notre exemple où une seule variable explicative est disponible et où l’on cherche à déterminer un bon modèle polynomial. Le Tableau 4.2 est le même que le Tableau 4.1 mais avec une colonne en plus, la dernière, \(\mathsf{VC} (K=10)\). Il s’agit des estimations de l’\(\mathsf{EQM}\) obtenues avec la validation croisée à 10 groupes. Notez que si vous exécutez le programme, vous n’obtiendrez pas les mêmes valeurs car il y a un élément aléatoire dans ce processus. La colonne représente la moyenne de 100 réplications.

Tableau 4.2: Mesures d’adéquation du modèle linéaire et estimés de l’erreur, incluant la validation croisée.
\(\mathsf{EQM}\) \(\widehat{\mathsf{EQM}}_a\) \(R^2\) \(R^2_a\) \(\mathsf{AIC}\) \(\mathsf{BIC}\) \(\mathsf{VC} (K=10)\)
1 3191 3674 0.65 0.65 1111 1119 3675
2 3133 2879 0.73 0.72 1088 1099 2898
3 2697 2620 0.75 0.74 1081 1094 2676
4 2767 2582 0.75 0.74 1081 1097 2666
5 2771 2581 0.75 0.74 1083 1102 2711
6 2780 2578 0.75 0.74 1085 1106 2757
7 2780 2577 0.75 0.74 1087 1111 2788
8 2797 2531 0.76 0.74 1087 1113 2846
9 2811 2528 0.76 0.73 1089 1118 2896
10 2849 2519 0.76 0.73 1091 1122 2976

Le modèle cubique (ordre 3) est aussi choisi par la validation croisée, en moyenne (comme il l’était par le \(\mathsf{AIC}\) et le \(\mathsf{BIC}\)). Le graphe qui suit trace les valeurs de l’estimation par validation croisée (courbe de validation croisée) et aussi le \(\mathsf{EQM}\). On voit que l’estimation par validation croisée suit assez bien la forme du \(\mathsf{EQM}\) (qu’il est supposé estimer). Les boîtes à moustache permettent d’apprécier la variabilité des estimés de l’erreur quadratique moyenne telles qu’estimée par validation croisée avec 10 groupes.

4.8 Cibler les clients pour l’envoi d’un catalogue

Nous allons présenter un exemple classique de commercialisation de bases de données qui nous servira à illustrer la sélection de modèles, la régression logistique et la gestion de données manquantes.

Le contexte est le suivant : une entreprise possède une grande base de données client. Elle désire envoyer un catalogue à ses clients mais souhaite maximiser les revenus d’une telle initiative. Il est évidemment possible d’envoyer le catalogue à tous les clients mais ce n’est possiblement pas optimal. La stratégie envisagée est la suivante :

  1. Envoyer le catalogue à un échantillon de clients et attendre les réponses. Le coût de l’envoi d’un catalogue est de 10$.
  2. Construire un modèle avec cet échantillon afin de décider à quels clients (parmi les autres) le catalogue devrait être envoyé, afin de maximiser les revenus.

Plus précisément, on s’intéresse aux clients de 18 ans et plus qui ont au moins un an d’historique avec l’entreprise et qui ont effectué au moins un achat au cours de la dernière année. Dans un premier lieu, on a envoyé le catalogue à un échantillon de 1000 clients. Un modèle sera construit avec ces 1000 clients afin de cibler lesquels des clients restants seront choisis pour recevoir le catalogue.

Pour les 1000 clients de l’échantillon d’apprentissage, les deux variables cibles suivantes sont disponibles :

  • yachat, une variable binaire qui indique si le client a acheté quelque chose dans le catalogue égale à 1 si oui et 0 sinon.
  • ymontant, le montant de l’achat si le client a acheté quelque chose.

Les 10 variables suivantes sont disponibles pour tous les clients et serviront de variables explicatives pour les deux variables cibles. Il s’agit de :

  • x1: sexe de l’individu, soit homme (0) ou femme (1);
  • x2: l’âge (en année);
  • x3: variable catégorielle indiquant le revenu, soit moins de 35 000$ (1), entre 35 000$ et 75 000$ (2) ou plus de 75 000$ (3);
  • x4: variable catégorielle indiquant la région où habite le client (de 1 à 5);
  • x5: conjoint : le client a-t-il un conjoint (0=non, 1=oui);
  • x6: nombre d’année depuis que le client est avec la compagnie;
  • x7: nombre de semaines depuis le dernier achat;
  • x8: montant (en dollars) du dernier achat;
  • x9: montant total (en dollars) dépensé depuis un an;
  • x10: nombre d’achats différents depuis un an.

Les données se trouvent dans le fichier DBM.sas7bdat. Voici d’abord des statistiques descriptives pour l’échantillon d’apprentissage.

Il y a 46,6% de femmes parmi les 1000 clients de l’échantillon. De plus, 39,7% ont un revenu de moins de 35 000$, 33,7% sont entre 35 000$ et 75 000$ et 26,6% ont plus de 75 000$. 42,5% de ces clients qui ont un conjoint. Le nombre d’achats différents depuis un an par ces clients varie entre 1 et 14. Un peu plus de la moitié (51,4%) ont fait cinq achats ou moins. Parmi les 1000 clients de l’échantillon d’apprentissage, 210 ont acheté quelque chose dans le catalogue. La variable yachat sera l’une des variables que nous allons chercher à modéliser en vue d’obtenir des prédictions.

L’âge des 1000 clients de l’échantillon d’apprentissage varie entre 20 et 70 avec une moyenne de 37,1 ans. En moyenne, ces clients ont acheté pour 229,30$ depuis un an. Le dernier achat de ces clients remonte, en moyenne, à 10 semaines. Nous chercherons également à modéliser la variable ymontant. Seuls 210 clients ont acheté quelque chose dans le catalogue et les statistiques rapportées correspondent seulement à ces derniers, car la variable ymontant est manquante si le client n’a rien acheté dans le catalogue. On pourrait également remplacer ces valeurs par des zéros et les modéliser, mais nous aborderons cet aspect ultérieurement. Les clients qui ont acheté quelque chose ont dépensé en moyenne 67,3$, et au minimum 25$. Les histogrammes de quelques unes de ces variables permet de mieux visualiser la répartition des observations.

Il y a plusieurs façons d’utiliser l’échantillon d’apprentissage afin de mieux cibler les clients à qui envoyer le catalogue et maximiser les revenus. En voici quelques unes.

  1. On pourrait développer un modèle afin d’estimer la probabilité qu’un client achète quelque chose si on lui envoie un catalogue. Plus précisément, on peut développer un modèle pour \(\Pr(\texttt{yachat}=1)\). Comme la variable yachat est binaire, un modèle possible est la régression logistique, que nous décrirons au chapitre suivant. Ainsi, en appliquant le modèle aux 100 000 clients restant, on pourra cibler les clients susceptibles d’acheter (ceux avec une probabilité élevée).
  2. Une autre façon serait de tenter de prévoir le montant d’argent dépensé. Nous venons de voir la distribution de la variable ymontant. Il y a deux situations, ceux qui ont acheté et ceux qui n’ont pas achetés. En conditionnant sur le fait d’avoir acheté quelque chose, il est possible de décomposer le problème de la manière suivante :

\[\begin{align*} {\mathsf E}\left(\texttt{ymontant}\right) &= {\mathsf E}\left(\texttt{ymontant} \mid \texttt{yachat}=1\right) {\mathsf P}\left(\texttt{yachat}=1\right) \\& \qquad \qquad \qquad + {\mathsf E}\left(\texttt{ymontant} \mid \texttt{yachat}=0\right) {\mathsf P}\left(\texttt{yachat}=0\right) \\ &= {\mathsf E}\left(\texttt{ymontant} \mid \texttt{yachat}=1\right) {\mathsf P}\left(\texttt{yachat}=1\right) \end{align*}\] En mots, la moyenne du montant dépensé est égale à la moyenne du montant dépensé étant donné qu’il y a eu achat, fois la probabilité qu’il ait eu achat.

On peut donc estimer \({\mathsf E}\left(\texttt{ymontant} \mid \texttt{yachat}=1\right)\) et \({\mathsf P}\left(\texttt{yachat}=1\right)\), pour ensuite les combiner et avoir une estimation de \({\mathsf E}\left(\texttt{ymontant}\right)\). Le développement du modèle pour \({\mathsf E}\left(\texttt{ymontant} \mid \texttt{yachat}=1\right)\) peut se faire avec la régression linéaire, en utilisant seulement les clients qui ont acheté dans l’échantillon d’apprentissage, car ymontant est une variable continue dans ce cas. Le développement du modèle pour \({\mathsf P}\left(\texttt{yachat}=1\right)\) peut se faire avec la régression logistique, tel que mentionné plus haut, en utilisant tous les 1000 clients de l’échantillon d’apprentissage. En fait, nous verrons plus loin qu’il est possible d’estimer conjointement les deux modèles avec un modèle Tobit. En appliquant le modèle aux 100 000 clients restants, on pourra cibler les clients qui risquent de dépenser un assez grand montant.

Comme nous n’avons pas encore vu la régression logistique, nous allons nous limiter à illustrer les méthodes qui restent à voir dans ce chapitre avec la régression linéaire en cherchant à développer un modèle pour \({\mathsf E}\left(\texttt{ymontant} \mid \texttt{yachat}=1\right)\), le montant d’argent dépensé par les clients qui ont acheté quelque chose.

La base de donnée contient deux variables explicatives catégorielles. Il s’agit de revenu (x3) et région (x4). Il faut coder d’une manière appropriée afin de pouvoir les incorporer dans les modèles. La manière habituelle est de créer des variables indicatrices (binaires) qui indiquent si la variable prend ou non une valeur particulière, dans SAS avec l’option class. En général, si une variable catégorielle possède \(K\) valeurs possibles, il est suffisant de créer \(K-1\) indicatrices, en laissant une modalité comme référence. Par exemple, pour x3, nous allons créer deux variables,

  • x31: variable binaire égale à 1 si x3 égale 1 et 0 sinon,
  • x32: variable binaire égale à 1 si x3 égale 2 et 0 sinon.

Ainsi, la valeur 3 est celle de référence. Ces deux indicatrices sont suffisantes pour récupérer toute l’information comme le démontre le tableau 4.3.

Tableau 4.3: Valeur des indicateurs en fonction du niveau de la variable catégorielle
x3 x31 x32
1 1 0
2 0 1
3 0 0

4.9 Recherche automatique du meilleur modèle

Lorsque nous voulons comparer un petit nombre de modèles, il est relativement aisé d’obtenir les critères (\(\mathsf{AIC}\), \(\mathsf{BIC}\) ou autre) pour tous les modèles et de choisir le meilleur. C’était le cas dans l’exemple du choix de l’ordre du polynôme où il y avait seulement 10 modèles en compétitions. Mais lorsqu’il y a plusieurs variables en jeu, le nombre de modèles potentiel augmente très rapidement.

En fait, supposons qu’on a \(p\) variables distinctes disponibles. Avant même de considérer les transformations des variables et les interactions entre elles, il y a déjà modèles possibles. En effet, chaque variable est soit incluse ou pas (deux possibilités) et donc il y a \(2^p=2\times 2 \times \cdots \times 2\) (\(p\) fois) modèles en tout à considérer. Ce nombre augmente très rapidement comme en témoigne le tableau 4.4.

Tableau 4.4: Nombres de modèles en fonction du nombre de paramètres \(p\).
\(p\) nombre de paramètres
5 32
10 1024
15 32768
20 1048576
25 33554432
30 1073741824

Ainsi, si le nombre de variables est restreint, il est possible de comparer tous les modèles potentiels et de choisir le meilleur (selon un critère). II existe même des algorithmes très efficaces qui permettent de trouver le meilleur modèle sans devoir examiner tous les modèles possibles. Le nombre de variables qu’il est possible d’avoir dépend de la puissance de calcul et augmente d’année en année. Par contre, dans plusieurs applications, il ne sera pas possible de comparer tous les modèles et il faudra effectuer une recherche limitée. Faire une recherche exhaustive parmi tous les modèles possibles s’appelle sélection de tous les sous-ensembles (best subsets). La procédure reg de SAS permet de faire cela pour la régression linéaire.

On veut trouver un bon modèle pour prévoir la valeur de ymontant des clients qui ont acheté quelque chose. On a vu qu’il y a 210 clients qui ont acheté dans l’échantillon d’apprentissage. Nous allons chercher à développer un « bon » modèle avec ces 210 clients. Dans ce premier exemple, nous allons seulement utiliser les 10 variables explicatives de base (14 variables avec les indicatrices). Le code suivant montre comment faire une sélection de variables selon le critère du \(R^2\) et demande à SAS de présenter le modèle à \(k\) variables (\(k=1, \ldots, 14\)) qui a le plus grand \(R^2\); voir selection2_methodes.sas pour plus de détails.

Pour un nombre de variables fixé, le meilleur modèle selon le \(R^2\) est aussi le meilleur selon les critères d’information \(\mathsf{AIC}\) et \(\mathsf{BIC}\), pour ce nombre fixé de variables. Pour vous convaincre de cette affirmation, fixons le nombre de variables et restreignons-nous seulement aux modèles avec ce nombre de variables. Comme = \(1 - \mathsf{SCE}/\mathsf{SCT}\) et que \(\mathsf{SCT}\) est une constante indépendante du modèle, le modèle avec le plus grand coefficient de détermination, \(R^2\), est aussi celui avec la plus petite somme du carré des erreurs (\(\mathsf{SCE}\)). Comme \(\mathsf{AIC}=n (\ln (\mathsf{SCE}/n)) + 2p\), ce sera aussi celui avec le plus petit \(\mathsf{AIC}\) car la pénalité \(2p\) est la même si on fixe le nombre de variables; la même remarque est valide pour le \(\mathsf{BIC}\).

Ainsi, pour trouver le meilleur modèle globalement (sans fixer le nombre de variables), il suffit de trouver le modèle à \(k\) variables explicatives ayant le coefficient de détermination le plus élevé pour tous les nombres de variables fixés et d’ensuite de trouver celui qui minimise le \(\mathsf{AIC}\) (ou le \(\mathsf{BIC}\)) parmi ces modèles. Cette astuce est utile dans la mesure où SAS ne permet pas de faire cette même recherche avec les critères d’information. Ainsi, le modèle linéaire simple qui a le plus grand \(R^2\) est celui qui inclut le conjoint (x5). Le meilleur modèle (selon le \(R^2\)) parmi tous les modèles avec deux variables est celui avec x5 et x6.

Dans SAS, seule la procédure reg permet de faire cette recherche exhaustive et nous guaranti de recouvrer le modèle avec le plus grand \(R^2\). Un algorithme par séparation et évaluation permet d’effectuer cette recherche de manière efficace sans essayer tous les candidats pour ces sous-ensembles. Dans l’exemple, on voit que le modèle avec les variables x1 x2 x31 x44 x5 x6 x7 x8 x9 et x10 est celui qui minimise le \(\mathsf{AIC}\) globalement (\(\mathsf{AIC}=660.15\)). Le modèle choisi par le \(\mathsf{BIC}\) contient seulement sept variables explicatives (plutôt que 10), soit x1 x31 x5 x6 x7 x8 x10.

Nous avons seulement inclus les variables de base pour ce premier essai. Il est possible qu’ajouter des variables supplémentaires améliore la performance du modèle. Pour cet exemple, nous allons considérer les variables suivantes:

  • les variables continues au carré, comme \(\texttt{age}^2\).
  • toutes les interactions d’ordre deux entre les variables de base, comme \(\texttt{sexe}\cdot\texttt{age}\).

Si on suppose que la vrai moyenne de la variable réponse ymontant est lisse, le modèle ajusté précédemment capture l’approximation de degré 1 (série de Taylor) du vrai modèle et le modèle avec termes quadratiques (incluant les interactions) l’approximation de degré 2. Aux variables de base (10 variables explicatives, mais 14 avec les indicatrices pour les variables catégorielles), s’ajoutent ainsi 90 autres variables. Il y a donc 104 variables explicatives potentielles si on inclut les interactions et les termes quadratiques. Notez qu’il y a des interactions entre chacune des variables indicatrices et chacune des autres variables, mais il ne sert à rien de calculer une interaction entre deux indicatrices d’une même variable (car une telle variable est zéro pour tous les individus). De même, il ne sert à rien de calculer le carré d’une variable binaire.

Dans la mesure ou on aura presque un ratio d’un paramètre pour deux observations, ajuster un modèle avec toutes les variables explicatives est profondément stupide ici. Pensez à la taille de votre échantillon comme à un budget et aux paramètres comme à un nombre d’items: plus vous achetez d’items, moins votre budget est élevé pour chacun et leur qualité en pâtira. Le modèle à 104 variables servira uniquement à illustrer le surajustement. Réalistement, un modèle avec plus d’une vingtaine de variables ici serait difficilement estimable de manière fiable et l’inclusion d’interactions et de termes quadratiques sert surtout à augmenter la flexibilité et les possibilités lors de la sélection de variables.

Lancer une sélection exhaustive de tous les sous-modèles avec 104 variables risque de prendre un temps énorme. Que faire alors? Il y a plusieurs possibilités. Nous pourrions faire une recherche limitée avec les méthodes que nous allons voir à partir de la section suivante. Nous pourrions aussi combiner les deux approches. Supposons que notre ordinateur permet de faire une recherche exhaustive de tous les sous-modèles avec 40 variables. Nous pourrions alors commencer avec une recherche limitée pour trouver un sous-ensemble de 40 « bonnes » variables et faire une recherche exhaustive, mais en se restraignant à ces 40 variables.

4.10 Méthodes classiques de sélection

Les méthodes de sélection ascendante, descendante et séquentielle sont des algorithmes gloutons qui permettent de choisir des variables. Elles ont été développées à une époque où la puissance de calcul était bien moindre, et où il était impossible de faire une recherche exhaustive des sous-modèles. Avec l’approche classique, ces méthodes font une recherche séquentielle guidée parmi un nombre limité de modèles, à l’aide des valeurs-p du test-t pour la significativité des paramètres individuels du modèle avec \(p\) prédicteurs potentiels \(X_1,\ldots, X_p\). Les procédures glmselect et reg permettent une sélection de modèle avec une approche séquentielle, ascendante ou descendante.

4.10.1 Sélection ascendante

L’idée de la sélection ascendante est de tester l’ajout de chaque variable individuellement et d’ajouter celle qui est la plus significative selon le test-t si elle a une valeur-p assez petite.

  • Initialisation: le modèle linéaire de départ est celui qui n’inclut que l’ordonnée à l’origine, \(Y=\beta_0+\varepsilon\), où \(\varepsilon\) est une erreur centrée.
  • Critère d’entrée: \(c\), valeur-p minimale à partir de laquelle une variable peut être incluse dans le modèle
  • Boucle soit \(X_{(1)}, \ldots, X_{(k)}\), les variables explicatives à l’étape \(k<p\).
    • pour chaque \(j\) (\(j=\{1,\ldots, p\}\setminus \{(1), \ldots (k)\}\)), on ajuste tour à tour le modèle \(Y=\beta_0+\sum_{i=1}^k \beta_i X_{(i)} + \beta_{k+1}X_{j}\) et on calcule la valeur-p du test-t pour les hypothèses \(\mathcal{H}_0: \beta_{k+1}=0\) contre l’alternative bilatérale \(\mathcal{H}_1: \beta_{k+1} \neq 0\).
    • Soit \(p_{\min}\) la plus petite des \(p-k\) valeurs-p qui correspond à \(X_{(k+1)}\), disons.
      • si \(p_{\min}<c\), continuer la procédure.
      • si \(p_{\min} \geq c\), retourner le modèle \(Y=\beta_0 + \sum_{i=1}^k \beta_iX_{(i)}+\varepsilon\).

On continue ainsi à ajouter des variables jusqu’à ce que le critère d’entrée ne soit pas satisfait. Si on se rend jusqu’au bout, on va terminer avec le modèle complet qui contient toutes les variables.

4.10.2 Sélection descendante

  • Initialisation: le modèle linéaire de départ est celui qui inclut toutes les variables explicatives, \(Y=\beta_0+\sum_{j=1}^p \beta_j X_{(j)}+\varepsilon\), où \(\varepsilon\) est une erreur centrée.
  • Critère de sortie: \(c\), valeur-p maximale à partir de laquelle une variable peut être excluse du modèle
  • Boucle soit \(X_{(1)}, \ldots, X_{(p-k)}\), les variables explicatives présentes dans le modèle à l’étape \(k<p\).
    • pour chaque \(j\) (\(j =1, \ldots, p-k\)), on calcule la valeur-p du test-t \(\mathcal{H}_0: \beta_{j}=0\) contre l’alternative bilatérale \(\mathcal{H}_1: \beta_{j} \neq 0\).
    • si toutes ces valeurs sont inférieures à \(c\), on retourne le modèle \(Y=\beta_0 + \sum_{j=1}^{p-k} \beta_j X_{(j)}\).
    • sinon, on enlève la variable qui a la plus grande valeur-p (disons \(X_{(p-k)}\)), on réajuste le modèle sans cette variable et on recommence la procédure.

L’idée est l’inverse de la méthode ascendante. On va tester le retrait de chaque variable individuellement et retirer celle qui est la moins significative, si sa valeur-p est assez grande. Si la procédure se termine après \(p\) itérations, aucune variable explicative n’est retenue.

4.10.3 Méthode séquentielle

Il s’agit d’une méthode hybride entre ascendante et descendante. On sélectionne un critère d’entrée et de sortie pour chacune des deux et on début la recherche à partir du modèle ne contenant que l’ordonnée à l’origine. À chaque étape, on fait une étape ascendante suivie de une (ou plusieurs) étapes descendantes. On continue ainsi tant que le modèle retourné par l’algorithme n’est pas identique à celui de l’étape précédente. Le dernier modèle est celui retenu.

Avec la méthode séquentielle, une fois qu’on entre une variable (étape ascendante), on fait autant d’étapes descendante afin de retirer toutes les variables qui satisfont le critère de sortie (il peut ne pas y en avoir). Une fois cela effectué, on refait une étape ascendante pour voir si on peut ajouter une nouvelle variable.

Remarques sur ces méthodes: avec la méthode ascendante, une fois qu’une variable est dans le modèle, elle y reste. Avec la méthode descendante, une fois qu’une variable est sortie du modèle, elle ne peut plus y entrer. Avec la méthode séquentielle, une variable peut entrer dans le modèle et sortir plus tard dans le processus. Par conséquent, parmi les trois, la méthode séquentielle est généralement préférable aux méthodes ascendante et descendante, car elle inspecte potentiellement un plus grand nombre de modèles.

On peut soi-même spécifier les critères d’entrée et de sortie. Plus le critère d’entrée est élevé, plus il y aura de variables dans le modèle final. De même, plus le critère de sortie est élevé, plus il y aura de variables dans le modèle. Il y également moyen avec glmselect de spécifier le nombre d’étapes de chaque procédure.

Utilisons la méthode de sélection séquentielle classique avec des critères d’entrée et de sortie de 0.15 et les 104 variables. La sortie SAS est assez volumineuse car elle retrace toutes étapes de la sélection séquentielle. L’historique montre qu’à l’étape 1, la variable d’interaction x5_x6 a été ajoutée, suivie de x31_x10 à l’étape 2. Un peu plus loin, à l’étape 6, x5_x6 est retirée et ainsi de suite. Il y a eu 40 étapes en tout et, à la fin, il reste 22 variables (parmi les 104) dans le modèle final. Le \(R^2\) du modèle final est 0,966.

On voit bien que toutes les valeurs-p (qui ne sont pas valides à cause de la sélection de modèles) sont toutes inférieures à 0.15. Le choix de 0.15 comme critère d’entrée et de sortie est complètement arbitraire. Il est fort possible que d’autres valeurs donnent de meilleurs résultats, mais il n’est pas évident de les choisir.

Une façon de contourner le problème de devoir spécifier les critères d’entrée et de sortie est de procéder en deux étapes. Supposons que notre ordinateur permet de faire une recherche exhaustive de tous les sous-modèles avec près de 60 variables. L’idée est alors de passer de 104 à un sous-ensemble d’environ 60 variables, avec une sélection séquentielle gloutonne, et d’ensuite utiliser une recherche exhaustive avec ce sous-ensemble de variables. Plus précisément:

  1. On fait une sélection séquentielle classique avec des valeurs élevées pour les critères d’entrée et de sortie afin que le modèle retenu contienne le nombre voulu de variables (par exemple, 60).
  2. En utilisant seulement ce sous-ensemble de variables, on choisit le meilleur modèle selon le \(\mathsf{AIC}\) ou le \(\mathsf{BIC}\) en faisant une recherche exhaustive de tous les sous-modèles.

En fixant, les critères d’entrée et de sortie à 0,6 pour la recherche séquentielle, le modèle retenu aura 56 variables. Il est possible de faire une recherche exhaustive avec 56 variables sur un ordinateur portable avec SAS. Le \(\mathsf{AIC}\) est mène à un modèle avec 38 de ces 56 variables, ce qui est probablement trop. Le \(\mathsf{BIC}\) est quant à lui beaucoup plus parcimonieux et choisit 15 de ces variables pour le modèle final. Nous verrons à la section suivante qu’il est possible de faire une recherche séquentielle en utilisant d’autres critères que la valeur-p du test-t pour faire ajouter ou enlever des variables.

4.11 Recherche séquentielle automatique limitée

L’idée de la procédure séquentielle classique est d’inclure ou d’exclure une variable à la fois sur la base des valeurs-p. La procédure glmselect permet de faire une sélection séquentielle en utilisant d’autres critères, comme le \(\mathsf{AIC}\) ou le \(\mathsf{BIC}\). Cette procédure permet de contrôler très finement le processus de sélection de variables. Le code qui suit fait une recherche séquentielle avec les particularités suivantes. À chaque étape ascendante de la procédure séquentielle, c’est la variable qui améliore le plus le \(\mathsf{AIC}\) (select=aic) qui est entrée. De plus, à chaque étape descendante de la procédure séquenctielle, c’est la (ou les) variable(s) qui détériore(nt) le plus le \(\mathsf{AIC}\) qui est (sont) retirée(s). À la toute fin du processus, c’est le modèle qui a le meilleur \(\mathsf{BIC}\) (choose=BIC) qui est retenu.

La procédure glmselect permet la déclarations de variables catégorielles (class); on évite ainsi de créer les variables binaires une par une. Voici les commandes pour faire une procédure séquentielle avec des critères plus généreux (entrée=sortie=0.6). Si le nombre de variables change en sortie, nous aurons ici 56 variables qui seront ensuite utilisées avec une recherche exhaustive: pour ce faire, on reprend la sortie, mais cette fois on choisit le modèle par la suite qui a le plus petit BIC (SBC) ou AIC.

 proc glmselect data=ymontant outdesign=glmselectoutput;
 partition role=train(train="1" validate="0");
 class x3(param=ref split) x4(param=ref split);
 model ymontant=x1|x2|x3|x4|x5|x6|x7|x8|x9|x10 @2
 x2*x2 x6*x6 x7*x7 x8*x8 x9*x9 x10*x10 /  
 selection=stepwise(slentry=0.6 slstay=0.6 select=60) hier=none;
 run;

 proc glmselect data=glmselectoutput;
 model ymontant= &_GLSMOD / 
 selection=backward(stop=1 choose=sbc) hier=none;
 ods output GLMSelect.Summary.SelectionSummary=historiqueSelection;
 run;

À l’étape 1, la variable x5_x6 est ajoutée au modèle de base car c’est celle qui fait diminuer le plus le \(\mathsf{AIC}\). À l’étape 2, la variable x31_x10 est ajoutée, À l’étape 6, la variable x5_x6 est retirée car cela fait baisser le \(\mathsf{AIC}\). Notez que le \(\mathsf{AIC}\) décroit toujours d’une étape à l’autre. SAS garde aussi la trace du \(\mathsf{BIC}\) car le modèle final sera choisi selon ce critère. Finalement le processus séquentiel se termine à l’étape 40, car il n’y a plus moyen de faire diminuer le \(\mathsf{AIC}\). Le modèle final retenu est celui de l’étape 18, car c’est celui qui a le \(\mathsf{BIC}\) le plus petit parmi tous ces modèles (\(\mathsf{BIC}=484.22\)).

Voici différentes statistiques ainsi que les estimations des paramètres de ce modèle qui contient 10 variables.

On peut voir sur le graphique @ref(fig:fig2p2_fig17) l’historique des valeurs de AIC et BIC à mesure qu’on diminue le nombre de variables dans le modèle: les mêmes variables sont enlevées à chaque étape, mais la valeur optimale du critère est différente pour la sélection finale. Sur l’axe des abscisses, j’ai ajouté l’erreur quadratique moyenne pour l’échantillon de validation contenant les 23 389 clients (sur 100 000) qui achèteraient si leur envoyait une copie du catalogue. Cet exemple n’est pas réaliste puisqu’on regarde la solution, mais il permet de nous comparer et de voir à quel point ici le critère d’information bayésien suit la même tendance que l’erreur moyenne quadratique de validation.

4.12 Méthodes de régression avec régularisation

Une façon d’éviter le surajustement est d’ajouter une pénalité sur les coefficients: ce faisant, on introduit un biais dans nos estimés, mais dans l’espoir de réduire leur variabilité et ainsi d’obtenir une meilleur erreur quadratique moyenne.

Les estimateurs des moindres carrés ordinaires pour la régression linéaire représentent la combinaison qui minimise la somme du carré des erreurs, \[\begin{align*} \mathsf{SCE} = \sum_{i=1}^n \left(y_i - \beta_0 - \sum_{j=1}^p\mathrm{X}_{ij}\beta_{j}\right)^2. \end{align*}\] On peut ajouter à cette fonction objective \(\mathsf{SCE}\) un terme additionnel de pénalité qui va contraindre les paramètres à ne pas être trop grand. On considère une pénalité additionnelle pour la valeur absolue des coefficients, \[ q_1(\lambda) = \lambda \sum_{j=1}^p |\beta_j|. \] Pour chaque valeur de \(\lambda\) donnée, on obtiendra une solution différente pour les estimés car on minimisera désormais \(\mathsf{SCE} + q_1(\lambda)\). On ne pénalise pas l’ordonnée à l’origine \(\beta_0\), parce que ce coefficient sert à recentrer les observations et a une signification particulière. Si on standardise les données, de manière à ce que leur moyenne empirique soit zéro et leur écart-type un, alors \(\widehat{\beta}_0 = \overline{y}\).

L’avantage des moindres carrés est que les valeurs ajustées et les prédictions ne changent pas si on fait une transformation affine (de type \(Z = aX+b\)). Peu importe le choix d’unités (par exemple, exprimer une distance en centimètres plutôt qu’en mètres, ou la température en Farenheit plutôt qu’en Celcius), on obtient le même modèle. En revanche, une fois qu’on introduit un terme de pénalité, notre solution dépendra de l’unité de mesure, d’où l’importance d’utiliser les données centrées et réduites pour que la solution reste la même.

La pénalité \(q_1(\lambda)\) a un rôle particulier parce qu’elle a deux effets: elle réduit la taille des paramètres, mais elle force également certains paramètres très proches de zéro à être exactement égaux à zéro, ce qui fait que la régression pénalité agit également comme outil de sélection de variables. Des algorithmes efficaces permettent de trouver la solution du problème d’optimisation \[\begin{align*} \min_{\boldsymbol{\beta}} \{\mathsf{SCE} + q_1(\lambda)\} = \min_{\boldsymbol{\beta}} \left\{\sum_{i=1}^n \left(y_i - \beta_0 - \sum_{j=1}^p\mathrm{X}_{ij}\beta_{j}\right)^2 + \lambda \sum_{j=1}^p |\beta_j|\right\} \end{align*}\] laquelle est appelée LASSO. La Figure 4.2 montre la fonction objective dans le cas où on a deux paramètres, \(\beta_1\) et \(\beta_2\). La solution des moindres carrés ordinaires, qui minimisent l’erreur quadratique moyenne, est au centre des ellipses de contour et correspond à la solution du modèle avec \(\lambda=0\). À mesure que l’on augmente la pénalité \(\lambda\), les coefficients rétrécissent vers \((0, 0)\). On peut interpréter la pénalité \(l_1\) comme une contraire budgétaire: les coefficients estimés pour une valeur de \(\lambda\) donnée sont ceux qui minimisent la somme du carré des erreurs, mais doivent être à l’intérieur d’un budget alloué (losange). La forme de la région fait en sorte que la solution, qui se trouve sur la bordure du losange, intervient dans un coin avec certaines coordonnées nulles.

Courbes de contour du critère de l'erreur quadratique moyenne (ellipses) et fonction de pénalité (losanges) pour différentes valeurs de $\lambda$. Les points dénotent des solutions différentes et intersectent les contours du losange.

Figure 4.2: Courbes de contour du critère de l’erreur quadratique moyenne (ellipses) et fonction de pénalité (losanges) pour différentes valeurs de \(\lambda\). Les points dénotent des solutions différentes et intersectent les contours du losange.

Plusieurs variantes existent dans la littérature qui généralisent le modèle à des contextes plus compliqués. Le choix des variables à inclure dans la sélection dépend du choix de la pénalité \(\lambda\), qui est règle générale estimée par validation croisée à 5 ou 10 groupes.

4.13 Moyenne de modèles

Une idée importante et moderne en statistique est qu’il est souvent préférable de combiner plusieurs modèles plutôt que d’en choisir un seul. La technique des forêts aléatoires (random forests) est une des meilleures techniques de prédiction disponibles de nos jours. Elle est basée sur cette idée, en combinant plusieurs arbres de classification (ou de régression) individuels. C’est une des techniques de base en exploitation de données.

Ici, nous allons voir comment cette idée peut être appliquée à notre contexte. Toutes les méthodes que nous avons vues jusqu’à maintenant font une sélection « rigide » de variables, dans le sens que chaque variable est soit sélectionnée pour faire partie du modèle, soit elle ne l’est pas. C’est donc tout ou rien pour chaque variable. Il y a beaucoup de variabilité associée à une telle forme de sélection. Une variable peut avoir été très près d’être choisie, mais elle ne l’a pas été et est éliminée complètement. Construire plusieurs modèles et en faire la moyenne permet d’adoucir le processus de sélection car une variable peut alors être partiellement sélectionnée.

Supposons qu’on dispose de deux échantillons et qu’on fasse une sélection de variables séparément pour les deux échantillons, avec l’une des approches que nous avons vues jusqu’à maintenant. Il est alors très probable qu’on ne va pas avoir exactement les mêmes variables sélectionnées pour les deux échantillons. Supposons ensuite qu’on fasse la moyenne des coefficients pour les deux modèles. Si une variable, disons \(X_1\), a été choisie les deux fois, alors la moyenne des deux coefficients devrait estimer en quelque sorte un effet global pour cette variable. Si une autre variable, disons \(X_2\), n’a pas été choisie du tout pour les deux échantillons, alors la moyenne de ses deux coefficients est nulle. Mais si une variable, disons, \(X_3\), a été choisie pour seulement l’un des deux échantillons, alors la moyenne de ses deux coefficients est la moitié du coefficient pour le modèle dans lequel elle a été choisie (car l’autre est zéro). Ainsi, cette variable est donc représentée par une « moitié » d’effet dans la moyenne des modèles. Donc au lieu d’être totalement là ou totalement absente, elle est présente en fonction de sa probabilité d’être sélectionnée. Ceci diminue de beaucoup la variabilité engendrée par une sélection « rigide » de variables et permet souvent de produire un modèle fort raisonnable.

Le problème est que l’on n’a pas plusieurs échantillons mais un seul. Une solution possible est de générer nous-mêmes des échantillons différents à partir de l‘échantillon original. Cela peut être fait avec l’autoamorçage (bootstrap). Un échantillon d’autoamorçage est tout simplement un échantillon choisi au hasard et avec remise dans l’échantillon original. Ainsi, une même observation peut être sélectionnée plus d’une fois tandis qu’une autre peut ne pas être sélectionnée du tout.

L’idée est alors la suivante :

  1. Générer plusieurs échantillons par autoamorçage nonparamétrique à partir de l‘échantillon original.
  2. Faire une sélection de variables pour chaque échantillon.
  3. Faire la moyenne des paramètres de ces modèles.

La procédure glmselect a une commande, modelaverage, qui permet de faire une moyenne de modèles. Le code suivant permet de faire une moyenne de modèles.

proc glmselect data=ymontant seed=57484765;
partition role=train(train="1" validate="0");
class x3(param=ref split) x4(param=ref split);
model ymontant=x1|x2|x3|x4|x5|x6|x7|x8|x9|x10 @2
x2*x2 x6*x6 x7*x7 x8*x8 x9*x9 x10*x10 / 
selection=stepwise(select=sbc choose=sbc) hier=none; 
score data=testymontant out=predaverage p=predymontant;
modelaverage nsamples=500 sampling=urs subset(best=500);
run;

Chaque modèle est construit à l’aide d’un échantillon aléatoire avec remise (sampling=urs). Il y aura 500 échantillons, et donc modèles, en tout (nsamples=500). L’option subset(best=500) indique à SAS de faire la moyenne des paramètres des 500 modèles. Notez l’option seed qui permet de reproduire les résultats, car elle fixe une valeur pour le générateur de nombre aléatoire (qui sera utilisé pour générer les échantillons d’autoamorçage). Cette fois-ci la sélection se fait avec le critère \(\mathsf{BIC}\) à tous les niveaux (select=aic choose=sbc).

Ce tableau présente les variables qui ont été choisies dans au moins 20% des modèles, c’est-à-dire, dans au moins 100 des 500 modèles ici. Il y a deux variables qui ont été retenues dans tous les modèles, x1_x6 et x31_x8. Le tableau rapporte aussi la moyenne des estimations pour ces paramètres.

Toutes les méthodes employées jusqu’à maintenant utilise une méthode de pénalisation pour déterminer le meilleur modèle. Une alternative avec glmselect serait de répéter la sélection en utilisant directement l’erreur moyenne quadratique estimée à l’aide de la validation croisée comme critère de sélection (en remplaçant l’option choose=cv avec par exemple cvmethod=random(10) pour de la sélection avec \(K=10\) groupes formées aléatoirement).

4.14 Évaluation de la performance

La direction de la compagnie a décidé de passer outre vos recommandations et d’envoyer le catalogue aux 100 000 clients restants; nous pouvons donc faire un post-mortem afin de voir ce que chaque modèle aurait donné comme profit, comparativement à la stratégie de référence. Les 100 000 autres clients serviront d’échantillon de validation pour évaluer la performance des modèles et, plus précisément, afin d’évaluer les revenus (ou d’autres mesures de performance) si ces modèles avaient été utilisés. L’échantillon de validation nous donnera donc l’heure juste quant aux mérites des différentes approches que nous allons comparer. En pratique, nous ne pourrions pas faire cela car la valeur de la variable cible ne serait pas connue pour ces clients et nous utiliserions plutôt les modèles pour obtenir des prédictions pour déterminer quels clients cibler avec l’envoi. Parmi, les 100 000 clients restants, il y en a 23 179 qui auraient acheté quelque chose si on leur avait envoyé le catalogue. Ces 23 179 observations vont nous servir pour estimer l’erreur quadratique moyenne (théorique) des modèles retenus par nos critères (voir le fichier selection2_methodes.sas pour les manipulations). Comme notre base de données contient les 101 000 observations avec une étiquette train, on peut spécifier la partition avec partition role=train(train="1" validate="0") pour obtenir directement l’erreur moyenne quadratique pour l’échantillon de validation; cette stratégie fonctionne avec toutes les méthodes, sauf la moyenne de modèles pour laquelle on devra faire le calcul manuellement en obtenant les prédictions.

Commençons par l’estimation de l’erreur quadratique moyenne (moyenne des carrés des erreurs) pour les deux modèles retenus par le \(\mathsf{AIC}\) et le \(\mathsf{BIC}\) avec les variables de base. Le tableau 4.5 contient aussi l’estimation de l’erreur quadratique moyenne si on utilise toutes les variables (14 en incluant les indicatrices) sans faire de sélection. On voit que le modèle choisi par le \(\mathsf{BIC}\) est le meilleur des trois, car l’erreur quadratique moyenne sur l’échantillon test est de 3,6% inférieure à celle du modèle choisi par le \(\mathsf{AIC}\). Ces deux méthodes font mieux que le modèle qui inclut toutes les variables sans faire de sélection, mais nous verrons que leur performance est exécrable: les variables de base ne permettent pas de capturer les effets présents dans les données et ce manque de flexibilité coûte cher.

Tableau 4.5: Estimation de l’erreur quadratique moyenne sur l’échantillon test avec les variables de base. Les meilleurs modèles selon les critères d’informations découlent d’une recherche exhaustive de tous les sous-ensembles.
nombre de variables \(\mathsf{EQM}\) méthode
14 25,69 toutes les variables
10 24,72 exhaustive - \(\mathsf{AIC}\)
7 23,83 exhaustive - \(\mathsf{BIC}\)
Tableau 4.6: Comparaison des méthodes selon l’erreur quadratique moyenne avec les variables de base, les interactions et les termes quadratiques.
nombre de variables \(\mathsf{EQM}\) méthode
104 19,63 toutes les variables
22 12,25 séquentielle classique
38 14,83 séquentielle classique, recherche exhaustive avec 56 variables \((\mathsf{AIC})\)
15 11,96 séquentielle classique, recherche exhaustive avec 56 variables \((\mathsf{BIC})\)
10 10,08 séquentielle avec critère \(\mathsf{AIC}\) (choix selon le \(\mathsf{BIC}\))
26 11,61 LASSO, validation croisée avec 10 groupes
10,57 moyenne de modèles

Le tableau 4.6 présente la performance de toutes les méthodes avec les autres variables. On voit d’abord qu’utiliser toutes les 104 variables sans faire de sélection fait mieux (\(\mathsf{EQM}=19.63\)) que les modèles précédents basés sur les 10 variables originales. Mais faire une sélection séquentielle classique permet une amélioration très importante de la performance (\(\mathsf{EQM}=12.25\)). Utiliser les 104 variables mène à du surajustement (over-fitting).

La stratégie consistant à sélectionner un sous-ensemble de 56 variables avec la méthode séquentielle classique pour ensuite faire une recherche exhaustive de tous les sous-modèles possibles avec ces 56 variables, selon le \(\mathsf{BIC}\), donne le meilleur résultat jusqu’à présent (\(\mathsf{EQM}=11.96\)). Le \(\mathsf{AIC}\) fait moins bien dans ce cas, avec une erreur quadratique moyenne estimée de \(14.83\). Les méthodes séquentielles avec un critère d’information (qui pénalisent davantage que les tests d’hypothèse classique) mènent à des modèles plus parcimonieux et qui réduisent encore l’erreur moyenne quadratique. La moyenne de modèles est compétitive, tandis que le LASSO performe moins bien (probablement parce que les coefficients sont tous rétrécis vers zéro, ce qui engendre du biais).

Dans cet exemple, la méthode séquentielle de glmselect avec les options select=aic et choose=bic aurait donné le meilleur résultat pour prévoir le montant acheté des clients restants (de ceux qui auraient acheté quelque chose). Le deuxième meilleur aurait été la moyenne des modèles. Il faut bien comprendre qu’il ne s’agit que d’un seul exemple: il ne faut surtout pas conclure que la méthode séquentielle de glmselect avec les options select=aic et choose=bic sera toujours la meilleure. En fait, il est impossible de prévoir quelle méthode donnera les meilleurs résultats.

Il y aurait plusieurs autres approches/combinaisons qui pourraient être testées. Le but de ce chapitre était simplement de présenter les principes de base en sélection de modèles et de variables ainsi que certaines approches pratiques. Il y a d’autres approches intéressantes, tels le filet élastique ou la régression LARS (least-angle regression) qui sont disponibles dans glmselect. Ces méthodes sont dans la même mouvance moderne que celle qui consiste à faire la moyenne de plusieurs modèles, en performant à la fois une sélection de variables et en permettant d’avoir des parties d’effet par le rétrécissement (shrinkage). De récents développements théoriques permettent aussi de corriger les valeurs-p pour faire de l’inférence post-sélection avec le LASSO.