08 - ANOVA e confronto multipli

Hai quattro varianti di una pagina web e vuoi sapere se i tempi di caricamento medi sono davvero diversi. La strada ovvia è confrontare le varianti a coppie con un t-test: sono sei confronti, ognuno al livello del 5%. Il problema è che, anche se le quattro medie sono identiche, la probabilità che almeno uno dei sei confronti risulti «significativo» è circa il 20%, non il 5% (la simulazione qui sotto dà 19,9%). L'analisi della varianza (ANOVA) risolve la prima metà del problema con un solo test: chiede se le medie sono tutte uguali, con livello 5% garantito.

La seconda metà, capire quali differiscono, richiede un confronto a coppie corretto per i confronti multipli. Questo articolo, l'ottavo della serie Probabilità e Statistica, copre entrambe, per il caso a un solo fattore, e ti mostra con simulazioni dove il test si rompe. Servono il t-test e il p-value dell'articolo precedente.

L'idea: due stime della stessa varianza

Immagina tre fertilizzanti, con cinque parcelle ciascuno. Le rese cambiano per due ragioni: perché i fertilizzanti sono diversi (variabilità fra i gruppi) e perché anche parcelle con lo stesso fertilizzante danno rese diverse (variabilità dentro i gruppi, il rumore). L'ANOVA misura le due variabilità con la stessa unità e ne fa il rapporto. Sia x̄_i la media del gruppo i con n_i osservazioni, x̄ la media di tutte le N osservazioni e k il numero di gruppi.

La somma dei quadrati fra i gruppi è SS_B = Σn_i(x̄_i − x̄)², la somma dei quadrati dentro i gruppi è SS_W = Σ(x_ij − x̄_i)², e la somma dei quadrati totale SS_T = Σ(x_ij − x̄)² è sempre uguale a SS_B + SS_W: un'identità algebrica, non un'approssimazione.

Dividendo ciascuna somma per i suoi gradi di libertà (k − 1 per SS_B, N − k per SS_W) si ottengono le medie dei quadrati MS_B e MS_W. MS_W stima la varianza σ² dei dati qualunque sia la verità sulle medie; MS_B la stima solo se le medie sono uguali, e se non lo sono tende a essere più grande. Quindi F = MS_B/MS_W vale circa 1 sotto H₀ e cresce con le differenze fra le medie.

Se i dati di ogni gruppo sono Normali con la stessa varianza e indipendenti, la distribuzione esatta di F sotto H₀ è una F di Fisher con (k − 1, N − k) gradi di libertà; la procedura e le sue assunzioni sono descritte nel manuale di statistica del NIST, sezione 7.4.3.

Il programma sviluppa un esercizio piccolo a mano e lo verifica: tre fertilizzanti con cinque parcelle, medie 23, 28 e 24, scelti perché i conti vengano tondi. La media generale è 25, SS_B = 5·(4 + 9 + 1) = 70, SS_W = 10 + 10 + 10 = 30, quindi MS_B = 70/2 = 35, MS_W = 30/12 = 2,5 e F = 14. Il programma ripete il calcolo, usa stats.f.sf per il p-value e confronta con scipy.stats.f_oneway; calcola anche tre indici di dimensione dell'effetto.

# Esercizio svolto: ANOVA a una via su tre gruppi da cinque misure, a mano e con scipy
import numpy as np
from scipy import stats

gruppi = [np.array([23, 25, 21, 24, 22.0]),    # media 23
          np.array([27, 29, 26, 28, 30.0]),    # media 28
          np.array([24, 22, 26, 25, 23.0])]    # media 24
k, N = len(gruppi), sum(len(g) for g in gruppi)
tutti = np.concatenate(gruppi)
gran_media = tutti.mean()

SS_B = sum(len(g) * (g.mean() - gran_media) ** 2 for g in gruppi)     # fra i gruppi
SS_W = sum(((g - g.mean()) ** 2).sum() for g in gruppi)               # dentro i gruppi
SS_T = ((tutti - gran_media) ** 2).sum()                              # totale
MS_B, MS_W = SS_B / (k - 1), SS_W / (N - k)
F = MS_B / MS_W
p = stats.f.sf(F, k - 1, N - k)                                       # sf = 1 - cdf, senza perdita di precisione

