Comment réaliser une MANOVA dans R

L’investigation empirique des sciences comportementales, cognitives et biomédicales se heurte de manière récurrente à la nature intrinsèquement multidimensionnelle de ses objets d’étude. Qu’il s’agisse de mesurer l’efficacité d’une psychothérapie sur une constellation de symptômes cliniques ou d’évaluer l’impact d’une lésion cérébrale sur différentes facettes des fonctions exécutives, les chercheurs ne peuvent que rarement résumer un phénomène complexe à un indicateur scalaire unique. L’approche traditionnelle, qui consiste à isoler chaque variable dépendante dans une analyse de variance univariée distincte, se heurte immédiatement à deux écueils majeurs : la prolifération incontrôlée de l’erreur de type I et l’occultation systématique des interdépendances et corrélations unissant les dimensions observées. Face à ces restrictions méthodologiques, l’analyse de variance multivariée (MANOVA) s’impose comme une extension rigoureuse, élégante et puissante du modèle linéaire généralisé.

Au sein de l’environnement de calcul statistique R, la conduite d’une MANOVA ne se résume pas à l’exécution d’une fonction élémentaire. Elle requiert une démarche analytique globale qui articule l’évaluation minutieuse des postulats probabilistes, la décomposition géométrique des matrices de dispersion, l’arbitrage entre les critères de décision canoniques et le déploiement d’analyses discriminantes ou univariées complémentaires. Ce guide propose une immersion exhaustive dans les fondements mathématiques, les arbitrages conceptuels et l’application programmatique de la MANOVA sous R. À travers un fil conducteur méthodologique ancré dans la recherche psychologique contemporaine, nous explorerons les mécanismes sous-tendant la modélisation vectorielle pour doter les chercheurs et analystes de données d’un cadre opératoire conforme aux standards de publication scientifique les plus exigeants.

En abordant les étapes successives de la préparation des matrices de données jusqu’aux visualisations de l’espace hypothèse-erreur (HE plots) et aux alternatives robustes ou non paramétriques, cet article vise à combler le fossé entre la théorie matricielle pure et l’interprétation concrète des résultats empiriques. L’adoption des normes de la septième édition de l’American Psychological Association (APA) servira d’ancrage tout au long de notre parcours pour assurer une communicabilité et une reproductibilité optimales des inférences statistiques obtenues.

1. Fondements théoriques de la MANOVA et distinction avec l’ANOVA

1.1 De l’analyse univariée à l’approche multivariée

L’analyse de variance univariée (ANOVA) constitue depuis plusieurs décennies l’outil par excellence pour tester l’égalité des moyennes entre plusieurs populations réparties sous l’effet d’une ou plusieurs variables indépendantes qualitatives. Pourtant, les construits psychologiques, psychiatriques et psychosociaux se caractérisent par une architecture profondément multifactorielle. Considérons un protocole destiné à tester l’efficacité d’une prise en charge de la dépression majeure : le clinicien évalue simultanément la détresse affective, l’anhédonie, les altérations neurovégétatives et les pensées intrusives. Réduire une telle constellation sémiologique à une note globale arbitrairement sommative implique l’hypothèse héroïque – et souvent fausse – que chaque sous-dimension contribue de manière uniforme au tableau clinique. À l’inverse, analyser chaque dimension de manière totalement disjointe via des tests ANOVA isolés méconnaît la réalité psychométrique sous-jacente.

Sur le plan mathématique, l’approche multivariée remplace la variable dépendante scalaire unique $Y$ par un vecteur aléatoire de réponses multidimensionnel $\mathbf{y}_i = [y_{i1}, y_{i2}, dots, y_{ip}]^top$, où $p$ correspond au nombre de variables dépendantes mesurées conjointement pour chaque unité statistique $i$. Cette transition implique que la dispersion n’est plus quantifiée par une simple variance univariée $\sigma^2$, mais par une matrice de variance-covariance de dimension $p \times p$, notée $boldsymbol{\Sigma}$. Les éléments diagonaux de $boldsymbol{\Sigma}$ contiennent les variances des $p$ indicateurs, tandis que les éléments hors-diagonale décrivent les covariances bilatérales entre les dimensions.

L’un des avantages les plus évidents de la modélisation vectorielle réside dans le contrôle rigoureux de l’inflation de l’erreur globale de type I, également désignée sous le vocable de taux d’erreur par famille (familywise error rate). Lorsqu’un expérimentateur conduit $p$ tests univariés au seuil nominal $\alpha = 0{,}05$, la probabilité cumulative d’engendrer au moins un faux positif s’élève théoriquement à $1 – (1 – \alpha)^p$, sous l’hypothèse d’indépendance statistique des mesures. Pour un vecteur de quatre variables dépendantes, ce risque global atteint déjà $1 – (0{,}95)^4 \approx 0{,}185$, soit près de 19 %, ce qui contrevient au principe fondamental de précaution inférentielle. Bien que des procédures d’ajustement univariées a posteriori, telles que les corrections de Bonferroni ou de Holm-Bonferroni, permettent d’endiguer cette dérive, elles s’accompagnent d’un coût sévère en matière de puissance statistique en ignorant délibérément l’espace de corrélation intraclasse partagé par les réponses.

Au-delà du seul contrôle de l’erreur de type I, la MANOVA offre une sensibilité analytique supérieure en capturant la structure des corrélations sous-jacentes. Deux groupes expérimentaux peuvent ne manifester aucune divergence statistiquement perceptible sur le plan univarié le long de l’axe $Y_1$ ou le long de l’axe $Y_2$, tout en présentant une séparation spatiale spectaculaire lorsqu’ils sont projetés dans le plan bivarié formé par l’intersection des deux mesures. La covariance entre les indicateurs modifie la géométrie des nuages de points ; en observant les centroïdes multivariés plutôt que les moyennes marginales, la MANOVA détecte des patrons de différenciation que l’ANOVA univariée est incapable de discriminer.

1.2 Formulation mathématique générale du modèle linéaire multivarié

Le modèle linéaire multivarié généralise directement le formalisme matriciel de la régression et de l’analyse de variance classique. Soit un échantillon de $N$ individus répartis dans un plan expérimental comportant $g$ groupes ou conditions. Les observations sont rassemblées dans une matrice de réponse $\mathbf{Y}$ de dimension $N \times p$, tandis que l’architecture du plan expérimental est encodée dans la matrice d’expérience $\mathbf{X}$ de dimension $N \times q$, où $q$ représente le nombre de paramètres de position estimés (incluant l’intercept et les contrastes factoriels). Le modèle s’écrit de manière synthétique :

$$\mathbf{Y} = \mathbf{X}\mathbf{B} + \mathbf{E}$$

Dans cette équation, $\mathbf{B}$ représente la matrice des coefficients de régression ou des effets de taille $q \times p$, liant chaque terme du prédicteur aux $p$ variables dépendantes. La matrice $\mathbf{E}$, de dimension $N \times p$, rassemble les résidus multivariés. On formule l’hypothèse que chaque ligne du vecteur d’erreur suit une loi normale multivariée d’espérance nulle et de matrice de covariance résiduelle homogène : $boldsymbol{\varepsilon}_i \sim \mathcal{N}_p(\mathbf{0}, boldsymbol{\Sigma})$.

L’inférence statistique au sein du modèle multivarié repose sur la décomposition spectrale de la dispersion globale en deux matrices fondamentales de sommes des carrés et des produits croisés (SSCP) : la matrice de l’hypothèse, notée $\mathbf{H}$ (pour Hypothesis), et la matrice de l’erreur résiduelle, notée $\mathbf{E}$ (pour Error). Ces deux matrices de dimension $p \times p$ constituent la généralisation multidimensionnelle des sommes des carrés inter-groupes ($SS_{\text{inter}}$) et intra-groupes ($SS_{\text{intra}}$) de l’ANOVA univariée :

$$\mathbf{H} = \sum_{j=1}^g n_j (\bar{\mathbf{y}}_j – \bar{\mathbf{y}}_\cdot)(\bar{\mathbf{y}}_j – \bar{\mathbf{y}}_\cdot)^top$$

$$\mathbf{E} = \sum_{j=1}^g \sum_{i=1}^{n_j} (\mathbf{y}_{ij} – \bar{\mathbf{y}}_j)(\mathbf{y}_{ij} – \bar{\mathbf{y}}_j)^top$$

où $n_j$ représente l’effectif du groupe $j$, $\bar{\mathbf{y}}_j$ désigne le centroïde (vecteur moyen) du groupe $j$, et $\bar{\mathbf{y}}_\cdot$ figure le grand centroïde global de l’échantillon. Pour évaluer la divergence spatiale entre les groupes relativement au bruit intra-groupe, on étudie l’opérateur linéaire $\mathbf{E}^{-1}\mathbf{H}$. Les valeurs propres non nulles $\lambda_1 ge \lambda_2 ge dots ge \lambda_s$ de cette matrice $\mathbf{E}^{-1}\mathbf{H}$ (avec $s = min(p, g – 1)$) quantifient les ratios maximaux de variance inter-groupes sur la variance intra-groupes le long des axes canoniques discriminants.

Il importe de distinguer clairement cette approche transversale inter-sujets de l’analyse multivariée pour mesures répétées. Dans une MANOVA inter-sujets standard, chaque ligne de $\mathbf{Y}$ correspond à des construits conceptuellement hétérogènes mesurés au même moment sur des participants indépendants. Dans un plan à mesures répétées analysé via l’optique multivariée (le profil d’analyse multivarié), les différentes colonnes de $\mathbf{Y}$ représentent la répétition temporelle d’un même indicateur scalaire dans des conditions successives. Cette distinction structurelle conditionne l’interprétation des contrastes matriciels et la nature des postulats de sphéricité.

1.3 Comparaison critique entre MANOVAs successives et séries d’ANOVAs

La question du choix entre le recours direct à une MANOVA ou l’enchaînement de séries d’ANOVAs univariées corrigées a suscité de vifs débats au sein de la communauté des statisticiens appliqués. Une première limite de l’approche univariée séquentielle concerne la redondance informationnelle. Si plusieurs variables dépendantes partagent une forte inter-corrélation, exécuter des tests d’hypothèse parallèles conduit à tester plusieurs fois la même parcelle de variance latente, sans apporter de clarification sur la spécificité des effets observés. La MANOVA traite ce faisceau de covariance de manière unifiée en le condensant sur des dimensions orthogonales discriminantes.

