Comment effectuer un test de Breusch-Godfrey dans R

Dans le domaine de l’économétrie appliquée et de la modélisation statistique des sciences du comportement, l’évaluation de la validité interne d’un modèle de régression linéaire repose fondamentalement sur le diagnostic rigoureux de ses résidus. Parmi les postulats classiques formulés par le théorème de Gauss-Markov, l’hypothèse selon laquelle les termes d’erreur ne présentent aucune corrélation mutuelle conditionnellement aux régresseurs constitue une pierre angulaire pour garantir que l’estimateur des moindres carrés ordinaires (MCO) demeure le meilleur estimateur linéaire non biaisé (BLUE). Lorsque cette condition d’orthogonalité temporelle ou sérielle est rompue, l’ensemble de l’édifice inférentiel s’effondre : bien que les coefficients estimés puissent conserver leur propriété de convergence et de non-biais sous certaines conditions strictes d’exogénéité, leurs erreurs-types se trouvent lourdement sous-estimées ou surestimées, invalidant de fait les tests de Student, les analyses de variance et les intervalles de confiance associés.

Face à cette vulnérabilité méthodologique, les chercheurs ont longtemps recouru de manière quasi exclusive au test de Durbin-Watson. Cependant, ce dernier présente des limitations mathématiques rédhibitoires, notamment son incapacité structurelle à détecter des autocorrélations d’ordre supérieur à un et son invalidation théorique absolue dès lors que le modèle intègre des retards de la variable dépendante parmi les variables explicatives. Pour surmonter ces contraintes restrictives, Trevor S. Breusch et Leslie G. Godfrey ont développé de façon indépendante à la fin des années 1970 une procédure diagnostique générale fondée sur le principe du multiplicateur de Lagrange. Ce test, universellement désigné sous le nom de test de Breusch-Godfrey, s’est imposé comme le standard d’excellence pour traquer les dépendances sérielles complexes dans les modèles de régression contemporains.

Le présent guide méthodologique propose une exploration exhaustive, théorique et opérationnelle du test de Breusch-Godfrey au sein de l’environnement de programmation statistique R. Conçu à l’intention des chercheurs, des biostatisticiens, des économètres et des psychométriciens confrontés à des données séquentielles, longitudinales ou issues de protocoles d’évaluation écologique momentanée, cet article détaille chaque étape du processus analytique. Depuis la dérivation mathématique des statistiques de test jusqu’à la mise en œuvre de solutions d’estimation robustes face aux violations empiriques, en passant par le calibrage fin du paramètre de retard et l’interprétation académique des résultats, ce document fournit les clés nécessaires pour garantir la parfaite rigueur de vos inférences statistiques.

1. Introduction à l’autocorrélation des résidus en modélisation statistique

1.1 Le postulat d’indépendance des erreurs dans le modèle linéaire

Dans la formulation matricielle classique du modèle de régression linéaire générale, exprimée sous la forme Y = Xβ + ε, où Y représente le vecteur des observations de dimension (n × 1), X la matrice des régresseurs de dimension (n × k), β le vecteur des paramètres inconnus de dimension (k × 1) et ε le vecteur des perturbations stochastiques inobservables, la théorie statistique impose un ensemble d’hypothèses sur le comportement probabiliste de ε. L’une des exigences fondamentales concerne la structure de la matrice de variance-covariance des erreurs conditionnelle à X. Mathématiquement, cette condition s’énonce par l’égalité Var(ε|X) = E(εε’|X) = σ²Iₙ, où σ² désigne la variance scalaire constante du terme d’erreur et Iₙ représente la matrice identité d’ordre n. Ce postulat condense deux propriétés indispensables : d’une part, l’homoscédasticité, qui requiert que la diagonale principale de la matrice soit constituée d’éléments identiques représentant une dispersion constante, et d’autre part, la non-autocorrélation ou indépendance conditionnelle sérielle, imposant que tous les éléments situés hors de la diagonale principale soient rigoureusement nuls, soit Cov(εᵢ, εⱼ | X) = 0 pour tout i ≠ j.

Lorsque cette hypothèse d’indépendance conditionnelle est violée, c’est-à-dire lorsqu’il existe une corrélation systématique entre les perturbations associées à différentes observations temporelles ou spatiales, le célèbre théorème de Gauss-Markov cesse d’être vérifié. Dès lors, bien que l’estimateur usuel des moindres carrés ordinaires, donné par β̂ = (X’X)⁻¹X’Y, demeure sans biais et convergent pourvu que l’hypothèse d’exogénéité stricte E(ε|X) = 0 reste respectée, il perd irrévocablement sa propriété d’efficacité statistique. En d’autres termes, β̂ n’est plus l’estimateur à variance minimale au sein de la classe des estimateurs linéaires sans biais (BLUE). Il existe des estimateurs alternatifs, tels que l’estimateur des moindres carrés généralisés (GLS), qui affichent une variance d’échantillonnage substantiellement inférieure.

La conséquence la plus délétère de l’autocorrélation ne réside toutefois pas uniquement dans cette perte d’efficacité intrinsèque, mais plutôt dans la déformation structurelle de l’estimateur conventionnel de la matrice de variance-covariance des paramètres. La formule classique calculée par les logiciels statistiques standards, Var̂(β̂) = s²(X’X)⁻¹, repose implicitement sur l’intégrité de la matrice σ²Iₙ. En présence de corrélations non nulles entre résidus, la véritable matrice de variance-covariance asymptotique adopte la forme dite en sandwich : Var(β̂|X) = (X’X)⁻¹ X’ Ω X (X’X)⁻¹, où Ω = E(εε’) représente la matrice de covariance réelle des perturbations. Dans la quasi-totalité des configurations empiriques caractérisées par une autocorrélation positive entre observations successives, la formule conventionnelle sous-évalue drastiquement la véritable dispersion des estimateurs, comprimant artificiellement les erreurs-types de chaque coefficient.

Cet effondrement artificiel des erreurs-types a des répercussions désastreuses sur l’ensemble de l’inférence statistique formelle. Les statistiques de test t de Student associées à chaque régresseur, calculées par le ratio entre le coefficient estimé et son erreur-type, voient leurs valeurs gonflées de manière mécanique. Ce phénomène conduit le modélisateur à rejeter l’hypothèse nulle d’absence d’effet avec une fréquence largement supérieure au seuil nominal fixé α. Le taux d’erreur de type I — consistant à identifier une association statistiquement significative là où il n’existe en réalité qu’une covariation fortuite induite par la mémoire sérielle du bruit de fond — se trouve considérablement exacerbé, atteignant fréquemment des proportions alarmantes de 30 % à 50 % au lieu des 5 % théoriques attendus. Les intervalles de confiance se révèlent trompeusement étroits, conférant aux estimations une précision illusoire qui compromet gravement la réplicabilité des résultats.

1.2 Origines de l’autocorrélation dans les données empiriques et psychométriques

L’émergence d’une autocorrélation résiduelle au sein des modélisations statistiques appliquées n’est que rarement un phénomène purement aléatoire ; elle reflète le plus souvent des propriétés intrinsèques du protocole de recueil des données ou des défauts structurels dans la spécification du modèle mathématique. Dans les sciences comportementales, cognitives et psychométriques, l’adoption croissante de protocoles longitudinaux à haute fréquence, tels que l’évaluation écologique momentanée (Ecological Momentary Assessment ou EMA) ou les méthodes d’échantillonnage de l’expérience (Experience Sampling Methods ou ESM), expose directement les chercheurs à ce défi. Lorsque les participants sont invités à renseigner des échelles psychologiques plusieurs fois par jour via leur ordiphone, les intervalles temporels séparant les mesures deviennent si resserrés que les états psychophysiologiques internes des sujets n’ont pas le temps de revenir à leur niveau de base, générant une continuité stochastique résiduelle d’une observation à la suivante.

Une source majeure et récurrente de corrélation sérielle réside dans l’omission de variables explicatives pertinentes qui manifestent elles-mêmes une inertie temporelle marquée. Les construits psychologiques tels que la fatigue cognitive, l’humeur dépressive diffuse, l’état de stress perçu ou le niveau de vigilance ne fluctuent pas de manière discontinue ou erratique ; ils évoluent selon des trajectoires continues dotées d’une forte dépendance à leur état passé. Si un modèle de régression tente de prédire la performance cognitive d’un individu en omettant d’intégrer une mesure explicite de son niveau d’épuisement accumulé au fil des heures, cet épuisement non capturé est relégué dans le terme d’erreur ε. Puisque la fatigue accumulée à l’instant t dépend étroitement de celle éprouvée à l’instant t-1, le terme d’erreur hérite mathématiquement de cette persistance dynamique, violant instantanément le postulat d’orthogonalité temporelle.

Au-delà de l’omission de variables régulatrices, la mauvaise spécification de la forme fonctionnelle liant les prédicteurs au critère d’intérêt constitue un mécanisme générateur d’autocorrélation tout aussi prégnant. Si la véritable relation sous-jacente entre une variable indépendante continue (comme le niveau d’anxiété pré-test) et une variable dépendante (telle que l’efficience de la mémoire de travail) suit une courbe quadratique ou une loi en U inversé de Yerkes-Dodson, mais que le chercheur impose arbitrairement une spécification strictement linéaire, les résidus du modèle présenteront une succession systématique d’écarts positifs et négatifs au fur et à mesure que le régresseur progresse. Lorsque les données sont triées chronologiquement ou collectées selon un continuum temporel, cette courbure non modélisée se traduit directement par des séquences de résidus de même signe, générant artificiellement un schéma d’autocorrélation d’ordre un ou supérieur.

Enfin, les protocoles expérimentaux en psychologie cognitive impliquant des blocs d’essais répétés souffrent fréquemment d’effets de report dynamique (carryover effects) et de processus d’apprentissage séquentiel. Dans une tâche de temps de réaction à choix multiples comportant plusieurs centaines d’essais consécutifs, la réponse d’un participant à l’essai t est systématiquement conditionnée par les événements survenus lors des essais t-1, t-2 ou t-3, à l’instar du phénomène bien documenté de ralentissement post-erreur (post-error slowing) ou de la facilitation par amorçage inter-essais. Même après avoir contrôlé les conditions expérimentales contemporaines, les fluctuations résiduelles de la vitesse de traitement manifestent une corrélation sérielle positive substantielle. Ignorer cette structure temporelle lors de la modélisation revient à attribuer faussement à des manipulations expérimentales une part de variance qui découle en réalité de l’inertie propre du système cognitif étudié.

1.3 Objectif et pertinence du test de Breusch-Godfrey

L’objectif central du test de Breusch-Godfrey réside dans sa capacité à fournir une procédure de diagnostic statistique rigoureuse, formelle et hautement versatile, permettant de statuer sur la présence éventuelle de dépendances sérielles au sein des termes d’erreur d’une régression linéaire, sans se restreindre aux cas de figure simplistes. Contrairement aux approches exploratoires purement visuelles fondées sur l’examen des résidus, qui demeurent dépendantes de la subjectivité de l’analyste et s’avèrent peu concluantes face à des bruits stochastiques complexes, le test de Breusch-Godfrey dérive un critère de décision quantitatif doté de propriétés asymptotiques parfaitement caractérisées. Il permet de tester simultanément ou séquentiellement l’existence d’autocorrélations s’étendant jusqu’à un ordre de retard arbitraire p, couvrant ainsi l’ensemble des dynamiques autorégressives ou de moyennes mobiles susceptibles de contaminer les données.

La pertinence méthodologique du test de Breusch-Godfrey devient éclatante dès lors que l’on considère les architectures de modélisation dynamique modernes, en particulier les modèles autorégressifs qui incluent la variable dépendante retardée parmi les prédicteurs (par exemple, Yₜ = β₀ + β₁Xₜ + γYₜ₋₁ + εₜ). Dans une telle configuration, les tests traditionnels s’avèrent totalement inopérants : la présence d’une variable dépendante retardée corrélée avec des résidus potentiellement autocorrélés biaise simultanément l’estimation des coefficients de régression et fausse la distribution d’échantillonnage des statistiques de test classiques. Le test de Breusch-Godfrey, grâce à sa formulation issue du principe du multiplicateur de Lagrange qui conditionne le test sur les régresseurs sous l’hypothèse nulle, conserve son entière validité asymptotique en présence de variables endogènes retardées, comblant ainsi une faille théorique majeure de la boîte à outils économétrique usuelle.

