Analisi esplorativa e data dredging: 3 esercizi svolti in Python

È un problema che riguarda chiunque lavori con i dati. Un team marketing confronta cinque landing page e proclama una vincitrice, un direttore commerciale individua il negozio con la crescita più alta e vuole replicarne il modello, mentre un reparto ospedaliero misura decine di biomarcatori e seleziona quelli con i p-value più piccoli per progettare un nuovo protocollo.

In tutti e tre i casi l’analisi esplorativa è utile, perché consente di individuare segnali, formulare spiegazioni plausibili e decidere quali ipotesi meritino ulteriori verifiche. Il rischio nasce quando un risultato emerso durante l’esplorazione viene presentato come se fosse stato confermato da un test indipendente.

Questo problema viene spesso chiamato data dredging, soprattutto quando si cercano sistematicamente risultati significativi fra molte analisi possibili. È collegato al problema dei confronti multipli, alla winner’s curse, alla regressione verso la media e alla selezione dei risultati più favorevoli. Può manifestarsi anche quando si osservano ripetutamente i dati durante un esperimento oppure quando si scelgono, dopo averli esaminati, le variabili, i segmenti e le specificazioni del modello da riportare.

I tre esercizi seguenti affrontano questi meccanismi con difficoltà crescente, attraverso esempi aziendali e clinici, formule statistiche, codice Python e interpretazioni operative. L’obiettivo non è scoraggiare l’esplorazione, ma distinguere un risultato che suggerisce una buona ipotesi da un risultato che offre evidenze sufficienti per prendere una decisione.

Esercizio 1 — Facile: le cinque landing page di Click&Co

Scenario

Click&Co. è un e-commerce di biciclette che vuole migliorare il tasso di conversione delle proprie pagine promozionali. Il team growth prepara cinque versioni della landing page, indicate con A, B, C, D ed E, assegnando casualmente 1.000 visitatori a ciascuna variante. La versione A è il controllo, mentre le altre quattro contengono modifiche al layout, ai messaggi commerciali e alla posizione del pulsante di acquisto.

Dopo una settimana, il marketer confronta ogni nuova versione con A attraverso quattro test z per due proporzioni. Il confronto con D restituisce un p-value di circa 0,036 e il team annuncia su Slack che D è la vincitrice, proponendo di trasferirle tutto il traffico.

Il risultato sembra convincente, perché la conversione di D è del 13%, contro il 10% del controllo, ma la decisione trascura un aspetto fondamentale: il team ha cercato un risultato favorevole fra quattro confronti diversi.

Dati osservati

Variante	Visitatori	Conversioni	Tasso di conversione
A — controllo	1.000	100	10,0%
B	1.000	120	12,0%
C	1.000	95	9,5%
D	1.000	130	13,0%
E	1.000	105	10,5%
Totale	5.000	550	11,0%

Calcolo del tasso complessivo e confronto D contro A

Il tasso di conversione complessivo è dato dal rapporto tra tutte le conversioni e tutti i visitatori:

[math]\displaystyle \begin{aligned}
\hat p=\frac{550}{5000}=0{,}11=11%.
\end{aligned}[/math]

Per confrontare D con A utilizziamo la proporzione pooled, cioè la proporzione stimata combinando i due gruppi sotto l’ipotesi nulla di uguaglianza:

[math]\displaystyle \begin{aligned}
\hat p_{DA}=\frac{130+100}{1000+1000}=0{,}115.
\end{aligned}[/math]

L’errore standard della differenza fra le proporzioni è:

[math]\displaystyle \begin{aligned}
SE=\sqrt{\hat p_{DA}(1-\hat p_{DA})
\left(\frac1{1000}+\frac1{1000}\right)}.
\end{aligned}[/math]

Sostituendo i valori:

[math]\displaystyle \begin{aligned}
SE=\sqrt{0{,}115\cdot0{,}885\cdot0{,}002}
\approx0{,}01427.
\end{aligned}[/math]

La statistica z vale quindi:

[math]\displaystyle \begin{aligned}
z=\frac{0{,}13-0{,}10}{0{,}01427}
\approx2{,}10.
\end{aligned}[/math]

Per un test bilaterale, il p-value è circa 0,036.

Se questo fosse stato l’unico confronto definito prima di osservare i dati, il risultato sarebbe significativo al livello nominale del 5%, assumendo che il test e le condizioni del disegno sperimentale siano appropriati. Il problema non è quindi il calcolo, ma il contesto nel quale il risultato viene interpretato.

Il rischio dei confronti multipli

Quando tutti e quattro i confronti vengono effettuati al livello [math]\alpha=0{,}05[/math], aumenta la probabilità di ottenere almeno un risultato apparentemente significativo anche se tutte le varianti hanno la stessa efficacia reale.

