Comment réaliser une ANOVA à mesures répétées sous R

L’analyse de variance à mesures répétées, communément désignée sous l’acronyme anglo-saxon d’ANOVA intra-sujets (Repeated Measures Analysis of Variance), constitue l’un des piliers méthodologiques et inférentiels les plus fondamentaux dans le domaine des sciences du comportement, des neurosciences cognitives et de la psychophysiologie expérimentale. Contrairement aux plans factoriels indépendants, où chaque unité expérimentale n’est exposée qu’à une seule modalité d’un facteur d’intérêt, les devis expérimentaux à mesures répétées évaluent les mêmes participants à travers une séquence temporelle définie ou sous une batterie complète de conditions expérimentales distinctes. Cette approche paradigmatique présente un avantage analytique exceptionnel : elle permet de s’affranchir de la variabilité interindividuelle statique, souvent considérable dans les populations humaines, afin de maximiser la sensibilité du test statistique face aux variations réelles induites par les manipulations expérimentales.

Dans la pratique contemporaine de la recherche universitaire et biomédicale, l’environnement logiciel statistique libre R s’est imposé comme la plateforme par excellence pour conduire ces investigations quantitatives complexes. Grâce à un écosystème modulaire particulièrement dynamique intégrant la philosophie du Tidyverse, ainsi que des bibliothèques hautement spécialisées telles que rstatix, ez ou emmeans, le chercheur dispose d’outils analytiques d’une flexibilité et d’une rigueur sans commune mesure avec les logiciels commerciaux traditionnels fermés. Cependant, la mise en œuvre de cette technique exige une maîtrise rigoureuse de ses postulats sous-jacents, en particulier le postulat de sphéricité de la matrice de covariance, dont la méconnaissance ou l’évaluation erronée peut gravement vicier les conclusions inférentielles en augmentant exponentiellement le taux d’erreur de première espèce.

Ce guide exhaustif a pour vocation d’offrir une référence méthodologique et computationnelle définitive pour les étudiants de cycle supérieur, les chercheurs en psychologie et les analystes de données. Nous y décortiquerons méticuleusement l’ensemble de la chaîne de traitement analytique : depuis la décomposition algébrique fondamentale des sommes de carrés et la formulation matricielle du modèle, jusqu’à la modélisation sous le langage R, le diagnostic visuel et inférentiel des postulats, l’ajustement automatique en cas de violations par les procédures d’epsilon de Greenhouse-Geisser et Huynh-Feldt, la décomposition fine des effets par contrastes post-hoc appariés, l’estimation des grandeurs d’effet généralisées et partielles, pour finir par la mise en forme graphique et textuelle rigoureuse conforme aux impératifs de publication de l’American Psychological Association (APA, 7e édition).

1. Fondements théoriques de l’ANOVA à mesures répétées en psychologie

1.1 Définition et principes des plans intra-sujets

La distinction conceptuelle entre un plan inter-sujets (ou plan à groupes indépendants) et un plan à mesures répétées (ou plan intra-sujets) repose sur l’assignation des unités d’observation aux conditions de l’expérience. Dans un plan inter-sujets classique, chaque individu sélectionné dans l’échantillon d’étude n’est assigné qu’à un et un seul niveau de la variable indépendante manipulée. Cette structure implique inévitablement que la variabilité observée entre les groupes résulte non seulement de l’impact de la variable indépendante d’intérêt, mais également de l’ensemble disparate des différences individuelles préexistantes qui séparent les participants composant les divers groupes, telles que leurs capacités cognitives basales, leurs prédispositions génétiques, leurs traits de personnalité stables ou leur historique développemental.

À l’opposé, les plans à mesures répétées soumettent l’ensemble des participants à l’intégralité des modalités de la variable indépendante, ou mesurent une variable dépendante donnée de manière sérielle à des temps successifs. Sur le plan structural, chaque participant devient ainsi son propre groupe témoin (own control). L’avantage méthodologique et statistique primordial réside dans l’élimination systématique de la variance interindividuelle stable du terme d’erreur statistique. Dans la mesure où les caractéristiques intrinsèques d’un individu demeurent globalement invariantes tout au long de la session de test, la variance associée à ces traits pérennes peut être statistiquement isolée et retirée de la variance résiduelle non expliquée. Il en découle une réduction substantielle du dénominateur du ratio de Fisher, ce qui accroît la sensibilité globale de l’analyse face aux micro-effets expérimentaux.

Ces dispositifs intra-sujets trouvent des applications prépondérantes au sein de la psychologie cognitive, de la neuropsychologie clinique et des neurosciences comportementales. À titre d’exemple, lors de l’étude des processus d’attention sélective au moyen de protocoles de chronométrie mentale (tels que la tâche de Stroop ou la tâche de Posner), il est capital d’évaluer le même individu dans les conditions congruentes, neutres et incongruentes pour déceler des différences de latence de l’ordre de quelques dizaines de millisecondes. En clinique, l’analyse de cinétiques de récupération post-accident vasculaire cérébral ou la mesure de l’efficacité d’un protocole de remédiation cognitive à travers des séances de mesure avant, pendant, immédiatement après et lors d’un suivi à six mois requièrent impérativement cette modélisation dépendante. Enfin, en neuropsychologie pharmacologique, où la variabilité métabolique individuelle est colossale, administrer successivement des agents psychotropes variés à un même panel de sujets restreint demeure l’approche la plus robuste scientifiquement.

L’autre atout décisif de cette architecture méthodologique est d’ordre pragmatique et concerne l’économie d’échantillonnage et la puissance statistique. Dans les domaines cliniques où le recrutement de patients porteurs de lésions cérébrales spécifiques ou de pathologies psychiatriques rares s’avère particulièrement ardu et chronophage, requérir quarante participants pour un plan inter-sujets à quatre modalités constitue une contrainte fréquemment rédhibitoire. Le plan intra-sujets permet d’obtenir une puissance statistique égale, voire supérieure, avec seulement dix participants évalués dans les quatre conditions, en maximisant le rendement informationnel par sujet tout en respectant les impératifs de rationalisation éthique des ressources humaines de recherche.

1.2 Décomposition mathématique de la variance intra-sujet

Pour appréhender l’essence computationnelle de l’ANOVA à mesures répétées, il convient d’examiner le partitionnement formel de la somme des carrés totale ($SS_{totale}$). Dans un modèle standard à mesures indépendantes à un facteur, la variation globale de la variable dépendante est simplement scindée en deux composantes orthogonales : la somme des carrés inter-groupes ($SS_{traitement}$) qui capte l’effet des modalités expérimentales, et la somme des carrés intra-groupes ($SS_{erreur}$), qui agrège indifféremment les idiosyncrasies individuelles et les erreurs de mesure fortuites. Le modèle intra-sujet raffine cette décomposition en fracturant le terme d’erreur généralisé en composantes distinctes et identifiables.

Dans un plan à un facteur intra-sujet, la somme des carrés totale s’exprime comme suit :
$$SS_{totale} = \sum_{i=1}^{n}\sum_{j=1}^{k}(Y_{ij} – \bar{Y}_{..})^2$$
Où $Y_{ij}$ désigne l’observation du participant $i$ sous la modalité $j$, $\bar{Y}_{..}$ représente la moyenne générale calculée sur l’intégralité des observations, $n$ est le nombre total de sujets uniques, et $k$ correspond au nombre total de conditions répétées. Cette quantité scalaire globale est d’abord divisée en deux sous-ensembles majeurs : la somme des carrés inter-sujets ($SS_{entre-sujets}$), qui quantifie l’écart des moyennes individuelles de chaque sujet par rapport à la moyenne globale, et la somme des carrés intra-sujets ($SS_{intra-sujets}$), qui mesure la dispersion des observations d’un sujet autour de sa propre moyenne individuelle :
$$SS_{entre-sujets} = k \sum_{i=1}^{n}(\bar{Y}_{i.} – \bar{Y}_{..})^2$$
$$SS_{intra-sujets} = \sum_{i=1}^{n}\sum_{j=1}^{k}(Y_{ij} – \bar{Y}_{i.})^2$$

La subtilité fondamentale du plan réside dans la partition ultérieure de cette somme des carrés intra-sujets. Cette dernière est segmentée en deux entités mutuellement exclusives : la somme des carrés due au traitement expérimental ($SS_{traitement}$ ou $SS_{condition}$), reflétant la variabilité imputable aux différentes modalités de l’intervention, et la somme des carrés de l’erreur résiduelle ($SS_{erreur}$ ou $SS_{résiduelle}$), qui correspond formellement à l’interaction entre les participants et les traitements ($Sujet \times Traitement$) :
$$SS_{traitement} = n \sum_{j=1}^{k}(\bar{Y}_{.j} – \bar{Y}_{..})^2$$
$$SS_{erreur} = \sum_{i=1}^{n}\sum_{j=1}^{k}(Y_{ij} – \bar{Y}_{i.} – \bar{Y}_{.j} + \bar{Y}_{..})^2$$
Cette interaction participant-traitement représente l’inconsistance de la réponse des sujets face aux conditions : si tous les sujets réagissent de façon rigoureusement proportionnelle et identique au changement de condition, cette quantité résiduelle tend vers zéro.