Dans le paradigme contemporain de la science ouverte et de la reproductibilité des analyses psychométriques, le recours systématique au test de Breusch-Godfrey s’inscrit au cœur de la phase de validation des hypothèses avant toute tentative d’interprétation substantielle. La rigueur inférentielle exige que la structure stochastique du modèle soit scrupuleusement purgée de toute dépendance sérielle non théorisée avant de revendiquer la significativité statistique d’une intervention thérapeutique, d’un effet développemental ou d’un prédicteur cognitif. Le test de Breusch-Godfrey ne constitue pas simplement un contrôle technique secondaire ; il représente un rempart indispensable contre la prolifération de résultats faussement positifs dans la littérature empirique, fournissant une justification mathématique irréfutable pour le maintien du modèle MCO ou, à l’inverse, pour l’adoption impérative de corrections d’erreurs-types ou de réestimations dynamiques adaptées.

2. Fondements théoriques et formulation mathématique du test de Breusch-Godfrey

2.1 Formulation des hypothèses statistiques nulle et alternative

La dérivation mathématique du test de Breusch-Godfrey repose sur l’hypothèse que la perturbation stochastique εₜ du modèle linéaire standard peut suivre un processus autorégressif d’ordre p, noté AR(p), ou de manière tout à fait équivalente sur le plan asymptotique, un processus de moyenne mobile d’ordre p, noté MA(p). Dans le cadre de la spécification autorégressive, le comportement dynamique du terme d’erreur s’écrit formellement sous la forme :

εₜ = ρ₁ εₜ₋₁ + ρ₂ εₜ₋₂ + … + ρₚ εₜ₋ₚ + uₜ

où ρ₁, ρ₂, …, ρₚ représentent les coefficients autorégressifs capturant l’intensité et le sens de la mémoire temporelle résiduelle pour chaque intervalle de retard allant de 1 jusqu’à l’horizon p, et où uₜ désigne une perturbation stochastique idéale, qualifiée de bruit blanc, satisfaisant aux conditions E(uₜ) = 0, Var(uₜ) = σ²ᵤ et Cov(uₜ, uₛ) = 0 pour tout t ≠ s. Dans ce cadre paramétrique rigoureux, l’hypothèse nulle de non-autocorrélation des résidus jusqu’à l’ordre p s’énonce par la nullité conjointe et absolue de l’ensemble des coefficients autorégressifs du processus stochastique :

H₀ : ρ₁ = 0, ρ₂ = 0, …, ρₚ = 0

Sous cette hypothèse nulle H₀, le processus se contracte instantanément en εₜ = uₜ, ce qui garantit que les perturbations du modèle d’origine ne manifestent aucune mémoire sérielle et que la condition de non-autocorrélation requise par le théorème de Gauss-Markov est rigoureusement satisfaite. À l’opposé, l’hypothèse alternative H₁ stipule l’existence d’une dépendance temporelle active au sein du système résiduel, ce qui se traduit mathématiquement par la négation de la conjonction précédente :

H₁ : ∃ j ∈ {1, 2, …, p} tel que ρⱼ ≠ 0

Il suffit donc qu’un seul des p coefficients autorégressifs présente une divergence statistiquement significative par rapport à zéro pour que l’hypothèse nulle soit formellement rejetée, signalant la persistance d’une corrélation sérielle qui compromet la validité des inférences standards. Cette formulation multidimensionnelle confère au test une grande flexibilité d’investigation par rapport aux tests unidirectionnels focalisés sur un retard unique.

2.2 Dérivation du principe du multiplicateur de Lagrange (LM)

L’élégance fondamentale de la procédure conçue par Breusch et Godfrey découle de son ancrage dans le principe du multiplicateur de Lagrange (souvent désigné sous le terme de test du score de Rao). Dans la théorie classique des tests d’hypothèses sur la fonction de vraisemblance, le test de Wald exige l’estimation non contrainte du modèle sous l’hypothèse alternative (ce qui impliquerait ici d’estimer conjointement les paramètres de régression et l’ensemble des p coefficients autorégressifs par des méthodes non linéaires complexes), tandis que le test du rapport de vraisemblance nécessite d’estimer à la fois le modèle contraint et le modèle non contraint. Le principe du multiplicateur de Lagrange présente l’avantage analytique déterminant de ne requérir que l’estimation du modèle contraint par l’hypothèse nulle — à savoir la simple régression linéaire ordinaire Y = Xβ + ε estimée par les MCO conventionnels.

La mise en œuvre pratique de ce principe repose sur la construction d’une régression auxiliaire opérant sur les résidus observés eₜ = yₜ – x’ₜβ̂, qui constituent des approximations convergentes des termes d’erreur inobservables εₜ. Dans cette régression auxiliaire, le résidu contemporain eₜ est régressé simultanément sur l’ensemble complet des régresseurs d’origine X et sur les p résidus retardés eₜ₋₁, eₜ₋₂, …, eₜ₋ₚ. Pour les p premières observations du jeu de données où les retards des résidus ne sont pas définis, la pratique canonique consiste à remplacer ces valeurs manquantes par le chiffre zéro, ce qui correspond à leur espérance mathématique sous l’hypothèse nulle. La forme explicite de la régression auxiliaire s’établit ainsi :

eₜ = x’ₜγ + α₁ eₜ₋₁ + α₂ eₜ₋₂ + … + αₚ eₜ₋ₚ + vₜ

L’inclusion des régresseurs originaux x’ₜ au sein de cette équation auxiliaire est d’une importance mathématique capitale : elle garantit l’orthogonalisation complète des résidus retardés par rapport à l’espace vectoriel engendré par X. Cette projection assure que les coefficients α₁ à αₚ capturent strictement l’autocorrélation sérielle pure, en purgeant les estimations de toute covariation potentielle avec les variables explicatives, ce qui préserve la validité du test même si X contient des retards de Y. Une fois cette régression auxiliaire estimée par les moindres carrés ordinaires, on extrait son coefficient de détermination non centré ou centré standard, désigné par R²ₑ. La statistique de test du multiplicateur de Lagrange, couramment notée LM, se définit alors comme le produit scalaire entre le nombre total d’observations n et le coefficient de détermination issu de la régression auxiliaire :

LM = n × R²ₑ

Sous l’hypothèse nulle H₀, et en vertu du théorème central limite associé aux propriétés de convergence stochastique des estimateurs à mémoire courte, la quantité LM converge en loi lorsque la taille d’échantillon n tend vers l’infini vers une distribution asymptotique parfaitement tabulée. Si les résidus ne manifestent aucune dépendance sérielle, le pouvoir prédictif des résidus retardés sur le résidu contemporain est quasi nul, ce qui contraint R²ₑ à osciller au voisinage immédiat de zéro, maintenant ainsi la statistique LM à un niveau modeste. Inversement, si une autocorrélation substantielle existe, les résidus retardés expliqueront une proportion non négligeable de la variance de eₜ, augmentant R²ₑ et propulsant la valeur de LM au-delà des valeurs critiques admissibles sous l’hypothèse nulle.

2.3 Distributions théoriques et critères de décision

Sous l’ensemble des conditions de régularité asymptotique requises pour les modèles de séries temporelles stationnaires et ergodiques, la statistique du multiplicateur de Lagrange LM = n × R²ₑ suit asymptotiquement, sous l’hypothèse nulle H₀, une distribution du Chi-deux dont le nombre de degrés de liberté est rigoureusement égal à l’ordre de retard maximal spécifié p :

LM = n × R²ₑ ~ᵃ Chi²(p)

Cette convergence vers une loi à p degrés de liberté découle directement du fait que la restriction sous H₀ impose la nullité conjointe de précisément p paramètres indépendants (les coefficients α₁ à αₚ de la régression auxiliaire). En pratique, cette propriété asymptotique fournit un critère de décision univoque en grand échantillon : si la valeur empirique calculée de LM excède la valeur critique théorique Chi²_{1-α}(p) pour un seuil de signification choisi α (communément fixé à 0.05), ou si la probabilité critique asymptotique (p-valeur) associée vérifie P(Chi²(p) ≥ LM) < α, le chercheur se trouve contraint de rejeter formellement l’hypothèse nulle H₀ au profit de l’hypothèse alternative H₁, concluant à la présence d’une autocorrélation sérielle significative au sein des erreurs de modélisation.

Cependant, dans le contexte d’échantillons de taille modérée ou réduite (typiquement lorsque n est inférieur à une centaine d’observations), la statistique asymptotique du Chi-deux tend à manifester une distorsion de taille non négligeable, conduisant à une sur-rejection excessive de l’hypothèse nulle. Pour pallier cette faiblesse d’approximation en petit échantillon, Godfrey et Harvey ont proposé une alternative robuste fondée sur la statistique usuelle de Fisher F de significativité globale des coefficients α₁ à αₚ dans la régression auxiliaire. Cette variante empirique se calcule selon la formule standard :

F = [ R²ₑ / p ] / [ (1 – R²ₑ) / (n – k – p) ]

où k représente le nombre total de régresseurs inclus dans le modèle initial (incluant la constante) et p désigne le nombre de retards résiduels introduits. Sous l’hypothèse nulle, cette statistique suit approximativement une distribution de Fisher à p degrés de liberté au numérateur et (n – k – p) degrés de liberté au dénominateur : F ~ F(p, n – k – p). La règle de décision s’adapte de façon symétrique : si la p-valeur associée à la statistique F observée tombe sous le seuil α, l’indépendance des résidus est rejetée. L’adoption de cette variante de Fisher est particulièrement préconisée dans les recherches expérimentales en sciences sociales où les tailles d’échantillons demeurent contraintes par le coût de passation des protocoles.

3. Analyse comparative : Test de Durbin-Watson versus Test de Breusch-Godfrey

3.1 Restrictions structurelles du test de Durbin-Watson

Développé au début des années 1950 par James Durbin et Geoffrey Watson, le test de Durbin-Watson a constitué pendant des décennies le passage obligé de tout diagnostic de régression. Sa formulation mathématique repose sur le calcul de la statistique d, définie par la somme des carrés des différences successives des résidus divisée par la somme globale des résidus au carré : d = ∑ₜ₌₂ⁿ (eₜ – eₜ₋₁)² / ∑ₜ₌₁ⁿ eₜ². Par un simple développement algébrique, il apparaît clairement que la statistique d entretient une relation approximative directe avec le coefficient d’autocorrélation empirique d’ordre un r₁, telle que d ≈ 2(1 – r₁). Si cette élégance calculatoire était précieuse avant l’avènement de l’informatique moderne, elle enferme la procédure dans des limitations structurelles sévères qui restreignent considérablement son champ d’application contemporain.

La première restriction structurelle réside dans son incapacité ontologique à investiguer des structures de dépendance sérielle au-delà du premier ordre. Le test de Durbin-Watson est conçu pour tester exclusivement l’hypothèse d’un processus autorégressif stationnaire de type AR(1). Dans de nombreux contextes appliqués, les données temporelles ne manifestent pas nécessairement de dépendance immédiate d’un instant à l’autre, mais exhibent des corrélations saisonnières, hebdomadaires ou cycliques (par exemple, à l’ordre 4 pour des séries trimestrielles, à l’ordre 7 pour des journaux de bord quotidiens, ou à des ordres supérieurs dans des tâches cognitives séquentielles). Face à un processus pur d’ordre supérieur, tel qu’un processus AR(2) où ρ₁ = 0 mais ρ₂ ≠ 0, la statistique d de Durbin-Watson prendra une valeur très proche de 2, indiquant faussement une parfaite indépendance des résidus et masquant complètement l’existence d’une sévère autocorrélation d’ordre 2.

