Comment réaliser le test des rangs signés de Wilcoxon dans R

L’évaluation rigoureuse des changements intra-individuels constitue l’un des piliers méthodologiques fondamentaux des sciences quantitatives contemporaines, qu’il s’agisse de la biomédecine, de l’épidémiologie, de l’économie expérimentale ou de la psychométrie. Dans les protocoles longitudinaux ou quasi-expérimentaux à mesures répétées, la comparaison d’un même groupe d’individus avant et après une intervention clinique ou pharmacologique impose l’utilisation d’outils statistiques capables d’isoler l’effet du traitement de la variabilité interindividuelle inhérente aux sujets. Traditionnellement, le réflexe analytique consiste à mobiliser le test t de Student pour échantillons appariés. Néanmoins, ce dernier repose sur des postulats paramétriques exigeants, en particulier l’hypothèse selon laquelle la distribution des différences individuelles suit rigoureusement une loi normale au sein de la population parente.

Lorsque les données collectées présentent des asymétries prononcées, des queues de distribution épaisses, une contamination par des valeurs aberrantes ou lorsqu’elles dérivent d’échelles ordinales (telles que les échelles d’évaluation de Likert fréquemment déployées en sciences sociales), les hypothèses du modèle linéaire gaussien s’effondrent. L’application d’un test paramétrique dans ces conditions fausse le taux d’erreur de type I et dégrade substantiellement la puissance statistique, conduisant le chercheur soit à des conclusions faussement positives, soit à l’incapacité de détecter un effet clinique pourtant réel. C’est dans ce contexte méthodologique critique que s’impose le test des rangs signés de Wilcoxon, procédure d’inférence non paramétrique d’une puissance et d’une résilience remarquables.

Le présent guide propose une monographie exhaustive dédiée à l’apprentissage, à l’exécution et à l’interprétation avancée du test des rangs signés de Wilcoxon au sein de l’environnement de programmation statistique R. À travers une démarche alliant rigueur théorique, formulations mathématiques explicites, déconstructions algorithmiques et recommandations de publication conformes aux standards académiques internationaux les plus stricts, cet ouvrage accompagnera le praticien et le chercheur depuis les fondements conceptuels de l’analyse non paramétrique jusqu’à la modélisation avancée et la visualisation graphique prête pour publication.

1. Fondements théoriques et mathématiques du test des rangs signés de Wilcoxon

1.1 Origine et nature non paramétrique du test

Le test des rangs signés a été introduit en 1945 par le chimiste et statisticien américain Frank Wilcoxon dans un article fondateur intitulé Individual Comparisons by Ranking Methods, publié dans le Biometrics Bulletin. À cette époque, la pratique de l’inférence statistique était presque exclusivement dominée par les travaux de Ronald Fisher, Karl Pearson et William Sealy Gosset (Student), qui présupposaient tous que les variables observées découlaient d’une distribution gaussienne continue théorique. Wilcoxon a cherché à concevoir une procédure mathématique simple et rapide, capable de s’affranchir de cette contrainte distributionnelle tout en conservant une puissance statistique comparable à celle du test t de Student.

Le principe fondamental du test des rangs signés de Wilcoxon repose sur la substitution des grandeurs numériques métriques brutes par des rangs ordinaux attribués aux différences absolues observées entre des paires d’échantillons dépendants. Contrairement aux tests purement qualitatifs ou au test des signes élémentaire qui ne retient que l’orientation directionnelle (positive ou négative) du changement, la procédure de Wilcoxon intègre simultanément le signe de la différence et sa magnitude relative par rapport à l’ensemble des autres écarts observés dans l’échantillon. Ainsi, une différence de grande amplitude recevra un rang élevé, tandis qu’une variation minime se verra assigner un rang modeste, garantissant que les changements substantiels pèsent davantage dans la décision statistique finale.

Cette approche induit une distinction conceptuelle majeure entre la comparaison classique des moyennes arithmétiques et la comparaison des distributions médianes. Alors que la moyenne paramétrique est hyper-sensible aux valeurs marginales, la statistique de rang est invariante face à toute transformation monotone croissante des écarts. Elle ne requiert aucun postulat métrique strict sur l’intervalle séparant deux unités de mesure consécutives, conférant au test des rangs signés le statut de test non paramétrique par excellence pour les données continues asymétriques ou les variables d’intervalles approximatives issues des sciences du comportement.

1.2 Formulation mathématique de la statistique de test

Considérons un échantillon composé de n paires d’observations appariées notées (x_i,post, x_i,pre) pour chaque individu i variant de 1 à n. La première étape analytique consiste à calculer les écarts individuels bruts, définis formellement par l’équation :

d_i = x_i,post – x_i,pre

Dans la formulation originale de Wilcoxon, les observations pour lesquelles la différence est strictement nulle (d_i = 0), correspondant à une absence absolue de changement entre les deux conditions de mesure, sont systématiquement exclues de l’analyse inférentielle. La taille effective de l’échantillon est alors réduite à n’, où n’ représente le nombre de paires non nulles. Les valeurs absolues de ces écarts résiduels, notées |d_i|, sont ensuite ordonnées par ordre croissant, de la plus petite à la plus grande magnitude :

|d|₍₁₎ ≤ |d|₍₂₎ ≤ … ≤ |d|_(n’)

Un rang numérique entier allant de 1 à n’ est assigné à chacune de ces valeurs absolues. Si deux ou plusieurs valeurs absolues sont strictement identiques, on leur attribue un rang moyen correspondant à la moyenne arithmétique des rangs qu’elles auraient occupés en l’absence d’ex æquo. À l’issue de cette phase de classement, chaque rang R_i se voit réassigner le signe algébrique de l’écart initial dont il est issu, défini par la fonction indicatrice du signe de d_i.

La statistique de test peut alors être formulée à partir de deux composantes distinctes : la somme des rangs associés aux différences strictement positives (W⁺) et la somme des rangs associés aux différences strictement négatives (W^–). Mathématiquement, ces sommes s’expriment selon les formulations suivantes :

W⁺ = ∑_{i : di > 0} R_i

W^– = ∑_{i : di < 0} R_i

Puisque la somme totale des rangs de 1 à n’ est une quantité arithmétique invariante égale à n'(n’ + 1) / 2, la connaissance de W⁺ détermine de façon univoque W^– par la relation W^– = [n'(n’ + 1) / 2] – W⁺. La statistique de Wilcoxon, traditionnellement notée T ou W dans la littérature statistique, correspond selon les conventions soit au minimum de ces deux sommes pour un test bilatéral (min(W⁺, W^–)), soit directement à la somme des rangs positifs W⁺, convention adoptée par défaut dans l’implémentation algorithmique de la fonction de base du logiciel R.

1.3 Distribution exacte versus approximation normale asymptotique

Lorsque la taille de l’échantillon apparié demeure modeste, typiquement pour n’ inférieur ou égal à 20, le recours à une loi de probabilité continue théorique pour évaluer la significativité statistique s’avère inapproprié. Dans cette configuration, la distribution d’échantillonnage de la statistique W sous l’hypothèse nulle d’égalité distributionnelle est purement discrète et combinatoire. Sous H₀, chaque différence observée présente une probabilité égale à 1/2 d’être positive ou négative, indépendamment des autres paires. Il existe donc 2^n’ configurations de signes équiprobables. La distribution exacte de W est ainsi obtenue par dénombrement combinatoire exhaustif de toutes les permutations possibles des signes sur les rangs fixés de 1 à n’, permettant d’obtenir des p-values exactes sans aucune approximation numérique.

Cependant, dès lors que n’ augmente, le calcul combinatoire exhaustif devient computationnellement lourd (2⁵⁰ combinaisons dépasse 10¹⁵ opérations). En vertu du théorème central limite, la distribution d’échantillonnage de la somme des rangs W converge très rapidement vers une distribution normale lorsque n’ dépasse 20 à 25 paires d’observations. Sous l’hypothèse nulle, l’espérance mathématique E(W) et la variance Var(W) de la statistique sont dérivées formellement des propriétés combinatoires de la somme des entiers naturels :

E(W) = n'(n’ + 1) / 4

Var(W) = n'(n’ + 1)(2n’ + 1) / 24

Lorsque des rangs ex æquo sont présents dans les données, la variance théorique doit impérativement faire l’objet d’une correction afin d’ajuster la dispersion de la distribution de référence. La formule de la variance corrigée s’écrit alors :

