15 - Analisi dati esplorativa (EDA) con Python

Ti consegnano un file di dati e ti chiedono un modello. Prima di scegliere l'algoritmo, devi sapere che cosa hai in mano: quante righe, quali colonne sono numeri e quali etichette, dove mancano valori, come sono fatte le distribuzioni, quali variabili si muovono insieme e se qualche valore è così lontano dagli altri da guidare ogni calcolo. L'analisi esplorativa dei dati (EDA, exploratory data analysis) è questo lavoro. Il manuale di statistica del NIST la descrive come un approccio, non un elenco di tecniche, che usa soprattutto grafici per far emergere la struttura dei dati, individuare valori anomali e controllare le ipotesi sottostanti, rimandando la scelta del modello; la sua opera fondativa è il libro di Tukey del 1977. In questo articolo, l'ultimo della serie Probabilità e Statistica, impari a farla in Python con pandas e matplotlib su un dataset reale, e vedi con i numeri perché le statistiche riassuntive da sole non bastano mai. È il passo che precede i test, gli intervalli e i modelli degli articoli precedenti.

Perché non basta riassumere: il quartetto di Anscombe

Nel 1973 Francis Anscombe pubblicò quattro piccoli dataset, di 11 coppie (x, y) ciascuno, costruiti per avere quasi le stesse statistiche riassuntive e forme completamente diverse. Il listato li ricostruisce dai valori pubblicati e calcola per ognuno media e varianza di x e y, correlazione di Pearson e retta dei minimi quadrati. Poi guarda ciò che quei numeri nascondono: nel secondo set misura quanto meglio una parabola descrive i dati rispetto alla retta; nel terzo confronta la pendenza dei minimi quadrati con quella robusta di Theil-Sen (la mediana delle pendenze di tutte le coppie di punti); nel quarto conta i valori distinti di x.

import numpy as np
from scipy import stats

# Il quartetto di Anscombe (1973): quattro dataset di 11 punti
x123 = np.array([10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5], float)
x4 = np.array([8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8], float)
dati = {
    "I":   (x123, np.array([8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68])),
    "II":  (x123, np.array([9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74])),
    "III": (x123, np.array([7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73])),
    "IV":  (x4,   np.array([6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89])),
}
print("set   media x  media y  var x  var y   r     retta y = a + b x")
for nome, (x, y) in dati.items():
    b, a = np.polyfit(x, y, 1)
    print(f"{nome:5s} {x.mean():6.2f} {y.mean():8.2f} {x.var(ddof=1):6.2f} {y.var(ddof=1):6.2f} "
          f"{np.corrcoef(x, y)[0, 1]:6.3f}   y = {a:.2f} + {b:.3f} x")

# Che cosa vede l'occhio e le statistiche no
x, y = dati["II"]
quad = np.polyfit(x, y, 2)
r2_lin = np.corrcoef(x, y)[0, 1] ** 2
r2_quad = 1 - np.sum((y - np.polyval(quad, x)) ** 2) / np.sum((y - y.mean()) ** 2)
print(f"II : R2 retta {r2_lin:.3f}  R2 parabola {r2_quad:.3f}  (la relazione e' curva)")
x, y = dati["III"]
ts = stats.theilslopes(y, x)
b_ols, a_ols = np.polyfit(x, y, 1)
print(f"III: pendenza con minimi quadrati {b_ols:.3f}  con Theil-Sen {ts.slope:.3f}  (un punto trascina la retta)")
x, y = dati["IV"]
print(f"IV : valori distinti di x = {sorted(set(x.tolist()))}  senza il punto x=19 tutti gli x sono uguali: pendenza non definita")