D’autre part, la MANOVA offre une sensibilité géométrique inaccessible à l’ANOVA univariée. Il est possible qu’aucun des tests univariés ne franchisse le seuil de significativité statistique, alors même que le test multivarié global s’avère hautement significatif. Cette situation émerge dès lors que l’écart entre les centroïdes s’exprime dans une direction oblique non alignée avec les axes initiaux des coordonnées univariées. Dans cette configuration, la covariance résiduelle masque l’effet lorsque les variables sont prises isolément, mais amplifie le signal de séparation dès lors que le vecteur de réponses est considéré comme une entité spatiale intégrée.

Néanmoins, la MANOVA n’est pas dépourvue d’inconvénients. L’incorporation inconsidérée d’indicateurs dépendants dépourvus de corrélation théorique ou empirique avec l’intervention expérimentale dilue la puissance statistique. Chaque variable dépendante additionnelle retranche un degré de liberté au niveau du dénominateur d’erreur pour l’estimation de la matrice de covariance. Si une variable n’apporte aucune variance explicative par rapport aux groupes, elle élève le seuil critique d’inférence sans contribuer à la valeur propre de $\mathbf{H}$, réduisant la puissance de détection du test global.

En conséquence, les critères méthodologiques recommandent d’opter pour la MANOVA uniquement lorsque les variables dépendantes reflètent les dimensions corrélées d’un même macro-construit théorique sous-jacent. Si les indicateurs empiriques sont totalement orthogonaux ou s’ils répondent à des questions scientifiques conceptuellement indépendantes, une série d’ANOVAs univariées, assortie d’un protocole d’ajustement du taux d’erreur par famille rigoureux, s’avère plus parcimonieuse et statistiquement préférable.

2. Justification méthodologique et applications en psychologie

2.1 Construits multidimensionnels en psychologie cognitive et clinique

La recherche en psychologie clinique contemporaine a largement abandonné les approches nosographiques catégorielles monolithiques au profit de modèles dimensionnels et transdiagnostiques, tels que le cadre Research Domain Criteria (RDoC) ou la taxonomie HiTOP (Hierarchical Taxonomy of Psychopathology). Dans ce contexte conceptuel, l’évaluation d’une pathologie affective ou d’un protocole d’intervention psychothérapeutique requiert d’appréhender des profils sémiologiques multiformes. Par exemple, une étude portant sur l’impact de la thérapie d’acceptation et d’engagement (ACT) face à un traitement cognitivo-comportemental classique évalue typiquement un vecteur associant la dépression clinique (mesurée par le BDI-II), l’anxiété état-trait (STAI), l’évitement expérientiel (AAQ-2) et la flexibilité psychologique. L’interdépendance structurelle de ces scores rend la MANOVA méthodologiquement incontournable pour statuer sur l’efficacité globale du traitement avant d’en isoler les dynamiques univariées.

En neuropsychologie et en psychologie cognitive expérimentale, les protocoles exploitent quasi-systématiquement des batteries composites pour mesurer les systèmes exécutifs centraux. L’évaluation de l’inhibition motrice, de la mémoire de travail et de la flexibilité cognitive implique la collecte synchronisée de temps de réaction en millisecondes, de taux d’erreurs de commission, et d’indices de sensibilité perceptuelle ($d’$ de la théorie de détection du signal). Isoler ces mesures appauvrirait considérablement l’analyse, car les sujets opèrent des compromis vitesse-précision (speed-accuracy tradeoff). Deux groupes expérimentaux peuvent sembler similaires en temps de réaction moyen, mais diverger drastiquement en tenant compte simultanément de l’exactitude de leurs réponses. La MANOVA intègre ces arbitrages comportementaux au cœur de la statistique inférentielle.

L’évaluation des interventions psychosociales et éducatives tire un parti similaire de l’analyse multivariée. Lorsque l’on évalue un programme pédagogique innovant destiné à des adolescents scolarisés, les retombées ne se bornent pas aux performances académiques pures ; elles englobent la motivation intrinsèque, l’auto-efficacité perçue et les conduites prosociales. Le recours à une MANOVA permet d’affirmer avec rigueur statistique si le programme modifie l’empreinte comportementale et cognitive globale des participants, tout en prévenant les conclusions artificielles fondées sur des sous-échelles déconnectées.

2.2 Gestion de la multicolinéarité et structure des corrélations

Le bon fonctionnement d’une MANOVA dépend étroitement de la structure des inter-corrélations observées entre les variables dépendantes. Sur le plan psychométrique, il existe une « zone optimale » de corrélation linéaire. Idéalement, les coefficients de corrélation de Pearson $r$ entre les variables du vecteur de réponse doivent se situer dans une plage comprise approximativement entre $0{,}30$ et $0{,}70$. Dans cette configuration, les variables partagent une quantité substantielle de variance commune relative au construit théorique ciblé, tout en conservant une part d’information résiduelle propre capable d’enrichir la discrimination spatiale des centroïdes.

Lorsque les corrélations bivariées franchissent le seuil critique de $0{,}80$ ou $0{,}90$, le modèle bascule dans une zone de multicolinéarité sévère ou de quasi-singularité matricielle. D’un point de vue algébrique, si deux variables dépendantes sont redondantes, l’une constitue une combinaison quasi-linéaire de l’autre. Le déterminant de la matrice de dispersion intra-groupe $\mathbf{E}$ s’approche dangereusement de zéro ($|\mathbf{E}| \approx 0$). L’opération d’inversion matricielle $\mathbf{E}^{-1}$ nécessaire au calcul des critères multivariés devient numériquement instable, amplifiant les erreurs d’arrondi, faussant le calcul des valeurs propres canoniques et provoquant une explosion artificielle des erreurs types des estimations.

À l’inverse, si les variables dépendantes affichent des corrélations proches de zéro ($r \approx 0{,}00$), l’apport mathématique de la MANOVA s’évanouit. Lorsque les indicateurs sont rigoureusement orthogonaux, la matrice $\mathbf{E}$ devient une matrice diagonale. Dans cette éventualité, la statistique multivariée globale ne fait que reproduire la somme des statistiques $F$ univariées, tout en subissant une pénalisation en degrés de liberté due à l’évaluation conjointe. L’analyste se retrouve alors dans une situation défavorable où la MANOVA perd de sa puissance statistique sans apporter de gain conceptuel par rapport à une batterie d’analyses univariées indépendantes corrigées.

2.3 Considérations relatives à la taille d’échantillon et à la puissance

La planification d’une MANOVA exige une réflexion rigoureuse quant à la puissance statistique prospective et à l’effectif requis par cellule expérimentale. Une règle empirique minimale, solidement ancrée dans la littérature psychométrique, stipule que le nombre d’observations par cellule (groupe) doit impérativement excéder le nombre total de variables dépendantes ($n_j > p$). Si cette condition fondamentale n’est pas satisfaite, la matrice de variance-covariance d’échantillon ne dispose pas d’un rang suffisant ; elle devient singulière, rendant l’exécution formelle du test impossible. Pour garantir une robustesse raisonnable des estimations des covariances, de nombreux auteurs préconisent un minimum absolu de $20$ à $30$ sujets par cellule expérimentale lorsque le modèle inclut 3 à 5 variables dépendantes.

Le calcul formel de puissance statistique a priori pour une MANOVA peut être conduit à l’aide de logiciels dédiés comme G*Power ou au moyen d’outils R spécialisés. Contrairement à l’ANOVA univariée où la taille d’effet est caractérisée par le $f$ de Cohen ou le $\eta^2$, la puissance multivariée dépend non seulement de l’écart spatial entre les centroïdes, mais également de la matrice de corrélation résiduelle et de la directionnalité de l’effet le long des vecteurs propres. En règle générale, la taille d’effet multivariée est définie par le $V$ de Cramer ou le $f^2$ multivarié dérivé de la trace de Pillai ou du Lambda de Wilks. Pour un effet moyen attendu ($f^2 = 0{,}15$), avec une puissance cible $1 – \beta = 0{,}80$ et un seuil $\alpha = 0{,}05$, le dimensionnement d’un plan à 3 groupes et 4 variables dépendantes requiert un effectif global d’au moins 120 à 150 participants.

L’impact du déséquilibre des effectifs entre les groupes (plans non orthogonaux) constitue un facteur aggravant qui compromet la puissance et la robustesse de la MANOVA. Lorsque les tailles de cellules diffèrent de manière substantielle ($n_{\max} / n_{\min} > 1{,}5$), la modélisation standard basée sur les sommes de carrés de Type I génère des résultats séquentiels biaisés dépendant de l’ordre d’introduction des variables. De surcroît, le déséquilibre des effectifs exacerbe considérablement la vulnérabilité du modèle face aux violations de l’hypothèse d’homogénéité des matrices de covariance, entraînant des taux de rejet empiriques imprévisibles.

3. Prérequis et hypothèses statistiques fondamentales

3.1 Indépendance des observations et plan d’échantillonnage

L’indépendance statistique des résidus constitue le postulat le plus déterminant de la modélisation linéaire multivariée. Cette hypothèse impose que l’erreur de mesure associée à l’observation $i$ n’entretienne aucune dépendance probabiliste avec l’erreur de mesure associée à l’observation $k$, que ce soit au sein d’une même cellule ou entre des cellules distinctes. La violation de cette condition altère l’estimation des erreurs types matricielles, ce qui se traduit presque invariablement par une inflation dramatique de l’erreur de type I, pouvant porter le taux de faux positifs effectif à 30 % ou 50 % au lieu des 5 % théoriques prévus par le seuil nominal.

Dans la pratique expérimentale, le respect de cette hypothèse repose sur la mise en œuvre de procédures d’échantillonnage probabiliste et sur l’assignation aléatoire des participants aux conditions expérimentales. Les biais d’interdépendance surgissent classiquement lors de l’évaluation de sujets regroupés par clusters spatiaux ou temporels – par exemple, des élèves au sein d’une même classe, des patients traités par le même thérapeute, ou des données recueillies en dyades relationnelles. Dans ces situations, la corrélation intra-classe attribuable au contexte partagé invalide l’hypothèse d’indépendance de la MANOVA standard. Si ces structures de dépendance existent, le chercheur doit se réorienter vers des modèles mixtes multivariés (multivariate linear mixed-effects models) ou la modélisation par équations structurelles multi-niveaux.