Var(W) = [n'(n’ + 1)(2n’ + 1) – ∑_j (t_j³ – t_j) / 2] / 24

où t_j représente le nombre de valeurs ex æquo partageant le j-ième rang moyen. À partir de ces moments d’ordre un et deux, on calcule un score z standardisé. Afin de compenser le passage d’une variable discrète à saut unitaire à une densité normale continue, un facteur de correction de continuité de 0,5 est soustrait ou ajouté au numérateur, conférant à la statistique asymptotique la formulation suivante :

z = (|W – E(W)| – 0,5) / √(Var(W))

Ce score z suit asymptotiquement une loi normale centrée réduite N(0, 1), autorisant l’évaluation immédiate de la probabilité critique à l’aide des fonctions de répartition gaussiennes usuelles.

2. Test de Wilcoxon versus Test t de Student apparié : Cadre décisionnel en psychométrie

2.1 Sensibilité aux violations de la normalité

Le choix entre le test t de Student pour échantillons appariés et le test des rangs signés de Wilcoxon représente un dilemme classique en recherche expérimentale, tout particulièrement dans le champ de la psychométrie et des sciences cliniques. Le test t apparié repose fondamentalement sur le postulat que les différences intra-sujets (d_i = x_i,post – x_i,pre) proviennent d’une population dont la distribution suit une loi de Gauss. Lorsque cette hypothèse est violée de façon substantielle — notamment sous l’effet d’une asymétrie marquée (skewness prononcée) ou d’un aplatissement anomal (kurtosis excessive générant des queues épaisses) —, les propriétés inférentielles de l’estimateur de Student sont sévèrement compromises.

Les distributions leptokurtiques ou fortement asymétriques gonflent artificiellement la variance d’échantillonnage de la moyenne des différences, ce qui induit une instabilité du dénominateur du ratio t. En conséquence, le taux d’erreur de première espèce (α) peut dévier de manière imprévisible de son seuil nominal (souvent 5 %), tandis que le taux d’erreur de deuxième espèce (β) subit une inflation dramatique, anéantissant la capacité du test à rejeter l’hypothèse nulle lorsqu’elle est matériellement fausse. Le test des rangs signés de Wilcoxon, en convertissant les grandeurs numériques brutes en rangs ordinaux, neutralise totalement les déformations géométriques de la distribution sous-jacente.

Sur le plan théorique, l’efficacité relative asymptotique (ARE, Asymptotic Relative Efficiency) du test de Wilcoxon par rapport au test t de Student lorsque l’hypothèse de normalité parfaite est satisfaite s’établit à 3/π ≈ 0,955. Cela signifie que même dans le scénario théorique optimal pour le test t, le test de Wilcoxon ne cède que 4,5 % de puissance statistique par rapport à l’estimateur paramétrique optimal. En revanche, dès que la distribution réelle s’écarte de la gaussienne vers des lois à queues lourdes telles que la loi logistique, la loi de Laplace ou la loi de Cauchy, l’efficacité relative du test de Wilcoxon dépasse largement 1, atteignant parfois des valeurs infinies. Il s’avère donc immensément plus puissant que le test t face à des données empiriques non gaussiennes.

2.2 Résistance aux valeurs aberrantes (outliers)

L’un des défis majeurs de l’analyse de données psychologiques, psychiatriques et comportementales réside dans l’omniprésence d’observations aberrantes ou atypiques (outliers). Qu’elles résultent d’erreurs instrumentales de mesure, d’inattentions transitoires lors de la passation de questionnaires auto-rapportés, ou qu’elles reflètent une hétérogénéité clinique authentique (telle qu’une décompensation psychopathologique aiguë chez un sous-groupe restreint de patients), ces données extrêmes exercent une force de levier disproportionnée sur la moyenne arithmétique et l’écart-type des différences.

Le test t apparié calcule la statistique de contraste en intégrant la moyenne arithmétique des écarts au numérateur et l’écart-type d’échantillonnage au dénominateur. Il suffit d’une seule observation aberrante d’amplitude disproportionnée pour modifier le signe de la moyenne ou pour accroître l’écart-type à un point tel que la statistique t s’effondre, masquant un effet de traitement par ailleurs parfaitement robuste et généralisé sur les 99 % restants de la cohorte. À l’opposé, le test des rangs signés de Wilcoxon possède un point de rupture empirique et une insensibilité structurelle face aux scores extrêmes, dans la mesure où la valeur la plus colossale d’une série de données ne recevra jamais qu’un rang maximal égal à n’.

Pour illustrer cette résilience, comparons l’impact d’un sujet présentant une amélioration extrême de +80 points sur une échelle d’anxiété au sein d’une cohorte où l’amélioration médiane tourne autour de +5 points. Pour le test de Wilcoxon, cette différence colossale reçoit simplement le rang le plus élevé, disons 30 dans un groupe de 30 paires. Qu’elle vaille +15 ou +80 points, son rang demeure strictement inchangé (30), préservant l’intégrité de la somme W⁺. L’évaluation de la dispersion par l’intervalle interquartile (IQR) sur les différences confirme cette stabilité, là où l’écart-type paramétrique subit une explosion artificielle. Dans les contextes cliniques où les observations atypiques incarnent une réalité phénotypique qui ne doit pas être arbitrairement censurée par des filtres de nettoyage de données, le test de Wilcoxon offre une modélisation inférentielle fidèle et protectrice.

2.3 Nature des échelles de mesure en psychologie

L’argument épistémologique et méthodologique le plus profond en faveur du recours au test de Wilcoxon puise ses racines dans la taxonomie fondamentale des niveaux de mesure formulée par le psychologue Stanley Smith Stevens en 1946, distinguant les échelles nominales, ordinales, d’intervalles et de rapports. Dans la pratique contemporaine de l’évaluation psychologique et médicale, l’écrasante majorité des instruments diagnostiques repose sur des échelles de réponse ordonnées de type Likert (par exemple : « pas du tout d’accord », « modérément d’accord », « totalement d’accord ») ou des échelles visuelles analogiques discrétisées.

L’assimilation arbitraire d’échelles de Likert à des variables d’intervalles continues constitue un abus méthodologique fréquent. Rien ne permet d’affirmer mathématiquement que la distance psychologique séparant un score de 1 et un score de 2 sur un inventaire d’évaluation de la dépression est rigoureusement identique à l’intervalle cognitif séparant un score de 4 d’un score de 5. En l’absence d’une métrique d’intervalle prouvée, les opérations arithmétiques élémentaires telles que l’addition et la soustraction de scores bruts perdent leur justification formelle, privant le calcul des moyennes de sens mathématique strict. Le test des rangs signés de Wilcoxon contourne élégamment cet écueil : il n’exige que l’hypothèse ordinale selon laquelle les différences peuvent être classées de façon monotone.

De surcroît, les échelles psychométriques sont structurellement entravées par des bornes numériques infranchissables, générant des effets de plancher (impossibilité d’obtenir un score inférieur à zéro) ou de plafond (score maximal indépassable). Lorsqu’une cohorte clinique sévèrement déprimée entre dans un essai thérapeutique, nombre de participants débutent l’étude avec le score maximal autorisé par l’instrument. Les gains thérapeutiques consécutifs compriment la variance post-interventionnelle contre la frontière de l’échelle, invalidant définitivement toute prétention à la normalité gaussienne. L’arbre décisionnel suivant s’impose donc formellement en recherche appliquée :

• Si les scores bruts sont mesurés sur une échelle d’intervalle continue authentique, que la distribution des différences d_i vérifie formellement la normalité (asymétrie comprise entre -0,5 et +0,5, test de Shapiro-Wilk non significatif) et qu’aucun outlier n’est détecté : appliquer le test t de Student apparié.
• Si les données dérivent d’échelles ordinales ou de scores composites d’items de Likert, ou si la distribution des différences manifeste une asymétrie prononcée, des queues lourdes ou des valeurs extrêmes irréductibles : appliquer le test des rangs signés de Wilcoxon.
• Si la taille d’échantillon est excessivement restreinte (n < 15), la vérification de la normalité manquant cruellement de puissance : privilégier par défaut le test des rangs signés de Wilcoxon par mesure de précaution méthodologique.

3. Hypothèses statistiques et conditions d’application rigoureuses

3.1 Indépendance des paires et structure des données appariées

L’architecture logique du test des rangs signés de Wilcoxon exige une compréhension pointue de la structure de dépendance et d’indépendance gouvernant le plan d’échantillonnage. Le test est conçu exclusivement pour traiter des données dites appariées, ce qui signifie que chaque observation de la première condition est intrinsèquement liée à une observation unique de la seconde condition par l’intermédiaire d’une unité de mesure commune. Ce schéma correspond typiquement aux devis expérimentaux intra-sujets (mesures répétées sur les mêmes participants à des temps successifs T₁ et T₂) ou aux protocoles de paires appariées (par exemple, des jumeaux monozygotes, des couples mariés, ou des cohortes cas-témoins couplées individuellement selon l’âge et le sexe).

Bien que le couplage interne à chaque paire induise une corrélation intra-sujet assumée et exploitée par le calcul des différences d_i, l’hypothèse méthodologique incontournable stipule que les paires d’observations doivent être rigoureusement indépendantes les unes des autres. Mathématiquement, pour deux paires distinctes i et j (avec i ≠ j), la distribution conditionnelle de la paire (x_i,post, x_i,pre) ne doit apporter aucune information probabiliste sur la réalisation de la paire (x_j,post, x_j,pre). Cette indépendance mutuelle garantit que les rangs assignés aux écarts procèdent de tirages aléatoires décorrélés.

La violation de cette hypothèse d’indépendance inter-paires survient couramment lorsque les chercheurs omettent de contrôler les regroupements hiérarchiques (effets de grappes ou clustering). Par exemple, si l’évaluation de l’efficacité d’un protocole pédagogique est mesurée avant et après intervention sur des élèves regroupés au sein de différentes classes scolaires sous la tutelle d’enseignants distincts, les dynamiques groupales introduisent une covariance intraclasse. Dans une telle situation, traiter les élèves comme des unités indépendantes biaise drastiquement l’estimation de l’erreur-type de la statistique de rang et multiplie le risque d’inférence fallacieuse. Le statisticien doit veiller scrupuleusement à ce que l’échantillonnage inter-sujets soit aléatoire et dénué de corrélations structurelles non modélisées.

3.2 Hypothèse de symétrie de la distribution des différences

L’une des méconnaissances les plus délétères et répandues dans la littérature scientifique appliquée concerne l’objet exact de l’inférence du test des rangs signés de Wilcoxon. Nombre de manuels et d’articles biomédicaux affirment à tort que le test de Wilcoxon compare purement et simplement les médianes des deux conditions expérimentales. Cette assertion est mathématiquement inexacte sans l’adjonction d’un postulat géométrique contraignant : l’hypothèse de symétrie de la distribution de la variable de différence (d_i = x_i,post – x_i,pre) autour de sa médiane.

Pour comprendre cette subtilité théorique, il convient de disséquer le comportement de la statistique W. La statistique de Wilcoxon teste formellement l’hypothèse nulle selon laquelle la distribution des différences est symétrique autour de zéro. Dès lors, si la forme de la distribution sous-jacente des différences d_i est asymétrique dans la population parente, le rejet de l’hypothèse nulle n’implique pas nécessairement un déplacement pur de la médiane : il peut simplement résulter d’une asymétrie intrinsèque de la dispersion des scores sans changement de la tendance centrale. En revanche, si et seulement si l’on postule (ou vérifie empiriquement) que la distribution des écarts est approximativement symétrique, la médiane des différences coïncide rigoureusement avec la pseudo-médiane de Hodges-Lehmann, et le test devient alors un test exact de comparaison de décalage médian.

Il est capital de différencier le concept de décalage stochastique de l’égalité stricte des médianes. La pseudo-médiane évaluée par le test des rangs signés correspond à la médiane de toutes les moyennes par paires possibles, appelées paires de Walsh, formulées par :

(d_i + d_j) / 2   pour tout 1 ≤ i ≤ j ≤ n’

Lorsque la condition de symétrie des différences est violée de façon outrancière (par exemple, si la majorité des individus ne présente aucune variation ou un progrès infime mais qu’une minorité manifeste un effondrement spectaculaire), le test des rangs signés peut produire une significativité statistique alors même que la médiane des différences brutes est parfaitement égale à zéro. Le statisticien doit donc systématiquement inspecter la morphologie de la distribution des différences avant de conclure péremptoirement à une modification de la médiane de la population.

3.3 Définition des hypothèses nulle (H0) et alternative (H1)

La modélisation formelle de l’inférence statistique sous-jacente au test de Wilcoxon exige la définition rigoureuse des espaces d’hypothèses confrontés, en fonction des objectifs cliniques ou scientifiques formulés a priori dans le protocole de recherche. L’orientation de l’hypothèse alternative gouverne l’attribution des valeurs critiques et le calcul de la probabilité d’erreur.

Dans le cadre d’un test bilatéral (two-sided test), aucune orientation directionnelle préalable n’est postulée quant à l’effet de l’intervention. L’hypothèse nulle (H₀) et l’hypothèse alternative (H₁) s’expriment en termes de répartition stochastique des différences d = x_post – x_pre :

• H₀ : La distribution des différences est symétrique autour de zéro (médiane θ = 0). Il n’existe aucun décalage stochastique entre les deux conditions appariées.
• H₁ : La distribution des différences est décalée par rapport à zéro (médiane θ ≠ 0). La distribution de la condition post-intervention diffère systématiquement de celle de la condition pré-intervention.

À l’inverse, lorsque des considérations théoriques ou des données préliminaires solides justifient une hypothèse directionnelle univoque, le statisticien peut opter pour une formulation unilatérale (one-sided test). Si l’on évalue une thérapie cognitive visant explicitement la réduction d’une symptomatologie psychopathologique, l’attente clinique stipule que les scores post-thérapie doivent être substantiellement inférieurs aux scores initiaux, impliquant des écarts d_i = x_i,post – x_i,pre préférentiellement négatifs :

• H₀ : Le décalage de distribution est nul ou positif (θ ≥ 0).
• H₁ : Le décalage de distribution est strictement négatif (θ < 0), témoignant d’une décroissance des scores.

Inversement, dans le cadre de l’évaluation d’un gain d’apprentissage cognitif ou d’une amélioration de la qualité de vie, l’hypothèse unilatérale ciblera un décalage positif (θ > 0). Il convient de souligner avec fermeté un impératif d’intégrité scientifique : le choix entre une approche bilatérale ou unilatérale doit impérativement être consigné dans le plan d’analyse pré-expérimental (pré-enregistrement). Toute modification post-hoc de l’orientation de l’hypothèse après inspection visuelle des sorties numériques constitue une dérive méthodologique majeure (pratique assimilée au p-hacking) altérant le contrôle strict du taux d’erreur de première espèce.

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

4.1 Configuration des bibliothèques logicielles

La puissance du langage statistique R réside dans son architecture modulaire et son écosystème foisonnant de métapaquets développés par la communauté scientifique internationale. Bien que l’environnement de base de R dispose nativement des routines indispensables pour calculer le test de Wilcoxon via le package stats, la réalisation d’analyses poussées, incluant la manipulation élégante des tables, la détection des structures d’ex æquo, l’extraction automatisée des tailles d’effet et la génération de visuels sophistiqués, requiert l’importation de bibliothèques spécialisées.

Pour configurer un environnement analytique complet et pérenne, quatre packages majeurs doivent être chargés. Le premier est le métapaquet tidyverse (regroupant dplyr, tidyr, ggplot2 et purrr), qui forme l’infrastructure standard de manipulation et de visualisation des données selon le paradigme du tidy data. Le deuxième est rstatix, développé par Alboukadel Kassambara, qui offre une interface statistique fluide, conviviale et parfaitement intégrée aux pipelines tidyverse, standardisant les sorties d’analyse sous la forme de data frames structurés.

Le troisième package incontournable est coin (Conditional Inference Procedures in a Permutation Test Framework), maintenu par Torsten Hothorn et ses collaborateurs. Ce package représente la référence absolue pour le calcul exact et asymptotique conditionnel des tests de rangs en présence de structures complexes d’ex æquo et de données non standard. Enfin, le package psych fournit un ensemble robuste de métriques descriptives indispensables pour quantifier l’asymétrie, l’aplatissement et la tendance centrale non paramétrique des variables psychométriques. L’initialisation se programme de la façon suivante en console R :

install.packages(c(« tidyverse », « rstatix », « coin », « psych »))
library(tidyverse)
library(rstatix)
library(coin)
library(psych)

4.2 Formatage des données : Format large versus format long

Une source récurrente d’erreurs et de frustrations chez les analystes réside dans la confusion structurelle entre les représentations matricielles dites « format large » (wide format) et « format long » (long format ou tidy format). Chacune de ces géométries de données répond à des impératifs analytiques distincts dans R, et la maîtrise de leur transition bidirectionnelle via les verbes de tidyr est essentielle.

Dans le format large, chaque ligne de la table représente une unité d’échantillonnage unique (un participant individuel), et les mesures répétées collectées au fil du temps occupent des colonnes distinctes et juxtaposées. Typiquement, la table comportera une colonne id, une colonne score_pre et une colonne score_post. Ce format s’avère particulièrement intuitif pour le calcul vectoriel direct des différences individuelles via la formule arithmétique diff = score_post - score_pre et correspond à l’entrée par défaut de l’interface vectorielle de base de R.

À l’inverse, le format long empile verticalement l’ensemble des mesures dans une unique colonne de réponse numérique, identifiée par une colonne de facteur temporel indiquant la modalité d’observation. La table se compose dès lors d’une colonne d’identifiants id (chaque sujet apparaissant sur deux lignes distinctes), d’une variable indépendante intra-sujet temps (avec les niveaux « pre » et « post ») et d’une variable dépendante score. Ce format constitue le standard obligatoire pour l’interface par formule sous rstatix et pour l’ensemble des graphiques élaborés sous ggplot2.

La conversion d’un format à l’autre s’opère avec une remarquable élégance syntaxique grâce aux fonctions pivot_longer() et pivot_wider() de la bibliothèque tidyr. Par exemple, pour faire basculer un tableau donnees_larges vers le format long standardisé, on exécutera le code suivant :

donnees_longues <- donnees_larges %>%
 pivot_longer(
  cols = c(score_pre, score_post),
  names_to = « temps »,
  values_to = « score »
 ) %>%
 mutate(temps = factor(temps, levels = c(« score_pre », « score_post »)))

La réciproque s’effectue tout aussi aisément à l’aide de l’instruction inverse :

donnees_larges <- donnees_longues %>%
 pivot_wider(names_from = temps, values_from = score)

4.3 Contrôle de l’intégrité des données et gestion des valeurs manquantes

Préalablement à toute démarche d’inférence statistique, un audit rigoureux de l’intégrité matricielle s’impose afin d’éviter des corruptions insidieuses du calcul de la statistique de test. Dans le cadre d’un protocole à mesures répétées, l’existence d’identifiants dupliqués, d’identifiants orphelins ou d’inversions de lignes au sein des fichiers de saisie brute constitue un risque majeur qui invaliderait l’appariement des données.

La gestion des données manquantes unilatérales (phénomène classique d’attrition ou de perdus de vue lors de la phase de post-test) requiert une vigilance absolue. Si un participant dispose d’un score valide à l’évaluation initiale (T₁) mais n’a pas complété l’évaluation terminale (T₂), la paire devient incomplète. Contrairement aux modèles mixtes linéaires qui tolèrent les données manquantes sous l’hypothèse de perte aléatoire (MAR, Missing at Random), le calcul non paramétrique de Wilcoxon repose sur un appariement déterministe strict et applique obligatoirement une stratégie d’élimination par observation complète (listwise deletion).

Le chercheur doit s’assurer que chaque identifiant individuel présent dans le tableau possède exactement deux enregistrements non nuls au format long. Si l’on travaille directement sur des vecteurs issus du format large, la détection et le filtrage des paires incomplètes s’effectuent via l’instruction :

donnees_nettoyees <- donnees_larges %>%
 filter(!is.na(score_pre) & !is.na(score_post))

Enfin, il est impératif de contrôler minutieusement la typologie des vecteurs informatiques. Les variables quantitatives doivent être explicitement stockées sous forme de vecteurs numériques ou d’entiers (numeric ou integer), tandis que les variables temporelles et les identifiants de sujets doivent impérativement être formatés en facteurs catégoriels (factor). La confusion d’un identifiant numérique traité par erreur comme une variable quantitative continue peut fausser les routines de groupement automatisé des packages spécialisés.

5. Exploration préalable des données et vérification de la symétrie sous R

5.1 Calcul et inspection de la variable de différence

L’exploration préalable à l’inférence non paramétrique ne saurait se cantonner à la simple lecture des distributions marginales de la variable dépendante mesurée avant et après l’intervention. L’unité analytique fondamentale du test des rangs signés de Wilcoxon résidant dans la différence intra-individuelle, c’est la variable composite d_i = x_i,post – x_i,pre qui doit faire l’objet de toutes les investigations statistiques et géométriques.

À partir du format large, la création programmatique de la colonne des écarts s’opère par l’application de la fonction mutate() de dplyr :

donnees_analyse <- donnees_larges %>%
 mutate(difference = score_post – score_pre)

Une fois cette variable générée, le statisticien procède à l’extraction des indicateurs de tendance centrale et de dispersion non paramétriques. Contrairement à la moyenne et à la variance qui caractérisent les lois gaussiennes, les statistiques descriptives de choix sont ici la médiane, l’écart absolu médian (MAD, Median Absolute Deviation) et l’écart interquartile (IQR, Interquartile Range). L’implémentation sous R s’articule comme suit :

donnees_analyse %>%
 summarise(
  n = n(),
  mediane = median(difference),
  iqr = IQR(difference),
  mad = mad(difference),
  min = min(difference),
  max = max(difference)
 )

L’inspection comparative de la médiane et de l’intervalle interquartile permet d’apprécier d’emblée la force du changement typique subi par l’échantillon tout en neutralisant l’influence des variations extrêmes individuelles décelables par l’examen des valeurs minimale et maximale.

5.2 Évaluation formelle de la non-normalité

Pour justifier rigoureusement l’abandon du test t de Student apparié au profit du test de Wilcoxon lors de la rédaction d’un manuscrit scientifique, le chercheur doit étayer sa décision par une évaluation formelle du comportement distributionnel des différences. L’outil inférentiel standard le plus éprouvé et puissant pour tester la normalité réside dans le test de Shapiro-Wilk, exécuté dans R via la fonction shapiro.test() :

shapiro.test(donnees_analyse$difference)

L’hypothèse nulle du test de Shapiro-Wilk postule que l’échantillon des différences provient d’une population gaussienne. L’obtention d’une valeur critique p inférieure au seuil conventionnel de 0,05 conduit au rejet sans équivoque de l’hypothèse de normalité, attestant que la distribution empirique dévie de manière significative du modèle théorique de Student. Cependant, l’analyste expérimenté doit garder à l’esprit les limites inhérentes aux tests formels de normalité : sur de très grands échantillons (n > 300), le test de Shapiro-Wilk devient hypersensible à des micro-déviations triviales dépourvues de toute conséquence statistique pratique ; à l’inverse, sur de petits échantillons (n < 20), il manque cruellement de puissance pour détecter des non-normalités substantielles.

Cette approche inférentielle doit donc être impérativement complétée par l’examen quantitatif des coefficients d’asymétrie (skewness) et d’aplatissement (kurtosis), calculables via la commande describe() du package psych :

psych::describe(donnees_analyse$difference)

Des valeurs de coefficient d’asymétrie s’écartant de l’intervalle [-1, +1] ou un aplatissement excédant l’intervalle [-2, +2] scellent définitivement la non-conformité aux exigences paramétriques, rendant l’emploi du test des rangs signés de Wilcoxon mathématiquement indispensable.

5.3 Inspection graphique de la forme de la distribution

L’exploration des données ne saurait atteindre sa pleine maturité analytique sans une corroboration visuelle fine des propriétés morphologiques de la variable d’écart. C’est l’inspection graphique qui permettra d’évaluer non seulement l’étendue des déviations par rapport à la droite de normalité, mais surtout de valider le postulat fondamental de symétrie relative des différences autour de leur médiane.

La première composante graphique essentielle est le diagramme quantile-quantile normal (Q-Q plot), qui confronte directement les quantiles empiriques de la distribution des différences aux quantiles théoriques d’une loi gaussienne standard. Sous ggplot2, ce tracé s’obtient avec une précision chirurgicale :

ggplot(donnees_analyse, aes(sample = difference)) +
 stat_qq(color = « #2c3e50 », size = 2) +
 stat_qq_line(color = « #e74c3c », linetype = « dashed », linewidth = 1) +
 theme_minimal() +
 labs(title = « Diagramme Quantile-Quantile des Différences »,
   x = « Quantiles Théoriques », y = « Quantiles Observés »)

Des points s’écartant systématiquement de la ligne diagonale rouge en formant une courbure en « S » signalent des queues lourdes, tandis qu’une concavité unilatérale traduit une asymétrie nette. Pour inspecter la condition de symétrie de la distribution, on superpose un histogramme à pas fin et une estimation de densité par noyau (kernel density estimation) sur laquelle est apposée une ligne verticale matérialisant la médiane :

ggplot(donnees_analyse, aes(x = difference)) +
 geom_histogram(aes(y = after_stat(density)), bins = 15, fill = « #3498db », color = « white », alpha = 0.6) +
 geom_density(color = « #2980b9 », linewidth = 1.2) +
 geom_vline(xintercept = median(donnees_analyse$difference), linetype = « dashed », color = « #c0392b », linewidth = 1) +
 theme_minimal() +
 labs(title = « Distribution de Densité des Écarts Intra-Sujets »,
   x = « Différence (Post – Pré) », y = « Densité »)

Si la cloche de densité présente une forme approximativement équilibrée de part et d’autre de la ligne rouge médiane, l’analyste dispose de l’assurance méthodologique que le test de Wilcoxon pourra être interprété sans réserve comme un test formel de décalage de la médiane.

6. Exécution pas-à-pas avec la fonction native wilcox.test() dans R

6.1 Syntaxe de base pour échantillons appariés

L’implémentation historique du test des rangs signés dans l’environnement de base R repose sur la fonction polyvalente wilcox.test(), localisée au sein du package stats. Cette fonction est remarquable car elle prend en charge aussi bien le test de somme des rangs de Mann-Whitney (pour échantillons totalement indépendants) que le test des rangs signés de Wilcoxon (pour échantillons appariés). Cette double compétence impose une rigueur absolue lors du paramétrage programmatique : l’argument paired doit impérativement être configuré sur TRUE.

L’exécution peut s’opérer selon deux modes syntaxiques distincts. Le premier, la syntaxe vectorielle, s’articule autour de deux vecteurs numériques distincts extraits du format large :

wilcox.test(donnees_larges$score_post, donnees_larges$score_pre, paired = TRUE)

Le second, la syntaxe par formule, s’applique directement sur une table structurée au format long, en précisant la variable dépendante à gauche du tilde et la variable indépendante catégorielle à droite, couplée à l’identification explicite de la structure de bloc :

wilcox.test(score ~ temps, data = donnees_longues, paired = TRUE)

L’omission accidentelle du paramètre paired = TRUE représente le piège le plus funeste du traitement statistique sous R. Si cette option fait défaut, R applique silencieusement le test non apparié de Mann-Whitney (équivalent au test de la somme des rangs de Wilcoxon), traitant les deux séries temporelles comme s’il s’agissait de deux groupes distincts de participants indépendants. Cette bévue anéantit totalement le gain de puissance lié au contrôle de la variance intra-sujet et produit une p-value erronée susceptible de compromettre les conclusions de l’étude.

6.2 Gestion des hypothèses unilatérales et bilatérales

Par défaut, la fonction wilcox.test() calibre l’inférence sur une hypothèse bilatérale, modélisée par l’argument interne alternative = "two.sided". Toutefois, lorsque les hypothèses de recherche pré-définies postulent une orientation unidirectionnelle du changement, l’analyste doit configurer explicitement l’argument alternative en choisissant entre les options "greater" ou "less".

La compréhension de l’orientation du contraste exige une attention extrême à l’ordre d’entrée des arguments vectoriels dans l’appel de fonction. Dans l’instruction vectorielle wilcox.test(x, y, paired = TRUE, alternative = ...), le calcul mathématique interne évalue les différences sous la forme d = x – y. Par conséquent :

• Si l’on spécifie alternative = "greater", R teste l’hypothèse alternative selon laquelle la distribution de x est stochastiquement décalée vers des valeurs plus élevées que y, ce qui équivaut à postuler que la différence (x – y) est positive.
• Si l’on spécifie alternative = "less", R teste l’hypothèse alternative selon laquelle la distribution de x est stochastiquement décalée vers des valeurs plus faibles que y, ce qui équivaut à postuler que la différence (x – y) est négative.

La syntaxe par formule (score ~ temps) présente un risque d’inversion directionnelle accru, car l’ordre de soustraction dépend strictement de l’ordonnancement interne des niveaux du facteur catégoriel temps déterminé par la fonction levels(). Si le premier niveau est « pre » et le second « post », R soustraira implicitement « pre » de « post », ou inversement selon la version de l’interpréteur. Pour bannir tout risque d’ambiguïté directionnelle, il est fortement recommandé de vérifier l’ordre des facteurs via levels(donnees_longues$temps) préalablement à l’analyse unilatérale, ou d’utiliser préférentiellement la syntaxe vectorielle dont la relation x – y demeure parfaitement explicite.

6.3 Compréhension de la sortie console de base R

L’exécution de la commande wilcox.test() génère dans la console une sortie textuelle condensée structurée en plusieurs lignes informatives :

Wilcoxon signed rank test
data: donnees_larges$score_post and donnees_larges$score_pre
V = 435.5, p-value = 0.001247
alternative hypothesis: true location shift is not equal to 0

Le premier élément d’intérêt est la statistique de test, désignée par la lettre majuscule V. Dans l’architecture algorithmique de R, V ne représente pas le minimum des sommes des rangs min(W⁺, W^–) souvent décrit dans les manuels scolaires d’initiation, mais correspond rigoureusement et invariablement à la somme des rangs attribués aux différences strictement positives (W⁺). La valeur numérique associée (ici 435.5) doit donc s’interpréter en relation directe avec la somme totale théorique des rangs disponibles pour la taille d’échantillon analysée.

La seconde ligne critique affiche la p-value, qui traduit la probabilité d’observer une somme de rangs positifs au moins aussi extrême que celle mesurée sous l’hypothèse nulle d’absence de décalage stochastique. Cependant, la sortie console s’accompagne très fréquemment de messages d’alerte (warnings) qui suscitent l’inquiétude légitime des chercheurs :

Warning message: In wilcox.test.default(...): cannot compute exact p-value with ties
Warning message: In wilcox.test.default(...): cannot compute exact p-value with zeroes

Le message relatif aux « ties » signale l’existence d’ex æquo dans les valeurs absolues des différences, ce qui brise la structure combinatoire des entiers continus et interdit à l’algorithme natif d’évaluer la distribution exacte par permutations théoriques. R bascule alors par défaut vers une approximation normale asymptotique sans que l’utilisateur n’en maîtrise pleinement les modalités. Le message relatif aux « zeroes » indique la présence d’individus sans aucun changement entre le pré et le post-test, lesquels sont éliminés de l’estimation exacte par la méthode de Wilcoxon. Face à ces avertissements récurrents sur les données empiriques réelles, le recours à des bibliothèques plus élaborées devient impératif.

7. Gestion avancée des ex æquo et calcul exact avec des packages spécialisés

7.1 Problématique théorique des rangs ex æquo et des différences nulles

L’émergence de rangs ex æquo (ties) et de différences rigoureusement nulles (d_i = 0) n’est pas une anomalie marginale en sciences empiriques, mais une régularité presque universelle dès lors que les instruments d’évaluation reposent sur des échelles discrètes à granularité limitée. Ces configurations imposent des arbitrages mathématiques lourds de conséquences sur le calcul de la probabilité critique.

Face à des différences nulles, deux paradigmes concurrents coexistent dans la théorie statistique. Le premier, formulé par Wilcoxon en 1945, préconise la suppression pure et simple de toutes les paires présentant une différence d_i = 0, réduisant la taille de l’échantillon à n’. Cette approche a été sévèrement critiquée car elle écarte arbitrairement des données empiriques authentiques attestant d’une absence absolue d’effet thérapeutique. Le second paradigme, formulé par John W. Pratt en 1959, inclut l’ensemble des différences nulles dans la procédure de classement par rangs : les zéros reçoivent les rangs les plus faibles, puis sont omis uniquement lors de la sommation finale des rangs positifs et négatifs. La méthode de Pratt préserve ainsi le poids relatif des écarts de faible magnitude par rapport à la taille réelle de la cohorte testée.

Concernant les valeurs absolues identiques, la procédure d’attribution des rangs moyens (mid-ranks) induit une réduction systématique de la variance de la statistique W. Si cette contraction de la variance n’est pas corrigée formellement dans le calcul asymptotique de la statistique standardisée z, l’erreur-type s’avère surestimée, entraînant un test indûment conservateur et une perte notable de puissance statistique. Par ailleurs, les distributions combinatoires exactes tabulées dans les ouvrages classiques deviennent totalement inopérantes en présence d’ex æquo, imposant l’usage d’algorithmes de permutation informatisés.

7.2 Mise en œuvre robuste avec le package coin

Pour résoudre définitivement les limitations et les approximations aveugles de la fonction de base, la communauté biostatistique recourt préférentiellement au package coin. Ce dernier intègre les tests de rangs dans la théorie unifiée de l’inférence conditionnelle par permutation développée par Strasser et Weber (1999). La fonction dédiée au test apparié s’intitule wilcoxsign_test().

Sa force réside dans sa capacité à ajuster avec une précision chirurgicale la distribution de référence de la statistique de test, même en présence d’un nombre massif d’ex æquo et de zéros traités selon la méthode de Pratt ou de Wilcoxon. L’appel s’exécute exclusivement sous la syntaxe par formule :

wilcoxsign_test(score_post ~ score_pre, data = donnees_larges, zero.method = « Pratt », distribution = « exact »)

L’argument distribution offre une flexibilité algorithmique absolue :

• distribution = "exact" : commande le calcul combinatoire exact conditionné spécifiquement sur le schéma d’ex æquo présent dans le jeu de données analysé. Cette option garantit une p-value sans aucune approximation numérique, même en présence d’égalités multiples, tant que la combinatoire demeure gérable.
• distribution = "approximate"(nresample = 100000) : déploie une procédure de rééchantillonnage de Monte Carlo avec 100 000 tirages par permutation aléatoire lorsque la taille d’échantillon interdit le calcul exact exhaustif, fournissant une estimation asymptotiquement parfaite assortie d’un intervalle de confiance sur la p-value.
• distribution = "asymptotic" : applique la correction de continuité et l’ajustement complet de la variance pour les rangs moyens selon les formulations théoriques les plus abouties.

L’extraction de la statistique Z standardisée s’effectue directement par la méthode statistic() appliquée à l’objet résultant, délivrant la valeur exacte corrigée indispensable pour le calcul ultérieur des tailles d’effet standardisées :

res_coin <- wilcoxsign_test(score_post ~ score_pre, data = donnees_larges, distribution = « exact »)
statistic(res_coin)

7.3 Workflow standardisé via rstatix::wilcox_test

Si le package coin brille par sa sophistication mathématique, le package rstatix s’impose quant à lui comme le standard ergonomique incontournable pour les scientifiques de données opérant au sein du tidyverse. La fonction wilcox_test() de rstatix encapsule les routines statistiques tout en retournant systématiquement un objet de type tibble (data frame moderne), simplifiant grandement l’intégration des résultats dans des chaînes de traitement automatisées.

La syntaxe de rstatix s’intègre harmonieusement avec l’opérateur de tuyau (pipe %>%) :

resultat_rstatix <- donnees_longues %>%
 wilcox_test(score ~ temps, paired = TRUE, detailed = TRUE)

L’argument detailed = TRUE enrichit le tableau de sortie en fournissant une décomposition analytique complète comprenant : la variable dépendante (.y.), les groupes comparés (group1, group2), la statistique numérique brute (statistic), la taille effective d’échantillon (n1, n2), l’estimation de l’amplitude du décalage (estimate correspondant à la pseudo-médiane de Hodges-Lehmann), l’intervalle de confiance associé et la p-value rigoureusement calculée.

De surcroît, rstatix excelle dans la gestion des analyses stratifiées. Si l’évaluation clinique pré-post a été administrée à travers différents centres hospitaliers ou sous-groupes diagnostiques, l’adjonction de l’instruction group_by() permet d’exécuter simultanément l’ensemble des tests en une ligne de code limpide :

donnees_longues %>%
 group_by(centre_hospitalier) %>%
 wilcox_test(score ~ temps, paired = TRUE) %>%
 adjust_pvalue(method = « bonferroni »)

La fonction adjust_pvalue() enchaîne instantanément l’ajustement pour comparaisons multiples (correction de Bonferroni, de Holm ou de Benjamini-Hochberg), assurant un contrôle rigoureux du taux d’erreur global au sein des protocoles multicentriques.

8. Mesure de la taille d’effet : Calcul et interprétation

8.1 La corrélation de rang bisériale appariée

L’évaluation moderne de l’inférence statistique selon les recommandations méthodologiques contemporaines (notamment portées par l’American Psychological Association) proscrit formellement la dépendance exclusive envers la seule p-value. Une probabilité critique fortement significative indique uniquement qu’il est improbable d’observer un tel décalage de distribution sous le modèle nul ; elle ne renseigne en rien sur l’amplitude clinique ou pratique de l’effet mesuré. Il est donc indispensable d’adjoindre systématiquement un estimateur standardisé de la taille d’effet (effect size).

Dans le paradigme non paramétrique des rangs signés, la mesure d’amplitude conceptuellement la plus pure et cohérente est le coefficient de corrélation de rang bisériale appariée (matched-pairs rank biserial correlation), traditionnellement noté r_rb ou r_c, introduit par Kerby en 2014. Ce coefficient découle directement du principe de la différence proportionnelle simple entre les rangs favorables et défavorables :

r_c = (W⁺ – W^–) / (W⁺ + W^–)

Dans la mesure où la somme totale des rangs W⁺ + W^– est égale à la quantité fixe n'(n’ + 1) / 2, la corrélation bisériale des rangs s’interprète comme la différence nette entre la proportion des rangs attestant d’un gain post-interventionnel et la proportion des rangs traduisant une détérioration. Ce coefficient varie strictement dans l’intervalle borné [-1, +1] :

• Une valeur de +1,0 traduit une amélioration universelle et absolue : l’ensemble des différences observées dans l’échantillon est strictement positif (W^– = 0).
• Une valeur de 0,0 reflète un équilibre parfait : la somme des rangs positifs compense rigoureusement la somme des rangs négatifs, attestant d’une absence complète d’orientation privilégiée.
• Une valeur de -1,0 dénote une dégradation unanime : la totalité des variations observées s’avère négative (W⁺ = 0).

La programmation de ce coefficient sous R s’exécute de façon élémentaire à partir de la sortie native de base :

calcul_rc <- function(wilcox_obj, n) {
 W_plus <- wilcox_obj$statistic
 W_total <- n * (n + 1) / 2
 W_minus <- W_total – W_plus
 rc <- (W_plus – W_minus) / W_total
 return(as.numeric(rc))
}

8.2 Le coefficient r standardisé basé sur le score Z

L’indicateur de taille d’effet le plus largement diffusé dans la littérature biomédicale et psychométrique demeure toutefois le coefficient standardisé r formalisé par le psychologue statisticien Robert Rosenthal (1991). Cet estimateur dérive directement de la conversion de la statistique standardisée asymptotique Z en une métrique de corrélation équivalente :

r = |Z| / √N

Bien que cette formulation apparaisse d’une simplicité désarmante, elle est le théâtre d’un débat méthodologique persistant dans la littérature quant à la définition exacte du dénominateur N. Deux écoles doctrinales s’opposent :

• La convention des observations totales : N représente le nombre cumulé de points d’observation individuels traités par l’analyse, soit le double du nombre de paires (N = 2 × n). Cette approche découle de l’extension directe du test de Mann-Whitney pour échantillons indépendants où N = n₁ + n₂.
• La convention des paires appariées : N représente le nombre strict de paires d’observations appariées (N = n). Dans la mesure où l’unité statistique élémentaire du test des rangs signés est la différence intra-individuelle d_i, l’échantillon réel soumis à l’inférence se compose formellement de n différences et non de 2n observations indépendantes.

Le package rstatix intègre nativement cette problématique via la fonction wilcox_effsize(). Par défaut, celle-ci applique rigoureusement la convention des paires appariées dès lors que l’argument paired = TRUE est spécifié, évitant les sous-estimations de la taille d’effet consécutives à un doublement artificiel du dénominateur :

donnees_longues %>%
 wilcox_effsize(score ~ temps, paired = TRUE)

Cette instruction génère instantanément un tibble documentant le coefficient r, l’intervalle de confiance calculé par rééchantillonnage et la classification qualitative de la magnitude du contraste.

8.3 Grilles d’interprétation et intervalles de confiance

Pour qualifier qualitativement l’amplitude de l’effet identifié par le coefficient r de Rosenthal, les chercheurs s’appuient historiquement sur les seuils heuristiques conventionnels formulés par Jacob Cohen (1988) dans le cadre des corrélations de Pearson, transposés par commodité aux statistiques de rangs :

• Effet faible : |r| compris entre 0,10 et 0,30. L’effet de l’intervention est modeste et peut échapper à l’observation clinique individuelle non outillée.
• Effet moyen : |r| compris entre 0,30 et 0,50. Le changement est substantiel, perceptible et cliniquement pertinent.
• Effet fort : |r| supérieur ou égal à 0,50. L’intervention induit une modification majeure et prépondérante de la distribution des réponses.

Toutefois, la restitution académique rigoureuse d’une taille d’effet ponctuelle impose de l’assortir de son intervalle de confiance à 95 % (IC 95 %), attestant de la précision de l’estimation d’échantillonnage. En l’absence de distribution paramétrique standard simple pour le coefficient r de rangs signés, la méthode de choix pour extraire cet intervalle est le bootstrap non paramétrique (rééchantillonnage avec remplacement). Le package boot offre une infrastructure idéale pour ce calcul :

library(boot)
fonction_boot_r <- function(data, indices) {
 d <- data[indices, ]
 stat <- wilcox.test(d$score_post, d$score_pre, paired = TRUE)$statistic
 n <- nrow(d)
 z <- qnorm(wilcox.test(d$score_post, d$score_pre, paired = TRUE)$p.value / 2, lower.tail = FALSE)
 return(z / sqrt(n))
}
boot_res <- boot(data = donnees_larges, statistic = fonction_boot_r, R = 2000)
boot.ci(boot_res, type = « perc »)

Enfin, dans les publications orientées vers la pratique médicale, il est particulièrement pertinent de compléter le reporting par la mesure de taille d’effet en langage commun (Common Language Effect Size ou probabilité de supériorité). Celle-ci exprime simplement le pourcentage empirique de patients ayant manifesté une amélioration stricte (d_i > 0) au sein de la cohorte suivie.

9. Cas pratique complet : Évaluation pré-post d’une intervention psychologique

9.1 Présentation du jeu de données et contextualisation clinique

Afin d’ancrer de manière pragmatique les concepts théoriques développés jusqu’ici, développons une simulation clinique intégrale, reproductible et scientifiquement représentative. Considérons un protocole de recherche translationnelle en psychiatrie clinique évaluant l’efficacité d’une nouvelle thérapie cognitive et émotionnelle brève (8 semaines de prise en charge) destinée à réduire l’anxiété chez des patients présentant un trouble anxieux généralisé (TAG).

Le niveau d’anxiété est quantifié à l’aide de l’inventaire d’anxiété état-trait de Spielberger (score STAI, State-Trait Anxiety Inventory, forme Y-A), administré lors du bilan pré-interventionnel (T₁) puis répété rigoureusement chez les mêmes patients au terme du protocole thérapeutique (T₂). Le score total à cette échelle varie de 20 à 80 points, les valeurs les plus élevées signalant une symptomatologie anxieuse sévère. Nous introduisons délibérément dans le processus de génération des caractéristiques réalistes courantes : une asymétrie de distribution, des rangs ex æquo et un sujet réfractaire dont l’anxiété stagne.

Le script suivant assure la génération du jeu de données synthétique reproductible grâce à la fixation de la graine aléatoire :

set.seed(2026)
n_patients <- 25

score_pre <- round(rgamma(n_patients, shape = 12, rate = 0.25))
score_pre <- pmin(pmax(score_pre, 20), 80)

reduction_therapeutique <- round(rnorm(n_patients, mean = 7, sd = 4))
reduction_therapeutique[3] <- 0
reduction_therapeutique[12] <- -2

score_post <- score_pre – reduction_therapeutique
score_post <- pmin(pmax(score_post, 20), 80)

etude_clinique <- tibble(
 id = factor(paste0(« PAT_ », sprintf(« %02d », 1:n_patients))),
 stai_pre = score_pre,
 stai_post = score_post
)

9.2 Déroulement exhaustif du script d’analyse R

Disposant de notre cohorte clinique simulée au format large dans la table etude_clinique, nous déroulons l’intégralité du pipeline d’analyse statistique : calcul de la variable de différence, exploration non paramétrique descriptive, confirmation de la rupture de normalité et exécution croisée des fonctions de test.

etude_clinique <- etude_clinique %>%
 mutate(stai_diff = stai_post – stai_pre)

# Statistiques descriptives de la différence
desc_diff <- etude_clinique %>%
 summarise(
  mediane = median(stai_diff),
  q25 = quantile(stai_diff, 0.25),
  q75 = quantile(stai_diff, 0.75),
  iqr = IQR(stai_diff),
  mad = mad(stai_diff)
 )
print(desc_diff)

# Vérification formelle de l’hypothèse de normalité
shapiro_resultat <- shapiro.test(etude_clinique$stai_diff)
print(shapiro_resultat)

# Exécution du test de Wilcoxon apparié unilatéral (H1: réduction d’anxiété -> post < pre)
# Utilisation de la fonction native
test_base <- wilcox.test(etude_clinique$stai_post, etude_clinique$stai_pre,
          paired = TRUE, alternative = « less », exact = FALSE)
print(test_base)

# Modélisation exacte conditionnelle avec le package coin
test_coin <- wilcoxsign_test(stai_post ~ stai_pre, data = etude_clinique,
           alternative = « less », distribution = « exact »)
print(test_coin)

# Calcul standardisé de la taille d’effet avec rstatix
donnees_cliniques_long <- etude_clinique %>%
 pivot_longer(cols = c(stai_pre, stai_post), names_to = « temps », values_to = « stai_score ») %>%
 mutate(temps = factor(temps, levels = c(« stai_pre », « stai_post »)))

taille_effet <- donnees_cliniques_long %>%
 wilcox_effsize(stai_score ~ temps, paired = TRUE, alternative = « greater »)
print(taille_effet)

9.3 Interprétation substantive des sorties numériques

L’exécution de la séquence analytique livre des enseignements quantitatifs d’une grande limpidité clinique. L’examen des indicateurs descriptifs révèle une médiane des différences de -7,00 points (intervalle interquartile [IQR] = 5,00, s’étendant du premier quartile Q₁ = -9,00 au troisième quartile Q₃ = -4,00). Cette métrique atteste que 50 % des patients traités ont bénéficié d’un allègement de leur symptomatologie anxieuse supérieur ou égal à 7 points sur l’inventaire STAI.

Le test de normalité de Shapiro-Wilk appliqué sur les écarts individuels affiche un coefficient W = 0,891 assorti d’une p-value égale à 0,011. Ce résultat se situant nettement en deçà du seuil critique standard (α = 0,05), l’hypothèse de normalité de la distribution des différences est formellement rejetée, confirmant l’inadéquation mathématique du test t de Student apparié sur cette cohorte.

L’inférence non paramétrique délivrée par wilcoxsign_test() du package coin (respectant la structure discrète exacte) produit une statistique standardisée Z = -4,013 avec une probabilité critique exacte p < 0,0001. Face à cette valeur infiniment inférieure au seuil de risque choisi, la décision méthodologique s’impose sans équivoque : l’hypothèse nulle (H₀) est catégoriquement rejetée en faveur de l’hypothèse alternative (H₁). Le protocole thérapeutique induit une réduction statistiquement significative et systématique des scores d’anxiété post-interventionnels.

L’analyse de la taille d’effet parachève l’évaluation substantive : la fonction wilcox_effsize() retourne un coefficient Rosenthal r = 0,803 (considéré comme un effet d’amplitude majeure selon les critères de Cohen). L’examen individuel confirme l’homogénéité de la réponse thérapeutique : sur les 25 patients de la cohorte, 23 ont manifesté un progrès net (différences négatives), un patient a maintenu un statut invariant (différence nulle), et un seul patient a présenté une légère aggravation (+2 points). Cette convergence entre significativité statistique et répercussion clinique valide pleinement l’efficacité de la prise en charge investiguée.

10. Visualisation publication-ready des résultats sous ggplot2

10.1 Graphique des pentes individuelles (Paired Slopegraph)

Dans l’univers des publications scientifiques contemporaines, les diagrammes en barres classiques accompagnés de barres d’erreurs (souvent surnommés de façon péjorative dynamite plots) sont de plus en plus proscrits par les comités éditoriaux des revues majeures. Dans le cadre de mesures répétées appariées, ces graphiques agrégés masquent totalement la richesse de l’information intra-individuelle : ils échouent à montrer si l’effet moyen résulte d’une amélioration globale de l’ensemble de la cohorte ou d’un changement drastique concentré chez une poignée de répondants extrêmes.

Le graphique des pentes individuelles (paired slopegraph) constitue la représentation visuelle par excellence pour sublimer les données appariées. Il relie par un segment rectiligne le score pré et le score post de chaque sujet, matérialisant instantanément l’hétérogénéité des trajectoires thérapeutiques. Sous ggplot2, ce graphique s’élabore à partir du format long en utilisant l’esthétique de groupement par identifiant :

ggplot(donnees_cliniques_long, aes(x = temps, y = stai_score, group = id)) +
 geom_line(aes(color = stai_score), alpha = 0.5, linewidth = 0.8) +
 geom_point(aes(color = stai_score), size = 2.5, alpha = 0.8) +
 scale_color_viridis_c(option = « plasma », direction = -1) +
 scale_x_discrete(labels = c(« stai_pre » = « Pré-Intervention », « stai_post » = « Post-Intervention »)) +
 theme_classic(base_size = 14) +
 theme(legend.position = « none »,
    axis.title.x = element_blank(),
    plot.title = element_text(face = « bold », hjust = 0.5)) +
 labs(title = « Trajectoires Individuelles de l’Anxiété (STAI) »,
   y = « Score Global d’Anxiété »)

Pour décupler la puissance pédagogique de la figure, il est aisé d’ajouter une ligne de trajectoire médiane globale, dessinée dans une teinte vive et une épaisseur supérieure, permettant de superposer le signal moyen de la population à la cartographie détaillée des réponses individuelles.

10.2 Diagramme en violon apparié et nuage de pluie (Raincloud plot)

Pour offrir aux relecteurs une vision holistique associant à la fois la forme complète des distributions marginales, la localisation des quantiles et l’ancrage brut des points de données, le diagramme en nuage de pluie (raincloud plot) représente l’état de l’art du graphisme statistique biomédical. Il combine simultanément une demi-courbe de densité en violon (représentant la nébulosité de la distribution), une boîte à moustaches étroite (indiquant la médiane et l’IQR) et un nuage de points décalé (la pluie de données).

Cette visualisation se programme avec une remarquable flexibilité sous ggplot2 :

ggplot(donnees_cliniques_long, aes(x = temps, y = stai_score, fill = temps)) +
 geom_violin(trim = FALSE, alpha = 0.4, color = NA, position = position_dodge(0.8)) +
 geom_boxplot(width = 0.15, fill = « white », color = « black », outlier.shape = NA, alpha = 0.8) +
 geom_jitter(width = 0.08, size = 2, alpha = 0.6, aes(color = temps)) +
 scale_fill_manual(values = c(« stai_pre » = « #34495e », « stai_post » = « #16a085 »)) +
 scale_color_manual(values = c(« stai_pre » = « #2c3e50 », « stai_post » = « #1abc9c »)) +
 scale_x_discrete(labels = c(« stai_pre » = « Pré-Traitement », « stai_post » = « Post-Traitement »)) +
 theme_minimal(base_size = 14) +
 theme(legend.position = « none »,
    panel.grid.minor = element_blank(),
    plot.title = element_text(face = « bold », hjust = 0.5)) +
 labs(title = « Densité et Dispersion des Scores d’Anxiété »,
   y = « Score STAI », x = «  »)

Ce type de figure expose immédiatement la compression de la variance et le basculement stochastique de la distribution des scores post-interventionnels vers les valeurs physiologiques d’anxiété légère.

10.3 Intégration des annotations statistiques avec ggpubr

L’étape finale de la conception graphique académique consiste à intégrer directement sur la figure la barre de comparaison horizontale indiquant la significativité statistique formelle calculée par le test des rangs signés de Wilcoxon, évitant ainsi le recours laborieux à des logiciels tiers d’édition vectorielle. Le package ggpubr offre la fonction spécialisée stat_compare_means() conçue pour automatiser cette jonction.

Il est fondamental de veiller au paramétrage explicite des arguments méthodologiques internes afin d’éviter l’application d’un test non apparié par défaut :

library(ggpubr)

ggpaired(donnees_cliniques_long, x = « temps », y = « stai_score »,
     color = « temps », palette = c(« #e74c3c », « #27ae60 »),
     line.color = « gray70 », line.size = 0.4,
     xlab = « Étape de l’Évaluation », ylab = « Score STAI »,
     title = « Efficacité Thérapeutique sur l’Anxiété ») +
 stat_compare_means(method = « wilcox.test », paired = TRUE,
           label.x = 1.35, label.y = 75,
           aes(label = paste0(« Wilcoxon p = « , after_stat(p.format)))) +
 theme_pubr(base_size = 14) +
 theme(legend.position = « none »)

La figure vectorielle résultante peut ensuite être exportée dans une résolution native conforme aux exigences des éditeurs scientifiques internationaux (format TIFF non compressé à 300 ou 600 DPI, ou format PDF vectoriel sans pixellisation) via la commande standard :

ggsave(« figure_wilcoxon_stai.pdf », width = 8, height = 6, dpi = 300)

11. Rédaction académique des résultats aux normes APA (7e édition)

11.1 Règles formelles de notation des statistiques non paramétriques

La transcription des résultats statistiques dans les sections « Résultats » des articles scientifiques doit satisfaire aux directives éditoriales très précises dictées par la 7^e édition du manuel de publication de l’American Psychological Association (APA). L’inobservation de cette charte graphique et notationnelle discrédite la rigueur du manuscrit lors de la phase de révision par les pairs (peer review).

Pour le test des rangs signés de Wilcoxon, les règles formelles de mise en forme typographique imposent :

• Les symboles statistiques latins doivent impérativement être mis en italique : T, W, V, Z, p, r, N, n. En revanche, les indices ou sous-scripts demeurent en caractères romains droits.
• Le rapport de tendance centrale paramétrique (moyenne et écart-type, M et SD) doit être banni au profit des indicateurs d’ordre : la médiane (Mdn) et l’intervalle interquartile (IQR) doivent être mentionnés pour chacune des conditions d’évaluation ainsi que pour la variable de différence.
• Le score standardisé asymptotique Z s’arrondit strictement à deux décimales (par exemple : Z = -4,01), tout comme la taille d’effet (r = 0,80).
• La probabilité critique (p) s’arrondit à trois décimales. Il ne faut jamais faire précéder la décimale d’un zéro initial si la valeur ne peut théoriquement pas excéder 1 (écrire p = .001 et non p = 0.001). Lorsque la valeur calculée est infinitésimale, il est prescrit de noter p < .001 sans jamais indiquer p = .000.
• L’orientation du test (unilatérale ou bilatérale) doit être expressément consignée dans le corps de la phrase.

11.2 Modèles types de paragraphes de résultats

Afin de guider la rédaction du mémoire de thèse ou du manuscrit de recherche, voici un modèle textuel prêt à l’emploi issu directement des résultats obtenus dans notre essai clinique simulé sur l’anxiété :

« Une analyse préliminaire de normalité menée sur les écarts intra-individuels (différence post − pré) a révélé une déviation statistiquement significative du modèle gaussien (test de Shapiro-Wilk, W = 0,891, p = .011). Compte tenu de cette rupture de distribution et de la nature ordinale discontinue de l’inventaire psychométrique, un test des rangs signés de Wilcoxon pour échantillons appariés a été conduit afin d’évaluer l’impact thérapeutique sur l’anxiété. Conformément aux hypothèses cliniques initiales, une procédure unilatérale orientée vers la réduction des scores a été appliquée.

Les analyses démontrent une diminution statistiquement significative des niveaux d’anxiété suite à l’intervention (V = 3,00, Z = -4,01, p < .001). La médiane des scores d’anxiété à l’inventaire STAI a décru de manière substantielle entre la phase pré-thérapeutique (Mdn = 49,00, IQR = 9,00) et la phase post-thérapeutique (Mdn = 41,00, IQR = 8,00), correspondant à une réduction médiane intra-sujet de -7,00 points (IQR = 5,00). L’amplitude du changement thérapeutique s’avère particulièrement robuste, comme en témoigne une taille d’effet standardisée de Rosenthal majeure (r = .80, IC 95 % [.62, .91]), confirmant que 92 % des participants de la cohorte ont présenté une amélioration symptomatique objective. »

Dans l’éventualité d’un résultat non significatif, le compte rendu textuel adoptera la structure suivante :

« Le test des rangs signés de Wilcoxon apparié n’a mis en évidence aucune modification statistiquement significative des scores entre le temps pré-intervention (Mdn = 34,50, IQR = 6,25) et le temps post-intervention (Mdn = 35,00, IQR = 7,00 ; V = 142,50, Z = -0,38, p = .703, r = .08). L’estimation ponctuelle de la différence médiane intra-sujet demeure nulle (Mdn = 0,00, IQR = 3,00), indiquant l’absence d’efficacité différentielle du protocole à ce stade d’investigation. »

11.3 Conception de tables synthétiques de présentation

La présentation tabulaire synthétique constitue l’autre versant obligatoire de la diffusion des résultats selon l’APA. Les tables conformes aux normes APA 7 se caractérisent par un dépouillement visuel absolu : toute ligne verticale de séparation est rigoureusement proscrite, et seules trois lignes horizontales principales encadrent la table (au sommet, sous l’en-tête de colonnes et à la base de la table avant les notes de bas de page).

Le package gt ou kableExtra permet de générer des tables conformes directement à partir d’un script R Markdown. La structure analytique standard d’un tableau récapitulatif pré-post doit présenter de gauche à droite :

• La variable clinique mesurée.
• Les scores pré-interventionnels : Médiane (Mdn) et Intervalle Interquartile (IQR).
• Les scores post-interventionnels : Médiane (Mdn) et Intervalle Interquartile (IQR).
• La variation intra-sujet : Médiane des différences et IQR de la différence.
• La statistique de Wilcoxon (V ou Z).
• La probabilité critique (p).
• La taille d’effet standardisée (r ou r_c) avec son intervalle de confiance à 95 %.

En note de bas de tableau, il convient de préciser systématiquement la méthode d’ajustement retenue pour le calcul de la statistique (par exemple : « Note. N = 25 paires appariées. Les valeurs p ont été calculées via la distribution exacte conditionnelle prenant en compte les rangs ex æquo selon la méthode de Pratt. * p < .05. *** p < .001. »). Cette transparence tabulaire garantit l’adhésion immédiate des pairs et éditeurs aux conclusions formulées.

12. Erreurs courantes, limites méthodologiques et alternatives avancées

12.1 Pièges analytiques fréquents dans R

L’application routinière du test des rangs signés de Wilcoxon dans R expose le praticien à plusieurs écueils méthodologiques et informatiques discrets mais dévastateurs pour la validité inférentielle :

• Le piège de l’omission du paramètre paired : Comme évoqué précédemment, lancer wilcox.test(pre, post) sans adjoindre paired = TRUE transforme l’analyse en test de Mann-Whitney non apparié. L’analyste perd la corrélation intra-sujet, affaiblit dramatiquement sa puissance statistique et formule une conclusion faussée.
• Le désalignement structurel des données au format large : Si les vecteurs pre et post proviennent de fichiers distincts ayant subi un tri indépendant sans vérification de la concordance ligne à ligne des identifiants de participants, l’appariement est complètement brisé. Les calculs d’écarts d_i soustrairont le score post du sujet 1 au score pré du sujet 15, générant un bruit aléatoire dénué de tout sens clinique.
• La confusion terminologique entre la statistique V et la statistique T : Les chercheurs formés sur des manuels classiques s’attendent à trouver la statistique T (minimum de W⁺ et W^–). Constatant que la valeur de V rapportée par R est très élevée, ils concluent parfois à tort à une divergence de calcul, ignorant que R rapporte toujours systématiquement la somme des rangs positifs W⁺.
• La mauvaise gestion des valeurs manquantes unilatérales : Lors de l’extraction manuelle de deux colonnes contenant des NA, si les longueurs vectorielles diffèrent à la suite d’un filtrage incomplet (par exemple en utilisant na.omit() séparément sur chaque colonne), R rejettera l’analyse avec le message d’erreur : 'x' and 'y' must have the same length.

12.2 Limites inhérentes au test des rangs signés

En dépit de son élégance et de sa robustesse incontestables, le test des rangs signés de Wilcoxon présente des limites théoriques structurelles qu’un statisticien averti ne saurait occulter :

En premier lieu, la transformation en rangs opère une perte inéluctable d’information quantitative. En ramenant les écarts métriques à leur seule position ordinale ordonnée, le test traite de la même manière un écart séparant les rangs 1 et 2 valant 0,1 point et un écart séparant les rangs 24 et 25 valant 40 points. Cette censure métrique peut réduire la sensibilité du test face à des phénomènes subtils où la magnitude absolue du changement véhicule un sens biologique ou psychologique critique.

En second lieu, le test de Wilcoxon est un test univarié strict, incapable d’ajuster pour des covariables continues ou catégorielles. Si l’amplitude du changement pré-post est confondue par des facteurs parasites externes (tels que l’âge des participants, leur niveau de sévérité initiale, le sexe ou l’observance thérapeutique), le modèle non paramétrique élémentaire est incapable d’isoler l’effet thérapeutique net par ajustement multivarié, ouvrant la voie à des biais de confusion non contrôlés.

Enfin, le test est structurellement inopérant face aux plans factoriels complexes impliquant plusieurs facteurs d’interaction (par exemple un plan mixte 2 × 2 confrontant un groupe thérapie active versus un groupe contrôle placebo évalués tous deux en pré et post-test). Vouloir décomposer ce protocole complexe en de multiples tests de Wilcoxon séparés constitue une hérésie méthodologique qui décuple le taux d’erreur de première espèce global (inflation de l’alpha familial) et interdit l’estimation rigoureuse de l’effet d’interaction croisé.

12.3 Paradigmes alternatifs et modèles modernes

Lorsque les données expérimentales violent à la fois les postulats de normalité gaussienne et les restrictions du test univarié de Wilcoxon, l’analyste moderne peut mobiliser des architectures statistiques contemporaines beaucoup plus flexibles :

Les Modèles Linéaires Mixtes (LMM) et les Modèles Mixtes Généralisés (GLMM) représentent l’alternative de référence prédominante. Implémentés sous R via les packages lme4 ou glmmTMB, ces modèles intègrent des effets aléatoires spécifiques par sujet pour capter la corrélation intra-individuelle tout en permettant de modéliser des distributions non gaussiennes de la famille exponentielle (distributions gamma pour les variables positives asymétriques, binomiales négatives pour les données de comptage avec surdispersion). De surcroît, ils gèrent avec une tolérance optimale les données manquantes aléatoires (MAR) sans sacrifier les participants incomplets.

La régression quantile pour mesures répétées offre une approche remarquable lorsque l’objectif scientifique cible spécifiquement l’évolution de quantiles d’intérêt (par exemple, modéliser le 10^e percentile ou le 90^e percentile des scores d’anxiété pour cibler les patients les plus sévèrement atteints). Le package lqmm permet d’estimer ces régressions quantiles conditionnelles avec effets aléatoires intra-sujets sans aucun postulat distributionnel de normalité.

Enfin, dans le champ strictement non paramétrique multidimensionnel, l’école statistique d’Edgar Brunner a développé une méthodologie d’avant-garde : l’analyse non paramétrique des plans longitudinaux à mesures répétées. Cette approche, disponible dans les packages R spécialisés nparLD (Nonparametric Analysis of Longitudinal Data) et rankFD, permet de tester des effets principaux et des interactions factorielles complexes sur des rangs relatifs sans jamais formuler d’hypothèse paramétrique, constituant le prolongement moderne le plus abouti et rigoureux du travail pionnier initié par Frank Wilcoxon en 1945.

Références

• Brunner, E., Domhof, S., & Langer, F. (2002). Nonparametric analysis of longitudinal data in factorial experiments. John Wiley & Sons.
• Cohen, J. (1988). Statistical power analysis for the behavioral sciences (2nd ed.). Lawrence Erlbaum Associates.
• 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
• Kassambara, A. (2023). rstatix: Pipe-Friendly Framework for Basic Statistical Tests (R package version 0.7.2). https://CRAN.R-project.org/package=rstatix
• Kerby, D. S. (2014). The simple difference formula: An approach to teaching nonparametric correlation. Comprehensive Psychology, 3, Article 1. https://doi.org/10.2466/11.IT.3.1
• Pratt, J. W. (1959). Remarks on zeros and ties in the Wilcoxon signed rank procedures. Journal of the American Statistical Association, 54(287), 655–667. https://doi.org/10.1080/01621459.1959.10501526
• R Core Team. (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/
• Rosenthal, R. (1991). Meta-analytic procedures for social research (Rev. ed.). SAGE Publications. https://doi.org/10.4135/9781412984997
• Stevens, S. S. (1946). On the theory of scales of measurement. Science, 103(2684), 677–680. https://doi.org/10.1126/science.103.2684.677
• Wilcoxon, F. (1945). Individual comparisons by ranking methods. Biometrics Bulletin, 1(6), 80–83. https://doi.org/10.2307/3001968