print(f"medie: {[float(g.mean()) for g in gruppi]}, gran media {float(gran_media)}")
print(f"SS_B = {SS_B:.1f}, SS_W = {SS_W:.1f}, SS_T = {SS_T:.1f}, SS_B + SS_W = {SS_B + SS_W:.1f}")
print(f"MS_B = {MS_B:.1f}, MS_W = {MS_W:.1f}, F = {F:.1f}, gradi di liberta' ({k - 1}, {N - k}), p = {p:.5f}")
r = stats.f_oneway(*gruppi)
print(f"scipy.stats.f_oneway: F = {r.statistic:.1f}, p = {r.pvalue:.5f}")
eta2 = SS_B / SS_T
omega2 = (SS_B - (k - 1) * MS_W) / (SS_T + MS_W)
print(f"eta^2 = {eta2:.4f}, omega^2 = {omega2:.4f}, f di Cohen = {np.sqrt(eta2 / (1 - eta2)):.4f}")

# Output:
# medie: [23.0, 28.0, 24.0], gran media 25.0
# SS_B = 70.0, SS_W = 30.0, SS_T = 100.0, SS_B + SS_W = 100.0
# MS_B = 35.0, MS_W = 2.5, F = 14.0, gradi di liberta' (2, 12), p = 0.00073
# scipy.stats.f_oneway: F = 14.0, p = 0.00073
# eta^2 = 0.7000, omega^2 = 0.6341, f di Cohen = 1.5275

I numeri coincidono con quelli a mano: SS_B = 70, SS_W = 30, SS_T = 100 = 70 + 30, F = 14 con (2, 12) gradi di libertà e p = 0,00073, identico a f_oneway. Un p-value così basso dice che differenze fra le medie come queste sono molto improbabili se i fertilizzanti fossero equivalenti, ma non dice quale sia diverso né quanto.

Gli indici rispondono alla seconda parte: η² = SS_B/SS_T = 0,70 è la quota della variabilità totale spiegata dal gruppo; ω² = (SS_B − (k − 1)MS_W)/(SS_T + MS_W) = 0,6341 è una versione meno gonfiata (η² tende a sovrastimare con pochi dati); f = √(η²/(1 − η²)) = 1,5275 è l'effetto espresso come scarto tipo delle medie rispetto allo scarto tipo dei dati.

Livello e potenza del test F

Che F segua davvero la legge F sotto H₀ si controlla simulando: 200.000 esperimenti con quattro gruppi da 30 osservazioni Normali con la stessa media. Sotto H₁, la statistica segue una F non centrale con parametro di non centralità λ = Σn_i(μ_i − μ̄)²/σ², e scipy.stats.ncf ne dà le probabilità: la documentazione di SciPy definisce la distribuzione. Le medie vere (10, 10,7, 10,7, 11,4 con σ = 2) sono scelte per ottenere un effetto «medio» secondo Cohen, f ≈ 0,25. Il programma confronta potenza esatta e simulata, e cerca quanti dati per gruppo servono per arrivare all'80%.

# Sotto H0 la statistica F ha la legge F(k-1, N-k); sotto H1 una F non centrale: livello e potenza
import numpy as np
from scipy import stats

def statistica_F(X):
    """X ha forma (R, k, n): R ripetizioni, k gruppi, n osservazioni per gruppo."""
    R, k, n = X.shape
    medie = X.mean(axis=2)
    SS_B = n * ((medie - medie.mean(axis=1, keepdims=True)) ** 2).sum(axis=1)
    SS_W = ((X - medie[:, :, None]) ** 2).sum(axis=(1, 2))
    return (SS_B / (k - 1)) / (SS_W / (k * (n - 1)))

k, n, sigma, R = 4, 30, 2.0, 200_000
N, soglia = k * n, stats.f.ppf(0.95, k - 1, k * n - k)
rng = np.random.default_rng(42)

F0 = statistica_F(rng.normal(10.0, sigma, size=(R, k, n)))             # H0: quattro medie uguali
print(f"H0: soglia F(0.95; {k - 1}, {N - k}) = {soglia:.4f}; rifiuti {(F0 > soglia).mean():.4f} (atteso 0.05)")
for q in (0.50, 0.90, 0.99):
    sotto = (F0 <= stats.f.ppf(q, k - 1, N - k)).mean()
    print(f"    quota di F sotto il quantile {q:.2f} della legge F: {sotto:.4f} +/- {np.sqrt(q * (1 - q) / R):.4f}")

