L’analyse statistique des données issues de protocoles expérimentaux en psychologie et en sciences comportementales se heurte fréquemment aux contraintes imposées par la distribution empirique des variables mesurées. Lorsque les chercheurs étudient l’évolution temporelle d’un construit psychologique, l’impact différentiel de plusieurs psychotropes ou l’efficacité relative d’interventions thérapeutiques successives auprès d’une cohorte identique, le plan expérimental à mesures répétées s’impose naturellement. Toutefois, l’application de l’analyse de variance paramétrique classique pour mesures répétées requiert le respect scrupuleux de postulats contraignants, notamment la normalité univariée et multivariée des distributions, l’absence de valeurs aberrantes influentes et la sphéricité de la matrice de variance-covariance. Dans la pratique clinique, ces hypothèses théoriques se révèlent souvent violées, particulièrement face à des échantillons de taille modeste ou à des échelles de mesure de nature ordinale.
Face à ces défaillances distributionnelles, l’approche non paramétrique offre des alternatives méthodologiques robustes. Le test de Friedman constitue précisément l’homologue non paramétrique fondamental de l’analyse de variance à mesures répétées à un facteur. En opérant une conversion des valeurs quantitatives ou ordinales en rangs attribués au sein de chaque individu, cette méthode élimine l’influence déstabilisatrice des asymétries de distribution et neutralise la variabilité interindividuelle initiale. L’environnement statistique open source R s’impose aujourd’hui comme l’écosystème de référence pour l’implémentation, l’évaluation critique et la restitution graphique de ces protocoles non paramétriques.
Ce guide exhaustif a pour vocation d’accompagner les chercheurs, psychologues cliniciens, neuropsychologues et analystes de données dans la maîtrise intégrale du test de Friedman sous R. De ses fondements mathématiques théoriques jusqu’aux préconisations de rédaction académique selon les normes de l’American Psychological Association, ce document détaille chaque étape opérationnelle. Nous aborderons la structuration des jeux de données, le diagnostic statistique préalable, la mise en œuvre du test omnibus, l’ajustement rigoureux des comparaisons post-hoc, l’estimation de la taille d’effet par le coefficient de concordance de Kendall, ainsi que la gestion experte des données manquantes et des ex-æquo.
- 1. Introduction théorique et fondements du test de Friedman en psychologie
- 2. Conditions d’application et hypothèses méthodologiques préalables
- 3. Préparation de l’environnement de travail et structuration des données sous R
- 4. Formulation des hypothèses statistiques et principe mathématique du test
- 5. Exécution du test standard avec la fonction native friedman.test()
- 6. Application pratique en psychologie : analyse du temps de réaction sous quatre traitements
- 7. Interprétation académique et décomposition des résultats statistiques
- 8. Analyses post-hoc et ajustement des comparaisons multiples
- 9. Calcul et interprétation de la taille d’effet : le W de Kendall
- 10. Visualisation graphique avancée des résultats sous R
- 11. Gestion des cas particuliers et écueils méthodologiques
- 12. Rédaction des résultats selon les normes de l’APA et bonnes pratiques
- Références
1. Introduction théorique et fondements du test de Friedman en psychologie
1.1 Origine historique et place dans la démarche non paramétrique
Le test de Friedman a été formulé pour la première fois en 1937 par l’économiste et statisticien américain Milton Friedman dans son article pionnier publié dans le Journal of the American Statistical Association. Conçu initialement pour évaluer des classements préférentiels et des séries temporelles économiques sans présumer de la normalité sous-jacente des variables, cet outil a très rapidement dépassé son champ disciplinaire originel pour s’ancrer au cœur des sciences du comportement, de la psychologie expérimentale et des neurosciences cognitives.
La nécessité méthodologique à l’origine de ce développement repose sur les limites intrinsèques des estimateurs gaussiens face à des phénomènes psychologiques asymétriques. En psychométrie et en clinique, les distributions tronquées, bimodales ou fortement étirées vers les valeurs extrêmes constituent la règle plutôt que l’exception. L’analyse de variance traditionnelle (ANOVA) à mesures répétées repose sur le postulat d’une distribution normale des résidus intra-sujets. Lorsque ce postulat s’effondre, l’estimation des carrés moyens et des rapports F s’avère biaisée, entraînant une distorsion incontrôlée de l’erreur de première espèce ou une perte dramatique de puissance statistique.
Le test de Friedman contourne cette difficulté conceptuelle en substituant aux observations brutes leurs rangs relatifs, calculés exclusivement au sein de chaque bloc expérimental horizontal, c’est-à-dire au sein de chaque sujet. Cette opération de mise en rang transforme l’espace métrique original en une métrique ordinale pure. Ainsi, le test s’affranchit totalement de la forme de la distribution d’origine tout en préservant l’ordonnancement intra-individuel des performances à travers les conditions expérimentales, neutralisant par là même l’hétérogénéité des moyennes individuelles.
1.2 Pertinence du test dans les protocoles expérimentaux en psychologie
Dans la recherche contemporaine en psychologie, le test de Friedman trouve une application privilégiée lors de l’évaluation longitudinale d’interventions thérapeutiques chez un groupe unique de patients. Considérons, à titre d’illustration, l’administration séquentielle de trois modules de thérapie cognitive et comportementale (gestion du stress, restructuration cognitive et exposition graduée) auprès de sujets présentant un trouble anxieux généralisé. Chaque individu agissant comme son propre contrôle, les caractéristiques idiosyncrasiques (telles que le niveau socio-économique, la personnalité de base ou la vulnérabilité génétique) sont mécaniquement contrôlées par le dispositif en blocs complets.
Un autre champ d’application critique concerne l’analyse de variables cognitives soumises à d’intenses effets de plafond ou de plancher. Lors de tâches neuropsychologiques évaluant les fonctions exécutives, la mémoire de travail ou l’attention sélective, il est fréquent qu’une proportion substantielle de participants atteigne le score maximal dans les conditions de faible charge cognitive, ou au contraire échoue systématiquement dans les conditions de surcharge. Ces saturations rendent les données non normales et non linéaires. L’utilisation des rangs intra-sujets permise par Friedman permet de capturer la dégradation ordinale des performances sans subir les distorsions d’échelle inhérentes aux scores bruts.
Enfin, les échelles psychométriques ordinales de type Likert, massivement déployées dans les questionnaires de personnalité, de bien-être ou d’évaluation de l’alliance thérapeutique, ne possèdent pas les propriétés métriques d’intervalles constants. Traiter la distance séparant « tout à fait d’accord » de « d’accord » comme strictement équivalente à celle séparant « neutre » de « pas d’accord » constitue une entorse méthodologique majeure. Le test de Friedman traite légitimement ces variables comme de véritables données ordinales, offrant une robustesse inférentielle inégalée face à la variabilité interindividuelle et aux artefacts de mesure.
2. Conditions d’application et hypothèses méthodologiques préalables
2.1 Structure du plan expérimental et dépendance des observations
La validité des conclusions issues du test de Friedman dépend de l’adéquation formelle du protocole d’échantillonnage à la structure théorique du modèle en blocs complets randomisés. Ce plan expérimental exige qu’un échantillon unique de n participants ou unités d’observation soit évalué de manière répétée à travers l’ensemble des k conditions expérimentales ou points temporels retenus, avec k supérieur ou égal à 3. Lorsque k est égal à 2, le test de Friedman devient mathématiquement équivalent au test des rangs signés de Wilcoxon pour échantillons appariés.
L’hypothèse maîtresse de ce modèle réside dans l’indépendance mutuelle stricte entre les différents blocs ou sujets de l’étude. Si les mesures réalisées au sein d’un même individu à travers les k conditions sont intrinsèquement appariées et dépendantes, le comportement ou la performance d’un sujet donné ne doit en aucun cas influencer ou contraindre les réponses observées chez un autre participant. Toute contamination interindividuelle, telle qu’une passation en groupe sans contrôle de l’influence sociale ou un effet de contagion relationnelle, vicie fondamentalement la structure de dépendance et invalide le calcul probabiliste.
Par ailleurs, le modèle standard de Friedman postule l’absence totale d’interaction statistique entre les blocs et les traitements. Cela signifie mathématiquement que la hiérarchie relative des effets de traitement est présumée constante ou homogène au sein de la population parente, sans qu’un sous-groupe latent d’individus ne réagisse de façon diamétralement inverse aux autres participants sous l’effet d’une variable modératrice non contrôlée.
2.2 Nature des variables et sensibilité de l’échelle de mesure
La variable dépendante analysée dans le cadre du test de Friedman doit impérativement posséder une échelle de mesure au moins ordinale. Elle peut être continue (exprimée par exemple en temps de réaction mesuré en millisecondes, en concentrations plasmatiques de cortisol salivaire ou en pourcentages d’erreurs) ou discrètement ordinale (scores agrégés d’inventaires d’anxiété, échelles de douleur visuelles analogiques ou réponses étalonnées sur des échelles de Likert). L’exigence essentielle réside dans la capacité de l’expérimentateur à établir un ordre strict ou partiel sans ambiguïté entre les réalisations de la variable.
La variable indépendante, quant à elle, représente le facteur intra-sujet ou traitement. Il s’agit obligatoirement d’une variable nominale ou ordinale à niveaux fixes, comportant un minimum de trois modalités distinctes. Ces modalités peuvent refléter une progression temporelle ordonnée (par exemple : ligne de base, post-intervention immédiate, suivi à six mois) ou des conditions qualitatives contrastées présentées dans un ordre contrebalancé (par exemple : condition placebo, psychotrope molécule A, psychotrope molécule B).
Il est en outre indispensable d’assurer l’homogénéité conceptuelle des mesures à travers les répétitions. L’instrument psychométrique ou le protocole de recueil doit mesurer scrupuleusement le même construit théorique sous chaque condition expérimentale. L’introduction subreptice d’un instrument alternatif ou d’une version révisée en cours d’expérimentation introduirait une erreur de mesure différentielle qui compromettrait l’interprétation des rangs intra-blocs.
2.3 Vérification des postulats statistiques préalables sous R
Avant d’opter formellement pour le test de Friedman, le chercheur doit soumettre ses données à une inspection diagnostique rigoureuse au sein de l’environnement R afin de déterminer si l’ANOVA paramétrique à mesures répétées doit être rejetée. La première étape consiste à examiner la normalité de la variable dépendante au sein de chaque condition expérimentale, ou préférentiellement la normalité des résidus du modèle linéaire mixte. L’évaluation s’effectue au moyen de diagrammes quantile-quantile (Q-Q plots) générés via les fonctions qqnorm() et qqline() du package de base, ou avec ggqqplot() du package spécialisé ggpubr, complétés par le test d’adéquation de Shapiro-Wilk appliqué par la fonction shapiro_test() du package rstatix.
La seconde étape diagnostique capitale concerne le postulat de sphéricité, également nommé condition de circularité de Huynh-Feldt et Box. La sphéricité stipule que les variances de l’ensemble des différences possibles entre chaque paire de conditions expérimentales doivent être rigoureusement égales. Cette condition est formellement éprouvée sous R à l’aide du test de sphéricité de Mauchly, implémenté par exemple via la fonction anova_test() du package rstatix.
Lorsque le test de Mauchly s’avère statistiquement significatif (p inférieur à 0.05), le postulat de sphéricité est formellement réfuté. Si cette violation s’accompagne d’une distribution asymétrique mise en évidence par les diagrammes Q-Q et d’un effectif d’échantillon modeste (n inférieur à 30), les corrections paramétriques usuelles de Greenhouse-Geisser ou de Huynh-Feldt peuvent présenter une puissance statistique médiocre ou des biais d’estimation. Dans cette configuration empirique précise, le basculement vers le test non paramétrique de Friedman s’impose comme la stratégie d’analyse la plus parcimonieuse et scientifiquement défendable.
3. Préparation de l’environnement de travail et structuration des données sous R
3.1 Installation et chargement des bibliothèques nécessaires
Pour mener à bien l’ensemble du pipeline analytique, depuis le reformatage des matrices brutes jusqu’aux calculs d’effets de taille et visualisations graphiques avancées, il convient de doter la session R des extensions méthodologiques appropriées. Le moteur statistique de base de R dispose nativement du package stats, qui intègre la fonction canonique de calcul du test. Cependant, l’enrichissement de l’espace de travail par des packages spécialisés garantit une fluidité syntaxique et une conformité aux standards contemporains de la science ouverte.
Le méta-package tidyverse constitue la pierre angulaire de la manipulation des structures de données et de l’ingénierie des variables grâce à ses composantes dplyr, tidyr et ggplot2. Pour l’inférence non paramétrique spécifique, la bibliothèque rstatix offre des fonctions conviviales compatibles avec les opérations de chaînage (pipe). Parallèlement, le package coin s’avère incontournable pour l’estimation de tests exacts et d’approches par permutation en présence de micro-échantillons. Enfin, la bibliothèque effectsize met à disposition des estimateurs standardisés de magnitude d’effet.
L’initialisation de l’environnement s’effectue classiquement par la commande d’installation install.packages(c("tidyverse", "rstatix", "coin", "effectsize", "PMCMRplus")), suivie de leur importation séquentielle dans la session active à l’aide de la primitive library(). Cette préparation méthodique prévient tout conflit d’espace de noms et arme le praticien pour l’intégralité des étapes computationnelles ultérieures.
3.2 Formatage des données : transition du format large au format long
L’une des sources majeures d’erreurs d’exécution dans R réside dans la confusion entre les formats structurels tabulaires. Historiquement, les fichiers de saisie clinique (issus de logiciels tels que SPSS ou de tableurs comme Excel) adoptent quasi systématiquement le format large, qualifié de wide format. Dans cette disposition horizontale, chaque ligne représente un participant unique, et les colonnes successives matérialisent les mesures répétées prélevées aux différentes conditions ou échéances temporelles.
Bien que le format large soit directement accepté par certaines signatures syntaxiques vectorielles de base, la grammaire unifiée du tidyverse et les fonctions analytiques modernes exigent impérativement une structure au format long, ou tidy format. Dans un tableau au format long, chaque ligne unitaire correspond à une observation singulière. Le jeu de données se compose alors d’au minimum trois colonnes fondamentales : une colonne d’identification unique du participant ou bloc (par exemple, patient_id), une colonne catégorielle codifiant la condition expérimentale ou le temps (par exemple, condition), et une colonne numérique consignant le score observé (par exemple, score).
La transition opératoire du format large vers le format long s’exécute avec une remarquable élégance sous R grâce à la fonction pivot_longer() issue du package tidyr. L’analyste y spécifie les colonnes sources à regrouper, le nom de la variable destination recevant les étiquettes de traitement via l’argument names_to, et le vecteur de stockage des valeurs numériques via values_to. Immédiatement après cette mutation, il est impératif d’assigner formellement aux variables d’identification et de traitement le type facteur (factor), condition sine qua non pour que les algorithmes de test reconnaissent adéquatement la partition en blocs.
4. Formulation des hypothèses statistiques et principe mathématique du test
4.1 Définition des hypothèses nulle et alternative
La conceptualisation probabiliste du test de Friedman obéit à une logique rigoureuse centrée sur l’homogénéité des positions ordinales. L’hypothèse nulle, formellement désignée sous le symbole H0, énonce qu’il n’existe aucune différence systématique entre les distributions sous-jacentes de la variable dépendante à travers les k conditions de traitement au sein de la population cible. Sous H0, pour chaque participant, l’ordonnancement des scores à travers les k conditions résulte uniquement des fluctuations aléatoires de l’échantillonnage ou du bruit de mesure. En d’autres termes, les k variables aléatoires continues appariées possèdent des fonctions de répartition marginales interchangeables, ce qui implique l’égalité des médianes relatives de rang.
À l’inverse, l’hypothèse alternative H1 stipule qu’au moins l’une des conditions expérimentales tend systématiquement à engendrer des valeurs plus élevées ou plus faibles que les autres, provoquant un décalage horizontal de sa distribution relative. Il est crucial de noter que le test de Friedman est un test omnibus directionnel de position : il ne recherche pas de simples divergences de dispersion ou de kurtosis, mais bien une asymétrie systématique dans la prédominance stochastique d’au moins un traitement vis-à-vis d’un autre.
Le rejet de l’hypothèse nulle ne permet en aucun cas d’isoler unilatéralement quelle modalité diffère de quelle autre. La découverte d’une significativité statistique globale à ce stade indique uniquement l’invraisemblance d’une répartition aléatoire uniforme des rangs au sein des blocs, imposant la mise en œuvre ultérieure de procédures exploratoires post-hoc.
4.2 Algorithme de calcul du rang et statistique du Chi-deux
Le mécanisme calculatoire interne du test de Friedman est remarquablement transparent. Soit un plan constitué de n sujets (blocs horizontaux) évalués sous k conditions de traitement (colonnes verticales). L’algorithme opère de manière purement intra-individuelle : pour chaque bloc i (allant de 1 à n), les k observations sont ordonnées par ordre croissant et reçoivent un rang r_ij variant de 1 (attribué à la valeur la plus faible du sujet) à k (attribué à sa valeur la plus élevée). Si des égalités surviennent au sein d’une ligne, on leur confère classiquement le rang moyen correspondant à leurs positions virtuelles.
Une fois l’intégralité de la matrice transformée en rangs intra-sujets, la somme des rangs attribués à chaque condition expérimentale j (allant de 1 à k) est calculée et notée R_j. Sous l’hypothèse nulle d’indépendance et d’équiprobabilité, la somme espérée des rangs pour n’importe quelle colonne correspond à la constante théorique n(k + 1) / 2. Toute déviation observée de R_j par rapport à cette valeur moyenne centrale témoigne d’un effet systématique du traitement.
La statistique de test standard de Friedman, couramment notée Q ou Chi-deux de Friedman, quantifie l’ampleur quadratique globale de ces écarts. Elle s’exprime par la formule analytique suivante :
Q = [12 / (n * k * (k + 1))] * Somme(R_j^2) – 3 * n * (k + 1)
Lorsque le produit n * k est suffisamment grand, la distribution d’échantillonnage de la statistique Q converge asymptotiquement vers une loi théorique du Chi-deux à k – 1 degrés de liberté. C’est sur cette comparaison entre la valeur empirique obtenue de Q et le quantile critique de la distribution théorique du Chi-deux que repose la dérivation de la p-valeur du test.
5. Exécution du test standard avec la fonction native friedman.test()
5.1 Syntaxe vectorielle directe : friedman.test(y, groups, blocks)
L’implémentation originelle du test au sein de l’environnement de base de R réside dans la fonction friedman.test() hébergée dans le module fondamental stats. Cette fonction offre deux interfaces syntaxiques distinctes. La première correspond à l’approche vectorielle explicite, dans laquelle l’utilisateur fournit les données décomposées en trois vecteurs indépendants de même longueur.
Dans cette signature, le premier paramètre nommé y désigne le vecteur numérique ou ordinal contenant l’ensemble des scores mesurés. Le deuxième paramètre, groups, prend la forme d’un vecteur ou facteur codifiant l’assignation de chaque observation à son groupe de traitement ou sa condition temporelle. Enfin, le troisième paramètre obligatoire, blocks, spécifie le vecteur ou facteur identifiant l’unité expérimentale ou le participant ayant fourni la réponse.
Bien que cette syntaxe vectorielle soit parfaitement rigoureuse sur le plan computationnel, elle présente l’inconvénient pratique d’exiger la manipulation permanente de vecteurs isolés, augmentant le risque d’erreurs d’alignement d’index lors des manipulations préliminaires en mémoire. Elle demeure néanmoins utile dans les environnements algorithmiques automatisés où les vecteurs sont générés dynamiquement sans transit par un conteneur tabulaire structuré.
5.2 Syntaxe par formule : la notation y ~ groups | blocks
La seconde interface offerte par friedman.test() repose sur le formalisme des formules symboliques de S et R, largement préféré par la communauté académique pour sa clarté conceptuelle et son élégance d’écriture. Cette syntaxe mobilise l’opérateur de conditionnement vertical conditionnel (la barre verticale |) pour spécifier la dépendance intra-bloc.
La formule s’articule sous l’expression canonique variable_reponse ~ facteur_traitement | identifiant_sujet. L’argument data permet d’adjoindre directement le tableau de données sous-jacent, qu’il s’agisse d’un data frame standard ou d’un tibble manipulé sous le tidyverse. Cette écriture isole impeccablement la structure du plan : à gauche du tilde figure la réponse continue ou ordinale ; à droite, le facteur expérimental principal, suivi de la spécification explicite de la variable de blocage agissant comme facteur de regroupement apparié.
L’utilisation de cette syntaxe par formule présente l’immense avantage préventif d’interdire implicitement toute altération de l’alignement des lignes et garantit la portabilité immédiate des scripts d’analyse au sein de rapports reproductibles sous R Markdown ou Quarto.
6. Application pratique en psychologie : analyse du temps de réaction sous quatre traitements
6.1 Contexte clinique et génération du jeu de données
Afin de concrétiser l’apprentissage par un cas empirique réaliste, imaginons une étude de psychopharmacologie cognitive menée auprès d’une cohorte clinique restreinte composée de dix patients souffrant de troubles déficitaires de l’attention avec hyperactivité (TDAH). L’équipe de recherche souhaite évaluer l’impact aigu sur la vitesse de traitement de l’information de quatre molécules psychotropes distinctes administrées selon un plan croisé contrebalancé : un placebo inerte (Placebo), une molécule dopaminergique standard (DopaStandard), une formulation novatrice à libération prolongée (ProlongEx) et une molécule sotoninergique modératrice (SotoMod).
La variable dépendante retenue est le temps moyen de réaction (exprimé en millisecondes) enregistré lors d’une tâche d’inhibition motrice de type Go/No-Go. Compte tenu de la petite taille de l’échantillon (n = 10) et de la présence notable d’asymétries dans les temps de latence neuropsychologiques, le recours à une modélisation non paramétrique s’impose impérativement.
Nous instancions ce jeu de données sous R en construisant un data frame structuré directement au format long. Pour garantir la reproductibilité absolue de nos calculs, nous définissons une graine pseudo-aléatoire via set.seed(123). Les dix sujets sont étiquetés séquentiellement de P01 à P10, et chacun traverse successivement les quatre conditions thérapeutiques, générant une matrice complète et parfaitement équilibrée de quarante observations individuelles exemptes de cellules orphelines.
6.2 Mise en œuvre du test pas à pas sur le jeu de données clinique
Une fois le tableau clinique généré dans l’environnement de travail sous le nom explicite de tdah_data, la première manœuvre analytique consiste à inspecter sa structure interne à l’aide de la commande str(tdah_data). Cette étape vérifie impérativement que la colonne patient est instanciée sous la forme d’un facteur à 10 niveaux, que la colonne traitement est un facteur à 4 niveaux et que la colonne temps_reaction correspond à un vecteur numérique de nombres réels.
L’exécution formelle du test de Friedman est alors déclenchée par l’instruction :
friedman_resultat <- friedman.test(temps_reaction ~ traitement | patient, data = tdah_data)
L’assignation explicite du résultat dans un objet dédié nommé friedman_resultat permet de capturer la structure de sortie sous forme de liste statistique enrichie de la classe htest. La simple saisie du nom de cet objet dans la console interactive de R provoque l’impression formatée du rapport d’analyse officiel, détaillant le nom du test exécuté, l’intitulé des données sources, la statistique numérique calculée, le nombre de degrés de liberté associés et la p-valeur exacte dérivée de l’asymptote.
7. Interprétation académique et décomposition des résultats statistiques
7.1 Compréhension de la sortie console de R
L’examen minutieux de la sortie textuelle générée par l’instruction précédente livre trois informations numériques fondamentales qui président à la prise de décision statistique. Le premier indicateur correspond au Friedman chi-squared, qui n’est autre que la réalisation empirique de la statistique Q détaillée dans notre section théorique. Une valeur élevée témoigne d’un écart substantiel entre les sommes de rangs observées sous chaque traitement et la moyenne théorique attendue sous l’hypothèse d’équiprobabilité.
Le second paramètre affiché est représenté par les degrés de liberté, abrégés df. Dans le cadre canonique du test de Friedman, cette quantité est strictement bornée par l’égalité df = k – 1, où k représente le nombre de conditions de traitement étudiées. Dans notre application neuropsychologique comprenant quatre psychotropes, le système statistique possède précisément 4 – 1 = 3 degrés de liberté.
Le troisième paramètre, décisif pour l’inférence, est la p-value. Cette valeur probabiliste quantifie le risque d’observer une statistique de test au moins aussi extrême que celle mesurée empiriquement sous la supposition formelle que l’hypothèse nulle soit strictement vraie. Si cette p-valeur s’avère inférieure au seuil conventionnel alpha fixé à 0.05 (ou a fortiori à 0.01 ou 0.001), le chercheur rejette formellement l’hypothèse nulle H0 et conclut à l’existence d’un effet statistiquement significatif du facteur traitement sur la vitesse de réaction cognitive des participants.
7.2 Extraction programmatique des éléments du test
Dans un contexte moderne de rédaction automatisée de manuscrits scientifiques, d’intégration continue sous Quarto ou de développement d’applications analytiques dynamiques avec Shiny, l’analyste ne saurait se satisfaire d’une simple lecture visuelle de la console. Il est essentiel d’extraire programmatiquement les composantes atomiques de l’objet de test.
L’objet de classe htest produit par friedman.test() est une liste ordonnée dont les éléments internes sont directement adressables à l’aide de l’opérateur dollar $. Ainsi, l’expression friedman_resultat$statistic extrait vectoriellement la valeur scalaire de la statistique du Chi-deux. De façon analogue, friedman_resultat$parameter restitue les degrés de liberté, tandis que friedman_resultat$p.value fournit la p-valeur brute sous forme d’un flottant double précision.
Cette modularité logicielle permet d’injecter ces scalaires directement au sein de fonctions de formatage textuel automatisées telles que sprintf() ou le package glue, éradiquant tout risque d’erreur humaine de transcription manuelle lors de la finalisation des sections de résultats empiriques.
8. Analyses post-hoc et ajustement des comparaisons multiples
8.1 Nécessité des tests a posteriori après un test global significatif
L’obtention d’un résultat omnibus statistiquement significatif lors de l’exécution du test de Friedman constitue une étape d’étape indispensable mais fondamentalement incomplète. En rejetant l’hypothèse nulle d’homogénéité globale des distributions, le chercheur sait qu’au moins une condition expérimentale diverge de l’une des autres, sans disposer de la moindre information sur l’architecture précise de ces contrastes. Est-ce la molécule ProlongEx qui accélère significativement les réponses par rapport au Placebo ? La molécule DopaStandard se distingue-t-elle de SotoMod ?
Toutefois, la conduite désordonnée de comparaisons par paires expose l’expérimentateur à un péril méthodologique redoutable : l’inflation exponentielle du risque d’erreur de première espèce (l’erreur de type I). Pour k conditions expérimentales, le nombre total de comparaisons bilatérales possibles s’élève à C = k * (k – 1) / 2. Dans notre protocole clinique à quatre traitements, six comparaisons par paires doivent être menées de front.
Si chaque test individuel est exécuté au seuil de signification standard alpha de 0.05 sans correction protectrice, le risque global cumulé d’observer au moins un faux positif (le Family-Wise Error Rate ou FWER) grimpe drastiquement selon la loi 1 – (1 – 0.05)^6, atteignant environ 26.5 %. La mise en œuvre de procédures de contrôle et d’ajustement statistique a posteriori constitue donc un impératif de rigueur déontologique et scientifique incontournable.
8.2 Test des rangs signés de Wilcoxon apparié avec correction
L’approche post-hoc la plus intuitive et couramment adoptée dans la littérature psychologique consiste à appliquer de manière séquentielle le test non paramétrique des rangs signés de Wilcoxon pour mesures appariées sur chacune des paires de conditions envisageables, en y adjoignant un algorithme d’ajustement de la p-valeur. Sous R, cette opération s’exécute avec fluidité via la commande pairwise.wilcox.test() en fixant obligatoirement l’argument booléen paired = TRUE.
Plusieurs stratégies de contrôle de l’erreur globale peuvent être spécifiées via l’argument p.adjust.method. La méthode de Bonferroni représente l’approche la plus classique et conservatrice : elle consiste à multiplier chaque p-valeur individuelle observée par le nombre total de contrastes C, ou à abaisser le seuil de décision critique à alpha / C. Bien qu’elle garantisse un contrôle hermétique du FWER, son extrême sévérité tend à laminer la puissance statistique, augmentant considérablement le risque d’erreur de seconde espèce (erreur de type II).
Pour parer à ce conservatisme excessif, la procédure séquentielle descendante de Holm (souvent appelée méthode de Holm-Bonferroni) est universellement recommandée par les méthodologistes. Elle ordonne hiérarchiquement les p-valeurs obtenues de la plus petite à la plus grande et ajuste le critère de rejet de façon dynamique étape par étape, offrant un contrôle tout aussi infaillible du risque de faux positif tout en préservant substantiellement la sensibilité de détection. Enfin, dans les contextes exploratoires à grand nombre de modalités, l’ajustement de Benjamini-Hochberg (contrôle du taux de fausse découverte ou False Discovery Rate, FDR) peut constituer un arbitrage judicieux.
8.3 Approches post-hoc spécifiques : tests de Conover et de Nemenyi
En dépit de sa popularité, l’enchaînement de tests de Wilcoxon appariés isolés souffre d’un défaut théorique : chaque test de Wilcoxon réassigne de nouveaux rangs sur la base exclusive de la paire de conditions examinée, ignorant la structure globale de rangs préalablement établie sur l’ensemble de l’expérience par le test omnibus de Friedman. Pour pallier cette rupture d’invariance, des procédures post-hoc dévolues spécifiquement aux plans en blocs ont été développées, au premier rang desquelles figurent les tests de Nemenyi et de Conover.
Le test post-hoc de Nemenyi, accessible sous R via la fonction frdManyOneDurbinTest() ou kwManyOneConoverTest() du package spécialisé PMCMRplus, évalue les différences absolues observées entre les sommes globales de rangs des traitements en s’appuyant sur la distribution de l’étendue studentisée. Très conservateur, il s’avère particulièrement indiqué lorsque l’analyste compare toutes les conditions deux à deux dans un plan rigoureusement équilibré.
À l’opposé, le test post-hoc de Conover (implémenté via la fonction frdAllPairsConoverTest() du package PMCMRplus) réutilise la variance résiduelle de la matrice globale des rangs de Friedman pour formuler une statistique distribuée selon une loi de Student. Il en résulte un gain notable de puissance statistique, autorisant la mise en lumière de contrastes fins que la rudesse de Bonferroni ou de Nemenyi aurait laissés dans l’ombre. L’analyste rigoureux prendra soin d’expliciter le choix de son test post-hoc dans son rapport de recherche en justifiant l’équilibre recherché entre puissance et conservatisme.
9. Calcul et interprétation de la taille d’effet : le W de Kendall
9.1 Importance de quantifier la magnitude de l’effet en psychologie
Dans la recherche contemporaine en psychologie et en sciences comportementales, la seule déclaration de la significativité statistique via la p-valeur est jugée notoirement insuffisante. La p-valeur informe uniquement sur la plausibilité de l’hypothèse nulle au regard des données recueillies, mais elle demeure viscéralement dépendante de la taille de l’échantillon. Un échantillon de très grande taille peut révéler une différence statistiquement significative pour un effet dont l’amplitude pratique ou clinique est rigoureusement négligeable. Inversement, une étude portant sur un effectif clinique réduit mais investiguant une intervention d’une efficacité spectaculaire peut échouer à atteindre le seuil de 0.05 par simple manque de puissance.
Face à cette aporie, les comités de lecture internationaux et le manuel éditorial de l’APA imposent la communication systématique d’une mesure standardisée de la magnitude de l’effet. L’estimation de la taille d’effet offre une mesure métrique universelle de l’impact substantiel du traitement, autorisant la comparaison directe des résultats à travers différents protocoles et fournissant le substrat indispensable aux futures méta-analyses quantitatives.
9.2 Calcul du coefficient de concordance W de Kendall sous R
Dans l’architecture non paramétrique des plans à mesures répétées, l’estimateur de taille d’effet standard dévolu au test de Friedman est le coefficient de concordance W de Kendall. Cet indice statistique capture l’harmonie ou le degré d’accord manifesté par l’ensemble des n blocs (les participants) quant au classement ordonné des k traitements expérimentaux.
Le coefficient W de Kendall se déduit de manière limpide de la statistique du Chi-deux de Friedman Q selon la formulation algébrique exacte suivante :
W = Q / (n * (k – 1))
L’indice W est ainsi naturellement borné au sein de l’intervalle continu allant de 0 à 1. Une valeur de 0 indique une absence absolue de concordance entre les participants (les classements individuels se répartissent de manière chaotique et aléatoire, chaque traitement ayant une probabilité identique d’être classé premier ou dernier). À l’extrême opposé, une valeur de 1 matérialise une unanimité parfaite et invariable (l’ensemble unanime des participants classe les traitements dans un ordre strictement identique).
Sous l’environnement R, le calcul de cet indice s’exécute de manière élégante et automatisée grâce à la fonction friedman_effsize() mise à disposition par le package rstatix. Cette fonction prend en charge les formules de même structure que le test omnibus et peut en outre générer un intervalle de confiance robuste à 95 % par rééchantillonnage non paramétrique (procédure de bootstrapping) en spécifiant l’argument ci = TRUE.
9.3 Grille d’interprétation clinique et psychométrique
La qualification qualitative de la magnitude du coefficient de concordance W de Kendall repose sur des repères conventionnels hérités des travaux de Jacob Cohen et adaptés par les statisticiens non paramétriques. Ces seuils de référence guident le jugement diagnostique du psychologue dans l’appréciation de son effet expérimental :
- Effet négligeable : W strictement inférieur à 0.10. Les classements intra-sujets divergent profondément et l’effet thérapeutique s’avère indétectable en pratique clinique courante.
- Effet de faible amplitude : W compris entre 0.10 et 0.29. Une tendance ordinale commence à émerger au sein de la cohorte, mais la variabilité des réponses individuelles reste prépondérante.
- Effet modéré : W compris entre 0.30 et 0.49. L’alignement des classements individuels est substantiel ; une proportion significative de la trajectoire clinique est directement gouvernée par la nature des traitements.
- Effet de forte amplitude : W supérieur ou égal à 0.50. Il existe une concordance majeure entre les participants quant à la hiérarchie des conditions expérimentales, attestant d’une efficacité ou d’une toxicité différentielle massive.
L’interprétation de ces seuils doit néanmoins toujours être modulée en fonction du contexte disciplinaire. Dans les recherches exploratoires sur des construits de personnalité stables ou des troubles psychiatriques complexes et réfractaires, l’observation d’un W modéré de 0.35 peut représenter une avancée clinique remarquable, tandis que lors de tests psychophysiques élémentaires hautement contrôlés en laboratoire, un W de 0.60 constituera un prérequis standard de validation interne.
10. Visualisation graphique avancée des résultats sous R
10.1 Boîtes à moustaches enrichies avec trajectoires individuelles
La communication scientifique des résultats d’un plan à mesures répétées analysé par le test de Friedman ne saurait se limiter à de simples tableaux d’effectifs ou à des graphiques à barres agrégés. Ces derniers occultent la spécificité nodale du protocole : le couplage intra-individuel des observations. Pour rendre compte de la distribution ordinale tout en révélant la dynamique longitudinale propre à chaque participant, la boîte à moustaches (boxplot) enrichie des trajectoires individuelles constitue le standard de visualisation d’excellence sous la bibliothèque ggplot2.
Dans ce dispositif graphique élaboré, la boîte à moustaches centrale matérialise pour chaque traitement les trois quartiles cardinaux (médiane, 25e et 75e centiles) ainsi que les limites de dispersion de Tukey. Par-dessus cette armature descriptive, les points de données brutes individuels sont superposés à l’aide de la couche geom_point(). L’élément déterminant réside dans l’adjonction de segments continus reliant les scores d’un même patient à travers les quatre conditions au moyen de l’instruction geom_line(aes(group = patient), alpha = 0.4, linetype = "dashed").
Cette visualisation permet au lecteur de constater immédiatement si le basculement d’une condition à l’autre induit une élévation ou un abaissement unanime de la trajectoire chez l’ensemble des sujets, ou si au contraire certains profils cliniques atypiques présentent une résistance paradoxale au traitement, éclairant visuellement la valeur obtenue pour le W de Kendall.
10.2 Graphiques de rangs moyens et profils de variations
Puisque le mécanisme décisionnel du test de Friedman opère exclusivement sur les rangs relatifs et non sur les valeurs métriques directes, il est méthodologiquement pertinent de construire un graphique représentant les rangs moyens assignés à chaque condition expérimentale. Ce type de figure offre une fidélité géométrique parfaite au fonctionnement mathématique interne du modèle non paramétrique.
Pour concrétiser cette restitution, les données sont préalablement transformées par bloc afin d’extraire le rang intra-sujet de chaque observation via la fonction rank() au sein d’une instruction groupée de dplyr. Les rangs moyens et leurs intervalles de dispersion associés (tels que l’écart-type des rangs ou l’écart interquartile) sont ensuite agrégés sous forme de coordonnées discrètes. L’utilisation conjointe des fonctions geom_point(size = 4) et geom_errorbar(width = 0.2) permet de matérialiser élégamment la position relative de chaque modalité.
Enfin, l’incorporation d’annotations statistiques directes sur le canevas graphique, synthétisant la statistique globale de Friedman (valeur de Q, degrés de liberté et p-valeur) ainsi que les accolades signalant les paires statistiquement discriminées par les analyses post-hoc (réalisables via le package complémentaire ggpubr ou ggsignif), confère à la figure une autonomie documentaire conforme aux exigences des publications scientifiques de premier plan.
11. Gestion des cas particuliers et écueils méthodologiques
11.1 Traitement des ex-æquo dans l’attribution des rangs
La dérivation mathématique fondamentale du test de Friedman postule la continuité absolue de la variable dépendante sous-jacente, ce qui devrait théoriquement exclure l’émergence de valeurs parfaitement identiques au sein d’un même bloc. Dans la réalité empirique des laboratoires de psychologie, la présence de scores ex-æquo (communément désignés sous le vocable anglais de ties) est extrêmement fréquente, que ce soit en raison de la discrétisation imposée par les échelles de Likert ou de l’imprécision inhérente aux appareils de mesure chronométrique.
Lorsque deux ou plusieurs observations d’un même participant présentent des valeurs numériques rigoureusement identiques, la procédure conventionnelle leur attribue le rang moyen des rangs qu’elles auraient occupés si elles avaient été distinctes. Par exemple, si les deux valeurs les plus basses d’un individu sont à égalité, elles occupent virtuellement les positions 1 et 2, recevant chacune le rang moyen de (1 + 2) / 2 = 1.5.
L’existence d’ex-æquo a pour conséquence mécanique de réduire la variance totale de la distribution d’échantillonnage des rangs. Si la proportion d’égalités est substantielle et qu’aucune correction n’est appliquée, la statistique de test standard Q sous-estime l’écart réel et conduit à une p-valeur artificiellement conservatrice. Heureusement, la fonction friedman.test() de R intègre automatiquement un coefficient d’ajustement qui divise l’expression de la statistique par un facteur réducteur dépendant de la somme des cubes des multiplicités d’égalités au sein de chaque bloc. Ce mécanisme automatique garantit la robustesse du seuil décisionnel tant que la proportion d’égalités ne sature pas la quasi-totalité de l’espace d’observation.
11.2 Gestion des données manquantes et plans incomplets
Une limitation fonctionnelle majeure de la fonction native friedman.test() réside dans son intransigeance face aux valeurs manquantes. L’algorithme exige impérativement un plan en blocs complets et équilibrés : si un seul participant omet de renseigner une condition sur les quatre proposées, la moindre valeur NA provoque le rejet complet du bloc ou l’avortement de l’exécution.
Face à des déperditions d’échantillon (phénomène de drop-out fréquent lors de suivis thérapeutiques longitudinaux), l’analyste ne doit pas recourir à l’amputation aveugle de tous les patients incomplets (analyse de cas complets), ce qui amputerait drastiquement la puissance statistique et introduirait de redoutables biais d’attrition. La solution méthodologique non paramétrique la plus élégante consiste à recourir au test de Skillings-Mack.
Le test de Skillings-Mack représente une généralisation mathématique du test de Friedman spécifiquement conçue pour accommoder les plans en blocs incomplets randomisés comportant des données manquantes complètement aléatoires (MCAR) ou aléatoires (MAR). Cette méthode applique une pondération différentielle aux blocs selon le nombre de mesures effectivement validées par chaque sujet. Sous R, cette procédure robuste s’exécute aisément via le package spécialisé Skillings.Mack à l’aide de la commande éponyme Ski.Mack(), fournissant une statistique de type Chi-deux parfaitement interprétable sans recourir à des procédures d’imputation artificielle hasardeuses.
11.3 Approches exactes versus asymptotiques pour très petits échantillons
Comme explicité dans la section mathématique, la conversion de la statistique Q en une p-valeur décisionnelle repose sur l’hypothèse d’une convergence asymptotique vers une loi théorique du Chi-deux. Cette approximation gaussienne n’est formellement valide que lorsque la taille de l’échantillon n et le nombre de modalités k sont suffisamment imposants. Lorsque les études portent sur des micro-échantillons cliniques (par exemple, n inférieur à 8 avec k égal à 3 ou 4), la distribution empirique exacte s’écarte notablement de la courbe théorique du Chi-deux, augmentant considérablement le risque d’erreurs inférentielles.
Dans ce contexte de rareté d’observations, l’analyste doit abandonner l’approximation asymptotique au profit d’un test exact fondé sur le calcul combinatoire exhaustif de toutes les permutations conditionnelles possibles des rangs au sein des blocs. L’environnement R permet cette transition computationnelle de haut niveau grâce à la bibliothèque d’inférence conditionnelle coin (abréviation de Conditional Inference Procedures in R).
Au travers de la fonction friedman_test() issue de ce package, l’utilisateur a la faculté de spécifier le paramètre distribution = "exact" (pour un dénombrement combinatoire complet) ou distribution = approximate(nresample = 100000) (pour une approximation par permutations de Monte Carlo à très haute résolution). Cette démarche garantit une rigueur méthodologique irréprochable et protège les conclusions de la publication contre toute remise en cause relative à la représentativité asymptotique de l’échantillon.
12. Rédaction des résultats selon les normes de l’APA et bonnes pratiques
12.1 Directives de rédaction selon la 7e édition de l’APA
La rédaction de la section des résultats statistiques au sein d’un manuscrit soumis aux revues de psychologie régies par l’American Psychological Association (APA, 7e édition) obéit à un protocole standardisé strict. La transparence et la complétude méthodologique exigent de ne jamais dissocier la significativité statistique de sa magnitude descriptive et de sa portée d’effet.
Pour un test de Friedman, le formalisme prescrit impose de rapporter explicitement :
- L’indicateur de tendance centrale approprié pour chaque modalité expérimentale (compte tenu de la nature ordinale ou asymétrique des données, il s’agit obligatoirement de la médiane accompagnée de l’écart interquartile, ou alternativement de la moyenne des rangs) ;
- La lettre grecque stylisée ou la désignation textuelle de la statistique de test, notée Chi-deux de Friedman ou classiquement Chi-deux avec en indice ses degrés de liberté entre parenthèses ;
- La taille de l’échantillon de blocs N inclus dans l’analyse finale ;
- La valeur numérique exacte de la statistique de test arrondie à deux décimales ;
- La p-valeur exacte arrondie à trois décimales (sauf lorsqu’elle est inférieure à 0.001, auquel cas elle est formulée sous la mention p < .001, sans zéro préliminaire conformément aux règles anglo-saxonnes si le texte est en anglais, ou avec zéro selon l’adaptation francophone) ;
- L’indice de taille d’effet standardisé (le W de Kendall), accompagné le cas échéant de son intervalle de confiance à 95 % ;
- La description explicite de la méthode d’ajustement employée pour neutraliser l’inflation de l’erreur de type I lors des comparaisons multiples post-hoc (par exemple, la correction de Holm ou de Bonferroni).
12.2 Exemple de paragraphe rédigé pour une publication scientifique
Pour matérialiser ces prescriptions normatives en situation réelle, voici un exemple clé en main de paragraphe rédigé, directement transposable dans la section Résultats d’un article académique ou d’une thèse de doctorat en psychologie :
« Afin d’évaluer l’impact différentiel des quatre conditions pharmacologiques sur la vitesse de traitement de l’information, un test de Friedman pour mesures répétées a été conduit sur les temps de réaction des participants. L’analyse révèle une variation statistiquement significative des performances cognitives entre les différents traitements, Chi-deux(3, N = 10) = 18.84, p < .001, attestant d’une magnitude d’effet très élevée avec un coefficient de concordance de Kendall W = 0.63, IC 95 % [0.38, 0.86]. Les temps médians de réaction s’élevaient respectivement à 485.5 ms (EIQ = 42.0) pour le Placebo, 412.0 ms (EIQ = 35.5) pour la molécule SotoMod, 360.5 ms (EIQ = 28.0) pour DopaStandard, et 315.0 ms (EIQ = 24.5) pour la formulation ProlongEx.
Les comparaisons post-hoc par paires, exécutées au moyen de tests des rangs signés de Wilcoxon pour échantillons appariés et assorties de la correction séquentielle de Holm pour comparaisons multiples, démontrent que le traitement ProlongEx induit une réduction statistiquement significative du temps de réaction par rapport au Placebo (p_ajusté = .012) ainsi que par rapport à SotoMod (p_ajusté = .036). De surcroît, la molécule DopaStandard se démarque significativement du Placebo (p_ajusté = .024). En revanche, la différence observée entre DopaStandard et ProlongEx n’atteint pas le seuil de significativité après ajustement conservateur (p_ajusté = .148). L’ensemble de ces profils est illustré à la Figure 1. »
12.3 Garantie de reproductibilité avec R Markdown et Quarto
La crise contemporaine de la reproductibilité scientifique a conduit la communauté des chercheurs en sciences humaines à réviser en profondeur leurs chaînes de traitement statistique. La juxtaposition manuelle d’un logiciel de statistiques, d’un tableur intermédiaire et d’un traitement de texte bureautique constitue le vecteur principal d’erreurs matérielles, d’incohérences de signalement et de perte de traçabilité des analyses.
L’adoption des systèmes de publication computationnelle reproductible tels que Quarto et R Markdown résout définitivement cette fragilité. En amalgamant au sein d’un document source unique le code exécutable R, la prose explicative et les appels directs aux variables dynamiques en ligne (inline R code), l’expérimentateur s’assure que chaque chiffre figurant dans son manuscrit découle directement et sans médiation de la base de données brute.
Pour sanctifier cette rigueur méthodologique, plusieurs bonnes pratiques s’imposent :
- Fixation explicite des graines aléatoires : Toute procédure computationnelle faisant intervenir un rééchantillonnage stochastique (notamment le bootstrap pour les intervalles de confiance ou les tests de permutations de Monte Carlo) doit impérativement être précédée de la commande
set.seed()afin de permettre la réplication rigoureusement identique des résultats à la décimale près par les réviseurs ; - Documentation de l’environnement matériel et logiciel : La clôture systématique de tout script d’analyse par l’instruction
sessionInfo()permet d’archiver la version exacte du compilateur R, du système d’exploitation sous-jacent et des versions incrémentales de chaque package invoqué ; - Formatage tabulaire normalisé : L’utilisation d’outils de mise en forme automatisée tels que les packages knitr::kable, flextable ou gtsummary permet d’exporter des tableaux d’analyses post-hoc respectant scrupuleusement les exigences typographiques de l’APA sans aucune retouche manuelle ultérieure, scellant ainsi l’intégrité de la démarche scientifique.
Références
- American Psychological Association. (2020). Publication manual of the American Psychological Association (7th ed.). American Psychological Association. https://doi.org/10.1037/0000165-000
- 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
- Conover, W. J. (1999). Practical nonparametric statistics (3rd ed.). John Wiley & Sons.
- Friedman, M. (1937). The use of ranks to avoid the assumption of normality implicit in the analysis of variance. Journal of the American Statistical Association, 32(200), 675–701. https://doi.org/10.1080/01621459.1937.10503522
- Friedman, M. (1940). A comparison of alternative tests of significance for the problem of m rankings. The Annals of Mathematical Statistics, 11(1), 86–92. https://doi.org/10.1214/aoms/1177731944
- Holm, S. (1979). A simple sequentially rejective multiple test procedure. Scandinavian Journal of Statistics, 6(2), 65–70. https://www.jstor.org/stable/4615733
- Hothorn, T., Hornik, K., van de Wiel, M. A., & Zeileis, A. (2006). A Lego system for conditional inference. The American Statistician, 60(3), 257–263. https://doi.org/10.1198/000313006X118430
- Kendall, M. G., & Babington Smith, B. (1939). The problem of m rankings. The Annals of Mathematical Statistics, 10(3), 275–287. https://doi.org/10.1214/aoms/1177732186
- Kassambara, A. (2023). rstatix: Pipe-friendly framework for basic statistical tests (R package version 0.7.2). Comprehensive R Archive Network. https://cran.r-project.org/package=rstatix
- Pohlert, T. (2014). The pairwise multiple comparison of mean ranks package (PMCMR). R package. https://cran.r-project.org/package=PMCMRplus
- R Core Team. (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/
- Skillings, J. H., & Mack, G. A. (1981). On the use of a Friedman-type statistic in balanced and unbalanced block designs. Technometrics, 23(2), 171–177. https://doi.org/10.1080/00401706.1981.10486261
- Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L. D., François, R., Grolemund, G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Pedersen, T. L., Miller, E., Bache, S. M., Müller, K., Ooms, J., Robinson, D., Seidel, D. P., Spinu, V., … Yutani, H. (2019). Welcome to the Tidyverse. Journal of Open Source Software, 4(43), 1686. https://doi.org/10.21105/joss.01686