13 - Statistica non parametrica: Wilcoxon, Mann-Whitney

Hai misurato lo stipendio di sette persone di un reparto e di sei di un altro, e uno dei valori è 9,5 mila euro quando tutti gli altri stanno fra 1,8 e 3,8. Con un t-test confronteresti due medie, e quel valore sposta la media del primo di circa un punto. Con un test non parametrico ti chiedi una cosa diversa: se, scegliendo a caso una persona per reparto, la prima guadagna più della seconda più o meno una volta su due. Questi test non assumono che i dati seguano una distribuzione particolare, come la normale: lavorano sull'ordine dei valori, non sulla loro grandezza. Questo articolo, il tredicesimo della serie Probabilità e Statistica, calcola a mano il test di Mann-Whitney, il test di Wilcoxon per dati appaiati, quello di Kruskal-Wallis e le correlazioni di rango, li confronta con SciPy, e misura con simulazioni quando conviene usarli e che cosa testano davvero. Se non conosci il p-value e i test t, parti da lì.

Dai valori ai ranghi: il test U di Mann-Whitney

Il rango di un valore è la sua posizione nell'elenco ordinato di tutti i dati insieme: il più piccolo ha rango 1, il successivo 2, e così via; ai valori uguali si dà la media delle posizioni che occuperebbero. La statistica U di Mann-Whitney (1947) si può leggere in due modi equivalenti: è il numero di coppie (a, b), con a del primo gruppo e b del secondo, in cui a supera b (con mezzo punto per le parità), oppure è la somma dei ranghi del primo gruppo meno il minimo possibile, n₁(n₁ + 1)/2. Il rapporto U / (n₁ n₂) è la probabilità che un valore estratto dal primo gruppo superi uno estratto dal secondo.

Il listato calcola U in entrambi i modi e poi ottiene il p-value esatto senza formule: sotto l'ipotesi nulla (i due gruppi vengono dalla stessa distribuzione) tutte le assegnazioni dei 13 ranghi ai due gruppi sono equiprobabili, e le assegnazioni possibili sono C(13, 6) = 1716, poche abbastanza da elencarle. Il p-value è la quota di quelle che producono un U lontano dal valore atteso almeno quanto quello osservato. Ci sono due modi di definire «altrettanto estremo» e il codice stampa entrambi, perché la distribuzione di U con le parità non è simmetrica. Il test di Welch (un t-test che non assume varianze uguali) serve da confronto.

import itertools
import numpy as np
from scipy import stats

# Stipendi mensili netti (migliaia di euro): reparto A (7 persone, un valore alto), reparto B (6 persone)
A = np.array([1.8, 2.1, 2.1, 2.4, 2.6, 3.0, 9.5])
B = np.array([2.2, 2.5, 2.9, 3.1, 3.4, 3.8])
n1, n2 = len(A), len(B)

# U = numero di coppie (a, b) con a > b, contando 1/2 per ogni parita'
U = sum((a > b) + 0.5 * (a == b) for a in A for b in B)
print(f"U(A) = {U}   su {n1 * n2} coppie   P(A > B) = {U / (n1 * n2):.4f}")

# Stessa cosa dai ranghi: U = R_A - n1 (n1 + 1) / 2, con i ranghi medi per le parita'
rangi = stats.rankdata(np.concatenate([A, B]))
R_A = rangi[:n1].sum()
print(f"somma dei ranghi di A = {R_A}   U da ranghi = {R_A - n1 * (n1 + 1) / 2}")

# p esatto per permutazione: provo tutte le C(13,6) = 1716 assegnazioni dei ranghi ai due gruppi
tutte = np.array([rangi[list(idx)].sum() - n1 * (n1 + 1) / 2
                  for idx in itertools.combinations(range(n1 + n2), n1)])
media_U = n1 * n2 / 2
p_distanza = np.mean(np.abs(tutte - media_U) >= abs(U - media_U) - 1e-12)   # "altrettanto lontano dalla media"
p_code = 2 * min(np.mean(tutte <= U + 1e-12), np.mean(tutte >= U - 1e-12))  # il doppio della coda piu' piccola
print(f"p esatto: {int(p_distanza * len(tutte))}/{len(tutte)} = {p_distanza:.4f} (distanza dalla media)   {p_code:.4f} (doppio della coda)")

