09 - Chi-quadrato: indipendenza e bontà di adattamento

Un A/B test ha dato 120 conversioni su 1000 visite per la versione A e 150 su 1000 per la B: un aumento relativo del 25%, o solo fortuna? I dati non sono misure ma conteggi in categorie (convertito sì o no, versione A o B), e il t-test dell'articolo precedente non si applica. Il test giusto è il chi-quadro (si scrive χ²), e sulla tua tabella dà p = 0,0496 senza correzione e p = 0,0577 con la correzione di continuità di Yates: un verdetto che dipende da una scelta che i manuali presentano come ovvia.

Questo articolo, il nono della serie Probabilità e Statistica, spiega che cosa calcola il test, quando vale l'approssimazione su cui poggia, e che cosa fare quando i dati sono pochi. Servono l'analisi della varianza e il p-value dell'articolo sui test di ipotesi.

Che cosa calcola la statistica χ²

Hai k categorie. Per ciascuna conti quante osservazioni ci sono (O_i, le frequenze osservate) e calcoli quante ne ti aspetteresti se l'ipotesi nulla fosse vera (E_i, le frequenze attese). La statistica di Pearson è χ² = Σ(O_i − E_i)²/E_i. Ogni addendo è lo scarto al quadrato, diviso per ciò che ti aspettavi: uno scarto di 10 su una cella con 1000 attesi pesa poco, lo stesso scarto su una cella con 20 attesi pesa molto. Se la statistica è piccola, i dati somigliano al modello; se è grande, non gli somigliano. La definizione è nel manuale di statistica del NIST (sezione 1.3.5.15).

La distribuzione di riferimento è la chi-quadro con d gradi di libertà: la legge della somma dei quadrati di d variabili Normali standard indipendenti, con media d e varianza 2d. Il numero di gradi di libertà conta quante quantità libere ci sono nel confronto.

Con k categorie e il totale fissato, solo k − 1 conteggi sono liberi: da qui d = k − 1 nel test di bontà di adattamento. Il motivo per cui la statistica ha questa legge è un'approssimazione, valida quando i conteggi sono grandi (il teorema centrale del limite applicato ai conteggi). Per questo il test è «asintotico»: il resto dell'articolo misura quanto bene regge a numeri piccoli.

Bontà di adattamento: un dado è equo?

Lanci un dado 60 volte e ottieni le facce 1-6 rispettivamente 8, 12, 9, 15, 6 e 10 volte. Se il dado è equo ti aspetti 10 per faccia, quindi χ² = ((−2)² + 2² + (−1)² + 5² + (−4)² + 0²)/10 = (4 + 4 + 1 + 25 + 16 + 0)/10 = 5,0, con 6 − 1 = 5 gradi di libertà. Il p-value è la probabilità che una chi-quadro con 5 gradi di libertà superi 5,0.

Il programma fa il calcolo e lo confronta con scipy.stats.chisquare. Poi misura due cose. La prima: il p-value ottenuto simulando un milione di dadi equi lanciati 60 volte, che non dipende dall'approssimazione. La seconda: quanto l'approssimazione sbaglia con pochi dati, su un caso con quattro categorie di probabilità 0,90, 0,05, 0,03 e 0,02. Con 4 categorie si possono enumerare tutti gli esiti possibili dei lanci e sommare le loro probabilità esatte, quindi la «dimensione vera» del test al 5% è esatta, non simulata.

# Bonta' di adattamento: un dado e' equo? Esercizio svolto, poi errore dell'approssimazione chi-quadro
import numpy as np
from scipy import stats

osservate = np.array([8, 12, 9, 15, 6, 10])          # 60 lanci
n, k = osservate.sum(), len(osservate)
attese = np.full(k, n / k)                           # 10 per faccia se il dado e' equo
chi2 = ((osservate - attese) ** 2 / attese).sum()
p = stats.chi2.sf(chi2, df=k - 1)
r = stats.chisquare(osservate, attese)
print(f"n = {n}, attese = {attese[0]:.0f}, chi2 = {chi2:.3f}, gradi di liberta' = {k - 1}, p = {p:.4f}")
print(f"scipy.stats.chisquare: chi2 = {r.statistic:.3f}, p = {r.pvalue:.4f}")

