La modélisation statistique contemporaine se trouve confrontée à un défi épistémologique permanent : concilier l’ajustement empirique aux données observées avec la parcimonie structurelle indispensable à la généralisation prédictive. Dans le cadre de la régression linéaire multiple, l’ajout systématique de covariables tend inévitablement à réduire ou maintenir constante la somme des carrés des résidus, procurant l’illusion fallacieuse d’une amélioration continue de la qualité du modèle. Cette inflation artificielle des performances sur l’échantillon d’apprentissage masque un phénomène délétère bien connu des biostatisticiens et des psychométriciens : le sur-apprentissage ou sur-ajustement (overfitting), caractérisé par une explosion de la variance des estimateurs et une dégradation dramatique du pouvoir prédictif sur des données indépendantes.
Pour arbitrer rigoureusement ce compromis fondamental entre biais et variance, la théorie statistique a vu émerger au début des années 1970 plusieurs critères d’information et de pénalisation. Parmi ces instruments de décision algorithmique, le critère Cp de Mallows occupe une place matricielle. Conçu spécifiquement pour comparer des sous-modèles dérivés d’un modèle complet présumé non biaisé, le Cp évalue l’erreur quadratique totale moyenne de prédiction normalisée en pénalisant explicitement la dimensionnalité de l’espace des prédicteurs. Contrairement aux coefficients de détermination classiques qui récompensent aveuglément la complexité paramétrique, le critère de Mallows offre une boussole analytique permettant d’isoler l’architecture optimale garantissant une estimation robuste et non biaisée des relations sous-jacentes.
L’implémentation de cette métrique au sein de l’environnement de calcul statistique R constitue aujourd’hui une compétence fondamentale pour tout chercheur ou analyste de données opérant en économétrie, en épidémiologie ou en psychométrie. Ce guide exhaustif explore les fondements mathématiques, les subtilités interprétatives et les protocoles logiciels nécessaires pour calculer, visualiser et contextualiser le critère Cp de Mallows à l’aide de routines R éprouvées, allant des bibliothèques spécialisées comme olsrr et leaps jusqu’à l’écriture d’algorithmes vectoriels autonomes.
- 1. Introduction théorique au critère Cp de Mallows en modélisation statistique
- 2. Fondements mathématiques et dérivation du Cp de Mallows
- 3. Le compromis biais-variance et l’interprétation des paramètres
- 4. Environnement informatique R et préparation des données psychométriques
- 5. Ajustement du modèle complet et des sous-modèles via la fonction lm()
- 6. Calcul du Cp de Mallows avec le package olsrr et la fonction ols_mallows_cp()
- 7. Approches alternatives sous R : calcul manuel et packages complémentaires
- 8. Procédure de sélection automatisée par sous-ensembles et comparaison graphique
- 9. Analyse comparative : Cp de Mallows versus AIC, BIC et R-carré ajusté
- 10. Cas pratique approfondi en recherche psychologique et comportementale
- 11. Diagnostics résiduels et validation post-sélection du modèle retenu
- 12. Bonnes pratiques, limites méthodologiques et considérations éditoriales en recherche quantitative
- Références
1. Introduction théorique au critère Cp de Mallows en modélisation statistique
1.1 Origine historique et place dans l’économétrie et la psychométrie
Proposé formellement par le statisticien britannique Colin Lingwood Mallows dans son article séminal publié en 1973 dans la revue Technometrics (« Some Comments on Cp »), le critère Cp est né des réflexions entourant la sélection de variables dans les contextes de régression linéaire gaussienne. À une époque où les capacités de calcul informatique commençaient à permettre l’estimation séquentielle de multiples équations de régression, la communauté scientifique manquait d’une métrique standardisée pour sanctionner l’introduction superflue de régresseurs non informatifs. Mallows a développé cette approche non pas comme un test d’hypothèse formel reposant sur des valeurs critiques rigides, mais comme un dispositif graphique et diagnostique destiné à guider le praticien dans l’exploration raisonnée des sous-espaces vectoriels de prédiction.
Au fil des décennies, cette méthodologie s’est profondément enracinée dans l’économétrie appliquée et la psychométrie quantitative. Dans les sciences du comportement, l’investigation des construits latents implique fréquemment l’administration de batteries de questionnaires générant des dizaines d’indicateurs psychologiques potentiellement redondants. Face à des matrices de données où le nombre de prédicteurs croît rapidement par rapport au nombre de participants, le risque d’incorporer du bruit stochastique dans les équations structurelles devient critique. Le critère de Mallows a offert aux chercheurs un cadre formel pour élaguer les inventaires de personnalité, les échelles cliniques et les indices cognitifs sans sacrifier la précision globale de mesure.
Dans le paysage contemporain de la science des données, marqué par l’hégémonie apparente des algorithmes d’apprentissage automatique de type boîte noire, la pertinence du Cp de Mallows demeure intacte. Alors que les macro-données psychologiques issues du pistage numérique et des enquêtes en ligne massives multiplient les covariables disponibles, la nécessité de préserver des modèles explicatifs hautement interprétables réhabilite les procédures de sélection parcimonieuse. Le Cp de Mallows s’affirme ainsi comme un pont conceptuel fondamental entre la rigueur de l’inférence statistique classique et les impératifs modernes de validation prédictive.
1.2 Définition conceptuelle et utilité en sélection de modèles
Sur le plan conceptuel, le critère Cp de Mallows mesure l’erreur quadratique totale moyenne standardisée des prédictions générées par un modèle restreint contenant un sous-ensemble de variables, comparativement à un modèle saturé ou complet englobant l’intégralité des prédicteurs disponibles. L’intuition fondamentale qui sous-tend ce critère repose sur le constat suivant : l’omission d’un prédicteur véritablement associé à la variable réponse introduit un biais systématique dans l’estimation des coefficients restants, tandis que l’inclusion d’un prédicteur superflu augmente inutilement la variance d’échantillonnage de toutes les prédictions. Le Cp opère ainsi un arbitrage objectif et transparent entre ces deux forces antagonistes.
Dans une démarche de modélisation par régression, l’analyste se trouve confronté à un ensemble de modèles candidats générés soit par des considérations théoriques a priori, soit par des algorithmes combinatoires exhaustifs. L’utilité du Cp réside dans sa capacité à projeter l’ensemble de ces spécifications compétitives sur une échelle unidimensionnelle standardisée, directement comparable au nombre de paramètres estimés. Cette métrique ne requiert pas que les modèles comparés soient strictement emboîtés au sens des tests de rapport de vraisemblance, pourvu qu’ils dérivent tous du même espace de variables défini par le modèle de référence complet.
Il convient néanmoins de souligner la distinction épistémologique cruciale entre la prédiction empirique et l’inférence causale lors de l’usage du Cp de Mallows. Ce critère est explicitement calibré pour minimiser l’erreur quadratique moyenne de prédiction (Mean Squared Error of Prediction). Par conséquent, un modèle sélectionné par le Cp optimise la précision avec laquelle de futures observations issues de la même population pourront être estimées. Cela n’implique nullement que les coefficients de régression retenus possèdent une interprétabilité causale non biaisée, particulièrement si des variables de confusion majeures ont été omises dès la formulation du modèle complet de référence.
1.3 La règle décisionnelle fondamentale du Cp de Mallows
La mise en œuvre pratique du critère de Mallows repose sur une règle de décision élégante et immédiatement visualisable. Soit un modèle candidat comportant p paramètres (incluant usuellement la constante d’ajustement ou intercept). L’espérance mathématique de la statistique Cp démontre que pour un modèle correctement spécifié, exempt de biais d’omission de variables pertinentes, la valeur attendue de Cp est approximativement égale à p. Cette égalité fondamentale, formalisée par l’expression mathématique E[Cp | biais nul] = p, constitue la pierre angulaire de l’interprétation diagnostique.
Dès lors, l’analyste évalue les modèles candidats en confrontant la valeur numérique de leur Cp à leur cardinalité paramétrique p. Deux trajectoires de déviation apparaissent fréquemment dans les applications concrètes :
- Cas où Cp est nettement supérieur à p (Cp >> p) : Cette situation traduit l’existence d’un biais substantiel dû à la sous-spécification du modèle. Le sous-ensemble de variables sélectionné échoue à capturer une part essentielle de la variance structurée du phénomène étudié, conduisant à des prédictions systématiquement décalées par rapport à la réalité.
- Cas où Cp est proche ou inférieur à p (Cp ≤ p) : Cette configuration indique que le biais résiduel est négligeable ou statistiquement indiscernable du bruit d’échantillonnage. Lorsque le Cp tombe légèrement en dessous de p, ce phénomène résulte le plus souvent de fluctuations stochastiques d’échantillonnage favorables dans l’estimation de la variance de l’erreur.
La règle heuristique universellement reconnue en modélisation consiste donc à identifier, parmi l’ensemble des spécifications candidates, les modèles pour lesquels la valeur de Cp est à la fois faible dans l’absolu et extrêmement proche de p. En présence de plusieurs modèles satisfaisant à cette proximité (Cp ≈ p), le principe de parcimonie scientifique enjoint de retenir le modèle présentant la valeur de p la plus faible, réduisant ainsi la complexité cognitive et opérationnelle du système explicatif.
2. Fondements mathématiques et dérivation du Cp de Mallows
2.1 Formulation algébrique standard du critère
Pour appréhender la mécanique interne du critère Cp de Mallows, il est nécessaire de formaliser le cadre de la régression linéaire sous forme matricielle. Considérons un échantillon composé de n observations indépendantes. Le modèle complet, supposé contenir l’ensemble des K régresseurs potentiels plus la constante (soit P = K + 1 paramètres), s’écrit sous la forme vectorielle standard Y = Xβ + ε, où Y est le vecteur de dimension (n × 1) des réponses observées, X est la matrice de dessein de dimension (n × P), β est le vecteur des paramètres de régression et ε est le vecteur des perturbations aléatoires suivant une loi normale multivariée d’espérance nulle et de matrice de covariance σ2In.
Soit désormais un sous-modèle restreint n’incorporant qu’un sous-ensemble de prédicteurs correspondant à p paramètres (avec p < P). La matrice de projection associée à ce sous-modèle projette le vecteur Y sur le sous-espace vectoriel engendré par les colonnes sélectionnées de X, produisant le vecteur des valeurs ajustées Y_chapeaup. La somme des carrés des résidus associée à ce modèle réduit, usuellement désignée par l’acronyme RSSp (pour Residual Sum of Squares), s’exprime formellement par :
RSSp = ∑i=1n (yi – y_chapeaui,p)2 = YT (In – Hp) Y
où Hp = Xp(XpTXp)-1XpT représente la matrice de projection orthogonale (ou matrice chapeau) de rang p. Le critère Cp de Mallows normalise cette quantité résiduelle en exploitant un estimateur sans biais de la variance de l’erreur issu exclusivement du modèle complet, noté s2 (ou σ_chapeau2). L’équation algébrique définitive du critère s’énonce comme suit :
Cp = (RSSp / s2) – n + 2p
Dans cette expression fondamentale, le terme (RSSp / s2) reflète l’adéquation brute du modèle réduit aux données, pondérée par l’incertitude intrinsèque du système. Le terme constant (- n) opère un recentrage dimensionnel, tandis que la composante (+ 2p) agit comme un mécanisme de pénalisation linéaire proportionnel au nombre de coefficients alloués dans le sous-modèle.
2.2 Hypothèses sous-jacentes du modèle complet
La validité mathématique et l’efficacité opérationnelle du critère Cp reposent sur une condition préalable impérative : le modèle complet à partir duquel s2 est dérivé doit être rigoureusement exempt de biais de spécification. En termes statistiques, cela signifie que la matrice de régresseurs du modèle saturé doit inclure toutes les variables explicatives structurales pertinentes pour décrire l’espérance de Y, garantissant ainsi que l’espérance mathématique de la somme des carrés résiduelle du modèle complet divisée par ses degrés de liberté (n – P) constitue un estimateur sans biais et convergent de la variance théorique σ2 :
E[s2] = E[RSSP / (n – P)] = σ2
Si cette hypothèse de base est violée — par exemple si le chercheur a omis des variables déterminantes dès la phase de formulation théorique globale —, l’estimateur s2 absorbera une fraction de la variance systématique non expliquée, conduisant à une surévaluation grossière de la variance de l’erreur résiduelle. Dans une telle conjoncture, le dénominateur de la fraction (RSSp / s2) devient artificiellement élevé, ce qui fausse l’ensemble des calculs de Cp, comprime les valeurs obtenues vers le bas et entraîne une sélection erronée de sous-modèles sous-dimensionnés.
Par ailleurs, les hypothèses canoniques de la régression linéaire doivent être satisfaites : homoscédasticité parfaite (variance constante des résidus à travers toutes les strates de prédiction), indépendance stochastique des perturbations résiduelles et normalité gaussienne des erreurs. De surcroît, la présence de valeurs aberrantes extrêmes (outliers) ou de points à fort effet de levier (high leverage points) dans la matrice X peut altérer drastiquement l’estimation de s2 et de RSSp, déstabilisant complètement la hiérarchie relative des valeurs de Cp calculées sur les différents sous-ensembles.
2.3 Composantes de l’espérance mathématique de Cp
L’élégance théorique de la métrique de Mallows se révèle pleinement à travers la décomposition analytique de son espérance mathématique. Soit l’erreur quadratique totale moyenne de prédiction normalisée, notée Γp (Gamma de p), définie par la sommation des espérances des écarts au carré entre la véritable fonction d’espérance conditionnelle μ = E[Y|X] et les prédictions vectorielles Y_chapeaup issues du sous-modèle :
Γp = (1 / σ2) ∑i=1n E[(y_chapeaui,p – μi)2]
En décomposant cet écart en termes de biais quadratique d’omission et de variance d’estimation par les propriétés des projecteurs matriciels, on démontre rigoureusement l’égalité suivante :
Γp = (Biaisp / σ2) + p
où Biaisp correspond à la quantité ∑i=1n (E[y_chapeaui,p] – μi)2. Lorsque le modèle réduit ne souffre d’aucun biais de spécification — c’est-à-dire que tous les régresseurs omis dans le passage de P à p possédaient des coefficients théoriques rigoureusement égaux à zéro —, la composante Biaisp s’annule identiquement. Dans ces conditions de spécification adéquate, Γp se réduit strictement à sa composante de variance, soit exactement la valeur p.
La statistique Cp a été précisément configurée par Colin Mallows pour constituer un estimateur non biaisé de cette quantité théorique Γp. En prenant l’espérance mathématique conditionnelle de Cp sous l’hypothèse que s2 est un estimateur cohérent de σ2 indépendant des résidus du sous-modèle, on obtient :
E[Cp] = E[(RSSp / σ2) – n + 2p] = ((n – p + Biaisp / σ2)) – n + 2p = p + (Biaisp / σ2)
Cette dérivation formelle démontre sans équivoque que si le sous-modèle souffre d’un biais d’omission, E[Cp] excédera systématiquement p d’une quantité strictement proportionnelle à l’ampleur du biais quadratique standardisé. Inversement, si le modèle n’introduit aucun biais d’omission, l’espérance s’établit exactement à p. Toute déviation empirique observable où Cp est inférieur à p s’explique par la nature aléatoire de l’échantillonnage, la variable aléatoire RSSp suivant une distribution du chi-deux non centrée.
3. Le compromis biais-variance et l’interprétation des paramètres
3.1 Dilemme classique entre complexité et généralisation
Le développement de modèles statistiques explicatifs ou prédictifs en psychologie expérimentale et dans les sciences sociales s’inscrit au cœur du dilemme classique entre biais et variance. Un modèle sur-paramétré présente une très faible erreur d’ajustement sur l’échantillon initial parce qu’il possède une flexibilité excessive qui lui permet d’épouser non seulement les tendances structurelles profondes des variables, mais également les fluctuations aléatoires idiosyncrasiques propres au recueil de données. Cette perméabilité au bruit statistique génère une variance d’échantillonnage élevée : si l’on répétait l’étude sur une autre cohorte de sujets, les coefficients estimés subiraient d’importantes oscillations.
À l’opposé du spectre, un modèle чрезмерement parcimonieux compresse l’espace explicatif en éliminant des covariables pertinentes. Si cette simplification drastique a pour effet de stabiliser les coefficients restants — conférant ainsi une faible variance d’échantillonnage aux estimateurs —, elle injecte un biais systématique non nul dans toutes les inférences subséquentes. Les prédictions émises par un tel modèle seront invariablement décalées par rapport aux véritables processus générateurs de données.
Le critère Cp de Mallows opérationnalise précisément ce compromis en formalisant le point d’inflexion où le coût marginal lié à l’augmentation de la variance (pénalisé par le terme + 2p) compense exactement le gain marginal obtenu en termes de réduction de la somme des carrés des résidus (RSSp). Il protège l’investigateur contre les mirages de l’optimisation intra-échantillon en sélectionnant les architectures de régression qui maximisent les capacités de réplication externe et de généralisation hors-échantillon.
3.2 Convention de notation : p versus p + 1
Une source récurrente de confusion et d’erreurs méthodologiques dans la littérature empirique et les manuels d’économétrie réside dans la dénomination formelle du nombre de paramètres, spécifiquement en ce qui concerne le statut de la constante d’ajustement (le terme d’ordonnée à l’origine ou intercept). Selon les disciplines et les traditions logicielles, le symbole p peut désigner soit le nombre total de prédicteurs stochastiques (hors constante), soit le nombre absolu de paramètres libres de l’équation de régression.
Dans la formulation théorique originale de Mallows ainsi que dans la majorité des progiciels statistiques sous R, p représente l’ensemble des paramètres estimés, ce qui inclut systématiquement la constante si celle-ci figure dans l’équation. Ainsi, pour une régression linéaire standard modélisant une variable réponse continue à partir de 4 variables indépendantes et d’une constante, le nombre de paramètres s’élève rigoureusement à :
p = 4 régresseurs + 1 intercept = 5 paramètres
Cette précision est fondamentale pour l’application des règles de décision. Si un chercheur confronte une valeur de Cp calculée à 4 au lieu de 5 en omettant de comptabiliser l’ordonnée à l’origine, il conclura à tort à l’existence d’une inadéquation ou d’un biais systématique. Les packages R modernes, tels que olsrr et leaps, adoptent heureusement cette convention intégrée : la ligne d’équivalence sur laquelle le Cp théorique idéal s’aligne correspond toujours au nombre total de colonnes linéairement indépendantes engagées dans la sous-matrice de projection Xp.
3.3 Typologie des scénarios de sélection
L’analyse des valeurs de Cp issues d’une procédure d’évaluation de modèles candidats permet de catégoriser les sous-ensembles de variables en plusieurs scénarios distincts, dictant la conduite diagnostique du statisticien :
- Le scénario d’adéquation optimale (Cp ≈ p avec p minimal) : Le modèle réduit présente une statistique Cp extrêmement proche de son cardinal de paramètres (par exemple un Cp de 4.12 pour p = 4). Cela signifie que le biais d’omission est virtuellement nul tout en ayant épuré l’équation de régression des dimensions parasitaires. Ce modèle constitue la cible privilégiée du chercheur.
- Le scénario de sous-paramétrisation sévère (Cp >> p) : Le modèle testé affiche une valeur Cp disproportionnellement élevée par rapport à sa dimension (par exemple un Cp de 38.6 pour p = 3). Un tel fossé numérique atteste qu’une ou plusieurs variables déterminantes ont été exclues de la spécification. Les estimations sont grevées d’un biais structurel substantiel qui rend le modèle inapte à la prédiction rigoureuse.
- Le scénario de sur-paramétrisation redondante (Cp ≈ p avec p élevé) : Lorsque l’on examine des modèles dont la taille approche celle du modèle complet P, les valeurs de Cp convergent nécessairement vers P (puisque pour le modèle complet, Cp est par définition égal à P). Cependant, conserver ces variables surnuméraires n’apporte aucune réduction du biais et dégrade la variance des estimations.
- Le scénario sous-unitaire ou d’échantillonnage favorable (Cp < p) : Dans certains cas empiriques, l’analyste observe un Cp inférieur à la dimension p (par exemple un Cp de 2.8 pour p = 4). Loin d’indiquer une anomalie, ce résultat traduit le fait que le sous-modèle bénéficie d’une variance d’échantillonnage résiduelle légèrement inférieure à la prédiction asymptotique moyenne sous distribution gaussienne.
Face à plusieurs spécifications présentant simultanément un écart |Cp – p| minimal, le chercheur devra arbitrer non seulement sur la base du critère numérique, mais également à la lumière de la cohérence théorique, de la faisabilité logistique de collecte des variables et de la parcimonie conceptuelle.
4. Environnement informatique R et préparation des données psychométriques
4.1 Installation et configuration des bibliothèques nécessaires
La mise en œuvre computationnelle du critère Cp de Mallows sous l’environnement CRAN (Comprehensive R Archive Network) bénéficie d’un riche écosystème logiciel. Afin de garantir un environnement de travail reproductible et fluide, plusieurs bibliothèques spécialisées doivent être déployées. Le package olsrr (Ordinary Least Squares Regression Reports) s’impose comme la référence contemporaine pour l’analyse pédagogique et visuelle des régressions linéaires ordinaires. Parallèlement, le package leaps, conçu originellement par Thomas Lumley, offre des routines algorithmiques compilées en Fortran permettant l’exploration exhaustive ultra-rapide des espaces combinatoires de prédicteurs.
Pour orchestrer la manipulation des structures de données, la production de sorties graphiques de haute qualité éditoriale et l’évaluation psychométrique préalable des échelles de mesure, nous mobiliserons conjointement l’écosystème tidyverse (notamment dplyr et ggplot2), le package car (Companion to Applied Regression) pour le calcul des diagnostics de colinéarité, ainsi que le package psych dédié à l’analyse descriptive des construits psychologiques. L’instruction d’installation et de chargement coordonné de ces extensions s’articule via le script R suivant :
install.packages(c("olsrr", "leaps", "car", "tidyverse", "psych"))
library(olsrr)
library(leaps)
library(car)
library(tidyverse)
library(psych)
Il est recommandé de s’assurer que la version active de R soit égale ou supérieure à la version 4.1.0, garantissant l’accès aux optimisations natives de l’opérateur de tuyau (pipe) et une gestion mémoire optimisée lors de la manipulation de matrices de régression de taille intermédiaire à grande.
4.2 Structure et prétraitement du jeu de données
En psychométrie et en sciences sociales appliquées, les données empiriques proviennent typiquement d’échelles d’auto-évaluation de type Likert agrégées en scores composites continus, combinées à des mesures biométriques, comportementales ou sociodémographiques. Considérons, à titre d’illustration opérationnelle tout au long de ce guide, une cohorte de recherche en psychologie de la santé évaluant le niveau d’épuisement professionnel (burnout) auprès d’un échantillon d’individus suivis en milieu organisationnel.
La variable d’intérêt continue à prédire est le score global d’épuisement émotionnel (mesuré par l’inventaire MBI). Les covariables candidates à l’inclusion dans le système prédictif comprennent : le stress perçu (PSS-10), le névrosisme (facette du Big Five), la charge horaire hebdomadaire moyenne, l’indice de soutien social perçu au travail, l’auto-efficacité générale, la qualité subjective du sommeil (indice PSQI), et l’ancienneté professionnelle en années. Avant toute entreprise de modélisation mathématique, une inspection rigoureuse de l’intégrité des vecteurs de données est impérative.
Le prétraitement doit adresser de manière systématique les valeurs manquantes (missing values). Dans le contexte spécifique de la sélection de variables via le critère Cp de Mallows, une règle méthodologique fondamentale s’impose : tous les sous-modèles candidats doivent impérativement être estimés sur un échantillon d’observations rigoureusement identique. Si des lignes sont écartées dynamiquement en raison de valeurs manquantes dispersées à travers les différents régresseurs, la variance de référence s2 et les sommes des carrés des résidus RSSp deviendront mathématiquement incommensurables, invalidant totalement la procédure de comparaison. Une imputation prudente ou l’exclusion cas-par-cas globale (listwise deletion) préalable est donc requise.
4.3 Définition formelle de l’échantillon d’analyse
Pour assurer la solidité inférentielle de l’exercice, nous générons synthétiquement un jeu de données psychométriques réaliste comportant 250 participants (n = 250), dimensionné de façon adéquate pour respecter un ratio sujets/variables d’au moins 25 observations par prédicteur, minimisant ainsi les risques de sur-ajustement accidentel dès le stade de la conception. Le script R ci-dessous configure la matrice de données avec une structure de covariance empirique maîtrisée :
set.seed(1973) # Hommage a l'annee de publication de Colin Mallows
n <- 250
stress <- rnorm(n, mean = 25, sd = 6)
nevrosisme <- rnorm(n, mean = 50, sd = 10) + 0.4 * stress
charge_travail <- rnorm(n, mean = 40, sd = 8)
soutien_social <- rnorm(n, mean = 60, sd = 12) - 0.3 * stress
auto_efficacite <- rnorm(n, mean = 30, sd = 5) - 0.2 * nevrosisme
sommeil_degrade <- rnorm(n, mean = 8, sd = 3) + 0.3 * stress
anciennete <- runif(n, min = 1, max = 25)
Dans ce modèle de simulation générative, concevons la variable réponse epuisement de telle sorte que seules certaines variables exercent un effet structural direct réel, les autres n’introduisant que des corrélations collatérales ou du bruit stochastique :
epuisement <- 15 + 0.8 * stress + 0.5 * nevrosisme - 0.4 * soutien_social + 1.2 * sommeil_degrade + rnorm(n, mean = 0, sd = 8)
df_psycho <- data.frame(epuisement, stress, nevrosisme, charge_travail, soutien_social, auto_efficacite, sommeil_degrade, anciennete)
Cette configuration matricielle servira de substrat expérimental transparent : nous savons a priori que la charge de travail, l’auto-efficacité et l’ancienneté sont structurellement redondantes dans l’équation de génération une fois contrôlés le stress, le névrosisme, le soutien social et le sommeil. L’objectif analytique consistera à vérifier si le critère Cp de Mallows sous R isole avec précision ce sous-ensemble théorique exact.
5. Ajustement du modèle complet et des sous-modèles via la fonction lm()
5.1 Spécification rigoureuse du modèle saturé
L’estimation sous R débute obligatoirement par la spécification formelle du modèle linéaire complet au moyen de la primitive fondamentale lm(). Ce modèle complet, parfois qualifié de modèle saturé ou de référence, doit incorporer simultanément l’ensemble des régresseurs identifiés comme scientifiquement plausibles. La syntaxe compacte de R permet de déclarer la régression de la variable réponse sur l’intégralité des colonnes résiduelles de la matrice de données grâce à la notation par point :
modele_complet <- lm(epuisement ~ ., data = df_psycho)
L’ajustement de cet objet de classe lm constitue la clé de voûte de toute l’évaluation subséquente. C’est à partir de cet objet que l’estimateur sans biais de la variance de l’erreur résiduelle sera extrait. Dans l’environnement R, cette variance correspond au carré de la statistique désignée sous le terme de residual standard error présente dans la sortie textuelle standard. Algébriquement, cette valeur est accessible directement via la structure interne de la liste R :
sigma2_complet <- summary(modele_complet)$sigma^2
Il est capital d’inspecter visuellement et numériquement ce modèle initial. Si le modèle complet souffre de singularités mathématiques — par exemple si le rang de la matrice X est inférieur au nombre de colonnes en raison d’une colinéarité parfaite —, la fonction lm() produira des coefficients non définis (marqués par la valeur NA). Une telle situation invaliderait immédiatement l’estimation de σ2 et doit être résolue avant de procéder au calcul du critère de Mallows.
5.2 Construction manuelle de modèles concurrents emboîtés et non emboîtés
Pour illustrer la mécanique comparée, procédons à l’ajustement explicite de plusieurs modèles concurrents représentant des hypothèses théoriques alternatives quant aux déterminants de l’épuisement professionnel. Cette approche manuelle reflète fidèlement la pratique de recherche confirmatoire où le statisticien teste un ensemble restreint de spécifications dictées par la littérature scientifique :
- Modèle 1 (Modèle purement biologique et contextuel) : Épuisement prédit par la charge de travail, l’ancienneté et le sommeil dégradé.
m1 <- lm(epuisement ~ charge_travail + anciennete + sommeil_degrade, data = df_psycho) - Modèle 2 (Modèle cognitivo-émotionnel) : Épuisement prédit par le stress perçu, le névrosisme et l’auto-efficacité.
m2 <- lm(epuisement ~ stress + nevrosisme + auto_efficacite, data = df_psycho) - Modèle 3 (Modèle structurel congruent) : Épuisement prédit par le stress perçu, le névrosisme, le soutien social et le sommeil dégradé.
m3 <- lm(epuisement ~ stress + nevrosisme + soutien_social + sommeil_degrade, data = df_psycho) - Modèle 4 (Modèle quasi-complet avec bruit) : Épuisement prédit par toutes les variables à l’exclusion unique de l’ancienneté.
m4 <- lm(epuisement ~ stress + nevrosisme + charge_travail + soutien_social + auto_efficacite + sommeil_degrade, data = df_psycho)
Chacun de ces objets encapsule une structure matricielle spécifique, caractérisée par une somme des carrés résiduelle distincte (accessible par la commande deviance(m)) et un nombre de degrés de liberté résiduels propre. Le tableau d’analyse comparée consistera à évaluer comment chacune de ces configurations théoriques se positionne par rapport au modèle complet de référence.
5.3 Contrôle des sorties intermédiaires
Avant d’engager le calcul formel du Cp, l’investigateur doit soumettre chaque modèle candidat à un examen rigoureux via la fonction summary(). Cet examen intermédiaire ne vise pas à sélectionner les variables sur la base exclusive de leur significativité marginale (valeur de p individuelle du test de Student), car cette pratique courante souffre notoirement du piège de l’omission conjointe : deux prédicteurs individuellement non significatifs en raison d’une corrélation mutuelle modérée peuvent s’avérer conjointement indispensables pour réduire la variance d’erreur.
L’inspection des sorties intermédiaires permet de relever les coefficients de détermination classiques R2 et ajustés R2adj. On observe systématiquement que le R2 ordinaire progresse de façon monotone à mesure que l’on élargit la spécification paramétrique : le Modèle 4 affichera invariablement un R2 supérieur au Modèle 3, qui lui-même surpassera les modèles 1 et 2. Cependant, cette progression brute est incapable d’indiquer si le gain explicatif compense la perte de parcimonie.
Le contrôle dimensionnel s’opère en extrayant le nombre effectif de paramètres pour chaque modèle. Sous R, la dimension p d’un modèle estimé par lm() s’obtient avec exactitude via la longueur du vecteur des coefficients : length(coef(m)). Il convient de s’assurer scrupuleusement que ce décompte inclut bien le terme (Intercept), conformément aux fondements mathématiques du critère explicités à la section 3.2.
6. Calcul du Cp de Mallows avec le package olsrr et la fonction ols_mallows_cp()
6.1 Mise en œuvre technique de ols_mallows_cp()
Le package olsrr a standardisé l’extraction du critère de Mallows au travers d’une fonction dédiée, au nom explicite et à la signature syntaxique intuitive : ols_mallows_cp(). Cette fonction matérialise l’interaction computationnelle directe entre un modèle candidat restreint et le modèle saturé de référence. Sa signature logicielle requiert deux arguments positionnels obligatoires : l’objet lm représentant le sous-modèle analysé et l’objet lm incarnant le modèle de référence complet.
L’appel canonique s’exécute selon le schéma syntaxique suivant :
ols_mallows_cp(model, full_model)
Lors de l’évaluation de cette expression, R exécute en arrière-plan une séquence d’opérations vectorielles hautement optimisées : il extrait la déviance résiduelle du sous-modèle, calcule la variance résiduelle non biaisée du modèle saturé, quantifie le nombre d’observations effectives (n) ainsi que la dimension paramétrique (p) de l’équation réduite, puis applique rigoureusement l’algorithme :
Cp = (deviance(model) / (deviance(full_model) / df.residual(full_model))) – nobs(model) + 2 * length(coef(model))
La sortie renvoyée par la fonction se présente sous la forme d’une valeur scalaire numérique continue à double précision. Cette simplicité d’utilisation autorise une intégration aisée au sein de boucles de traitement, de fonctions personnalisées de la famille apply ou de flux de calcul construits avec le package purrr.
6.2 Exemple pas à pas basé sur un jeu de données de référence
Appliquons immédiatement cette fonction opérationnelle à notre ensemble de modèles psychométriques concurrents initialisés à la section 5.2, en utilisant l’objet modele_complet comme étalon de référence pour la variance de l’erreur :
cp_m1 <- ols_mallows_cp(m1, modele_complet)
cp_m2 <- ols_mallows_cp(m2, modele_complet)
cp_m3 <- ols_mallows_cp(m3, modele_complet)
cp_m4 <- ols_mallows_cp(m4, modele_complet)
Pour structurer l’analyse comparative et faciliter l’interprétation substantielle des résultats, synthétisons ces sorties numériques au sein d’un tableau récapitulatif structuré (tibble ou data.frame) intégrant la dimension p correspondante ainsi que l’écart absolu |Cp – p| :
resultats_cp <- data.frame(
Modele = c("M1: Bio-contexte", "M2: Cognitivo-emotif", "M3: Structurel congruent", "M4: Quasi-complet"),
Parametres_p = c(length(coef(m1)), length(coef(m2)), length(coef(m3)), length(coef(m4))),
Mallows_Cp = c(cp_m1, cp_m2, cp_m3, cp_m4)
)
resultats_cp$Ecart_Absolu <- \abs(resultats_cp$Mallows_Cp - resultats_cp$Parametres_p)
print(resultats_cp)
L’exécution de cette séquence met en relief la puissance discriminative du critère. Le Modèle 1 (p = 4) produit une valeur de Cp massivement élevée (dépassant typiquement le seuil de 80), traduisant un biais d’omission sévère consécutif à l’exclusion du stress et du névrosisme. Le Modèle 2 (p = 4), bien que conceptuellement séduisant, exhibe également un Cp nettement supérieur à 4 en raison de l’absence du soutien social et du sommeil. Le Modèle 4 (p = 7) affiche quant à lui un Cp très proche de sa dimension (approchant 7.0), mais au prix d’une sur-paramétrisation inutile. Enfin, le Modèle 3 (p = 5) se singularise par une valeur de Cp gravitant rigoureusement autour de 4.8 à 5.1, minimisant simultanément le critère global et l’écart |Cp – p|, confirmant empiriquement la structure génératrice sous-jacente.
6.3 Erreurs classiques de manipulation et résolution
Lors de l’exploitation de la fonction ols_mallows_cp() en contexte appliqué, plusieurs erreurs de manipulation fréquentes sont susceptibles d’interrompre l’exécution du script ou de corrompre silencieusement la validité des inférences :
- L’inversion involontaire des arguments fonctionnels : Une méprise classique consiste à spécifier le modèle complet en premier argument et le sous-modèle en second argument (
ols_mallows_cp(modele_complet, m3)). Dans ce cas, R calcule un Cp mathématiquement absurde où le modèle complet est pénalisé par rapport à la variance d’un modèle sous-spécifié. L’analyste doit mémoriser la séquence impérative : premier argument = modèle réduit à tester ; second argument = modèle complet saturé. - La divergence du nombre d’observations (taille d’échantillon inégale) : Si le jeu de données initial contient des valeurs manquantes non traitées et que le modèle complet intègre une variable présentant des données manquantes que le sous-modèle ne contient pas, la fonction
lm()appliquera des suppressions d’observations asymétriques. R renverra alors un message d’erreur explicite signalant une non-concordance de dimensionnalité vectorielle : « models were not fitted to the same size of dataset ». La solution consiste à créer un sous-ensemble préalablement filtré parna.omit()sur l’ensemble des variables du modèle saturé. - L’omission de termes d’interaction dans le modèle de référence : Si le chercheur évalue un sous-modèle incorporant une interaction multiplicative (ex.
stress * soutien_social), mais que cette interaction a été omise de la déclaration du modèle complet, l’hypothèse de non-biais du modèle saturé s’effondre. Le modèle de référence doit obligatoirement inclure l’espace vectoriel le plus vaste envisagé dans l’étude.
7. Approches alternatives sous R : calcul manuel et packages complémentaires
7.1 Implémentation d’une fonction R personnalisée pour le calcul de Cp
Bien que les bibliothèques spécialisées offrent un confort opératoire appréciable, la maîtrise conceptuelle approfondie d’un critère statistique s’enrichit considérablement de sa programmation sous forme d’algorithme matriciel ou vectoriel autonome (from scratch). L’implémentation d’une fonction personnalisée sous R permet de s’affranchir des dépendances logicielles externes, d’accélérer les temps de calcul lors de simulations Monte Carlo intensives et d’intégrer des routines de contrôle d’intégrité sur mesure.
Considérons l’écriture d’une fonction robuste baptisée custom_mallows_cp(), conçue pour extraire les scalaires analytiques indispensables directement depuis les environnements d’évaluation des objets lm :
custom_mallows_cp <- function(sub_model, full_model) {
# Verification formelle de la classe des objets en entree
if (!inherits(sub_model, "lm") || !inherits(full_model, "lm")) {
stop("Les deux arguments doivent obligatoirement etre des objets de classe 'lm'.")
}
# Controle de l'identite rigoureuse du nombre d'observations
n_sub <- stats::nobs(sub_model)
n_full <- stats::nobs(full_model)
if (n_sub != n_full) {
stop(paste("Incompatibilite de taille d'echantillon : modele reduit n =", n_sub,
"| modele complet n =", n_full))
}
# Extraction de la somme des carres des residus (RSS) et des degres de liberte
rss_p <- stats::deviance(sub_model)
rss_full <- stats::deviance(full_model)
df_full <- stats::df.residual(full_model)
# Estimation non biaisee de la variance de l'erreur residuelle du modele complet
sigma2_est <- rss_full / df_full
# Nombre total de parametres du sous-modele (incluant l'intercept)
p <- length(stats::coef(sub_model))
# Application formelle de la formule algebrique de Mallows
cp_val <- (rss_p / sigma2_est) - n_sub + (2 * p)
return(cp_val)
}
Cette implémentation native garantit une transparence totale sur le flux de calcul. La vérification croisée entre custom_mallows_cp(m3, modele_complet) et ols_mallows_cp(m3, modele_complet) produit des résultats rigoureusement identiques au niveau de précision de la virgule flottante, attestant de la concordance théorique absolue du processus.
7.2 Utilisation de la bibliothèque leaps pour la régression exhaustive
Lorsque le nombre de variables prédictives potentielles excède quatre ou cinq régresseurs, l’ajustement manuel individuel de toutes les combinaisons possibles devient rapidement irréalisable sur le plan combinatoire. Pour un ensemble de K prédicteurs, le nombre total d’équations candidates envisageables s’élève à 2K (soit 128 modèles pour K = 7, mais 1 048 576 pour K = 20). Pour relever ce défi computationnel, la bibliothèque leaps implémente l’algorithme hautement efficace de Furnival et Wilson (1974), reposant sur une exploration arborescente par séparation et évaluation progressive (branch-and-bound).
La fonction maîtresse regsubsets() permet d’exécuter cette recherche combinatoire intégrale en une seule ligne de code, en isolant instantanément la meilleure combinaison de prédicteurs pour chaque taille de sous-modèle donnée :
selection_exhaustive <- regsubsets(epuisement ~ ., data = df_psycho, nvmax = 7, method = "exhaustive")
resume_leaps <- summary(selection_exhaustive)
L’objet de synthèse produit par summary(selection_exhaustive) recèle une matrice logique détaillant les prédicteurs engagés, ainsi qu’un vecteur numérique contenant les valeurs de Cp associées à la meilleure combinaison identifiée pour chaque strate de dimensionnalité :
resume_leaps$cp
Pour convertir ces sorties vectorielles brutes en un tableau analytique contemporain et manipulable, nous pouvons structurer les résultats à l’aide des fonctions de manipulation de données modernes :
synthese_leaps <- tibble(
Nombre_Predicteurs = 1:length(resume_leaps$cp),
Parametres_p = (1:length(resume_leaps$cp)) + 1, # Prise en compte de la constante
Cp_Mallows = resume_leaps$cp,
R2 = resume_leaps$rsq,
R2_Ajuste = resume_leaps$adjr2,
BIC = resume_leaps$bic
)
print(synthese_leaps)
Cette approche systématique met en évidence la trajectoire d’évolution du Cp à travers l’espace des dimensions, identifiant instantanément à quel niveau de cardinalité le critère plonge pour la première fois à proximité immédiate de la ligne d’équipartition p.
7.3 Sélection multimorphique avec MuMIn et car
L’écosystème logiciel de R dispose d’autres extensions sophistiquées dédiées à l’inférence multimorphique et à la sélection d’équations structurales. Parmi celles-ci, le package MuMIn (Multi-Model Inference) fournit un cadre de travail avancé pour le classement, le partitionnement et la combinaison bayésienne ou fréquentiste de modèles de régression. Bien que MuMIn soit principalement orienté vers le calcul de l’AICc (critère d’Akaike corrigé pour petits échantillons), il permet d’intégrer des critères de pénalisation personnalisés au sein de ses tables d’ordonnancement global.
La fonction dredge() de MuMIn permet de générer automatiquement la totalité des sous-modèles dérivés d’un modèle global préalablement ajusté. Elle requiert toutefois que l’argument d’ajustement global soit configuré de façon à préserver l’accès intègre aux données d’origine (via l’instruction options(na.action = "na.fail")) :
library(MuMIn)
options(na.action = "na.fail")
modele_global <- lm(epuisement ~ ., data = df_psycho)
table_selection <- dredge(modele_global)
head(table_selection, 5)
Parallèlement, la bibliothèque car apporte un support méthodologique déterminant en permettant de tester si l’omission conjointe d’un ensemble de prédicteurs écartés par l’arbitrage du Cp induit une dégradation de variance statistiquement significative. L’emploi conjoint de la fonction linearHypothesis() ou de l’analyse de variance comparative anova(modele_reduit, modele_complet) permet ainsi d’établir un pont diagnostique formel entre la règle décisionnelle heuristique de Colin Mallows et les tests de rapport de variance de Fisher-Snedecor traditionnels.
8. Procédure de sélection automatisée par sous-ensembles et comparaison graphique
8.1 Exécution de la sélection globale avec ols_step_best_subset()
La bibliothèque olsrr propose une interface d’évaluation automatisée de niveau supérieur grâce à la fonction ols_step_best_subset(). Contrairement aux approches séquentielles pas à pas (stepwise regression) qui souffrent de limitations analytiques sévères en présence de colinéarité, cette routine examine rigoureusement l’ensemble des 2K configurations possibles et synthétise les métriques d’évaluation décisionnelles au sein d’une structure de données unifiée.
L’exécution de la procédure s’opère directement à partir du modèle complet initial :
selection_sous_ensembles <- ols_step_best_subset(modele_complet)
print(selection_sous_ensembles)
La console affiche un rapport tabulaire d’une richesse exceptionnelle. Pour chaque niveau de dimension paramétrique (de 1 à K prédicteurs), la routine consigne la combinaison gagnante, son coefficient de détermination R2, son R2 ajusté, le critère d’information d’Akaike (AIC), le critère de Schwarz-Bayes (BIC), le critère de Sawa, et de façon centrale, le Cp de Mallows. L’analyste peut ainsi retracer visuellement le cheminement de la pénalisation : à mesure que le modèle intègre les variables structurellement déterminantes (le stress, puis le névrosisme, puis le soutien social, puis le sommeil), le Cp s’effondre de valeurs stratosphériques (supérieures à 200 pour le modèle à une variable) jusqu’à atteindre un point d’équilibre stabilisé.
8.2 Visualisation graphique : le diagramme de Mallows
L’apport méthodologique fondamental de Colin Mallows en 1973 ne résidait pas uniquement dans la dérivation algébrique de sa formule, mais dans la proposition d’un dispositif d’exploration visuelle désormais célèbre : le diagramme de Cp (ou Cp plot). Ce graphique cartographie sur l’axe des ordonnées la valeur numérique de la statistique Cp calculée pour chaque modèle candidat, et sur l’axe des abscisses le nombre de paramètres totaux p engagés dans l’équation.
La fonction générique intégrée à olsrr génère automatiquement ce type de rendu :
plot(selection_sous_ensembles)
Cependant, afin de répondre aux exigences de clarté et de personnalisation des publications scientifiques à fort facteur d’impact, il est hautement préférable de concevoir ce tracé directement sous l’environnement ggplot2. Ce contrôle programmatique permet d’incorporer explicitement la ligne théorique d’équipartition y = x symbolisant l’absence de biais :
df_graph <- as.data.frame(selection_sous_ensembles$metrics)
df_graph$p <- df_graph$n + 1 # n representant ici le nombre de predicteurs hors intercept
ggplot(df_graph, aes(x = p, y = mallows_cp)) +
geom_point(color = "#1f77b4", size = 3.5, alpha = 0.8) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "#d62728", linewidth = 0.9) +
geom_text(aes(label = paste0("p=", p, "nCp=", round(mallows_cp, 1))),
vjust = -0.7, size = 3, check_overlap = TRUE) +
scale_x_continuous(breaks = 2:(ncol(df_psycho))) +
labs(
title = "Diagramme de selection de Mallows (Cp Plot)",
subtitle = "La droite en pointille rouge figure la reference d'absence theorique de biais (y = p)",
x = "Nombre total de parametres estimes p (constante incluse)",
y = "Statistique Cp de Mallows"
) +
theme_minimal(base_size = 12) +
theme(
panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5)
)
8.3 Interprétation sémiologique du graphique de sélection
L’interprétation sémiologique de cette cartographie visuelle repose sur la position relative des points empiriques par rapport à la droite diagonale de référence y = x :
- Zone supérieure gauche (Cp >> p) : Les modèles positionnés haut au-dessus de la ligne diagonale souffrent d’un biais d’omission sévère. Bien que parcimonieux en termes de paramètres, leur erreur quadratique totale de prédiction est inacceptable. Ces modèles doivent être catégoriquement rejetés.
- Trajectoire d’inflexion (Le point de coude) : En examinant le tracé de gauche à droite, l’analyste observe généralement une chute abrupte de la courbe du Cp jusqu’à atteindre pour la première fois la proximité immédiate de la droite y = x. Ce point d’inflexion signale la saturation de l’information explicative majeure.
- Zone de concordance optimale (Cp ≈ p) : Le premier modèle situé sur la ligne diagonale ou immédiatement en dessous de celle-ci avec le plus petit nombre de paramètres identifie la spécification idéale. Dans notre simulation, le point correspondant à p = 5 (stress, nevrosisme, soutien_social, sommeil_degrade + intercept) s’ancre précisément sur la droite.
- Plateau asymptotique à droite (Cp ≈ p pour p croissant) : À mesure que p progresse au-delà du point optimal vers le modèle saturé P, les points continuent de coller à la diagonale. Toutefois, cette trajectoire n’apporte aucune réduction substantielle du biais : elle ne fait qu’ajouter des paramètres redondants dont le coût en variance compense stérilement l’infime gain résiduel.
9. Analyse comparative : Cp de Mallows versus AIC, BIC et R-carré ajusté
9.1 Propriétés asymptotiques face à l’Akaike Information Criterion (AIC)
Le critère Cp de Mallows entretient des affinités formelles profondes avec le célèbre critère d’information d’Akaike (AIC). Sous l’hypothèse de résidus gaussiens distribués selon une loi normale d’espérance nulle et de variance constante σ2, il est mathématiquement démontrable que le critère Cp de Mallows et l’AIC sont asymptotiquement équivalents. En effet, l’expression de la log-vraisemblance d’un modèle de régression linéaire conduit à la formulation suivante de l’AIC :
AIC = n × ln(RSSp / n) + 2p + Cte
En effectuant un développement limité au premier ordre de la fonction logarithmique au voisinage de la variance de référence s2, la pénalité de l’AIC converge vers la fonction de pénalisation linéaire 2p caractéristique du Cp de Mallows. Dans les deux cas, le coût marginal alloué à l’introduction d’un paramètre supplémentaire correspond exactement à deux unités de variance résiduelle normalisée.
Néanmoins, des nuances opérationnelles notables subsistent à petite échelle d’échantillonnage. Le Cp de Mallows utilise explicitement l’estimateur sans biais de la variance résiduelle dérivé du modèle complet de référence au dénominateur, tandis que l’AIC classique évalue la log-vraisemblance interne de chaque sous-modèle de manière autonome (en estimant sa variance par le maximum de vraisemblance RSSp / n, qui constitue un estimateur biaisé pour les faibles valeurs de n). Par conséquent, sur des échantillons restreints, le Cp de Mallows offre fréquemment une évaluation plus stable et moins sensible aux fluctuations locales d’ajustement que l’AIC non corrigé.
9.2 Cp de Mallows face au Bayesian Information Criterion (BIC)
La confrontation entre le Cp de Mallows et le critère d’information bayésien (BIC ou critère de Schwarz) illustre deux philosophies statistiques distinctes de la sélection de modèles. L’équation canonique du BIC s’énonce selon :
BIC = n × ln(RSSp / n) + ln(n) × p
La différence fondamentale réside dans la sévérité de la pénalisation paramétrique. Alors que le Cp de Mallows et l’AIC appliquent une pénalité fixe égale à 2 unités par paramètre quel que soit le volume d’observations, le BIC applique un coefficient multiplicateur logarithmique ln(n). Dès lors que la taille d’échantillon n excède e2 ≈ 7.39 observations (ce qui est systématiquement le cas en pratique), ln(n) devient strictement supérieur à 2. Le BIC pénalise donc la complexité structurelle avec une rigueur exponentiellement croissante par rapport au Cp.
Sur le plan des propriétés fondamentales, le BIC est asymptotiquement convergent (consistent) : si le « vrai » modèle générateur de données appartient formellement à la famille des sous-modèles examinés, la probabilité que le BIC sélectionne ce modèle exact tend vers 1 lorsque n tend vers l’infini. En revanche, le Cp de Mallows est efficace pour la prédiction (asymptotically efficient) : il minimise l’erreur quadratique moyenne de prédiction sur de nouvelles observations, quitte à conserver une légère propension à sur-sélectionner des prédicteurs marginaux sur les très grands échantillons. En psychologie et en épidémiologie, où les phénomènes comportementaux sont multifactoriels et où la notion de « vrai » modèle parcimonieux est une abstraction théorique, l’orientation prédictive du Cp s’avère souvent hautement pragmatique.
9.3 Limites critiques du R-carré ajusté
Le coefficient de détermination ajusté (R2adj ou adjusted R-squared), proposé initialement par Henri Theil, constitue traditionnellement la métrique la plus familière des étudiants et des praticiens de la régression linéaire. Conçu pour pallier l’accroissement mécanique du R2 standard, il pondère le rapport des sommes des carrés par leurs degrés de liberté respectifs :
R2adj = 1 – [ (RSSp / (n – p)) / (TSS / (n – 1)) ]
Malgré sa popularité pédagogique, le R2 ajusté présente des limites critiques structurelles qui le rendent très inférieur au Cp de Mallows pour la sélection rigoureuse de modèles :
- Pénalisation marginale insuffisante : On démontre analytiquement que le R2 ajusté augmente dès lors qu’un prédicteur additionnel présente une statistique de Student en valeur absolue strictement supérieure à 1 (|t| > 1, soit une valeur de p approximativement inférieure à 0.32). Cette condition d’inclusion est particulièrement laxiste et conduit fréquemment le chercheur à retenir des prédicteurs purement parasitaires qui détériorent la variance de prédiction.
- Absence de quantification explicite du biais : Le R2 ajusté normalise la variance intra-modèle sans jamais la confronter à un étalon de référence extérieur impartial. Il est incapable d’indiquer si le sous-modèle retenu génère un biais systématique d’omission.
- Divergences fréquentes en présence de bruit : Dans des contextes de colinéarité modérée, le R2 ajusté tend à sélectionner le modèle contenant le plus grand nombre de variables, là où le diagramme de Mallows met en exergue la stabilisation du Cp bien en amont de cette saturation artificielle.
10. Cas pratique approfondi en recherche psychologique et comportementale
10.1 Présentation du cas d’étude et hypothèses scientifiques
Afin d’ancrer ces concepts dans la réalité d’un laboratoire de recherche quantitative, développons une étude de cas détaillée portant sur l’évaluation psychométrique de la détresse psychologique en contexte post-traumatique ou de stress chronique. Une équipe de chercheurs souhaite concevoir un outil de dépistage rapide (screening) permettant de prédire le score d’anxiété clinique généralisée (mesuré par l’échelle continue GAD-7) auprès d’une cohorte de patients admis en consultation de soutien.
Le protocole d’admission initial administre un protocole exhaustif documentant dix prédicteurs psychométriques et physiologiques potentiels :
- Névrosisme (inventaire NEO-PI-R)
- Réactivité émotionnelle (score composite d’intensité affective)
- Soutien familial perçu (échelle MSPSS)
- Soutien amical perçu (échelle MSPSS)
- Stratégies de coping actif (échelle Brief-COPE)
- Stratégies de coping par évitement (échelle Brief-COPE)
- Fréquence cardiaque au repos (biomarqueur continu)
- Durée moyenne de sommeil paradoxal (mesure polysomnographique)
- Niveau de revenu annuel déclaré (variable socio-économique)
- Âge chronologique du participant
L’hypothèse scientifique fondamentale postule que la détresse anxieuse clinique peut être modélisée avec une précision prédictive optimale au moyen d’un sous-ensemble restreint combinant le névrosisme, le coping par évitement et la perturbation du sommeil, sans qu’il soit nécessaire de mobiliser les mesures biométriques coûteuses ou les indicateurs socio-économiques périphériques. L’objectif méthodologique consiste à appliquer la procédure de Mallows pour arbitrer rigoureusement cette hypothèse.
10.2 Pipeline de traitement R complet pas à pas
Déployons l’intégralité du pipeline computationnel sous R, depuis la simulation contrôlée du jeu de données cliniques (n = 300) jusqu’à la sélection automatisée guidée par le Cp :
# 1. Generation de la matrice clinique simulee
set.seed(2024)
n_clinique <- 300
nevrosisme <- rnorm(n_clinique, mean = 55, sd = 12)
reactivite <- rnorm(n_clinique, mean = 50, sd = 10) + 0.35 * nevrosisme
soutien_famille <- rnorm(n_clinique, mean = 28, sd = 6)
soutien_amis <- rnorm(n_clinique, mean = 26, sd = 5) + 0.25 * soutien_famille
coping_actif <- rnorm(n_clinique, mean = 32, sd = 7)
coping_evitement <- rnorm(n_clinique, mean = 24, sd = 6) + 0.3 * nevrosisme
freq_cardiaque <- rnorm(n_clinique, mean = 72, sd = 9)
sommeil_rem <- rnorm(n_clinique, mean = 90, sd = 20) - 0.4 * nevrosisme
revenu <- runif(n_clinique, min = 15, max = 80)
age <- rnorm(n_clinique, mean = 42, sd = 11)
# 2. Generation de la variable reponse : Anxiete GAD-7
anxiete_gad7 <- 10 + 0.45 * nevrosisme + 0.65 * coping_evitement - 0.15 * sommeil_rem + rnorm(n_clinique, mean = 0, sd = 6)
# 3. Assemblage du dataframe d'etude
df_clinique <- data.frame(
anxiete_gad7, nevrosisme, reactivite, soutien_famille, soutien_amis,
coping_actif, coping_evitement, freq_cardiaque, sommeil_rem, revenu, age
)
# 4. Ajustement du modele complet sature
mod_complet_clinique <- lm(anxiete_gad7 ~ ., data = df_clinique)
# 5. Selection exhaustive des sous-ensembles optimaux
selection_clinique <- ols_step_best_subset(mod_complet_clinique)
# 6. Extraction et structuration des metriques comparatives
metrics_tbl <- as_tibble(selection_clinique$metrics) %>%
mutate(p_total = n + 1,
ecart_mallows = abs(mallows_cp - p_total)) %>%
select(n, p_total, predictors, rsquare, adjr, mallows_cp, ecart_mallows, aic, bic)
# 7. Isolement de la solution minimisant le compromis de Mallows
modele_selectionne <- metrics_tbl %>%
filter(mallows_cp <= p_total + 1.5) %>%
arrange(p_total) %>%
slice(1)
print(modele_selectionne)
10.3 Synthèse et discussion des résultats psychométriques
L’inspection analytique de la table générée met en évidence des constats psychométriques d’une grande limpidité. Pour n = 1 régresseur (le névrosisme seul, p = 2), la statistique Cp s’établit à une valeur massivement prohibitive (Cp > 120), démontrant sans équivoque que la prise en compte exclusive du trait de personnalité n’épuise nullement la variance clinique de l’anxiété. Pour n = 2 régresseurs (névrosisme et coping par évitement, p = 3), le Cp régresse spectaculairement autour de 22.4, mais demeure substantiellement au-dessus du seuil de non-biais.
La bascule décisive s’opère rigoureusement à n = 3 prédicteurs (p = 4 paramètres), correspondant à l’intégration conjointe de nevrosisme, coping_evitement et sommeil_rem. À cet échelon précis, le critère Cp de Mallows chute à 3.78. Cette valeur est à la fois remarquablement basse et strictement inférieure à la dimension paramétrique nominale (p = 4), attestant d’une absence complète de biais de spécification couplée à une variance minimale.
L’examen des strates ultérieures (de n = 4 à n = 10 prédicteurs) confirme la pertinence de la règle de décision. L’adjonction de la fréquence cardiaque, du revenu ou des facettes de soutien social ne provoque qu’une réduction marginale et stochastique de la somme des carrés résiduelle : le Cp oscille stérilement entre 4.8 et 11.0 en collant à la droite d’équipartition sans jamais améliorer le pouvoir explicatif net. Grâce au critère de Mallows, le protocole clinique peut valablement élaguer 7 des 10 questionnaires initiaux, allégeant drastiquement le fardeau cognitif imposé aux patients sans aucune dégradation de la validité de dépistage.
11. Diagnostics résiduels et validation post-sélection du modèle retenu
11.1 Validation formelle des postulats de régression gaussienne
L’isolement d’un modèle parcimonieux par le critère Cp de Mallows ne dispense aucunement le statisticien de l’obligation méthodologique de soumettre la spécification retenue à une batterie complète de diagnostics résiduels. L’ensemble de la dérivation mathématique du Cp reposant sur l’hypothèse de perturbations résiduelles sphériques et gaussiennes, toute entorse sévère à ces principes fragiliserait les conclusions tirées.
Ajustons de façon autonome le modèle clinique optimisé isolé à la section précédente :
modele_optimal <- lm(anxiete_gad7 ~ nevrosisme + coping_evitement + sommeil_rem, data = df_clinique)
La validation formelle implique l’exécution des tests diagnostiques canoniques :
- Normalité des résidus : Évaluée conjointement par l’inspection du diagramme quantile-quantile (QQ-plot) et par le test formel de Shapiro-Wilk :
shapiro.test(residuals(modele_optimal))
Une valeur de p supérieure au seuil usuel de 0.05 confirme que la distribution des erreurs ne s’écarte pas substantiellement d’une loi normale univariée. - Homoscédasticité des perturbations : Contrôlée par le test d’homogénéité des variances de Breusch-Pagan implémenté dans le package car :
car::ncvTest(modele_optimal)
L’absence de significativité statistique valide l’hypothèse de variance constante de l’erreur σ2 à travers tout le spectre prédictif. - Linéarité structurelle : Vérifiée par la confrontation graphique des résidus ordinaires standardisés aux valeurs ajustées (fitted values), garantissant l’absence de courbures polynomiales non modélisées.
11.2 Évaluation de la colinéarité résiduelle (VIF)
L’un des bénéfices collatéraux majeurs de l’élagage parcimonieux orchestré par le critère de Mallows réside dans l’assainissement de la matrice de covariance des prédicteurs. Dans le modèle complet initial, la présence simultanée de multiples indicateurs psychologiques corrélés engendre fréquemment une instabilité numérique délétère désignée sous le nom de multicolinéarité.
Cette instabilité se quantifie au moyen du Facteur d’Inflation de la Variance (Variance Inflation Factor ou VIF). Pour un prédicteur j donné, le VIF mesure à quel point la variance de l’estimateur β_chapeauj est gonflée par sa corrélation linéaire avec l’ensemble des autres régresseurs :
VIFj = 1 / (1 – Rj2)
Comparons les indices VIF calculés avant et après sélection à l’aide de la primitive vif() :
vif_complet <- car::vif(mod_complet_clinique)
vif_optimal <- car::vif(modele_optimal)
print("VIF du modele complet :")
print(vif_complet)
print("VIF du modele optimal apres arbitrage par Mallows :")
print(vif_optimal)
Alors que le modèle complet présentait des facteurs VIF élevés sur les échelles de coping et de réactivité en raison de redondances conceptuelles manifestes, le modèle restreint affiche des indices VIF strictement inférieurs au seuil critique de 2.0. Cela confirme que le critère Cp a naturellement éliminé les prédicteurs redondants, restituant une base vectorielle quasi-orthogonale garantissant des estimations de coefficients stables et aisément interprétables.
11.3 Validation croisée et stabilité de l’échantillonnage
Bien que le Cp de Mallows constitue un estimateur asymptotique sans biais de l’erreur quadratique de prédiction, son comportement effectif hors-échantillon gagne à être corroboré par une procédure de validation croisée à K blocs (K-fold cross-validation). Cette technique empirique consiste à partitionner l’échantillon global en K sous-groupes indépendants (usuellement K = 10), à entraîner itérativement le modèle sur 9 blocs et à mesurer son erreur quadratique moyenne de prédiction (RMSEP) sur le dixième bloc non vu.
Moyennant l’utilisation des structures vectorielles de base de R, la validation croisée s’implémente en quelques instructions limpides :
set.seed(42)
k <- 10
folds <- sample(rep(1:k, length.out = nrow(df_clinique)))
cv_errors <- numeric(k)
for (i in 1:k) {
train_data <- df_clinique[folds != i, ]
test_data <- df_clinique[folds == i, ]
fit_cv <- lm(anxiete_gad7 ~ nevrosisme + coping_evitement + sommeil_rem, data = train_data)
preds <- predict(fit_cv, newdata = test_data)
cv_errors[i] <- mean((test_data$anxiete_gad7 - preds)^2)
}
rmsep_cv <- sqrt(mean(cv_errors))
cat("Erreur quadratique moyenne de prediction hors-echantillon (RMSEP) :", round(rmsep_cv, 3))
Parallèlement, une démarche de rééchantillonnage non paramétrique (bootstrap) peut être engagée pour éprouver la stabilité structurelle du choix de variables opéré par le Cp. En générant 1000 réplicats de l’échantillon initial et en ré-exécutant l’algorithme ols_step_best_subset() sur chaque tirage avec remise, le statisticien peut quantifier la fréquence relative à laquelle la combinaison ternaire {névrosisme, coping par évitement, sommeil} est formellement isolée. Une fréquence de réplication élevée (par exemple > 85%) confirme la robustesse écologique du modèle face aux fluctuations aléatoires du recrutement clinique.
12. Bonnes pratiques, limites méthodologiques et considérations éditoriales en recherche quantitative
12.1 Le biais induit par la sélection de variables (post-selection inference)
L’utilisation d’algorithmes de sélection de modèles — qu’ils reposent sur le critère de Mallows, sur l’AIC ou sur des procédures pas à pas — soulève une préoccupation épistémologique majeure souvent passée sous silence dans les manuels appliqués : le problème de l’inférence post-sélection (post-selection inference ou biais d’optimisation par les données). Lorsque la même matrice empirique est exploitée simultanément pour identifier le sous-ensemble optimal de covariables et pour tester la significativité statistique des coefficients retenus, les postulats fondamentaux des tests d’hypothèse classiques sont violés.
En sélectionnant les régresseurs qui maximisent l’ajustement empirique tout en minimisant RSSp, l’algorithme capitalise inévitablement sur le hasard d’échantillonnage. En conséquence directe :
- Les valeurs de p issues du tableau
summary()du modèle final sont artificiellement trop optimistes (inflation du taux d’erreur de type I, assimilable à une forme involontaire de p-hacking). - Les amplitudes des coefficients de régression standardisés sont systématiquement surestimées en valeur absolue (phénomène connu sous le nom d’illusion du gagnant ou winner’s curse).
- Les intervalles de confiance construits selon la loi de Student sont trop étroits, produisant un faux sentiment de précision métrique.
Conformément aux recommandations éditoriales de l’American Psychological Association (APA) et de l’Open Science Framework, les chercheurs doivent expliciter sans ambiguïté dans leurs manuscrits si la spécification finale résulte d’une démarche purement exploratoire assistée par le Cp de Mallows ou d’un schéma d’analyse confirmatoire pré-enregistré. Idéalement, une partition initiale de l’échantillon (data splitting) réservant une fraction des données (ex. 30%) à la confirmation indépendante des inférences post-sélection permet de restaurer la validité nominale des tests d’hypothèse.
12.2 Limites structurelles du critère de Mallows
En dépit de son élégance matricielle et de son universalité dans le cadre de la régression linéaire standard, le critère Cp de Mallows présente des frontières structurelles d’applicabilité que le modélisateur doit impérativement respecter :
- Inapplicabilité aux Modèles Linéaires Généralisés (GLM) non gaussiens : La dérivation algébrique de l’espérance mathématique E[Cp] = p repose explicitement sur l’hypothèse d’additivité des résidus gaussiens homoscédastiques. Dans le contexte de la régression logistique binaire, de la régression de Poisson pour données de comptage ou des modèles de survie à risques proportionnels, le critère de Mallows n’est pas transposable directement sous sa forme classique. Le chercheur devra lui substituer l’AIC, le critère de déviance d’information (DIC) ou la déviance résiduelle généralisée.
- Vulnérabilité face aux non-linéarités complexes : Le critère de Mallows évalue des sous-espaces vectoriels linéaires. Si la relation véritable entre un régresseur et la variable réponse est hautement non linéaire (ex. fonctions sinusoïdales, seuils critiques, effets de saturation non linéaires) et que celle-ci n’a pas été préalablement linéarisée via des splines ou des transformations polynomiales explicites dans le modèle complet, le Cp échouera à quantifier adéquatement le biais de spécification.
- Sensibilité extrême à la spécification du modèle complet initial : Le critère Cp n’évalue jamais la validité intrinsèque absolue d’un modèle dans l’absolu ; il n’évalue qu’une distance relative par rapport à l’étalon constitué par le modèle complet. Si ce modèle complet est lui-même déficient, mal mesuré ou amputé de régresseurs structurels indispensables, l’ensemble de l’édifice comparatif est corrompu dès sa base.
- Concurrence des méthodes de régularisation moderne : Face à des matrices de très haute dimension où le nombre de prédicteurs P approche ou dépasse la taille d’échantillon n (problématique P ≥ n courante en neuro-imagerie ou en génomique comportementale), la régression complète ne peut plus être estimée par moindres carrés ordinaires en raison de la singularité de XTX. Le critère de Mallows devient mathématiquement inopérant. Dans de telles configurations, les techniques de régularisation pénalisée continue — telles que la régression Lasso (L1), la régression Ridge (L2) ou l’Elastic Net — s’imposent comme les alternatives contemporaines incontournables.
12.3 Recommandations de rédaction pour publications académiques
La présentation des résultats issus d’une sélection de modèles par le critère Cp de Mallows au sein d’une revue internationale à comité de lecture requiert une transparence méthodologique totale. L’espace éditorial ne doit pas se limiter à mentionner laconiquement le modèle sélectionné, mais doit documenter le processus d’arbitrage ayant guidé l’investigateur.
Les standards de communication scientifique imposent de respecter les consignes suivantes :
- Spécification exhaustive du modèle de référence : Déclarer formellement l’ensemble des variables incluses dans le modèle complet initial, leur échelle de mesure, le nombre d’observations effectives (après traitement contrôlé des données manquantes), ainsi que la variance résiduelle de référence s2 (avec ses degrés de liberté résiduels n – P).
- Présentation d’un tableau synthétique comparatif : Structurer les résultats sous forme d’une table normalisée alignant, pour les meilleurs sous-modèles compétitifs de chaque strate dimensionnelle : la liste des prédicteurs engagés, le nombre total de paramètres estimés p (incluant explicitement l’intercept), le R2, le R2 ajusté, la statistique Cp de Mallows, et l’écart absolu |Cp – p|.
- Justification formelle du choix final : Énoncer explicitement la règle heuristique mobilisée pour départager les spécifications concurrentes (ex. sélection du modèle présentant le Cp le plus bas sous la condition de proximité Cp ≤ p). Argumenter substantiellement la plausibilité théorique des variables écartées par l’arbitrage mathématique.
- Mise à disposition des scripts et des données : Conformément aux principes de la science ouverte (Open Science), déposer les scripts de calcul R rigoureusement commentés ainsi que les jeux de données anonymisés sur des dépôts institutionnels publics pérennes tels que l’Open Science Framework (OSF), Zenodo ou GitHub, garantissant la reproductibilité computationnelle immédiate des analyses par les pairs.
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
- American Psychological Association. (2020). Publication manual of the American Psychological Association (7th ed.). American Psychological Association. https://doi.org/10.1037/0000165-000
- Burnham, K. P., & Anderson, D. R. (2002). Model selection and multimodel inference: A practical information-theoretic approach (2nd ed.). Springer-Verlag. https://doi.org/10.1007/b97636
- Fox, J., & Weisberg, S. (2019). An R companion to applied regression (3rd ed.). SAGE Publications. https://socialsciences.mcmaster.ca/jfox/Books/Companion/
- Furnival, G. M., & Wilson, R. W. (1974). Regressions by leaps and bounds. Technometrics, 16(4), 499–511. https://doi.org/10.1080/00401706.1974.10489231
- Hebbali, A. (2020). olsrr: Tools for building OLS regression models (R package version 0.5.3). Comprehensive R Archive Network. https://cran.r-project.org/package=olsrr
- 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. https://cran.r-project.org/package=leaps
- Mallows, C. L. (1973). Some comments on Cp. Technometrics, 15(4), 661–675. https://doi.org/10.1080/00401706.1973.10489103
- R Core Team. (2023). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/
- Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6(2), 461–464. https://doi.org/10.1214/aos/1176344136
- Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L. D., François, R., Grolemund, G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Pedersen, T. L., Miller, E., Bache, S. M., Müller, K., Ooms, J., Robinson, D., Seidel, D. P., Spinu, V., … Yutani, H. (2019). Welcome to the Tidyverse. Journal of Open Source Software, 4(43), Article 1686. https://doi.org/10.21105/joss.01686