Dès lors, le calcul du ratio $F$ de Fisher s’effectue en rapportant le carré moyen du traitement expérimental ($MS_{traitement}$) non pas à l’ensemble de la variabilité résiduelle, mais strictement au carré moyen de l’erreur résiduelle ($MS_{erreur}$), où chaque carré moyen est obtenu en divisant la somme des carrés par ses degrés de liberté ($df$) respectifs :
$$df_{traitement} = k – 1$$
$$df_{erreur} = (n – 1)(k – 1)$$
$$F = \frac{MS_{traitement}}{MS_{erreur}} = \frac{SS_{traitement} / (k – 1)}{SS_{erreur} / ((n – 1)(k – 1))}$$
Cette construction mathématique met en lumière le mécanisme d’amplification de puissance de l’ANOVA intra-sujets : en soustrayant $SS_{entre-sujets}$ du dénominateur, la statistique de test $F$ subit une élévation substantielle, à condition que les scores des participants soient positivement corrélés à travers les conditions répétées, ce qui est la règle absolue en psychologie empirique.

1.3 Présentation du paradigme expérimental d’illustration

Afin d’ancrer de manière pragmatique les démonstrations théoriques, computationnelles et d’analyse de code de cet article, nous formalisons un paradigme expérimental standard issu du champ de la neuropsychopharmacologie cognitive. L’investigation porte sur l’évaluation de l’impact différentiel de quatre substances psychotropes distinctes sur la vitesse de traitement de l’information et la vigilance psychomotrice, quantifiées au moyen d’un chronomètre numérique de haute précision mesurant le temps de réaction (exprimé en millisecondes) lors d’une épreuve d’inhibition motrice de type Go/No-Go informatisée.

L’échantillon d’étude comprend une cohorte homogène de cinq patients diagnostiqués avec des troubles légers de l’attention (identifiés de façon anonyme par les codes P1, P2, P3, P4 et P5). Chaque participant a été soumis à l’ensemble des quatre conditions pharmacologiques suivantes, administrées selon une méthodologie en double aveugle croisé avec une période de sevrage et de clairance biologique (washout) adéquate entre chaque session pour prévenir tout effet de report ou d’accumulation métabolique rémanente :

• Placebo : Une solution saline neutre sans principe actif faisant office de ligne de base intra-individuelle.
• M1 (Molécule 1) : Un anxiolytique léger à action sédative présumée susceptible de ralentir les vitesses d’amorce motrice.
• M2 (Molécule 2) : Un psychostimulant dopaminergique théoriquement capable d’optimiser l’efficience psychomotrice.
• M3 (Molécule 3) : Un régulateur glutamatergique expérimental dont les répercussions comportementales demeurent indéterminées.

La structure formelle de ce dispositif implique une variable indépendante catégorielle intra-sujet comportant quatre niveaux mutuellement exclusifs (nommée drug), et une variable dépendante continue quantitative (nommée rt, pour reaction time). Le tableau expérimental intégral est ainsi un plan factoriel équilibré à mesures répétées à un facteur (désigné dans la nomenclature expérimentale par la notation $S_5 \times T_4$), totalisant vingt observations expérimentales couplées. Cette modélisation permettra de mettre en exergue de manière transparente la manipulation fine des données, les pièges syntaxiques sous R et l’arbitrage méthodologique entre les diverses fonctions statistiques disponibles.

2. Postulats statistiques préalables et gestion des violations

2.1 Le postulat de normalité des résidus de distribution

À l’instar de tous les modèles linéaires généraux paramétriques, l’ANOVA à mesures répétées repose sur l’hypothèse distributionnelle que les erreurs du modèle (les résidus empiriques) sont réparties selon une distribution normale gaussienne à chaque niveau de la variable indépendante, ou globalement au sein de la matrice résiduelle. Formellement, il est postulé que pour chaque condition $j$, l’erreur résiduelle $\epsilon_{ij}$ suit une loi normale d’espérance nulle et de variance conditionnelle : $\epsilon_{ij} \sim \mathcal{N}(0, \sigma^2_j)$.

Dans les situations empiriques de recherche, une distribution asymétrique (présentant une forte obliquité ou skewness) ou affectée par une anomalie d’aplatissement (kurtosis léptokurtique ou platykurtique) peut altérer la validité de l’inférence. Les conséquences directes d’une violation substantielle de la normalité se traduisent par une distorsion du taux d’erreur de première espèce ($\alpha$) réel, qui peut s’écarter significativement du seuil nominal de 0,05, ainsi que par une dégradation dramatique de la puissance statistique du test de Fisher. Les estimateurs des moyennes deviennent alors vulnérables aux valeurs extrêmes, biaisant le ratio des carrés moyens vers des conclusions fallacieuses.

Toutefois, la littérature méthodologique contemporaine démontre que le test d’ANOVA demeure remarquablement robuste face à des déviations modérées de la normalité, en vertu du théorème central limite, pour autant que la taille d’échantillon soit substantielle ($n > 30$) et que les distributions ne présentent pas de valeurs aberrantes outrancières. En revanche, dans le cadre d’échantillons cliniques de taille modeste, comme dans notre paradigme illustratif comportant cinq individus, cette robustesse asymptotique ne peut être invoquée d’office. L’évaluation de la normalité requiert alors une analyse combinée de tests d’adéquation statistiques, tels que le test de Shapiro-Wilk, et d’inspections graphiques diagnostiques via des diagrammes quantile-quantile (Q-Q plots).

2.2 Le postulat crucial de sphéricité de la matrice de covariance

Si la normalité constitue une exigence commune à l’ensemble des techniques paramétriques, l’hypothèse de sphéricité (souvent désignée sous le terme de condition de circularité de Huynh-Feldt) représente le postulat idiosyncrasique et le plus restrictif de l’analyse de variance à mesures répétées. Dans un devis inter-sujets, le statisticien postule l’homogénéité des variances entre les groupes indépendants (postulat d’homoscédasticité testé via l’épreuve de Levene). Dans un modèle intra-sujets, la structure d’interdépendance des mesures complexifie radicalement cette hypothèse.

Le postulat de sphéricité stipule formellement que les variances de toutes les différences possibles calculées entre chaque paire de conditions expérimentales doivent être rigoureusement équivalentes dans la population. Mathématiquement, pour un facteur comportant $k$ niveaux, si l’on calcule pour chaque sujet la différence entre son score sous la condition $A$ et son score sous la condition $B$ ($D_{A-B} = Y_{iA} – Y_{iB}$), puis la différence entre les conditions $A$ et $C$ ($D_{A-C} = Y_{iA} – Y_{iC}$), et ainsi de suite pour toutes les paires possibles, la variance de ces variables différentielles doit être constante :
$$Var(Y_A – Y_B) = Var(Y_A – Y_C) = Var(Y_B – Y_C) = dots = Cste$$

Une condition suffisante (quoique non strictement nécessaire) pour que la sphéricité soit satisfaite est l’hypothèse de symétrie composée, qui postule à la fois que toutes les variances des conditions soient identiques (homogénéité des variances) et que toutes les covariances entre les paires de conditions soient rigoureusement égales. Cette configuration est quasiment systématiquement violée dans les plans à mesures répétées longitudinaux ou temporels en psychologie : deux mesures prises à court intervalle temporel sont invariablement plus corrélées entre elles que deux mesures séparées par de longues périodes. Lorsque la condition de sphéricité n’est pas satisfaite, les distributions théoriques du ratio $F$ tabulées ne s’appliquent plus. Le ratio de Fisher calculé empiriquement suit une distribution théorique décalée, conduisant à une inflation gravissime de l’erreur de type I (rejet à tort de l’hypothèse nulle d’absence d’effet, le taux réel pouvant atteindre 0,10 voire 0,15 pour un seuil nominal fixé à 0,05).

2.3 Procédures de correction en cas d’asphéricité

Lorsque la matrice de variance-covariance dévie de la sphéricité parfaite, il n’est heureusement pas nécessaire d’abandonner le cadre de l’ANOVA au profit de modèles non paramétriques moins puissants. Les statisticiens ont élaboré des facteurs correctifs univariés élégants visant à ajuster à la baisse les degrés de liberté du numérateur et du dénominateur du test $F$, compensant ainsi précisément l’inflation de l’erreur $\alpha$. Ces facteurs sont synthétisés par l’indice $epsilon$ (epsilon) d’asphéricité, dont la valeur est bornée supérieurement par 1,0 (représentant une sphéricité absolue et parfaite) et inférieurement par $1/(k – 1)$ (correspondant au niveau d’asphéricité maximal envisageable théoriquement pour $k$ conditions).