# p-value "quasi esatto": frazione di 1.000.000 di dadi equi, lanciati 60 volte, con chi2 almeno cosi' grande
rng = np.random.default_rng(7)
sim = rng.multinomial(n, np.full(k, 1 / k), size=1_000_000)
chi2_sim = ((sim - attese) ** 2 / attese).sum(axis=1)
print(f"p-value per simulazione = {(chi2_sim >= chi2 - 1e-9).mean():.4f} +/- 0.0005 (approssimazione chi-quadro: {p:.4f})")

# Dimensione vera del test al 5% con 4 categorie di probabilita' (0.90, 0.05, 0.03, 0.02):
# enumero TUTTI gli esiti multinomiali, quindi il risultato e' esatto
prob = np.array([0.90, 0.05, 0.03, 0.02])
crit = stats.chi2.ppf(0.95, df=3)
print("   n  attesa minima  dimensione esatta del test al 5%")
for n in (10, 20, 50, 100, 200):
    a, b, c = np.meshgrid(np.arange(n + 1), np.arange(n + 1), np.arange(n + 1), indexing="ij")
    ok = a + b + c <= n
    x = np.stack([a[ok], b[ok], c[ok], n - a[ok] - b[ok] - c[ok]], axis=1)
    pmf = stats.multinomial.pmf(x, n, prob)
    att = n * prob
    stat = ((x - att) ** 2 / att).sum(axis=1)
    print(f"{n:4d}  {att.min():13.2f}  {pmf[stat >= crit].sum():.4f}   (somma delle probabilita' = {pmf.sum():.6f})")

# Output:
# n = 60, attese = 10, chi2 = 5.000, gradi di liberta' = 5, p = 0.4159
# scipy.stats.chisquare: chi2 = 5.000, p = 0.4159
# p-value per simulazione = 0.4306 +/- 0.0005 (approssimazione chi-quadro: 0.4159)
#    n  attesa minima  dimensione esatta del test al 5%
#   10           0.20  0.0721   (somma delle probabilita' = 1.000000)
#   20           0.40  0.0750   (somma delle probabilita' = 1.000000)
#   50           1.00  0.0556   (somma delle probabilita' = 1.000000)
#  100           2.00  0.0515   (somma delle probabilita' = 1.000000)
#  200           4.00  0.0491   (somma delle probabilita' = 1.000000)

Il calcolo a mano coincide con SciPy: χ² = 5,000, p = 0,4159, e quindi il dado non dà alcuna evidenza di essere truccato. La simulazione dà 0,4306 (errore standard 0,0005): più dell'approssimazione, perché la statistica assume pochi valori distinti e «almeno 5,0» conta anche i lanci che danno esattamente 5,0.

Nella tabella della dimensione vera, con attesa minima 0,20 il test «al 5%» rifiuta in realtà il 7,21% delle volte, con 0,40 il 7,50%; con 1,00 il 5,56%, con 2,00 il 5,15% e con 4,00 il 4,91%. La soglia pratica che trovi nei manuali (frequenze attese di almeno 5; il manuale NIST e la documentazione di SciPy la riportano) è una regola prudente, non un confine netto: il danno cresce gradualmente man mano che le attese scendono sotto 2.

Parametri stimati dai dati: i gradi di libertà scendono

Spesso l'ipotesi nulla non è «il dado è equo» ma «i dati sono Poisson», senza dire con quale parametro. Devi stimare il parametro dai dati, e questo ha un costo: il modello si è già adattato ai conteggi, quindi gli scarti sono più piccoli di quelli di un modello fissato. Il teorema corretto, riportato nella lezione sul test di bontà di adattamento del corso MIT 18.650, dice che la statistica tende a una chi-quadro con K − d − 1 gradi di libertà, dove d è il numero di parametri stimati. Ogni parametro stimato toglie un grado di libertà.