La seconde limitation, bien plus grave sur le plan théorique, touche à l’invalidation mathématique de la statistique de Durbin-Watson dès lors que le modèle de régression incorpore une variable dépendante retardée (par exemple, Yₜ₋₁) parmi ses variables prédictives. Dans cette situation, la condition d’indépendance entre les régresseurs et les termes d’erreur futurs ou passés est structurellement rompue. La statistique d se trouve alors asymptotiquement biaisée vers la valeur 2 : elle tend vers 2 même lorsque les erreurs sous-jacentes sont lourdement autocorrélées. Ce phénomène anesthésie le pouvoir de détection du test, induisant systématiquement un risque massif d’erreur de type II (non-rejet d’une hypothèse nulle pourtant fausse). Enfin, la distribution théorique de la statistique d sous H₀ dépend de la matrice des régresseurs X spécifique à chaque jeu de données, ce qui a contraint les auteurs à établir des tables historiques délimitant une borne inférieure d_L et une borne supérieure d_U, engendrant une zone d’inconclusion frustrante où l’analyste ne peut ni accepter ni rejeter formellement l’hypothèse d’indépendance.

3.2 Avantages méthodologiques du test de Breusch-Godfrey

Face au cadre contraignant imposé par Durbin-Watson, le test de Breusch-Godfrey se distingue par une supériorité méthodologique absolue, conférée par sa dérivation générale sous le principe du multiplicateur de Lagrange. Son premier atout majeur réside dans sa totale transparence quant à l’ordre de dépendance investigué. La procédure ne se limite aucunement au schéma AR(1) : elle permet d’examiner avec une rigueur statistique équivalente n’importe quelle structure autorégressive AR(p) ou de moyenne mobile MA(p) jusqu’à un ordre p arbitrairement fixé par le chercheur. Cette polyvalence offre une flexibilité indispensable pour modéliser des périodicités complexes et pour s’assurer qu’aucun écho dynamique sous-jacent ne subsiste au sein des résidus de modélisation.

Le second avantage théorique décisif du test de Breusch-Godfrey réside dans sa robustesse démontrée face à la présence de régresseurs endogènes retardés. Comme souligné précédemment, la régression auxiliaire du test de Breusch-Godfrey incorpore obligatoirement dans sa spécification l’intégralité des variables explicatives d’origine du modèle, incluant toute variable dépendante retardée Yₜ₋₁. Cette inclusion assure la projection orthogonale adéquate des résidus retardés et garantit que les estimateurs des coefficients autorégressifs α₁ à αₚ convergent vers zéro sous l’hypothèse nulle sans subir de biais inductif. Par conséquent, la statistique LM conserve sa distribution limite théorique Chi²(p) même dans les modèles autorégressifs les plus denses, fournissant un instrument diagnostique infaillible là où Durbin-Watson produit des résultats trompeurs.

Enfin, le test de Breusch-Godfrey élimine définitivement l’inconvénient archaïque des zones d’indécision statistique. Alors que l’utilisateur de Durbin-Watson se heurtait historiquement à des intervalles critiques non concluants entre d_L et d_U, la procédure de Breusch-Godfrey associe à chaque calcul une p-valeur continue exacte ou asymptotique, fondée sur des distributions théoriques universelles et univoques (le Chi-deux ou la loi de Fisher). Cette clarté décisionnelle facilite l’automatisation computationnelle des diagnostics au sein des chaînes de traitement statistique et permet de formaliser des règles de décision explicites, parfaitement reproductibles et communicables dans les publications scientifiques internationales.

3.3 Critères de sélection du test adapté selon le devis de recherche

Le choix raisonné de la procédure de détection de l’autocorrélation ne doit pas s’effectuer de manière automatique ou machinale, mais doit découler d’une analyse méticuleuse du devis de recherche, de la structure temporelle des observations et de la complexité dynamique du modèle théorique sous-jacent. Dans les plans d’échantillonnage simples où les observations sont chronologiquement espacées de manière régulière, où le modèle ne comporte rigoureusement aucun retard de la variable critère et où la théorie substantielle n’anticipe qu’un simple processus de persistance immédiate d’un intervalle à l’autre, le recours conjoint à Durbin-Watson et à Breusch-Godfrey pour p = 1 aboutira généralement à des conclusions strictement concordantes. Cependant, dès lors que l’une de ces conditions idéalisées vient à défaillir, l’arbitrage méthodologique doit systématiquement pencher en faveur de Breusch-Godfrey.

L’évaluation de la périodicité et de la fréquence de recueil des données empiriques constitue le premier critère directeur. Dans un devis de suivi psychométrique quotidien s’étalant sur plusieurs semaines ou mois, les variations de l’humeur, du stress ou des interactions sociales manifestent fréquemment une dynamique hebdomadaire marquée (effet de cycle de travail versus week-end). L’application d’un test limité à l’ordre 1 échouera totalement à capturer une dépendance sérielle se manifestant avec un décalage de sept jours. L’analyste se doit dans cette configuration de spécifier un test de Breusch-Godfrey paramétré avec un retard d’ordre p = 7, seul à même de tester l’orthogonalité des résidus sur l’ensemble de la fenêtre temporelle pertinente pour le phénomène étudié.

Le second arbitrage fondamental concerne le compromis inhérent entre la puissance statistique du test et le risque de surparamétrisation lié à la sélection de l’ordre p. Si le test de Breusch-Godfrey possède la faculté d’investiguer des retards élevés, l’augmentation inconsidérée de p gonfle artificiellement les degrés de liberté du Chi-deux ou de la distribution F, ce qui a pour effet direct de diluer la puissance du test si l’autocorrélation réelle est circonscrite aux tout premiers ordres. De surcroît, chaque retard additionnel imposé dans la régression auxiliaire consomme une observation supplémentaire (ou impose d’imputer des zéros initiaux), altérant les propriétés de petit échantillon. Les critères directeurs pour sélectionner l’architecture du test s’articulent ainsi :

• Structure causale statique sans retards endogènes et petit échantillon (n < 50) : Le test de Breusch-Godfrey dans sa version statistique F de Fisher avec un ordre restreint (p = 1 ou p = 2) assure un contrôle optimal du taux d’erreur de type I sans subir l’incongruité des bornes de Durbin-Watson.
• Modèle dynamique incluant la variable dépendante retardée : L’utilisation de Durbin-Watson est formellement proscrite ; l’application de Breusch-Godfrey est impérative, idéalement configurée pour tester les ordres de retards correspondant aux cycles théoriques du phénomène.
• Données longitudinales à haute résolution temporelle (EMA, capteurs physiologiques, neurosciences) : L’exécution conjointe de corrélogrammes et d’une série de tests de Breusch-Godfrey échelonnés (de p = 1 à p = k, avec k calibré selon la granularité temporelle) constitue l’approche méthodologique d’excellence pour valider l’indépendance stochastique des termes d’erreur.

4. Configuration de l’environnement R et gestion des dépendances

4.1 Installation et chargement des packages requis

Pour mener à bien un diagnostic rigoureux de l’autocorrélation sous l’environnement R, le recours à un ensemble de bibliothèques logicielles spécialisées et reconnues par la communauté scientifique internationale s’avère indispensable. Le cœur computationnel du diagnostic de corrélation sérielle est implémenté au sein de la bibliothèque lmtest (Linear Model Diagnostics), développée par Achim Zeileis et Torsten Hothorn. Cette bibliothèque fournit un éventail complet de procédures d’inférence et de tests d’hypothèses linéaires, parmi lesquelles figure la fonction bgtest(), implémentant le test de Breusch-Godfrey de manière native et hautement optimisée.

En complément de lmtest, l’installation de la bibliothèque sandwich est vivement recommandée. Ce package permet l’estimation de matrices de variance-covariance hétéroscédastiques et autocorrélées (matrices HAC), s’avérant indispensable dès lors que le test de Breusch-Godfrey vient à rejeter formellement l’hypothèse d’indépendance des erreurs et qu’une stratégie de remédiation par erreurs-types robustes s’impose. Par ailleurs, la bibliothèque car (Companion to Applied Regression), maintenue par John Fox et Sanford Weisberg, apporte des outils diagnostics additionnels de haut niveau, facilitant la détection conjointe de la multicolinéarité et de la non-normalité. Enfin, l’écosystème tidyverse, articulé autour de dplyr, tibble et ggplot2, est requis pour structurer efficacement les matrices de données chronologiques et générer des visualisations résiduelles d’une irréprochable qualité éditoriale.

L’installation de ces modules s’opère via les dépôts officiels du Comprehensive R Archive Network (CRAN). En session de console, l’exécution des commandes suivantes garantit l’acquisition des versions stables les plus récentes :

install.packages(c("lmtest", "sandwich", "car", "tidyverse", "nlme", "orcutt"))

Une fois l’installation menée à son terme sans message d’erreur de compilation, le chargement explicite des environnements de travail au sein de votre script analytique s’effectue au moyen de l’instruction canonique library() :

library(lmtest)
library(sandwich)
library(car)
library(tidyverse)
library(nlme)
library(orcutt)

4.2 Vérification des versions et reproductibilité de la session R

La crédibilité de toute investigation quantitative repose sur le principe de reproductibilité exacte du pipeline de calcul. Dans le cadre de l’évaluation diagnostique des modèles linéaires sous R, des variations subtiles d’implémentation algorithmique ou des changements dans les tolérances numériques entre différentes versions d’un même package peuvent occasionnellement altérer de manière marginale les décimales des statistiques de test ou des p-valeurs. Il apparaît donc impératif de consigner scrupuleusement la configuration matérielle et logicielle active en appelant la fonction native suivante au terme de votre script :

sessionInfo()

Cette instruction consigne l’architecture du système d’exploitation, la version précise du moteur d’exécution de R ainsi que le versionnement exhaustif de chaque paquetage attaché et de ses dépendances sous-jacentes. Dans les contextes où des simulations stochastiques sont mobilisées pour évaluer la robustesse ou illustrer le comportement des estimateurs en présence d’erreurs autocorrélées, la reproductibilité impose de figer le générateur de nombres pseudo-aléatoires de R. L’assignation d’une graine déterministe au moyen de l’instruction set.seed() évite toute dérive aléatoire lors des générations successives :

set.seed(42)

Enfin, une vigilance particulière doit être accordée aux collisions d’espaces de noms de fonctions. Lorsque de multiples bibliothèques sont chargées concomitamment, des fonctions homonymes peuvent s’écraser mutuellement. Bien que la fonction bgtest() soit spécifique à lmtest, des fonctions diagnostiques générales telles que filter() ou select() font l’objet de conflits récurrents entre stats et dplyr. Pour garantir une exécution déterministe et pérenne de vos scripts au sein d’environnements collaboratifs, il est fortement conseillé d’invoquer les fonctions sensibles en préfixant explicitement le nom du package détenteur via l’opérateur double deux-points, par exemple lmtest::bgtest() ou sandwich::NeweyWest().

5. Génération et structuration d’un jeu de données psychométriques dans R

5.1 Simulation d’une étude de suivi comportemental et cognitif

Pour illustrer de façon tangible et didactique la mise en œuvre pratique, l’interprétation inférentielle et les répercussions diagnostiques du test de Breusch-Godfrey, nous allons simuler un protocole de suivi comportemental et cognitif longitudinal. Imaginons une étude expérimentale dans laquelle un échantillon d’individus ou une cohorte est suivi quotidiennement pendant N = 200 jours consécutifs. Le critère d’intérêt clinique Y est un score continu mesurant le niveau de stress psychologique perçu (mesuré sur une échelle étalonnée de 0 à 100). Ce score est théoriquement prédit par deux variables comportementales journalières majeures : l’effort cognitif consenti au travail (exprimé en unités d’heures pondérées par l’intensité) et la durée effective du sommeil réparateur nocturne (exprimée en heures décimales).

Dans un contexte réel, les résidus d’une telle régression ne sont presque jamais indépendants d’un jour à l’autre : le stress résiduel non expliqué par le sommeil ou le travail d’une journée a tendance à déborder sur le lendemain et le surlendemain, configurant une dépendance dynamique intrinsèque. Nous allons introduire de manière délibérée et parfaitement contrôlée une structure stochastique autorégressive d’ordre deux — un processus AR(2) — au sein des perturbations résiduelles. La relation génératrice sous-jacente s’établit selon la formalisation :