La correction canonique la plus célèbre est celle de Greenhouse-Geisser ($\hat{\epsilon}_{GG}$). Développée dès 1959, cette méthode estime le degré d’asphéricité directement à partir des termes de la matrice de covariance empirique. Les degrés de liberté ajustés deviennent :
$$df_{traitement}^{*} = \hat{\epsilon}_{GG} \times (k – 1)$$
$$df_{erreur}^{*} = \hat{\epsilon}_{GG} \times (n – 1)(k – 1)$$
La statistique $F$ observée demeure strictement identique, mais la p-valeur associée est recalculée sous la loi de Fisher paramétrée par ces nouveaux degrés de liberté réduits (qui deviennent généralement fractionnaires). La correction de Greenhouse-Geisser est reconnue pour son caractère particulièrement conservateur, tendant à sous-estimer légèrement $epsilon$ lorsque celui-ci est proche de l’unité, ce qui majore le risque d’erreur de type II (non-détection d’un effet pourtant réel).

Pour parer à ce conservatisme excessif, Huynh-Feldt ont développé en 1976 une alternative moins pénalisante ($\tilde{\epsilon}_{HF}$), qui recalcule l’indice epsilon en corrigeant le biais d’échantillonnage de Greenhouse-Geisser, tout particulièrement lorsque les échantillons sont limités en taille. La formulation de Huynh-Feldt permet théoriquement à l’indice d’excéder la borne maximale de 1,0 ; lorsque cela se produit, il est coutumier de le tronquer manuellement à 1,00. Selon les recommandations actuelles en psychométrie (notamment formalisées par Maxwell et Delaney), il est préconisé d’utiliser la correction de Greenhouse-Geisser lorsque $hat{epsilon} < 0,75$, car elle protège infailliblement contre l’erreur de type I, et de basculer sur la correction de Huynh-Feldt lorsque$hat{epsilon} ge 0,75$, afin de préserver une puissance d’inférence optimale.

3. Configuration de l’environnement de travail et des packages sous R

3.1 Sélection des bibliothèques R spécialisées

L’exécution fluide, moderne et reproductible d’une ANOVA à mesures répétées au sein de la plateforme R nécessite la mise en place d’un environnement computationnel robuste. Bien que le moteur de base de R (base R) fournisse nativement des fonctions pour conduire de telles analyses (notamment par le biais de la fonction aov()), son ergonomie syntaxique s’avère complexe et ses sorties brutes manquent cruellement d’automatisations indispensables telles que le test de Mauchly ou l’application instantanée des corrections d’asphéricité. C’est pourquoi nous mobilisons un ensemble d’extensions contemporaines incontournables.

En premier lieu, le métapackage tidyverse constitue l’infrastructure indispensable pour l’ingestion, le nettoyage, la manipulation et la restructuration vectorielle des données grâce aux bibliothèques dplyr et tidyr, tout en offrant avec ggplot2 le standard absolu de la datavisualisation scientifique. En second lieu, la bibliothèque rstatix a été spécialement conçue pour convertir les routines statistiques conventionnelles en fonctions compatibles avec les principes du pipe (%>% ou le natif |>), permettant de générer des tableaux de résultats parfaitement agencés sous forme de data frames typés (tibbles).

Pour des analyses factorielles expérimentales avancées, la bibliothèque ez développée par Michael Lawrence demeure un standard historiquement éprouvé, particulièrement réputé pour sa fonction ezANOVA() qui exécute simultanément l’analyse factorielle, le test de sphéricité de Mauchly, les calculs des epsilons correctifs et les tailles d’effet généralisées. Enfin, pour l’exploration post-hoc et la décomposition analytique des interactions ou des effets principaux, le package emmeans (Estimated Marginal Means) constitue la référence épistémologique la plus solide, permettant de modéliser les moyennes marginales estimées et d’effectuer des tests de contrastes linéaires avec ajustement strict des seuils d’erreur de multiplicité.

3.2 Organisation architecturale des jeux de données

La structuration matricielle des jeux de données constitue l’une des sources d’erreurs les plus récurrentes rencontrées par les analystes lors de l’exécution d’analyses intra-sujets sous R. Deux formats d’agencement structurel s’opposent fondamentalement : le format dit large (wide format) et le format dit long (long format ou tidy data). Dans le format large, chaque unité expérimentale (chaque participant) est assignée à une ligne unique, et les différentes mesures répétées sont alignées horizontalement dans des colonnes dédiées distinctes (par exemple une colonne rt_placebo, une colonne rt_m1, une colonne rt_m2 et une colonne rt_m3).

Ce format large, bien qu’intuitif lors de la phase de recueil de données sous tableur de type Excel ou lors de l’encodage sur le logiciel SPSS, s’avère totalement inexploitable pour la grande majorité des algorithmes de modélisation du langage R, y compris les fonctions de régression à effets mixtes (lme4) et les moteurs de visualisation moderne (ggplot2). Pour être traitées adéquatement, les données doivent impérativement être converties au format long. Dans cette architecture conforme aux principes de Hadley Wickham, chaque ligne représente une mesure univariée singulière : une colonne identifie formellement le sujet, une seconde colonne précise l’étiquette de la condition expérimentale sous laquelle la mesure a été acquise, et une troisième colonne héberge la valeur numérique de la variable dépendante mesurée.

La manipulation de ces structures morphologiques au sein de R s’opère de façon déterministe au moyen des fonctions duales pivot_longer() et pivot_wider() issues de l’extension tidyr. La fonction pivot_longer() compresse les multiples colonnes de conditions pour créer un vecteur de variables indicatrices et un vecteur de scores unifié. À l’inverse, si des diagnostics multivariés exigent une décomposition matricielle (comme l’inspection des matrices de covariance pour le calcul manuel de la sphéricité), pivot_wider() permet de reconstituer instantanément la représentation tabulaire horizontale. La maîtrise parfaite de ces fonctions structurelles conditionne directement le bon déroulement de tout le pipeline d’analyse quantitative sous R.

4. Saisie, structuration et inspection du jeu de données expérimental

4.1 Création programmatique du tableau de données sous R

Afin de garantir l’autonomie et la parfaite reproductibilité de l’analyse, nous allons instancier directement dans la console R le jeu de données issu de notre paradigme expérimental. Cette démarche méthodologique consiste à structurer un tableau de données comportant exactement vingt lignes, correspondant aux cinq patients consécutivement testés sous les quatre modalités pharmacologiques (Placebo, M1, M2, M3). L’implémentation repose sur l’utilisation stratégique des fonctions de répétition vectorielle rep() associées à la concaténation standard.

Dans cette architecture, la colonne assignée aux identifiants individuels doit associer cinq identifiants de patients (de P1 à P5) répétés à travers les quatre conditions pharmacologiques. Corrélativement, le vecteur assigné aux traitements doit spécifier chaque étiquette de psychotrope répétée pour l’ensemble des cinq individus. Enfin, le vecteur contenant les temps de latence psychomotrice (rt) recueille les latences précises observées au laboratoire. Cette initialisation permet de construire un objet natif de classe data.frame, garantissant que chaque observation univariée est strictement indexée sur le couple sujet-condition approprié.

L’exécution de cette structure matricielle permet de vérifier la parfaite intégrité des valeurs numériques saisies. Il est primordial d’observer que les vingt observations sont complètes : aucune donnée manquante (valeur NA) ne doit entacher le protocole, car l’ANOVA à mesures répétées standard repose impérativement sur un plan factoriel équilibré et complet. En présence de données manquantes, les fonctions traditionnelles comme aov() procèdent à l’éviction totale du sujet concerné (suppression listwise), amputant drastiquement la représentativité de l’échantillon d’étude.

4.2 Déclaration explicite des variables facteurs

L’une des défaillances techniques les plus insidieuses sous R découle de la gestion implicite des types de données. Par défaut, lorsqu’un vecteur comporte des étiquettes textuelles ou des identifiants numériques (tels que 1, 2, 3... pour numéroter les participants), R est susceptible de les stocker sous la classe générique de chaînes de caractères (character) ou d’entiers numériques continus (numeric ou integer). Si l’on applique directement un algorithme d’ANOVA sur un tableau de données sans typage préalable explicite, le modèle traitera l’identifiant du patient ou le traitement pharmacologique comme une covariable linéaire continue, faussant l’intégralité du calcul des sommes de carrés et produisant une régression multivariée aberrante à un degré de liberté au lieu du test factoriel approprié.

Pour immuniser l’analyse contre cette dérive computationnelle, il est impératif d’utiliser la fonction factor() pour coercer formellement les colonnes qualitatives en variables nominales catégorielles. La variable identifiant le sujet (patient) doit être explicitement déclarée comme un facteur sans ordre intrinsèque. De la même manière, la variable pharmacologique (drug) doit être convertie en facteur nominal, dont les niveaux peuvent être formellement ordonnés si l’on souhaite établir le Placebo comme condition de contrôle basale de référence.