Se i test fossero indipendenti e tutte le ipotesi nulle fossero vere, il family-wise error rate, o FWER, sarebbe:

[math]\displaystyle \begin{aligned}
FWER&=1-(1-\alpha)^k\
&=1-0{,}95^4\
&\approx0{,}1855.
\end{aligned}[/math]

Si tratta di circa il 18,5% di probabilità di almeno un falso positivo, anziché del 5% associato a un singolo test.

Questa formula è esatta sotto l’indipendenza, ma i quattro confronti di Click&Co. condividono il controllo A e non sono indipendenti. Il valore del 18,5% è quindi un riferimento teorico, non la probabilità esatta di falso positivo in questo esperimento. Il principio resta valido: eseguire più test senza tenere conto della molteplicità può aumentare il rischio di conclusioni errate.

Correzione di Holm e Bonferroni

Un metodo semplice per controllare il FWER è la correzione di Bonferroni, che divide il livello di significatività per il numero di confronti:

[math]\displaystyle \begin{aligned}
\alpha’=\frac{0{,}05}{4}=0{,}0125.
\end{aligned}[/math]

Il confronto D contro A non supera questa soglia, perché [math]0{,}036>0{,}0125[/math].

Un’alternativa spesso preferibile è il metodo di Holm, che ordina i p-value dal più piccolo al più grande e applica soglie progressivamente meno restrittive. Holm controlla il FWER senza richiedere indipendenza e risulta almeno altrettanto potente di Bonferroni per la stessa famiglia di test.

Con i dati di Click&Co., il p-value di D corretto con Holm è circa 0,14, quindi la variante non risulta significativa dopo la correzione.

Questo non dimostra che D sia uguale alle altre pagine, ma indica che i dati disponibili non consentono di dichiararla vincitrice con il livello di controllo degli errori scelto.

Il test omnibus: esiste una differenza complessiva?

Prima di interpretare i singoli confronti, può essere utile valutare l’ipotesi globale secondo cui tutte le varianti hanno la stessa probabilità di conversione. Il test del chi-quadrato di Pearson utilizza una tabella di contingenza con cinque righe, una per variante, e due colonne, una per le conversioni e l’altra per le non-conversioni.

Sotto l’ipotesi nulla, ogni variante ha 110 conversioni attese e 890 non-conversioni attese. La statistica è:

[math]\displaystyle \begin{aligned}
\chi^2=\sum_{i=1}^{5}
\left[
\frac{(O_{i,1}-E_{i,1})^2}{E_{i,1}}
+
\frac{(O_{i,2}-E_{i,2})^2}{E_{i,2}}
\right].
\end{aligned}[/math]

Poiché il totale delle conversioni è 550, la somma dei quadrati degli scarti delle conversioni è 850; anche la somma dei quadrati degli scarti delle non-conversioni è 850. Pertanto:

[math]\displaystyle \begin{aligned}
\chi^2
&=850\left(\frac1{110}+\frac1{890}\right)\
&\approx8{,}60.
\end{aligned}[/math]

I gradi di libertà sono [math](5-1)(2-1)=4[/math], e il p-value è circa 0,072. Il test globale non è significativo al 5%.

Il risultato non è in contraddizione con il p-value grezzo di D contro A, perché i due test rispondono a domande diverse. Il test omnibus valuta se esiste una differenza complessiva fra le cinque varianti, mentre il test z considera uno specifico confronto. Inoltre, il confronto grezzo è stato selezionato dopo aver esaminato quattro risultati, circostanza che ne modifica l’interpretazione inferenziale.

Un test omnibus non significativo non dimostra che tutte le varianti siano uguali e non è una condizione universale che debba precedere ogni confronto pianificato. La strategia inferenziale va definita in base alla domanda iniziale, al disegno sperimentale e alla famiglia di ipotesi che si intende controllare.

Implementazione Python

Il codice seguente calcola il tasso complessivo, i quattro test z, i p-value corretti con Holm e il test globale di Pearson, includendo anche le frequenze attese per verificare che l’approssimazione sia ragionevole.

import numpy as np
import pandas as pd
from scipy.stats import chi2_contingency
from statsmodels.stats.proportion import proportions_ztest
from statsmodels.stats.multitest import multipletests

# Dati osservati per ciascuna variante.
conversioni = {"A": 100, "B": 120, "C": 95, "D": 130, "E": 105}
visitatori = {nome: 1000 for nome in conversioni}

alpha = 0.05
controllo = "A"

# Tasso complessivo.
tot_conv = sum(conversioni.values())
tot_vis = sum(visitatori.values())
print(f"Tasso complessivo: {tot_conv}/{tot_vis} = {tot_conv/tot_vis:.2%}")