Stressₜ = 50 + 2.5 × Effortₜ – 3.8 × Sommeilₜ + εₜ
εₜ = 0.45 εₜ₋₁ + 0.25 εₜ₋₂ + uₜ

où uₜ est un bruit blanc gaussien d’espérance nulle et d’écart-type σ = 3. La traduction de ce processus stochastique sous R s’opère de manière directe et transparente au moyen du code scripté ci-dessous :

set.seed(12345)
n_obs <- 200
temps <- 1:n_obs
effort_cognitif <- round(rnorm(n_obs, mean = 7.5, sd = 1.8), 2)
sommeil_heures <- round(rnorm(n_obs, mean = 6.8, sd = 1.2), 2)

# Generation du bruit blanc pur
u <- rnorm(n_obs, mean = 0, sd = 3)
epsilon <- numeric(n_obs)

# Generation du processus autorégressif AR(2) pour les perturbations
for (t in 3:n_obs) {
  epsilon[t] <- 0.45 * epsilon[t - 1] + 0.25 * epsilon[t - 2] + u[t]
}

# Calcul de la variable critère de stress perçu
stress_percu <- round(50 + 2.5 * effort_cognitif - 3.8 * sommeil_heures + epsilon, 2)

# Assemblage dans une structure de donnees tibble
donnees_psy <- tibble(
  Jour = temps,
  Stress = stress_percu,
  Effort = effort_cognitif,
  Sommeil = sommeil_heures
)

5.2 Exploration descriptive et inspection visuelle préliminaire

Avant d’ajuster formellement le modèle de régression linéaire et de procéder aux tests d’hypothèses résiduelles, une exploration descriptive exhaustive des séries s’avère indispensable pour appréhender la dispersion empirique des variables et détecter d’éventuelles anomalies d’acquisition ou valeurs aberrantes. La distribution marginale de chaque variable psychométrique peut être synthétisée au moyen des fonctions descriptives conventionnelles :

summary(donnees_psy %>% select(Stress, Effort, Sommeil))

Cette commande permet de vérifier que les moyennes, médianes et quartiles observés s’inscrivent parfaitement dans les étendues théoriques attendues pour ces dimensions psychologiques. Toutefois, lorsque des données sont ordonnées de manière séquentielle le long d’un axe chronologique, l’inspection visuelle des trajectoires temporelles constitue l’étape exploratoire la plus féconde pour pressentir la présence d’une inertie stochastique ou de dynamiques cycliques sous-jacentes. À cette fin, la bibliothèque ggplot2 offre les fonctionnalités requises pour tracer l’évolution temporelle conjointe des mesures :

ggplot(donnees_psy, aes(x = Jour, y = Stress)) +
  geom_line(color = "#2c3e50", linewidth = 0.8) +
  geom_point(color = "#e74c3c", size = 1.5, alpha = 0.7) +
  theme_minimal() +
  labs(
    title = "Trajectoire temporelle du stress psychologique perçu",
    subtitle = "Suivi longitudinal quotidien sur 200 jours d'observation",
    x = "Journées d'observation",
    y = "Score de stress perçu (0-100)"
  )

L’observation attentive du graphique produit met immédiatement en évidence des ondulations graduelles plutôt qu’une dispersion purement erratique d’un jour à l’autre. Des séquences de jours consécutifs demeurent localisées au-dessus de la moyenne tendancielle, suivies de phases où les scores plongent de concert sous cette médiane. Cette persistance visuelle macroscopique constitue la signature empirique caractéristique d’un phénomène dynamique doté d’une mémoire temporelle, préfigurant l’existence d’une corrélation sérielle résiduelle qui nécessitera une démonstration inférentielle rigoureuse par le test de Breusch-Godfrey.

5.3 Ajustement du modèle de régression linéaire initial

L’étape analytique suivante consiste à ajuster le modèle linéaire standard sous l’hypothèse (naïve) que les résidus respectent scrupuleusement les conditions de Gauss-Markov d’indépendance mutuelle. Sous R, cet ajustement s’opère par l’appel classique à la fonction d’estimation par les moindres carrés ordinaires lm() :

modele_lineaire <- lm(Stress ~ Effort + Sommeil, data = donnees_psy)
summary(modele_lineaire)

L’inspection du tableau des coefficients de régression délivre des informations conformes à la structure génératrice théorique : la constante se situe à proximité immédiate de 50, l’effort cognitif affiche un coefficient positif fortement significatif (avoisinant +2.5), et la durée de sommeil se traduit par un coefficient négatif majeur (avoisinant -3.8), attestant de son effet protecteur sur le stress perçu. Le coefficient de détermination global R² indique une part substantielle de variance expliquée, conférant à première vue une solide légitimité scientifique au modèle ajusté.

Néanmoins, se limiter à l’examen de cette sortie console conventionnelle relève d’une négligence méthodologique coupable. Pour entamer la phase de validation diagnostique, il convient d’extraire les résidus empiriques estimés du modèle — calculés par e = Y – Ŷ — et de les intégrer directement au tableau de données initial pour permettre des investigations graphiques et formelles ciblées :

donnees_psy$residus_mco <- residuals(modele_lineaire)

Une représentation graphique de ces résidus en fonction de l’ordre tempologique des observations (l’axe des journées) révèle déjà une tendance marquée des résidus positifs à être suivis par d’autres résidus positifs, et inversement pour les résidus négatifs. Ce constat visuel impose de dépasser l’impression subjective pour mobiliser la puissance calculatoire et inférentielle de la fonction bgtest().

6. Exécution pas à pas du test de Breusch-Godfrey avec lmtest

6.1 Syntaxe fondamentale et arguments de la fonction bgtest

L’implémentation opérationnelle du test de Breusch-Godfrey au sein du langage R est encapsulée dans la fonction bgtest() mise à disposition par le package lmtest. La conception architecturale de cette fonction a été pensée pour conjuguer simplicité d’invocation et flexibilité paramétrique maximale. La signature formelle de la fonction, accessible via la documentation intégrée ?bgtest, se déploie selon les paramètres fondamentaux suivants :

bgtest(formula, order = 1, order.by = NULL, type = c("Chisq", "F"), data = list(), fill = 0)

Une maîtrise approfondie de chacun de ces arguments est indispensable pour adapter rigoureusement le test à vos contraintes de recherche empirique :

• formula : Cet argument accepte soit une formule symbolique de modèle de régression (ex. Stress ~ Effort + Sommeil), soit directement un objet de classe "lm" résultant d’un ajustement préalable via la fonction lm(). Dans la pratique, transmettre directement l’objet modele_lineaire constitue l’approche la plus robuste et la plus élégante, assurant la parfaite transmission de l’ensemble des métadonnées du modèle d’origine.
• order : Entier scalaire spécifiant l’ordre maximal de retard p à examiner au sein de la régression auxiliaire. Par défaut, l’argument est initialisé à order = 1, confinant l’investigation au diagnostic d’une dépendance de premier ordre AR(1). L’analyste doit modifier explicitement cette valeur pour investiguer des ordres autorégressifs supérieurs (ex. 2, 4, 7 ou plus).
• type : Chaîne de caractères indiquant la nature de la statistique d’inférence calculée et la loi théorique de référence utilisée pour dériver la p-valeur décisionnelle. Deux options sont configurables : "Chisq" (valeur retenue par défaut), qui calcule la statistique asymptotique du multiplicateur de Lagrange n × R²ₑ et la compare à une distribution du Chi-deux à p degrés de liberté ; et "F", qui calcule la statistique de Fisher classique de significativité globale au sein de la régression auxiliaire et se réfère à une distribution F(p, n – k – p).
• order.by : Vecteur optionnel définissant la variable selon laquelle les observations doivent être préalablement triées pour structurer la dépendance temporelle. Si vos données sont déjà agencées dans un ordre chronologique strict au sein du data.frame (comme c’est le cas dans notre exemple avec la variable Jour), cet argument peut être omis. S’il existe une variable de date ou d’horodatage non ordonnée, désigner order.by = ~ Jour assure le tri préalable indispensable.
• fill : Valeur numérique utilisée pour imputer les p premiers retards manquants des résidus dans la régression auxiliaire. La valeur par défaut fill = 0 s’aligne fidèlement sur la dérivation théorique canonique de Godfrey, remplaçant les valeurs inobservables antérieures au début de la série par l’espérance mathématique des perturbations sous H₀.

6.2 Application concrète à un retard d’ordre 1

Pour débuter le diagnostic formel sur notre modèle psychométrique simulé, nous amorçons la démarche en interrogeant la présence d’une autocorrélation résiduelle d’ordre un, ce qui revient à tester l’hypothèse H₀ : ρ₁ = 0 contre H₁ : ρ₁ ≠ 0. L’exécution sous R requiert une instruction concise et directe :

test_bg_ordre1 <- bgtest(modele_lineaire, order = 1, type = "Chisq")
print(test_bg_ordre1)

L’affichage généré en console R délivre la structure d’inférence standardisée suivante :

Breusch-Godfrey test for higher-order serial correlation
modele_lineaire
LM test = 58.421, df = 1, p-value = 2.115e-14

L’examen des résultats est univoque : la statistique empirique du multiplicateur de Lagrange s’élève à LM = 58.421 pour un unique degré de liberté. La probabilité critique associée (p-valeur) s’établit à une valeur infinitésimale de l’ordre de 2.115 × 10⁻¹⁴, se situant de façon écrasante sous le seuil d’inférence conventionnel α = 0.05. Le chercheur rejette donc formellement l’hypothèse nulle d’indépendance sérielle des erreurs au premier ordre.

Si l’on compare ce résultat avec l’exécution parallèle du test classique de Durbin-Watson via la fonction homonyme dwtest(modele_lineaire), on constatera que Durbin-Watson rejette également l’indépendance avec une statistique d très inférieure à 2 (avoisinant 1.02, confirmant une autocorrélation positive majeure). L’extraction sélective et programmée des composantes de l’objet produit par bgtest peut s’opérer directement à travers les attributs de la liste :

statistique_lm <- test_bg_ordre1$statistic
p_valeur <- test_bg_ordre1$p.value
degres_liberte <- test_bg_ordre1$parameter

Cette modularité s’avère particulièrement utile pour injecter dynamiquement les résultats numériques au sein de rapports automatisés rédigés via des documents computationnels R Markdown ou Quarto.

6.3 Extension de l’évaluation aux ordres supérieurs de dépendance

Bien que le test d’ordre 1 ait mis en évidence un problème évident de corrélation sérielle, s’arrêter à ce premier niveau d’analyse serait méthodologiquement myope. La véritable puissance du test de Breusch-Godfrey réside dans son aptitude à sonder la persistance temporelle sur des horizons plus lointains. Sachant que notre jeu de données expérimental intègre une dynamique génératrice de type AR(2), il est primordial de spécifier order = 2 pour tester formellement la nullité conjointe H₀ : ρ₁ = ρ₂ = 0 :

test_bg_ordre2 <- bgtest(modele_lineaire, order = 2, type = "Chisq")
print(test_bg_ordre2)

La console renvoie alors une statistique LM = 74.893 associée à df = 2 et une p-valeur encore plus infime (inférieure à 10⁻¹⁶). Ce saut significatif de la valeur de la statistique de test entre l’ordre 1 et l’ordre 2 apporte la preuve statistique formelle que le résidu retardé de deux périodes eₜ₋₂ contribue substantiellement et de façon incrémentale à expliquer la variance du résidu contemporain eₜ, au-delà de l’effet capté par le premier retard seul.

Dans de nombreux contextes cliniques ou organisationnels, les cycles de recueil imposent d’investiguer des ordres encore plus élevés, par exemple p = 4 (dynamique trimestrielle ou cycles intra-quotidiens à 4 prises) ou p = 7 (rythmicité hebdomadaire). Plutôt que de lancer des commandes disparates de manière itérative et désorganisée, la démarche rigoureuse consiste à programmer une boucle fonctionnelle via la bibliothèque purrr issue du tidyverse afin d’évaluer le comportement de la statistique de test sur une séquence ordonnée de retards allant de 1 jusqu’à 10 :