# Output:
# set   media x  media y  var x  var y   r     retta y = a + b x
# I       9.00     7.50  11.00   4.13  0.816   y = 3.00 + 0.500 x
# II      9.00     7.50  11.00   4.13  0.816   y = 3.00 + 0.500 x
# III     9.00     7.50  11.00   4.12  0.816   y = 3.00 + 0.500 x
# IV      9.00     7.50  11.00   4.12  0.817   y = 3.00 + 0.500 x
# II : R2 retta 0.666  R2 parabola 1.000  (la relazione e' curva)
# III: pendenza con minimi quadrati 0.500  con Theil-Sen 0.346  (un punto trascina la retta)
# IV : valori distinti di x = [8.0, 19.0]  senza il punto x=19 tutti gli x sono uguali: pendenza non definita

I quattro set hanno medie 9,00 e 7,50, varianze di x pari a 11,00 e di y fra 4,12 e 4,13, correlazione 0,816-0,817 e la stessa retta, y = 3,00 + 0,500 x. Eppure il primo è una nuvola lineare con rumore, nel secondo la relazione è una curva perfetta (R² della retta 0,666, della parabola 1,000), nel terzo un solo punto anomalo trascina la pendenza a 0,500 mentre la stima robusta dà 0,346, nel quarto tutti gli x valgono 8 tranne uno a 19, e senza quel punto la pendenza non esisterebbe. Nessun numero di quella tabella avrebbe fatto sospettare tutto questo; un grafico a dispersione sì. È la ragione per cui l'EDA comincia guardando i dati prima di riassumerli.

Il primo sguardo a un dataset reale

Come esempio uso il dataset wine, incluso in scikit-learn: 178 vini di una stessa regione italiana, 13 misure chimiche per ciascuno e una classe fra tre produttori (la descrizione è nella documentazione di load_wine). Il primo gesto è sempre lo stesso: la forma della tabella, i tipi delle colonne e i valori mancanti. Poi una sintesi numerica delle variabili, con la media e la mediana affiancate, perché una distanza fra le due è già un indizio di asimmetria. Infine la recinzione di Tukey: si considera sospetto ogni valore che cade a più di 1,5 volte lo scarto interquartile (IQR = Q3 − Q1) sotto il primo quartile o sopra il terzo, la stessa regola con cui matplotlib disegna i baffi dei boxplot. L'asimmetria (skewness) è positiva se la coda destra è più lunga.

import numpy as np
import pandas as pd
from sklearn.datasets import load_wine

wine = load_wine(as_frame=True)
df = wine.frame.rename(columns={"target": "classe"})
print("forma:", df.shape, "  valori mancanti:", int(df.isna().sum().sum()))
print("tipi:", df.dtypes.astype(str).value_counts().to_dict())
print("righe per classe:", df["classe"].value_counts().sort_index().to_dict())

cols = ["alcohol", "malic_acid", "ash", "magnesium", "proline"]
riassunto = df[cols].describe().loc[["mean", "std", "min", "50%", "max"]].round(2)
print(riassunto.to_string())

# Forma delle distribuzioni: asimmetria e valori fuori dalle "recinzioni" di Tukey (Q1 - 1.5 IQR, Q3 + 1.5 IQR)
q1, q3 = df[cols].quantile(.25), df[cols].quantile(.75)
iqr = q3 - q1
fuori = ((df[cols] < q1 - 1.5 * iqr) | (df[cols] > q3 + 1.5 * iqr)).sum()
print(pd.DataFrame({"asimmetria": df[cols].skew().round(2), "fuori recinzione": fuori}).to_string())

# Output:
# forma: (178, 14)   valori mancanti: 0
# tipi: {'float64': 13, 'int64': 1}
# righe per classe: {0: 59, 1: 71, 2: 48}
#       alcohol  malic_acid   ash  magnesium  proline
# mean    13.00        2.34  2.37      99.74   746.89
# std      0.81        1.12  0.27      14.28   314.91
# min     11.03        0.74  1.36      70.00   278.00
# 50%     13.05        1.87  2.36      98.00   673.50
# max     14.83        5.80  3.23     162.00  1680.00
#             asimmetria  fuori recinzione
# alcohol          -0.05                 0
# malic_acid        1.04                 3
# ash              -0.18                 3
# magnesium         1.10                 4
# proline           0.77                 0