# Confronti di ciascuna variante con il controllo.
risultati = []
pvals = []

for nome in conversioni:
    if nome == controllo:
        continue

    z, p = proportions_ztest(
        count=[conversioni[nome], conversioni[controllo]],
        nobs=[visitatori[nome], visitatori[controllo]],
        alternative="two-sided"
    )

    risultati.append({
        "Variante": nome,
        "Tasso": conversioni[nome] / visitatori[nome],
        "z": z,
        "p grezzo": p
    })
    pvals.append(p)

# Correzione di Holm per il controllo del FWER.
rifiuta, p_holm, _, _ = multipletests(
    pvals, alpha=alpha, method="holm"
)

for i, r in enumerate(risultati):
    r["p Holm"] = p_holm[i]
    r["Significativo Holm"] = rifiuta[i]

tabella = pd.DataFrame(risultati)
print("\nConfronti contro A:")
print(tabella.round(4).to_string(index=False))

# Tabella 5 x 2: conversioni e non-conversioni.
osservate = np.array([
    [conversioni[nome], visitatori[nome] - conversioni[nome]]
    for nome in conversioni
])

chi2, p_globale, df, attese = chi2_contingency(
    osservate, correction=False
)

print("\nTest omnibus di Pearson:")
print(f"Chi-quadrato = {chi2:.4f}")
print(f"Gradi di libertà = {df}")
print(f"p-value = {p_globale:.4f}")
print(f"Frequenza attesa minima = {attese.min():.1f}")
print(f"Significativo al 5%: {p_globale < alpha}")

L’output mostra un tasso complessivo dell’11%, un p-value grezzo di circa 0,036 per D contro A e un p-value corretto con Holm di circa 0,14. Il test omnibus restituisce un chi-quadrato vicino a 8,60, con quattro gradi di libertà e un p-value di circa 0,072.

Dal punto di vista gestionale, la decisione corretta non consiste nel trasferire automaticamente tutto il traffico a D, ma nel riconoscere che il risultato è promettente e merita una verifica con un disegno confermativo adeguato. Se il confronto D-A era stato pianificato prima dell’esperimento, va valutato secondo quel piano; se D è stata scelta dopo aver confrontato tutte le varianti, occorre considerare la selezione effettuata.

Domande di riflessione

Perché Bonferroni funziona anche senza indipendenza?

La disuguaglianza di Boole stabilisce che la probabilità dell’unione di eventi non può superare la somma delle loro probabilità. Se ogni test ha probabilità di falso positivo non superiore a [math]\alpha/k[/math], la probabilità di almeno un falso positivo nella famiglia non supera [math]\alpha[/math], qualunque sia la struttura di dipendenza.

Cosa cambia se D contro A era l’unico confronto pianificato?

Cambia la famiglia di ipotesi che si intende sottoporre a inferenza. Un confronto realmente definito prima di osservare i risultati può essere valutato al livello nominale previsto, mentre scegliere D fra quattro confronti dopo aver visto i p-value introduce una selezione che deve essere gestita.

Cosa succede se si confrontano tutte le varianti fra loro?

Con cinque varianti si ottengono dieci coppie. Sotto indipendenza e nullo globale, il FWER teorico non corretto diventa [math]1-0{,}95^{10}\approx40{,}1\%[/math]. I confronti sono però dipendenti, quindi questa formula non è esatta per il disegno reale. Per confronti tutti contro tutti si possono usare procedure specifiche, come Tukey HSD nei contesti in cui sono appropriate, oppure altre correzioni coerenti con il disegno.

Esercizio 2 — Facile/medio: il negozio star di FrescoMarket

Scenario

FrescoMarket gestisce venti supermercati e registra una crescita media annua del 2%. Il direttore commerciale osserva che il punto vendita di Via Verdi è cresciuto del 18% e conclude che il negozio abbia sviluppato una strategia particolarmente efficace, proponendo di replicarne il modello nell’intera catena.

Il dato è reale nel senso che rappresenta la crescita osservata, ma non è sufficiente per stabilire quanto del risultato dipenda da un vantaggio persistente e quanto dalla variabilità delle misure. Quando si seleziona il negozio con il risultato più alto, infatti, si sceglie anche quello che ha maggiori probabilità di aver beneficiato di un errore casuale favorevole.

Il modello statistico

Assumiamo che la crescita vera del negozio [math]i[/math] sia descritta da:

[math]\displaystyle \begin{aligned}
\theta_i\sim N(\mu,\tau^2),
\end{aligned}[/math]

con media di catena [math]\mu=2\%[/math] e deviazione standard fra negozi [math]\tau=3[/math] punti percentuali.

La crescita osservata è:

[math]\displaystyle \begin{aligned}
y_i=\theta_i+\varepsilon_i,
\qquad
\varepsilon_i\sim N(0,\sigma^2),
\end{aligned}[/math]

dove [math]\sigma=6[/math] punti percentuali rappresenta il rumore anno-su-anno. In questo modello, la crescita vera varia fra negozi e la crescita osservata contiene anche una componente casuale.

Quanto può essere alto il massimo per puro caso?

Consideriamo inizialmente uno scenario nel quale tutti i negozi abbiano la stessa crescita vera, pari al 2%, mentre le osservazioni differiscano per effetto del rumore. In questo caso le crescite osservate seguono una distribuzione normale con media 2% e deviazione standard 6 punti percentuali.

Il valore atteso del massimo di venti variabili normali standard indipendenti è circa 1,87. Di conseguenza, il massimo osservato atteso in questo scenario è:

[math]\displaystyle \begin{aligned}
E[\max(y_1,\ldots,y_{20})]
&\approx2%+1{,}87\cdot6%\
&\approx13{,}2%.
\end{aligned}[/math]

Una crescita osservata del 18% è superiore a questo massimo atteso, ma non è sufficiente per concludere che Via Verdi possieda un vantaggio strutturale. Il valore 13,2% è una media teorica del massimo sotto un’ipotesi specifica, non una soglia oltre la quale un risultato diventa automaticamente significativo.

La winner’s curse nasce proprio da questo meccanismo: selezionando il massimo fra molte stime rumorose, si tende a sovrastimare il valore che si osserverebbe di nuovo. Il negozio può avere un vantaggio autentico, ma il risultato osservato tende a essere più estremo della sua prestazione attesa.

Shrinkage: quanto è plausibile la crescita vera di Via Verdi?

Nel modello normal-normal, se i parametri della distribuzione iniziale e del rumore sono noti, la media condizionale della crescita vera è:

[math]\displaystyle \begin{aligned}
E[\theta_i\mid y_i]
=\mu+B(y_i-\mu),
\end{aligned}[/math]

dove il fattore di shrinkage è:

[math]\displaystyle \begin{aligned}
B=\frac{\tau^2}{\tau^2+\sigma^2}.
\end{aligned}[/math]

Con [math]\tau=3[/math] e [math]\sigma=6[/math], si ottiene:

[math]\displaystyle \begin{aligned}
B=\frac{3^2}{3^2+6^2}
=\frac9{45}=0{,}20.
\end{aligned}[/math]

Per Via Verdi, la stima diventa:

[math]\displaystyle \begin{aligned}
\hat\theta
&=2%+0{,}20(18%-2%)\
&=5{,}2%.
\end{aligned}[/math]

La stima condizionale della crescita vera è quindi 5,2%, molto inferiore al 18% osservato, perché il modello attribuisce un peso limitato allo scarto individuale quando il rumore è elevato rispetto all’eterogeneità fra negozi.

Questo risultato non significa che sia stato dimostrato che l’80% della crescita osservata è rumore. Significa che, date le ipotesi del modello e i parametri forniti, la stima della crescita vera combina il dato individuale con la media della popolazione, riducendo l’estremità delle osservazioni più rumorose.

 

Che cosa aspettarsi l’anno successivo?

Nel modello specificato, il valore atteso della crescita vera di Via Verdi, condizionatamente al dato osservato, è 5,2%. Anche la crescita osservata del prossimo anno ha questo valore come previsione media, ma sarà soggetta a ulteriore variabilità.

La distribuzione predittiva per una nuova osservazione include sia l’incertezza sulla crescita vera sia il nuovo rumore annuale. Con i parametri noti dell’esempio, la varianza predittiva è:

[math]\displaystyle \begin{aligned}
\operatorname{Var}(y_{\text{nuovo}}\mid y)
=\sigma^2+\frac{\tau^2\sigma^2}{\tau^2+\sigma^2}
=36+7{,}2=43{,}2.
\end{aligned}[/math]

La deviazione standard predittiva è quindi circa 6,57 punti percentuali. Questa dispersione è importante: il 5,2% non è un obiettivo garantito, ma una previsione centrale sotto il modello adottato.

Il direttore commerciale dovrebbe quindi verificare quali fattori possano spiegare la performance di Via Verdi, come la composizione della clientela, le promozioni, la disponibilità dei prodotti, la concorrenza locale e le condizioni iniziali del negozio. La replica del modello dovrebbe essere valutata attraverso un piano che consenta di distinguere gli effetti persistenti dalle oscillazioni casuali.

Implementazione Python

Il codice applica direttamente il modello normal-normal con i parametri dichiarati, senza stimare il rumore dalla stessa variabile che dovrebbe essere corretta. Questa distinzione è essenziale: se si sottrae la varianza campionaria di un vettore dalla stessa varianza campionaria per stimare l’eterogeneità, si rischia di ottenere una stima nulla per costruzione.