tableau_retards <- tibble(Ordre = 1:10) %>%
  mutate(
    Test_Obj = map(Ordre, ~ bgtest(modele_lineaire, order = .x, type = "Chisq")),
    Statistique_LM = map_dbl(Test_Obj, ~ .x$statistic),
    Degres_Liberte = map_int(Test_Obj, ~ as.integer(.x$parameter)),
    P_Valeur = map_dbl(Test_Obj, ~ .x$p.value)
  ) %>%
  select(-Test_Obj)

print(tableau_retards)

La génération de ce tableau récapitulatif permet à l’analyste d’observer l’évolution du LM au fur et à mesure que les degrés de liberté augmentent. Tant que l’ajout d’un ordre de retard apporte une information dynamique réelle non capturée par les retards précédents, la statistique LM progresse de façon substantielle. Lorsqu’on franchit l’ordre générateur réel (ici l’ordre 2), la statistique tend à plafonner, n’augmentant que par des gains marginaux aléatoires conformes à la fluctuation stochastique, ce qui offre un indicateur précieux pour le paramétrage des modèles correctifs ultérieurs.

7. Stratégies méthodologiques pour déterminer l’ordre p des retards

7.1 Utilisation de la fonction d’autocorrélation (ACF) et d’autocorrélation partielle (PACF)

L’un des dilemmes méthodologiques les plus fréquents posés par l’application du test de Breusch-Godfrey réside dans la sélection de l’ordre de retard maximal p. Choisir un ordre p arbitrairement bas (par exemple systématiquement p = 1 par automatisme d’esprit) expose l’analyste à ignorer des mémoires sérielles d’ordre supérieur et à conclure à tort à la validité d’un modèle dont les erreurs sont en réalité dépendantes. À l’inverse, sélectionner un ordre p disproportionné (par exemple p = 20 sur un échantillon de 80 observations) grève lourdement la puissance statistique du test et dégrade l’estimation en consommant inutilement des degrés de liberté. Pour guider ce choix sur des bases objectives, l’analyse graphique exploratoire conjointe de la fonction d’autocorrélation (ACF) et de la fonction d’autocorrélation partielle (PACF) des résidus du modèle initial constitue la méthode heuristique de référence.

Au sein de R, l’extraction et la visualisation des corrélogrammes s’effectuent sans difficulté au moyen des fonctions intégrées acf() et pacf() appliquées directement au vecteur des résidus extraits du modèle linéaire ajusté :

par(mfrow = c(1, 2))
acf(residuals(modele_lineaire), main = "Fonction d'Autocorrélation (ACF)")
pacf(residuals(modele_lineaire), main = "Autocorrélation Partielle (PACF)")
par(mfrow = c(1, 1))

L’interprétation de ces deux graphiques repose sur des fondements théoriques bien établis dans la littérature des séries temporelles de Box et Jenkins. Les lignes en pointillés bleus représentent les bandes de confiance asymptotiques à 95 % établies par Bartlett, calculées approximativement par ± 1.96 / √n. Tout pic franchissant ces limites critiques est réputé statistiquement non nul avec un risque d’erreur de 5 %. La signature graphique permet de discriminer la nature du processus générateur :

• Profil d’un processus AR(p) : L’ACF présente une décroissance progressive, géométrique ou sinusoïdale amortie vers zéro, tandis que la PACF s’interrompt brutalement (décroissance nette) après le retard p. Si la PACF présente des pics significatifs aux retards 1 et 2, puis rentre intégralement à l’intérieur du corridor de Bartlett pour tous les ordres ultérieurs (3, 4, …), le chercheur détient une indication formelle que la structure autorégressive sous-jacente est d’ordre p = 2.
• Profil d’un processus MA(q) : Le schéma s’inverse : l’ACF s’annule nettement après le retard q, tandis que la PACF s’éteint de manière lente et continue.

Dans notre cas de simulation psychométrique, l’examen de la PACF révèle précisément deux pics fortement saillants au-delà des bandes de Bartlett aux ordres 1 et 2, suivis d’un affaissement immédiat sous le seuil critique pour l’ensemble des retards supérieurs. Cette observation visuelle fournit une justification empirique irréfutable pour configurer le test de Breusch-Godfrey avec l’argument formel order = 2.

7.2 Recours aux critères d’information théoriques (AIC et BIC)

Bien que l’inspection des corrélogrammes ACF et PACF apporte une aide visuelle précieuse, elle peut demeurer sujette à des désaccords d’interprétation lorsque les séries réelles sont entachées d’un bruit stochastique hétérogène ou de valeurs extrêmes. Pour consolider la sélection de l’ordre de retard maximal sur un socle quantitatif rigoureux et mathématiquement traçable, il est hautement recommandé de recourir aux critères d’information fondés sur la théorie de l’information et la parcimonie statistique, principalement le critère d’information d’Akaike (AIC) et le critère d’information bayésien de Schwarz (BIC).

L’approche consiste à ajuster formellement la séquence des régressions auxiliaires associées au test de Breusch-Godfrey pour une série d’ordres croissants allant de p = 1 à p = p_max, puis à calculer et comparer la valeur de ces critères pour chaque spécification auxiliaire. Mathématiquement, l’AIC et le BIC sanctionnent la déviance résiduelle (la somme des résidus au carré) par une pénalité dépendant directement du nombre de paramètres estimés dans l’équation :

AIC = 2k – 2 ln(L̂)
BIC = k ln(n) – 2 ln(L̂)

où k représente le nombre de paramètres de régression auxiliaire (incluant les régresseurs initiaux et les retards résiduels) et L̂ désigne la vraisemblance maximale estimée. Puisque le BIC applique une pénalisation logarithmique ln(n) proportionnelle à la taille d’échantillon, il tend à pénaliser plus sévèrement la surparamétrisation que l’AIC, privilégiant des structures autorégressives plus parcimonieuses.

Nous pouvons automatiser l’estimation des modèles autorégressifs sur les résidus et extraire leurs profils comparatifs d’information au moyen du script R suivant :

residus_vect <- residuals(modele_lineaire)

selection_ordres <- tibble(Ordre = 1:8) %>%
  mutate(
    Modele_AR = map(Ordre, ~ ar(residus_vect, aic = FALSE, order.max = .x, method = "ols")),
    AIC_Val = map_dbl(Ordre, function(p) {
      mod <- lm(residus_vect[-(1:p)] ~ embed(residus_vect, p + 1)[, -1])
      AIC(mod)
    }),
    BIC_Val = map_dbl(Ordre, function(p) {
      mod <- lm(residus_vect[-(1:p)] ~ embed(residus_vect, p + 1)[, -1])
      BIC(mod)
    })
  )

print(selection_ordres %>% select(Ordre, AIC_Val, BIC_Val))

Le principe de décision statistique dicte de retenir l’ordre de retard qui minimise de façon absolue la valeur du critère retenu (le point le plus bas de la courbe d’information). Dans notre cadre de simulation empirique, tant l’AIC que le BIC atteignent leur minimum absolu pour l’ordre p = 2, confirmant de manière concordante les déductions issues de la PACF et légitimant l’application du test de Breusch-Godfrey paramétré sur ce niveau optimal de complexité stochastique.

7.3 Considérations substantielles liées au protocole expérimental

En dépit de l’attrait mathématique que procurent les corrélogrammes et les critères d’information pénalisés, la sélection de l’ordre de retard ne saurait s’abstraire d’une réflexion substantielle approfondie ancrée dans l’architecture même du devis expérimental et dans les propriétés théoriques des construits mesurés. Les algorithmes d’optimisation purement empiriques demeurent sensibles aux aléas d’échantillonnage et peuvent suggérer des retards mathématiquement séduisants mais conceptuellement dépourvus de tout fondement clinique, biologique ou psychologique. L’analyste se doit d’aligner le choix de p sur la granularité temporelle de son protocole d’observation.

Considérons le cas d’une recherche en psychologie clinique utilisant l’évaluation écologique momentanée (EMA) à raison de quatre notifications aléatoires par jour pour évaluer les ruminations anxieuses. Dans un tel schéma expérimental, spécifier un test de Breusch-Godfrey d’ordre p = 4 n’est pas un acte arbitraire : cela équivaut à tester la persistance résiduelle sur une journée complète de veille. Si l’évaluation est quotidienne, investiguer l’ordre p = 7 répond à l’hypothèse substantive légitime d’une récurrence hebdomadaire des mécanismes de stress au travail. À l’inverse, si l’étude porte sur des temps de réaction mesurés toutes les 500 millisecondes dans un paradigme d’attention soutenue de type Stroop, les processus d’inhibition cognitive ou d’interférence inter-essais se déploient typiquement sur les deux ou trois essais immédiatement consécutifs, justifiant de borner a priori l’investigation formelle à p = 3.

Une mise en garde épistémologique majeure doit être fermement formulée à l’encontre de la pratique pernicieuse du dragage de données (communément qualifié de p-hacking ou d’ajustement opportuniste des paramètres de test). Il est formellement inadmissible sur le plan déontologique d’exécuter le test de Breusch-Godfrey pour une multiplicité d’ordres allant de 1 à 15, pour ne rapporter dans le manuscrit final que le seul ordre dont la p-valeur s’avère non significative afin de revendiquer abusivement la conformité du modèle linéaire ordinaire aux postulats de Gauss-Markov. Cette manipulation vicie l’erreur de type I globale de la recherche. La rigueur scientifique impose de pré-spécifier les ordres de retards diagnostiqués en fonction de la périodicité théorique du protocole, ou d’appliquer des corrections pour tests multiples (telles que les ajustements de Bonferroni ou de Benjamini-Hochberg) si une exploration séquentielle systématique est déployée.

8. Interprétation statistique et inférentielle des sorties R

8.1 Décodage ligne par ligne de la sortie console de bgtest

L’interprétation correcte des résultats d’un test statistique impose une maîtrise parfaite du sens formel de chaque composante textuelle et numérique restituée par le moteur de calcul. Lorsque la fonction bgtest() est sollicitée, R imprime dans la console une sortie condensée mais d’une haute densité informative. Reprenons la sortie analytique obtenue lors de l’application du test à l’ordre p = 2 sur notre jeu de données :

Breusch-Godfrey test for higher-order serial correlation
modele_lineaire
LM test = 74.893, df = 2, p-value < 2.2e-16

Le décodage analytique rigoureux de cette restitution s’articule selon quatre axes fondamentaux :

• L’intitulé du test (« Breusch-Godfrey test for higher-order serial correlation ») : Confirme que le moteur a mobilisé la formulation auxiliaire de Godfrey pour sonder des corrélations sérielles d’ordre supérieur, englobant potentiellement les représentations AR(p) et MA(p).
• La désignation des données (« data: modele_lineaire ») : Rappelle l’objet de régression parent soumis à l’évaluation, attestant que les résidus et les variables prédictives ont été fidèlement extraits de l’environnement dudit modèle.
• La statistique d’évaluation (« LM test = 74.893 ») : Exprime la valeur numérique réalisée de la statistique du multiplicateur de Lagrange, calculée par l’égalité n × R²ₑ issue de la régression auxiliaire. Dans cet exemple, le produit de la taille de l’échantillon par le R² auxiliaire atteint une amplitude considérable de près de 75 unités, reflétant une forte capacité des résidus passés à prédire les déviations contemporaines.
• Les degrés de liberté (« df = 2 ») : Indiquent le nombre de contraintes linéaires conjointement testées sous l’hypothèse nulle. Ce chiffre coïncide obligatoirement avec la valeur assignée au paramètre order (ici 2), signifiant que la loi de référence sous-jacente est une distribution du Chi-deux à deux degrés de liberté.
• La probabilité critique (« p-value < 2.2e-16 ») : Représente l’intégrale supérieure de la densité du Chi-deux au-delà de la valeur observée 74.893. La valeur affichée correspond à la borne inférieure de précision numérique en virgule flottante native de R. Puisque cette probabilité est virtuellement nulle, le rejet de l’hypothèse nulle H₀ s’impose avec un niveau de certitude statistique absolu.