Il est également capital de préserver une distinction claire entre les variables indépendantes catégorielles, qui définissent les partitions d’échantillonnage, et les covariables continues potentielles. Les variables indépendantes doivent être assignées ou constatées sans rétroaction issue de la mesure du vecteur dépendant, sous peine d’induire une endogénéité structurelle qui compromettrait l’interprétabilité causale des inférences multivariées.

3.2 Normalité multivariée et détection des valeurs aberrantes

La MANOVA postule que le vecteur conjoint des résidus suit une distribution normale multivariée au sein de chaque groupe expérimental. La normalité multivariée impose non seulement que chaque variable dépendante affiche isolément une distribution normale univariée (avec une asymétrie et un aplatissement conformes aux attentes gaussiennes), mais requiert également que toute combinaison linéaire des variables dépendantes suive une loi normale univariée. Il s’ensuit que la vérification de la normalité sur chaque axe de manière disjointe est une condition nécessaire mais rigoureusement insuffisante pour attester de la normalité multivariée conjointe.

L’existence de valeurs aberrantes multivariées (multivariate outliers) représente une menace directe pour la validité du modèle. Un sujet peut présenter des valeurs parfaitement acceptables et plausibles sur chacune des échelles considérées isolément, tout en manifestant une configuration bivariée ou multivariée hautement aberrante compte tenu de la structure de corrélation générale. La métrique de référence pour quantifier cet éloignement géométrique est la distance de Mahalanobis ($D^2$), définie pour un individu $i$ par :

$$D_i^2 = (\mathbf{y}_i – \bar{\mathbf{y}})^top \mathbf{S}^{-1} (\mathbf{y}_i – \bar{\mathbf{y}})$$

où $\bar{\mathbf{y}}$ représente le vecteur moyen et $\mathbf{S}^{-1}$ est l’inverse de la matrice de variance-covariance de l’échantillon. En l’absence de déviation par rapport à la normalité multivariée, les distances de Mahalanobis se distribuent selon une loi du Chi-deux avec un nombre de degrés de liberté égal au nombre de variables dépendantes ($\chi^2_p$). Une observation présentant une distance de Mahalanobis associée à une probabilité d’occurrence inférieure à $p < 0{,}001$ est généralement qualifiée de valeur aberrante multivariée.

Sur le plan formel, l’évaluation de la normalité multivariée repose sur des tests probabilistes robustes, notamment le test d’asymétrie et d’aplatissement de Mardia. L’asymétrie multivariée de Mardia quantifie la projection tridimensionnelle des moments d’ordre 3, tandis que le coefficient d’aplatissement évalue le comportement des queues de distribution multidimensionnelles. Des déviations substantielles de la normalité multivariée affectent la taille réelle des tests d’hypothèse multivariés, en particulier lorsque la non-normalité découle de la présence de queues lourdes ou d’observations influentes extrêmes.

3.3 Homogénéité des matrices de variance-covariance

L’homogénéité des matrices de variance-covariance à travers les différentes strates de la variable indépendante constitue l’équivalent multivarié du postulat d’homoscédasticité de l’ANOVA classique. Cette exigence impose que chaque groupe $j$ partage la même matrice de dispersion théorique :

$$boldsymbol{\Sigma}_1 = boldsymbol{\Sigma}_2 = dots = boldsymbol{\Sigma}_g = boldsymbol{\Sigma}$$

Cette égalité signifie non seulement que les variances univariées de chaque variable dépendante doivent être identiques d’un groupe à l’autre, mais également que les coefficients de covariance (et donc les corrélations) reliant les variables deux à deux doivent être strictement constants à travers l’ensemble des conditions expérimentales.

Le test statistique standard utilisé pour évaluer cette égalité est le test $M$ de Box. Ce test exploite un ratio fondé sur les déterminants des matrices de variance-covariance des sous-groupes comparés au déterminant de la matrice de covariance résiduelle poolée :

$$M = (N – g) \ln |\mathbf{S}_p| – \sum_{j=1}^g (n_j – 1) \ln |\mathbf{S}_j|$$

où $\mathbf{S}_p$ est la matrice de covariance groupée et $\mathbf{S}_j$ celle du groupe $j$. Toutefois, le test $M$ de Box se caractérise par une hypersensibilité notoire aux violations de la normalité multivariée. Sous des conditions de non-normalité modérée, le test $M$ rejette fréquemment l’hypothèse nulle d’homogénéité à tort, identifiant une fausse hétéroscédasticité qui reflète simplement l’asymétrie des distributions.

Les conséquences d’une violation de cette hypothèse dépendent directement de l’équilibre des tailles de groupes. Lorsque les effectifs sont rigoureusement équilibrés ($n_1 = n_2 = dots = n_g$), les critères de décision de la MANOVA manifestent une robustesse remarquable, limitant la distorsion du taux d’erreur de type I. À l’opposé, en présence d’un plan déséquilibré, la violation du postulat induit des distorsions asymétriques préjudiciables : si les groupes aux effectifs les plus réduits affichent les variances et covariances les plus larges, le test devient excessivement libéral (taux d’erreur de type I réel très supérieur au seuil nominal). Si les groupes à grand effectif présentent les variances les plus importantes, le test devient excessivement conservateur, causant un effondrement de la puissance statistique.

3.4 Absence de multicolinéarité parfaite et linéarité bivariée

La validité des solutions géométriques générées par la MANOVA suppose que les relations entretenues par les variables dépendantes deux à deux soient approximativement linéaires au sein de chaque groupe expérimental. La modélisation linéaire multivariée cherche à optimiser la séparation des centroïdes via des combinaisons linéaires des indicateurs observés. Si la trajectoire bivalente unissant deux variables présente une courbure curvilinéaire prononcée (par exemple, une fonction quadratique ou exponentielle), la projection matricielle échouera à capturer la totalité de l’interdépendance partagée, affaiblissant la capacité du modèle à isoler l’effet de groupe.

De plus, l’absence de multicolinéarité parfaite ou quasi-parfaite demeure un impératif technique absolu. La multicolinéarité parfaite se produit lorsque l’une des variables dépendantes n’est rien d’autre qu’une combinaison arithmétique directe d’autres variables présentes dans l’analyse (par exemple, inclure simultanément les scores à trois sous-échelles d’un test psychométrique et le score global calculé par leur somme directe). Dans un tel scénario, le rang mathématique de la matrice $\mathbf{Y}$ devient déficient ($rg(\mathbf{Y}) < p$), et le déterminant de la matrice résiduelle s'annule rigoureusement ($|\mathbf{E}| = 0$). L'inversion matricielle devient strictement impossible.

Pour prévenir cet écueil, l’analyste doit contrôler le déterminant de la matrice de corrélation résiduelle et examiner les facteurs d’inflation de la variance (VIF – Variance Inflation Factor) adaptés au cadre multivarié. Des valeurs de VIF excédant le seuil conventionnel de 5 ou 10 signalent une redondance structurelle incompatible avec l’exploitation sécurisée des algorithmes de décomposition spectrale.

4. Préparation de l’environnement de travail et importation des données dans R

4.1 Sélection et installation des packages spécialisés

L’écosystème open-source R propose une architecture modulaire particulièrement performante pour le calcul matriciel, l’estimation des modèles multivariés et les représentations graphiques diagnostiques. Pour mener à bien l’ensemble du pipeline analytique d’une MANOVA, la mise en place d’un environnement intégrant plusieurs packages complémentaires est nécessaire. Le socle moderne de manipulation de données s’articule autour du métapackage tidyverse, dont les composants dplyr et tidyr simplifient le prétraitement matriciel, et ggplot2 qui gère les représentations graphiques personnalisées.

Pour la modélisation statistique avancée, le package car (Companion to Applied Regression) développé par John Fox constitue un outil indispensable. Contrairement à la fonction native de R qui privilégie les sommes de carrés de type I, le package car intègre la fonction Manova() qui permet d’extraire les sommes de carrés de Type II et de Type III, cruciales pour traiter les plans factoriels déséquilibrés. Le package MASS fournit quant à lui les fonctions requises pour l’analyse discriminante linéaire complémentaire via lda().

L’évaluation des hypothèses probabilistes et la visualisation diagnostique nécessitent des bibliothèques dédiées aux données multidimensionnelles. Le package MVN constitue l’outil le plus complet pour exécuter les tests de normalité multivariée (Mardia, Royston, Henze-Zirkler). Le package biotools renferme l’algorithme standard du test $M$ de Box avec la fonction boxM(). Enfin, le package heplots fournit une suite de fonctions graphiques projetant visuellement les ellipsoïdes de variance d’hypothèse et d’erreur dans les espaces à deux et trois dimensions. Les extensions psych et broom facilitent quant à elles l’extraction des statistiques descriptives et la conversion des sorties d’analyse en structures de données ordonnées.

4.2 Structuration et typage du jeu de données

L’exécution d’une MANOVA dans R exige que le jeu de données d’origine soit formaté selon une structure rectangulaire stricte en format large (wide format). Dans cette disposition tabulaire, chaque ligne doit représenter de manière univoque une unité expérimentale unique (un participant indépendant), tandis que chaque colonne encode soit un prédicteur factoriel, soit une covariable continue, soit l’une des composantes numériques du vecteur dépendant. Cette organisation tranche avec le format long (long format) traditionnellement requis pour les modèles mixtes ou les ANOVAs univariées à mesures répétées.

Le typage explicite des colonnes au sein du data.frame ou du tibble R constitue une étape préliminaire déterminante. Les variables indépendantes doivent impérativement être instanciées en tant que facteurs (factor), et non comme des vecteurs numériques ou des chaînes de caractères brutes (character). L’encodage explicite via les commandes factor() ou as.factor() garantit que le moteur linéaire construira correctement la matrice de contrastes expérimentaux $\mathbf{X}$ (par exemple via contr.sum ou contr.treatment). Les variables dépendantes doivent obligatoirement présenter le type numérique à virgule flottante (numeric ou double).