asin = stats.mannwhitneyu(A, B, alternative="two-sided", method="asymptotic")
perm = stats.mannwhitneyu(A, B, alternative="two-sided", method=stats.PermutationMethod(n_resamples=np.inf))
print(f"scipy asintotico (correzione per le parita') p = {asin.pvalue:.4f}   scipy permutazione esatta p = {perm.pvalue:.4f}")
t = stats.ttest_ind(A, B, equal_var=False)
print(f"t di Welch p = {t.pvalue:.4f}   medie {A.mean():.2f} contro {B.mean():.2f}   mediane {np.median(A):.2f} contro {np.median(B):.2f}")

# Output:
# U(A) = 12.0   su 42 coppie   P(A > B) = 0.2857
# somma dei ranghi di A = 40.0   U da ranghi = 12.0
# p esatto: 377/1716 = 0.2197 (distanza dalla media)   0.2214 (doppio della coda)
# scipy asintotico (correzione per le parita') p = 0.2240   scipy permutazione esatta p = 0.2214
# t di Welch p = 0.7357   medie 3.36 contro 2.98   mediane 2.40 contro 3.00

U vale 12 su 42 coppie: nel 28,6% dei confronti il reparto A guadagna più del B. Le due vie di calcolo danno lo stesso U. Il p-value esatto è 0,2197 con la regola «distanza dalla media» e 0,2214 con la regola «doppio della coda più piccola», quest'ultima usata da SciPy; l'approssimazione normale con correzione per le parità dà 0,2240. Tre cifre diverse e nessuna conclusione diversa: con 13 persone non c'è evidenza di una differenza. Il t di Welch dà 0,7357, perché la media del reparto A (3,36) supera quella del B (2,98) solo per effetto del valore 9,5, mentre la mediana di A (2,40) è inferiore a quella di B (3,00). Il caso mostra che l'esito dei due test non è lo stesso quando un valore estremo guida la media, ma non che il test sui ranghi sia «più potente»: qui nessuno dei due rifiuta. Il metodo exact di mannwhitneyu non corregge le parità, come dichiara la documentazione, e per questo il listato usa PermutationMethod.

Che cosa testa davvero Mann-Whitney

Nei manuali si legge spesso che Mann-Whitney «confronta le mediane». È vero solo sotto un'ipotesi in più: che le due distribuzioni abbiano la stessa forma e differiscano per uno spostamento. In generale il test rileva se P(X > Y) si discosta da 1/2, cioè se una variabile tende a superare l'altra, e due popolazioni con mediane identiche possono violare questa condizione. La simulazione seguente lo mostra con X esponenziale e Y normale, entrambe con mediana ln 2 = 0,6931. La probabilità esatta P(X > Y), calcolata per integrazione numerica, vale 0,5218, non 0,5.

import numpy as np
from scipy import stats, integrate

rng = np.random.default_rng(2027)
alfa = 0.05

# Falsa credenza: "Mann-Whitney confronta le mediane". Due popolazioni con la STESSA mediana:
# X ~ esponenziale (mediana ln 2), Y ~ normale con media ln 2 e dev. standard 0.3
med = np.log(2)
X, Y = stats.expon(), stats.norm(med, 0.3)
p_sup = integrate.quad(lambda y: Y.pdf(y) * X.sf(y), -np.inf, np.inf)[0]
print(f"mediane: {X.median():.4f} e {Y.median():.4f}   P(X > Y) esatta = {p_sup:.4f}")
for m in (50, 200, 1000):
    x = X.rvs((2000, m), random_state=rng); y = Y.rvs((2000, m), random_state=rng)
    rif = np.mean(stats.mannwhitneyu(x, y, axis=1, method="asymptotic").pvalue < alfa)
    print(f"n = {m:4d} per gruppo: rifiuta H0 nel {rif:.3f} dei casi")

# Output:
# mediane: 0.6931 e 0.6931   P(X > Y) esatta = 0.5218
# n =   50 per gruppo: rifiuta H0 nel 0.081 dei casi
# n =  200 per gruppo: rifiuta H0 nel 0.122 dei casi
# n = 1000 per gruppo: rifiuta H0 nel 0.386 dei casi

