Dans le champ contemporain de l’apprentissage automatique et de la modélisation statistique avancée, l’abondance des données constitue simultanément une opportunité remarquable et un défi analytique considérable. L’accès banalisé à des vecteurs de caractéristiques de très grande dimension expose le statisticien et l’ingénieur de données à une pathologie bien connue : le fléau de la dimensionnalité. Lorsqu’un modèle intègre un nombre disproportionné de variables prédictives par rapport au nombre d’observations empiriques disponibles, les estimations perdent leur stabilité, la variance des paramètres explose et l’interprétabilité phénoménologique du système s’effondre. Face à cette prolifération de prédicteurs potentiels, la sélection de variables s’impose non pas comme une simple étape technique d’ingénierie préliminaire, mais comme un impératif épistémologique visant à isoler le signal fondamental du bruit aléatoire inhérent au monde physique et humain.
Parmi l’arsenal méthodologique dédié à la réduction dimensionnelle, la sélection du meilleur sous-ensemble — universellement désignée sous l’appellation anglo-saxonne de Best Subset Selection — occupe une position canonique et théoriquement pure. Contrairement aux approches heuristiques descendantes, ascendantes ou aux méthodes de régularisation continue fondées sur des pénalités douces, la méthode du meilleur sous-ensemble s’attache à explorer l’intégralité combinatoire des configurations possibles de prédicteurs. Elle résout directement le problème d’optimisation mathématique sous contrainte de parcimonie stricte, formulé au moyen de la norme L0, garantissant ainsi l’identification du modèle optimal pour chaque degré de complexité fixé a priori. Cette rigueur absolue confère au procédé le statut de standard d’or (gold standard), servant d’étalon théorique à l’aune duquel toutes les méthodes approchées sont invariablement évaluées.
Ce guide exhaustif se propose d’explorer en profondeur l’ensemble des dimensions conceptuelles, mathématiques, algorithmiques et appliquées du Best Subset Selection. De l’élucidation du compromis biais-variance à la formalisation des critères d’information pénalisés tels que l’AIC, le BIC et le Cp de Mallows, en passant par la gestion des coûts combinatoires exponentiels et l’implémentation industrielle sous les écosystèmes R et Python, ce traité offre une analyse rigoureuse destinée aux chercheurs, psychomètres et praticiens de l’apprentissage statistique. Une attention toute particulière sera portée à la modélisation des sciences cognitives et comportementales, domaine où la nécessité d’une parcimonie intelligible rencontre la complexité subtile des mesures psychométriques.
1. Introduction conceptuelle à la sélection de variables et au principe de parcimonie
1.1 Le rôle fondamental de la sélection d’attributs en modélisation statistique
La formulation théorique de l’apprentissage supervisé repose sur l’hypothèse de l’existence d’une relation fonctionnelle sous-jacente reliant un vecteur de prédicteurs, noté X appartenant à un espace de dimension p, à une variable réponse continue ou catégorielle notée Y, évoluant dans l’espace unidimensionnel des réels. L’objectif fondamental de la démarche statistique consiste à estimer cette fonction au moyen d’un opérateur mathématique paramétrique ou non paramétrique, capable d’extraire la structure génératrice des données. Dans cette optique, la définition formelle de l’espace des prédicteurs détermine le périmètre de l’information accessible à l’algorithme. Toutefois, l’inclusion indiscriminée de toutes les covariables disponibles au sein du modèle engendre une tension immédiate entre la représentativité descriptive sur l’échantillon observé et la capacité de généralisation prédictive sur des données nouvelles non observées.
L’introduction de covariables non informatives — c’est-à-dire de variables aléatoires orthogonales au signal théorique ou ne présentant que des corrélations résiduelles fortuites avec la cible — dégrade gravement les propriétés de l’estimateur. En augmentant artificiellement le nombre de degrés de liberté consommés par la procédure d’ajustement, ces variables parasites gonflent la variance de l’estimation sans apporter de réduction correspondante du biais. D’un point de vue géométrique, l’espace vectoriel engendré par les colonnes de la matrice des données s’élargit inutilement, forçant la surface de régression à s’adapter aux fluctuations stochastiques propres à l’échantillon d’apprentissage, un phénomène désigné sous le terme de surajustement structurel.
Sur le plan épistémologique, cette quête d’économie conceptuelle s’inscrit en droite ligne dans la tradition philosophique du rasoir d’Ockham (lex parsimoniae), stipulant que la pluralité ne doit pas être postulée sans nécessité. Transposé à l’apprentissage supervisé, ce principe dicte qu’entre deux modèles offrant une capacité prédictive rigoureusement comparable, l’analyste doit impérativement privilégier celui reposant sur le sous-ensemble de prédicteurs le plus restreint. La parcimonie statistique ne constitue donc pas une simple commodité algorithmique, mais une condition sine qua non de la découvrabilité des lois scientifiques sous-jacentes et de la robustesse opérationnelle des inférences.
1.2 Le compromis biais-variance face à la dimensionnalité
L’arbitrage entre la complexité structurelle d’un modèle et sa performance empirique est rigoureusement formalisé par la décomposition de l’erreur quadratique moyenne de prédiction (Mean Squared Error ou MSE). Pour une observation non incluse dans l’ensemble d’apprentissage, l’erreur quadratique attendue se décompose analytiquement en trois termes distincts : le carré du biais de l’estimateur, la variance intrinsèque de l’estimateur induite par les fluctuations d’échantillonnage, et la variance du terme d’erreur stochastique, souvent qualifiée de bruit irréductible. Mathématiquement, cette égalité fondamentale s’écrit : E[(Y – f̂(X))²] = (Biais[f̂(X)])² + Var(f̂(X)) + σ², où σ² représente la variance résiduelle incompressible de la nature.
Lorsqu’un modèle linéaire incorpore l’ensemble des p covariables disponibles sans discernement, le biais d’estimation tend vers son minimum théorique si la spécification fonctionnelle est correcte, car le système intègre l’ensemble des médiateurs et modérateurs potentiels du phénomène. En contrepartie, la variance de l’estimateur croît proportionnellement au nombre de paramètres libres à calibrer. Chaque nouveau coefficient estimé agit comme un capteur supplémentaire capturant l’aléa spécifique de l’échantillon. La sélection de variables opère ainsi une réduction délibérée de la complexité structurelle en contraignant certains coefficients à zéro, acceptant potentiellement une fraction infinitésimale de biais supplémentaire en échange d’une réduction drastique de la variance d’estimation.
Cette dynamique devient critique en présence de multicollinéarité, situation où les colonnes de la matrice d’expérience présentent de fortes interdépendances linéaires. Dans ce régime, la matrice de Gram s’approche de la singularité, provoquant une instabilité numérique sévère de l’inversion matricielle standard des moindres carrés. Les variances des coefficients estimés deviennent alors exorbitantes, conduisant à des intervalles de confiance démesurément larges et à des inversions erratiques de signes de paramètres pourtant établis par la littérature. L’élimination systématique des coefficients redondants restaure le conditionnement de la matrice d’information et garantit la reproductibilité des inférences prédictives.
1.3 Importance de la parcimonie en sciences cognitives et comportementales
Dans le domaine des sciences cognitives, de la neuropsychologie et de la psychiatrie computationnelle, l’enjeu de la sélection de variables se double d’une responsabilité éthique et clinique prépondérante. Les protocoles expérimentaux contemporains recueillent couramment des centaines d’indicateurs issus de batteries de tests psychométriques, de biomarqueurs neurophysiologiques et de questionnaires auto-rapportés évaluant des dimensions latentes complexes telles que l’anxiété, la mémoire de travail ou l’inhibition exécutive. Tenter d’expliquer une altération fonctionnelle ou un profil psychopathologique à l’aide d’un modèle contenant des dizaines de dimensions corrélées compromet l’interprétabilité clinique, privant le praticien de leviers d’intervention thérapeutiques intelligibles.
L’utilisation de modèles parcimonieux permet de dissocier les véritables prédicteurs causaux des corrélations fallacieuses (spurious correlations). Dans les grandes cohortes d’études comportementales, des associations statistiques significatives peuvent émerger sous l’effet de biais de confusion non contrôlés ou d’erreurs de mesure partagées propres aux instruments psychométriques. En isolant le sous-ensemble minimal de prédicteurs assurant la quasi-totalité du pouvoir prédictif, l’analyste isole les construits psychologiques qui portent la véritable valeur étiologique, limitant les interprétations pseudoscientifiques dérivées de covariations accidentelles propres à un échantillon spécifique.
Enfin, l’optimisation des architectures prédictives répond à un impératif pragmatique de réduction du fardeau d’évaluation cognitive imposé aux participants ou aux patients. Soumettre un individu souffrant de troubles dépressifs sévères ou d’un déclin cognitif débutant à des batteries d’évaluation exhaustives durant plusieurs heures engendre fatigue, inattention et désengagement, dégradant la fidélité même des mesures recueillies. L’identification, par des démarches de type meilleur sous-ensemble, d’une combinaison restreinte de deux ou trois indices clés permettant un diagnostic prédictif de haute précision autorise la conception de protocoles de dépistage rapides, éthiques et parfaitement calibrés pour la pratique hospitalière de routine.
2. Fondements mathématiques et algorithmiques du Best Subset Selection
2.1 Formulation formelle de l’espace des modèles candidats
La recherche du meilleur sous-ensemble s’énonce rigoureusement comme un problème d’optimisation mathématique combinatoire sous contrainte de parcimonie, formalisée par la pseudo-norme L0 du vecteur des coefficients de régression. Considérons le modèle linéaire classique Y = Xβ + ε, où Y est le vecteur de dimension (n × 1) des observations cibles, X est la matrice d’expérience de dimension (n × p), β est le vecteur (p × 1) des paramètres réels à estimer, et ε désigne le vecteur des erreurs résiduelles indépendantes et identiquement distribuées. La pseudo-norme ||β||₀ compte de façon binaire le nombre de composantes non nulles de β, représentant la cardinalité exacte du sous-ensemble de prédicteurs mobilisés.
Le problème canonique du meilleur sous-ensemble pour une taille k donnée (avec 0 ≤ k ≤ p) s’écrit de la manière suivante : minimiser la somme des carrés des résidus ||Y – Xβ||₂² sous la contrainte structurelle ||β||₀ ≤ k. L’espace des modèles candidats explore l’ensemble des sous-ensembles de l’ensemble des indices {1, 2, …, p}. La collection totale des configurations possibles pour une taille fixée k correspond au coefficient binomial C(p, k) = p! / (k! (p – k)!). Lorsque l’on somme ces combinaisons pour toutes les cardinalités possibles de k = 0 jusqu’à k = p, la dimension globale de l’espace de recherche atteint précisément 2ᵖ modèles concurrents.
À l’origine de cette séquence combinatoire trône le modèle nul, traditionnellement désigné par M₀. Ce modèle de référence ne comporte aucun prédicteur dans sa matrice d’expérience et se limite à l’ajustement d’un terme constant, l’interception, estimée par la moyenne empirique de la variable réponse Ȳ. Pour toute dimension supérieure, chaque sous-ensemble d’indices définit un sous-espace vectoriel linéaire propre au sein de l’espace euclidien de dimension n, dans lequel la projection orthogonale du vecteur réponse sera minutieusement calculée et évaluée.
2.2 Le protocole séquentiel d’optimisation en deux étapes
La résolution opératoire du Best Subset Selection procède selon une méthodologie rigoureusement scindée en deux phases algorithmiques distinctes, conjuguant l’optimisation intra-dimensionnelle locale et l’arbitrage inter-dimensionnel global. Cette séparation conceptuelle est fondamentale pour éviter la contamination entre l’ajustement mécanique aux données d’apprentissage et l’évaluation de la généralisation statistique :
- Étape 1 : Optimisation intra-dimensionnelle. Pour chaque entier k appartenant à l’intervalle {1, 2, …, p}, l’algorithme génère exhaustivement les C(p, k) combinaisons distinctes de k prédicteurs. Pour chaque combinaison spécifique, le modèle linéaire associé est ajusté par la méthode des moindres carrés ordinaires. On calcule alors la somme des carrés des résidus (RSS) ou le coefficient de détermination R². Le modèle présentant la variance inexpliquée minimale parmi les C(p, k) configurations est formellement désigné comme le représentant optimal de dimension k, noté canoniquement Mₖ. À l’issue de cette étape, l’analyste dispose d’une famille réduite de (p + 1) modèles optimaux : {M₀, M₁, M₂, …, Mₚ}.
- Étape 2 : Sélection inter-dimensionnelle globale. La seconde phase consiste à désigner le modèle définitif au sein de la séquence {M₀, M₁, …, Mₚ}. Cette décision ne peut en aucun cas reposer sur la minimisation brute de la RSS d’apprentissage, cette dernière étant une fonction nécessairement monotone décroissante de la dimension k. L’arbitrage final requiert l’application de critères pénalisés explicitement pour la complexité structurelle (tels que le Critère d’Information d’Akaike, le critère bayésien ou le Cₚ de Mallows) ou l’estimation directe de l’erreur de généralisation via des protocoles de rééchantillonnage empirique par validation croisée.
Cette distinction méthodologique garantit que la comparaison finale ne souffre pas de l’optimisme trompeur associé à l’augmentation mécanique de la dimensionnalité, protégeant l’analyste contre l’illusion de performance prédictive au sein du corpus d’apprentissage.
2.3 Hypothèses sous-jacentes du modèle linéaire classique
L’optimalité théorique conférée aux solutions identifiées à chaque niveau k repose intégralement sur la validité des postulats classiques du modèle de régression linéaire gaussien. Ces hypothèses conditionnent la distribution d’échantillonnage des estimateurs des moindres carrés et la légitimité probabiliste des critères d’arbitrage. En premier lieu, les perturbations résiduelles εᵢ doivent suivre une distribution normale de moyenne nulle et posséder une variance conditionnelle constante σ² sur l’ensemble du domaine des prédicteurs, satisfaisant ainsi l’exigence formelle d’homoscédasticité. Les écarts substantiels à la normalité ou la présence d’hétéroscédasticité déforment la projection géométrique et compromettent l’estimation non biaisée de la variance globale.
En second lieu, l’hypothèse d’indépendance conditionnelle des observations garantit que la matrice de covariance du terme d’erreur est purement diagonale : Var(ε | X) = σ² Iₙ. En cas d’autocorrélation temporelle ou de corrélation intra-classe (dans les plans factoriels hiérarchiques ou longitudinaux), la somme des carrés résiduels sous-évalue systématiquement l’incertitude réelle, ce qui biaise l’identification du sous-ensemble optimal en favorisant artificiellement des combinaisons de covariables surspécifiées.
Enfin, la validité computationnelle de l’algorithme exige que pour chaque sous-ensemble évalué de taille k, la matrice sous-jacente des prédicteurs Xₖ soit de rang plein colonne, c’est-à-dire rang(Xₖ) = k, avec impérativement la contrainte géométrique n > k. Si cette condition est respectée, la matrice de Gram restreinte (Xₖᵀ Xₖ) est strictement définie positive, garantissant l’existence et l’unicité de l’estimateur des moindres carrés β̂ₖ = (Xₖᵀ Xₖ)⁻¹ Xₖᵀ Y. La robustesse statistique de cet estimateur au sein des modèles restreints découle directement du théorème de Gauss-Markov, qui affirme que sous ces hypothèses, l’estimateur ainsi construit demeure le meilleur estimateur linéaire non biaisé (BLUE) pour l’espace vectoriel circonscrit.
3. Complexité combinatoire et limites de calculabilité
3.1 Démonstration du comportement exponentiel 2 puissance p
L’élégance conceptuelle du Best Subset Selection se heurte immédiatement à une barrière computationnelle insurmontable dès que la dimension de l’espace des prédicteurs franchit les standards des données contemporaines. Pour quantifier rigoureusement cette explosion, appliquons le théorème du binôme de Newton à la sommation exhaustive de toutes les combinaisons possibles de taille k prélevées parmi p covariables candidates : ∑ₖ₌₀ᵖ C(p, k) = (1 + 1)ᵖ = 2ᵖ. Le cardinal de l’espace de recherche obéit ainsi à une trajectoire géométrique exponentielle pure de base deux.
Cette croissance se traduit par une multiplication vertigineuse des inversions matricielles requises. Si pour p = 10, l’espace se limite à 1 024 modèles — charge de travail triviale résolue en quelques millisecondes par n’importe quel processeur contemporain —, pour p = 20, le nombre d’ajustements s’élève à 1 048 576 modèles. À p = 30, le compteur atteint 1 073 741 824 équations matricielles à résoudre. Pour p = 50, le volume d’opérations dépasse 1,12 × 10¹⁵ modèles candidats, ce qui nécessiterait des mois d’exécution continue sur des architectures de calcul distribué à haute performance.
Le point de rupture computationnel pour une station de travail moderne effectuant des inversions via des décompositions QR optimisées se situe classiquement autour de p = 30 à p = 40 variables, à condition de recourir à des implémentations en langage compilé (C ou Fortran). Au-delà de ce seuil, le temps nécessaire pour simplement énumérer les adresses mémoire des structures d’indices dépasse la capacité pratique des systèmes d’exploitation conventionnels, rendant la recherche exhaustive matériellement inaccessible.
3.2 La classification NP-difficile de l’optimisation globale
La barrière algorithmique du Best Subset Selection ne découle pas seulement d’une limitation technologique passagère relative à la vitesse des microprocesseurs, mais prend racine dans la théorie de la complexité computationnelle. Le problème de la régression sous contrainte de parcimonie L0 a été formellement démontré comme appartenant à la classe des problèmes NP-difficiles (NP-hard), notamment dans les travaux séminales de Natarajan (1995). Cette classification implique qu’il n’existe aucun algorithme déterministe connu capable d’identifier la solution globale exacte en un temps polynomial par rapport à la dimension p dans le cas général.
Cette réalité théorique a des conséquences fondamentales sur l’analyse de données modernes à haute dimension, telles que les matrices de transcriptomique génomique (comportant fréquemment p > 20 000 gènes pour quelques centaines d’échantillons) ou l’imagerie par résonance magnétique fonctionnelle (contenant des centaines de milliers de voxels). Dans ces contextes où la dimensionnalité excède de plusieurs ordres de grandeur la taille de l’échantillon (régime p ≫ n), l’application du Best Subset Selection brut relève de l’impossibilité physique absolue.
Cette impasse calculatoire constitue la motivation théorique directe ayant conduit la communauté statistique à concevoir des alternatives opérationnelles. Ces dernières se scindent en deux familles majeures : d’une part, les algorithmes gloutons à complexité polynomiale qui tronquent l’exploration combinatoire par des choix locaux séquentiels, et d’autre part, les méthodes de pénalisation continue convexe, au premier rang desquelles figure le Lasso, qui remplace la contrainte de norme L0 discontinue par une enveloppe convexe sous forme de norme L1.
3.3 Stratégies d’élagage d’arbre et branch-and-bound
Afin de repousser les limites de la faisabilité calculatoire sans renoncer à l’exactitude de la solution globale, des informaticiens et statisticiens ont développé des méthodes d’optimisation combinatoire avancées, fondées sur le paradigme de séparation et évaluation (Branch-and-Bound). L’algorithme le plus emblématique de cette catégorie demeure sans conteste celui formulé par Furnival et Wilson (1974), connu sous le nom de routine Leaps and Bounds. Cette approche permet de découvrir le sous-ensemble minimisant la RSS sans avoir à ajuster explicitement la totalité des 2ᵖ configurations.
Le principe fondamental repose sur l’organisation des sous-ensembles sous la forme d’un arbre d’exploration arborescent et sur une propriété géométrique stricte de la somme des carrés des résidus : l’adjonction d’un prédicteur à un sous-ensemble existant ne peut en aucun cas accroître sa RSS. Par conséquent, la RSS d’un sous-ensemble parent constitue une borne inférieure infranchissable pour l’ensemble des sous-modèles dérivés par expansion de ce nœud. Si, au cours de la descente de l’arbre, la RSS d’un nœud intermédiaire s’avère déjà supérieure à la meilleure RSS identifiée pour un modèle de taille k préalablement visité, il devient inutile d’explorer les branches filles. La coupure (ou élagage) de ce sous-arbre entier permet d’économiser l’évaluation de millions de combinaisons superflues.
L’efficacité de cette stratégie est magnifiée par la vectorisation des opérations sur les matrices de covariances et l’usage de balayages d’élimination de Gauss-Jordan (opérations de sweep). Néanmoins, cette optimisation rencontre des limites empiriques sévères en présence de prédicteurs fortement collinéaires. Dans une telle configuration, la décroissance de la RSS d’un niveau à l’autre devient erratique et marginale, empêchant l’algorithme d’activer précocement ses critères d’élagage. L’arborescence doit alors être visitée dans sa quasi-totalité, ramenant le coût de calcul à sa borne exponentielle originelle.
4. Critères d’évaluation au sein d’une dimension fixée
4.1 La somme des carrés résiduels (RSS)
Au sein d’une dimensionnalité k strictement prédéterminée, le critère universel permettant d’ordonner les C(p, k) modèles candidats et d’extraire le champion local Mₖ est la somme des carrés des résidus (Residual Sum of Squares ou RSS). Analytiquement, si l’on note y le vecteur des observations réelles et ŷ = Xₖ β̂ₖ le vecteur des prédictions générées par la sous-matrice Xₖ, le résidu d’ajustement est défini par le vecteur d’erreur e = y – ŷ. La RSS correspond à la norme euclidienne au carré de ce vecteur d’écarts : RSS = eᵀ e = ∑ᵢ₌₁ⁿ (yᵢ – ŷᵢ)².
Sur le plan géométrique, l’ajustement par les moindres carrés équivaut rigoureusement à projeter orthogonalement le vecteur de données y sur le sous-espace vectoriel de dimension k engendré par les colonnes de Xₖ. L’opérateur de projection orthogonale, communément qualifié de matrice chapeau (hat matrix), s’exprime par la relation matricielle Hₖ = Xₖ(Xₖᵀ Xₖ)⁻¹ Xₖᵀ. Le résidu s’écrit alors e = (Iₙ – Hₖ) y, où (Iₙ – Hₖ) est le projecteur orthogonal sur l’espace orthogonal à l’image de Xₖ. La minimisation de la RSS consiste géométriquement à trouver la distance minimale entre le vecteur réponse et le sous-espace considéré.
Une propriété fondamentale régit l’évolution de la RSS : la monotonie décroissante par rapport à l’inclusion de variables. Tout agrandissement d’un sous-espace vectoriel englobe nécessairement l’espace précédent, garantissant que la distance résiduelle ne peut qu’être réduite ou, au pire, demeurer strictement identique. De cette propriété découle une conséquence capitale : la RSS brute constitue un critère absolu d’arbitrage intra-dimensionnel (pour comparer des modèles partageant une même taille k), mais s’avère intrinsèquement incapable d’opérer un arbitrage inter-dimensionnel loyal entre modèles de dimensions distinctes, favorisant invariablement le modèle surparamétré complet Mₚ.
4.2 Le coefficient de détermination empirique R-deux
Directement corrélé à la somme des carrés résiduels, le coefficient de détermination empirique, universellement désigné par R², représente la proportion de la variance totale de la variable cible expliquée par le modèle de régression. Sa formulation algébrique s’appuie sur la décomposition de la somme totale des carrés corrigée de la moyenne (Total Sum of Squares ou TSS) en une composante expliquée par la régression (Explained Sum of Squares ou ESS) et la composante résiduelle : TSS = ESS + RSS, avec TSS = ∑ᵢ₌₁ⁿ (yᵢ – Ȳ)². Le coefficient se formalise alors par le ratio : R² = 1 – (RSS / TSS).
Dans l’espace euclidien centré des observations, le R² possède une interprétation trigonométrique élégante : il correspond exactement au carré du cosinus de l’angle géométrique formé entre le vecteur centré de la réponse empirique (y – Ȳ) et le vecteur des prédictions ajustées (ŷ – Ȳ). Un R² tendant vers l’unité traduit un alignement vectoriel quasi-parfait, signifiant que la surface de régression capture fidèlement la variabilité spatiale des points expérimentaux.
Cependant, l’utilisation naïve du R² comme boussole de sélection de modèles constitue l’un des pièges les plus documentés de l’inférence statistique. Puisque la TSS est une constante invariante déterminée exclusivement par la distribution brute de la cible y, et que la RSS décroît mécaniquement à chaque introduction de variable supplémentaire, le R² augmente de manière inéluctable avec la dimension k. L’adjonction d’une variable générée par un processus de bruit blanc gaussien sans aucun lien causal avec le phénomène élèvera systématiquement le R² de l’échantillon d’apprentissage, créant une illusion mathématique d’amélioration prédictive qui masque une dégradation sévère des capacités de généralisation du système.
4.3 Propriétés statistiques de l’estimateur des moindres carrés restreint
Lorsque le statisticien restreint délibérément l’ajustement à un sous-ensemble k de prédicteurs issus de la matrice globale, l’estimateur restreint obtenu β̂ₖ hérite de propriétés mathématiques rigoureusement décrites par le formalisme de Gauss-Markov. Si le véritable processus générateur de données ne mobilise effectivement que ce sous-ensemble de variables, l’estimateur des moindres carrés restreint est sans biais, c’est-à-dire E[β̂ₖ] = βₖ, et sa matrice de variance-covariance s’établit à Var(β̂ₖ) = σ² (Xₖᵀ Xₖ)⁻¹. L’espérance de cette matrice est étroitement liée au spectre des valeurs propres de la matrice d’information de Fisher restreinte Xₖᵀ Xₖ.
Le conditionnement de cette matrice joue un rôle déterminant dans la précision des estimations. En éliminant les prédicteurs redondants ou quasi-collinéaires du jeu global, l’inversion matricielle restreinte devient nettement mieux conditionnée : le ratio entre la plus grande et la plus petite valeur propre se stabilise, évitant l’amplification artificielle des erreurs de calcul et stabilisant l’estimation des erreurs-types des coefficients résiduels.
Il est toutefois impératif d’analyser le comportement de β̂ₖ dans l’éventualité où des prédicteurs véritablement informatifs auraient été exclus de la sélection, phénomène connu sous le nom d’effet de masquage ou de biais de variable omise. Si la matrice des variables omises Xₒ présente une covariance non nulle avec les variables retenues Xₖ, l’espérance de l’estimateur restreint devient : E[β̂ₖ] = βₖ + (Xₖᵀ Xₖ)⁻¹ Xₖᵀ Xₒ βₒ. Les coefficients des variables retenues absorbent alors mécaniquement l’effet marginal des variables manquantes, entraînant un biais structurel non réductible qui altère l’interprétation causale des estimateurs au profit exclusif de leur capacité d’approximation descriptive locale.
5. Critères de sélection globale pénalisée pour arbitrer entre tailles de modèles
5.1 Le critère d’information d’Akaike (AIC)
Pour dépasser l’incapacité structurelle de la RSS à arbitrer entre des modèles de complexités distinctes, la théorie de l’information appliquée à la statistique propose des métriques d’évaluation pénalisée. Parmi elles, le Critère d’Information d’Akaike (AIC), introduit par Hirotugu Akaike en 1974, repose sur une fondation épistémologique particulièrement féconde : la minimisation de la divergence relative de Kullback-Leibler entre le processus générateur de données véritable (généralement inaccessible) et le modèle probabiliste proposé.
Dans le cadre de la modélisation linéaire sous hypothèse de perturbations gaussiennes indépendantes et homoscédastiques, l’expression générale de l’AIC, formulée à partir de la fonction de log-vraisemblance maximale ln(L), prend la forme analytique suivante : AIC = -2 ln(L) + 2k*, où k* représente le nombre total de paramètres libres estimés par le modèle (incluant les k coefficients de régression, le terme d’interception et l’estimateur de la variance résiduelle σ², soit k* = k + 2). En exprimant la log-vraisemblance en fonction de la somme des carrés résiduels, la formule devient : AIC = n ln(RSS / n) + 2k + C, où C est une constante d’échelle additive indépendante des modèles en compétition.
L’AIC possède la propriété théorique d’efficacité asymptotique : lorsque la taille de l’échantillon tend vers l’infini, l’AIC sélectionne le modèle dont l’erreur quadratique moyenne de prédiction hors échantillon est minimale, en supposant que le modèle véritable ne figure pas strictement dans l’ensemble candidat ou soit de dimension infinie. Cependant, en contexte d’échantillons finis, le terme de pénalité linéaire 2k s’avère relativement modéré, induisant une tendance documentée à la sursélection (overfitting), ce qui conduit fréquemment l’algorithme à retenir des prédicteurs superficiels apportant des gains résiduels marginaux.
5.2 Le critère d’information bayésien (BIC / SC)
Formulé par Gideon Schwarz en 1978, le Critère d’Information Bayésien (Bayesian Information Criterion ou BIC), également désigné sous le nom de critère de Schwarz, adopte une posture méthodologique radicalement différente, enracinée dans le paradigme probabiliste de l’inférence bayésienne. L’objectif mathématique du BIC est de maximiser la probabilité a posteriori conditionnelle qu’un modèle candidat soit le modèle rigoureusement « vrai », sous l’hypothèse d’une distribution a priori non informative uniforme sur l’ensemble de l’espace des modèles concurrents.
L’équation canonique du BIC s’énonce selon une structure formellement analogue à celle de l’AIC, mais introduit une pondération de pénalité dépendante de la volumétrie de l’échantillon d’apprentissage : BIC = -2 ln(L) + k ln(n), ce qui, sous l’approximation par la régression linéaire aux moindres carrés ordinaires, s’exprime par : BIC = n ln(RSS / n) + k ln(n). Le multiplicateur de complexité passe ainsi de la constante 2 dans l’AIC au logarithme naturel de la taille de l’échantillon ln(n) dans le BIC.
Dès lors que n ≥ 8, la valeur de ln(n) surpasse strictement la pénalité de 2 imposée par l’AIC (par exemple, pour n = 1 000, ln(1 000) ≈ 6,91). Le BIC inflige par conséquent un coût computationnel bien plus sévère à chaque paramètre additionnel. Cette sévérité accrue lui confère une propriété asymptotique fondamentale : la consistance forte. Si le modèle générateur réel se trouve inclus dans la famille des sous-ensembles candidats, la probabilité que le BIC sélectionne ce modèle exact converge vers 1 lorsque n tend vers l’infini. En pratique clinique et appliquée, le BIC est systématiquement préféré lorsque l’interprétabilité et l’élimination des faux positifs priment sur la sensibilité prédictive pure.
5.3 Le Cp de Mallows et ses relations structurelles
Formulé par Colin Lingwood Mallows en 1973, le statistique Cp de Mallows a été spécialement conçu pour évaluer la qualité des régressions linéaires par la quantification de l’erreur quadratique moyenne standardisée totale de prédiction sur l’ensemble des données d’apprentissage. Contrairement aux démarches fondées sur la vraisemblance pure, le Cₚ évalue directement le compromis entre la distorsion induite par les variables omises et l’inflation de variance causée par les prédicteurs surabondants.
L’expression analytique du Cₚ pour un sous-modèle de cardinalité k s’établit selon la formule : Cₚ = (RSSₖ / σ̂²) – n + 2k, où RSSₖ désigne la somme des carrés résiduels du modèle restreint évalué, et σ̂² représente l’estimateur sans biais de la variance résiduelle issu du modèle saturé complet intégrant la totalité des p covariables : σ̂² = RSSₚ / (n – p – 1). L’examen géométrique et analytique de cette relation révèle que si un sous-modèle ne souffre d’aucun biais systématique notable (c’est-à-dire si les variables omises sont dénuées d’influence réelle), l’espérance mathématique de sa RSSₖ est égale à (n – k)σ². En injectant cette égalité dans la formule, on constate que l’espérance théorique du critère devient E[Cₚ | biais = 0] = k.
Le protocole de décision graphique associé au Cₚ consiste à tracer les valeurs de Cₚ calculées en fonction de la taille k et à superposer la ligne d’identité bissectrice Cₚ = k. Les modèles se positionnant nettement au-dessus de cette droite souffrent d’un biais d’omission substantiel, tandis que les modèles se situant à proximité immédiate ou en dessous de la droite représentent des architectures sans biais notable. L’analyste retiendra le modèle affichant la plus petite valeur de Cₚ tout en maintenant une proximité étroite avec la droite d’égalité. D’un point de vue structurel, il est aisé de démontrer qu’en contexte gaussien, la minimisation du Cₚ de Mallows est rigoureusement équivalente à la minimisation de l’AIC d’Akaike.
6. Validation croisée et ajustements empiriques de la performance
6.1 Le coefficient de détermination ajusté
Afin de corriger l’inflation artificielle inhérente au R² classique sans devoir recourir à la formulation théorique des critères d’information, le statisticien Robert Wherry et ses contemporains ont introduit le coefficient de détermination ajusté, conventionnellement noté R̄² ou R²_adj. Ce critère introduit une compensation mécanique en pondérant la somme des carrés résiduels et la somme totale des carrés par leurs degrés de liberté respectifs, rééquilibrant ainsi l’analyse statistique.
L’expression algébrique formelle de cet indicateur s’articule comme suit : R̄² = 1 – [ (RSS / (n – k – 1)) / (TSS / (n – 1)) ] = 1 – (1 – R²) [ (n – 1) / (n – k – 1) ], où n est la taille d’échantillon et k le nombre de prédicteurs inclus dans le sous-ensemble. Le ratio (n – 1) / (n – k – 1) agit comme un coefficient multiplicateur strictement supérieur à l’unité, qui augmente à mesure que de nouveaux paramètres sont intégrés au modèle. Dès lors, l’adjonction d’une covariable additionnelle n’augmente le R̄² que si la réduction relative de la RSS compense largement la perte du degré de liberté résiduel.
Contrairement au R² brut, le profil d’évolution du R̄² en fonction de la cardinalité k est non monotone : il présente une trajectoire concave caractérisée par une ascension rapide, un sommet identifiant la parcimonie optimale, puis une phase descendante marquant l’entrée dans la zone de surparamétrisation inutile. Cependant, sur le plan de la théorie de la décision statistique, le R̄² présente des limites notables : il ne dérive pas d’un principe d’optimisation formel de la fonction de perte prédictive et possède une propension à retenir un nombre de variables systématiquement supérieur aux modèles désignés par le BIC ou la validation croisée.
6.2 Protocole de validation croisée k-fold pour la sélection finale
L’évaluation par validation croisée à K plis (K-fold Cross-Validation) constitue l’alternative non paramétrique par excellence aux critères analytiques d’information pénalisée. Son principe opérationnel repose sur le partitionnement stochastique de l’échantillon empirique initial D en K blocs (ou plis) disjoints de tailles approximativement équivalentes. La robustesse de cette approche réside dans son affranchissement total de toute hypothèse distributionnelle quant à la gaussianité, l’homoscédasticité ou l’indépendance paramétrique des résidus.
Le protocole de sélection s’opère selon une rigueur méthodologique impérative pour éviter la fuite d’information (data leakage) :
- Pour chaque sous-ensemble de taille k ∈ {1, …, p}, le modèle champion Mₖ est évalué à travers K itérations successives.
- À chaque itération j ∈ {1, …, K}, le pli j est préservé comme ensemble de validation indépendant, tandis que les (K – 1) plis restants sont agrégés pour former l’ensemble d’entraînement temporaire.
- L’algorithme de régression linéaire est calibré sur les données d’entraînement, puis les paramètres obtenus sont appliqués aux prédicteurs du pli de test pour générer les prédictions ŷᵢ.
- L’erreur quadratique moyenne de validation hors échantillon est alors calculée : MSE_test(j) = (1 / nⱼ) ∑ᵢ (yᵢ – ŷᵢ)².
L’estimation globale de l’erreur de généralisation pour chaque taille de modèle k s’obtient par la moyenne arithmétique : CV₍ₖ₎ = (1 / K) ∑ⱼ₌₁ᴷ MSE_test(j). Pour arbitrer définitivement entre les tailles de modèles, les praticiens aguerris appliquent la célèbre règle du décalage standard (one-standard-error rule), popularisée par Leo Breiman. Cette heuristique consiste à calculer l’erreur-type empirique de la distribution des scores de validation croisée et à retenir le modèle le plus parcimonieux (c’est-à-dire affichant la plus faible valeur de k) dont l’erreur moyenne se situe à l’intérieur d’un intervalle d’un écart-type au-dessus de l’erreur minimale absolue observée.
6.3 Validation croisée répétée et Leave-One-Out (LOOCV)
Le cas limite où le nombre de plis K est fixé rigoureusement égal à l’effectif d’échantillonnage n est désigné sous le terme de validation croisée avec exclusion unitaire (Leave-One-Out Cross-Validation ou LOOCV). Dans ce cadre analytique, le modèle est ajusté n fois sur (n – 1) observations et testé sur l’unique individu exclu. L’estimateur LOOCV offre l’avantage théorique majeur de ne présenter quasiment aucun biais d’estimation de la fonction d’erreur globale, puisque chaque ensemble d’entraînement mobilise une taille d’échantillon quasi-identique au jeu de données complet.
Bien que cette procédure paraisse computationnellement rédhibitoire (exigeant n réajustements matriciels pour chacun des modèles candidats), la modélisation linéaire bénéficie d’un raccourci analytique spectaculaire découlant des propriétés de la matrice chapeau H. La somme des carrés des résidus de prédiction (statistique PRESS) peut être calculée à partir d’un unique ajustement du modèle sur l’ensemble complet, grâce à l’identité suivante : yᵢ – ŷ₋ᵢ = eᵢ / (1 – hᵢᵢ), où eᵢ = yᵢ – ŷᵢ désigne le résidu d’ajustement standard et hᵢᵢ représente le i-ème élément diagonal de la matrice chapeau H = X(Xᵀ X)⁻¹ Xᵀ, correspondant au levier statistique de l’observation i. L’erreur de LOOCV s’obtient alors instantanément par : LOOCV₍ₖ₎ = (1 / n) ∑ᵢ₌₁ⁿ [ eᵢ / (1 – hᵢᵢ) ]².
Cependant, le LOOCV souffre d’une variance d’estimation souvent substantielle : dans la mesure où les n ensembles d’apprentissage sont quasi identiques (ne différant que d’une seule observation), les modèles ajustés à chaque itération sont très fortement corrélés entre eux. La moyenne de variables corrélées affichant une variance supérieure à la moyenne de variables décorrélées, la validation croisée à 10 plis (10-fold CV) offre généralement un bien meilleur compromis biais-variance que le LOOCV pour l’arbitrage final entre sous-ensembles optimaux.
7. Démonstration manuelle et analytique avec un système à trois prédicteurs
7.1 Définition du problème expérimental tridimensionnel
Pour illustrer de manière tangible les mécanismes intimes de la sélection du meilleur sous-ensemble, développons une démonstration analytique complète à partir d’un système factoriel à trois covariables continues X = {x₁, x₂, x₃} et une variable cible y, mesurées sur un échantillon contrôlé de n = 50 participants. Cet environnement contrôlé génère un espace combinatoire complet de 2³ = 8 configurations de modèles distinctes, autorisant un examen exhaustif pas à pas sans opacité algorithmique.
La structure des relations de covariance au sein de notre jeu d’essai est délibérément paramétrée pour reproduire les écueils réels de la modélisation empirique : la variable x₁ possède une corrélation substantielle directe avec la cible (r(x₁, y) = 0,65) ; la variable x₂ présente une corrélation modérée avec la cible (r(x₂, y) = 0,45), mais se trouve partiellement collinéaire avec x₁ (r(x₁, x₂) = 0,50) ; enfin, la variable x₃ a été synthétisée comme une variable de bruit stochastique pur, indépendante de y (r(x₃, y) = 0,04) et orthogonale aux deux autres régresseurs.
Le tableau des huit modèles candidats qui s’affrontent s’établit ainsi selon leur complexité structurelle croissante :
- Taille 0 : Le modèle nul M₀ : ŷ = β₀ (constante d’interception seule).
- Taille 1 : Trois modèles univariés : M_{1a} : ŷ = β₀ + β₁ x₁ ; M_{1b} : ŷ = β₀ + β₂ x₂ ; M_{1c} : ŷ = β₀ + β₃ x₃.
- Taille 2 : Trois modèles bivariés : M_{2a} : ŷ = β₀ + β₁ x₁ + β₂ x₂ ; M_{2b} : ŷ = β₀ + β₁ x₁ + β₃ x₃ ; M_{2c} : ŷ = β₀ + β₂ x₂ + β₃ x₃.
- Taille 3 : Le modèle saturé complet M₃ : ŷ = β₀ + β₁ x₁ + β₂ x₂ + β₃ x₃.
La somme totale des carrés de l’expérience, calculée pour la variable réponse centrée, s’élève à TSS = 150,00.
7.2 Ajustement et compétition des modèles de taille 1 et de taille 2
La première phase d’optimisation locale requiert le calcul systématique des régressions pour chaque modèle univarié de taille 1. L’ajustement du modèle M_{1a} utilisant le prédicteur x₁ aboutit à une somme des carrés résiduels RSS = 86,63, ce qui correspond à un R² = 1 – (86,63 / 150,00) = 0,422. Le modèle M_{1b} basé sur x₂ génère une RSS = 119,63 (R² = 0,202). Enfin, le modèle M_{1c} fondé sur le bruit x₃ affiche une RSS = 149,76 (R² = 0,002). L’application du critère de minimisation intra-dimensionnelle de la RSS conduit sans ambiguïté à sélectionner M_{1a} comme le champion absolu de dimension un, désormais noté officiellement M₁ (avec RSS(M₁) = 86,63).
La procédure examine ensuite de façon exhaustive les trois combinaisons de cardinalité deux. L’ajustement de M_{2a}, combinant x₁ et x₂, exploite la variance complémentaire non redondante des deux facteurs et ramène la RSS à 74,10, portant le coefficient de détermination à R² = 0,506. L’évaluation de M_{2b} combinant x₁ et la variable bruitée x₃ produit une RSS = 86,52 (R² = 0,423), montrant une diminution marginale de la RSS purement accidentelle par rapport à x₁ seul. Enfin, M_{2c} associant x₂ et x₃ délivre une RSS = 119,41 (R² = 0,204). Le champion indiscutable de dimension deux est donc M_{2a}, que l’on consigne comme le modèle M₂ (avec RSS(M₂) = 74,10).
L’observation analytique des sauts de performance met en évidence le phénomène des rendements marginaux décroissants : le passage de la taille 0 (RSS(M₀) = 150,00) à la taille 1 engendre une réduction massive de la RSS de ΔRSS = 63,37. L’adjonction d’un deuxième paramètre informatif au sein de M₂ permet une diminution résiduelle additionnelle substantielle de ΔRSS = 12,53. En revanche, l’intégration finale du troisième prédicteur dans le modèle complet M₃ (qui intègre x₁, x₂, x₃) n’aboutit qu’à une RSS(M₃) = 73,95, soit une réduction résiduelle dérisoire de ΔRSS = 0,15, caractéristique typique de l’absorption d’une composante d’erreur stochastique.
7.3 Arbitrage final entre M0, M1, M2 et M3
L’étape ultime du protocole consiste à confronter les quatre champions issus de la première phase : la séquence {M₀, M₁, M₂, M₃}. Dressons le tableau analytique comparatif de leurs indicateurs de performance pénalisée pour n = 50 observations :
- Modèle M₀ (k=0) : RSS = 150,00 | R² = 0,000 | R̄² = 0,000 | AIC = 56,93 | BIC = 58,85 | Cp = 51,46
- Modèle M₁ (k=1) : RSS = 86,63 | R² = 0,422 | R̄² = 0,410 | AIC = 29,48 | BIC = 33,30 | Cp = 10,36
- Modèle M₂ (k=2) : RSS = 74,10 | R² = 0,506 | R̄² = 0,485 | AIC = 23,61 | BIC = 29,35 | Cp = 2,10
- Modèle M₃ (k=3) : RSS = 73,95 | R² = 0,507 | R̄² = 0,475 | AIC = 25,51 | BIC = 33,16 | Cp = 4,00
L’analyse des critères pénalisés permet de trancher avec clarté. L’examen du R̄² ajusté montre qu’il atteint son maximum structurel pour M₂ (0,485), avant de régresser pour M₃ (0,475) en raison de la perte non compensée d’un degré de liberté. Le critère Cₚ de Mallows affiche sa valeur minimale pour M₂ (Cₚ = 2,10), valeur extrêmement proche de sa cardinalité théorique sans biais k = 2, tandis que M₃ subit une pénalisation mécanique qui le porte à Cₚ = 4,00. L’AIC et le BIC convergent tous deux vers le même diagnostic en désignant sans équivoque M₂ comme la configuration optimale, affichant des valeurs minimales respectives de 23,61 et 29,35.
Le diagnostic de multicollinéarité résiduelle au sein de la sélection M₂ confirme la pertinence de ce choix : le facteur d’inflation de la variance (VIF) des prédicteurs x₁ et x₂ s’établit à VIF = 1 / (1 – 0,50²) = 1,33, une valeur largement inférieure au seuil critique conventionnel de 5 ou 10. La méthode a donc su extraire l’architecture d’explication optimale en éliminant précisément la variable parasite x₃, tout en préservant l’équilibre statistique fondamental des estimateurs retenus.
8. Application empirique en modélisation psychologique et psychométrique
8.1 Étude de cas : Prédiction des niveaux d’épuisement professionnel (Burnout)
Pour illustrer la portée concrète du Best Subset Selection dans le champ de la recherche appliquée, étudions un protocole en psychologie organisationnelle et ergonomie cognitive dédié à la modélisation de l’épuisement professionnel (syndrome de burnout). L’investigation mobilise une cohorte clinique de n = 450 personnels soignants hospitaliers. La variable dépendante continue à prédire correspond au score standardisé de la dimension « Épuisement Émotionnel » de l’inventaire canonique de Maslach (Maslach Burnout Inventory ou MBI), mesuré sur une échelle d’intensité continue de 0 à 54 points.
L’espace initial des prédicteurs rassemble sept indicateurs psychologiques, comportementaux et environnementaux majeurs, standardisés sous forme de scores z :
- x₁ : La charge mentale perçue (évaluée via l’échelle subjective multidimensionnelle NASA-TLX).
- x₂ : Le soutien social et managérial organisationnel (score agrégé de l’échelle de Karasek).
- x₃ : Le névrosisme / instabilité émotionnelle (trait de personnalité extrait de l’inventaire NEO-PI-R).
- x₄ : L’indice objectif de fragmentation du sommeil (mesuré par actigraphie connectée nocturne).
- x₅ : L’extraversion (dimension complémentaire issue du modèle en cinq facteurs).
- x₆ : L’ancienneté professionnelle cumulée au sein de l’unité de soins (exprimée en années).
- x₇ : Le sentiment d’auto-efficacité professionnelle (mesuré par l’échelle de Bandura).
Le traitement préliminaire des données a consisté en une imputation multiple par équations enchaînées (MICE) pour traiter un taux résiduel de 2,3 % de données manquantes sur les mesures actigraphiques, suivie d’un centrage et d’une réduction systématique des variables. L’algorithme Best Subset Selection a ensuite été appliqué pour explorer l’ensemble des 2⁷ = 128 régressions possibles afin d’extraire la structure d’évaluation psychométrique la plus parcimonieuse possible, minimisant l’erreur de détection clinique sans imposer une lourdeur diagnostique inutile.
8.2 Interprétabilité cognitive et validité de construit des sous-ensembles
L’application du critère BIC et de la validation croisée 10-fold identifie de manière convergente un sous-ensemble optimal constitué de quatre prédicteurs majeurs (M₄) : la charge mentale (x₁), le soutien organisationnel (x₂), le névrosisme (x₃) et la fragmentation du sommeil (x₄). L’adjonction de l’extraversion (x₅), de l’ancienneté (x₆) et de l’auto-efficacité (x₇) n’apporte aucun gain d’information statistiquement discernable sur l’échantillon de validation, leurs coefficients affichant des intervalles de confiance englobant zéro une fois les quatre premiers facteurs contrôlés.
L’évaluation épistémologique de cette sélection confirme sa haute validité de construit au regard des modèles contemporains de la psychologie de la santé au travail, notamment le modèle Exigences-Ressources du Travail (Job Demands-Resources model de Demerouti et Bakker). La sélection automatique isole avec une cohérence remarquable deux « exigences professionnelles » majeures (la charge cognitive et la dégradation de la récupération physiologique via le sommeil) et les met en balance avec une « ressource externe » fondamentale (le soutien managérial protecteur) et une « vulnérabilité psychologique interne » reconnue (le névrosisme prédisposant à la réactivité négative face au stress).
Cette concordance théorique démontre la supériorité de l’interprétabilité des modèles parcimonieux extraits par sélection de variables sur les approches de réduction de dimension par projection géométrique pure, telles que l’Analyse en Composantes Principales (ACP). L’ACP aurait combiné les sept variables au sein d’axes composites artificiels mêlant charge de travail, traits de personnalité et années de pratique, rendant les leviers d’action managériale flous et conceptuellement indéchiffrables pour les directions des ressources humaines et la médecine du travail.
8.3 Résolution des artefacts de collinéarité dans les batteries de tests
L’un des apports méthodologiques majeurs de la sélection exhaustive dans les batteries de tests psychométriques réside dans sa capacité à résoudre les phénomènes de quasi-isomorphisme conceptuel. En psychométrie différentielle, différentes échelles construites par des laboratoires distincts partagent fréquemment des formulations sémantiques redondantes. Il n’est pas rare de constater des intercorrélations empiriques excédant r = 0,70 entre des sous-facettes de l’anxiété, de l’affectivité négative et de l’auto-efficacité altérée. L’inclusion simultanée de ces construits au sein d’une régression ordinaire engendre des artefacts de suppression statistique et déstabilise l’estimation des poids bêta.
En confrontant rigoureusement les combinaisons de prédicteurs, l’algorithme Best Subset Selection arbitre automatiquement entre ces dimensions redondantes en ne retenant que le représentant démontrant la liaison fonctionnelle la plus pure avec la variable dépendante, écartant les duplications structurales sans intervention subjective de l’opérateur. La stabilité de cette sélection a été éprouvée en subdivisant la cohorte selon des variables modératrices potentielles (groupes d’âges, départements de soins intensifs versus services de rééducation) : la même configuration à quatre variables s’est systématiquement révélée optimale, confirmant la robustesse transculturelle et environnementale du modèle extrait.
D’un point de vue éthique et opérationnel, ce résultat a permis de condenser le protocole de dépistage préventif hospitalier. En supprimant les inventaires d’extraversion et d’auto-efficacité, la durée de passation du protocole d’évaluation est passée de quarante-cinq minutes à moins de douze minutes par collaborateur, éliminant l’attrition psychologique des répondants et favorisant un suivi longitudinal à grande échelle.
9. Implémentation technique et calculatoire sous R
9.1 Exploitation du package leaps et de la fonction regsubsets
Au sein de l’environnement de calcul statistique R, l’outil de référence pour exécuter la sélection du meilleur sous-ensemble est le package hautement optimisé leaps, développé par Thomas Lumley. Ce composant logiciel implémente au niveau de sa couche système l’algorithme compilé en Fortran de Furnival et Wilson, autorisant l’exploration rapide d’espaces combinatoires pouvant atteindre une trentaine de variables prédictives.
La fonction maîtresse du package s’intitule regsubsets(). Sa configuration élémentaire s’appuie sur la syntaxe canonique des formules statistiques sous R, complétée par des hyperparamètres critiques régissant le périmètre de l’exploration :
formula: la relation symbolique définissant la variable dépendante et les prédicteurs (par exemple :burnout ~ .pour mobiliser l’ensemble des colonnes disponibles).data: le data.frame contenant le tableau rectangulaire des données expérimentales.nvmax: la cardinalité maximale autorisée pour les sous-ensembles (fixée par défaut à 8, devant impérativement être augmentée pour correspondre au nombre total p de variables candidates).method: la stratégie algorithmique de parcours (le paramètre"exhaustive"active formellement le Best Subset Selection exact, tandis que"forward"ou"backward"déclenchent les approximations gloutonnes séquentielles).
L’objet résultant de l’exécution encapsule la totalité des matrices d’information associées aux modèles optimaux locaux Mₖ. L’extraction de son résumé via summary() donne accès à une matrice booléenne d’inclusion logique identifiant les variables sélectionnées pour chaque dimension, ainsi qu’aux vecteurs d’indicateurs pénalisés : rsq, rss, adjr2, cp et bic. L’écosystème propose une fonction de restitution graphique spécialisée, plot.regsubsets(), qui projette les grilles de coefficients matriciels associées à leurs scores d’évaluation ordonnés, visualisant instantanément l’architecture du modèle optimal dominant l’espace de recherche.
9.2 Automatisation du choix du modèle optimal par validation croisée externe
Une lacune structurelle du package leaps réside dans l’absence de méthode intégrée predict() pour les objets de classe regsubsets, empêchant le calcul automatique direct de prédictions sur de nouveaux jeux de données. Pour orchestrer une sélection irréprochable fondée sur la validation croisée à 10 plis, l’analyste doit impérativement concevoir une routine de prédiction personnalisée et encapsuler l’évaluation au sein d’une double boucle programmatique vectorisée.
La fonction de prédiction matricielle requise s’articule autour de l’extraction de la formule du modèle et de la multiplication directe des covariables par le vecteur des paramètres calibrés :
predict.regsubsets <- function(object, newdata, id) {
form <- as.formula(object$call[[2]])
mat <- model.matrix(form, newdata)
coefi <- coef(object, id = id)
xvars <- names(coefi)
mat[, xvars] %*% coefi
}
Une fois cet opérateur établi, la boucle de validation croisée procède à l’affectation aléatoire d’un identifiant de pli (1 à 10) à chaque ligne du tableau de données. À chaque itération, la fonction regsubsets() est exécutée exclusivement sur les observations d’entraînement (constituées de 9 plis sur 10), garantissant l’étanchéité absolue de l’évaluation. Les prédictions sont ensuite générées pour chaque taille k ∈ {1, …, p} sur les données du pli de test isolé, et les résidus quadratiques sont stockés au sein d’une matrice rectangulaire d’erreurs.
La moyenne des colonnes de cette matrice fournit la trajectoire empirique de l’erreur quadratique moyenne de prédiction (MSE) en fonction de la taille du modèle. L’adjonction des barres d’écart-type d’échantillonnage permet d’appliquer graphiquement la règle du one-standard-error, identifiant formellement la dimension optimale qui concilie stabilité prédictive et compacité structurelle.
9.3 Diagnostic approfondi et validation des résidus du modèle extrait
L’extraction algorithmique du meilleur sous-ensemble ne dispense en aucun cas l’analyste d’un diagnostic approfondi des résidus sur le modèle final retenu. Une régression sélectionnée par les procédures d’optimisation peut s’avérer artéfactuelle si son ajustement apparent repose sur un sous-groupe d’observations extrêmes ou si elle viole grossièrement les postulats gaussiens. Une fois le modèle champion de taille k extrait, celui-ci doit être réajusté explicitement via la fonction canonique lm() pour générer l’ensemble des statistiques d’influence et de diagnostic.
Le protocole de validation résiduelle requiert la vérification des points suivants :
- Détection des points leviers et des observations influentes : Le calcul systématique des distances de Cook permet d’isoler toute observation affichant une valeur supérieure au seuil usuel 4 / (n – k – 1), signalant un individu déformant à lui seul la géométrie de la surface d’ajustement.
- Examen des résidus studentisés : L’analyse de la distribution des résidus studentisés externes permet de vérifier la normalité des erreurs et de détecter la présence éventuelle de valeurs aberrantes au-delà de l’intervalle [-3, +3].
- Test d’homoscédasticité : L’application du test de Breusch-Pagan (disponible via la fonction
bptest()du packagelmtest) confirme ou infirme formellement l’hypothèse de variance résiduelle constante. - Facteur d’inflation de la variance (VIF) : L’utilisation de la fonction
vif()du packagecaratteste de l’absence de collinéarité dommageable au sein du sous-ensemble sélectionné, une valeur inférieure à 2,5 garantissant l’orthogonalité quasi-parfaite des prédicteurs.
L’ensemble de ces métriques diagnostiques est consigné au sein d’un tableau récapitulatif formaté rigoureusement selon les standards de publication de l’American Psychological Association (APA 7th edition), incluant les coefficients de régression non standardisés B, les erreurs-types SE, les coefficients standardisés β, les valeurs du test t de Student et les probabilités critiques p associées.
10. Implémentation technique sous Python (Scikit-Learn, Statsmodels, MLxtend)
10.1 Exploitation du module ExhaustiveFeatureSelector de MLxtend
Dans l’écosystème Python orienté vers la science des données, la bibliothèque native Scikit-Learn ne propose pas de classe d’estimation combinatoire exhaustive directe, privilégiant les algorithmes gloutons tels que RFE (Recursive Feature Elimination) ou SequentialFeatureSelector. Pour réaliser un Best Subset Selection strict conforme aux canons statistiques, la bibliothèque spécialisée MLxtend (Machine Learning Extensions), développée par Sebastian Raschka, propose le module wrapper de référence : ExhaustiveFeatureSelector.
Cette classe s’intègre parfaitement aux API de Scikit-Learn en encapsulant n’importe quel estimateur de base (comme LinearRegression()) et en exécutant l’exploration intégrale sous contrainte combinatoire. Les hyperparamètres de paramétrage fondamental comprennent :
estimator: l’algorithme d’apprentissage sous-jacent servant à l’évaluation des sous-espaces.min_featuresetmax_features: les bornes entières définissant l’intervalle des cardinalités explorées (pour un Best Subset complet, fixer min = 1 et max = p).scoring: la fonction de perte ou métrique d’évaluation optimisée à travers les plis (couramment'neg_mean_squared_error'pour la MSE ou'r2'pour le coefficient de détermination).cv: la stratégie de validation croisée interne, acceptant un entier K (par exemple cv = 10) ou un générateur de partitionnement stratifié (KFold).n_jobs: le paramètre de parallélisation computationnelle, permettant de distribuer l’ajustement des milliers de combinaisons matricielles sur l’ensemble des cœurs physiques du processeur (fixé à-1pour saturer les capacités disponibles de la machine).
L’appel de la méthode fit(X, y) déclenche la recherche exhaustive. Une fois l’exploration achevée, l’attribut subsets_ restitue un dictionnaire Pandas exhaustif documentant, pour chaque combinaison évaluée, les indices des variables, les scores individuels obtenus sur chaque pli de validation croisée, le score moyen agrégé et l’erreur-type d’échantillonnage correspondante.
10.2 Construction d’une routine algorithmique native en Python pur
Pour les contextes où l’analyste statistique recherche la traçabilité intégrale de l’inférence théorique et la restitution explicite des critères analytiques formels (AIC, BIC, log-vraisemblance) sans recourir à des bibliothèques tierces, la construction d’une routine algorithmique modulaire exploitant la puissance conjointe de itertools.combinations et de la bibliothèque économétrique Statsmodels représente une démarche pédagogique et opérationnelle exemplaire.
L’architecture algorithmique s’articule autour d’une boucle énumérative générant les sous-ensembles de cardinalité k via le générateur itératif de bas niveau itertools.combinations(colonnes_disponibles, k). Pour chaque sous-ensemble d’attributs, la sous-matrice correspondante est transmise à l’estimateur des moindres carrés ordinaires sm.OLS(y, sm.add_constant(X_subset)). L’ajustement vectorisé permet d’extraire instantanément la log-vraisemblance maximale, le R² centré, le R̄² ajusté ainsi que les valeurs exactes de l’AIC et du BIC calculées au sens de la théorie de la vraisemblance.
L’ensemble de ces informations est agrégé au sein d’une structure de données structurée de type pandas.DataFrame. L’extraction du modèle champion pour chaque strate de dimension s’opère par une simple opération de groupement et de minimisation conditionnelle : df.loc[df.groupby('k')['rss'].idxmin()]. Cette routine native permet de concevoir une classe orientée objet élégante, dont les performances de vitesse s’avèrent hautement compétitives pour p ≤ 20 prédicteurs, bénéficiant des bibliothèques numériques compilées BLAS/LAPACK sous-jacentes à Numpy et Scipy.
10.3 Intégration dans un pipeline de production Scikit-Learn
Dans un contexte de déploiement applicatif ou industriel d’apprentissage automatique, la sélection de variables ne doit jamais opérer comme une étape déconnectée du cycle d’ingénierie des données. L’erreur méthodologique majeure consiste à sélectionner le sous-ensemble de variables sur l’intégralité du jeu de données avant d’opérer la standardisation ou l’évaluation, introduisant un biais optimiste sévère par fuite d’information statistique.
Pour prévenir cette contamination, l’algorithme de sélection exhaustive doit être impérativement inséré au sein d’un Pipeline Scikit-Learn unifié. L’architecture de production connecte séquentiellement un module de centrage-réduction des données (StandardScaler()), le sélecteur combinatoire ExhaustiveFeatureSelector, et enfin le régresseur terminal d’inférence. Dans cette configuration géométrique, lors de chaque découpage de validation croisée, les moyennes et écarts-types de standardisation sont calculés exclusivement sur les observations d’entraînement du pli, garantissant que les variables sélectionnées ne bénéficient d’aucune information sur la structure distributionnelle du pli d’évaluation.
Le pipeline ainsi constitué peut faire l’objet d’une optimisation conjointe de ses hyperparamètres via GridSearchCV, permettant d’arbitrer simultanément la cardinalité des variables et les seuils de pré-conditionnement. Une fois calibré sur l’intégralité du corpus d’entraînement selon les meilleures pratiques, le pipeline complet est sérialisé et figé dans un artefact binaire persistant à l’aide de la bibliothèque joblib, garantissant sa déployabilité immédiate sur des infrastructures de production en microservices ou en traitement par lots.
11. Analyse comparative : Best Subset face aux heuristiques concurrentes
11.1 Confrontation avec la sélection pas à pas ascendante (Forward Selection)
Face à la lourdeur exponentielle 2ᵖ du Best Subset Selection, la communauté statistique a historiquement privilégié des heuristiques polynomiales de type glouton (greedy algorithms). La plus célèbre d’entre elles est la sélection pas à pas ascendante (Forward Stepwise Selection). Cette méthode débute invariablement par le modèle nul M₀, puis évalue à chaque itération l’adjonction individuelle de chacune des variables non encore intégrées. La variable induisant le plus grand gain marginal de variance expliquée (ou la plus faible p-value au test de Fisher partiel) est incorporée de manière irréversible. Le processus se répète jusqu’à l’incorporation séquentielle de la totalité des p variables.
L’avantage computationnel de la sélection ascendante est spectaculaire : au lieu d’explorer 2ᵖ modèles, l’algorithme n’ajuste qu’une fraction restreinte de configurations, égale à 1 + ∑ₖ₌₀ᵖ⁻¹ (p – k) = 1 + p(p + 1) / 2 ≈ O(p²) modèles. Pour une dimension de p = 50, le nombre d’équations à résoudre passe ainsi de 1,12 × 10¹⁵ pour le meilleur sous-ensemble à seulement 1 276 ajustements pour la sélection ascendante, ramenant l’exécution à une fraction de seconde sur un ordinateur portable.
Toutefois, cette accélération radicale s’accompagne d’une pathologie structurelle majeure : la myopie algorithmique. La sélection pas à pas effectue des choix locaux irrévocables sans jamais réévaluer les décisions antérieures. Or, il est mathématiquement démontré qu’un sous-ensemble composé de deux variables conjointement très explicatives peut comporter deux prédicteurs qui, pris individuellement de manière univariée, s’avèrent de piètres compétiteurs. Dans un tel cas, la sélection ascendante ignorera totalement cette combinaison optimale bivariée car aucun de ses composants n’aura réussi à franchir le premier filtre univarié de l’étape 1. La sélection pas à pas ne garantit donc aucunement l’identification de l’optimum global, contrairement au Best Subset.
11.2 Confrontation avec l’élimination pas à pas descendante (Backward Elimination)
Symétrique de la démarche précédente, l’élimination pas à pas descendante (Backward Stepwise Selection ou élagage rétrograde) amorce sa trajectoire au sommet de la complexité, c’est-à-dire au niveau du modèle complet Mₚ intégrant simultanément l’ensemble des p covariables. À chaque itération séquentielle, l’algorithme identifie et retranche le prédicteur dont la contribution marginale est la plus faible, attestée par la statistique t de Student la plus insignifiante ou la plus faible augmentation de la RSS. Ce processus de délestage s’interrompt lorsque tous les prédicteurs résiduels franchissent un seuil de significativité statistique préalablement prescrit.
À l’instar de l’approche ascendante, l’élimination descendante bénéficie d’une complexité polynomiale en O(p²), offrant une alternative computationnelle attractive à l’exploration exhaustive. De surcroît, elle possède une supériorité conceptuelle reconnue sur la sélection ascendante : en démarrant du modèle complet, elle évalue chaque variable conditionnellement à la présence de toutes les autres. Cette formulation permet de préserver des couples ou des triplets de covariables qui n’expriment leur puissance explicative que par des mécanismes d’interaction ou de suppression mutuelle du bruit.
Cependant, l’élimination descendante se heurte à une contrainte mathématique absolue : elle exige impérativement que la taille de l’échantillon soit strictement supérieure au nombre de prédicteurs (n > p). Dans le cas contraire, la matrice de Gram globale n’est pas inversible, empêchant tout ajustement initial du modèle de départ Mₚ. De surcroît, tout comme la sélection ascendante, elle reste une heuristique gloutonne exposée au risque d’éliminations prématurées non récupérables, ce qui la rend incapable de garantir l’optimalité globale que seul le Best Subset Selection délivre avec certitude.
11.3 Confrontation avec la régularisation Lasso (norme L1) et Ridge (norme L2)
L’alternative théorique contemporaine la plus puissante au Best Subset Selection réside dans les méthodes de régression pénalisée, incarnées par le Lasso (Least Absolute Shrinkage and Selection Operator) de Robert Tibshirani (1996) et la régression Ridge d’Hoerl et Kennard (1970). Alors que le Best Subset résout un problème d’optimisation non convexe discontinu basé sur la norme L0, le Lasso réalise une relaxation convexe naturelle en remplaçant la contrainte discontinue par une pénalité fondée sur la norme L1 de la somme des valeurs absolues des coefficients : Minimise ||Y – Xβ||₂² + λ ||β||₁.
Sur le plan géométrique, l’espace des contraintes de la norme L1 forme un polytope à sommets aigus (un hyperoctaèdre dans l’espace des paramètres). Lors de l’expansion quadratique des contours elliptiques de la fonction de perte des moindres carrés, le point de contact tangentiel s’établit préférentiellement sur les coins et les arêtes du polytope, contraignant mécaniquement certains coefficients à s’annuler exactement. Le Lasso opère ainsi simultanément l’estimation continue des paramètres et la sélection parcimonieuse d’attributs. Grâce à la convexité du problème, il est résolu en temps polynomial par des algorithmes de descente par coordonnées cycliques (Coordinate Descent), gérant avec aisance des matrices comprenant des centaines de milliers de variables.
Néanmoins, cette relaxation continue n’est pas sans coût méthodologique. La pénalité L1 applique une force de contraction uniforme qui biaise systématiquement vers zéro l’estimation des coefficients des véritables prédicteurs dominants (biais d’atténuation ou de rétrécissement). Le Best Subset Selection ne souffre d’aucun biais de ce type : une fois le sous-ensemble identifié, les coefficients sont ajustés par moindres carrés purs sans contraction artificielle. Des travaux récents en optimisation conique et programmation mixte en nombres entiers (notamment menés par Dimitris Bertsimas au MIT) ont d’ailleurs conçu des solveurs modernes (tels que L0Learn) qui combinent la rigueur mathématique du Best Subset pur et des relaxations continues pour résoudre des problèmes de sélection L0 exacte pour des dimensions atteignant p = 100 000, ouvrant la voie à une renaissance spectaculaire du meilleur sous-ensemble en haute dimension.
12. Considérations méthodologiques critiques, pièges et inférence post-sélection
12.1 L’illusion de la significativité et l’invalidation des p-values classiques
L’utilisation de la sélection de variables par meilleur sous-ensemble engendre une pathologie inférentielle subtile mais dévastatrice pour la rigueur scientifique : l’invalidation pure et simple de la théorie classique des tests d’hypothèses statistiques. Lorsqu’un chercheur extrait le meilleur sous-ensemble Mₖ et examine la table de résultats standard rapportant les statistiques t de Student et les probabilités critiques associées (p-values), ces indicateurs sont frappés d’un biais de sélection massif qui les rend mathématiquement fallacieux.
La théorie asymptotique standard de la distribution t de Student repose sur l’hypothèse incontournable que le modèle évalué a été spécifié a priori, de manière totalement indépendante des données empiriques ayant servi à l’ajustement. Or, dans le cas du Best Subset Selection, le modèle Mₖ a été extrait précisément parce qu’il maximisait la variance expliquée parmi des milliers ou des millions de sous-ensembles concurrents testés sur le même échantillon. Les estimateurs des moindres carrés ont ainsi été mécaniquement sélectionnés sur la base de leurs résidus les plus faibles.
Ce processus d’exploration exhaustive entraîne une inflation astronomique du taux d’erreur de type I (faux positifs) : la probabilité de découvrir un prédicteur présentant une pseudo-significativité apparente (p < 0,05) en l’absence de tout lien théorique réel approche 100 % dès que l’espace des covariables s’élargit. Les distributions réelles des estimateurs sous sélection conditionnelle ne suivent plus des lois de Student ou de Fisher-Snedecor, mais des distributions asymétriques tronquées à queues lourdes. L’interprétation naïve des p-values post-sélection sans recours à la théorie moderne de l’inférence sélective (Selective Inference, initiée par Taylor, Tibshirani et Berk) relève d’une illusion méthodologique qui compromet la reproductibilité de la science empirique.
12.2 Instabilité structurelle de l’échantillonnage et sensibilité aux données aberrantes
Un défaut intrinsèque du Best Subset Selection réside dans son instabilité structurelle aiguë face aux fluctuations aléatoires de l’échantillonnage. En statistique computationnelle, une procédure est dite instable si une perturbation minime apportée aux données d’apprentissage (comme la suppression, la modification ou l’ajout de quelques observations) aboutit à la sélection d’un sous-ensemble de covariables radicalement différent.
Cette instabilité structurelle est particulièrement prononcée lorsque l’espace des prédicteurs comprend des variables mutuellement corrélées. Dans cette configuration, des combinaisons distinctes de variables peuvent afficher des valeurs de RSS quasiment indiscernables (ne différant que de quelques millièmes). Une infime fluctuation stochastique dans l’échantillon suffit à faire basculer l’algorithme vers une combinaison alternative, provoquant l’éviction brutale d’un ensemble de variables au profit d’un autre. Ce phénomène complique considérablement l’établissement de conclusions phénoménologiques définitives quant à l’importance théorique réelle de chaque covariable.
Pour mesurer et mitiger cette vulnérabilité, il est impératif de coupler la sélection du meilleur sous-ensemble à des protocoles de rééchantillonnage par bootstrap non paramétrique (générant par exemple 1 000 répliques avec remise du jeu de données initial). L’analyste évalue alors l’indice de stabilité en quantifiant la fréquence relative d’inclusion de chaque variable au fil des itérations. Une variable qui n’émerge comme sélectionnée que dans 45 % des répliques bootstrap ne peut être considérée comme un prédicteur structurellement fiable, indépendamment de ses performances au sein du jeu complet d’origine.
Par ailleurs, la sélection combinatoire est hypersensible à l’influence des valeurs aberrantes (outliers) et des observations à fort levier spatial. Un point unique mal conditionné situé aux confins de l’espace des covariables peut modifier artificiellement l’orientation du vecteur résiduel et contraindre l’algorithme à incorporer un prédicteur inutile dans le seul but de compenser géométriquement l’effet de levier de cet individu marginal. L’examen des distances de levier préliminaire constitue donc un prérequis absolu avant tout lancement d’un balayage combinatoire.
12.3 Guide de recommandations pratiques pour les chercheurs et praticiens
Au terme de cette analyse approfondie, l’analyste, le chercheur ou le modélisateur doit naviguer avec discernement pour choisir la méthodologie de sélection la plus adaptée à ses objectifs scientifiques et aux contraintes computationnelles de son jeu de données. L’arbre décisionnel suivant formalise les bonnes pratiques méthodologiques :
- Régime de faible à moyenne dimension (p ≤ 25 – 30) avec n > p : Le Best Subset Selection exact représente le choix méthodologique optimal absolu. L’analyste exploitera le package
leapssous R ouMLxtend/ routines Statsmodels sous Python. L’arbitrage inter-dimensionnel final devra impérativement s’appuyer sur la validation croisée à 10 plis (en mobilisant la règle du one-standard-error) ou sur le critère BIC si l’objectif principal est la parcimonie et l’explication causale. - Régime de moyenne dimension (30 < p ≤ 60) avec n > p : L’exploration combinatoire exhaustive devenant intraitable, le chercheur privilégiera des algorithmes de Branch-and-Bound optimisés ou se tournera vers les variantes pas à pas bidirectionnelles (Stepwise) avec réévaluation systématique des prédicteurs, tout en gardant à l’esprit les risques de myopie locale.
- Régime de haute à très haute dimension (p > n ou p ≫ 1 000) : L’abandon du Best Subset classique est obligatoire. La stratégie de référence consiste à déployer une régularisation continue de type Lasso ou Elastic Net, ou à exploiter les solveurs modernes d’optimisation en nombres entiers pour la norme L0 développés en recherche opérationnelle (tels que la boîte à outils
L0Learn). - Sanctuarisation de la validation : Quelle que soit la méthode mobilisée, un tiers de l’échantillon empirique initial doit être sanctuarisé sous forme d’ensemble de test indépendant et étanche (holdout test set). Aucun paramètre, aucun choix d’hyperparamètre et aucune phase de sélection ne doivent avoir accès à cet ensemble de test. L’erreur de généralisation définitive et les conclusions scientifiques finales seront mesurées exclusivement sur cet échantillon préservé.
- Intégration du savoir disciplinaire : Les algorithmes combinatoires ne doivent jamais remplacer l’expertise théorique du domaine d’investigation. Si la littérature clinique ou cognitive documente solidement l’influence d’un facteur modérateur ou d’un covariable de confusion, ce dernier doit être maintenu dans la structure du modèle comme régresseur forcé, la sélection combinatoire ne s’exerçant que sur les dimensions exploratoires additionnelles.
En respectant ces canons méthodologiques et déontologiques, la sélection du meilleur sous-ensemble s’affirme comme un levier analytique remarquable, capable d’extraire la substantifique moelle des architectures de données empiriques tout en préservant l’élégance de la parcimonie scientifique.
Références
Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6), 716–723. https://doi.org/10.1109/TAC.1974.1100705
Berk, R., Brown, L., Buja, A., Zhang, K., & Zhao, L. (2013). Valid post-selection inference. The Annals of Statistics, 41(2), 802–837. https://doi.org/10.1214/12-AOS1077
Bertsimas, D., King, A., & Mazumder, R. (2016). Best subset selection via a modern optimization lens. The Annals of Statistics, 44(2), 813–852. https://doi.org/10.1214/15-AOS1388
Breiman, L., Friedman, J. H., Olshen, R. A., & Stone, C. J. (1984). Classification and Regression Trees. Wadsworth & Brooks/Cole Advanced Books & Software.
Furnival, G. M., & Wilson, R. W. (1974). Regressions by leaps and bounds. Technometrics, 16(4), 499–511. https://doi.org/10.2307/1267608
Hastie, T., Tibshirani, R., & Friedman, J. (2009). The Elements of Statistical Learning: Data Mining, Inference, and Prediction (2nd ed.). Springer. https://doi.org/10.1007/978-0-387-84858-7
James, G., Witten, D., Hastie, T., & Tibshirani, R. (2021). An Introduction to Statistical Learning with Applications in R (2nd ed.). Springer. https://doi.org/10.1007/978-1-0716-1418-1
Lumley, T. (2020). leaps: Regression Subset Selection (R package version 3.1). Comprehensive R Archive Network (CRAN). https://cran.r-project.org/package=leaps
Mallows, C. L. (1973). Some comments on Cp. Technometrics, 15(4), 661–675. https://doi.org/10.2307/1267380
Natarajan, B. K. (1995). Sparse approximate solutions to linear systems. SIAM Journal on Computing, 24(2), 227–234. https://doi.org/10.1137/S0097539792240406
Raschka, S. (2018). MLxtend: Providing machine learning and data science utilities and extensions to Python’s scientific computing stack. The Journal of Open Source Software, 3(24), 638. https://doi.org/10.21105/joss.00638
Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6(2), 461–464. https://doi.org/10.1214/aos/1176344136
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Methodological), 58(1), 267–288. https://doi.org/10.1111/j.2517-6161.1996.tb02080.x