Il programma lo verifica: genera 20.000 campioni di 200 osservazioni Poisson con parametro 3, li raggruppa in 8 categorie (0, 1, …, 6 e «7 o più»), stima il parametro con la massima verosimiglianza sui conteggi raggruppati (che è la stima per cui vale il teorema), calcola χ² e conta quante volte supera la soglia al 5% con K − 1 = 7 gradi di libertà e con K − 2 = 6.

# Se stimi un parametro dai dati, i gradi di liberta' scendono: k - 1 - (parametri stimati)
import numpy as np
from scipy import stats, optimize

lam_vero, n, R = 3.0, 200, 20_000
k = 8                                          # categorie: 0, 1, ..., 6 e "7 o piu'"
rng = np.random.default_rng(12)

def probabilita(lam):
    p = stats.poisson.pmf(np.arange(k - 1), lam)
    return np.append(p, 1 - p.sum())           # l'ultima categoria raccoglie la coda

def stima_ml_raggruppata(conteggi):
    """Massima verosimiglianza di lambda sui conteggi raggruppati (verosimiglianza multinomiale)."""
    f = lambda lam: -(conteggi * np.log(probabilita(lam))).sum()
    return optimize.minimize_scalar(f, bounds=(0.1, 10), method="bounded", options={"xatol": 1e-9}).x

chi2 = np.empty(R)
for r in range(R):
    x = np.minimum(rng.poisson(lam_vero, size=n), k - 1)       # dati gia' raggruppati
    o = np.bincount(x, minlength=k)
    lam = stima_ml_raggruppata(o)
    e = n * probabilita(lam)
    chi2[r] = ((o - e) ** 2 / e).sum()

print(f"attesa minima con lambda vero = {n * probabilita(lam_vero).min():.2f}; ripetizioni = {R}")
print("gradi di liberta'  soglia al 5%  rifiuti con quella soglia (H0 vera)")
for df in (k - 1, k - 2):
    soglia = stats.chi2.ppf(0.95, df)
    print(f"{df:17d}  {soglia:11.3f}  {(chi2 > soglia).mean():.4f}")
print(f"media simulata di chi2 = {chi2.mean():.3f} (k-1 = {k - 1}, k-2 = {k - 2}); errore standard del livello: {np.sqrt(0.05 * 0.95 / R):.4f}")

# Output:
# attesa minima con lambda vero = 6.70; ripetizioni = 20000
# gradi di liberta'  soglia al 5%  rifiuti con quella soglia (H0 vera)
#                 7       14.067  0.0303
#                 6       12.592  0.0478
# media simulata di chi2 = 5.987 (k-1 = 7, k-2 = 6); errore standard del livello: 0.0015

Con 7 gradi di libertà il test rifiuta il 3,03% delle volte, cioè è troppo prudente: perde potenza senza dirtelo. Con 6 rifiuta il 4,78%, vicino al 5% dichiarato (errore standard 0,0015). La media simulata della statistica, 5,987, è praticamente 6, il valore atteso di una chi-quadro con 6 gradi di libertà, e non 7. L'errore più comune è calcolare il parametro dai dati e poi usare comunque K − 1: il test è troppo prudente e accetta più facilmente modelli sbagliati. La stima va fatta sui conteggi raggruppati; con la stima sui dati grezzi la legge limite non è più esattamente quella.

Indipendenza in una tabella di contingenza

Una tabella di contingenza con r righe e c colonne contiene i conteggi per ogni combinazione di due variabili categoriche. L'ipotesi nulla è che le variabili siano indipendenti: nel tuo A/B test, che la probabilità di convertire non dipenda dalla versione. Sotto indipendenza, la frequenza attesa della cella (i, j) è E_ij = (totale riga i) × (totale colonna j)/n, perché le probabilità di riga e di colonna si moltiplicano.

I gradi di libertà sono (r − 1)(c − 1): le rc celle meno il totale, meno le r − 1 probabilità di riga e le c − 1 di colonna che hai stimato dai margini. Nel tuo caso, (2 − 1)(2 − 1) = 1. Con due campioni già separati (le due versioni) si parla di test di omogeneità, ma la statistica è la stessa.