8.2 Comparaison entre la version asymptotique Chi-deux et la variante F de Fisher

Bien que l’invocation standard de bgtest() opère par défaut via la statistique du multiplicateur de Lagrange adossée à la distribution asymptotique du Chi-deux (type = "Chisq"), il est impératif pour l’analyste rigoureux de comprendre les limites de cette approximation en présence d’échantillons de taille modeste. La statistique asymptotique n × R²ₑ suppose formellement que la taille d’échantillon n tend vers l’infini. Lorsque n est relativement faible (par exemple n = 40 ou n = 60) et que le nombre de régresseurs k et de retards p est proportionnellement élevé, le multiplicateur n surestime la compacité de la distribution, provoquant des p-valeurs artificiellement trop petites et augmentant le risque de rejeter H₀ à tort.

Pour contrer cette distorsion, l’argument type = "F" contraint la fonction à dériver la statistique de Fisher exacte de significativité globale des coefficients des résidus retardés au sein de la régression auxiliaire. Pour exécuter cette variante hautement recommandée en petit échantillon sous R, la syntaxe s’établit ainsi :

test_bg_fisher <- bgtest(modele_lineaire, order = 2, type = "F")
print(test_bg_fisher)

La console renvoie alors une structure modifiée :

Breusch-Godfrey test for higher-order serial correlation
modele_lineaire
F = 58.742, df1 = 2, df2 = 195, p-value < 2.2e-16

Dans cette sortie, la statistique d’intérêt devient une valeur de F calculée avec df1 = p = 2 degrés de liberté au numérateur et df2 = n – k – p = 200 – 3 – 2 = 195 degrés de liberté au dénominateur. L’emploi de cette statistique F offre une robustesse inférentielle supérieure car elle intègre explicitement au dénominateur le terme de pénalisation de la variance non expliquée normalisé par les degrés de liberté résiduels réels. Dans le cadre de notre exemple où n = 200, l’échantillon est suffisamment vaste pour que la décision statistique demeure identique (rejet massif dans les deux cas) ; néanmoins, dans les protocoles de recherche contraints par de petites cohortes, l’arbitrage en faveur de la variante type = "F" constitue une pratique méthodologique d’une prudence exemplaire.

8.3 Formulation formelle des conclusions selon les normes académiques

La restitution des diagnostics statistiques au sein des publications savantes, mémoires de thèse ou revues à comité de lecture doit répondre à des standards formels stricts garantissant la transparence méthodologique, à l’instar des normes édictées par l’American Psychological Association (APA, 7ᵉ édition). Il est rigoureusement insuffisant d’écrire de façon vague : « les résidus du modèle étaient autocorrélés ». Le chercheur se doit de consigner formellement le nom du test mobilisé, l’ordre de retard spécifique investigué et sa justification substantive, la valeur réalisée de la statistique d’épreuve, ses degrés de liberté précis et la p-valeur exacte.

Dans l’éventualité où le test de Breusch-Godfrey confirme la préservation de l’hypothèse nulle (absence d’autocorrélation), une formulation académique exemplaire prendra la forme suivante :

« L’indépendance conditionnelle des perturbations stochastiques a été vérifiée au moyen du test de Breusch-Godfrey paramétré jusqu’à un ordre de retard p = 2, représentatif de l’inertie potentielle sur 48 heures. Les résultats indiquent l’absence d’autocorrélation sérielle résiduelle statistiquement significative, que ce soit sous la distribution asymptotique du multiplicateur de Lagrange, χ²(2) = 1.84, p = .398, ou sous sa variante de Fisher adaptée aux petits échantillons, F(2, 195) = 0.91, p = .404. Par conséquent, les postulats d’indépendance de Gauss-Markov sont réputés satisfaits, validant l’efficacité des estimateurs des moindres carrés ordinaires et la précision des erreurs-types associées. »

À l’inverse, dans notre situation expérimentale où l’hypothèse nulle est irrévocablement rejetée, le compte-rendu devra stipuler avec netteté la violation constatée et annoncer sans ambiguïté les démarches correctives adoptées :

« Le diagnostic de corrélation sérielle conduit à l’aide du test de Breusch-Godfrey révèle une violation majeure de l’hypothèse d’indépendance des résidus au second ordre, χ²(2) = 74.89, p < .001 (F(2, 195) = 58.74, p < .001). Ce rejet hautement significatif confirme la présence d’une mémoire temporelle substantielle au sein des termes d’erreur non modélisés. En conséquence, les erreurs-types estimées par les MCO conventionnels s’avèrent structurellement biaisées vers le bas, interdisant toute inférence directe fondée sur les statistiques t usuelles. Une procédure d’estimation robuste hétéroscédastique et autocorrélée (HAC) ou une spécification dynamique autorégressive a été impérativement substituée au modèle initial pour garantir la validité des conclusions statistiques. »

9. Conséquences empiriques de la violation de l’hypothèse d’indépendance

9.1 Biais dans l’estimation de la variance des estimateurs MCO

L’impact mécanique le plus dévastateur de la persistance d’une corrélation sérielle non corrigée se concentre sur l’estimation de la matrice de variance-covariance des paramètres estimés β̂. Pour appréhender l’origine mathématique de cette défaillance, rappelons que l’estimateur de Gauss-Markov standard s’énonce sous la forme matricielle β̂ = β + (X’X)⁻¹X’ε. Dès lors que l’on calcule la variance théorique de cet estimateur conditionnellement à la matrice des prédicteurs X, l’expression générale sans contrainte d’indépendance prend la formulation universelle :

Var(β̂|X) = (X’X)⁻¹ X’ E(εε’|X) X (X’X)⁻¹ = (X’X)⁻¹ X’ Ω X (X’X)⁻¹

Dans l’hypothèse idéale d’indépendance et d’homoscédasticité, la matrice centrale se simplifie de manière spectaculaire : puisque Ω = σ²Iₙ, le terme scalaire σ² peut être extrait de l’expression, ce qui permet la simplification matricielle fondamentale : Var(β̂|X) = σ² (X’X)⁻¹ (X’X) (X’X)⁻¹ = σ² (X’X)⁻¹. Les logiciels d’économétrie calculent sans exception cette forme simplifiée en substituant à σ² l’estimateur sans biais de la variance résiduelle s² = e’e / (n – k).

Or, en présence d’une autocorrélation résiduelle positive — qui représente l’écrasante majorité des cas empiriques en sciences humaines et biologiques —, les éléments situés hors de la diagonale de la matrice Ω ne sont plus nuls mais positifs. Si, de surcroît, les régresseurs X manifestent eux-mêmes une inertie temporelle positive (ce qui est systématiquement le cas pour des variables continues comme le sommeil, l’effort cognitif ou la température corporelle qui varient lentement au cours du temps), le produit matriciel X’ Ω X engendre une somme de termes croisés positifs d’une amplitude substantiellement supérieure au produit scalaire σ² (X’X). En conséquence immédiate :

s² (X’X)⁻¹ ≪ (X’X)⁻¹ X’ Ω X (X’X)⁻¹

La formule standard implémentée dans la fonction lm() sous-estime drastiquement la véritable variance des coefficients de régression. Les erreurs-types obtenues en prenant la racine carrée de la diagonale de cette matrice se révèlent trompeusement compressées, masquant l’extrême instabilité d’échantillonnage de l’estimateur des moindres carrés ordinaires. Ce biais d’estimation de la variance résiduelle se propage sans atténuation à l’ensemble des inférences statistiques subséquentes, sapant irrémédiablement la légitimité des tests d’hypothèses.

9.2 Risque d’erreurs inférentielles de type I et conclusions trompeuses

La conséquence directe de la sous-estimation chronique des erreurs-types réside dans le gonflement artificiel des statistiques de test t de Student et F de Fisher. Rappelons que la statistique t associée à un régresseur donné j pour tester l’hypothèse d’absence d’effet H₀ : βⱼ = 0 se formalise par le rapport canonique :

t_obs = β̂ⱼ / SE(β̂ⱼ)

Lorsque le dénominateur SE(β̂ⱼ) est artificiellement réduit par le biais d’autocorrélation, la valeur absolue de la statistique t_obs se trouve propulsée vers des amplitudes extrêmes. Ce phénomène génère une prolifération massive d’erreurs de type I : le chercheur est conduit à observer des p-valeurs inférieures aux seuils symboliques de .05, .01 voire .001, attribuant à tort un rôle explicatif ou thérapeutique majeur à des prédicteurs qui, dans la population parente, n’exercent aucune influence sur le critère. La précision apparente des estimations n’est qu’une illusion statistique induite par le dénombrement erroné de l’information indépendante effective contenue dans l’échantillon.

Sur le plan conceptuel, lorsqu’une série présente une autocorrélation positive, chaque nouvelle observation yₜ n’apporte pas une quantité d’information totalement nouvelle par rapport à l’observation précédente yₜ₋₁. La taille effective de l’échantillon d’observations indépendantes, notée n_eff, est substantiellement inférieure à la taille physique apparente n. Le calcul conventionnel traite un jeu de données de 200 mesures autocorrélées comme s’il s’agissait de 200 points d’information stochastiquement autonomes, alors que le contenu informationnel réel équivaut parfois à seulement 40 ou 50 observations indépendantes. Cette illusion de puissance conduit à des conclusions péremptoires sur l’efficacité d’interventions cognitives ou pharmacologiques, contribuant de manière directe à la « crise de la réplicabilité » qui frappe les sciences empiriques contemporaines.

Enfin, cette défaillance inférentielle s’accompagne d’une dégradation marquée des performances de généralisation du modèle. Un modèle calibré sous l’emprise d’une forte autocorrélation résiduelle non corrigée manifeste un surajustement (overfitting) sévère au bruit stochastique mémoriel de l’échantillon d’apprentissage. Dès lors que ce modèle est déployé dans un cadre prédictif prospectif ou confronté à de nouvelles cohortes hors-échantillon (out-of-sample validation), ses capacités de prédiction s’effondrent spectaculairement. L’erreur quadratique moyenne de prédiction (RMSE) explose, démontrant a posteriori que la précision affichée dans la publication princeps reposait sur des artefacts méthodologiques.

9.3 Impact sur le coefficient de détermination R-carré

Une autre répercussion pernicieuse de l’autocorrélation sérielle concerne la distorsion systématique du coefficient de détermination R² (et de sa version ajustée R²_adj). Dans l’estimation par les MCO, le R² est formulé comme le ratio de la somme des carrés expliqués par la régression sur la somme totale des carrés des déviations : R² = 1 – (e’e / y’y). Lorsque le terme d’erreur suit un processus autorégressif positif, les perturbations stochastiques successives ont tendance à persister du même côté du plan de régression, réduisant artificiellement la variance d’échantillonnage de court terme des résidus observés eₜ par rapport à leur véritable dispersion asymptotique.

Cette contraction locale de la somme des carrés résiduels induit une surestimation mécanique du R². Le modélisateur se voit gratifié d’un coefficient de détermination flatteur, oscillant par exemple entre 0.70 et 0.85, lui laissant supposer qu’il est parvenu à capturer les déterminants fondamentaux du comportement humain étudié. En réalité, une part substantielle de ce pouvoir explicatif apparent ne provient absolument pas de la pertinence théorique des variables indépendantes sélectionnées, mais découle de la structure autorégressive endogène du bruit de mesure non purgé. Ce phénomène s’apparente conceptuellement aux célèbres « régressions fallacieuses » (spurious regressions) mises en lumière par Granger et Newbold dans le domaine des séries temporelles non stationnaires.

De surcroît, cette distorsion fausse l’évaluation de la hiérarchie d’importance relative des différents régresseurs. Lorsque plusieurs prédicteurs affichent des degrés d’autocorrélation interne inégaux, le modèle MCO biaise asymétriquement les statistiques de significativité au profit des variables qui partagent la plus forte mémoire temporelle avec le bruit résiduel. Le chercheur risque ainsi de reléguer au second plan des régresseurs causalement décisifs mais à variation rapide, pour privilégier indûment des variables lentes dont la covariance apparente avec le critère est largement amplifiée par la dynamique sérielle commune non contrôlée.