La validation de cette étape structurante s’effectue au moyen de fonctions d’interrogation du schéma de données, telles que str() ou glimpse() de dplyr. La vérification visuelle doit confirmer sans équivoque que la variable patient affiche l’attribut Factor w/ 5 levels et que la variable drug affiche l’attribut Factor w/ 4 levels. Dès que cette précaution architecturale est rigoureusement scellée, l’environnement R est optimalement configuré pour engager la phase d’exploration statistique descriptive des données.

5. Exploration visuelle et caractérisation descriptive des mesures

5.1 Synthèse descriptive conditionnelle

Avant d’initier toute procédure de modélisation inférentielle, une démarche empirique saine exige de procéder à une caractérisation descriptive méticuleuse de la distribution centrale et de la dispersion des scores à l’intérieur de chaque modalité expérimentale. Cette synthèse quantitative s’articule conventionnellement autour du calcul des moyennes empiriques ($\bar{Y}_{.j}$) et des écarts-types échantillonnaux ($SD_j$) pour chacune des conditions pharmacologiques (Placebo, M1, M2, M3).

Cependant, dans le cadre spécifique des devis à mesures répétées, l’utilisation de l’écart-type conventionnel ou de l’erreur type classique ($SE = SD / \sqrt{n}$) pour illustrer les intervalles d’incertitude induit un biais d’interprétation méthodologique sévère. En effet, l’écart-type brut intègre conjointement la variabilité intra-sujet liée au traitement et la variance interindividuelle globale (les différences basales entre les sujets lents et les sujets rapides). Dès lors, afficher des barres d’erreur basées sur les variances inter-sujets non corrigées masque l’homogénéité réelle des trajectoires individuelles et fausse l’appréciation visuelle de la significativité statistique.

Pour surmonter cette aporie descriptive, il est fortement recommandé d’implémenter les intervalles de confiance intra-sujets corrigés, basés sur la méthodologie de normalisation des scores développée par Loftus et Masson (1994) et affinée par Morey (2008). Cette procédure consiste à soustraire de chaque score individuel la moyenne générale du sujet concerné, puis à réinjecter la moyenne globale de l’expérience, avant d’appliquer un coefficient multiplicateur correctif dépendant du nombre de conditions ($k / (k – 1)$). La synthèse descriptive résultante, produite de façon élégante via group_by() et summarise() sous dplyr, fournit ainsi une estimation non polluée de la précision des moyennes conditionnelles, propice à une intégration directe dans les manuscrits de publication.

5.2 Visualisation des trajectoires individuelles

Si la synthèse tabulaire des moyennes est indispensable, l’examen graphique de la variabilité interindividuelle constitue l’étape la plus révélatrice d’une analyse intra-sujets réussie. L’outil graphique par excellence pour appréhender ces dynamiques est le graphique dit « spaghetti » (spaghetti plot), généré avec une grande finesse plastique sous ggplot2. Contrairement à un simple diagramme à barres qui agrège et opacifie l’information, ce tracé représente simultanément chaque participant sous la forme d’une ligne polyséquentielle reliant ses scores à travers l’ensemble des conditions expérimentales.

La contemplation d’un graphique de trajectoires individuelles permet d’isoler visuellement deux composantes orthogonales de l’architecture des données. D’une part, l’écartement vertical global séparant les différentes lignes individuelles illustre avec netteté l’ampleur de la variance interindividuelle de base : on discerne immédiatement les sujets structurellement lents (dont l’ensemble de la trajectoire se situe en haut de l’axe des ordonnées) des sujets globalement véloces. Cette variabilité correspond précisément à la somme des carrés inter-sujets ($SS_{entre-sujets}$) qui sera purgée de l’erreur résiduelle lors du test d’ANOVA.

D’autre part, le parallélisme relatif des trajectoires fournit un indice intuitif direct sur la régularité de l’effet expérimental et sur le respect du postulat de sphéricité. Si toutes les trajectoires individuelles adoptent une pente rigoureusement descendante ou ascendante similaire d’une condition à l’autre, cela indique une consistance parfaite de la réaction des sujets face aux drogues, minimisant le terme d’interaction sujet-traitement ($SS_{erreur}$). En superposant sur ce faisceau de lignes ténues la courbe épaisse de la moyenne générale d’échantillon, le chercheur dispose d’une vue synoptique exhaustive qui conjugue transparence des données brutes et lisibilité des tendances collectives.

5.3 Détection exploratoire des valeurs aberrantes

La présence d’observations aberrantes univariées (outliers) constitue l’une des menaces les plus pernicieuses pesant sur la validité de l’ANOVA paramétrique, en particulier au sein de petits effectifs où une unique mesure déviante peut modifier artificiellement la moyenne d’une condition et gonfler démesurément la somme des carrés résiduelle. L’identification rigoureuse de ces points nécessite un recours combiné à l’inspection géométrique par diagrammes en boîtes à moustaches (boxplots) et à des critères algorithmiques formalisés.

L’algorithme canonique de John Tukey, implémenté sous R via la fonction identify_outliers() du package rstatix, s’avère particulièrement efficace pour cette tâche d’audit qualité. Cet algorithme calcule pour chaque modalité expérimentale l’intervalle interquartile ($IQR = Q_3 – Q_1$) de la distribution des données. Les seuils de démarcation sont formellement établis selon deux frontières critiques :

• Valeurs extrêmes modérées : Toute observation située au-delà de la limite $Q_1 – 1,5 \times IQR$ ou excédant la limite $Q_3 + 1,5 \times IQR$. Dans la sortie logicielle, ces points sont qualifiés d’aberrants simples (attribut is.outlier = TRUE).
• Valeurs aberrantes sévères : Toute observation franchissant les bornes drastiques de $Q_1 – 3 \times IQR$ ou de $Q_3 + 3 \times IQR$ (attribut is.extreme = TRUE). Ces anomalies extrêmes résultent fréquemment d’erreurs de mesure techniques, de défaillances de capteurs ou de ruptures attentionnelles majeures chez le sujet.

La découverte d’une valeur aberrante sévère impose un arbitrage méthodologique délicat : le chercheur doit investiguer l’origine empirique de la déviance, évaluer l’opportunité d’une modélisation robuste, envisager une transformation mathématique des données (telle que la transformation logarithmique ou inverse, fréquemment requise pour les temps de réaction chronométriques), ou justifier formellement l’exclusion de l’observation incriminée afin de préserver l’invariance des estimateurs statistiques.

6. Vérification empirique des postulats sous R

6.1 Évaluation de la normalité univariée et distributionnelle

L’évaluation de la normalité au sein d’un modèle à mesures répétées s’opère soit sur la distribution des scores conditionnels pour chaque niveau de la variable indépendante, soit directement sur la série unifiée des résidus du modèle d’analyse de variance. Dans le cadre de protocoles impliquant un nombre limité d’observations, l’analyse condition par condition demeure la méthode diagnostique la plus commune et didactique.

La boîte à outils rstatix facilite grandement cette vérification grâce à la combinaison des verbes de commande group_by(drug) et de la fonction shapiro_test(rt). Le test de Shapiro-Wilk constitue à ce jour le test inférentiel d’adéquation gaussienne le plus puissant pour les échantillons restreints ($n < 50$). Il confronte l’hypothèse nulle ($H_0$) d’une distribution normale des données dans la population à l’hypothèse alternative ($H_1$) d’une déviation non gaussienne. L’obtention d’une statistique$W$ associée à une valeur de probabilité supérieure au seuil canonique de 0,05 conduit à la non-rejection de l’hypothèse nulle, validant formellement le postulat de normalité.

Néanmoins, les statisticiens rappellent avec insistance qu’au sein de très petits échantillons (tels que notre cohorte de cinq patients), le test de Shapiro-Wilk souffre d’un manque intrinsèque de puissance statistique : il s’avère fréquemment incapable de rejeter $H_0$ même en présence de déviations substantielles. C’est pourquoi le diagnostic inférentiel doit impérativement être corroboré par une inspection graphique par diagrammes quantile-quantile (Q-Q plots) générés via ggqqplot(). Ces graphiques confrontent les quantiles empiriques observés aux quantiles théoriques attendus sous une loi normale parfaite. L’alignement étroit des points de mesure le long de la bissectrice linéaire confirme graphiquement la viabilité de l’hypothèse de normalité résiduelle.

6.2 Évaluation diagnostique du postulat de sphéricité

L’évaluation du postulat de sphéricité représente le jalon décisionnel le plus déterminant dans la conduite de l’ANOVA intra-sujets. L’outil inférentiel standardisé consacré pour tester cette propriété de la matrice de covariance est le test de sphéricité de Mauchly (Mauchly, 1940). Ce test confronte l’hypothèse nulle de sphéricité parfaite ($H_0 : \Sigma \text{ est sphérique}$) à l’hypothèse alternative d’une violation substantielle de l’égalité des variances des différences ($H_1 : \Sigma \text{ est asphérique}$).