Le mediane coincidono e la probabilità di superamento è 0,5218, ma il test rifiuta l'ipotesi nulla con frequenza crescente al crescere del campione: 8,1% con 50 casi per gruppo, 12,2% con 200 e 38,6% con 1000 (errore di simulazione circa 0,01). Non è un difetto del test, risponde a una domanda diversa da quella che si credeva di porre. Se la domanda è «le mediane sono uguali?», Mann-Whitney non è lo strumento giusto senza l'ipotesi di stessa forma; se è «una popolazione tende a dare valori più alti dell'altra?», lo è. Per questo conviene riportare accanto al p-value la stima di P(X > Y), cioè U / (n₁ n₂), che dice anche di quanto.

Quanto si perde rinunciando al t-test, e quando si guadagna

Se i dati sono davvero normali il t-test è il più potente tra quelli possibili; usare i ranghi costa poco. L'efficienza relativa asintotica (ARE), cioè il rapporto fra le dimensioni campionarie di due test che raggiungono la stessa potenza, vale 3/π ≈ 0,955 per il test di Wilcoxon contro il t-test sotto normalità: perdi circa il 4,5% di efficienza. Hodges e Lehmann (1956) hanno dimostrato che, per qualsiasi distribuzione, l'ARE non scende sotto 0,864. Il guadagno invece può essere grande quando le code sono pesanti, perché pochi valori estremi gonfiano la varianza usata dal t-test, mentre i ranghi li limitano.

Il listato lo verifica con 4000 simulazioni per riga e con 30 osservazioni per gruppo: stessa varianza (uguale a 1) e stesso spostamento di 0,6 deviazioni standard fra i due gruppi, cambia soltanto la forma. Si conta la frazione di simulazioni in cui ciascun test rifiuta l'ipotesi nulla al livello 5%, cioè la potenza.

import numpy as np
from scipy import stats, integrate

rng = np.random.default_rng(2026)
n, R, alfa = 30, 4000, 0.05

def potenza(estrai, spostamento):
    x = estrai((R, n))
    y = estrai((R, n)) + spostamento
    p_t = stats.ttest_ind(x, y, axis=1, equal_var=False).pvalue
    p_u = stats.mannwhitneyu(x, y, axis=1, method="asymptotic").pvalue
    return np.mean(p_t < alfa), np.mean(p_u < alfa)

casi = [("normale",           lambda s: rng.standard_normal(s),          0.6),
        ("Laplace (code pesanti)", lambda s: rng.laplace(size=s) / np.sqrt(2), 0.6),   # varianza 1
        ("t con 3 gdl",       lambda s: rng.standard_t(3, s) / np.sqrt(3),  0.6)]      # varianza 1
print("distribuzione            potenza t   potenza Mann-Whitney")
for nome, f, d in casi:
    a, b = potenza(f, d)
    print(f"{nome:24s} {a:8.3f}   {b:8.3f}    (errore ~0.008)")

# Output:
# distribuzione            potenza t   potenza Mann-Whitney
# normale                     0.628      0.594    (errore ~0.008)
# Laplace (code pesanti)      0.636      0.749    (errore ~0.008)
# t con 3 gdl                 0.696      0.856    (errore ~0.008)

Con dati normali la potenza è 0,628 per il t-test e 0,594 per Mann-Whitney: una perdita di 3-4 punti percentuali, coerente con 0,955. Con la distribuzione di Laplace (code più pesanti) i ruoli si invertono, 0,636 contro 0,749, e con la t di Student a 3 gradi di libertà, che ha code ancora più pesanti, la distanza arriva a 0,696 contro 0,856. L'errore delle simulazioni è circa 0,008. La regola pratica che ne esce è asimmetrica: il costo di usare i ranghi quando non servivano è piccolo, il costo di usare il t-test su code pesanti arriva a 16 punti percentuali di potenza (0,696 contro 0,856). Con pochi dati, inoltre, la normalità non si riesce a verificare bene, e il test sui ranghi è una scelta prudente.

Dati appaiati: il test dei ranghi con segno di Wilcoxon