La tabella ha 178 righe e 14 colonne (13 numeriche reali più la classe intera), nessun valore mancante e classi da 59, 71 e 48 vini: i conteggi coincidono con la descrizione ufficiale del dataset. Per la proline, la media (746,89) supera la mediana (673,50) di oltre 70 unità, e l'asimmetria è 0,77, segno di una coda destra; per l'alcol, media e mediana quasi coincidono (13,00 e 13,05) e l'asimmetria è −0,05. La regola di Tukey segnala 3 valori fuori recinzione per l'acido malico, 3 per la cenere, 4 per il magnesio e nessuno per alcol e proline. Segnalare non significa eliminare: un valore fuori recinzione può essere un errore di battitura, una misura di un'unità diversa o un caso reale e raro, e le tre ipotesi chiedono azioni opposte. Prima di toglierlo si controlla da dove viene.

Relazioni fra variabili, e che cosa non dicono

Il passo successivo guarda le variabili a coppie. La correlazione di Pearson misura la tendenza lineare e quella di Spearman (articolo 13) quella monotona: se divergono molto, la relazione non è una retta o c'è qualche valore anomalo. Il listato elenca le tre coppie più correlate, mostra le medie di alcune variabili per ciascuna classe e calcola il rapporto F dell'ANOVA a una via per la proline, cioè il rapporto fra la variabilità fra le classi e quella dentro le classi, a mano e con scipy.stats.f_oneway (vedi ANOVA e confronti multipli).

import numpy as np
import pandas as pd
from scipy import stats
from sklearn.datasets import load_wine

wine = load_wine(as_frame=True)
df = wine.frame.rename(columns={"target": "classe"})

# Relazioni: le coppie di variabili piu' correlate (Pearson e Spearman)
feat = df.drop(columns="classe")
pearson = feat.corr().where(np.triu(np.ones((13, 13), bool), 1)).stack()
top = pearson.reindex(pearson.abs().sort_values(ascending=False).index).head(3)
for (a, b), r in top.items():
    print(f"{a} ~ {b}: Pearson {r:.3f}  Spearman {feat[a].corr(feat[b], method='spearman'):.3f}")

# Una variabile separa le classi? media per classe e rapporto fra varianza tra classi e dentro le classi
medie = df.groupby("classe")[["alcohol", "proline", "ash"]].mean().round(2)
print(medie.to_string())
tra = df.groupby("classe")["proline"].agg(["mean", "count"])
F = (tra["count"] * (tra["mean"] - df["proline"].mean()) ** 2).sum() / 2 / (df.groupby("classe")["proline"].var() * (tra["count"] - 1)).sum() * (len(df) - 3)
gruppi = [g["proline"].to_numpy() for _, g in df.groupby("classe")]
print(f"proline: F di Fisher a mano = {F:.1f}   scipy.stats.f_oneway = {stats.f_oneway(*gruppi).statistic:.1f}")

# Output:
# total_phenols ~ flavanoids: Pearson 0.865  Spearman 0.879
# flavanoids ~ od280/od315_of_diluted_wines: Pearson 0.787  Spearman 0.742
# total_phenols ~ od280/od315_of_diluted_wines: Pearson 0.700  Spearman 0.687
#         alcohol  proline   ash
# classe                        
# 0         13.74  1115.71  2.46
# 1         12.28   519.51  2.24
# 2         13.15   629.90  2.44
# proline: F di Fisher a mano = 207.9   scipy.stats.f_oneway = 207.9