mu = np.array([10.0, 10.7, 10.7, 11.4])                                  # H1: quattro medie diverse
lam = n * ((mu - mu.mean()) ** 2).sum() / sigma ** 2                     # parametro di non centralita'
F1 = statistica_F(rng.normal(mu[None, :, None], sigma, size=(R, k, n)))
esatta = stats.ncf.sf(soglia, k - 1, N - k, lam)
sim = (F1 > soglia).mean()
print(f"H1: lambda = {lam:.3f}, f di Cohen = {np.sqrt(lam / N):.4f}")
print(f"potenza esatta (F non centrale) = {esatta:.4f}; simulata = {sim:.4f} +/- {np.sqrt(sim * (1 - sim) / R):.4f}")

def potenza(n_gruppo, f):
    """Potenza esatta di F con k gruppi da n_gruppo osservazioni e effetto f di Cohen."""
    Ntot = k * n_gruppo
    return stats.ncf.sf(stats.f.ppf(0.95, k - 1, Ntot - k), k - 1, Ntot - k, Ntot * f ** 2)

print(f"f = 0.25, {n} per gruppo: potenza esatta = {potenza(n, 0.25):.4f}")
n80 = next(m for m in range(3, 500) if potenza(m, 0.25) >= 0.80)
print(f"f = 0.25: primo n per gruppo con potenza >= 0.80 = {n80} (potenza {potenza(n80, 0.25):.4f})")

# Output:
# H0: soglia F(0.95; 3, 116) = 2.6828; rifiuti 0.0498 (atteso 0.05)
#     quota di F sotto il quantile 0.50 della legge F: 0.5003 +/- 0.0011
#     quota di F sotto il quantile 0.90 della legge F: 0.9002 +/- 0.0007
#     quota di F sotto il quantile 0.99 della legge F: 0.9901 +/- 0.0002
# H1: lambda = 7.350, f di Cohen = 0.2475
# potenza esatta (F non centrale) = 0.5968; simulata = 0.5952 +/- 0.0011
# f = 0.25, 30 per gruppo: potenza esatta = 0.6065
# f = 0.25: primo n per gruppo con potenza >= 0.80 = 45 (potenza 0.8040)

Sotto H₀ la soglia è 2,6828, rifiuti il 4,98% delle volte, e le frazioni di F sotto i quantili 50%, 90% e 99% della legge F sono 0,5003, 0,9002 e 0,9901, tutte entro l'errore standard stampato: la distribuzione teorica regge.

Sotto H₁ il parametro di non centralità è 7,350 (f = 0,2475) e la potenza esatta è 0,5968 contro 0,5952 simulata, con errore standard 0,0011. Il risultato che pesa di più per chi progetta l'esperimento è l'ultimo: con 30 osservazioni per gruppo e un effetto medio la potenza è 0,6065, cioè sbaglierai a non vedere l'effetto quattro volte su dieci, e per l'80% servono 45 osservazioni per gruppo.

Quali medie differiscono: Bonferroni e Tukey

Il test F risponde «almeno una media differisce». Per sapere quali, si confrontano le coppie, controllando l'errore familywise (FWER): la probabilità di avere almeno un falso positivo fra tutti i confronti. Con k = 4 ci sono k(k − 1)/2 = 6 coppie. Bonferroni usa per ciascuna il livello 0,05/6: semplice e valido con qualunque test, ma conservativo.