La gestion rigoureuse des données manquantes (NA) s’avère critique en analyse multivariée. Par défaut, les fonctions R exécutent une suppression cas par cas (listwise deletion). Dès lors qu’une seule valeur dépendante manque sur un vecteur de 5 dimensions, le participant est intégralement exclu de l’estimation de l’ensemble des matrices $\mathbf{H}$ et $\mathbf{E}$. Cette attrition peut réduire l’échantillon utile de façon démesurée. Le chercheur doit évaluer le patron des données manquantes à l’aide de packages dédiés (tels que mice ou VIM). Si le mécanisme de perte respecte le postulat de données manquantes aléatoires (MAR – Missing At Random), le recours à l’imputation multiple par équations chaînées avant la modélisation multivariée doit être systématiquement envisagé pour prévenir l’érosion de la puissance statistique.

4.3 Constitution de la matrice des variables dépendantes

Dans la syntaxe du langage R, l’évaluation conjointe de plusieurs variables dépendantes au sein de la fonction modèle requiert la création d’un objet matriciel à l’aide de l’opérateur de concaténation de colonnes cbind(). Cet opérateur assemble les colonnes sélectionnées au sein d’une matrice bidimensionnelle de dimensions $N \times p$. Par exemple, si l’on évalue l’inhibition motrice (Inhibition), la flexibilité attentionnelle (Flexibility) et la mise à jour mnésique (Updating), la matrice de réponse se déclare sous la forme :

# Définition formelle de la matrice de variables dépendantes
Y <- cbind(donnees$Inhibition, donnees$Flexibility, donnees$Updating)

La question de la standardisation préalable des échelles de mesure constitue un arbitrage méthodologique important. La MANOVA repose théoriquement sur le principe d’invariance par transformation affine : l’application d’une multiplication par une constante ou l’addition d’un scalaire sur l’une des variables dépendantes ne modifie aucunement la valeur propre $lambda$ de $\mathbf{E}^{-1}\mathbf{H}$, ni les statistiques de test canoniques (Wilks, Pillai). Par conséquent, transformer les variables en scores centrés-réduits ($z$-scores) au moyen de la fonction scale() ne modifie pas les valeurs de $p$ multivariées.

Néanmoins, la standardisation préalable offre des avantages substantiels lors des étapes d’analyse subséquentes. Elle harmonise l’échelle d’affichage des graphiques bivariés, élimine les distorsions visuelles causées par des amplitudes d’échelles divergentes (par exemple, un temps de réaction oscillant entre 400 et 900 ms comparé à un taux d’exactitude variant de 0 à 1), et simplifie grandement l’interprétation des coefficients discriminants standardisés lors de l’analyse canonique post-hoc.

5. Exploration visuelle et analyse préliminaire des données multivariées

5.1 Statistiques descriptives groupées et bivariées

L’exploration analytique des données multivariées débute systématiquement par l’estimation des métriques de tendance centrale et de dispersion séparément pour chaque cellule du plan expérimental. À l’aide du package dplyr, l’analyste calcule de manière vectorielle les moyennes empiriques, les écarts-types, ainsi que les intervalles de confiance à 95 % pour chaque dimension dépendante à travers les niveaux du facteur inter-sujets. Cette première synthèse numérique permet de vérifier l’absence d’incohérences de saisie, de vérifier les bornes théoriques de chaque test psychométrique et d’identifier de prime abord les variables qui semblent manifester des contrastes univariés perceptibles.

Dans un second temps, il est indispensable de générer la matrice des corrélations bivariées de Pearson sur l’ensemble de l’échantillon, puis au sein de chaque groupe expérimental séparément. L’examen des matrices de corrélation permet de vérifier si les associations empiriques s’inscrivent bien dans l’intervalle méthodologiquement pertinent (typiquement entre $r = 0{,}30$ et $r = 0{,}70$). Une corrélation intra-groupe stable entre deux mesures confirme la pertinence de leur agrégation dans un vecteur dépendant conjoint. La fonction corr.test() du package psych fournit simultanément la matrice des coefficients de Pearson et la matrice associée des niveaux de significativité corrigés.

L’évaluation des statistiques descriptives groupées offre également une première fenêtre diagnostique sur l’homogénéité des variances univariées. Si le rapport entre la variance maximale et la variance minimale pour une variable dépendante donnée excède un ratio de 4:1 (soit un ratio d’écarts-types supérieur à 2:1) entre les cellules, l’analyste doit redoubler de prudence, particulièrement si ce déséquilibre de variance coïncide avec des effectifs de groupes inégaux.

5.2 Visualisations bivariées et projections multidimensionnelles

La représentation graphique préliminaire en contexte multivarié ne peut se cantonner aux traditionnels diagrammes en barres ou boîtes à moustaches univariées. L’outil diagnostique de prédilection est la matrice de nuages de points bivariés (pairs plot), générée par des fonctions telles que pairs() du système de base de R, ggpairs() du package GGally, ou la fonction spécialisée scatterplotMatrix() du package car. Cet agencement affiche sur sa diagonale les densités de distribution univariées pour chaque groupe, tandis que les cellules hors-diagonale projettent les nuages de dispersion bivariés croisant chaque paire de variables dépendantes avec un code couleur identifiant les modalités expérimentales.

L’incorporation d’ellipses de confiance à 95 % superposées sur les diagrammes de dispersion bivariés constitue un moyen efficace d’appréhender la forme et l’orientation des nuages de points. La géométrie de ces ellipses reflète directement la matrice de covariance intra-groupe : la longueur des axes principaux traduit l’amplitude des variances, tandis que l’angle d’inclinaison de l’ellipse matérialise la direction et l’intensité de la corrélation bivariée. Une superposition quasi-totale des ellipses de différents groupes signale une absence d’effet de séparation, tandis qu’un décalage spatial prononcé de leurs centres matérialise une divergence des centroïdes.

L’inspection attentive de ces nuages de points permet d’identifier visuellement les ruptures de linéarité bivariée ou la présence d’hétéroscédasticité multivariée. Si l’ellipse d’un groupe expérimental apparaît circulaire tandis que celle d’un autre groupe présente un profil hyper-elliptique fortement étiré et oblique, cela indique que la structure de corrélation fluctue sous l’effet du traitement expérimental, présageant un rejet potentiel du postulat d’homogénéité des matrices de covariance.

6. Vérification empirique des postulats sous R

6.1 Contrôle des valeurs aberrantes multivariées sous R

La recherche de valeurs aberrantes multivariées requiert le calcul formel des distances de Mahalanobis pour chaque observation par rapport au centroïde multivarié moyen. Bien que la fonction native mahalanobis() de R permette cette opération, elle doit être calculée de préférence sur les matrices résiduelles ou au sein de chaque groupe expérimental pris individuellement, afin que la distance d’un individu soit évaluée par rapport au centroïde de son propre groupe et non par rapport à la moyenne générale de l’échantillon global.

La syntaxe de programmation pour isoler et identifier les observations aberrantes au sein d’un jeu de données s’articule comme suit :

# Calcul de la distance de Mahalanobis groupe par groupe
calculer_mahalanobis <- function(df, variables, groupe) {
 df %>%
 dplyr::group_by(across(all_of(groupe))) %>%
 dplyr::mutate(
 d2 = mahalanobis(
 x = pick(all_of(variables)),
 center = colMeans(pick(all_of(variables))),
 cov = cov(pick(all_of(variables)))
 ),
 p_val_d2 = pchisq(d2, df = length(variables), lower.tail = FALSE)
 ) %>%
 dplyr::ungroup()
}

Le seuil décisionnel canonique d’exclusion ou d’investigation approfondie repose sur un niveau de significativité très conservateur de $\alpha = 0{,}001$ appliqué à la distribution $\chi^2$ à $p$ degrés de liberté. Toute observation affichant une probabilité $p_{\text{val_d2}} < 0{,}001$ doit faire l'objet d'une analyse clinique et méthodologique détaillée. L'exclusion d'un participant ne doit jamais être automatique ; le chercheur doit consigner rigoureusement s'il s'agit d'une erreur d'enregistrement, d'une anomalie procédurale ou d'une variabilité phénotypique réelle mais extrême.

6.2 Test de la normalité multivariée avec le package MVN

L’évaluation inférentielle de la normalité multivariée dans l’environnement R est grandement optimisée par l’utilisation du package MVN. La fonction centrale mvn() centralise les algorithmes d’estimation de Mardia, de Royston et de Henze-Zirkler, tout en produisant simultanément les tests univariés de Shapiro-Wilk pour chaque composante vectorielle :

library(MVN)
# Exécution des tests de normalité multivariée par groupe expérimental
resultats_mvn <- mvn(
 data = donnees[, c("Inhibition", "Flexibility", "Updating")],
 subset = donnees$Groupe,
 mvnTest = "mardia",
 univariateTest = "SW",
 multivariatePlot = "qq"
)

Le test d’asymétrie multivariée de Mardia évalue si les déformations obliques du nuage multidimensionnel dépassent les fluctuations d’échantillonnage attendues sous la loi normale, tandis que le test d’aplatissement de Mardia évalue la lourdeur des queues de distribution conjointe. L’obtention d’une valeur de $p > 0{,}05$ pour ces deux indices indique la conformité des données avec le postulat de normalité multivariée.

En complément des indicateurs numériques, la fonction génère un graphique quantile-quantile multivarié (Chi-square Q-Q plot). Ce diagnostic visuel projette les distances de Mahalanobis empiriques en fonction des quantiles théoriques issus de la distribution du Chi-deux. La normalité multivariée est validée lorsque les points empiriques s’alignent rigoureusement le long de la diagonale de référence à 45 degrés. Si les données s’écartent modérément de la linéarité mais que l’échantillon comporte plus de 30 observations par cellule sans asymétrie univariée sévère, le théorème central limite multivarié assure une robustesse relative des statistiques de test, notamment pour la trace de Pillai.

6.3 Évaluation de l’homogénéité des covariances avec le test M de Box

La vérification de l’égalité structurelle des matrices de variance-covariance s’effectue au moyen de la fonction boxM() accessible via le package biotools. Cette fonction prend en arguments la matrice des variables dépendantes et le facteur définissant les groupes expérimentaux :