Le tre coppie più correlate sono tutte fenoli: total_phenols e flavanoids (Pearson 0,865, Spearman 0,879), poi flavanoids e od280/od315 (0,787 e 0,742) e total_phenols con od280/od315 (0,700 e 0,687). Sono misure della stessa famiglia chimica, quindi non indipendenti: un modello che le usa tutte insieme dovrà tenerne conto. La proline varia molto con la classe: media 1115,71 nella classe 0 contro 519,51 e 629,90 nelle altre, con F = 207,9 sia a mano sia con SciPy. Questo dice che la variabile separa bene le classi; non dice che ne sia la causa, e non dice quale classe di vino abbia «più proline» per un motivo che il dataset spieghi.

Il paradosso di Simpson: quando il totale mente

Una trappola dell'EDA è aggregare troppo presto. Nel paradosso di Simpson (descritto nel 1951) una tendenza che compare in ogni sottogruppo si inverte, o sparisce, nei dati aggregati, perché i gruppi non hanno la stessa composizione. L'esempio seguente è costruito, e i numeri sono inventati: due varianti di una pagina, A e B, e due tipi di visitatore, esperti e principianti, i secondi molto meno propensi a convertire.

import pandas as pd

# Dati INVENTATI per l'esempio: conversioni di due varianti, divise per tipo di visitatore
righe = [("A", "esperto",   90,  100), ("A", "principiante", 160, 400),
         ("B", "esperto",  800, 1000), ("B", "principiante",  10, 100)]
t = pd.DataFrame(righe, columns=["variante", "gruppo", "conversioni", "visite"])
t["tasso"] = (t["conversioni"] / t["visite"]).round(3)
print(t.to_string(index=False))

tot = t.groupby("variante")[["conversioni", "visite"]].sum()
tot["tasso"] = (tot["conversioni"] / tot["visite"]).round(3)
print(tot.to_string())

# Perche': le due varianti non vedono lo stesso mix di visitatori
mix = t.pivot(index="variante", columns="gruppo", values="visite")
print((mix.div(mix.sum(axis=1), axis=0) * 100).round(1).to_string())

# Output:
# variante       gruppo  conversioni  visite  tasso
#        A      esperto           90     100    0.9
#        A principiante          160     400    0.4
#        B      esperto          800    1000    0.8
#        B principiante           10     100    0.1
#           conversioni  visite  tasso
# variante                            
# A                 250     500  0.500
# B                 810    1100  0.736
# gruppo    esperto  principiante
# variante                       
# A            20.0          80.0
# B            90.9           9.1

La variante A converte meglio di B in entrambi i gruppi (90% contro 80% tra gli esperti, 40% contro 10% tra i principianti) ma peggio sul totale (50,0% contro 73,6%). La spiegazione sta nell'ultima tabella: ad A arriva l'80% di principianti, a B il 9,1%. Il confronto aggregato misura la variante e il mix di visitatori insieme, e non permette di separarli. Quale lettura è giusta dipende da come sono stati assegnati i visitatori: se casualmente, i due mix dovrebbero somigliarsi e uno così diverso indica un difetto nell'assegnazione; se no, il confronto per gruppo è più affidabile. Il modo di accorgersene è esplorare i dati spezzati per le variabili che ne spiegano la composizione. Nessun test dopo questo passo rimedia a un'aggregazione sbagliata (vedi anche Test and learn: esperimenti controllati e A/B test).

I grafici: quale per quale domanda

I numeri riassumono, i grafici mostrano la forma. La scelta dipende dalla domanda che stai ponendo ai dati, e la tabella segue questo criterio:

Il listato produce tre pannelli sul dataset wine: l'istogramma della proline, i boxplot della stessa variabile per classe e il diagramma a dispersione di flavanoids contro total_phenols colorato per classe. Usa il backend Agg di matplotlib, che disegna su file senza aprire finestre: è ciò che serve su un server o in uno script. L'immagine non è riprodotta nell'articolo; l'uscita del programma conferma che i tre pannelli esistono e che l'ultimo contiene i 178 punti.

import matplotlib
matplotlib.use("Agg")                       # nessuna finestra: salva su file
import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import load_wine