import numpy as np

# Parametri del modello, espressi in punti percentuali.
mu = 2.0
tau = 3.0
sigma = 6.0
crescita_osservata = 18.0

# Shrinkage normal-normal.
var_eterogeneita = tau**2
var_rumore = sigma**2

B = var_eterogeneita / (var_eterogeneita + var_rumore)
stima_vera = mu + B * (crescita_osservata - mu)

# Massimo atteso di 20 osservazioni normali standard.
n_negozi = 20
massimo_standard_atteso = 1.87
massimo_atteso = mu + sigma * massimo_standard_atteso

# Varianza e deviazione standard della previsione di una nuova osservazione.
var_posteriore = (
    var_eterogeneita * var_rumore
    / (var_eterogeneita + var_rumore)
)
var_predittiva = var_posteriore + var_rumore
sd_predittiva = np.sqrt(var_predittiva)

print(f"Massimo atteso sotto crescita vera comune: {massimo_atteso:.2f}%")
print(f"Fattore di shrinkage: {B:.2f}")
print(f"Crescita osservata: {crescita_osservata:.2f}%")
print(f"Stima della crescita vera: {stima_vera:.2f}%")
print(f"Deviazione standard predittiva: {sd_predittiva:.2f} punti percentuali")


def shrink_normal_normal(y, mu, sigma, tau):
    """
    Calcola le stime normal-normal quando mu, sigma e tau
    sono noti o stimati separatamente con un metodo appropriato.

    y: osservazioni, in punti percentuali
    mu: media della popolazione
    sigma: deviazione standard del rumore osservativo
    tau: deviazione standard dell'eterogeneità reale
    """
    y = np.asarray(y, dtype=float)

    if sigma < 0 or tau < 0:
        raise ValueError("sigma e tau devono essere non negativi.")

    denominatore = tau**2 + sigma**2
    B = tau**2 / denominatore if denominatore > 0 else 0.0

    stime = mu + B * (y - mu)
    return stime, B


stima_verificata, _ = shrink_normal_normal(
    [crescita_osservata], mu=mu, sigma=sigma, tau=tau
)

print(f"Verifica funzione: {stima_verificata[0]:.2f}%")

L’output atteso, con gli arrotondamenti indicati, è:

Massimo atteso sotto crescita vera comune: 13.22%
Fattore di shrinkage: 0.20
Crescita osservata: 18.00%
Stima della crescita vera: 5.20%
Deviazione standard predittiva: 6.57 punti percentuali
Verifica funzione: 5.20%

La funzione è riutilizzabile, ma non sostituisce la stima di un modello gerarchico completo quando i parametri della popolazione sono ignoti. Per un’analisi aziendale reale sarebbe opportuno stimare l’eterogeneità sui dati di tutti i negozi, verificare le ipotesi del modello e valutare la previsione su un periodo successivo.

Domande di riflessione

Da dove deriva la formula dello shrinkage?

Deriva dalla coniugatezza della distribuzione normale nel modello normal-normal, nel quale la distribuzione iniziale della crescita vera e la distribuzione dell’osservazione condizionata alla crescita vera sono entrambe normali. La media condizionale risultante è una combinazione pesata della media della popolazione e del valore osservato.

Che cosa succede al fattore di shrinkage se aumenta il rumore?

Se [math]\sigma[/math] aumenta, il fattore [math]B[/math] diminuisce e la stima si avvicina alla media della popolazione, perché l’osservazione individuale diventa meno informativa. Se [math]\tau[/math] tende a zero, il modello assume che i negozi siano sostanzialmente identici nella crescita vera e la stima converge alla media comune.

Perché essere il negozio con la crescita più alta non dimostra di avere la strategia migliore?

Essere il massimo è una proprietà della classifica osservata, mentre possedere il massimo valore vero è una proprietà diversa, che dipende anche dal rumore e dall’incertezza sulle stime. Selezionare il risultato più alto garantisce che qualcuno sia primo, ma non garantisce che il suo vantaggio sia persistente o replicabile.

Esercizio 3 — Medio: i 30 biomarcatori dell’Ospedale San Luca

Scenario

Un reparto di medicina interna analizza 30 biomarcatori su 400 pazienti dimessi, con l’obiettivo esplorativo di individuare possibili associazioni con il rientro in ospedale entro 30 giorni. Il data analyst esegue 30 regressioni logistiche univariate e presenta al primario i cinque p-value più piccoli, proponendo di utilizzare due biomarcatori per costruire un nuovo protocollo.