library(biotools)
# Exécution formelle du test M de Box
test_box <- boxM(
 data = donnees[, c("Inhibition", "Flexibility", "Updating")],
 grouping = donnees$Groupe
)
print(test_box)

L’interprétation de la sortie du test $M$ de Box requiert une rigueur méthodologique particulière. En raison de sa sensibilité extrême à d’infimes déviations de la normalité multivariée sans impact pratique sur les taux d’erreur du modèle linéaire, les statisticiens recommandent de manière unanime d’adopter un seuil alpha ajusté hautement restrictif, fixé à $\alpha = 0{,}001$. Par conséquent, l’hypothèse nulle d’homogénéité des matrices de covariance n’est rejetée que si la $p$-valeur associée à la statistique $M$ est strictement inférieure à $0{,}001$.

Dans l’éventualité où le test $M$ de Box s’avère statistiquement significatif au seuil de $0{,}001$ :

• Si le plan d’expérience est strictement équilibré ($n_1 = n_2 = dots = n_g$), l’effet de cette violation sur les taux d’erreur de première espèce reste marginal. Il convient alors de privilégier systématiquement la Trace de Pillai-Bartlett, reconnue pour sa robustesse intrinsèque aux hétérogénéités de variance.
• Si le plan d’expérience présente un déséquilibre marqué des effectifs, l’interprétation de la MANOVA paramétrique classique devient hasardeuse. Il est alors impératif de se tourner vers des alternatives robustes, des méthodes de rééchantillonnage bootstrap ou la MANOVA non paramétrique basée sur des permutations (PERMANOVA).

7. Spécification et exécution du modèle MANOVA avec la fonction manova()

7.1 Syntaxe générale de la modélisation sous R

Au sein du système de base R, la déclaration formelle d’une analyse de variance multivariée repose sur la fonction manova(). La formule de modélisation exploite l’opérateur d’assignation tilde (~), en plaçant à sa gauche la matrice multidimensionnelle assemblée via cbind() et à sa droite la structure des prédicteurs catégoriels. Dans le cadre d’un plan factoriel à un facteur inter-sujets (one-way MANOVA), la spécification s’écrit :

# Spécification d'un modèle MANOVA à un facteur
modele_manova_simple <- manova(
 cbind(Inhibition, Flexibility, Updating) ~ Groupe,
 data = donnees
)

Cette architecture syntaxique s’étend aisément aux plans factoriels à plusieurs variables indépendantes (two-way MANOVA, three-way MANOVA). Lorsqu’un protocole évalue simultanément les effets de deux facteurs inter-sujets (par exemple, le type d’intervention thérapeutique Traitement et la sévérité initiale du trouble Severite), le modèle doit intégrer les effets principaux ainsi que leur terme d’interaction multivarié :

# Modèle factoriel complet avec interaction multivariée
modele_manova_factoriel <- manova(
 cbind(Inhibition, Flexibility, Updating) ~ Traitement * Severite,
 data = donnees
)

L’opérateur astérisque (*) instancie automatiquement les termes simples Traitement et Severite, ainsi que leur produit d’interaction Traitement:Severite. Cette interaction multivariée évalue formellement si le vecteur de divergence spatiale induit par les différentes modalités de traitement varie en magnitude ou en direction selon le niveau de sévérité initiale des participants.

7.2 Choix des types de sommes de carrés (Type I, II ou III)

La fonction native manova() de R procède par défaut à une décomposition séquentielle des sommes de carrés et de produits croisés, correspondant aux sommes de carrés de Type I. Dans une décomposition de Type I, l’ordre de déclaration des variables indépendantes dans la formule modifie directement les matrices $\mathbf{H}$ associées à chaque terme. Le premier terme absorbe toute la variance explicative qu’il partage potentiellement avec les facteurs suivants, ces derniers n’étant testés que sur la variance résiduelle restante. Dans les plans d’expérience parfaitement orthogonaux (caractérisés par des effectifs rigoureusement identiques dans chaque cellule), les décompositions de Type I, II et III conduisent à des résultats numériques identiques.

Cependant, dans la quasi-totalité des recherches quasi-expérimentales ou cliniques, les effectifs des cellules s’avèrent inégaux. L’application des sommes de carrés de Type I devient alors inappropriée, car les conclusions statistiques dépendent arbitrairement de l’ordre d’écriture choisi par l’analyste. Il est donc indispensable d’utiliser la fonction Manova() (avec une majuscule initiale) fournie par le package car, qui implémente les décompositions de Type II et de Type III :

library(car)
# Spécification des contrastes orthogonaux indispensables au Type III
options(contrasts = c("contr.sum", "contr.poly"))
# Exécution de la MANOVA factorielle avec sommes des carrés de Type III
modele_manova_type3 <- Manova(
 lm(cbind(Inhibition, Flexibility, Updating) ~ Traitement * Severite, data = donnees),
 type = "III"
)

Dans cette syntaxe, les sommes de carrés de Type II testent chaque effet principal en contrôlant pour les autres effets principaux, mais sans ajuster pour l’interaction d’ordre supérieur. Les sommes de carrés de Type III testent chaque terme (effet principal ou interaction) conditionnellement à l’ensemble des autres termes du modèle complet. Le Type III requiert impérativement de paramétrer les contrastes sous forme orthogonale centrée (contr.sum) avant l’estimation du modèle linéaire sous-jacent lm(), sous peine de biaiser l’interprétation des intercepts et des effets marginaux multivariés.

7.3 Gestion des covariables continues (MANCOVA)

L’analyse de covariance multivariée (MANCOVA) représente une extension mathématique directe de la MANOVA, permettant d’intégrer une ou plusieurs variables prédictives continues au sein du modèle expérimental. Ces variables, qualifiées de covariables, ont pour vocation d’ajuster statistiquement les centroïdes de groupes en éliminant la fraction de variance résiduelle attribuable à des facteurs de confusion préexistants non contrôlés (tels que l’âge chronologique, le quotient intellectuel global ou le niveau symptomatique pré-test). Cette réduction du bruit intra-groupe accroît la puissance statistique du test de séparation spatiale.

L’implémentation d’une MANCOVA exige au préalable la validation d’une hypothèse fondamentale : l’homogénéité des pentes de régression multivariées. Ce postulat stipule que la relation linéaire multivariée unissant les covariables continues au vecteur de variables dépendantes doit demeurer rigoureusement constante à travers l’ensemble des groupes expérimentaux. La violation de cette condition indique une interaction entre la covariable et le traitement expérimental, invalidant le principe d’ajustement global.

Sous R, l’évaluation de cette hypothèse s’effectue en testant le terme d’interaction entre le facteur de groupe et la covariable continue :

# Évaluation de l'homogénéité des pentes de régression multivariées
test_homogeneite_pentes <- manova(
 cbind(Inhibition, Flexibility, Updating) ~ Groupe * Age,
 data = donnees
)
summary(test_homogeneite_pentes, test = "Wilks")

Si le terme d’interaction multivarié Groupe:Age s’avère non significatif ($p > 0{,}05$), l’hypothèse de parallélisme des pentes multivariées est retenue. L’analyste spécifie alors le modèle MANCOVA définitif en retirant le terme d’interaction, en plaçant la covariable continue en tête de formule :

# Modèle MANCOVA ajusté définitif
modele_mancova <- Manova(
 lm(cbind(Inhibition, Flexibility, Updating) ~ Age + Groupe, data = donnees),
 type = "II"
)
summary(modele_mancova, test = "Pillai")

8. Interprétation approfondie des statistiques de test multivariées

8.1 Les quatre critères multivariés canoniques

La statistique univariée $F$ de Fisher repose sur le simple quotient de deux variances scalaires ($MS_{\text{inter}} / MS_{\text{intra}}$). En revanche, la comparaison des matrices multidimensionnelles $\mathbf{H}$ et $\mathbf{E}$ ne peut se résumer à un ratio univoque. Les statisticiens ont élaboré quatre critères multivariés canoniques fondés sur les valeurs propres non nulles $\lambda_i$ de l’opérateur linéaire $\mathbf{E}^{-1}\mathbf{H}$, où $i = 1, dots, s$ et $s = min(p, g – 1)$.

1. La Trace de Pillai-Bartlett ($V$) : Elle est définie comme la somme des variances expliquées le long de l’ensemble des axes canoniques discriminants :

$$V = operatorname{tr}\left( \mathbf{H}(\mathbf{H} + \mathbf{E})^{-1} \right) = \sum_{i=1}^s \frac{\lambda_i}{1 + \lambda_i}$$

Sa valeur oscille entre $0$ et $s$. Plus $V$ est élevé, plus la divergence spatiale globale entre les centroïdes de groupes est marquée relativement au bruit résiduel. En raison de sa structure de sommation pondérée, la trace de Pillai est reconnue comme la statistique la plus robuste face aux déviations modérées de la normalité multivariée et aux hétérogénéités des matrices de variance-covariance.

2. Le Lambda de Wilks ($Lambda$) : Reposant sur le principe du rapport de vraisemblance (likelihood-ratio criterion), le Lambda de Wilks mesure la proportion de variance non expliquée par l’effet expérimental :

$$\Lambda = \frac{|\mathbf{E}|}{|\mathbf{H} + \mathbf{E}|} = \prod_{i=1}^s \frac{1}{1 + \lambda_i}$$

La valeur de $Lambda$ est strictement bornée dans l’intervalle $[0, 1]$. Une valeur proche de $1$ traduit une absence totale d’effet de groupe (matrice $\mathbf{H}$ négligeable devant $\mathbf{E}$), tandis qu’une valeur tendant vers $0$ indique que la dispersion inter-groupes supplante largement la dispersion intra-groupe. Historiquement, le Lambda de Wilks constitue le critère le plus largement publié dans les revues internationales de psychologie.

3. La Trace de Hotelling-Lawley ($T$) : Cette statistique représente la somme arithmétique directe des valeurs propres de la matrice $\mathbf{E}^{-1}\mathbf{H}$ :

$$T = operatorname{tr}(\mathbf{E}^{-1}\mathbf{H}) = \sum_{i=1}^s \lambda_i$$

Directement proportionnelle à la variance totale inter-groupes projetée sur le sous-espace résiduel, cette mesure s’avère particulièrement sensible lorsque la séparation des groupes se déploie de manière homogène le long de plusieurs dimensions canoniques distinctes.