L’algorithme de Mauchly procède en extrayant le déterminant et la trace de la matrice orthogonale des contrastes normalisés construite sur les données expérimentales. Il condense cette relation sous une statistique scalaire notée $W$, dont la grandeur est mathématiquement comprise entre 0 et 1 :
$$W = \frac{|\mathbf{S}^*|}{\left[\frac{1}{k-1}\text{tr}(\mathbf{S}^*)\right]^{k-1}}$$
Où $\mathbf{S}^*$ représente la matrice de variance-covariance des différences transformées. Lorsque la matrice empirique converge vers une sphéricité absolue, la valeur de $W$ approche rigoureusement 1,0. À mesure que l’hétérogénéité des covariances s’accroît, $W$ s’amenuise vers zéro. Une transformation asymptotique du $W$ de Mauchly suit une distribution du Chi-carré ($\chi^2$), permettant de dériver une p-valeur formelle.

L’interprétation de ce test suit une logique rigoureusement dichotomique :

• Si la probabilité associée à la statistique de Mauchly est strictement supérieure à 0,05 ($p > 0,05$), le postulat de sphéricité ne peut être rejeté : le chercheur est formellement autorisé à interpréter la statistique $F$ standard univariée non corrigée avec ses degrés de liberté entiers nominaux.
• Si la probabilité s’avère inférieure ou égale à 0,05 ($p le 0,05$), l’hypothèse de sphéricité est formellement rejetée. L’asphéricité avérée de la matrice impose d’emblée l’application de facteurs correctifs d’ajustement des degrés de liberté, sous peine de commettre des erreurs d’inférence inacceptables.

Il importe de souligner une limitation structurelle : lorsque le facteur intra-sujet ne comprend que deux modalités ($k = 2$), le test de Mauchly ne peut être calculé mathématiquement (la fonction renverra une valeur NA). Cette propriété est parfaitement normale : n’existant qu’une unique différence possible entre deux conditions, le concept même d’égalité des variances de multiples différences devient sans objet, et la sphéricité est alors automatiquement garantie par définition.

7. Exécution de l’ANOVA à mesures répétées avec la fonction de base aov()

7.1 Paramétrage du terme d’erreur intra-sujet

Historiquement, l’exécution de l’analyse de variance au sein du système R natif repose sur l’appel à la fonction universelle aov(). Toutefois, l’ajustement d’un modèle à mesures répétées requiert une formulation syntaxique très spécifique, souvent déconcertante pour les praticiens habitués aux modèles de régression simple. La syntaxe canonique s’articule impérativement sous la forme suivante : aov(rt ~ drug + Error(patient/drug), data = mon_dataframe).

Le fondement théorique sous-jacent à cette formulation réside dans la spécification du terme d’erreur composite via l’opérateur Error(). Dans un modèle inter-sujets standard formulé par rt ~ drug, R assume l’indépendance statistique absolue de chaque ligne d’observation et alloue la totalité de la variabilité non expliquée à une strate résiduelle unique globale. L’inclusion de la clause Error(patient/drug) instruit explicitement le compilateur mathématique de R d’opérer une décomposition hiérarchique de l’espace des résidus en différentes strates de variance emboîtées (multistrata design).

L’expression patient/drug correspond formellement à une relation d’emboîtement (nesting), qui est automatiquement expansée par le parseur de R sous la forme patient + patient:drug. Dès lors, l’objet résultant (qui n’est pas un modèle linéaire simple mais une structure composite de classe aovlist) sépare distinctement la strate de variabilité inter-sujets (associée au terme patient) de la strate de variabilité intra-sujets (associée à l’interaction patient:drug). C’est précisément au sein de cette seconde strate que réside l’effet du traitement pharmacologique, testé rigoureusement contre cette interaction sujet-traitement faisant office de terme d’erreur spécifique.

7.2 Lecture et décomposition du tableau d’analyse de variance

L’interrogation de l’objet produit par aov() au moyen de la fonction générique summary() génère une restitution tabulaire segmentée en plusieurs sections analytiques distinctes correspondant aux différentes strates de projection orthogonale des moindres carrés. Pour interpréter correctement ces résultats, il convient d’ignorer la première section intitulée Error: patient, qui ne fait que quantifier la somme des carrés interindividuelle brute à $n – 1 = 4$ degrés de liberté sans formuler de test $F$ (puisqu’aucun facteur inter-sujets n’est modélisé ici pour confronter cette variabilité).

L’attention analytique doit se focaliser intégralement sur la section subséquente, expressément désignée sous l’étiquette Error: patient:drug. C’est à cet emplacement précis que le modèle déploie le test de Fisher de la manipulation expérimentale. Le tableau décline méthodiquement :

• La somme des carrés du facteur drug ($SS_{traitement}$) associée à ses degrés de liberté nominaux ($k – 1 = 3$).
• La somme des carrés résiduelle (notée Residuals, correspondant à $SS_{erreur}$) associée à ses degrés de liberté conjoints ($(n – 1)(k – 1) = 4 \times 3 = 12$).
• Les carrés moyens respectifs ($MS$), obtenus par le ratio des sommes de carrés sur leurs degrés de liberté.
• La statistique de test observée $F$, dérivée du quotient $MS_{drug} / MS_{Residuals}$.
• La p-valeur asymptotique (notée Pr(>F)), quantifiant la probabilité d’observer un tel ratio de variance sous l’hypothèse nulle d’égalité stricte des moyennes pharmacologiques.

Bien que la fonction aov() offre une rigueur computationnelle irréprochable et demeure un standard historique, elle souffre d’un défaut critique rédhibitoire dans la pratique de recherche appliquée : elle n’intègre aucun mécanisme natif d’évaluation de la sphéricité de Mauchly et ne fournit aucune correction d’ajustement d’epsilon automatique. Si la matrice de covariance viole le postulat de sphéricité, l’utilisateur d’aov() est contraint d’extraire manuellement les résidus matriciels pour calculer laborieusement les facteurs correctifs de Greenhouse-Geisser. Cette lacune justifie pleinement l’adoption de bibliothèques modernes automatisées.

8. Implémentation moderne et automatisée avec le package rstatix

8.1 Syntaxe et exécution de la fonction anova_test()

Pour répondre aux impératifs d’automatisation, d’ergonomie et d’intégration dans les flux de données modernes, la bibliothèque rstatix propose la fonction anova_test(). Cette fonction constitue une avancée méthodologique remarquable, car elle encapsule les calculs mathématiques sous-jacents les plus sophistiqués au sein d’une interface syntaxique limpide, épurée et entièrement orientée vers le paradigme du tidy data.

La formulation d’une ANOVA à mesures répétées via anova_test() ne requiert que la spécification intuitive de trois arguments essentiels :

• dv (dependent variable) : Spécifie la variable quantitative dépendante mesurée (dans notre cas, dv = rt).
• wid (within-subject identification) : Indique formellement la variable identifiant de manière univoque les sujets ou unités d’échantillonnage répétées (ici, wid = patient).
• within : Définit la ou les variables indépendantes manipulées de manière répétée au sein des mêmes individus (ici, within = drug).

L’appel s’exécute avec une élégance structurelle remarquable : resultat_anova <- mon_dataframe %>% anova_test(dv = rt, wid = patient, within = drug).

L’exécution de cet appel génère instantanément un objet composite de classe anova_test. Lorsque l’on affiche cet objet dans la console, la fonction ne renvoie pas une simple chaîne de texte brute non structurée, mais un tableau de données typé (un tibble) regroupant exhaustivement l’ensemble des grandeurs d’intérêt : l’effet testé, les degrés de liberté du numérateur (DFn), les degrés de liberté du dénominateur (DFd), la valeur numérique exacte de la statistique $F$, la p-valeur non ajustée, une colonne d’astérisques d’évaluation de la significativité statistique, ainsi que la taille d’effet mesurée par l’êta-carré généralisé (ges).

8.2 Application dynamique des corrections d’asphéricité

La supériorité opérationnelle du package rstatix se manifeste avec le plus d’éclat dans sa prise en charge native et totalement automatisée du postulat de sphéricité et des corrections associées. Lors de l’invocation de anova_test(), la fonction procède en coulisses à l’estimation immédiate du test de sphéricité de Mauchly ainsi qu’au calcul rigoureux des coefficients d’asphéricité $\hat{\epsilon}_{GG}$ et $\tilde{\epsilon}_{HF}$.

L’analyste peut interroger directement ces diagnostics sous-jacents en accédant aux attributs de l’objet d’ANOVA ou en exploitant les fonctions dédiées :

• get_mauchly(resultat_anova) : Extrait instantanément la statistique $W$ de Mauchly, la p-valeur associée et un indicateur logique stipulant si la sphéricité est violée au seuil nominal usuel de 5 %.
• get_anova_table(resultat_anova, correction = "auto") : Génère le tableau final d’ANOVA en appliquant de manière dynamique et intelligente la correction adéquate.

