07 - Test di ipotesi: t-test, z-test, p-value

Una macchina dovrebbe riempire confezioni da 500 g. Ne peschi 25 dalla linea, le pesi e trovi una media di 497,6 g con uno scarto tipo campionario di 6 g: la macchina è starata, oppure è solo la normale oscillazione del campione? Con i numeri dell'esempio la statistica vale esattamente −2,00, ma la conclusione dipende da quale distribuzione le dai in pasto: con la Normale il p-value è 0,0455 e rifiuti l'ipotesi che la macchina sia a posto, con la t di Student è 0,0569 e non la rifiuti.

Gli stessi dati, due verdetti opposti al livello del 5%. Questo articolo spiega perché, e ti fa passare da «ho una formula» a «so che cosa garantisce». È il settimo della serie Probabilità e Statistica; ti servono gli stimatori dell'articolo precedente e la Normale del teorema centrale del limite.

Come ragiona un test di ipotesi

Un test parte da due affermazioni sul parametro che non conosci. L'ipotesi nulla H₀ è quella da scartare solo davanti a prove forti (nell'esempio, «la media delle confezioni è 500 g»); l'alternativa H₁ è ciò che accetti se scarti la nulla («la media è diversa da 500 g»).

Poi scegli una statistica test, una funzione dei dati la cui distribuzione, se H₀ è vera, è nota: nell'esempio è T = (X̄ − 500)/(s/√n), cioè lo scarto fra la media osservata e quella dichiarata misurato in errori standard (qui 6/√25 = 1,2 g, dunque (497,6 − 500)/1,2 = −2). Infine fissi una regola: rifiuti H₀ se la statistica cade in una zona così lontana da essere poco probabile sotto H₀.