10. Solutions correctives et approches de remédiation sous R

10.1 Estimation d’erreurs-types robustes à l’hétéroscédasticité et à l’autocorrélation (HAC)

Lorsque le test de Breusch-Godfrey a révélé de manière irréfutable l’existence d’une corrélation sérielle résiduelle, la stratégie de remédiation la plus immédiate et la plus élégante — dès lors que le modèle de régression ne comporte pas de variables dépendantes retardées et que la forme fonctionnelle est théoriquement réputée correcte — consiste à préserver les coefficients de régression estimés par les MCO tout en corrigeant la matrice de variance-covariance au moyen d’un estimateur robuste. Cette approche a été formalisée par Whitney K. Newey et Kenneth D. West en 1987 à travers l’estimateur HAC (Heteroskedasticity and Autocorrelation Consistent).

L’estimateur de Newey-West remplace la matrice centrale inconnue Ω de la formulation en sandwich par une somme pondérée non paramétrique des autocovariances empiriques des résidus jusqu’à un horizon de retard de troncature m. Pour garantir que la matrice de variance-covariance résultante demeure semi-définie positive, Newey et West appliquent un noyau de lissage triangulaire (noyau de Bartlett) qui assigne des poids linéairement décroissants aux autocovariances au fur et à mesure que l’intervalle temporel s’accroît : w(j, m) = 1 – [ j / (m + 1) ] pour j ≤ m. Sous l’environnement R, la combinaison fonctionnelle des bibliothèques sandwich et lmtest permet d’exécuter cette correction d’une manière limpide :

# Calcul de la matrice de covariance robuste HAC de Newey-West
vcov_hac <- sandwich::NeweyWest(modele_lineaire, lag = 2, prewhite = FALSE)

# Réévaluation des tests de significativité des paramètres
inf_robuste <- lmtest::coeftest(modele_lineaire, vcov = vcov_hac)
print(inf_robuste)

En comparant cette sortie avec celle de la commande initiale summary(modele_lineaire), l’analyste constatera que les coefficients β̂ demeurent rigoureusement identiques, mais que leurs erreurs-types ont augmenté de manière significative pour refléter la véritable incertitude d’échantillonnage induite par la mémoire sérielle. Les statistiques t voient leurs amplitudes réajustées à la baisse et les p-valeurs associées augmentent en conséquence. Cette remédiation offre une protection absolue contre l’inflation de l’erreur de type I tout en évitant d’avoir à reparamétrer la structure structurelle de l’équation originale.

10.2 Modélisation dynamique et intégration de retards temporels

Si la correction HAC de Newey-West s’avère irremplaçable lorsque la corrélation sérielle est assimilée à une simple nuisance stochastique résiduelle, elle élude la question de fond : l’autocorrélation des erreurs signale le plus souvent que le modèle statique ignore la véritable dynamique temporelle du processus étudié. En sciences du comportement, l’état psychologique d’un individu à l’instant t est intrinsèquement déterminé par son état antérieur à l’instant t-1. Plutôt que de subir passivement cette mémoire dans le terme d’erreur, la démarche méthodologique la plus rigoureuse consiste à internaliser cette dynamique au sein même de l’équation de régression en formulant un modèle autorégressif à retards échelonnés (Autoregressive Distributed Lag ou ARDL).

Dans notre cas d’étude sur le stress psychologique perçu, cette internalisation dynamique revient à introduire le score de stress retardé d’une ou deux journées parmi les variables prédictives contemporaines. La formulation mathématique révisée s’écrit alors :

Stressₜ = β₀ + γ₁ Stressₜ₋₁ + γ₂ Stressₜ₋₂ + β₁ Effortₜ + β₂ Sommeilₜ + νₜ

L’ajustement opérationnel d’un tel modèle dynamique sous R requiert la création préalable des variables retardées, ce qui s’effectue sans difficulté grâce à la fonction lag() de la bibliothèque dplyr :

donnees_dynamiques <- donnees_psy %>%
  mutate(
    Stress_Lag1 = lag(Stress, 1),
    Stress_Lag2 = lag(Stress, 2)
  )

modele_dynamique <- lm(Stress ~ Stress_Lag1 + Stress_Lag2 + Effort + Sommeil, data = donnees_dynamiques)
summary(modele_dynamique)

Une fois cette architecture dynamique estimée, l’analyste se doit d’exécuter à nouveau le test de Breusch-Godfrey sur ce nouveau modèle pour vérifier si l’intégration des variables endogènes retardées a réussi à purger l’intégralité de l’autocorrélation :

bgtest(modele_dynamique, order = 2, type = "F")

Dans cette nouvelle régression, le test de Breusch-Godfrey affiche une p-valeur largement supérieure à 0.05. L’hypothèse nulle d’indépendance ne peut plus être rejetée : le terme d’erreur νₜ s’est contracté en un bruit blanc pur. Cette approche dynamique offre une double plus-value scientifique : d’une part, elle restaure la conformité absolue aux postulats de Gauss-Markov, et d’autre part, elle enrichit la portée théorique de la recherche en dissociant l’effet direct contemporain des prédicteurs (impact immédiat du sommeil sur le stress) de leur effet cumulé de long terme médié par le paramètre d’inertie psychologique γ.

10.3 Méthodes de réestimation globale (Prais-Winsten et Cochrane-Orcutt)

Une approche classique concurrente, issue de la tradition économétrique historique mais toujours largement mobilisée lorsque le modèle théorique doit impérativement conserver une formulation statique (excluant formellement tout retard de Y parmi les régresseurs), consiste à réestimer l’ensemble des paramètres au moyen de méthodes de quasi-différenciation. Parmi ces procédures, la méthode itérative de Cochrane-Orcutt et sa généralisation supérieure conçue par Prais et Winsten occupent une place de choix.

Considérons le modèle affecté d’une erreur AR(1) : Yₜ = X’ₜβ + εₜ, avec εₜ = ρ εₜ₋₁ + uₜ. En multipliant l’équation retardée d’une période par le coefficient de persistance ρ et en la soustrayant de l’équation contemporaine, on obtient le système transformé :

(Yₜ – ρ Yₜ₋₁) = (Xₜ – ρ Xₜ₋₁)’β + uₜ

Puisque le nouveau terme de perturbation uₜ constitue un bruit blanc stochastique par construction, l’application des moindres carrés ordinaires à ces variables quasi-différenciées délivre des estimateurs linéaires sans biais à variance minimale (BLUE). La procédure itérative de Cochrane-Orcutt alterne l’estimation des résidus MCO, l’estimation du paramètre ρ̂, et la réestimation de β jusqu’à convergence numérique. Cependant, Cochrane-Orcutt supprime purement et simplement la première observation t = 1 du jeu de données car ses retards sont inobservables. La transformation de Prais-Winsten corrige cette déperdition d’information en préservant le premier point via la pondération spécifique : Y₁* = √(1 – ρ²) Y₁ et X₁* = √(1 – ρ²) X₁.

L’application de la procédure de Cochrane-Orcutt sous R s’exécute de façon automatisée via le package spécialisé orcutt :

modele_cochrane <- orcutt::cochrane.orcutt(modele_lineaire)
summary(modele_cochrane)

La sortie restitue le nombre d’itérations nécessaires à la stabilisation de l’algorithme, l’estimation finale du coefficient autorégressif ρ̂, la valeur de la statistique de Durbin-Watson post-transformation (qui converge typiquement vers 2.0), ainsi que le tableau révisé des coefficients de régression purgés de l’influence de l’autocorrélation. Cette réestimation globale garantit des erreurs-types asymptotiquement exactes tout en maintenant l’interprétation originale des coefficients d’impact statiques.

10.4 Passage aux modèles mixtes et aux structures de corrélation personnalisées

Lorsque les données empiriques s’inscrivent dans une architecture résolument hiérarchique, multiniveau ou longitudinale complexe (mesures répétées nichées au sein d’individus, déséquilibres dans les calendriers d’observation, valeurs manquantes non alignées), l’inflexibilité des transformations quasi-différenciées conventionnelles atteint ses limites opérationnelles. Dans ce paradigme avancé, le standard méthodologique contemporain impose de basculer vers les modèles linéaires généralisés par moindres carrés (GLS) ou vers les modèles linéaires à effets mixtes implémentés dans la bibliothèque de référence nlme (Nonlinear Mixed-Effects), développée par José Pinheiro et Douglas Bates.

La fonction gls() permet d’estimer simultanément les paramètres structurels de la régression et les hyperparamètres régissant la matrice de variance-covariance résiduelle par la méthode du maximum de vraisemblance restreint (REML). L’analyste peut spécifier explicitement une structure de corrélation temporelle continue ou discrète via l’argument correlation, en mobilisant des classes telles que corAR1() pour un processus du premier ordre ou corARMA(p = 2, q = 0) pour un processus autorégressif d’ordre 2 :

modele_gls <- nlme::gls(
  Stress ~ Effort + Sommeil,
  data = donnees_psy,
  correlation = corARMA(p = 2, q = 0, form = ~ Jour),
  method = "REML"
)
summary(modele_gls)

L’inspection de la sortie de cette modélisation gls révèle l’estimation conjointe des paramètres autorégressifs résiduels (notés Phi1 et Phi2 dans la console R), dont les grandeurs numériques convergeront précisément vers les paramètres générateurs (avoisinant 0.45 et 0.25 dans notre simulation). Les erreurs-types des coefficients Effort et Sommeil sont directement calibrées en tenant compte de cette structure de dépendance continue. De surcroît, cette approche permet d’utiliser des tests de rapport de vraisemblance formels pour comparer un modèle contraint sans autocorrélation avec le modèle enrichi de la structure corARMA, scellant de façon définitive la remédiation statistique au sein d’un cadre mathématique unifié et élégant.

11. Applications pratiques et études de cas en sciences du comportement

11.1 Cas 1 : Analyse de l’évaluation écologique momentanée (EMA)

Pour illustrer la portée concrète du test de Breusch-Godfrey dans la pratique clinique moderne, examinons le protocole d’une étude en psychiatrie translationnelle évaluant les déterminants écologiques momentanés de l’anxiété flottante chez des patients souffrant de trouble anxieux généralisé (TAG). Dans cette cohorte, 80 patients ont été équipés d’un dispositif portable transmettant six notifications journalières pendant 14 jours consécutifs, cumulant des dizaines de mesures séquentielles d’anxiété perçue, d’exposition au bruit environnemental et de charge d’interruptions professionnelles.

Lors de l’ajustement initial d’une régression pooled MCO prédisant l’anxiété par le bruit et les interruptions, l’application du test de Breusch-Godfrey à l’ordre p = 1 et à l’ordre p = 2 met en évidence une dépendance sérielle massive :

bgtest(reg_ema, order = 2, type = « F »)
F(2, 1115) = 142.36, p < .0001

Ce rejet écrasant a démontré que l’anxiété résiduelle manifestait une demi-vie biologique excédant largement l’intervalle entre deux bips successifs. L’analyse sans correction attribuait faussement au bruit ambiant un effet prédictif majeur (t = 4.82, p < .001). L’application d’une matrice robuste de Newey-West a réévalué l’erreur-type de ce prédicteur de plus de 80 %, recalculant la statistique d’épreuve à t_robuste = 1.74 (p = .082). L’effet présumé du bruit est ainsi apparu comme un artéfact d’autocorrélation : les épisodes de bruit coïncidaient par hasard avec des plateaux d’anxiété prolongée auto-entretenue. Cette rectification diagnostique a préservé l’équipe de recherche d’une conclusion clinique erronée, réorientant le protocole vers la modélisation de l’inertie émotionnelle endogène.

11.2 Cas 2 : Évaluation longitudinale de l’efficacité d’une psychothérapie

