Dans le paysage contemporain de la modélisation statistique, de l’économétrie et de la psychométrie computationnelle, le problème du choix de modèle constitue l’un des défis méthodologiques les plus critiques auxquels font face les chercheurs et analystes de données. L’accroissement de la puissance de calcul et la disponibilité croissante de données hautement multidimensionnelles incitent fréquemment à formuler des modèles d’une complexité algorithmique démesurée. Cependant, l’augmentation du nombre de paramètres libres dans une spécification statistique s’accompagne inévitablement d’une vulnérabilité accrue au surapprentissage, réduisant dramatiquement la validité externe et la portée reproductible des inférences scientifiques. Face à cet écueil, les critères d’information pénalisés se sont imposés comme des outils d’arbitrage indispensables entre la qualité d’ajustement empirique et la parcimonie structurelle des modèles concurrents.
Le critère d’information bayésien, désigné sous l’acronyme BIC pour Bayesian Information Criterion ou critère de Schwarz, occupe une place prépondérante dans cette taxonomie décisionnelle. Ancré dans les fondements de la théorie de la décision bayésienne et de l’analyse asymptotique, le BIC propose une formalisation rigoureuse du célèbre rasoir d’Ockham en introduisant une pénalisation explicite et proportionnelle à la fois au nombre de coefficients estimés et au logarithme népérien de la taille de l’échantillon. Cette caractéristique lui confère une propriété asymptotique singulière, la consistance, qui le distingue fondamentalement d’autres métriques telles que le critère d’information d’Akaike (AIC).
L’environnement statistique R offre un écosystème d’une souplesse remarquable pour estimer, analyser et comparer le BIC sur une multitude d’architectures de modèles, allant des régressions linéaires par moindres carrés ordinaires jusqu’aux modèles linéaires généralisés et aux structures hiérarchiques multiniveaux. Ce guide exhaustif explore l’ensemble des ramifications théoriques, mathématiques et informatiques relatives au calcul du BIC dans R. En combinant dérivation analytique des équations sous-jacentes, décomposition manuelle d’algorithmes et manipulation approfondie des fonctions natives et spécialisées du langage, cette ressource ambitionne de doter la communauté scientifique d’une compréhension doctrinale et pratique sans compromis.
1. Introduction conceptuelle au critère d’information bayésien (BIC)
1.1 Définition et origines du critère de Schwarz
Le critère d’information bayésien a été formulé par le mathématicien et statisticien Gideon Schwarz dans un article séminal publié en 1978 au sein des Annals of Statistics, intitulé « Estimating the Dimension of a Model ». Dans ce travail fondateur, Schwarz cherchait à résoudre la problématique récurrente de la sélection de dimension dans les modèles linéaires emboîtés sans dépendre d’un seuil de signification conventionnel ou arbitraire, tel que la valeur p issue d’un test d’hypothèse de rapport de vraisemblance.
Sur le plan théorique, le BIC découle d’une approximation par la méthode du point-selle (ou approximation de Laplace) de l’intégrale définissant la vraisemblance marginale d’un modèle sous un a priori diffus ou faiblement informatif. Lorsque la taille d’échantillon n tend vers l’infini, le logarithme de la probabilité a posteriori intégrée d’un modèle statistique se décompose en un terme principal dominé par le maximum de la fonction de log-vraisemblance et un terme d’ordre inférieur reflétant l’incertitude paramétrique, formalisé par le produit de la dimension du modèle par le logarithme de la taille de l’échantillon.
Cette approche consacre une rupture paradigmatique par rapport à l’inférence statistique fréquentiste classique. Alors que la tradition de Neyman-Pearson s’attache à contrôler l’erreur de type I via des tests d’hypothèses nulles souvent rigides, la perspective initiée par Schwarz positionne le statisticien dans une logique de comparaison bayésienne où chaque modèle candidat se voit attribuer une probabilité relative d’avoir généré les données observées. Le principe de parcimonie ne relève dès lors plus d’une préférence philosophique subjective mais émerge naturellement de l’intégration analytique sur l’espace des paramètres.
1.2 Utilité méthodologique en modélisation psychologique et statistique
Dans les sciences comportementales, cognitives et sociales, les chercheurs sont confrontés à des systèmes complexes caractérisés par une variabilité résiduelle substantielle et de multiples sources potentielles de confusion. L’utilité méthodologique du BIC s’avère particulièrement critique dans la gestion du compromis universel entre le biais d’omission de variables et la variance d’estimation paramétrique. Ajouter des variables explicatives au sein d’une équation structurelle améliore inéluctablement le coefficient de détermination R-carré au sein de l’échantillon d’apprentissage, créant l’illusion d’une explication théorique plus robuste alors même que le modèle peut capturer un simple bruit stochastique.
La prévention du surapprentissage revêt une acuité majeure dans les devis de recherche expérimentaux et quasi-expérimentaux où la taille de l’échantillon est contrainte par des limites logistiques, financières ou éthiques. Un modèle saturé qui épouse fidèlement les singularités idiosynchratiques d’un échantillon restreint échouera inévitablement à reproduire ses prédictions sur de nouveaux participants ou dans d’autres laboratoires. Le BIC impose un coût numérique sévère pour chaque degré de liberté consommé, agissant comme un garde-fou méthodologique contre l’inflation non contrôlée des composantes latentes ou des termes d’interaction superflus.
En outre, l’application du critère de Schwarz s’avère déterminante pour la validation de théories psychométriques concurrentes. Qu’il s’agisse de trancher entre une structure factorielle unidimensionnelle et un modèle bifactoriel hiérarchique, ou de sélectionner le schéma de covariance le plus vraisemblable dans un modèle de croissance latente, le BIC fournit une métrique scalaire objective et standardisée. Cette quantification rigoureuse de l’évidence garantit une comparabilité transdisciplinaire indispensable à la reproductibilité des inférences scientifiques.
1.3 Critères de décision pour le meilleur modèle selon le BIC
La règle décisionnelle fondamentale présidant à l’évaluation des spécifications concurrentes repose sur un principe de minimisation : le modèle statistiquement optimal, au sein de l’ensemble fini des candidats soumis à l’évaluation, est invariablement celui qui présente la valeur numérique de BIC la plus faible. Il est essentiel de souligner que la valeur absolue brute d’un score BIC ne possède que peu de signification intrinsèque isolée ; seule la différence arithmétique relative entre deux scores renseigne le modélisateur sur l’intensité du soutien empirique en faveur de l’une ou l’autre architecture.
Une variation négative de score BIC lors de l’adjonction d’un prédicteur indique que le gain en adéquation empirique mesuré par la log-vraisemblance compense largement la pénalité logarithmique induite par le paramètre supplémentaire. Inversement, une augmentation du score BIC atteste que la diminution de l’erreur résiduelle demeure insuffisante face à la complexité structurelle nouvellement concédée, justifiant le rejet de cette spécification plus lourde.
Il importe enfin de distinguer conceptuellement la qualité d’ajustement globale de l’utilité prédictive. Le BIC ne valide pas l’exactitude ontologique d’un modèle statistique au sens absolu, mais classe un ensemble restreint de propositions théoriques quantitatives. Un modèle affichant le BIC minimal peut demeurer structurellement insatisfaisant si l’espace des modèles explorés souffre d’une mauvaise formulation théorique initiale. Le critère doit donc être mobilisé conjointement avec des inspections graphiques diagnostiques des résidus et des évaluations de pertinence théorique de fond.
2. Fondements mathématiques et dérivation analytique du BIC
2.1 Formulation générale basée sur la log-vraisemblance
La formulation analytique générale du critère d’information bayésien s’énonce selon l’équation canonique suivante :
BIC = -2 * ln(L) + k * ln(n)
Dans cette expression fondamentale, la composante L désigne la valeur maximale atteinte par la fonction de vraisemblance conjointe des données conditionnellement aux paramètres optimaux du modèle, calculée selon la méthode du maximum de vraisemblance. Le terme -2 * ln(L), souvent désigné sous le vocable de déviance non calibrée, quantifie l’adéquation empirique brute entre la distribution observée des variables dépendantes et la distribution théorique postulée par la famille paramétrique sélectionnée. Plus cette valeur est faible, plus le modèle théorique est capable de reproduire les régularités présentes dans la matrice empirique.
Le second terme, formalisé par le produit k * ln(n), constitue la pénalité paramétrique explicite. La variable scalaire k représente la dimension dimensionnelle intégrale du modèle, c’est-à-dire le nombre total de paramètres libres devant être estimés empiriquement à partir des données. Cela englobe non seulement les coefficients de régression associés aux variables explicatives, mais également l’ordonnée à l’origine (l’intercept) ainsi que tout paramètre d’échelle auxiliaire, tel que la variance de l’erreur résiduelle dans les modèles linéaires gaussiens. La grandeur n correspond quant à elle au nombre total d’observations indépendantes constituant l’échantillon d’analyse.
La validité mathématique de cette approximation asymptotique repose sur des hypothèses de régularité statistique strictes issues de la théorie du maximum de vraisemblance : l’espace des paramètres doit être compact, le véritable vecteur de paramètres ne doit pas se situer sur la frontière de cet espace, et la fonction de log-vraisemblance doit être deux fois continûment différentiable au voisinage du point d’estimation maximale, admettant une matrice d’information de Fisher définie positive.
2.2 Formulation alternative basée sur la somme des carrés des résidus (RSS)
Dans le cadre particulier et largement répandu de la régression linéaire multiple univariée, où l’on formule l’hypothèse canonique selon laquelle les résidus sont indépendants et identiquement distribués selon une loi normale d’espérance nulle et de variance homoscédastique constante, la fonction de log-vraisemblance s’exprime directement en fonction de la somme des carrés des résidus (RSS, pour Residual Sum of Squares). En substituant l’estimateur du maximum de vraisemblance de la variance résiduelle au sein de l’équation de vraisemblance gaussienne multivariée, on obtient une dérivation analytique spécifique du BIC :
BIC = n * ln(RSS / n) + k * ln(n) + C
Dans cette équation, la constante additive C est une fonction exclusive de la taille de l’échantillon n et de constantes universelles telles que pi (C = n * (ln(2 * pi) + 1)). Étant donné que cette constante demeure rigoureusement invariante d’un modèle à l’autre dès lors que la base de données analysée est strictement constante, elle est fréquemment omise ou neutralisée dans les comparaisons algorithmiques, conduisant à la formule opérationnelle suivante :
BIC = n * ln(RSS / n) + k * ln(n)
Il existe également dans la littérature économétrique une formulation normalisée par l’échantillon d’observation : BIC = (RSS + ln(n) * d * sigma_chapeau^2) / n, où d représente le nombre de régresseurs et sigma_chapeau^2 un estimateur sans biais de la variance résiduelle issue d’un modèle complet de référence. Néanmoins, l’équivalence formelle avec la formulation basée sur la log-vraisemblance maximale est préservée à une transformation linéaire près, garantissant une cohérence absolue dans l’ordonnancement final des spécifications concurrentes.
2.3 Comportement asymptotique et consistance du critère
La singularité théorique majeure du BIC réside dans sa propriété de consistance de sélection (model selection consistency). Un critère de sélection est qualifié de consistant si la probabilité qu’il sélectionne le véritable modèle générateur des données tend asymptotiquement vers 1 lorsque la taille de l’échantillon n croît vers l’infini, sous réserve formelle que ce modèle générateur fasse effectivement partie de l’ensemble fermé des spécifications concurrentes évaluées par le chercheur.
Cette consistance découle directement de la vitesse de croissance logarithmique de la pénalité paramétrique k * ln(n). Alors que la pénalité d’Akaike demeure fixée à 2k quelle que soit la masse de données accumulée, la pénalité du critère de Schwarz croît de manière monotone et non bornée avec n. Par conséquent, pour de très grands volumes d’échantillonnage, la pénalité par paramètre devient suffisamment dominante pour éliminer impitoyablement toute variable auxiliaire présentant un apport marginal insignifiant, écrasant ainsi la probabilité asymptotique de surparamétrisation.
Cependant, cette propriété remarquable présente un corollaire théorique qu’il convient de ne jamais occulter : si le « vrai » modèle sous-jacent est d’une dimensionnalité infinie ou s’il n’appartient pas à la famille des spécifications postulées (ce qui représente la quasi-totalité des contextes empiriques réels en psychologie et en sociologie), la sévérité de la pénalisation logarithmique peut contraindre le chercheur à sélectionner une approximation excessivement simpliste au détriment de l’efficacité prédictive globale.
3. Configuration de l’environnement R et préparation des données
3.1 Initialisation de la session et gestion des dépendances
Avant d’entreprendre l’implémentation algorithmique ou le calcul automatisé du BIC, il est impératif de configurer l’environnement de calcul R selon les préceptes de la science ouverte et de la recherche computationnelle reproductible. Une session d’analyse doit être initialisée en réinitialisant si nécessaire l’environnement global et en configurant des options globales d’affichage numérique claires afin de prévenir les troncatures scientifiques prématurées qui pourraient masquer de faibles deltas d’information.
L’utilisation de la commande options(scipen = 999, digits = 8) force l’interpréteur R à afficher les résultats sous forme décimale standard plutôt qu’en notation exponentielle, facilitant ainsi la détection des écarts subtils de vraisemblance entre modèles à fortes dimensions. Par ailleurs, la fixation explicite d’une graine de reproductibilité pseudo-aléatoire à l’aide de l’instruction set.seed(42) assure que toute procédure stochastique ultérieure, telle qu’une validation croisée ou une rééchantillonnage par bootstrap, produira des métriques strictement identiques lors des exécutions subséquentes.
Bien que les fonctions fondamentales de calcul du BIC soient intégrées de manière native au sein de la distribution standard de R via le package stats, l’importation de bibliothèques spécialisées s’avère indispensable pour manipuler des structures complexes. L’environnement d’analyse s’articulera avantageusement autour des packages tidyverse pour l’ingénierie des données, MASS pour les distributions avancées, lme4 pour les architectures hiérarchiques, ainsi que MuMIn et modelsummary pour l’automatisation tabulaire des diagnostics de sélection.
3.2 Importation et inspection structurelle du jeu de données
Pour illustrer les calculs et procédures décrits tout au long de cet article, nous mobilisons le jeu de données classique mtcars, disponible nativement dans R, tout en l’adaptant conceptuellement aux problématiques d’évaluation comportementale et de modélisation multivariée. Ce corpus compile des métriques de performance et de conception mécanique pour 32 véhicules automobiles, incluant des variables telles que la consommation de carburant (mpg, miles per gallon), la puissance du moteur (hp, horsepower), le poids du véhicule (wt, weight) et le type de transmission (am).
L’inspection structurelle préalable débute par l’interrogation systématique des dimensions globales de l’objet à l’aide des primitives dim(mtcars) et nrow(mtcars). La confirmation rigoureuse du nombre de lignes observées s’avère névralgique, car le scalaire n conditionne directement la magnitude de la pénalité logarithmique du BIC. Une altération involontaire de cette grandeur invaliderait irrémédiablement l’ensemble des dérivations théoriques subséquentes.
L’application de la fonction str(mtcars) et des utilitaires de diagnostic de complétude permet de vérifier l’absence de discordances dans le typage des prédicteurs (variables continues, discrètes ou catégorielles codées sous forme de facteurs) et d’identifier la présence éventuelle de données manquantes (NA). Dans les jeux de données observationnels ou psychométriques plus volumineux, ces vérifications initiales constituent la première ligne de défense contre les distorsions computationnelles de vraisemblance.
3.3 Nettoyage préalable et respect des prérequis d’échantillonnage
L’une des fautes méthodologiques les plus communes et les plus dévastatrices dans l’application des critères d’information réside dans la comparaison de modèles ajustés sur des échantillons de tailles différentes. Si un modèle A inclut des covariables dépourvues de valeurs manquantes tandis qu’un modèle B incorpore un prédicteur additionnel présentant des observations non renseignées, l’interpréteur R procèdera par défaut à une suppression par liste (listwise deletion) lors de l’estimation du modèle B. La taille d’échantillon effective n diminuera alors substantiellement pour cette seconde équation.
Une telle disparité rend la comparaison directe des valeurs de BIC strictement illégitime sur le plan mathématique. La fonction de vraisemblance d’un échantillon de taille n = 100 n’est en aucun cas comparable à la fonction de vraisemblance calculée sur n = 85, car la distribution de probabilité conjointe sous-jacente ne porte pas sur le même espace d’échantillonnage. En conséquence, les métriques d’information deviennent totalement incommensurables.
Pour immuniser l’analyse contre cet écueil, il est impératif d’établir un sous-ensemble d’apprentissage parfaitement complet en amont de toute modélisation en utilisant des instructions de filtrage explicites telles que na.omit() ou tidyr::drop_na() sur l’ensemble des variables participant à l’espace global des modèles candidats. Dans des contextes de recherche fondamentale où l’attrition des données est substantielle, des stratégies d’imputation multiple par équations chaînées doivent être privilégiées, le calcul du BIC devant alors être adapté pour agréger l’incertitude entre les multiples jeux de données imputés.
4. Calcul algorithmique et manuel du BIC à partir de la formule statistique
4.1 Extraction programmatique des composantes de la régression linéaire
Afin de démystifier le fonctionnement interne des fonctions de sélection de modèles sous R, il convient de programmer manuellement le calcul du BIC à partir de ses composantes scalaires fondamentales. Ajustons à cette fin un modèle de régression linéaire multiple de base cherchant à modéliser la consommation (mpg) à partir du poids (wt) et de la puissance (hp) au sein du jeu de données mtcars :
Le modèle est estimé via l’instruction modele_base <- lm(mpg ~ wt + hp, data = mtcars). À partir de cet objet ajusté de classe lm, nous pouvons extraire la somme des carrés des résidus (RSS) par le biais de la primitive residuals() combinée à l’opérateur de sommation : rss_manuel <- sum(residuals(modele_base)^2).
Simultanément, la taille effective de l’échantillon d’analyse est isolée de manière fiable à l’aide de la fonction nobs(modele_base), qui extrait directement l’attribut du nombre d’observations mobilisées dans l’estimation sans dépendre des dimensions de l’objet de données parent. Enfin, le dénombrement scrupuleux des paramètres k requiert une vigilance absolue : dans une régression linéaire classique, k comprend les coefficients associés aux prédicteurs, l’ordonnée à l’origine (intercept), ainsi que le paramètre d’échelle résiduel correspondant à la variance de l’erreur (sigma carré). Pour un modèle à deux prédicteurs, k s’établit ainsi à 2 + 1 + 1 = 4 paramètres libres.
4.2 Implémentation scriptée de la formule basée sur les résidus
Nous pouvons désormais formaliser cette logique au sein d’une fonction R personnalisée, démontrant la traduction directe des équations mathématiques en langage de programmation statistique :
L’écriture d’une telle routine prend la forme suivante :
calcul_bic_rss <- function(modele) {
rss <- sum(residuals(modele)^2)
n <- nobs(modele)
k <- length(coef(modele)) + 1
bic_val <- n * log(rss / n) + k * log(n)
return(bic_val)
}
Il importe de noter que cette fonction calcule la variante opérationnelle sans les constantes universelles d’intégration. Dans la littérature académique, certains auteurs omettent délibérément le terme additionnel correspondant au paramètre d’échelle de variance résiduelle dans le décompte de k, en considérant k uniquement comme le nombre de colonnes de la matrice du plan d’expérience. Bien que cette convention n’altère en rien les deltas comparatifs entre modèles partageant la même distribution d’erreur, elle introduit un décalage d’ordonnée dans la valeur absolue finale qu’il est capital de documenter scrupuleusement pour assurer la traçabilité des scripts d’analyse.
4.3 Implémentation scriptée basée sur la fonction logLik()
L’approche la plus universelle et rigoureuse sur le plan computationnel consiste à extraire la log-vraisemblance maximale maximisée de l’objet de modèle via la méthode générique logLik(). Cette primitive renvoie un objet scalaire enrichi d’un attribut précieux nommé df (degrés de liberté), lequel correspond exactement au nombre total k de paramètres libres estimés par la méthode du maximum de vraisemblance :
L’implémentation scriptée procède selon la logique suivante :
ll_obj <- logLik(modele_base)
valeur_ll <- as.numeric(ll_obj)
k_parametres <- attr(ll_obj, "df")
taille_n <- nobs(modele_base)
bic_loglik_manuel <- -2 * valeur_ll + k_parametres * log(taille_n)
La confrontation systématique de la valeur obtenue par cette dernière équation avec celle retournée par l’équation basée sur la somme des carrés des résidus met en lumière l’impact exact des constantes d’intégration normalisées. La méthode reposant sur logLik() incorpore le terme complet comprenant le facteur logarithmique de 2 * pi. Cette dérivation exacte constitue la référence canonique alignée sur la définition originale de Schwarz, garantissant une cohérence transversale absolue lorsque l’on compare des régressions linéaires avec des familles de modèles non linéaires.
5. Utilisation des fonctions natives de R pour le calcul du BIC
5.1 Calcul unitaire avec la fonction générique BIC() du package stats
Dans la pratique courante de l’analyse de données, il n’est heureusement pas nécessaire de reprogrammer systématiquement les équations fondamentales. Le package de base stats fournit une fonction générique polymorphe extrêmement puissante baptisée BIC(). Cette méthode prend en charge nativement un large éventail de classes d’objets statistiques produits par l’environnement R, allant des objets lm et glm jusqu’aux modèles de séries temporelles arima.
L’exécution de la syntaxe unitaire est d’une concision exemplaire : BIC(modele_base). Dès son appel, R procède à l’interrogation de la méthode logLik() appropriée pour la classe de l’objet passé en argument, détermine automatiquement le nombre total de paramètres via les attributs de degrés de liberté, interroge la taille de l’échantillon par la méthode nobs() correspondante, et restitue un scalaire numérique rigoureusement exact.
L’analyse comparative de la sortie produite par BIC(modele_base) face à notre implémentation manuelle basée sur logLik() démontre une identité arithmétique absolue jusqu’à la huitième décimale. Cette concordance totale valide le fait que la fonction native intègre scrupuleusement l’ensemble des constantes de normalisation de la vraisemblance gaussienne, y compris le terme de variance d’erreur.
5.2 Calcul simultané pour une collection de modèles emboîtés
L’un des atouts opérationnels majeurs de la primitive BIC() réside dans sa capacité à recevoir simultanément une pluralité d’objets de modèles au sein d’une instruction unique, générant instantanément une matrice récapitulative indispensable à l’analyse comparative. Supposons la définition de trois spécifications distinctes :
Un modèle univarié simple : mod1 <- lm(mpg ~ wt, data = mtcars)
Un modèle additif à deux prédicteurs : mod2 <- lm(mpg ~ wt + hp, data = mtcars)
Un modèle étendu avec interaction : mod3 <- lm(mpg ~ wt * hp, data = mtcars)
L’appel simultané tableau_comparatif <- BIC(mod1, mod2, mod3) produit un tableau structuré sous forme de data.frame comportant deux colonnes indispensables : df, explicitant le nombre de paramètres k associés à chaque spécification (ici 3 pour mod1, 4 pour mod2, et 5 pour mod3), et BIC, consignant les scores calculés respectifs. Cette structure tabulaire native permet un repérage immédiat de la dynamique d’évolution du critère face à l’accroissement des degrés de liberté consommés.
5.3 Alternative avec la fonction extractAIC() paramétrée
Il existe au sein du package stats une fonction alternative fréquemment mobilisée par les algorithmes de sélection automatique internes de R : la fonction extractAIC(). Contrairement à la fonction générique BIC(), extractAIC() a été conçue pour retourner le critère d’Akaike généralisé défini par la formule -2 * ln(L) + k * p, où p représente la pénalité paramétrique assignée.
Pour émuler rigoureusement le comportement du BIC à travers cette interface, il convient de modifier explicitement l’argument optionnel k en lui affectant la valeur du logarithme népérien de la taille d’échantillon, selon la syntaxe suivante :
taille_n <- nobs(mod1)
extractAIC(mod1, k = log(taille_n))
À l’exécution, le modélisateur observera une divergence numérique apparente entre le score retourné par extractAIC() et celui fourni par BIC(). Cette divergence ne découle aucunement d’une erreur mathématique, mais d’un choix algorithmique historique : extractAIC.lm() élimine l’ensemble des constantes résiduelles universelles indépendantes des paramètres de régression et ne comptabilise pas la variance de l’erreur dans la dimension k affichée. Les écarts relatifs (deltas) entre modèles concurrents évalués via extractAIC(..., k = log(n)) demeurent rigoureusement identiques à ceux calculés via BIC(), préservant intégralement l’ordre de préférence des modèles.
6. Application empirique : comparaison de modèles de régression linéaire multiple
6.1 Définition théorique des modèles concurrents
Pour illustrer la mise en œuvre pratique de la sélection de modèles guidée par le BIC, déployons un protocole expérimental complet visant à modéliser la consommation énergétique des véhicules au sein du jeu de données mtcars. Considérons la formulation de quatre hypothèses de recherche formalisées par des modèles économétriques concurrents :
L’hypothèse 1 postule un déterminisme purement mécanique lié à la masse du véhicule (Modèle M1 : mpg ~ wt). L’hypothèse 2 propose un modèle bivarié additif associant le poids et le rendement thermique du moteur exprimé par sa puissance brute (Modèle M2 : mpg ~ wt + hp). L’hypothèse 3 enrichit cette structure par l’ajout de la cylindrée du moteur (disp) et du rapport de pont arrière (drat), traduisant une architecture plus détaillée de la chaîne cinématique (Modèle M3 : mpg ~ wt + hp + disp + drat). Enfin, l’hypothèse 4 postule l’existence d’une dynamique non linéaire complexe en intégrant un terme quadratique sur le poids ainsi qu’une interaction de second ordre entre le poids et la puissance (Modèle M4 : mpg ~ wt + I(wt^2) + hp + wt:hp).
6.2 Exécution séquentielle et stockage des modèles
Afin de garantir une programmation modulaire et élégante, les modèles statistiques concurrents sont ajustés séquentiellement et archivés au sein d’une structure de liste nommée en langage R. Cette approche facilite l’application d’itérateurs fonctionnels et standardise les traitements ultérieurs :
L’ajustement procède selon les commandes suivantes :
liste_modeles <- list(
M1_Simple = lm(mpg ~ wt, data = mtcars),
M2_Additif = lm(mpg ~ wt + hp, data = mtcars),
M3_Etendu = lm(mpg ~ wt + hp + disp + drat, data = mtcars),
M4_Complexe = lm(mpg ~ wt + I(wt^2) + hp + wt:hp, data = mtcars)
)
L’extraction comparative des métriques d’ajustement conventionnelles met immédiatement en relief le piège classique des métriques non pénalisées :
sapply(liste_modeles, function(m) summary(m)$r.squared)
sapply(liste_modeles, function(m) summary(m)$adj.r.squared)
Le coefficient de détermination R-carré croît de manière monotone au fil de la complexité des modèles, atteignant son zénith absolu avec le modèle M4. Si un analyste non averti se fiait exclusivement à cette mesure, il opterait pour la spécification la plus lourde sans mesurer le risque de surparamétrisation encouru.
6.3 Calcul et mise en tableau des valeurs de BIC comparatives
Procédons à présent à l’extraction ordonnée des scores de BIC et à la dérivation des écarts d’information relatifs. L’application d’un script structuré permet de générer une table d’aide à la décision méthodologique :
tableau_bic <- data.frame(
Modele = names(liste_modeles),
k = sapply(liste_modeles, function(m) attr(logLik(m), "df")),
R2_Ajuste = sapply(liste_modeles, function(m) summary(m)$adj.r.squared),
BIC = sapply(liste_modeles, BIC)
)
Pour affiner l’analyse, calculons le différentiel par rapport au modèle le plus performant, défini formellement par Delta BIC = BIC_i – BIC_min, et ordonnons la structure tabulaire selon cette grandeur :
tableau_bic$\Delta_BIC <- tableau_bic$BIC - min(tableau_bic$BIC)
tableau_bic <- tableau_bic[order(tableau_bic$Delta_BIC), ]
L’examen des résultats révèle un enseignement capital : le modèle M2 (additif, poids et puissance) affiche le score de BIC le plus faible, s’imposant comme la spécification optimale. Bien que les modèles M3 et M4 affichent un R-carré brut supérieur, l’adjonction de variables supplémentaires ou de termes non linéaires complexes entraîne une pénalisation logarithmique k * ln(n) qui excède largement le gain modeste de réduction de l’erreur résiduelle. Le critère de Schwarz protège ainsi l’analyste contre l’illusion de complexité.
7. Calcul du BIC pour les modèles linéaires généralisés (GLM)
7.1 Particularités du calcul dans les modèles de régression logistique binaire
L’extension du calcul du critère d’information bayésien aux modèles linéaires généralisés (GLM) requiert la prise en compte des spécificités fonctionnelles de la famille exponentielle et de l’estimation par réestimation pondérée des moindres carrés (IRLS). Dans le cas d’une régression logistique binaire visant à modéliser la probabilité d’occurrence d’un événement dichotomique (tel que le type de transmission mécanique ou la réussite à un item psychométrique), l’équation d’ajustement est spécifiée via la primitive glm() combinée à l’argument family = binomial(link = "logit").
Dans cette architecture, la métrique fondamentale d’écart n’est plus la somme des carrés des résidus, mais la déviance résiduelle, définie formellement comme le double de la différence de log-vraisemblance entre le modèle saturé et le modèle estimé : Déviance = 2 * (ln(L_sature) – ln(L_modele)). La log-vraisemblance du modèle maximisé L_modele est directement exploitée par la fonction logLik.glm().
Lors de l’appel BIC(modele_glm), R applique avec exactitude la formule fondamentale -2 * ln(L) + k * ln(n). La composante k correspond ici au nombre de coefficients de régression logistique estimés (incluant l’ordonnée à l’origine). Contrairement à la régression linéaire sous hypothèse gaussienne, aucun paramètre additionnel de variance résiduelle n’est estimé, car la variance d’une loi de Bernoulli est intrinsèquement liée à son espérance mathématique p * (1 – p), ce dont la primitive native tient rigoureusement compte dans son décompte des degrés de liberté.
7.2 Régression de Poisson et distribution binomiale négative pour données de comptage
La modélisation de données discrètes de comptage, récurrentes dans l’étude des fréquences comportementales ou des récidives cliniques, repose conventionnellement sur la régression de Poisson (family = poisson(link = "log")). Une contrainte restrictive de la distribution de Poisson réside dans l’équidispersion, stipulant l’égalité théorique entre la moyenne conditionnelle et la variance. Dans les contextes empiriques réels, les données présentent fréquemment une surdispersion substantielle (la variance excédant la moyenne), ce qui biaise l’estimation des erreurs-types et fausse les tests de vraisemblance.
Face à ce phénomène, l’analyste peut hésiter entre une spécification de quasi-Poisson ou une régression binomiale négative. Il convient de souligner que les modèles de quasi-Poisson ne définissant pas une véritable fonction de vraisemblance conjointe (ils reposent sur une quasi-vraisemblance sans intégration de densité), la méthode logLik() renvoie une valeur indéterminée (NA), rendant le calcul formel du BIC impossible dans ce cadre strict.
Pour arbitrer rigoureusement entre une structure équidispersée et une distribution prenant en charge la surdispersion, il est nécessaire d’ajuster un modèle binomial négatif via la fonction glm.nb() du package MASS. La distribution binomiale négative intégrant un paramètre de dispersion additionnel (thêta), ce paramètre libre supplémentaire est automatiquement comptabilisé dans k par logLik(). Dès lors, la comparaison directe BIC(modele_poisson, modele_binomial_negatif) permet de déterminer objectivement si la prise en compte de la surdispersion compense la perte d’un degré de liberté additionnel.
7.3 Détermination rigoureuse de n dans les modèles GLM
Une complexité subtile mais majeure relative à l’application du critère de Schwarz aux modèles linéaires généralisés réside dans la définition précise de la grandeur n intervenant dans le terme de pénalité logarithmique k * ln(n). Lorsque les données binomiales sont structurées sous un format agrégé (c’est-à-dire sous la forme d’un tableau synthétisant le nombre de succès et le nombre d’échecs pour des covariables groupées, manipulé via une variable réponse à deux colonnes cbind(succes, echecs)), l’interpréteur R attribue par défaut à n le nombre de lignes du tableau agrégé, et non le nombre total d’essais individuels de Bernoulli sous-jacents.
Une telle disparité induit une sous-estimation systématique et sévère de la taille réelle de l’échantillon, conduisant à une pénalisation paramétrique anormalement clémente. Pour garantir une rigueur méthodologique absolue lors du traitement de formats groupés, il est impératif d’ajuster manuellement la valeur de n en extrayant la somme totale des effectifs observés.
Dans un tel scénario, le calcul du BIC rectifié doit être exécuté par une formulation personnalisée :
taille_reelle <- sum(donnees_groupees$succes + donnees_groupees$echecs)
k_params <- attr(logLik(modele_glm), "df")
bic_corrige <- -2 * as.numeric(logLik(modele_glm)) + k_params * log(taille_reelle)
Cette rectification évite le risque de retenir des spécifications surparamétrées résultant d’une illusion d’échantillon restreint.
8. Comparaison approfondie : BIC vs AIC et critères dérivés
8.1 Divergences philosophiques et structurelles entre BIC et AIC
Le débat relatif à l’arbitrage méthodologique entre le critère d’Akaike (AIC) et le critère de Schwarz (BIC) constitue l’une des controverses fondamentales de la statistique appliquée. Cette divergence prend racine dans des fondements épistémologiques et mathématiques radicalement contrastés. L’AIC, formulé par Hirotugu Akaike en 1974, vise à minimiser la divergence de Kullback-Leibler entre le modèle postulé et la véritable réalité génératrice des données. Sa formulation générale s’écrit AIC = -2 * ln(L) + 2k.
La distinction mathématique centrale entre ces deux métriques réside exclusivement dans la structure de leur pénalisation paramétrique : tandis que l’AIC impose une pénalité fixe de 2 unités pour chaque paramètre additionnel quelle que soit la masse de données analysées, le BIC applique une pénalité dynamique égale à ln(n) par paramètre. Une simple résolution algébrique de l’inéquation ln(n) > 2 permet d’identifier le seuil critique à partir duquel le comportement des deux critères diverge :
e^2 ≈ 7,389
Dès lors que la taille de l’échantillon n est supérieure ou égale à 8 observations, ln(n) excède strictement 2. Par conséquent, pour toute analyse standard comportant plusieurs dizaines, centaines ou milliers d’observations, le critère de Schwarz sanctionne la complexité paramétrique avec une sévérité considérablement supérieure à celle imposée par l’AIC.
8.2 Risques comparés : surparamétrisation contre sous-paramétrisation
Cette divergence structurelle de pénalisation engendre des profils de risque asymétriques qu’il est capital de pondérer en fonction des objectifs scientifiques poursuivis. En vertu de sa pénalité invariante face à l’accroissement d’échantillon, l’AIC n’est pas un critère consistant. Lorsque n croît de manière illimitée, la probabilité que l’AIC sélectionne un modèle contenant des variables superflues ne tend pas vers zéro ; elle se stabilise autour d’un taux d’erreur de surparamétrisation documenté dans la littérature asymptotique comme pouvant atteindre approximativement 15 % à 20 % dans les régressions univariées.
À l’inverse, l’AIC optimise l’efficacité prédictive à échantillon fini (propriété d’optimalité minimax). Si la réalité empirique est caractérisée par une multitude d’effets faibles infiniment imbriqués sans frontière nette, l’AIC fournira un modèle doté d’une meilleure capacité de généralisation prédictive hors échantillon.
Le BIC, en revanche, privilégie résolument la parcimonie et la consistance. Face à un échantillon de taille importante, il identifiera avec une quasi-certitude le sous-ensemble minimal de variables véritablement déterminantes. Cependant, cette intransigeance algorithmique comporte un revers : dans des contextes d’échantillons réduits à modérés (n < 100), la sévérité de ln(n) peut conduire à une sous-paramétrisation excessive, éliminant des covariables porteuses d’un signal théorique réel mais affichant une puissance statistique marginalement insuffisante.
8.3 Présentation des critères alternatifs : AICc, QAIC, WAIC
Pour surmonter les limites intrinsèques des critères classiques dans des configurations empiriques complexes, plusieurs métriques alternatives ont été développées. Dans les contextes où le ratio entre le nombre d’observations et le nombre de paramètres est faible (n / k < 40), l’AIC standard sous-pénalise sévèrement les modèles. Sugiura (1978) et Hurvich & Tsai (1989) ont ainsi introduit l’AIC corrigé (AICc), qui intègre un terme de pénalité quadratique de second ordre :
AICc = AIC + (2k * (k + 1)) / (n – k – 1)
Dans l’écosystème R, le package MuMIn et la bibliothèque AICcmodavg permettent de calculer simultanément l’AICc et le BIC à travers une interface standardisée :
library(MuMIn)
model.sel(mod1, mod2, mod3, mod4, rank = "BIC")
Face à des structures de comptage affichant une surdispersion majeure, le critère quasi-AIC (QAIC) incorpore une estimation empirique du paramètre de surdispersion c-chapeau au sein de la fonction de vraisemblance. Enfin, dans le paradigme de la modélisation bayésienne computationnelle par chaînes de Markov Monte Carlo (MCMC), les critères traditionnels sont avantageusement supplantés par le WAIC (Widely Applicable Information Criterion) ou critère de Watanabe, accessible sous R via le package loo. Le WAIC présente la particularité remarquable d’évaluer la log-vraisemblance ponctuelle sur l’ensemble de la distribution a posteriori sans s’appuyer sur une simple approximation asymptotique de point-selle.
9. Sélection automatisée de variables basée sur le BIC dans R
9.1 Approches pas à pas avec la fonction step()
La recherche d’une spécification parcimonieuse au sein d’un espace combinatoire étendu de régresseurs peut être automatisée sous R à l’aide de la primitive step(). Par défaut, cette fonction implémente une procédure séquentielle guidée par l’AIC classique (la pénalité par défaut étant fixée à k = 2). Toutefois, il est tout à fait aisé de contraindre cet algorithme à naviguer selon le critère d’information bayésien en reparamétrant explicitement l’argument de pénalité par le logarithme népérien de la taille de l’échantillon d’analyse : k = log(nobs(modele_sature)).
Considérons un modèle saturé incorporant la totalité des régresseurs disponibles au sein du jeu de données :
mod_complet <- lm(mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb, data = mtcars)
n_taille <- nobs(mod_complet)
selection_bic <- step(mod_complet, direction = "both", k = log(n_taille), trace = TRUE)
L’argument direction = "both" active un processus bidirectionnel progressif (alternant adjonctions et suppressions de termes). À chaque itération de la boucle de calcul, l’algorithme évalue la valeur du BIC résultant du retrait ou de l’adjonction de chaque variable candidate. L’affichage textuel détaillé (trace = TRUE) consigne scrupuleusement la trajectoire descendante de la déviance pénalisée, arrêtant la procédure lorsqu’aucune modification structurelle ne permet d’obtenir un score BIC inférieur à celui du modèle courant.
9.2 Sélection exhaustive de sous-ensembles avec leaps
Bien que les procédures pas à pas (stepwise) soient computationnellement efficientes, elles souffrent d’une limite heuristique notoire : elles explorent des trajectoires locales et ne garantissent aucunement l’identification de l’optimum global au sein de l’espace des modèles. Pour pallier cette faiblesse, une approche par recherche exhaustive de tous les sous-ensembles possibles (all-subsets regression) s’avère nettement supérieure lorsque le nombre total de régresseurs demeure inférieur à une trentaine.
Le package leaps met à disposition la fonction hautement optimisée regsubsets(), qui déploie l’algorithme de branchement et coupure (branch-and-bound) de Furnival et Wilson pour évaluer l’intégralité des combinaisons sans ajuster explicitement chaque régression :
library(leaps)
selection_exhaustive <- regsubsets(mpg ~ ., data = mtcars, nvmax = 10, method = "exhaustive")
synthese_leaps <- summary(selection_exhaustive)
L’extraction de la composante synthese_leaps$bic permet d’isoler avec une rigueur absolue la combinaison de prédicteurs minimisant globalement le critère de Schwarz pour chaque taille de modèle k donnée. La fonction graphique associée plot(selection_exhaustive, scale = "bic") génère une matrice visuelle cartographiant les covariables retenues dans les architectures les plus performantes, offrant une vision synoptique immédiate de la structure de sélection.
9.3 Limites épistémologiques et écueils de la sélection automatisée
En dépit de leur séduction pratique évidente, les algorithmes de sélection automatique de variables guidés par les critères d’information soulèvent des réserves méthodologiques majeures documentées par les statisticiens théoriciens. Le dragage systématique de données (data dredging) génère une inflation drastique des risques d’erreurs d’inférence statistique. Lorsqu’un algorithme passe au crible des dizaines de combinaisons, les estimations finales des coefficients de régression ainsi que leurs erreurs-types se trouvent substantiellement biaisées en raison de la sélection conditionnelle.
En présence d’une forte multicolinéarité entre les variables explicatives, les procédures de sélection basées sur le BIC affichent une instabilité extrême. Une modification infime d’une observation ou la suppression d’une seule ligne d’échantillon peut modifier radicalement l’identité des régresseurs sélectionnés par l’algorithme, remplaçant arbitrairement une variable par une autre avec laquelle elle est étroitement corrélée.
Par conséquent, les critères de décision algorithmiques ne sauraient en aucun cas se substituer à une conceptualisation théorique préalable rigoureuse. La sélection de modèles doit demeurer guidée par des hypothèses scientifiques déductives étayées par la littérature du domaine, le calcul du BIC servant alors d’outil d’arbitrage entre un nombre restreint de spécifications théoriquement justifiées plutôt que de moteur de génération aveugle de formulations quantitatives.
10. Calcul du BIC pour les modèles multiniveaux et mixtes avec lme4
10.1 Particularités des modèles à effets mixtes (LMM et GLMM)
L’analyse de données imbriquées, hiérarchiques ou longitudinales (telles que des mesures répétées sur des participants ou des élèves nichés au sein de classes scolaires) repose sur les modèles linéaires à effets mixtes (LMM) et leurs homologues généralisés (GLMM). Dans l’environnement R, le package de référence pour l’estimation de ces structures est lme4, articulé autour des commandes lmer() et glmer().
Dans un modèle mixte, la variance de la variable réponse est décomposée entre des coefficients structurels fixes décrivant les tendances moyennes de la population et des composantes de variance aléatoires (ordonnées à l’origine et pentes aléatoires) modélisant l’hétérogénéité spécifique aux sous-unités d’agrégation. Cette complexité structurelle modifie profondément le calcul de la log-vraisemblance maximale ainsi que le dénombrement effectif des degrés de liberté libres k.
L’appel de la méthode unitaire BIC() sur un objet de classe merMod (produit par lmer) applique une formule qui englobe dans k non seulement les coefficients de régression fixe, mais également l’ensemble des éléments uniques de la matrice de covariance des effets aléatoires ainsi que la variance du terme d’erreur résiduel de niveau 1.
10.2 Le dilemme de la taille d’échantillon n dans les structures hiérarchiques
L’application du critère de Schwarz aux architectures multiniveaux se heurte à un problème théorique d’envergure, fréquemment occulté dans les applications empiriques : la définition ambiguë de la taille d’échantillon effective n. Dans une structure où J groupes de niveau 2 regroupent un total cumulé de N observations élémentaires de niveau 1, quelle grandeur doit être injectée dans la pénalité logarithmique k * ln(n) ?
Par défaut, l’implémentation native de BIC.merMod() dans le package lme4 utilise la taille d’échantillon totale de niveau 1 (n = N, c’est-à-dire le nombre total d’observations individuelles). Si cette stratégie s’avère parfaitement adéquate pour évaluer l’impact de prédicteurs de niveau 1, elle impose une pénalisation excessivement sévère lors de l’évaluation de paramètres ou de prédicteurs mesurés exclusivement au niveau agrégé (niveau 2).
La littérature statistique méthodologique contemporaine (notamment les travaux de Delattre et al., 2014) préconise dans certains contextes la partition du terme de pénalité en fonction de l’échelle d’action des covariables, en appliquant ln(N) pour les composantes de niveau 1 et ln(J) pour les composantes de groupe de niveau 2. Les chercheurs utilisant lme4 doivent être pleinement conscients de ce comportement par défaut, documentant explicitement dans leurs protocoles que la pénalité appliquée par R repose sur la taille globale de niveau micro.
10.3 Comparaison de structures d’effets aléatoires
Lors de l’arbitrage entre différentes configurations d’effets aléatoires au sein d’un modèle mixte, une règle d’estimation absolue doit être impérativement respectée pour que le calcul du BIC demeure valide : les modèles concurrents doivent impérativement être ajustés par la méthode du Maximum de Vraisemblance ordinaire (ML, spécifiée via REML = FALSE au sein de l’instruction lmer).
Par défaut, la fonction lmer() recourt au Maximum de Vraisemblance Restreint (REML), car ce dernier produit des estimations sans biais des composantes de variance dans les échantillons finis. Cependant, la log-vraisemblance restreinte sous REML est conditionnée par la structure spécifique des effets fixes. Par conséquent, comparer les valeurs de BIC de modèles présentant des structures d’effets fixes différentes estimées sous REML constitue une erreur computationnelle grave conduisant à des inférences totalement erronées.
Illustrons une séquence rigoureuse d’évaluation de la structure aléatoire :
library(lme4)
lmm_intercept_seul <- lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy, REML = FALSE)
lmm_pente_aleatoire <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy, REML = FALSE)
BIC(lmm_intercept_seul, lmm_pente_aleatoire)
Cette comparaison systématique permet d’établir objectivement si la corrélation et l’hétérogénéité inter-individuelle des trajectoires de réponse temporelle justifient l’adjonction des paramètres de covariance aléatoires supplémentaires face à la pénalité logarithmique du critère de Schwarz.
11. Interprétation avancée : poids de Schwarz, probabilités et facteurs de Bayes
11.1 Calcul des deltas BIC et seuils interprétatifs de Raftery
Une fois les scores de BIC extraits pour l’ensemble des modèles concurrents, l’analyse ne saurait se restreindre à la simple désignation mécanique de la valeur numérique minimale. Pour quantifier la force de l’évidence relative apportée par les données en faveur du meilleur modèle par rapport à chaque compétiteur alternatif, l’analyste doit impérativement calculer l’écart relatif noté Delta BIC (ou ΔBIC) :
Delta_i = BIC_i – BIC_min
Par définition, le modèle le plus performant affiche un Delta de 0, tandis que tous les modèles alternatifs présentent des valeurs positives. Pour structurer l’interprétation qualitative de ces grandeurs scalaires, le statisticien Adrian Raftery a établi en 1995 une grille d’évaluation empirique universellement reconnue dans la communauté scientifique internationale :
- Un écart compris entre 0 et 2 : Évidence très faible ou négligeable en faveur du modèle préféré. Les données ne permettent pas de discriminer significativement les deux spécifications, qui doivent être considérées comme empiriquement équivalentes.
- Un écart compris entre 2 et 6 : Évidence positive modérée. Le modèle présentant le BIC le plus bas bénéficie d’un soutien substantiel face à son concurrent.
- Un écart compris entre 6 et 10 : Évidence forte. Le soutien empirique en faveur de la spécification optimale est marqué, justifiant le rejet de l’alternative.
- Un écart strictement supérieur à 10 : Évidence très forte ou décisive. Le modèle alternatif présente une probabilité conditionnelle quasi-nulle d’adéquation relative.
11.2 Transformation en poids de BIC (poids de Schwarz) et probabilités a posteriori
L’une des propriétés les plus élégantes du critère d’information bayésien réside dans sa convertibilité directe en probabilités a posteriori de modèles à travers la formulation des « poids de BIC » (ou poids de Schwarz), notés w_i. Sous l’hypothèse canonique d’indifférence a priori assignant une probabilité uniforme équiprobable à chaque modèle avant observation des données (P(M_i) = 1 / K), le poids de Schwarz formalise la probabilité conditionnelle que le modèle M_i soit effectivement le meilleur modèle au sein de l’espace fermé des spécifications examinées :
w_i = exp(-0.5 * Delta_i) / sum_{j=1}^K exp(-0.5 * Delta_j)
Ces poids sont rigoureusement bornés au sein de l’intervalle unitaire [0, 1] et leur sommation sur l’ensemble des K modèles concurrents est rigoureusement égale à l’unité. Leur implémentation sous forme d’algorithme R vectorisé s’exécute avec une clarté remarquable :
deltas <- tableau_bic$Delta_BIC
poids_schwarz <- exp(-0.5 * deltas) / sum(exp(-0.5 * deltas))
tableau_bic$Poids <- round(poids_schwarz, 4)
Cette formalisation probabiliste offre un support de communication remarquablement intuitif pour le rapportage scientifique. Elle permet d’affirmer, par exemple, qu’un modèle particulier concentre 82 % de la probabilité a posteriori totale d’adéquation empirique parmi l’ensemble des architectures théoriques postulées par l’équipe de recherche.
11.3 Approximation des facteurs de Bayes à partir du BIC
Dans le cadre de l’inférence bayésienne rigoureuse, le test comparatif direct de deux hypothèses ou modèles concurrents M1 et M0 s’exprime par le facteur de Bayes (BF_10), qui quantifie l’actualisation du rapport de cotes opérée par l’observation des données empiriques : Cote a posteriori = BF_10 * Cote a priori.
Cependant, le calcul exact du facteur de Bayes requiert l’évaluation analytique d’intégrales hautement complexes et multidimensionnelles sur l’ensemble de l’espace des paramètres, nécessitant des techniques computationnelles lourdes telles que l’échantillonnage de Gibbs ou l’intégration par ponts (bridge sampling). Schwarz a démontré que le différentiel de score BIC fournit une approximation asymptotique directe de premier ordre du logarithme du facteur de Bayes :
ln(BF_10) ≈ -0.5 * (BIC_1 – BIC_0) = -0.5 * Delta BIC_{1-0}
En transformant cette égalité par la fonction exponentielle, on obtient la formule opérationnelle fondamentale :
BF_10 ≈ exp(-0.5 * (BIC_1 – BIC_0))
Si la différence BIC_1 – BIC_0 est égale à -6 (indiquant que le modèle 1 surpasse le modèle 0 de 6 unités de BIC), le facteur de Bayes BF_10 approximatif s’établit à exp(3) ≈ 20,08. Ce résultat indique concrètement que les données observées sont approximativement 20 fois plus probables sous la formulation du modèle 1 que sous celle du modèle 0, offrant ainsi une quantification bayésienne limpide sans nécessiter la spécification laborieuse de lois a priori conjuguées ou d’algorithmes MCMC complexes.
12. Diagnostic, pièges classiques et bonnes pratiques sous R
12.1 Erreurs fréquentes d’implémentation et de comparaison
L’utilisation non critique du critère d’information bayésien au sein des environnements d’analyse statistique automatisés expose les chercheurs à des erreurs d’interprétation méthodologique destructrices. La première d’entre elles concerne la non-invariance stricte de l’échantillon d’entraînement. Comme développé précédemment, il est impératif d’auditer systématiquement la conformité de l’attribut nobs() sur chacun des objets de modèles comparés. Toute divergence, ne serait-ce que d’une seule observation, invalide formellement l’ordonnancement des BIC.
Un deuxième piège majeur concerne la comparaison de modèles appliquant des transformations non linéaires différentes à la variable réponse dépendante. Considérons une analyse comparant un modèle linéaire direct lm(y ~ x) avec un modèle log-transformé lm(log(y) ~ x). L’application directe de la fonction BIC() sur ces deux objets retournera des scores numériques que le chercheur ne doit en aucun cas comparer directement.
La raison théorique de cette interdiction réside dans le théorème de changement de variable des densités de probabilité continues : la densité conjointe d’une variable transformée log(y) intègre le jacobien de la transformation (1 / y). La log-vraisemblance retournée par logLik() pour le modèle log-transformé omet le facteur correctif d’échelle -sum(log(y)). Pour comparer légitimement ces modèles via le BIC, il est impératif d’ajouter manuellement ce terme jacobien à la log-vraisemblance du modèle log-transformé avant d’appliquer l’équation de Schwarz.
12.2 Conflits entre modèles non emboîtés et validation croisée
Contrairement aux tests statistiques de rapport de vraisemblance classiques (tels que la fonction anova(m1, m2, test = "Chisq")), qui imposent une relation formelle d’emboîtement structurel (l’espace paramétrique du modèle restreint devant être un sous-espace strict de celui du modèle général), le critère d’information bayésien s’applique en toute rigueur à la comparaison de modèles non emboîtés. Cette propriété lui confère une flexibilité considérable pour trancher entre des théories scientifiques incompatibles formulées par des systèmes d’équations distincts.
Néanmoins, les conclusions formulées sur la base exclusive du BIC peuvent parfois entrer en conflit avec les métriques empiriques obtenues par validation croisée (k-fold cross-validation). En vertu de sa pénalité asymptotiquement lourde k * ln(n), le BIC manifeste une tendance conservatrice prononcée qui le conduit fréquemment à sélectionner des modèles plus compacts que ceux favorisés par une validation croisée à 10 blocs visant la minimisation stricte de l’erreur quadratique moyenne de prédiction (RMSEP).
Face à une telle divergence empirique, le chercheur doit réconcilier ses diagnostics en explicitant son objectif d’analyse principal : si la visée ultime du projet réside dans la prédiction pure hors échantillon sans préoccupation majeure pour l’explicabilité théorique, les arbitrages issus de la validation croisée doivent primer. Si, en revanche, l’investigation ambitionne l’identification structurelle des véritables déterminants causaux au sens de l’explication théorique parcimonieuse, la consistance du modèle désigné par le BIC offre une garantie théorique plus robuste.
12.3 Protocole de rapportage scientifique selon les standards académiques
La transcription des résultats d’une sélection de modèles guidée par le critère d’information bayésien au sein d’une publication scientifique évaluée par les pairs doit obéir à des standards formels de transparence et de reproductibilité computationnelle. Les normes de l’American Psychological Association (APA 7th edition) et des principales revues d’économétrie préconisent la présentation des architectures concurrentes au sein d’un tableau récapitulatif synoptique et normalisé.
Ce tableau méthodologique doit obligatoirement expliciter pour chaque modèle examiné : la désignation conceptuelle claire de la spécification théorique, la log-vraisemblance maximale maximisée (ln(L)), le nombre de paramètres effectifs estimés (k), le score BIC résultant, le différentiel d’information relatif (Delta BIC), ainsi que le poids de probabilité a posteriori normalisé (w_i). L’inclusion du coefficient R-carré ajusté fournit une information complémentaire sur la fraction de variance expliquée.
Dans l’écosystème logiciel R, la formalisation automatisée de ces tableaux éditoriaux est facilitée par l’utilisation de la bibliothèque modelsummary, qui permet d’exporter directement les structures d’objets vers des formats LaTeX, HTML ou Microsoft Word sans réencodage manuel :
library(modelsummary)
modelsummary(liste_modeles, output = "table_bic.html", metrics = c("BIC", "AIC", "RMSE", "R2"), stars = TRUE)
Enfin, en accord avec les exigences éthiques de la science ouverte (Open Science), l’intégralité du script R ayant présidé au nettoyage des données, à la fixation de la taille d’échantillon effective n, au dénombrement méticuleux des paramètres et au calcul des métriques de Schwarz doit être mise à disposition publique au sein d’un dépôt institutionnel pérenne (tel que l’Open Science Framework ou Zenodo), garantissant l’auditabilité computationnelle absolue des inférences publiées.
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
- Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01
- 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
- Burnham, K. P., & Anderson, D. R. (2004). Multimodel inference: Understanding AIC and BIC in model selection. Sociological Methods & Research, 33(2), 261–304. https://doi.org/10.1177/0049124104268644
- Delattre, M., Lavielle, M., & Poursat, M. A. (2014). A note on BIC in mixed-effects models. Electronic Journal of Statistics, 8(1), 456–475. https://doi.org/10.1214/14-EJS890
- 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
- Hurvich, C. M., & Tsai, C. L. (1989). Regression and time series model selection in small samples. Biometrika, 76(2), 297–307. https://doi.org/10.1093/biomet/76.2.297
- Kass, R. E., & Raftery, A. E. (1995). Bayes factors. Journal of the American Statistical Association, 90(430), 773–795. https://doi.org/10.1080/01621459.1995.10476572
- R Core Team. (2023). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/
- Raftery, A. E. (1995). Bayesian model selection in social research. Sociological Methodology, 25, 111–163. https://doi.org/10.2307/271063
- Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6(2), 461–464. https://doi.org/10.1214/aos/1176344136
- Sugiura, N. (1978). Further analysis of the data by Akaike’s information criterion and the finite corrections. Communications in Statistics – Theory and Methods, 7(1), 13–26. https://doi.org/10.1080/03610927808827599
- Vehtari, A., Gelman, A., & Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27(5), 1413–1432. https://doi.org/10.1007/s11222-016-9696-4
- Wagenmakers, E. J. (2007). A practical solution to the pervasive problems of p values. Psychonomic Bulletin & Review, 14(5), 779–804. https://doi.org/10.3758/BF03194105