Le comportement de l’argument correction mérite une attention méthodologique singulière. Si l’argument est paramétré sur "auto" (ce qui constitue le réglage par défaut hautement recommandé), l’algorithme évalue d’abord le test de Mauchly : si la sphéricité est respectée ($p > 0,05$), aucune correction n’est appliquée ; en revanche, si le test s’avère statistiquement significatif ($p le 0,05$), la fonction substitue automatiquement les degrés de liberté nominaux par les degrés de liberté corrigés par le facteur d’epsilon de Greenhouse-Geisser. L’analyste peut également forcer manuellement l’application inconditionnelle de la correction de Greenhouse-Geisser (correction = "GG") ou de Huynh-Feldt (correction = "HF"), garantissant ainsi une conformité totale avec les préconisations éditoriales les plus rigoureuses des revues scientifiques internationales.

9. Modélisation robuste via la fonction ezANOVA du package ez

9.1 Configuration technique de l’appel ezANOVA()

Parallèlement à rstatix, le package ez représente un autre standard informatique historique extrêmement respecté par la communauté des psychologues quantitatifs et des neuroscientifiques. Sa fonction vedette, ezANOVA(), a été expressément conçue pour reproduire fidèlement, tout en les enrichissant, les fonctionnalités analytiques autrefois exclusives aux logiciels propriétaires tels que SPSS ou SAS, en instaurant un protocole d’interrogation particulièrement rigoureux.

L’interface d’appel de ezANOVA() impose une paramétrisation explicite de ses arguments clés, requérant la transmission directe des noms de colonnes sous forme de symboles non évalués ou d’expressions littérales :

• data : Le tableau de données d’entrée, impérativement structuré au format long.
• dv : La variable dépendante continue observée (ex. : .(rt) ou dv = rt).
• wid : Le vecteur identificateur des participants (ex. : .(patient)).
• within : Le ou les facteurs de mesures répétées intra-sujets (ex. : .(drug)).
• type : La méthode de calcul de la somme des carrés (Type I, Type II ou Type III).

L’argument type revêt une importance méthodologique fondamentale. Dans un plan intra-sujets équilibré sans données manquantes, les sommes des carrés de Type I (séquentielle), Type II (hiérarchique) et Type III (orthogonale marginale) convergent vers des résultats strictement identiques. Cependant, dès lors que l’expérience présente un déséquilibre structurel, la spécification type = 3 (qui est la norme préconisée par l’APA) devient indispensable pour estimer l’effet principal d’un facteur après avoir contrôlé marginalement l’ensemble des autres termes factoriels et interactions. La fonction ezANOVA() valide systématiquement l’équilibre du plan d’expérience avant de lancer les calculs et avertit immédiatement l’analyste par un message d’alerte explicite si le dispositif présente une quelconque asymétrie.

9.2 Interprétation des composantes de la sortie ezANOVA

L’un des atouts distinctifs majeurs de la fonction ezANOVA() réside dans le fait qu’un unique appel produit une liste composite structurée retournant simultanément l’ensemble des strates du diagnostic inférentiel. La sortie console s’articule méthodiquement en trois sections hiérarchisées :

La première section, simplement intitulée $ANOVA, présente la table standard d’analyse de variance univariée. On y retrouve l’effet factoriel testé (Effect), les degrés de liberté du numérateur (DFn) et du dénominateur (DFd), la valeur empirique du ratio $F$, la p-valeur non ajustée (p), la taille d’effet standardisée quantifiée par l’êta-carré généralisé (ges), ainsi qu’une colonne indicatrice de significativité binaire (p<.05).

La deuxième section est dédiée au test diagnostique : $Mauchly's Test for Sphericity. Elle détaille le résultat formel de l'épreuve de circularité de Mauchly appliquée au facteur intra-sujet concerné, en exposant la statistique $W$, la valeur de probabilité $p$ et une évaluation binaire signalant si l'hypothèse de sphéricité doit être formellement rejetée.

Enfin, la troisième section, dénommée $Sphericity Corrections, offre les estimations numériques des coefficients correctifs d'asphéricité : l'epsilon de Greenhouse-Geisser (GGe) et l'epsilon de Huynh-Feldt (HFe), accompagnés de leurs p-valeurs recalculées respectives (p[GG] et p[HF]). Cette disposition tabulaire unifiée permet au chercheur d'embrasser d'un seul coup d'œil l'ensemble de la séquence inférentielle : si le test de Mauchly de la section deux est significatif, il suffit de glisser le regard vers la section trois pour consigner directement la p-valeur pénalisée par Greenhouse-Geisser ou Huynh-Feldt dans le rapport de recherche final, éliminant tout risque d'erreur humaine de re-calcul.

10. Analyses post-hoc et décomposition des contrastes multiples

10.1 Tests de comparaisons par paires pour mesures dépendantes