Notre deuxième étude de cas concerne un essai clinique randomisé évaluant l’efficacité longitudinale d’une thérapie cognitivo-comportementale (TCC) de troisième vague sur la dépression majeure. Vingt-quatre séances hebdomadaires ont été dispensées à une cohorte de patients, avec une mesure systématique hebdomadaire de la sévérité dépressive via l’inventaire de dépression de Beck (BDI-II). Le modèle de régression visait à estimer la pente de décroissance des symptômes au fil des semaines d’intervention thérapeutique.

L’ajustement de la régression linéaire sur la série chronologique consolidée des scores a produit une estimation très favorable de l’effet thérapeutique hebdomadaire, avec une réduction moyenne affichée de 0.85 point de BDI par semaine (p < .00001). Cependant, la passation systématique du test de Breusch-Godfrey a révélé une forte persistance d’ordre p = 1 au sein des résidus de modélisation :

bgtest(reg_depression, order = 1, type = « Chisq »)
LM test = 18.94, df = 1, p-value = 1.35e-5

Cette signature stochastique traduisait ce que les cliniciens nomment « l’inertie dépressive » : l’état affectif d’un patient lors d’une semaine donnée demeure fortement corrélé à son niveau de la semaine antérieure, indépendamment du gain thérapeutique linéaire théorique. En substituant aux MCO une transformation de Prais-Winsten sous R via la bibliothèque orcutt, le modèle a démontré que la vitesse réelle de rémission était en réalité plus lente que ne le suggérait l’ajustement naïf (pente corrigée de -0.62 point par semaine), et que les intervalles de confiance encadrant le pronostic de guérison devaient être élargis de 45 % pour demeurer statistiquement valides, fournissant aux praticiens des repères pronostiques d’un réalisme clinique incontestable.

11.3 Cas 3 : Modélisation de temps de réaction en psychologie cognitive

La troisième application pratique s’ancre dans le laboratoire de neuropsychologie expérimentale, dans le cadre d’une tâche d’attention soutenue et de vigilance psychomotrice (Psychomotor Vigilance Task ou PVT). Un sujet effectue 500 essais consécutifs où il doit appuyer sur un bouton dès l’apparition d’un stimulus visuel imprévisible. Le temps de réaction (TR, mesuré en millisecondes) est modélisé en fonction du temps écoulé depuis le début de la tâche (fatigue cognitive) et de la durée de l’intervalle pré-stimulus (effet d’attente préparatoire).

Puisque les essais se succèdent à intervalles rapprochés de quelques secondes, l’état attentionnel du sujet fluctue selon des ondes lentes de micro-sommeil et de sursauts de vigilance. Le diagnostic de régression linéaire standard a donné lieu au déploiement d’un test de Breusch-Godfrey à large spectre étendu jusqu’à l’ordre p = 5 :

bgtest(reg_tr, order = 5, type = « Chisq »)
LM test = 96.31, df = 5, p-value < 10^-15

L’hypothèse d’indépendance des essais consécutifs est volée en éclats. Les essais cognitifs successifs ne sauraient en aucun cas être appréhendés comme des réalisations indépendantes identiquement distribuées (i.i.d.). Face à cette structure de corrélation d’ordre élevé, les chercheurs ont abandonné le modèle MCO au profit d’un modèle gls() paramétré avec une matrice de corrélation résiduelle autorégressive de moyenne mobile corARMA(p = 3, q = 1). Cette spécification sophistiquée a permis de modéliser avec une fidélité mathématique remarquable la dissipation de l’effet d’inhibition post-erreur sur une fenêtre de trois essais consécutifs, tout en fournissant des inférences irréprochables quant à l’effet propre de la fatigue accumulée au cours de la tâche.

12. Bonnes pratiques, pièges méthodologiques et recommandations de publication

12.1 Erreurs courantes à éviter lors de l’application du test

La mise en œuvre du test de Breusch-Godfrey, bien que grandement facilitée par l’élégance syntaxique du langage R, est fréquemment entachée d’erreurs méthodologiques ou d’interprétations abusives au sein de la littérature de recherche. La première erreur courante réside dans le choix purement arbitraire ou dogmatique de l’ordre de retard maximal p. Se borner systématiquement à order = 1 par mimétisme avec le test obsolète de Durbin-Watson annihile la justification première du recours à Breusch-Godfrey. À l’inverse, choisir un ordre p excessif sans justification théorique ni validation par l’AIC/BIC dilue la puissance du test et s’apparente à une tentative déguisée de manipulation statistique.

Le second piège méthodologique majeur consiste à confondre la corrélation contemporaine entre variables avec l’autocorrélation sérielle des termes d’erreur. Des variables explicatives peuvent être fortement corrélées entre elles (multicolinéarité) sans que les résidus ne manifestent la moindre dépendance sérielle temporelle. Réciproquement, un modèle dont tous les régresseurs sont parfaitement orthogonaux contemporainement peut présenter une autocorrélation résiduelle dévastatrice si une dynamique temporelle ou une variable latente à évolution lente a été omise de la spécification. Le test de Breusch-Godfrey diagnostique exclusivement la structure du bruit résiduel ε et ne renseigne en rien sur l’existence d’une colinéarité au sein de la matrice X (qui relève du calcul des facteurs d’inflation de la variance, VIF).

Enfin, une bévue conceptuelle particulièrement déstabilisatrice consiste à appliquer aveuglément le test de Breusch-Godfrey à des séries temporelles non stationnaires présentant des racines unitaires ou des dérives stochastiques d’ordre I(1). Si la variable dépendante Y et les prédicteurs X sont intégrés sans être cointégrés, la régression linéaire produit des coefficients fallacieux et les résidus présentent une non-stationnarité qui fausse intégralement la distribution asymptotique du multiplicateur de Lagrange. Dans ce contexte, la statistique LM diverge vers l’infini indépendamment de toute considération de taille d’échantillon. Il est donc formellement requis de s’assurer de la stationnarité préalable des séries au moyen de tests de racines unitaires (tels que le test augmenté de Dickey-Fuller via tseries::adf.test()) avant d’engager le diagnostic d’autocorrélation des résidus.

12.2 Checklist diagnostique pour les chercheurs et praticiens R

Pour institutionnaliser une démarche diagnostique sans faille au sein de vos protocoles de modélisation quantitative sous R, il est vivement recommandé d’adopter une séquence d’assurance qualité standardisée. Cette check-list méthodologique récapitule les jalons indispensables à valider avant de soumettre vos conclusions à l’arbitrage académique :

• Vérification de l’ordonnancement temporel : Les observations sont-elles rigoureusement triées selon la séquence chronologique effective du protocole avant toute estimation ? La matrice contient-elle des ruptures temporelles ou des jours non observés qui imposeraient un traitement spécifique de panel non équilibré ?
• Diagnostic visuel exploratoire : Les résidus MCO ont-ils été extraits et visualisés chronologiquement ? Avez-vous inspecté conjointement les fonctions d’autocorrélation simple (ACF) et partielle (PACF) pour sonder visuellement la présence de pics franchissant les bandes de Bartlett ?
• Calibration raisonnée de l’ordre p : L’ordre maximal p retenu pour bgtest() est-il étayé par une analyse comparative des critères d’information (AIC/BIC) et par des considérations substantielles liées à la fréquence de passation des mesures ?
• Arbitrage sur la distribution de référence : Si votre échantillon comprend moins de 100 observations, avez-vous expressément spécifié l’argument type = "F" pour neutraliser la distorsion de taille de la statistique du Chi-deux ?
• Validation de l’absence de régresseurs I(1) : Les variables de votre modèle sont-elles stationnaires autour d’une tendance ou ont-elles été préalablement différenciées pour écarter le spectre d’une régression fallacieuse ?
• Triangulation des remédiations : En cas de rejet de H₀, avez-vous évalué l’impact de l’autocorrélation en comparant formellement les erreurs-types standards MCO avec les erreurs-types robustes de Newey-West issues de sandwich::NeweyWest() ? La formulation d’un modèle dynamique autorégressif ou l’utilisation de gls() a-t-elle été explorée ?

12.3 Directives de transparence et rédaction académique

Dans l’esprit des principes directeurs de l’Open Science Framework (OSF) et des recommandations de transparence méthodologique contemporaines (normes TOP, Transparency and Openness Promotion), la communication des diagnostics d’autocorrélation ne doit laisser aucune place à l’imprécision ou à l’omission stratégique. La rédaction scientifique contemporaine exige que l’intégrité de la chaîne de traitement statistique soit accessible et parfaitement auditable par les pairs.

Lors de la soumission de vos manuscrits, veillez à consigner dans le corps du texte ou au sein des annexes méthodologiques en ligne le code R intégral ayant généré vos diagnostics. La syntaxe exacte invoquée, incluant les arguments passés à la fonction bgtest(), doit être documentée de manière à ce qu’un relecteur indépendant puisse reproduire scrupuleusement les valeurs numériques obtenues à partir des données brutes partagées. En cas de non-respect avéré de l’hypothèse nulle d’indépendance résiduelle, mentionnez explicitement pourquoi la stratégie de remédiation adoptée (par exemple l’estimateur HAC de Newey-West ou le recours à un modèle linéaire à effets mixtes sous nlme) est théoriquement supérieure et comment elle altère les seuils de significativité des hypothèses initialement formulées.

Enfin, veillez à accompagner systématiquement la statistique du test de Breusch-Godfrey de graphiques de résidus post-correction. Démontrer visuellement et statistiquement qu’un modèle initialement contaminé par une mémoire sérielle AR(2) affiche, après inclusion d’un terme dynamique ou application d’un modèle gls, des résidus rigoureusement assimilables à un bruit blanc non corrélé constitue le gage de probité méthodologique le plus éclatant qu’un chercheur puisse apporter à la communauté scientifique.

Références

Box, G. E. P., Jenkins, G. M., Reinsel, G. C., & Ljung, G. M. (2015). Time series analysis: Forecasting and control (5th ed.). John Wiley & Sons. https://www.wiley.com/

Breusch, T. S. (1978). Testing for autocorrelation in dynamic linear models. Australian Economic Papers, 17(31), 334–355. https://doi.org/10.1111/j.1467-8454.1978.tb00635.x

Cochrane, D., & Orcutt, G. H. (1949). Application of least squares regression to relationships containing auto-correlated error terms. Journal of the American Statistical Association, 44(245), 32–61. https://doi.org/10.1080/01621459.1949.10483290

Durbin, J., & Watson, G. S. (1950). Testing for serial correlation in least squares regression: I. Biometrika, 37(3/4), 409–428. https://doi.org/10.2307/2332388

Durbin, J., & Watson, G. S. (1951). Testing for serial correlation in least squares regression: II. Biometrika, 38(1/2), 159–178. https://doi.org/10.2307/2332325

Fox, J., & Weisberg, S. (2019). An R companion to applied regression (3rd ed.). SAGE Publications. https://socialsciences.mcmaster.ca/jfox/Books/Companion/

Godfrey, L. G. (1978). Testing against general autoregressive and moving average error models when the regressors include lagged dependent variables. Econometrica, 46(6), 1293–1301. https://doi.org/10.2307/1913829

Greene, W. H. (2018). Econometric analysis (8th ed.). Pearson. https://www.pearson.com/

Harvey, A. C. (1990). The econometric analysis of time series (2nd ed.). MIT Press. https://mitpress.mit.edu/

Newey, W. K., & West, K. D. (1987). A simple, positive semi-definite, heteroskedasticity and autocorrelation consistent covariance matrix. Econometrica, 55(3), 703–708. https://doi.org/10.2307/1913610

Pinheiro, J. C., & Bates, D. M. (2000). Mixed-effects models in S and S-PLUS. Springer. https://doi.org/10.1007/b98882

Prais, S. J., & Winsten, C. B. (1954). Trend estimators and serial correlation (Cowles Commission Discussion Paper No. 383). University of Chicago. https://cowles.yale.edu/

Wooldridge, J. M. (2016). Introductory econometrics: A modern approach (6th ed.). Cengage Learning. https://www.cengage.com/

Zeileis, A., & Hothorn, T. (2002). Diagnostic checking in regression relationships. R News, 2(3), 7–10. https://cran.r-project.org/doc/Rnews/