df = load_wine(as_frame=True).frame
fig, ax = plt.subplots(1, 3, figsize=(12, 3.6))
ax[0].hist(df["proline"], bins=20, edgecolor="white")          # forma: una gobba, coda a destra
ax[0].set(title="proline: distribuzione", xlabel="proline", ylabel="conteggio")
df.boxplot(column="proline", by="target", ax=ax[1])             # confronto fra gruppi
ax[1].set(title="proline per classe", xlabel="classe", ylabel="proline")
for c, g in df.groupby("target"):                               # relazione, colorata per gruppo
    ax[2].scatter(g["flavanoids"], g["total_phenols"], s=12, label=f"classe {c}")
ax[2].set(title="fenoli e flavonoidi", xlabel="flavanoids", ylabel="total_phenols"); ax[2].legend()
fig.suptitle("")
fig.tight_layout()
fig.savefig("eda_wine.png", dpi=80)
print("pannelli:", len(fig.axes), " punti nell'ultimo grafico:", sum(len(c.get_offsets()) for c in ax[2].collections))

# Output:
# pannelli: 3  punti nell'ultimo grafico: 178

La riga df.boxplot(column="proline", by="target", ax=ax[1]) fa il confronto fra gruppi senza ciclo; ax[2].scatter è dentro un ciclo sulle classi per assegnare un colore a ciascuna e una voce alla legenda. Se togli matplotlib.use("Agg") su una macchina senza schermo, il programma può fallire alla creazione della figura. Aprendo l'immagine dovresti ritrovare ciò che i numeri hanno già detto: una coda a destra nella proline (asimmetria 0,77), la classe 0 con valori mediamente più alti (media 1115,71) e, nel diagramma a dispersione, punti che seguono approssimativamente una retta (correlazione 0,865). I grafici non aggiungono fatti, li rendono più facili da cogliere e da mostrare a chi non legge tabelle.

Esercizio svolto e prova da fare da solo

Un'azienda misura i tempi di consegna in giorni: 2, 3, 3, 4, 4, 4, 5, 5, 6, 21. La media è 57/10 = 5,7; la mediana, media dei due valori centrali (4 e 4), è 4,0. Il primo quartile è 3,25 e il terzo 5,0 (interpolazione lineare, come in NumPy), quindi IQR = 1,75 e la recinzione superiore è 5,0 + 1,5 · 1,75 = 7,625: il 21 è fuori, gli altri no. Senza il 21 la media scende a 36/9 = 4,0 e la mediana resta 4,0. Un solo valore sposta la media di 1,7 giorni e la mediana di zero: la mediana è robusta, la media no. Il 21 va verificato (una consegna davvero lunghissima, o un errore?) prima di scegliere fra media e mediana per descrivere il ritardo tipico. I numeri sono stati controllati con np.percentile.

Che cosa sai fare adesso

Sai cominciare da forma, tipi e valori mancanti, descrivere una variabile con media e mediana affiancate e segnalare i valori fuori recinzione senza eliminarli per riflesso. Sai leggere due correlazioni diverse insieme, ricordando che nessuna delle due dice la causa, sai riconoscere un paradosso di Simpson spezzando i dati per gruppo e sai scegliere il grafico in base alla domanda. Con l'EDA fatta, il resto della serie diventa lo strumento giusto al momento giusto: gli stimatori per riassumere, i test di ipotesi per decidere, la regressione per modellare, il bootstrap per misurare l'incertezza. Per l'uso di pandas e matplotlib in dettaglio vedi Pandas: tabelle di dati senza sorprese e Visualizzare i dati con Matplotlib.

La serie «Probabilità e Statistica» è completa. I quindici articoli, nell'ordine: spazi di probabilità e assiomi di Kolmogorov, variabili discrete, variabili continue, distribuzioni multivariate e covarianza, teorema centrale del limite, stimatori, test di ipotesi, ANOVA, chi-quadro, regressione lineare, stima bayesiana, catene di Markov, statistica non parametrica, bootstrap e questo.