4. La Plus Grande Racine de Roy ($\theta$) : Contrairement aux trois métriques précédentes qui intègrent l’ensemble du spectre des valeurs propres, le critère de Roy retient exclusivement la valeur propre dominante $\lambda_{\max}$ :

$$\theta = \frac{\lambda_1}{1 + \lambda_1}$$

Le critère de Roy offre la puissance statistique la plus élevée parmi les quatre indices lorsque la séparation expérimentale réelle s’aligne rigoureusement sur une unique dimension canonique ($s = 1$). Cependant, elle perd considérablement en puissance si les centroïdes divergent selon une structure multidirectionnelle diffuse, et se révèle excessivement sensible aux violations de normalité.

8.2 Critères de sélection de la statistique adaptée au contexte de recherche

Le choix d’un critère multivarié parmi les quatre options disponibles dépend de la configuration structurelle des données et des violations éventuelles des postulats statistiques fondamentaux. Dans les conditions de validité idéales (normalité multivariée parfaite, matrices de variance-covariance rigoureusement homogènes et effectifs de groupes équilibrés), les quatre statistiques convergent asymptotiquement vers des conclusions inférentielles rigoureusement identiques.

Néanmoins, la réalité des données psychologiques et médicales confrontant fréquemment l’analyste à des imperfections d’échantillonnage, des règles d’arbitrage méthodologiques claires doivent être appliquées :

• Lorsque le jeu de données présente des tailles d’échantillons inégales, des matrices de variance-covariance modérément hétérogènes (test $M$ de Box significatif), ou de légères asymétries multivariées, la Trace de Pillai-Bartlett doit être systématiquement retenue. Des études de simulation de Monte-Carlo démontrent que Pillai maintient son taux d’erreur de première espèce au plus proche du seuil théorique de 5 % sous des conditions adverses où les autres critères deviennent soit trop libéraux, soit excessivement conservateurs.
• En situation de parfaite conformité aux postulats théoriques et d’effectifs équilibrés, le Lambda de Wilks demeure le standard académique conventionnel. Il offre une comparabilité directe avec les travaux antérieurs de la littérature scientifique et dispose d’une interprétation intuitive en tant que ratio de variance non expliquée.
• La Plus Grande Racine de Roy est à réserver aux configurations théoriques strictes où le chercheur postule à l’avance que la divergence entre les groupes s’articulera exclusivement le long d’un axe discriminant unique. En dehors de ce cas très spécifique, son utilisation comme preuve inférentielle est déconseillée en raison de sa fragilité face aux violations des postulats.

Toutes ces métriques multivariées sont converties par des approximations mathématiques – principalement l’approximation de Rao pour le Lambda de Wilks et l’approximation de Pillai pour la trace – en statistiques de test distribuées selon la loi $F$ de Fisher, permettant l’extraction d’une $p$-valeur opérationnelle.

8.3 Extraction et lecture critique de la sortie de summary.manova()

L’extraction des résultats inférentiels sous R s’opère au moyen de la méthode générique summary() appliquée à l’objet retourné par la modélisation multivariée. Par défaut, R affiche la Trace de Pillai, mais l’analyste peut spécifier explicitement le critère canonique souhaité via l’argument test :

# Extraction des résultats multivariés avec les différents critères
summary(modele_manova_simple, test = "Pillai")
summary(modele_manova_simple, test = "Wilks")
summary(modele_manova_simple, test = "Hotelling-Lawley")
summary(modele_manova_simple, test = "Roy")

L’exécution de cette commande génère une table inférentielle structurée comportant plusieurs colonnes clés. Prenons l’exemple d’une sortie typique obtenue avec le Lambda de Wilks :

 Df Wilks approx F num Df den Df Pr(>F) 
Groupe 2 0.58412 8.4125 6 190 2.145e-08 ***
Residuals 97 

L’évaluation rigoureuse de cette table implique une lecture attentive des paramètres :

• La colonne Df relative au terme inter-sujets indique les degrés de liberté factoriels ($g – 1$).
• La valeur numérique du critère (ici $\Lambda = 0{,}58412$) indique qu’environ 58 % de la variance multivariée globale n’est pas expliquée par le facteur de groupe, ce qui implique que l’effet expérimental capture plus de 41 % de la variance généralisée.
• La colonne approx F indique la valeur de la statistique $F$ de Fisher calculée par l’approximation de Rao (ici $F = 8{,}41$).
• Les degrés de liberté du numérateur (num Df) et du dénominateur (den Df) reflètent la dimensionnalité du problème mathématique. Le numérateur ($6$) résulte du produit du nombre de dimensions dépendantes par les degrés de liberté du groupe ($p \times (g – 1) = 3 \times 2$). Le dénominateur ($190$) correspond au pool de degrés de liberté résiduels matriciels.
• La probabilité de dépassement critique (Pr(>F)) est très inférieure à $0{,}001$. L’hypothèse nulle d’égalité des centroïdes multivariés dans l’espace tridimensionnel est donc rejetée avec une grande confiance statistique.

9. Analyses post-hoc et décomposition univariée après MANOVA

9.1 Décomposition univariée via summary.aov()

Le rejet de l’hypothèse nulle multivariée atteste qu’il existe une séparation globale statistiquement significative entre les centroïdes des groupes au sein du sous-espace multidimensionnel. Toutefois, ce résultat global ne renseigne pas sur les dimensions spécifiques qui portent cet effet. L’étape analytique subséquente classique consiste à décomposer le modèle multivarié en une série d’analyses de variance univariées d’accompagnement (follow-up ANOVAs).

Sous R, cette décomposition s’exécute de manière automatisée en appliquant la fonction summary.aov() à l’objet MANOVA préexistant :

# Décomposition de la MANOVA en ANOVAs univariées disjointes
deballage_univarie <- summary.aov(modele_manova_simple)
print(deballage_univarie)

Cette commande génère un bloc de sortie distinct pour chacune des variables dépendantes initialement intégrées dans l’objet cbind(). Chaque section présente la table d’ANOVA univariée correspondante, incluant la somme des carrés inter-groupes, la somme des carrés résiduelle, la valeur univariée de $F$, et sa $p$-valeur associée.

Cette démarche en deux temps (MANOVA globale suivie d’ANOVAs univariées) a historiquement été désignée sous le nom d’« approche de protection de Hummel-Sligo ». Cette protection n’élimine toutefois pas complètement le risque d’inflation de l’erreur de Type I lors de l’examen individuel des multiples variables dépendantes. Pour préserver l’intégrité de la décision inférentielle globale, l’analyste doit ajuster le seuil de significativité $\alpha$ appliqué à chaque test univarié d’accompagnement. La méthode de Bonferroni préconise d’utiliser un seuil ajusté $alpha’ = \alpha / p$. Pour 3 variables dépendantes et un seuil initial de $0{,}05$, une ANOVA univariée ne sera ainsi considérée comme significative que si sa $p$-valeur est strictement inférieure à $0{,}0167$. L’alternative séquentielle de Holm-Bonferroni peut également être mobilisée pour réduire la sévérité de cette correction sans sacrifier le contrôle du taux d’erreur par famille.

9.2 Comparaisons par paires et contrastes planifiés

Lorsqu’une variable dépendante franchit avec succès le filtre du seuil univarié ajusté et que la variable indépendante comprend au moins trois modalités expérimentales ($g ge 3$), il demeure impossible de statuer sur les paires de groupes à l’origine de la divergence sans procéder à des comparaisons multiples par paires (tests post-hoc) ou à des contrastes planifiés orthogonaux.

Le package emmeans (Estimated Marginal Means) constitue l’environnement le plus puissant et flexible sous R pour estimer les moyennes marginales ajustées et appliquer les tests de comparaison multiple post-hoc :

library(emmeans)
# Extraction des moyennes ajustées et comparaisons de Tukey pour chaque variable
emms_inhibition <- emmeans(modele_manova_simple, ~ Groupe, yvar = "Inhibition")
pairs(emms_inhibition, adjust = "tukey")
emms_flexibility <- emmeans(modele_manova_simple, ~ Groupe, yvar = "Flexibility")
pairs(emms_flexibility, adjust = "tukey")

L’argument adjust = "tukey" applique la méthode de Tukey HSD (Honestly Significant Difference), adaptée pour tester toutes les comparaisons deux à deux tout en maintenant rigoureusement le risque global d’erreur de première espèce au niveau fixé par l’expérimentateur. Si le chercheur dispose d’hypothèses théoriques a priori directes (par exemple, postuler qu’un groupe clinique expérimental surpasse conjointement deux groupes contrôles distincts), il est recommandé de spécifier des contrastes planifiés personnalisés plutôt que d’explorer l’ensemble des paires possibles, ce qui optimise la puissance statistique de l’analyse.

9.3 Analyse discriminante descriptive complémentaire

Bien que la décomposition univariée via summary.aov() demeure la stratégie la plus fréquente dans les publications, elle présente la limite conceptuelle de fractionner une structure multivariée en fragments univariés indépendants, occultant à nouveau les covariances croisées. L’analyse discriminante linéaire (LDA – Linear Discriminant Analysis) descriptive constitue le complément analytique le plus cohérent après une MANOVA significative. Elle cherche à identifier les combinaisons linéaires optimales des variables dépendantes qui maximisent la séparation statistique entre les groupes.

À l’aide de la fonction lda() du package MASS, l’analyste décompose directement l’espace vectoriel :

library(MASS)
# Modélisation discriminante descriptive sur le vecteur de variables dépendantes
modele_lda <- lda(Groupe ~ Inhibition + Flexibility + Updating, data = donnees)
print(modele_lda)

La sortie de cette analyse fournit les coefficients des fonctions discriminantes linéaires canoniques (LD1, LD2). Ces coefficients indiquent le poids relatif et la directionnalité de chaque indicateur dans la séparation spatiale maximale des centroïdes. L’examen des corrélations structurelles (ou charges canoniques), obtenues en calculant la corrélation de Pearson entre chaque variable dépendante brute et les scores discriminants individuels, renseigne précisément sur la dimensionnalité du phénomène sous-jacent :