Il programma applica la formula alla tabella dell'A/B test, la confronta con chi2_contingency con e senza correzione di Yates (che la funzione applica solo con un grado di libertà), mostra che per una 2 × 2 il test coincide con il confronto fra due proporzioni (z² = χ²), e misura la dimensione dell'effetto con tre indici: la V di Cramér (√(χ²/(n(min(r, c) − 1))), fra 0 e 1), l'odds ratio con l'intervallo di Woolf e la differenza fra le due proporzioni con il suo intervallo. Infine verifica i gradi di libertà su una tabella 3 × 4 con indipendenza vera, simulando 20.000 tabelle.

# Test di indipendenza su una tabella 2 x 2 (esito di un A/B test), a mano e con scipy
import numpy as np
from scipy import stats

tabella = np.array([[120, 880],      # versione A: 120 conversioni su 1000
                    [150, 850]])     # versione B: 150 conversioni su 1000
n = tabella.sum()
righe, colonne = tabella.sum(axis=1), tabella.sum(axis=0)
attese = np.outer(righe, colonne) / n                 # E_ij = (totale riga) * (totale colonna) / n
chi2 = ((tabella - attese) ** 2 / attese).sum()
gdl = (tabella.shape[0] - 1) * (tabella.shape[1] - 1)
print(f"attese = {attese.tolist()}")
print(f"chi2 = {chi2:.4f}, gradi di liberta' = {gdl}, p = {stats.chi2.sf(chi2, gdl):.4f}")
c, p, d, _ = stats.chi2_contingency(tabella, correction=False)
print(f"scipy (senza Yates): chi2 = {c:.4f}, p = {p:.4f}")
cy, py_, _, _ = stats.chi2_contingency(tabella, correction=True)
print(f"scipy (con Yates):   chi2 = {cy:.4f}, p = {py_:.4f}")

# Lo stesso test come confronto di due proporzioni: z^2 = chi2
pA, pB, pc = 120 / 1000, 150 / 1000, 270 / 2000
z = (pB - pA) / np.sqrt(pc * (1 - pc) * (1 / 1000 + 1 / 1000))
print(f"z = {z:.4f}, z^2 = {z ** 2:.4f}")

# Quanto e' grande l'effetto: V di Cramer, rapporto fra odds e intervallo di Woolf, differenza fra proporzioni
V = np.sqrt(chi2 / (n * (min(tabella.shape) - 1)))
OR = tabella[1, 0] * tabella[0, 1] / (tabella[0, 0] * tabella[1, 1])
se_log = np.sqrt((1 / tabella).sum())
lo, hi = np.exp(np.log(OR) + np.array([-1.96, 1.96]) * se_log)
se_d = np.sqrt(pA * (1 - pA) / 1000 + pB * (1 - pB) / 1000)
print(f"V di Cramer = {V:.4f}")
print(f"odds ratio B contro A = {OR:.4f}, IC 95% di Woolf = [{lo:.4f}, {hi:.4f}]")
print(f"differenza fra proporzioni = {pB - pA:.3f}, IC 95% = [{pB - pA - 1.96 * se_d:.4f}, {pB - pA + 1.96 * se_d:.4f}]")

# Verifica dei gradi di liberta' con una tabella 3 x 4 e indipendenza vera: 20.000 tabelle simulate
rng = np.random.default_rng(3)
pr, pc_ = np.array([0.5, 0.3, 0.2]), np.array([0.4, 0.3, 0.2, 0.1])
celle = np.outer(pr, pc_).ravel()
R, m = 20_000, 300
stat = np.empty(R)
for r in range(R):
    t = rng.multinomial(m, celle).reshape(3, 4)
    stat[r] = stats.chi2_contingency(t, correction=False)[0]
soglia = stats.chi2.ppf(0.95, 6)
print(f"tabella 3 x 4, n = {m}: gradi di liberta' 6, media di chi2 = {stat.mean():.3f}, rifiuti al 5% = {(stat > soglia).mean():.4f}")