L'obtention d'un ratio de Fisher statistiquement significatif au terme de l'ANOVA intra-sujets globale valide formellement l'hypothèse générale selon laquelle les moyennes des différentes modalités expérimentales ne sont pas toutes équivalentes dans la population d'étude ($H_1 : \exists (j, j') \text{ tel que } \mu_j \neq \mu_{j'}$). Néanmoins, cette décision inférentielle globale demeure d'une nature strictement omnibus : elle s'avère fondamentalement incapable d'indiquer précisément quelles paires spécifiques de conditions diffèrent significativement les unes des autres. Les temps de réaction sous Molécule 2 sont-ils réellement plus véloces que sous Placebo ? La Molécule 1 induit-elle un ralentissement significatif par rapport à la Molécule 3 ?

Pour résoudre cette interrogation substantielle, le chercheur doit obligatoirement procéder à des comparaisons multiples par paires (post-hoc pairwise comparisons). Sous R, l'approche computationnelle la plus limpide et contemporaine repose sur l'usage de la fonction pairwise_t_test() issue de la bibliothèque rstatix. Il est impératif de souligner une exigence méthodologique absolue : dans le cadre d'un plan à mesures répétées, l'argument de couplage doit impérativement être déclaré actif : paired = TRUE.

L'omission de cet argument paired = TRUE constitue une faute méthodologique majeure : la fonction exécuterait alors des tests t de Student pour échantillons indépendants, en utilisant comme terme d'erreur la variance interindividuelle globale au lieu de la variance de la différence intra-individuelle ($Var(D)$). La puissance statistique s'en trouverait dramatiquement effondrée, masquant de fait les contrastes significatifs réels existant entre les agents pharmacologiques.

10.2 Contrôle de l'inflation de l'erreur globale de type I

L'exécution itérative de multiples tests d'hypothèses sur un même jeu de données engendre un péril statistique redoutable connu sous le concept d'inflation de l'erreur de type I par famille d'hypothèses (family-wise error rate, FWER). Si l'on réalise $c$ comparaisons indépendantes au seuil nominal unitaire de $\alpha = 0,05$, la probabilité globale de commettre au moins une erreur de première espèce (déclarer à tort une différence significative) s'élève à :
$$\alpha_{FWER} = 1 - (1 - \alpha)^c$$
Pour notre paradigme comportant $k = 4$ molécules, le nombre total de combinaisons bilatérales possibles est de :
$$c = \frac{k(k - 1)}{2} = \frac{4 \times 3}{2} = 6 \text{ comparaisons}$$
En l'absence de procédure corrective d'ajustement, le risque d'erreur de première espèce cumulé atteint ainsi :
$$\alpha_{FWER} = 1 - (1 - 0,05)^6 \approx 1 - 0,735 = 0,265$$
Soit plus de 26 % de chances de publier une découverte fallacieuse.

Pour neutraliser ce risque inacceptable, la fonction pairwise_t_test() intègre nativement l'argument p.adjust.method, permettant d'appliquer divers algorithmes correctifs :

• Bonferroni ("bonferroni") : La correction la plus classique et intransigeante. Elle multiplie chaque p-valeur observée par le nombre total de comparaisons ($p_{ajustée} = \min(1, p \times c)$). Bien qu'infaillible pour juguler l'erreur $\alpha$, elle est notoirement hyper-conservatrice et pénalise excessivement la puissance statistique lorsque $c$ est élevé.
• Holm ("holm") : Procédure séquentielle par paliers (step-down) introduite par Sture Holm (1979). Elle ordonne les p-valeurs brutes de la plus faible à la plus forte et ajuste le seuil critique de manière progressive. Cette méthode domine strictement celle de Bonferroni : elle garantit un contrôle tout aussi infaillible du FWER tout en préservant une puissance de détection nettement supérieure. C'est l'ajustement par défaut recommandé dans la recherche contemporaine.
• Benjamini-Hochberg ("fdr") : Contrôle le taux de fausses découvertes (False Discovery Rate). Moins restrictif que le contrôle du FWER, il est principalement réservé aux plans exploratoires comportant des dizaines ou centaines de comparaisons multiples.

10.3 Approche alternative basée sur les moyennes marginales estimées (emmeans)

Si la méthode des tests t appariés successifs est très intuitive, l'approche la plus formelle et mathématiquement unifiée dans la modélisation statistique moderne repose sur le cadre des moyennes marginales estimées via la bibliothèque emmeans. Développée par Russell Lenth, cette méthode extrait les prédictions linéaires marginales du modèle préalablement ajusté pour exécuter les contrastes souhaités.

Pour mettre en œuvre cette approche sur un modèle intra-sujet sous R, on commence par convertir le modèle aov en un objet de référence marginale en exécutant la commande : emm <- emmeans(modele_aov, ~ drug). À partir de cette structure centrale, l'analyste déploie l'ensemble des comparaisons par paires au moyen de la fonction pairs(emm, adjust = "tukey"). La méthode de Tukey HSD (Honestly Significant Difference) est ici adaptée aux structures de mesures répétées en intégrant le carré moyen de l'erreur intra-sujet approprié.

L'immense valeur ajoutée méthodologique d'emmeans réside dans sa capacité à fournir non seulement des p-valeurs ajustées irréprochables, mais également les intervalles de confiance ajustés pour chaque différence de moyenne estimée (au moyen de la commande confint(pairs(emm))). L'analyste ne se borne plus à affirmer qu'un traitement diffère d'un autre ; il peut désormais déclarer avec précision que la molécule stimulante accélère la réponse motrice de 45 millisecondes en moyenne, avec un intervalle de confiance à 95 % compris par exemple entre 28 et 62 millisecondes, ancrant ainsi l'interprétation dans l'exactitude de l'estimation métrologique plutôt que dans le simple culte de la p-valeur binaire.

11. Quantification de la taille d'effet et puissance statistique

11.1 Calcul de l'êta-carré partiel (eta2 partiel) et généralisé

L'obtention d'une significativité statistique ($p < 0,05$) atteste de l'existence vraisemblable d'un effet non nul dans la population, mais ne renseigne en rien sur l'ampleur clinique, pratique ou cognitive de cet effet. Il est désormais impératif, conformément aux prescriptions internationales de publication scientifique, d'adjoindre systématiquement à chaque test$F$ un indice standardisé de taille d'effet (effect size). Dans les plans à mesures répétées, l'indice prédominant historique est l'êta-carré partiel ($\eta_p^2$).

Mathématiquement, l'êta-carré partiel exprime la proportion de la variance attribuable au facteur expérimental après avoir formellement retiré de l'équation la variance attribuable aux différences interindividuelles stables ($SS_{entre-sujets}$) :
$$\eta_p^2 = \frac{SS_{traitement}}{SS_{traitement} + SS_{erreur}}$$
Dans cette formulation, $SS_{erreur}$ correspond strictement à l'interaction $Sujet \times Traitement$. La bibliothèque rstatix ou la fonction eta_squared() du package effectsize permet d'extraire directement cette valeur. Selon les grilles d'interprétation conventionnelles établies par Jacob Cohen (1988) :

• $\eta_p^2 \approx 0,01$ : Effet d'amplitude modeste (faible).
• $\eta_p^2 \approx 0,06$ : Effet d'amplitude intermédiaire (moyen).
• $\eta_p^2 ge 0,14$ : Effet de forte amplitude (majeur).

Toutefois, de nombreux méthodologistes actuels (notamment Olejnik et Algina, 2003) soulignent les sérieuses limites comparatives de l'êta-carré partiel. En effet, la magnitude numérique de $\eta_p^2$ dépend intimement du design de recherche employé : un même effet expérimental sous-jacent produira une valeur de $\eta_p^2$ artificiellement plus élevée dans un plan intra-sujets que dans un plan inter-sujets indépendant, rendant toute comparaison méta-analytique périlleuse. C'est pourquoi ces auteurs recommandent l'usage privilégié de l'êta-carré généralisé ($\eta_G^2$ ou ges), qui incorpore l'ensemble des sources de variations systématiques afin d'assurer une parfaite comparabilité inter-protocoles.

11.2 L'oméga-carré partiel comme alternative non biaisée

Bien que l'êta-carré partiel demeure la mesure la plus abondamment rapportée dans les revues empiriques, il souffre d'un biais d'estimation structurel intrinsèque : il s'agit d'un estimateur purement descriptif qui présente une tendance systématique à surestimer la part de variance expliquée dans la population parente, ce phénomène de surajustement (overfitting) devenant particulièrement aigu au sein d'échantillons de recherche restreints ($n < 30$).

Pour pallier ce biais mécanique, les statisticiens préconisent le calcul de l'oméga-carré partiel ($\omega_p^2$). Dérivé des estimateurs non biaisés de composantes de variance de Cornfield et Tukey, l'oméga-carré partiel pénalise la somme des carrés de l'effet en lui soustrayant le carré moyen de l'erreur résiduelle, pondéré par le nombre d'observations et les degrés de liberté :
$$\omega_p^2 = \frac{df_{traitement} \times (MS_{traitement} - MS_{erreur})}{SS_{traitement} + (N_{total} - df_{traitement}) \times MS_{erreur} + SS_{entre-sujets}}$$
Où $N_{total}$ représente le nombre total de mesures agrégées ($n \times k$). Si la quantité résultante au numérateur s'avère négative (ce qui survient lorsque le $F$ observé est inférieur à 1), $\omega_p^2$ est conventionnellement fixé à 0.

L'obtention de l'oméga-carré partiel sous R s'opère sans difficulté au moyen de la bibliothèque effectsize via l'instruction omega_squared(modele, partial = TRUE). D'un point de vue épistémologique et éthique dans l'administration de la preuve empirique, rapporter $\omega_p^2$ plutôt que $\eta_p^2$ garantit une estimation d'une honnêteté intellectuelle maximale, prévenant la prolifération de revendications d'effets disproportionnellement gonflés dans la littérature scientifique.

11.3 Évaluation de la puissance a posteriori et planification future

L'évaluation quantitative d'un modèle statistique ne saurait être complète sans une analyse rigoureuse de sa puissance statistique ($1 - \beta$), c'est-à-dire la probabilité de rejeter formellement l'hypothèse nulle lorsque celle-ci est effectivement fausse. Dans le cadre d'un plan à mesures répétées, le calcul de la puissance se distingue fondamentalement des formules standard pour groupes indépendants en raison de l'intervention d'un paramètre crucial : le coefficient de corrélation moyen entre les mesures répétées ($r$).

En effet, plus la corrélation intra-individuelle entre les conditions est élevée, plus le dénominateur d'erreur résiduelle s'amenuise, et plus la puissance statistique s'envole. Le calcul de la puissance rétrospective (bien que critiqué si conduit de façon circulaire sur le $F$ observé) ou, plus fondamentalement, la planification prospective d'une réplication future s'effectue sous R à l'aide des bibliothèques pwr ou WebPower via la fonction spécialisée wp.rmanova().

Pour dimensionner de manière prospective un futur échantillon d'étude répliquant notre protocole neuropsychologique, le chercheur spécifie les paramètres cibles suivants :

• Le niveau de puissance désiré (conventionnellement fixé à $1 - \beta = 0,80$, voire 0,90).
• Le seuil nominal de significativité ($\alpha = 0,05$).
• Le nombre de conditions expérimentales répétées ($k = 4$).
• La taille d'effet prévisionnelle standardisée de Cohen ($f = \sqrt{\eta_p^2 / (1 - \eta_p^2)}$).
• Le coefficient de corrélation empirique observé entre les mesures appariées (par exemple, $r = 0,70$).
• Le facteur d'asphéricité anticipé $epsilon$ (en intégrant la valeur de Greenhouse-Geisser observée lors de l'étude pilote).

L'intégration de ces paramètres au sein de l'équation computationnelle de wp.rmanova() permet de déterminer l'effectif optimal requis avec une précision mathématique totale, évitant à la fois le sous-dimensionnement expérimental (qui expose à l'erreur de type II) et le sur-dimensionnement inutile (contraire aux normes éthiques de rationalisation de la recherche biomédicale).

12. Restitution académique et publication des résultats selon les normes APA

12.1 Règles formelles de rédaction de la 7e édition de l'APA

La transcription des analyses de variance au sein d'articles scientifiques ou de thèses de doctorat doit se conformer de façon intransigeante aux directives éditoriales édictées par l'American Psychological Association (APA 7th Edition). La norme prescrit une formalisation typographique extrêmement standardisée pour la déclaration des tests inférentiels paramétriques.

Le format textuel canonique pour rapporter une statistique $F$ d'ANOVA se structure sous la séquence suivante :
$$F(df_{numérateur}, df_{dénominateur}) = \text{valeur_F}, p = \text{valeur_p}, \eta_p^2 = \text{valeur_taille_effet}$$
Les impératifs typographiques suivants doivent être impérativement respectés :