# Extraction des scores discriminants individuels
scores_lda <- predict(modele_lda)$x
# Calcul des corrélations canoniques structurelles
cor_structurelles <- cor(donnees[, c("Inhibition", "Flexibility", "Updating")], scores_lda)
print(cor_structurelles)

Cette approche permet fréquemment de découvrir qu’une dimension univariée considérée comme non significative lors de l’étape summary.aov() contribue de manière substantielle à la fonction discriminante globale en ajustant ou en réprimant le bruit partagé avec d’autres échelles (effet de suppression multivarié).

10. Calcul et rapport de la taille d’effet multivariée et puissance statistique

10.1 Estimation de l’Eta-carré partiel multivarié

L’obtention d’une $p$-valeur hautement significative atteste de la reproductibilité probabiliste d’un effet, mais n’informe aucunement sur sa magnitude concrète ou sa pertinence clinique. Les standards méthodologiques internationaux imposent donc la communication systématique de la taille d’effet multivariée. L’indice standard le plus couramment mobilisé est l’Eta-carré partiel multivarié ($\eta_p^2$).

Lorsque l’inférence repose sur le Lambda de Wilks, l’Eta-carré partiel multivarié se calcule analytiquement au moyen de la formule de conversion suivante :

$$\eta_p^2 = 1 – \Lambda^{1 / s}$$

où $Lambda$ est la statistique observée de Wilks et $s = min(p, g – 1)$ représente le nombre maximal de dimensions canoniques indépendantes du modèle. Cet indice quantifie la proportion de variance généralisée totale attribuable au facteur inter-sujets, après élimination de la variance absorbée par les éventuels autres prédicteurs ou covariables du modèle.

Le package effectsize permet d’automatiser l’extraction standardisée de ces indices directement à partir des modèles multivariés sous R :

library(effectsize)
# Calcul automatisé des tailles d'effet multivariées partielles
tailles_effet <- eta_squared(modele_manova_simple, alternative = "two.sided")
print(tailles_effet)

L’interprétation qualitative de la magnitude de l’effet s’appuie généralement sur les repères conventionnels établis par Jacob Cohen pour le modèle linéaire, adaptés au contexte multidimensionnel :

• $\eta_p^2 \approx 0{,}01$ correspond à un effet de faible amplitude (environ 1 % de variance multivariée expliquée).
• $\eta_p^2 \approx 0{,}06$ correspond à un effet de magnitude moyenne (environ 6 % de variance multivariée expliquée).
• $\eta_p^2 ge 0{,}14$ correspond à un effet de forte amplitude (14 % ou plus de la variance multidimensionnelle totale absorbée par les groupes).

10.2 Intervalles de confiance des tailles d’effet

La communication d’une estimation ponctuelle de la taille d’effet ($\eta_p^2$) doit impérativement être consolidée par le rapport d’un intervalle de confiance probabiliste (typiquement à 90 % ou 95 %). En raison de la complexité algébrique de la distribution non centrale de Fisher et de la distribution de Wilks, la dérivation analytique exacte de ces intervalles est parfois délicate. Le recours aux procédures de rééchantillonnage non paramétrique (bootstrap) constitue l’approche méthodologique la plus robuste et la plus reproductible.

Le package boot permet d’estimer ces intervalles empiriques par rééchantillonnage itératif :

library(boot)
# Fonction d'estimation de l'Eta-carré partiel multivarié pour le bootstrap
boot_eta2 <- function(data, indices) {
 echantillon <- data[indices, ]
 mod <- manova(cbind(Inhibition, Flexibility, Updating) ~ Groupe, data = echantillon)
 wilks <- summary(mod, test = "Wilks")$stats[1, "Wilks"]
 s <- min(3, length(unique(data$Groupe)) - 1)
 eta2_p <- 1 - (wilks)^(1 / s)
 return(eta2_p)
}
# Lancement du rééchantillonnage sur 2000 itérations
set.seed(12345)
resultats_boot <- boot(data = donnees, statistic = boot_eta2, R = 2000)
ci_eta2 <- boot.ci(resultats_boot, type = "bca")
print(ci_eta2)

La précision de cet intervalle de confiance renseigne directement sur la stabilité des estimations. Un intervalle excessivement large (par exemple, $[0{,}02 ; 0{,}35]$) signale un manque de puissance statistique d’échantillonnage, invitant à la prudence quant à la robustesse des conclusions tirées de la taille d’effet ponctuelle.

10.3 Puissance observée versus puissance prospective

Une pratique encore répandue dans certaines revues empiriques consiste à calculer et rapporter la « puissance observée » (observed power ou post-hoc power) suite à l’obtention d’un résultat statistique, calculée à partir de la taille d’effet directement estimée sur l’échantillon. De nombreux statisticiens ont démontré que le concept de puissance post-hoc repose sur une tautologie mathématique sans valeur diagnostique : la puissance rétrospective n’est rien d’autre qu’une transformation monotone de la $p$-valeur observée. Un test affichant une valeur de $p = 0{,}049$ indiquera systématiquement une puissance rétrospective avoisinant 50 %, sans apporter d’information indépendante sur la qualité intrinsèque du protocole expérimental.

À l’inverse, l’analyse de puissance prospective planifiée constitue la seule démarche valide. Dans l’optique de concevoir de futures études de réplication ou de dimensionner adéquatement des extensions empiriques, le chercheur extrait les paramètres de non-centralité estimés à partir de la MANOVA actuelle pour paramétrer des simulations formelles. En exploitant les packages spécialisés comme pwr ou les scripts personnalisés de Monte-Carlo sous R, l’expérimentateur simule des matrices de dispersion semblables à celles observées empiriquement afin de déterminer le nombre exact de participants requis pour atteindre une puissance minimale de 80 % ou 90 % lors des travaux ultérieurs.

11. Diagnostic avancé des résidus et robustesse face aux violations d’hypothèses

11.1 Visualisation diagnostique avec les graphiques HE (Hypothesis-Error)

L’exploration diagnostique des modèles linéaires multivariés a été considérablement enrichie par le développement théorique et logiciel des graphiques Hypothèse-Erreur (HE plots), développés par Michael Friendly et intégrés dans le package R heplots. Contrairement aux visualisations marginales univariées classiques, les graphiques HE projettent simultanément les ellipsoïdes de covariance associés à l’effet du modèle (l’ellipsoïde de l’hypothèse $\mathbf{H}$) et les ellipsoïdes de covariance résiduelle (l’ellipsoïde d’erreur $\mathbf{E}$) pour des paires d’indicateurs dépendants.

La syntaxe de modélisation sous heplots s’interface directement avec les objets de modèle linéaire :

library(heplots)
# Modélisation linéaire multivariée sous-jacente
mod_lm <- lm(cbind(Inhibition, Flexibility, Updating) ~ Groupe, data = donnees)
# Génération du graphique HE bivarié
heplot(mod_lm, variables = c(1, 2),
 fill = TRUE, fill.alpha = 0.1,
 col = c("coral", "steelblue"),
 main = "Graphique HE : Inhibition vs Flexibilité")

L’interprétation géométrique de ce graphique est d’une grande clarté : l’ellipsoïde d’erreur $\mathbf{E}$ (généralement représenté en bleu) matérialise la dispersion résiduelle des observations après ajustement du modèle, son orientation diagonale reflétant la covariance résiduelle entre les variables. L’ellipsoïde d’hypothèse $\mathbf{H}$ (représenté en rouge ou orange) matérialise la dispersion entre les centroïdes des groupes. Dès lors que l’ellipsoïde de l’hypothèse $\mathbf{H}$ déborde de l’ellipsoïde résiduel $\mathbf{E}$ le long d’une direction spatiale donnée, l’effet multivarié est statistiquement significatif le long de cette orientation vectorielle.

Le package propose également une extension tridimensionnelle interactive via la fonction heplot3d(), qui permet de manipuler en temps réel dans l’espace les hyper-ellipsoïdes formés par trois variables dépendantes simultanées, facilitant la communication de configurations multivariées complexes auprès de la communauté scientifique.

11.2 Alternatives non paramétriques et approches par permutation

Lorsque les données empiriques violent simultanément la normalité multivariée et l’homogénéité des matrices de covariance, et que les effectifs par cellule sont inégaux ou limités, le cadre inférentiel paramétrique de la MANOVA standard perd de sa fiabilité. Dans ces circonstances, le chercheur doit mobiliser des alternatives non paramétriques, au premier rang desquelles figure l’analyse de variance multivariée par permutations (PERMANOVA), formalisée par Marti J. Anderson.

La PERMANOVA ne repose sur aucun postulat de normalité multivariée. Elle opère directement sur une matrice de distances géométriques calculée entre toutes les paires d’individus (distance euclidienne, distance de Bray-Curtis, distance de Mahalanobis ou distance de Manhattan). L’inférence statistique émerge de la permutation itérative (typiquement 999 à 9999 permutations aléatoires) des vecteurs d’observations entre les modalités expérimentales. Cette procédure est exécutée sous R via la fonction adonis2() du package écologique et statistique vegan :

library(vegan)
# Matrice de réponses multivariées
Y_mat <- as.matrix(donnees[, c("Inhibition", "Flexibility", "Updating")])
# Exécution d'une PERMANOVA basée sur la distance euclidienne (9999 permutations)
set.seed(42)
modele_permanova <- adonis2(
 Y_mat ~ Groupe,
 data = donnees,
 method = "euclidean",
 permutations = 9999
)
print(modele_permanova)

La sortie de adonis2() fournit un tableau d’analyse de variance basé sur les permutations, délivrant une statistique de pseudo-$F$ et une valeur de $p$ empirique hautement robuste, non sensible aux asymétries distributionnelles sous-jacentes. Si l’objectif principal est de s’affranchir de l’influence de données psychométriques asymétriques ou ordonnées sans recourir à des transformations arbitraires, la PERMANOVA constitue l’alternative non paramétrique de premier ordre.

11.3 Modélisation linéaire multivariée robuste

Une seconde alternative aux distorsions induites par la présence d’hétéroscédasticité ou de valeurs aberrantes non modélisées réside dans l’utilisation de la modélisation multivariée robuste fondée sur des estimateurs de position et de dispersion résistants aux contaminations d’échantillon. L’approche la plus reconnue exploite le déterminant minimal de covariance (MCD – Minimum Covariance Determinant) ou les $M$-estimateurs et $S$-estimateurs multivariés.