# Output:
# attese = [[135.0, 865.0], [135.0, 865.0]]
# chi2 = 3.8536, gradi di liberta' = 1, p = 0.0496
# scipy (senza Yates): chi2 = 3.8536, p = 0.0496
# scipy (con Yates):   chi2 = 3.6009, p = 0.0577
# z = 1.9630, z^2 = 3.8536
# V di Cramer = 0.0439
# odds ratio B contro A = 1.2941, IC 95% di Woolf = [0.9999, 1.6748]
# differenza fra proporzioni = 0.030, IC 95% = [0.0001, 0.0599]
# tabella 3 x 4, n = 300: gradi di liberta' 6, media di chi2 = 5.958, rifiuti al 5% = 0.0488

Le frequenze attese sono 135 e 865 in entrambe le righe, e χ² = 3,8536, p = 0,0496: con la correzione di Yates, 3,6009 e p = 0,0577. La statistica z del confronto fra proporzioni è 1,9630 e il suo quadrato è identico al χ². Il verdetto «significativo» o «non significativo» dipende da una correzione, e questo è un segnale: il risultato è al limite. Gli indici di effetto lo confermano.

La V di Cramér è 0,0439, cioè quasi nulla; l'odds ratio è 1,2941 con intervallo [0,9999; 1,6748], che sfiora 1; la differenza fra proporzioni è 3 punti percentuali con intervallo [0,0001; 0,0599]. L'effetto probabile è piccolo e incerto: i dati non bastano per decidere, e un p-value di 0,0496 senza questi numeri è una cattiva sintesi. Sulla tabella 3 × 4, la media di χ² è 5,958, vicina ai 6 gradi di libertà (3 − 1)(4 − 1), e rifiuta il 4,88% delle volte.

Pochi dati: Yates e Fisher

Quando le frequenze attese sono piccole, le alternative sono due. La correzione di Yates riduce ogni scarto |O − E| di 0,5 prima di elevarlo al quadrato, per compensare il fatto che i conteggi sono interi.

Il test esatto di Fisher non usa l'approssimazione: condiziona ai totali di riga e di colonna e calcola la probabilità esatta di ogni tabella con quei totali, che è una distribuzione ipergeometrica; scipy.stats.fisher_exact lo implementa per tabelle 2 × 2 (dalla versione 1.15 accetta anche tabelle più grandi, con il parametro method). Nessuno dei due è «il migliore»: sono due risposte diverse alla stessa difficoltà, e conviene vedere che cosa fanno.

Il programma prende la tabella [[8, 2], [1, 5]], con frequenza attesa minima 2,62, e calcola i tre p-value. Poi, per due gruppi da 10 osservazioni con probabilità di successo p₁ e p₂, enumera tutte le 121 tabelle possibili, somma le probabilità di quelle in cui ciascun test rifiuta al 5%, e ottiene le probabilità esatte di rifiutare: il livello vero quando p₁ = p₂ e la potenza quando differiscono. Le tabelle con una colonna vuota (nessun successo o nessun fallimento in totale) non danno informazione e non rifiutano.

# Tabelle 2 x 2 con pochi dati: dimensione vera (esatta) di chi-quadro, chi-quadro con Yates e test di Fisher
import numpy as np
from scipy import stats

print("tabella di esempio [[8, 2], [1, 5]]:")
t = np.array([[8, 2], [1, 5]])
print(f"  attese minime = {stats.contingency.expected_freq(t).min():.2f}")
print(f"  chi2 senza correzione: p = {stats.chi2_contingency(t, correction=False)[1]:.4f}")
print(f"  chi2 con Yates:        p = {stats.chi2_contingency(t, correction=True)[1]:.4f}")
print(f"  Fisher:                p = {stats.fisher_exact(t)[1]:.4f}")

# Due gruppi da 10, probabilita' di successo p1 e p2: enumero tutte le 121 tabelle possibili,
# quindi la probabilita' di rifiutare al 5% e' ESATTA (livello se p1 = p2, potenza se p1 != p2)
n1 = n2 = 10