Lo screening è un modo legittimo per generare ipotesi, ma la selezione dei risultati più favorevoli richiede un metodo statistico coerente con lo scopo dell’analisi. Inoltre, una regressione univariata non dimostra che un biomarcatore abbia un’associazione indipendente dagli altri fattori clinici, né stabilisce che il biomarcatore sia utile per una decisione terapeutica.

I p-value osservati

Biomarcatore	p-value
M1	0,0008
M2	0,0025
M3	0,0120
M4	0,0300
M5	0,0400

Questi sono i cinque valori più piccoli dei 30 test. I restanti 25 p-value non sono disponibili nello scenario e non devono essere inventati per completare artificialmente il calcolo.

Quanti falsi positivi ci si può aspettare?

Se tutte le 30 ipotesi nulle fossero vere e ogni test fosse condotto al livello [math]\alpha=0{,}05[/math], il numero atteso di falsi positivi sarebbe:

[math]\displaystyle \begin{aligned}
E[V]=m\alpha=30\cdot0{,}05=1{,}5.
\end{aligned}[/math]

Sotto l’ulteriore ipotesi di indipendenza dei test, la probabilità di ottenere almeno un falso positivo sarebbe:

[math]\displaystyle \begin{aligned}
P(V\ge1)
&=1-(1-\alpha)^{30}\
&=1-0{,}95^{30}\
&\approx0{,}785.
\end{aligned}[/math]

Il valore, circa il 78,5%, è valido sotto indipendenza e nullo globale. Se i test sono correlati, la formula non è in generale esatta: la probabilità dipende dalla struttura della dipendenza. Il numero atteso di falsi positivi resta invece pari a [math]m\alpha[/math] quando ciascun test ha probabilità di falso positivo [math]\alpha[/math], anche senza indipendenza.

Il messaggio operativo è chiaro: quando si esaminano molti biomarcatori, è plausibile trovare associazioni apparentemente interessanti anche in assenza di effetti reali. La correzione statistica non elimina la necessità di una validazione indipendente.

Bonferroni: controllare il rischio di almeno un falso positivo

La correzione di Bonferroni controlla il family-wise error rate dividendo il livello di significatività per il numero totale di ipotesi:

[math]\displaystyle \begin{aligned}
\alpha’=\frac{0{,}05}{30}\approx0{,}001667.
\end{aligned}[/math]

Fra i cinque p-value disponibili, soltanto M1 supera questo criterio, perché [math]0{,}0008<0{,}001667[/math]. M2, con p-value pari a 0,0025, non raggiunge la soglia Bonferroni.

Il vantaggio di questo approccio è che limita al 5% la probabilità di almeno un falso rifiuto nella famiglia, sotto la validità dei test individuali, senza richiedere indipendenza. Il costo può essere una minore potenza, soprattutto quando si analizzano molti marcatori correlati e gli effetti reali sono modesti.

Benjamini–Hochberg: controllare il false discovery rate

Il metodo Benjamini–Hochberg, indicato con BH, affronta una domanda diversa. Invece di controllare la probabilità di almeno un falso positivo, controlla la proporzione attesa di falsi positivi fra le ipotesi rifiutate, nelle condizioni di validità del metodo.

Ordinando i p-value dal più piccolo al più grande, il criterio confronta ogni valore [math]p_{(k)}[/math] con la soglia:

[math]\displaystyle \begin{aligned}
p_{(k)}\le\frac{k}{m}q,
\end{aligned}[/math]

dove [math]q=0{,}05[/math] è il livello desiderato per il controllo del FDR.

Per i primi cinque ranghi si ottiene:

Rango [math]k[/math]	p-value	Soglia BH
1	0,0008	0,001667
2	0,0025	0,003333
3	0,0120	0,005000
4	0,0300	0,006667
5	0,0400	0,008333

M1 e M2 superano le rispettive soglie. Questo permette di concludere che entrambi rientrano nell’insieme dei risultati rifiutati dal procedimento BH, purché i valori riportati siano davvero i cinque p-value più piccoli dell’intera famiglia di 30 test.

Non è invece possibile concludere che siano gli unici due risultati rifiutati, perché la decisione definitiva dipende anche dagli altri 25 p-value. Se uno dei ranghi successivi soddisfacesse il criterio BH, potrebbero essere selezionati anche altri biomarcatori.

La distinzione è importante: la tabella permette una conclusione parziale sui primi due marcatori, ma non una ricostruzione completa dei risultati aggiustati né una decisione definitiva sull’intera famiglia.

Implementazione Python

Il codice seguente mostra la soglia Bonferroni e verifica le soglie BH associate ai primi cinque ranghi. La funzione per il risultato definitivo richiede tutti i 30 p-value, in modo da evitare di correggere una famiglia di cinque test quando l’analisi originaria ne comprende trenta.