Una regola così può sbagliare in due modi. L'errore di primo tipo è rifiutare H₀ quando è vera (un falso allarme); la sua probabilità, che scegli tu prima di guardare i dati, si chiama livello α (di solito 0,05). L'errore di secondo tipo è non rifiutare H₀ quando è falsa (un'occasione persa); la probabilità di evitarlo è la potenza, e dipende da quanto è grande l'effetto vero, da quanti dati hai e da α. Le definizioni di errore, livello e potenza sono quelle della lezione sul test di ipotesi parametrico del corso MIT 18.650.

Il p-value è il più piccolo livello α al quale i tuoi dati avrebbero portato a rifiutare H₀ (stessa fonte): equivale alla probabilità, calcolata supponendo vera H₀, di ottenere una statistica estrema almeno quanto quella osservata. Nell'esempio, con la Normale, la probabilità di cadere oltre ±2 è 0,0455. Il p-value è una variabile aleatoria, perché dipende dal campione. Resta quindi il legame con gli intervalli di confidenza: il test bilaterale di livello 5% rifiuta H₀ esattamente quando il valore dichiarato (500 g) è fuori dall'intervallo di confidenza al 95%, come mostra il prossimo programma.

z o t: dipende da che cosa sai di σ

Il z-test vale quando la deviazione standard σ della popolazione è nota: la statistica Z = (X̄ − μ₀)/(σ/√n) è Normale standard sotto H₀ se i dati sono Normali.

Quando σ è ignota la sostituisci con lo scarto tipo campionario S, e la statistica diventa il rapporto fra una Normale e una quantità aleatoria; ha code più pesanti, ed è una t di Student con n − 1 gradi di libertà (il numero di scarti indipendenti che entrano in S). Per n piccolo la differenza è grande, per n grande sparisce. Il programma legge i dati dell'esempio in entrambi i modi: la stessa statistica −2, due code diverse.

# Le stesse 25 confezioni lette con z-test (sigma "nota" = 6) e con t-test (sigma stimata)
import numpy as np
from scipy import stats

n, media, s, mu0 = 25, 497.6, 6.0, 500.0
se = s / np.sqrt(n)                       # errore standard = 6/5 = 1.2 g
stat = (media - mu0) / se                 # stessa statistica numerica nei due test

p_z = 2 * stats.norm.sf(abs(stat))        # coda bilaterale della Normale standard
p_t = 2 * stats.t.sf(abs(stat), df=n - 1) # coda bilaterale della t con 24 gradi di liberta'
print(f"statistica = {stat:.3f}")
print(f"z-test: p = {p_z:.4f}  -> al 5% {'rifiuta' if p_z < 0.05 else 'non rifiuta'} H0")
print(f"t-test: p = {p_t:.4f}  -> al 5% {'rifiuta' if p_t < 0.05 else 'non rifiuta'} H0")

# Intervallo di confidenza al 95% costruito con la t: contiene 500 sse il t-test non rifiuta
q = stats.t.ppf(0.975, df=n - 1)
print(f"quantile t(24, 0.975) = {q:.4f} contro z(0.975) = {stats.norm.ppf(0.975):.4f}")
print(f"IC 95% (t): [{media - q * se:.2f}, {media + q * se:.2f}]")

# Output:
# statistica = -2.000
# z-test: p = 0.0455  -> al 5% rifiuta H0
# t-test: p = 0.0569  -> al 5% non rifiuta H0
# quantile t(24, 0.975) = 2.0639 contro z(0.975) = 1.9600
# IC 95% (t): [495.12, 500.08]

La statistica è identica, −2,000, ma la coda della t con 24 gradi di libertà vale 0,0569 contro 0,0455 della Normale, perché il quantile al 97,5% è 2,0639 e non 1,96. L'intervallo di confidenza al 95% costruito con la t è [495,12; 500,08] e contiene 500: coerente con il t-test che non rifiuta. In questo esempio σ non è nota (6 g è la stima del campione), quindi il test corretto è il t-test; il verdetto «macchina starata» del z-test è un falso allarme che nasce dall'aver usato la distribuzione sbagliata. Non è un'osservazione sul caso: il programma successivo misura quanto spesso succede.

Per misurarlo, simulo un milione di campioni da una Normale con media 0 (H₀ vera) per ciascun n, e conto con quale frequenza rifiuta il t-test con la sua soglia (quantile della t) e il «z-test con S» con la soglia 1,96. Sotto H₀ la statistica è esattamente una t con n − 1 gradi di libertà, quindi la frequenza esatta del secondo è 2 · P(t_n−1 > 1,96), che si calcola con la funzione di sopravvivenza della t.

# Tasso di falsi positivi sotto H0: z-test con sigma stimata contro t-test, per piccoli campioni
import numpy as np
from scipy import stats

R, BLOCCO, z = 1_000_000, 100_000, stats.norm.ppf(0.975)
rng = np.random.default_rng(99)
print("  n  t-test (sim.)  z con S (sim.)  z con S (esatto)")
for n in (5, 10, 30, 100):
    crit_t = stats.t.ppf(0.975, df=n - 1)              # soglia corretta: quantile della t
    rif_t = rif_z = 0
    for _ in range(R // BLOCCO):
        X = rng.normal(0.0, 2.0, size=(BLOCCO, n))     # H0 vera: media 0 (sigma = 2, ma il test non la usa)
        stat = np.sqrt(n) * X.mean(axis=1) / X.std(axis=1, ddof=1)
        rif_t += (np.abs(stat) > crit_t).sum()
        rif_z += (np.abs(stat) > z).sum()              # soglia sbagliata: quantile della Normale
    esatto = 2 * stats.t.sf(z, df=n - 1)               # sotto H0 la statistica e' una t(n-1)
    print(f"{n:3d}  {rif_t / R:12.4f}  {rif_z / R:14.4f}  {esatto:16.4f}")
print(f"errore standard di ogni valore simulato, con R = {R}: circa {np.sqrt(0.05 * 0.95 / R):.4f}")

# Output:
#   n  t-test (sim.)  z con S (sim.)  z con S (esatto)
#   5        0.0503          0.1222            0.1216
#  10        0.0501          0.0816            0.0816
#  30        0.0503          0.0599            0.0597
# 100        0.0501          0.0529            0.0528
# errore standard di ogni valore simulato, con R = 1000000: circa 0.0002

Il t-test rifiuta il 5,0% delle volte per ogni n, come promesso. Il «z-test con S» rifiuta il 12,2% delle volte con 5 osservazioni (esatto 12,16%), l'8,2% con 10, il 6,0% con 30 e il 5,3% con 100: la simulazione coincide con il calcolo esatto entro circa 0,0006. Ci vogliono un centinaio di osservazioni per rendere trascurabile l'errore di usare la Normale al posto della t, e questo è il motivo per cui la regola «con più di 30 dati puoi usare la Normale» è una scorciatoia accettabile solo se ti basta un livello vero di circa 6% al posto di 5%.

Due campioni: Student o Welch

Per confrontare le medie di due gruppi indipendenti, per esempio il tempo di risposta di due versioni di una pagina, esistono due t-test. Quello di Student assume che i due gruppi abbiano la stessa varianza, la stima con un valore comune (pooled) e usa n₁ + n₂ − 2 gradi di libertà. Quello di Welch stima le due varianze separatamente, usa come errore standard √(s₁²/n₁ + s₂²/n₂) e i gradi di libertà approssimati di Welch-Satterthwaite. In SciPy scipy.stats.ttest_ind ha equal_var=True come valore predefinito, cioè esegue Student: lo dice la documentazione della funzione. Per ottenere Welch devi scrivere equal_var=False.

Il programma scrive i due test in forma vettoriale, controlla che i p-value coincidano con quelli di SciPy su un campione, poi simula 200.000 coppie di campioni con la stessa media (H₀ vera) ma varianze diverse (deviazioni standard 3 e 1), in tre disposizioni: gruppo piccolo con varianza maggiore, gruppo piccolo con varianza minore, stesse numerosità. Per questo confronto non esiste una formula chiusa semplice, quindi la verifica è solo la simulazione, con il suo errore standard.

# Due campioni, medie uguali ma varianze diverse: t di Student (varianze uguali) contro Welch
import numpy as np
from scipy import stats

R = 200_000
rng = np.random.default_rng(8)

def p_values(a, b):
    """p-value bilaterali di Student (varianza comune) e di Welch; ogni riga e' un campione."""
    n1, n2 = a.shape[1], b.shape[1]
    m1, m2 = a.mean(axis=1), b.mean(axis=1)
    v1, v2 = a.var(axis=1, ddof=1), b.var(axis=1, ddof=1)
    # Student: varianza comune stimata (pooled), n1 + n2 - 2 gradi di liberta'
    sp2 = ((n1 - 1) * v1 + (n2 - 1) * v2) / (n1 + n2 - 2)
    t_st = (m1 - m2) / np.sqrt(sp2 * (1 / n1 + 1 / n2))
    p_st = 2 * stats.t.sf(np.abs(t_st), df=n1 + n2 - 2)
    # Welch: errori standard separati, gradi di liberta' di Welch-Satterthwaite
    e1, e2 = v1 / n1, v2 / n2
    t_w = (m1 - m2) / np.sqrt(e1 + e2)
    df_w = (e1 + e2) ** 2 / (e1 ** 2 / (n1 - 1) + e2 ** 2 / (n2 - 1))
    p_w = 2 * stats.t.sf(np.abs(t_w), df=df_w)
    return p_st, p_w

# controllo: le formule coincidono con scipy su un campione
a0, b0 = rng.normal(0, 3, (1, 10)), rng.normal(0, 1, (1, 40))
mio = p_values(a0, b0)
sci = (stats.ttest_ind(a0[0], b0[0], equal_var=True).pvalue, stats.ttest_ind(a0[0], b0[0], equal_var=False).pvalue)
print(f"p Student: mio {mio[0][0]:.6f}, scipy {sci[0]:.6f}; p Welch: mio {mio[1][0]:.6f}, scipy {sci[1]:.6f}")
print("caso                                      Student   Welch")
casi = {"campione piccolo ha la varianza maggiore": (10, 3.0, 40, 1.0),
        "campione piccolo ha la varianza minore": (10, 1.0, 40, 3.0),
        "numerosita' uguali, varianze diverse": (25, 3.0, 25, 1.0)}
for nome, (n1, s1, n2, s2) in casi.items():
    p_st, p_w = p_values(rng.normal(0, s1, (R, n1)), rng.normal(0, s2, (R, n2)))   # stessa media: H0 vera
    print(f"{nome:42s} {(p_st < 0.05).mean():6.4f}  {(p_w < 0.05).mean():6.4f}")
print(f"errore standard con R = {R}: circa {np.sqrt(0.05 * 0.95 / R):.4f}")

# Output:
# p Student: mio 0.040557, scipy 0.040557; p Welch: mio 0.268662, scipy 0.268662
# caso                                      Student   Welch
# campione piccolo ha la varianza maggiore   0.2568  0.0501
# campione piccolo ha la varianza minore     0.0014  0.0496
# numerosita' uguali, varianze diverse       0.0528  0.0493
# errore standard con R = 200000: circa 0.0005

Le prime due righe sono i verdetti importanti. Se il gruppo piccolo (10 osservazioni) è anche quello più disperso, Student rifiuta una H₀ vera il 25,7% delle volte invece del 5%, perché la varianza comune stimata è dominata dal gruppo grande e sottostima l'incertezza del piccolo. Nel caso opposto rifiuta quasi mai (0,14%), con una perdita di potenza che non si vede dal livello.

Welch resta intorno al 5% in tutti e tre i casi (0,0493-0,0501), e anche con numerosità uguali Student tiene (0,0528). Conclusione pratica: Welch non costa quasi nulla quando le varianze sono uguali e protegge quando non lo sono, quindi con SciPy conviene scrivere equal_var=False per abitudine.

Quanti dati servono: la potenza

La potenza risponde a una domanda da porsi prima di raccogliere i dati: se la differenza vera è quella che mi interessa, con quanti dati la scoprirò con probabilità alta? Si misura l'effetto in unità di deviazione standard, δ = (μ₁ − μ₀)/σ. Sotto l'alternativa la statistica del t-test segue una t non centrale con parametro di non centralità δ√n; scipy.stats.nct ne calcola le probabilità. La formula che si trova nei libri, n ≥ ((z_1−α/2 + z_0,80)/δ)², usa la Normale e quindi sottostima n di poco. Il programma confronta la formula con la potenza esatta e la conferma con una simulazione.

# Potenza del t-test a un campione (bilaterale, alfa = 0.05) e dimensione del campione
import numpy as np
from scipy import stats

alfa, delta = 0.05, 0.5            # delta = (mu1 - mu0)/sigma, effetto "medio" secondo Cohen
def potenza(n, d):
    """Esatta: sotto H1 la statistica e' una t non centrale con parametro d*sqrt(n)."""
    crit = stats.t.ppf(1 - alfa / 2, df=n - 1)
    nc = d * np.sqrt(n)
    return stats.nct.sf(crit, n - 1, nc) + stats.nct.cdf(-crit, n - 1, nc)

z = stats.norm.ppf
n_formula = ((z(1 - alfa / 2) + z(0.80)) / delta) ** 2     # formula con la Normale
print(f"formula Normale: n >= {n_formula:.2f} -> {int(np.ceil(n_formula))}")
for n in (32, 33, 34):
    print(f"n = {n}: potenza esatta del t-test = {potenza(n, delta):.4f}")
n_t = next(n for n in range(2, 1000) if potenza(n, delta) >= 0.80)
print(f"primo n con potenza >= 0.80 (t-test): {n_t}")

# Controllo per simulazione con n = 34
R, n = 200_000, 34
rng = np.random.default_rng(2026)
X = rng.normal(delta, 1.0, size=(R, n))
t = np.sqrt(n) * X.mean(axis=1) / X.std(axis=1, ddof=1)
sim = (np.abs(t) > stats.t.ppf(1 - alfa / 2, df=n - 1)).mean()
print(f"n = {n}: potenza simulata = {sim:.4f} +/- {np.sqrt(sim * (1 - sim) / R):.4f}")
for d in (0.4, 0.2):
    primo = next(n for n in range(2, 2000) if potenza(n, d) >= 0.80)
    print(f"delta = {d}: primo n con potenza >= 0.80 = {primo} (formula Normale: "
          f"{int(np.ceil(((z(1 - alfa / 2) + z(0.80)) / d) ** 2))})")
print(f"macchina: con n = 25 e delta = 0.4 la potenza esatta e' {potenza(25, 0.4):.4f}")

# Output:
# formula Normale: n >= 31.40 -> 32
# n = 32: potenza esatta del t-test = 0.7823
# n = 33: potenza esatta del t-test = 0.7954
# n = 34: potenza esatta del t-test = 0.8078
# primo n con potenza >= 0.80 (t-test): 34
# n = 34: potenza simulata = 0.8086 +/- 0.0009
# delta = 0.4: primo n con potenza >= 0.80 = 52 (formula Normale: 50)
# delta = 0.2: primo n con potenza >= 0.80 = 199 (formula Normale: 197)
# macchina: con n = 25 e delta = 0.4 la potenza esatta e' 0.4840

Per un effetto medio (δ = 0,5) la formula dice n ≥ 32, ma con 32 osservazioni la potenza esatta del t-test è 0,7823, con 33 è 0,7954 e solo con 34 arriva a 0,8078. La simulazione con 34 osservazioni dà 0,8086 con errore standard 0,0009.

Per l'esercizio della macchina, se vuoi accorgerti di uno scarto di 2,4 g con σ = 6 g (δ = 0,4) ti servono 52 confezioni, non 25: con 25 la potenza esatta è 0,4840, cioè meno di una volta su due, e un risultato «non significativo» come quello dell'esempio non dimostra che la macchina sia a posto. Con δ = 0,2 servono 199 osservazioni: dimezzare l'effetto da rilevare fa servire circa quattro volte i dati.

Che cosa dice un p-value, e che cosa no

Il comunicato dell'American Statistical Association del 7 marzo 2016 elenca sei principi sull'uso del p-value, riportati nel comunicato stampa dell'ASA; la dichiarazione completa, di Wasserstein e Lazar, è nel The American Statistician, 70(2), con DOI 10.1080/00031305.2016.1154108. Due principi riguardano le letture sbagliate più diffuse: un p-value indica quanto i dati sono incompatibili con un modello specificato, ma non misura la probabilità che l'ipotesi studiata sia vera, né la probabilità che i dati siano stati prodotti solo dal caso; e non misura la dimensione di un effetto né la sua importanza. Un terzo dice che le conclusioni non dovrebbero basarsi solo sul superamento di una soglia.

Il primo programma misura tre fatti sul p-value, con dati che non contengono alcun effetto. Primo, che sotto H₀ è uniforme su [0, 1], cioè che «p < 0,05» accade nel 5% dei casi, e anche «p < 0,50» nel 50%. Secondo, che 20 test indipendenti danno probabilità 1 − 0,95²⁰ = 0,6415 di avere almeno un p < 0,05, mentre la soglia di Bonferroni 0,05/20 riporta quella probabilità sotto il 5%. Terzo, che guardare i dati ogni 10 osservazioni fino a 100 e fermarsi al primo p < 0,05 gonfia i falsi positivi; qui non c'è formula chiusa, vale la simulazione.

# Sotto H0 il p-value e' uniforme su [0, 1]; che cosa succede con molti test o guardando i dati strada facendo
import numpy as np
from scipy import stats

def pvalue_t(X):
    """p-value bilaterale del t-test a un campione (H0: media 0) sulle righe di X."""
    n = X.shape[-1]
    t = np.sqrt(n) * X.mean(axis=-1) / X.std(axis=-1, ddof=1)
    return 2 * stats.t.sf(np.abs(t), df=n - 1)

rng = np.random.default_rng(20)

# 1) uniformita' del p-value quando H0 e' vera
R, n = 100_000, 30
p = pvalue_t(rng.normal(0, 1, size=(R, n)))
ks = stats.kstest(p, "uniform")
print(f"1) p < 0.05: {(p < 0.05).mean():.4f}; p < 0.50: {(p < 0.50).mean():.4f}; KS D = {ks.statistic:.4f}")

# 2) venti test indipendenti, nessun effetto vero: probabilita' di almeno un p < soglia
R, m = 20_000, 20
P = pvalue_t(rng.normal(0, 1, size=(R, m, n)))            # forma (R, m)
for nome, soglia in (("soglia 0.05      ", 0.05), ("Bonferroni 0.05/20", 0.05 / m)):
    sim = (P < soglia).any(axis=1).mean()
    print(f"2) {nome}: almeno un falso positivo, simulato {sim:.4f}, esatto {1 - (1 - soglia) ** m:.4f}")

# 3) guardare i dati ogni 10 osservazioni, fino a 100, e fermarsi al primo p < 0.05
R = 100_000
X = rng.normal(0, 1, size=(R, 100))
S1, S2 = X.cumsum(axis=1), (X ** 2).cumsum(axis=1)
fermati = np.zeros(R, dtype=bool)
for k in range(10, 101, 10):
    media = S1[:, k - 1] / k
    var = (S2[:, k - 1] - k * media ** 2) / (k - 1)
    pk = 2 * stats.t.sf(np.abs(np.sqrt(k) * media / np.sqrt(var)), df=k - 1)
    fermati |= pk < 0.05
print(f"3) 10 controlli intermedi: falsi positivi {fermati.mean():.4f} (nominale 0.05); solo alla fine: "
      f"{(pk < 0.05).mean():.4f}")

# Output:
# 1) p < 0.05: 0.0491; p < 0.50: 0.4995; KS D = 0.0015
# 2) soglia 0.05      : almeno un falso positivo, simulato 0.6459, esatto 0.6415
# 2) Bonferroni 0.05/20: almeno un falso positivo, simulato 0.0500, esatto 0.0488
# 3) 10 controlli intermedi: falsi positivi 0.1987 (nominale 0.05); solo alla fine: 0.0504

I tre risultati sono coerenti con la teoria. Il p-value ha frequenza 0,0491 sotto 0,05 e 0,4995 sotto 0,50, e la statistica di Kolmogorov-Smirnov contro l'uniforme è 0,0015, quindi non c'è scostamento osservabile. Con venti test la probabilità di almeno un falso positivo è 0,6459 simulata contro 0,6415 esatta (errore standard 0,0034 con 20.000 simulazioni), e 0,0500 contro 0,0488 con Bonferroni.

Con dieci controlli intermedi i falsi positivi salgono a 0,1987, quasi quattro volte il livello dichiarato, mentre il test unico alla fine resta al 5,04%. Un test di livello 5% vale solo se il piano (quanti test, quando fermarsi) è deciso prima di guardare i dati.

Resta il fraintendimento più dannoso: pensare che «p < 0,05» significhi «H₀ ha il 5% di probabilità di essere vera». La probabilità che H₀ sia vera dato un risultato significativo non si calcola dal p-value: dipende anche dalla frazione di ipotesi vere che stai mettendo alla prova e dalla potenza.

Il programma costruisce un mondo in cui una frazione π₀ delle ipotesi è nulla e le altre hanno un effetto δ = 0,5, con 34 osservazioni per test (potenza 0,8078), e conta fra i risultati significativi quanti sono nulle vere. La formula esatta è π₀α/(π₀α + (1 − π₀)·potenza), e si ottiene con la regola di Bayes.

# Fra i risultati "significativi", quanti sono ipotesi nulle vere? Dipende da quante ipotesi vere si mettono alla prova.
import numpy as np
from scipy import stats

alfa, n, delta = 0.05, 34, 0.5
crit = stats.t.ppf(1 - alfa / 2, df=n - 1)
nc = delta * np.sqrt(n)
potenza = stats.nct.sf(crit, n - 1, nc) + stats.nct.cdf(-crit, n - 1, nc)
print(f"alfa = {alfa}, potenza (n = {n}, effetto {delta}) = {potenza:.4f}")

rng = np.random.default_rng(5)
R, BLOCCO = 2_000_000, 100_000
print("quota di ipotesi nulle vere   esatto   simulato  (fra i risultati con p < 0.05)")
for pi0 in (0.5, 0.9, 0.99):
    significativi = nulle_fra_significativi = 0
    for _ in range(R // BLOCCO):
        vera_nulla = rng.random(BLOCCO) < pi0              # H0 vera con probabilita' pi0
        media = np.where(vera_nulla, 0.0, delta)           # effetto 0 oppure delta
        X = rng.normal(media[:, None], 1.0, size=(BLOCCO, n))
        t = np.sqrt(n) * X.mean(axis=1) / X.std(axis=1, ddof=1)
        signif = np.abs(t) > crit
        significativi += signif.sum()
        nulle_fra_significativi += (vera_nulla & signif).sum()
    esatto = pi0 * alfa / (pi0 * alfa + (1 - pi0) * potenza)
    print(f"{pi0:27.2f}  {esatto:8.4f}  {nulle_fra_significativi / significativi:8.4f}")

# Output:
# alfa = 0.05, potenza (n = 34, effetto 0.5) = 0.8078
# quota di ipotesi nulle vere   esatto   simulato  (fra i risultati con p < 0.05)
#                        0.50    0.0583    0.0577
#                        0.90    0.3578    0.3589
#                        0.99    0.8597    0.8592

Se metà delle ipotesi che metti alla prova è vera, il 5,8% dei risultati significativi sono falsi allarmi. Se ne è vera una su dieci, sono il 35,8%, e se una su cento, l'86%: la simulazione con due milioni di ipotesi coincide con la formula entro 0,11 punti percentuali. Il p-value è lo stesso (sotto 0,05) in tutti e tre i mondi, ma la probabilità che il risultato sia un falso allarme cambia di quindici volte. Per questo un solo p-value non basta: servono l'effetto stimato con il suo intervallo, la plausibilità a priori dell'ipotesi, e una conferma indipendente.

Cinque errori che si commettono davvero

• Usare la Normale al posto della t con campioni piccoli e σ stimata: il livello vero sale dal 5% al 12% con cinque osservazioni.
• Eseguire Student con ttest_ind senza accorgersi che è il valore predefinito, quando i gruppi hanno varianze e numerosità diverse: il livello vero arriva al 25,7%.
• Leggere un risultato non significativo come prova che non c'è differenza: con 25 osservazioni e un effetto di 0,4 deviazioni standard la potenza è 0,4840.
• Fermare la raccolta dati appena p scende sotto 0,05 o eseguire molti test e riportare il migliore: il livello nominale non vale più.
• Leggere il p-value come probabilità che l'ipotesi nulla sia vera: senza la frazione di ipotesi vere e senza la potenza, quella probabilità non si calcola.

I cinque errori nascono tutti dallo stesso scambio: trattare la garanzia di un test, che riguarda la procedura ripetuta su molti campioni e sotto condizioni precise, come una garanzia sul singolo risultato. Per ciascuno hai visto sopra un numero che mostra di quanto la garanzia si rompe, e ciascuno si controlla con un breve programma.

Adesso sai impostare un test (ipotesi, statistica, regola, livello), scegliere fra z, t, Student e Welch in base a ciò che sai dei dati, calcolare con scipy.stats.nct quanti dati servono per una potenza data, e dire che cosa un p-value afferma e che cosa no. Nell'articolo sull'ANOVA estenderai il confronto a più di due gruppi, dove il problema dei test multipli visto qui diventa il tema centrale, e in quello sul chi-quadro testerai dati categorici.

Nella stessa serie trovi anche Regressione lineare semplice e multipla, Statistica non parametrica: Wilcoxon e Mann-Whitney e Analisi dati esplorativa (EDA) con Python.