def prob_rifiuto(p1, p2):
    rif = np.zeros(3)
    for x1 in range(n1 + 1):
        for x2 in range(n2 + 1):
            if x1 + x2 in (0, n1 + n2):                    # una colonna vuota: nessuna evidenza, non si rifiuta
                continue
            tab = np.array([[x1, n1 - x1], [x2, n2 - x2]])
            w = stats.binom.pmf(x1, n1, p1) * stats.binom.pmf(x2, n2, p2)
            rif += w * np.array([stats.chi2_contingency(tab, correction=False)[1] < 0.05,
                                 stats.chi2_contingency(tab, correction=True)[1] < 0.05,
                                 stats.fisher_exact(tab)[1] < 0.05])
    return rif

print("situazione           chi2    Yates   Fisher   (probabilita' esatta di rifiutare al 5%)")
for etichetta, p1, p2 in (("H0, p = 0.1", 0.1, 0.1), ("H0, p = 0.3", 0.3, 0.3), ("H0, p = 0.5", 0.5, 0.5),
                          ("H1, 0.2 contro 0.8", 0.2, 0.8), ("H1, 0.1 contro 0.6", 0.1, 0.6)):
    r = prob_rifiuto(p1, p2)
    print(f"{etichetta:19s} {r[0]:7.4f} {r[1]:7.4f} {r[2]:7.4f}")

# Output:
# tabella di esempio [[8, 2], [1, 5]]:
#   attese minime = 2.62
#   chi2 senza correzione: p = 0.0134
#   chi2 con Yates:        p = 0.0510
#   Fisher:                p = 0.0350
# situazione           chi2    Yates   Fisher   (probabilita' esatta di rifiutare al 5%)
# H0, p = 0.1          0.0090  0.0011  0.0011
# H0, p = 0.3          0.0371  0.0119  0.0119
# H0, p = 0.5          0.0422  0.0128  0.0128
# H1, 0.2 contro 0.8   0.8054  0.6353  0.6353
# H1, 0.1 contro 0.6   0.6590  0.4740  0.4740

Sulla tabella di esempio i tre test non concordano al 5%: χ² senza correzione dà p = 0,0134, Yates 0,0510, Fisher 0,0350. Con attese di 2,62 nessuno dei tre è «quello vero» per definizione, ma Fisher è esatto condizionatamente ai totali. Il confronto esatto smentisce due idee comuni. Primo, a questa taglia il χ² non inflaziona i falsi positivi: rifiuta lo 0,90%, il 3,71% e il 4,22% con p uguale a 0,1, 0,3 e 0,5, quindi è prudente.

Secondo, Fisher e Yates sono ancora più prudenti, 0,11%-1,28%, e la prudenza si paga in potenza: contro un'alternativa 0,2 contro 0,8, il χ² rifiuta nell'80,5% dei casi, Yates e Fisher nel 63,5%; con 0,1 contro 0,6, 65,9% contro 47,4%. La scelta è quindi fra un test che sbaglia poco ma sfrutta poco i dati, e uno più potente che dichiara un livello non garantito.

Quattro errori che si commettono davvero

• Usare K − 1 gradi di libertà dopo aver stimato un parametro dai dati: il livello vero è 3,0% invece di 5%.
• Dichiarare un risultato «significativo» con p = 0,0496 senza guardare l'effetto: nell'A/B test V = 0,0439 e l'intervallo della differenza parte da 0,0001.
• Applicare il χ² con attese sotto 2 su categorie rare: il livello vero sale al 7,5%.
• Scegliere Yates o Fisher «per sicurezza» senza accorgersi che si perde potenza: dall'80,5% al 63,5% in uno dei casi esatti calcolati sopra.

Gli errori hanno in comune la fiducia in una regola non verificata: la soglia di 5, la correzione automatica, i gradi di libertà a memoria. Ciascuna regola ha una ragione e un costo, e il costo si misura con un programma che enumera le tabelle o simula i conteggi; se i tuoi dati stanno vicino a una soglia, rifallo con i tuoi numeri.

Adesso sai calcolare a mano e con SciPy la statistica χ² per una bontà di adattamento e per una tabella di contingenza, scegliere i gradi di libertà (anche con parametri stimati), riconoscere quando l'approssimazione è debole, confrontare χ², Yates e Fisher sui casi con pochi dati, e riportare l'effetto insieme al p-value. Nel prossimo articolo passi a variabili continue legate fra loro, con la regressione lineare.

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