import numpy as np
import pandas as pd
from statsmodels.stats.multitest import multipletests

# I cinque p-value più piccoli riportati nello scenario.
p_noti = {
    "M1": 0.0008,
    "M2": 0.0025,
    "M3": 0.0120,
    "M4": 0.0300,
    "M5": 0.0400
}

m = 30
alpha = 0.05

# Bonferroni: il denominatore è il numero totale di test.
soglia_bonferroni = alpha / m
print(f"Soglia Bonferroni: {soglia_bonferroni:.6f}")

for nome, p in p_noti.items():
    significativo = p <= soglia_bonferroni
    print(
        f"{nome}: p={p:.4f} -> "
        f"{'significativo' if significativo else 'non significativo'}"
    )

# Verifica BH parziale, limitata ai cinque valori noti.
tabella = pd.DataFrame({
    "Marcatore": list(p_noti.keys()),
    "p": list(p_noti.values())
})

tabella["Rango"] = np.arange(1, len(tabella) + 1)
tabella["Soglia BH"] = tabella["Rango"] / m * alpha
tabella["Supera la soglia al rango"] = (
    tabella["p"] <= tabella["Soglia BH"]
)

print("\nVerifica parziale BH:")
print(tabella.round(6).to_string(index=False))


def analizza_fdr(pvals, alpha=0.05):
    """
    Calcola p-value aggiustati BH e decisioni per l'intera famiglia.
    pvals deve contenere tutti i test previsti dal piano di analisi.
    """
    pvals = np.asarray(pvals, dtype=float)

    if pvals.ndim != 1 or len(pvals) == 0:
        raise ValueError("Fornire un vettore non vuoto di p-value.")

    if not np.all(np.isfinite(pvals)):
        raise ValueError("Tutti i p-value devono essere finiti.")

    if np.any((pvals < 0) | (pvals > 1)):
        raise ValueError("I p-value devono appartenere all'intervallo [0, 1].")

    rifiuta, p_aggiustati, _, _ = multipletests(
        pvals, alpha=alpha, method="fdr_bh"
    )

    return pd.DataFrame({
        "p grezzo": pvals,
        "p aggiustato BH": p_aggiustati,
        "Rifiuta H0": rifiuta
    }).sort_values("p grezzo").reset_index(drop=True)


# Per ottenere il risultato definitivo, utilizzare i 30 p-value originali:
#
# pvals_completi = [p1, p2, ..., p30]
# risultati_bh = analizza_fdr(pvals_completi, alpha=0.05)
# print(risultati_bh.to_string(index=False))

L’output parziale mostra una soglia Bonferroni di circa 0,001667, con M1 significativo e M2–M5 non significativi secondo questo criterio. La tabella BH evidenzia invece che M1 e M2 superano le rispettive soglie ai ranghi 1 e 2.

Per il risultato definitivo occorre eseguire la funzione sui 30 p-value originari. Il codice non attribuisce valori fittizi ai test mancanti e non confonde una verifica parziale con la decisione finale sull’intera famiglia.

Come spiegare FWER e FDR al primario

Il FWER riguarda il rischio di commettere almeno un errore nella famiglia di test, mentre il FDR riguarda la proporzione attesa di errori fra i risultati dichiarati significativi. Il primo è spesso preferibile quando un singolo falso positivo può avere conseguenze rilevanti; il secondo è utile negli screening esplorativi, quando si vuole costruire una lista di candidati da verificare in seguito.

Nessuno dei due metodi dimostra da solo che un biomarcatore sia clinicamente utile. I risultati devono essere valutati considerando la dimensione dell’effetto, gli intervalli d’incertezza, la qualità delle misure, i fattori confondenti e la capacità di predire il rientro ospedaliero in una coorte indipendente.

Domande di riflessione

Perché la formula della probabilità di almeno un falso positivo richiede indipendenza?

La formula [math]1-(1-\alpha)^m[/math] deriva dal calcolo della probabilità che nessun test produca un falso positivo, che sotto indipendenza è il prodotto delle probabilità dei singoli eventi. Con test correlati il prodotto non è generalmente valido, perché la probabilità di errori simultanei dipende dalla struttura della dipendenza.

Se i biomarcatori sono fortemente correlati, il numero effettivo di test è sempre 30?

Il numero di ipotesi testate è 30, ma la quantità di informazione indipendente contenuta nei marcatori può essere inferiore. La correlazione non autorizza però a ridurre arbitrariamente il numero di test nel calcolo della correzione. Eventuali metodi basati sul numero effettivo di test devono essere giustificati e adeguati al disegno; in alternativa, si possono utilizzare procedure di resampling o modelli che tengano conto della struttura dei dati.

Perché BH può rifiutare M2 mentre Bonferroni no?