• Les symboles statistiques alphabétiques latins majeurs ($F$, $p$, $t$, $N$) ainsi que les symboles d'indices grecs translittérés ($\eta_p^2$, $\omega_p^2$) doivent être obligatoirement composés en italique. Les chiffres numériques et les parenthèses demeurent quant à eux en caractères droits romains.
• En cas de correction pour violation du postulat de sphéricité, les degrés de liberté ajustés par Greenhouse-Geisser ou Huynh-Feldt doivent être rapportés sous leur forme fractionnaire exacte (avec deux décimales, par exemple : $F(1,82, 7,29) = 14,52$). De plus, la mention de la valeur numérique de l'indice epsilon utilisé ($epsilon = 0,61$) doit obligatoirement précéder ou accompagner le rapport du test.
• Concernant les valeurs de probabilité ($p$), si la p-valeur affichée par la console de R est infinitésimale (ex. : p < 0.001 ou sous notation scientifique 2.4e-05), il est strictement proscrit d'écrire $p = 0,000$. La rédaction formelle impose de formuler : $p < 0,001$. Pour toute autre valeur supérieure, la décimale exacte doit être précisée à deux ou trois chiffres après la virgule, sans zéro initial précédant la virgule ou le p\oint selon les normes anglophones (ex. :$p = ,024$).
• Chaque effet statistique inférentiel doit obligatoirement être articulé avec les statistiques descriptives complètes : les moyennes empiriques ($M$) et leurs écarts-types ($SD$) associés pour chacune des conditions doivent figurer explicitement soit dans le corps narratif du texte, soit au sein d'un tableau synthétique d'accompagnement.

12.2 Exemple de paragraphe de résultats pour publication

Afin de procurer un gabarit rédactionnel directement exploitable pour la communication scientifique, nous proposons ci-après l'agencement narratif exhaustif consignant l'intégralité de la chaîne d'analyse de notre paradigme expérimental neuropsychologique, articulant harmonieusement test omnibus, diagnostic de sphéricité, ajustement d'asphéricité et décomposition post-hoc :

« Une analyse de variance à mesures répétées à un facteur intra-sujet (Placebo vs. Molécule 1 vs. Molécule 2 vs. Molécule 3) a été conduite pour évaluer l'impact des différentes substances psychotropes sur la latence du temps de réaction psychomoteur. L'évaluation préalable du postulat de normalité résiduelle via le test de Shapiro-Wilk n'a révélé aucune déviation critique de la distribution gaussienne au sein des quatre conditions ($p > ,05$). Le test de circularité de Mauchly a toutefois mis en évidence une violation statistiquement significative du postulat de sphéricité de la matrice de variance-covariance, $W = 0,084, \chi^2(5) = 12,38, p = ,029$. En conséquence, les degrés de liberté du terme de traitement et de l'erreur ont été ajustés selon la procédure de correction conservatrice de Greenhouse-Geisser ($\hat{\epsilon} = 0,542$).

L'effet principal de la condition pharmacologique sur le temps de réaction s'avère hautement significatif et d'une magnitude considérable, $F(1,63, 6,51) = 24,87, p = ,001, \eta_p^2 = ,862, \eta_G^2 = ,415$. Les latences moyennes de réponse motrice s'échelonnent respectivement de $M = 345,2\text{ ms } (SD = 28,4)$ pour le Placebo, $M = 392,6\text{ ms } (SD = 31,1)$ sous Molécule 1, $M = 298,4\text{ ms } (SD = 22,7)$ sous Molécule 2, et $M = 338,1\text{ ms } (SD = 26,9)$ sous Molécule 3.

La décomposition analytique a posteriori de cet effet global par comparaisons multiples appariées deux à deux (tests t pour échantillons dépendants avec contrôle du taux d'erreur de première espèce selon la méthode séquentielle descendante de Holm) révèle que l'administration de la Molécule 2 entraîne une accélération statistiquement significative du temps de réaction comparativement à la condition Placebo, $t(4) = -4,89, p_{ajustée} = ,016, d = -1,82$, ainsi qu'une vélocité supérieure face à la Molécule 1, $t(4) = -6,72, p_{ajustée} = ,007$. À l'inverse, la Molécule 1 induit un allongement significatif de la latence de traitement vis-à-vis du Placebo, $t(4) = 3,98, p_{ajustée} = ,031$. Enfin, aucune différence statistiquement discernable n'émerge entre la condition Placebo et la Molécule 3, $t(4) = -0,81, p_{ajustée} = ,464$. Ces résultats corroborent pleinement l'hypothèse d'une facilitation psychomotrice sélective induite par le stimulant dopaminergique (M2), contrastant avec les propriétés sédatives inhibitrices de la Molécule 1. »

12.3 Production d'une figure de synthèse aux normes éditoriales

L'aboutissement culminant d'une communication empirique de haut niveau réside dans la transmission visuelle synthétique des données au moyen d'une figure de qualité éditoriale prête pour l'impression (camera-ready figure). Dans l'environnement R, la bibliothèque ggplot2 offre la plasticité nécessaire pour satisfaire aux standards esthétiques les plus exigeants.

Une figure conforme aux normes APA se caractérise par une épuration graphique stricte : élimination du quadrillage gris arrière-plan par l'application du thème minimaliste theme_classic(), choix d'une typographie proportionnelle neutre sans empattement, contrastes chromatiques distincts compatibles avec la reproduction en niveaux de gris, et étiquetage explicite et autonome des axes dimensionnels incluant les unités métrologiques standardisées (ex. : « Temps de réaction moyen (ms) »).

Au-delà de la représentation des moyennes conditionnelles sous forme de points massifs reliés par une trajectoire globale, l'élément technique prépondérant réside dans l'intégration rigoureuse des barres d'erreur basées sur les intervalles de confiance intra-sujets de Morey décrits à la section 5.1. L'utilisation d'intervalles inter-sujets non corrigés sur un graphique d'ANOVA à mesures répétées constitue une incohérence sémiotique majeure qui déroute l'œil du relecteur averti. Enfin, la figure doit idéalement intégrer des barres de contrastes superposées horizontales munies d'annotations textuelles conventionnelles (astérisques d'évaluation de la significativité : $*p < ,05$,$**p < ,01$,$***p < ,001$, ou mention explicite$ns$ pour les écarts non significatifs). Une telle construction iconographique offre une lisibilité instantanée et confère à la communication scientifique un professionnalisme méthodologique absolu.

Références

• Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate: A practical and powerful approach to multiple testing. Journal of the Royal Statistical Society: Series B (Methodological), 57(1), 289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x
• Cohen, J. (1988). Statistical power analysis for the behavioral sciences (2nd ed.). Lawrence Erlbaum Associates.
• Greenhouse, S. W., & Geisser, S. (1959). On methods in the analysis of profile data. Psychometrika, 24(2), 95–112. https://doi.org/10.1007/BF02289823
• Holm, S. (1979). A simple sequentially rejective multiple test procedure. Scandinavian Journal of Statistics, 6(2), 65–70. https://www.jstor.org/stable/4615733
• Huynh, H., & Feldt, L. S. (1976). Estimation of the Box correction for degrees of freedom from sample data in randomized block and split-plot designs. Journal of Educational Statistics, 1(1), 69–82. https://doi.org/10.3102/10769986001001069
• Kassambara, A. (2020). rstatix: Pipe-friendly framework for basic statistical tests (R package version 0.7.0). https://cran.r-project.org/package=rstatix
• Lawrence, M. A. (2016). ez: Easy analysis and visualization of factorial experiments (R package version 4.4-0). https://cran.r-project.org/package=ez
• Lenth, R. V. (2023). emmeans: Estimated marginal means, aka least-squares means (R package version 1.8.5). https://cran.r-project.org/package=emmeans
• Loftus, G. R., & Masson, M. E. (1994). Using confidence intervals in within-subject designs. Psychonomic Bulletin & Review, 1(4), 476–490. https://doi.org/10.3758/BF03210951
• Mauchly, J. W. (1940). Significance test for sphericity of a normal $n$-variate distribution. The Annals of Mathematical Statistics, 11(2), 204–209. https://doi.org/10.1214/aoms/1177731909
• Maxwell, S. E., & Delaney, H. D. (2004). Designing experiments and analyzing data: A model comparison perspective (2nd ed.). Lawrence Erlbaum Associates.
• Morey, R. D. (2008). Confidence intervals from normalized data: A correction to Cousineau (2005). Tutorials in Quantitative Methods for Psychology, 4(2), 61–64. https://doi.org/10.20982/tqmp.04.2.p061
• Olejnik, S., & Algina, J. (2003). Generalized eta and omega squared statistics: Measures of effect size for some common research designs. Psychological Methods, 8(4), 434–447. https://doi.org/10.1037/1082-989X.8.4.434
• R Core Team. (2023). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.r-project.org/
• Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L. D., François, R., Grolemund, G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Pedersen, T. L., Miller, E., Bache, S. M., Müller, K., Ooms, J., Robinson, D., Seidel, D. P., Spinu, V., Takahashi, K., Vaughan, D., Wilke, C., Woo, K., & Yutani, H. (2019). Welcome to the Tidyverse. Journal of Open Source Software, 4(43), 1686. https://doi.org/10.21105/joss.01686