Quando le misure sono accoppiate, come lo stesso gruppo di pagine prima e dopo un intervento, si lavora sulle differenze d = prima − dopo. Il test di Wilcoxon (1945) ordina i valori assoluti |d| e somma i loro ranghi separatamente per le differenze positive (W⁺) e per quelle negative (W⁻); la loro somma è sempre n(n + 1)/2. Sotto l'ipotesi nulla, ogni differenza ha segno positivo o negativo con probabilità 1/2, quindi con 8 coppie le assegnazioni possibili sono 2⁸ = 256 e si possono elencare. Il test presuppone che la distribuzione delle differenze sia simmetrica attorno a zero sotto l'ipotesi nulla.

import itertools
import numpy as np
from scipy import stats

# Tempo di caricamento di 8 pagine (secondi) prima e dopo un intervento
prima = np.array([4.1, 3.8, 5.0, 2.9, 6.2, 3.5, 4.7, 5.5])
dopo  = np.array([3.5, 3.9, 4.1, 2.6, 5.0, 3.0, 4.5, 4.2])
d = prima - dopo                       # differenze appaiate: il test lavora su queste
print("differenze:", np.round(d, 2))

rangi = stats.rankdata(np.abs(d))      # ranghi dei valori assoluti
W_piu = rangi[d > 0].sum()             # somma dei ranghi delle differenze positive
W_meno = rangi[d < 0].sum()
print("ranghi |d| :", rangi, "  W+ =", W_piu, " W- =", W_meno, " (W+ + W- =", rangi.sum(), ")")

# p esatto: sotto H0 ogni differenza ha segno + o - con prob. 1/2, quindi 2^8 = 256 assegnazioni
somme = np.array([sum(r for r, s in zip(rangi, segni) if s)
                  for segni in itertools.product([False, True], repeat=len(d))])
W = min(W_piu, W_meno)
p_esatto = np.mean(somme <= W) * 2     # la distribuzione e' simmetrica: raddoppio la coda
print(f"p esatto per enumerazione = {int(np.sum(somme <= W))}/{len(somme)} x 2 = {p_esatto:.4f}")

ris = stats.wilcoxon(prima, dopo, method="exact")
print(f"scipy.stats.wilcoxon: statistica = {ris.statistic}, p = {ris.pvalue:.4f}")
print(f"t appaiato: p = {stats.ttest_rel(prima, dopo).pvalue:.4f}")

# Con una differenza nulla la coppia viene scartata (zero_method='wilcox'): restano 7 differenze
dopo2 = dopo.copy(); dopo2[1] = prima[1]
r2 = stats.wilcoxon(prima, dopo2)
print(f"con una differenza nulla: statistica = {r2.statistic}, p = {r2.pvalue:.4f}")

# Output:
# differenze: [ 0.6 -0.1  0.9  0.3  1.2  0.5  0.2  1.3]
# ranghi |d| : [5. 1. 6. 3. 7. 4. 2. 8.]   W+ = 35.0  W- = 1.0  (W+ + W- = 36.0 )
# p esatto per enumerazione = 2/256 x 2 = 0.0156
# scipy.stats.wilcoxon: statistica = 1.0, p = 0.0156
# t appaiato: p = 0.0096
# con una differenza nulla: statistica = 0.0, p = 0.0156

Delle otto differenze, solo una è negativa (−0,1) e ha rango 1, quindi W⁻ = 1 e W⁺ = 35. Solo 2 delle 256 assegnazioni di segni producono una somma di ranghi negativi pari o inferiore a 1, e il p-value bilaterale è 2 · 2/256 = 0,0156, identico a quello di scipy.stats.wilcoxon. Il t appaiato dà 0,0096: stessa conclusione, p-value diverso. Una trappola vera: due differenze che a mano sono uguali, come 4,1 − 3,6 e 3,5 − 3,0 (entrambe 0,5), in virgola mobile valgono 0,49999999999999956 e 0,5, e i due valori ricevono ranghi diversi invece di un rango medio comune. La documentazione di SciPy avverte che conviene calcolare la differenza, arrotondarla e passarla al test come unico argomento. Le coppie con differenza nulla vengono scartate (zero_method='wilcox', il valore predefinito): l'ultima riga lo mostra con 7 differenze utili.

Più di due gruppi: Kruskal-Wallis e confronti a coppie