Tukey HSD usa invece la distribuzione del range studentizzato (la differenza massima fra medie in unità dell'errore standard) e confronta ogni coppia con una soglia comune. Il manuale NIST (sezione 7.4.7.1) precisa che con numerosità uguali il coefficiente di confidenza dell'intero insieme è esattamente 1 − α, e che con numerosità diverse il metodo è conservativo.

Il programma implementa le tre soglie sulla differenza fra medie, usando per tutte la stessa stima comune MS_W (così il confronto è omogeneo), e controlla il calcolo di Tukey confrontando i p-value con scipy.stats.tukey_hsd su un campione. Poi simula 100.000 esperimenti da quattro gruppi di 25: prima con tutte le medie uguali (per misurare l'errore familywise), poi con il quarto gruppo spostato di 1,5 (0,75 scarti tipo) per misurare la probabilità di scoprire la coppia (1, 4).

# Confronti a coppie fra k = 4 gruppi: errore familywise e potenza di t senza correzione, Bonferroni, Tukey
import numpy as np
from scipy import stats

k, n, sigma, R = 4, 25, 2.0, 100_000
N, m = k * n, k * (k - 1) // 2                     # m = 6 confronti
rng = np.random.default_rng(50)
i, j = np.triu_indices(k, 1)                       # le 6 coppie (i, j)

t_unc = stats.t.ppf(1 - 0.05 / 2, N - k)           # soglia t senza correzione (alfa = 0.05 per confronto)
t_bon = stats.t.ppf(1 - 0.05 / (2 * m), N - k)     # Bonferroni: alfa / 6 per confronto
q_tuk = stats.studentized_range.ppf(0.95, k, N - k)
# Soglie sulla differenza delle medie; MS_W e' la stima comune della varianza
def rifiuti(medie_vere):
    X = rng.normal(np.asarray(medie_vere)[None, :, None], sigma, size=(R, k, n))
    medie = X.mean(axis=2)
    MS_W = ((X - medie[:, :, None]) ** 2).sum(axis=(1, 2)) / (N - k)
    diff = np.abs(medie[:, i] - medie[:, j])                       # forma (R, 6)
    se = np.sqrt(MS_W * 2 / n)[:, None]
    return diff > t_unc * se, diff > t_bon * se, diff > q_tuk * se / np.sqrt(2)

# Controllo della implementazione di Tukey con scipy su un campione
x = rng.normal(0, sigma, size=(k, n))
mio = np.abs(x.mean(axis=1)[i] - x.mean(axis=1)[j]) / np.sqrt(((x - x.mean(axis=1, keepdims=True)) ** 2).sum() / (N - k) / n)
p_mio = stats.studentized_range.sf(mio, k, N - k)
print("controllo p di Tukey (mio, scipy):", np.round(p_mio, 5), np.round(stats.tukey_hsd(*x).pvalue[i, j], 5))

print(f"soglie: t = {t_unc:.3f}, t Bonferroni = {t_bon:.3f}, q/sqrt(2) di Tukey = {q_tuk / np.sqrt(2):.3f}")
unc, bon, tuk = rifiuti([10, 10, 10, 10])                          # H0: tutte le medie uguali
print("H0 vera, almeno un rifiuto fra i 6 confronti (errore familywise):")
for nome, r in (("senza correzione", unc), ("Bonferroni", bon), ("Tukey", tuk)):
    print(f"  {nome:17s} {r.any(axis=1).mean():.4f}")
print(f"  riferimento: 1 - 0.95^6 = {1 - 0.95 ** m:.4f} varrebbe per 6 confronti indipendenti; Tukey vale 0.05 per costruzione (gruppi bilanciati)")

unc, bon, tuk = rifiuti([10, 10, 10, 11.5])                        # H1: solo il gruppo 4 e' diverso (1,5 = 0,75 sigma)
coppia = [(a, b) for a, b in zip(i, j)].index((0, 3))
print("H1: il gruppo 4 ha media +1,5; frazione di campioni in cui la coppia (1, 4) risulta significativa:")
for nome, r in (("senza correzione", unc), ("Bonferroni", bon), ("Tukey", tuk)):
    print(f"  {nome:17s} {r[:, coppia].mean():.4f}")

# Output:
# controllo p di Tukey (mio, scipy): [0.33901 0.15543 0.28332 0.97242 0.99951 0.98823] [0.33901 0.15543 0.28332 0.97242 0.99951 0.98823]
# soglie: t = 1.985, t Bonferroni = 2.694, q/sqrt(2) di Tukey = 2.615
# H0 vera, almeno un rifiuto fra i 6 confronti (errore familywise):
#   senza correzione  0.1992
#   Bonferroni        0.0405
#   Tukey             0.0495
#   riferimento: 1 - 0.95^6 = 0.2649 varrebbe per 6 confronti indipendenti; Tukey vale 0.05 per costruzione (gruppi bilanciati)
# H1: il gruppo 4 ha media +1,5; frazione di campioni in cui la coppia (1, 4) risulta significativa:
#   senza correzione  0.7460
#   Bonferroni        0.4853
#   Tukey             0.5174

Con tutte le medie uguali, il t-test a coppie senza correzione dà almeno un falso positivo nel 19,92% degli esperimenti. Il riferimento 1 − 0,95⁶ = 0,2649 varrebbe per sei confronti indipendenti, ma le coppie condividono i gruppi e il tasso reale è più basso. Bonferroni è a 4,05%, Tukey a 4,95%, vicinissimo al 5% che garantisce per costruzione a gruppi bilanciati.

Il prezzo della correzione si vede sotto H₁: la coppia con differenza vera è scoperta nel 74,6% dei casi senza correzione, nel 48,5% con Bonferroni e nel 51,7% con Tukey. Il 74,6% è un guadagno illusorio, comprato con falsi positivi a un tasso che non puoi dichiarare; fra i due metodi onesti, Tukey è più potente di Bonferroni perché sfrutta la struttura dei confronti. Se i confronti sono solo contro un gruppo di controllo, esiste il test di Dunnett, disponibile come scipy.stats.dunnett.

Quando le assunzioni non reggono

Il test F richiede indipendenza fra le osservazioni, Normalità dentro i gruppi e varianze uguali (omoschedasticità); la documentazione di f_oneway elenca le tre e indica come alternative Kruskal-Wallis e il test di Alexander-Govern, con un po' di potenza persa. Di tre, la più insidiosa è quella sulle varianze, soprattutto con gruppi di numerosità diversa: come per il t-test a due campioni dell'articolo precedente, il guaio nasce quando il gruppo più piccolo è anche il più disperso. La correzione classica è l'ANOVA di Welch, che pesa ogni gruppo con la precisione della sua media (w_i = n_i/s_i²) e corregge i gradi di libertà.

Il programma implementa la formula di Welch (1951) e la controlla in modo indiretto: con due gruppi la statistica deve coincidere con il quadrato del t di Welch e dare lo stesso p-value, cosa che il confronto con ttest_ind(equal_var=False) verifica. Poi simula tre gruppi con medie uguali in tre disposizioni di numerosità e deviazioni standard, e misura la frequenza con cui rifiutano il F classico, il F di Welch e Kruskal-Wallis. Con 10.000 ripetizioni per caso l'errore standard è circa 0,002.

# Varianze e numerosita' diverse: F classico, F di Welch e Kruskal-Wallis, sotto H0 (medie uguali)
import numpy as np
from scipy import stats

def welch_F(X_list):
    """ANOVA di Welch (1951) per gruppi di numerosita' diversa: restituisce statistica e p-value."""
    k = len(X_list)
    n = np.array([len(x) for x in X_list], dtype=float)
    medie = np.array([x.mean() for x in X_list])
    varianze = np.array([x.var(ddof=1) for x in X_list])
    w = n / varianze                                   # peso = precisione della media di ogni gruppo
    media_pesata = (w * medie).sum() / w.sum()
    lam = (((1 - w / w.sum()) ** 2) / (n - 1)).sum()
    F = ((w * (medie - media_pesata) ** 2).sum() / (k - 1)) / (1 + 2 * (k - 2) * lam / (k ** 2 - 1))
    df2 = (k ** 2 - 1) / (3 * lam)
    return F, stats.f.sf(F, k - 1, df2)

rng = np.random.default_rng(15)
# controllo: con due gruppi l'ANOVA di Welch coincide con il t-test di Welch (F = t^2, stesso p-value)
a, b = rng.normal(2, 3, 10), rng.normal(0, 1, 40)    # medie diverse, cosi' F non e' vicino a zero
F2, p2 = welch_F([a, b])
tw = stats.ttest_ind(a, b, equal_var=False)
print(f"due gruppi: F di Welch {F2:.4f} contro t^2 {tw.statistic ** 2:.4f}; p {p2:.6f} contro {tw.pvalue:.6f}")

R = 10_000
casi = {"gruppo piccolo = varianza maggiore": ((10, 20, 40), (3.0, 2.0, 1.0)),
        "gruppo piccolo = varianza minore   ": ((10, 20, 40), (1.0, 2.0, 3.0)),
        "stesse numerosita', varianze diverse": ((25, 25, 25), (3.0, 2.0, 1.0))}
print("caso (3 gruppi, medie uguali)             F classico  F di Welch  Kruskal-Wallis")
for nome, (ns, sig) in casi.items():
    r_f = r_w = r_k = 0
    for _ in range(R):
        g = [rng.normal(0, s, m) for m, s in zip(ns, sig)]
        r_f += stats.f_oneway(*g).pvalue < 0.05
        r_w += welch_F(g)[1] < 0.05
        r_k += stats.kruskal(*g).pvalue < 0.05
    print(f"{nome:38s} {r_f / R:10.4f}  {r_w / R:10.4f}  {r_k / R:13.4f}")
print(f"errore standard con R = {R}: circa {np.sqrt(0.05 * 0.95 / R):.4f}")

# Output:
# due gruppi: F di Welch 3.0844 contro t^2 3.0844; p 0.111590 contro 0.111590
# caso (3 gruppi, medie uguali)             F classico  F di Welch  Kruskal-Wallis
# gruppo piccolo = varianza maggiore         0.1931      0.0485         0.1121
# gruppo piccolo = varianza minore           0.0120      0.0485         0.0209
# stesse numerosita', varianze diverse       0.0618      0.0489         0.0592
# errore standard con R = 10000: circa 0.0022

Il controllo con due gruppi è perfetto: F di Welch = t² = 3,0844 e p = 0,111590 in entrambi. Quando il gruppo più piccolo (10 osservazioni) ha la deviazione standard maggiore (3, contro 2 e 1), il F classico rifiuta una H₀ vera il 19,3% delle volte invece del 5%, Welch il 4,85%. Nel caso opposto il classico rifiuta l'1,2%, quindi perde potenza senza avvisare.

Con numerosità uguali l'errore è piccolo (6,2%), e Welch resta a 4,9%. Kruskal-Wallis non è la soluzione: confronta i ranghi, non le medie, e con dispersioni diverse rifiuta l'11,2% delle volte anche con medie uguali. La regola che ne esce è quella già vista per il t-test: Welch non costa quasi nulla con varianze uguali e protegge con varianze diverse.

Questo articolo tratta la sola ANOVA a un fattore. L'analisi a due fattori, con le interazioni, e i modelli a effetti misti per misure ripetute sugli stessi soggetti usano la stessa scomposizione, ma richiedono un'impostazione propria e non sono qui; per le misure ripetute tieni presente che la prima assunzione, l'indipendenza, non vale e il F a una via non è applicabile.

Quattro errori che si commettono davvero

• Confrontare tutte le coppie con t-test non corretti: con quattro gruppi e medie uguali, almeno un falso positivo nel 19,9% dei casi.
• Leggere un F significativo come «tutte le medie sono diverse»: dice solo che non sono tutte uguali.
• Usare il F classico con gruppi di numerosità e varianze diverse: il livello vero è 19,3% quando il gruppo piccolo è il più disperso.
• Progettare con «30 per gruppo» e «effetto medio» senza calcolare la potenza: è 0,61, non 0,80.

Gli errori sono tutti di calibrazione: il test funziona come dichiarato solo se rispetti i suoi presupposti e dichiari il numero di confronti. Ogni numero della lista viene da una simulazione di questo articolo, e nessuno richiede più di qualche decina di secondi per essere rifatto; per questo conviene rifarlo sui tuoi dati, invece di fidarti di «30 per gruppo bastano».

Adesso sai scomporre la variabilità in SS_B e SS_W e calcolare F, p-value ed effetto a mano e con SciPy, calcolare la potenza con la F non centrale e il numero di dati per gruppo, confrontare le coppie con Tukey o Bonferroni controllando l'errore familywise, e scegliere Welch quando le varianze non sono uguali. Nel prossimo articolo passi dai dati continui ai dati in categorie, con il test chi-quadro, e nell'articolo sulla regressione lineare vedrai che l'ANOVA è un caso speciale di modello lineare.

Nella stessa serie trovi anche Inferenza statistica: stimatori e proprietà, Statistica non parametrica: Wilcoxon e Mann-Whitney e Analisi dati esplorativa (EDA) con Python.