I due metodi controllano errori diversi. Bonferroni limita la probabilità di almeno un falso positivo nella famiglia, mentre BH controlla il FDR nelle condizioni previste dal metodo. Per questo BH può individuare più candidati, circostanza utile nella ricerca esplorativa, purché la selezione sia seguita da una validazione indipendente.

Peeking, garden of forking paths e selezione fra migliaia di strategie

I confronti multipli non sono l’unica forma di data dredging. Un rischio analogo compare quando si osservano ripetutamente i risultati di un esperimento e si decide di interromperlo nel momento in cui il p-value scende sotto 0,05, senza aver previsto una procedura di analisi sequenziale. Questa pratica, chiamata peeking, altera le proprietà del test e può aumentare il rischio di falsi positivi.

Il garden of forking paths descrive invece una situazione nella quale il ricercatore prende molte decisioni analitiche possibili, come scegliere le variabili da includere, escludere osservazioni, provare trasformazioni differenti, selezionare segmenti o modificare le soglie, e alla fine riporta soltanto il percorso che produce il risultato più favorevole. Anche se ogni singola scelta può apparire ragionevole, la combinazione delle alternative esplorate rende fuorviante interpretare il p-value finale come se l’intero percorso fosse stato definito in anticipo.

Un problema simile riguarda la selezione fra migliaia di strategie di trading, campagne commerciali o modelli predittivi. Quando si sceglie la strategia con la migliore performance storica, il risultato può essere influenzato dalla selezione e non riflettere un vantaggio stabile. Una valutazione affidabile richiede una separazione rigorosa fra dati utilizzati per sviluppare le strategie e dati utilizzati per valutarle, oltre a procedure che tengano conto della molteplicità delle scelte effettuate.

In questi contesti non basta dividere meccanicamente il livello di significatività per il numero dei test dichiarati, perché la famiglia di analisi può essere difficile da ricostruire dopo il fatto. È più utile documentare le decisioni analitiche, definire in anticipo gli obiettivi confermativi, utilizzare procedure appropriate per gli esperimenti sequenziali e verificare i risultati su dati che non hanno partecipato alla selezione.

Dall’esplorazione alla conferma: la checklist dell’analista

Prima di trasformare un risultato esplorativo in una decisione, è opportuno verificare cinque aspetti.

• Quante analisi sono state provate? Occorre considerare variabili, gruppi, segmenti, finestre temporali, modelli e soglie esaminate, ricordando che il numero dei tentativi non coincide sempre con il numero dei test riportati.
• Quale domanda si sta testando? Un’ipotesi definita prima di osservare i dati ha uno status diverso da una spiegazione selezionata dopo aver esaminato i risultati, e questa distinzione deve essere esplicita.
• Quale errore si vuole controllare? Il FWER è adatto quando interessa limitare il rischio di almeno un falso positivo nella famiglia, mentre il FDR è utile quando l’obiettivo è selezionare candidati per ulteriori verifiche.
• Quanto è stabile la stima? È necessario valutare la dimensione dell’effetto, gli intervalli d’incertezza, il rumore di misura, la selezione dei risultati e la possibile regressione verso la media.
• Come verrà verificata l’ipotesi? Prima di procedere è opportuno definire la metrica, il criterio decisionale, la finestra temporale e il metodo di validazione su dati indipendenti, quando disponibili.

L’analisi esplorativa non è un errore da evitare, ma una componente essenziale del lavoro sui dati. Il problema nasce quando un segnale individuato durante la ricerca viene presentato come una conferma, senza considerare quanti tentativi sono stati effettuati e quale incertezza accompagna il risultato.

Un buon analista non è quello che trova sempre un p-value significativo, ma quello che sa spiegare quanto è solida l’evidenza, quali alternative sono state considerate e quali verifiche servono prima di trasformare una scoperta in una decisione aziendale o clinica.

Conclusione

I tre casi mostrano come il data dredging possa assumere forme diverse: un confronto selezionato fra molte varianti, una performance estrema scelta da una classifica rumorosa e uno screening di numerose ipotesi cliniche. In ciascuna situazione, l’analisi esplorativa ha un ruolo legittimo, purché i risultati siano interpretati in modo coerente con il processo che li ha generati.

La correzione per i confronti multipli, lo shrinkage, la documentazione delle scelte analitiche e la validazione indipendente rispondono a problemi diversi e non sono strumenti intercambiabili. La scelta del metodo deve dipendere dalla domanda, dal disegno dei dati e dalle conseguenze di una decisione errata.

La regola pratica è semplice: un risultato esplorativo può indicare dove cercare, ma la sua capacità di sostenere una decisione deve essere valutata con procedure che tengano conto della selezione, dell’incertezza e della possibilità di replicare l’effetto.