La generalizzazione di Mann-Whitney a k gruppi è il test di Kruskal-Wallis (1952), l'analogo non parametrico dell'ANOVA a una via vista nell'articolo su ANOVA e confronti multipli. Si ordinano tutti gli N dati insieme e si calcola H = 12 / (N(N + 1)) · Σ R_i² / n_i − 3(N + 1), dove R_i è la somma dei ranghi del gruppo i e n_i la sua numerosità; per campioni non troppo piccoli H si confronta con una distribuzione chi-quadro con k − 1 gradi di libertà. Se il test rifiuta, dice che almeno un gruppo tende a dare valori diversi, non quale: per saperlo si fanno confronti a coppie, e qui serve correggere i p-value per i confronti multipli, altrimenti la probabilità di almeno un falso positivo supera il livello dichiarato. La correzione di Holm li ordina dal più piccolo, moltiplica il primo per il numero di confronti m, il secondo per m − 1 e così via, e rende il risultato non decrescente.

import numpy as np
from scipy import stats

# Tempi di risposta (ms) di tre versioni di un servizio, 5 misure ciascuna
g = [np.array([212, 230, 245, 251, 300.0]),
     np.array([260, 275, 290, 310, 335.0]),
     np.array([280, 320, 340, 360, 410.0])]
tutti = np.concatenate(g)
N = len(tutti)
rangi = stats.rankdata(tutti)
R = [rangi[:5], rangi[5:10], rangi[10:]]
H = 12 / (N * (N + 1)) * sum(r.sum() ** 2 / len(r) for r in R) - 3 * (N + 1)   # senza parita'
p_chi2 = stats.chi2.sf(H, df=len(g) - 1)
kw = stats.kruskal(*g)
print(f"somme dei ranghi: {[float(r.sum()) for r in R]}   H a mano = {H:.4f}  p (chi2, 2 gdl) = {p_chi2:.4f}")
print(f"scipy.stats.kruskal: H = {kw.statistic:.4f}  p = {kw.pvalue:.4f}")

# Confronti a coppie con la correzione di Holm: ordino i p, moltiplico per m, m-1, ..., 1 e rendo monotono
coppie = [(0, 1), (0, 2), (1, 2)]
p = np.array([stats.mannwhitneyu(g[i], g[j], method="exact").pvalue for i, j in coppie])
ordine = np.argsort(p)
m = len(p)
holm = np.empty(m); corrente = 0.0
for rango, k in enumerate(ordine):
    corrente = max(corrente, min(1.0, (m - rango) * p[k])); holm[k] = corrente
for (i, j), pg, ph in zip(coppie, p, holm):
    print(f"versione {i + 1} contro {j + 1}: p = {pg:.4f}  p corretto (Holm) = {ph:.4f}")

# Output:
# somme dei ranghi: [19.0, 41.0, 60.0]   H a mano = 8.4200  p (chi2, 2 gdl) = 0.0148
# scipy.stats.kruskal: H = 8.4200  p = 0.0148
# versione 1 contro 2: p = 0.0556  p corretto (Holm) = 0.1111
# versione 1 contro 3: p = 0.0159  p corretto (Holm) = 0.0476
# versione 2 contro 3: p = 0.0952  p corretto (Holm) = 0.1111

I ranghi dei tre gruppi sommano 19, 41 e 60, e H vale 8,42 sia con la formula a mano sia con SciPy (che applica la correzione per le parità, assente qui), con p = 0,0148 sotto l'approssimazione chi-quadro. Nei confronti a coppie solo versione 1 contro versione 3 resta sotto 0,05 dopo la correzione di Holm (0,0159 diventa 0,0476), mentre 1 contro 2 passa da 0,0556 a 0,1111. Il test globale e i confronti a coppie non devono per forza accordarsi: Kruskal-Wallis usa tutti i dati insieme, i confronti a coppie ne usano una parte e pagano la correzione. Esiste anche il test post-hoc di Dunn, che usa i ranghi globali, ma non è nella libreria standard di SciPy e non è eseguito qui.

Correlazioni di rango: Spearman e Kendall