Le package spécialisé rrcov (Robust Robust Covariance) fournit une implémentation de ces techniques pour l’analyse multivariée, notamment via la fonction Wilks.test() robuste :

library(rrcov)
# Test multivarié robuste de Wilks fondé sur l'estimateur MCD
test_robuste <- Wilks.test(
 x = as.matrix(donnees[, c("Inhibition", "Flexibility", "Updating")]),
 grouping = as.factor(donnees$Groupe),
 method = "mcd"
)
print(test_robuste)

Cette méthode permet d’attribuer automatiquement un poids réduit aux observations marginales situées en périphérie de l’espace multivarié sans exiger leur suppression pure et simple de l’échantillon empirique. Elle évite ainsi la censure artificielle de la variabilité psychologique individuelle tout en assurant que les statistiques de test et les estimations des matrices de dispersion ne soient pas indûment influencées par une poignée d’individus atypiques.

12. Rédaction académique des résultats aux normes APA et bonnes pratiques

12.1 Format de transcription textuelle selon les normes APA (7e édition)

La communication scientifique des résultats d’une analyse de variance multivariée au sein des périodiques régis par l’American Psychological Association (APA 7e édition) exige une rigueur formelle stricte. Le texte principal doit rapporter de manière exhaustive le critère multivarié retenu (généralement la trace de Pillai ou le Lambda de Wilks), la valeur de l’approximation de Fisher, les degrés de liberté du numérateur et du dénominateur, le seuil de probabilité exact ($p$-valeur rapportée avec trois décimales sans zéro initial, sauf pour $p < 0{,}001$), ainsi que la taille d'effet multivariée ($\eta_p^2$).

Voici un gabarit canonique de rédaction textuelle intégrant l’ensemble de ces exigences académiques :

« Une analyse de variance multivariée à un facteur (MANOVA) a été conduite afin d’évaluer l’impact du type d’intervention thérapeutique (Groupe Contrôle, Thérapie Cognitive, Thérapie d’Acceptation et d’Engagement) sur l’ensemble vectoriel des trois fonctions exécutives mesurées : l’inhibition motrice, la flexibilité attentionnelle et la mise à jour mnésique. Les analyses préliminaires confirment l’absence de valeurs aberrantes multivariées significatives (distances de Mahalanobis, $p > 0{,}001$). L’hypothèse de normalité multivariée a été jugée recevable par le test de Mardia ($p > 0{,}05$). Bien que le test $M$ de Box ait révélé une violation mineure du postulat d’homogénéité des covariances ($M = 18{,}42, p = 0{,}008$), la parfaite équivalence de la taille des cellules ($n = 50$ par groupe) garantit la robustesse du modèle. En vertu de ces constats diagnostiques, la trace de Pillai-Bartlett a été retenue comme critère d’inférence décisionnel.

L’effet global du type d’intervention s’avère statistiquement significatif, $V = 0{,}42, F(6, 292) = 12{,}84, p < 0{,}001, \eta_p^2 = 0{,}21$, révélant une taille d'effet forte selon les critères conventionnels de Cohen. Afin de préciser la nature de cette différenciation, des analyses de variance univariées d'accompagnement ont été exécutées en appliquant une correction du seuil alpha de Bonferroni fixée à $alpha' = 0{,}05 / 3 = 0{,}017$. Une divergence inter-groupes significative émerge sur la flexibilité attentionnelle, $F(2, 147) = 18{,}73, p < 0{,}001, \eta_p^2 = 0{,}20$, ainsi que sur la capacité d'inhibition, $F(2, 147) = 8{,}45, p < 0{,}001, \eta_p^2 = 0{,}10$. En revanche, aucune modulation n'a été détectée pour la mise à jour mnésique, $F(2, 147) = 1{,}12, p = 0{,}328, \eta_p^2 = 0{,}01$. Les comparaisons par paires post-hoc de Tukey HSD indiquent que les participants du groupe ACT manifestent des scores de flexibilité significativement supérieurs à ceux du groupe Contrôle ($p < 0{,}001$) et du groupe Thérapie Cognitive ($p = 0{,}004$). »

12.2 Conception de tableaux synthétiques de résultats

L’intégration d’un tableau récapitulatif standardisé est vivement recommandée pour fluidifier la lecture et éviter l’alourdissement du corps de texte. Ce tableau doit synthétiser les moyennes observées, les écarts-types, les sommes des carrés univariées, les statistiques $F$, les seuils de significativité ajustés et les tailles d’effet pour chacune des dimensions étudiées.

Le package flextable ou modelsummary sous R permet de concevoir et d’exporter directement ces tableaux au format Word ou HTML en respectant scrupuleusement les exigences typographiques de l’APA (absence de bordures verticales, lignes horizontales d’en-tête et de bas de tableau uniquement) :

library(flextable)
# Exemple de construction d'un tableau récapitulatif APA
tableau_synthese <- data.frame(
 Variable = c("Inhibition", "Flexibilité", "Mise à jour"),
 Controle_M_SD = c("45.2 (8.1)", "52.1 (7.4)", "12.3 (2.8)"),
 TCC_M_SD = c("50.4 (7.9)", "58.3 (6.9)", "12.8 (3.1)"),
 ACT_M_SD = c("52.8 (8.5)", "67.4 (8.2)", "13.1 (2.9)"),
 F_univarie = c("8.45", "18.73", "1.12"),
 p_valeur = c("< .001", "< .001", ".328"),
 Eta2_partiel = c("0.10", "0.20", "0.01")
)
ft <- flextable(tableau_synthese) %>%
 set_header_labels(
 Variable = "Variable Dépendante",
 Controle_M_SD = "Contrôle (n=50)",
 TCC_M_SD = "TCC (n=50)",
 ACT_M_SD = "ACT (n=50)",
 F_univarie = "F (2, 147)",
 p_valeur = "p",
 Eta2_partiel = "η²p"
 ) %>%
 theme_apa() %>%
 autofit()

Ce format tabulaire assure que le lecteur dispose instantanément des métriques descriptives indispensables pour appréhender l’amplitude concrète des décalages de moyennes, tout en intégrant les inférences statistiques validées lors des tests post-hoc.

12.3 Bonnes pratiques de science ouverte et de reproductibilité

L’application rigoureuse des principes de la science ouverte (Open Science) constitue désormais une condition sine qua non de la crédibilité académique. Pour assurer la reproductibilité computationnelle intégrale d’une MANOVA, la simple description des étapes dans la section méthodologique d’un article ne suffit plus. L’analyste doit archiver l’intégralité du pipeline d’analyse au sein d’un script structuré, auto-suffisant et entièrement reproductible, idéalement documenté sous forme de document dynamique Quarto ou RMarkdown.

Plusieurs bonnes pratiques méthodologiques et de programmation doivent être rigoureusement appliquées :

• Fixation explicite des graines aléatoires : Dès lors qu’une étape d’imputation multiple, de bootstrap ou de permutation (comme dans le cas de adonis2()) intervient dans le script, la commande set.seed() doit être initialisée en amont afin que les calculs produisent des résultats numériques identiques sur toute machine.
• Gestion des environnements de dépendances : L’utilisation de packages de gestion d’environnement tels que renv permet de figer les versions exactes des packages R utilisés lors de l’analyse, prévenant toute rupture de compatibilité ascendante lors des mises à jour futures des bibliothèques.
• Dépôt public des données et du code : Les scripts de nettoyage, les scripts de modélisation et, dans la mesure où les considérations éthiques et médico-légales de confidentialité le permettent, les matrices de données brutes anonymisées doivent être déposés sur des registres ouverts pérennes bénéficiant d’identifiants numériques d’objets (DOI), tels que l’Open Science Framework (OSF) ou Zenodo.

En respectant scrupuleusement l’ensemble de ces préconisations techniques, théoriques et de communication, l’expérimentateur garantit que son recours à l’analyse de variance multivariée sous R constituera une démarche robuste, transparente, conforme à l’éthique de la recherche moderne et capable d’apporter un éclairage décisif sur la structure complexe des construits psychologiques et comportementaux.

Références

• Anderson, M. J. (2001). A new method for non-parametric multivariate analysis of variance. Austral Ecology, 26(1), 32–46. https://doi.org/10.1111/j.1442-9993.2001.01070.pp.x
• Box, G. E. P. (1949). A general distribution theory for a class of likelihood criteria. Biometrika, 36(3/4), 317–346. https://doi.org/10.2307/2332671
• Cohen, J. (1988). Statistical power analysis for the behavioral sciences (2nd ed.). Lawrence Erlbaum Associates.
• Fox, J., & Weisberg, S. (2019). An R companion to applied regression (3rd ed.). Sage Publications. https://socialsciences.mcmaster.ca/jfox/Books/Companion/
• Friendly, M. (2007). HE plots for multivariate linear models. Journal of Computational and Graphical Statistics, 16(2), 421–444. https://doi.org/10.1198/106186007X208407
• Huberty, C. J., & Olejnik, S. (2006). Applied MANOVA and discriminant analysis (2nd ed.). John Wiley & Sons. https://doi.org/10.1002/047178947X
• Korkmaz, S., Goksuluk, D., & Zararsiz, G. (2014). MVN: An R package for assessing multivariate normality. The R Journal, 6(2), 151–162. https://doi.org/10.32614/RJ-2014-031
• Mardia, K. V. (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika, 57(3), 519–530. https://doi.org/10.1093/biomet/57.3.519
• Pillai, K. C. S. (1955). Some new test criteria in multivariate analysis. The Annals of Mathematical Statistics, 26(1), 117–121. https://doi.org/10.1214/aoms/1177728599
• R Core Team. (2023). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/
• Tabachnick, B. G., & Fidell, L. S. (2019). Using multivariate statistics (7th ed.). Pearson.
• Todorov, V., & Filzmoser, P. (2009). An object-oriented framework for robust multivariate analysis. Journal of Statistical Software, 32(3), 1–47. https://doi.org/10.18637/jss.v032.i03
• Wilks, S. S. (1932). Certain generalizations in the analysis of variance. Biometrika, 24(3/4), 471–494. https://doi.org/10.2307/2331979