La correlazione di Pearson misura l'associazione lineare ed è sensibile ai valori anomali. Le correlazioni di rango misurano l'associazione monotòna, cioè quanto una variabile cresca insieme all'altra senza chiedere che lo faccia lungo una retta. La ρ di Spearman è la correlazione di Pearson calcolata sui ranghi. La τ di Kendall (1938) è (C − D) / (numero di coppie), dove C conta le coppie concordanti (le due variabili si muovono nello stesso verso) e D quelle discordanti.

import numpy as np
from scipy import stats

# Correlazioni di rango su un legame monotono ma non lineare, con un punto anomalo
x = np.arange(1.0, 11.0)
y = np.exp(0.5 * x)                          # crescente, curvo
y_out = y.copy(); y_out[-1] = 1.0            # l'ultimo punto crolla
def kendall(a, b):                           # (concordanti - discordanti) / coppie, senza parita'
    c = d = 0
    for i in range(len(a)):
        for j in range(i + 1, len(a)):
            s = (a[i] - a[j]) * (b[i] - b[j]); c += s > 0; d += s < 0
    return (c - d) / (len(a) * (len(a) - 1) / 2)
for nome, yy in (("curva esponenziale", y), ("con punto anomalo", y_out)):
    sp = np.corrcoef(stats.rankdata(x), stats.rankdata(yy))[0, 1]
    print(f"{nome:19s} Pearson {stats.pearsonr(x, yy)[0]:.4f}  Spearman {sp:.4f} (scipy {stats.spearmanr(x, yy)[0]:.4f})  "
          f"Kendall {kendall(x, yy):.4f} (scipy {stats.kendalltau(x, yy)[0]:.4f})")

# Output:
# curva esponenziale  Pearson 0.8654  Spearman 1.0000 (scipy 1.0000)  Kendall 1.0000 (scipy 1.0000)
# con punto anomalo   Pearson 0.5920  Spearman 0.4545 (scipy 0.4545)  Kendall 0.6000 (scipy 0.6000)

Su una curva esponenziale crescente Pearson è 0,8654 perché la relazione non è una retta, mentre Spearman e Kendall valgono 1,0000: l'ordine è perfettamente rispettato. Se l'ultimo punto crolla a 1, Pearson scende a 0,5920, Spearman a 0,4545 e Kendall a 0,6000: tutti e tre risentono dell'anomalia, ma la risentono in misura diversa, e nessun coefficiente è «robusto» in assoluto. Il calcolo a mano coincide con SciPy fino alla quarta cifra. Quale scegliere? Per i dati senza parità e con poche osservazioni τ ha un'interpretazione diretta come differenza fra la probabilità di concordanza e di discordanza di una coppia; ρ è più usata per tradizione.

Esercizio svolto e prova da fare da solo

Cinque pagine, differenze di tempo (secondi) +1,5; −0,5; +2,0; +0,8; +3,0. I valori assoluti ordinati sono 0,5 / 0,8 / 1,5 / 2,0 / 3,0 con ranghi 1-5, quindi W⁻ = 1 (il rango di 0,5). Delle 2⁵ = 32 assegnazioni dei segni, due hanno somma negativa ≤ 1 (nessun segno negativo, oppure solo il rango 1), quindi il p-value bilaterale esatto è 2 · 2/32 = 0,125. Con cinque coppie il valore più piccolo possibile è 2/32 = 0,0625: nessun risultato può rifiutare al 5%, per quanto netta sia la differenza. Il test è stato verificato con scipy.stats.wilcoxon(..., method='exact'), che restituisce 0,125.

Che cosa sai fare adesso

Sai calcolare i ranghi e la statistica U, ottenere il p-value esatto per enumerazione e leggere U come probabilità che un gruppo superi l'altro. Sai che Mann-Whitney testa il superamento stocastico e non le mediane, che costa circa il 4,5% di efficienza sotto normalità e può guadagnare molto con code pesanti, e sai applicare Wilcoxon, Kruskal-Wallis con la correzione di Holm e le correlazioni di rango. Il prossimo articolo, il bootstrap, costruisce intervalli e test senza ipotesi sulla distribuzione, ricampionando i dati. Il chi-quadro per le tabelle di frequenze è nell'articolo 9; la catena di Markov dell'articolo 12 precede questo.

Nella stessa serie trovi anche Analisi dati esplorativa (EDA) con Python.