Test multipli e analisi di sottogruppo (cosa c’è di sbagliato nel trovare sempre qualcosa di significativo)

Introduzione

Nella precedente rubrica avevamo discusso come riportare alcune stime con gli intervalli di confidenza nel riassunto di un lavoro o nella sintesi di un rapporto (in inglese l’abstract di un articolo e l’executive summary di un rapporto tecnico).¹ Rispetto alla procedura di calcolo degli intervalli di confidenza abbiamo mostrato come sia necessaria una correzione e come gli intervalli di confidenza in queste situazioni saranno più ampi al fine di tener conto adeguatamente dell’incertezza che è introdotta dalla selezione di alcune stime e del loro intervallo di confidenza dall’insieme di stime di effetto che sono state calcolate e riportate nell’articolo o nel rapporto.

Selezionare un certo numero di stime da un insieme più grande è appunto il compito che ci si prefigge ogni qualvolta si esegue un test d’ipotesi: scegliamo (dichiariamo significativi) alcuni risultati dall’insieme dei risultati ottenuti. La procedura originariamente proposta da Fisher si basa sulla valutazione della probabilità di errore di primo tipo, o valore di p (p-valore, in inglese p-value), assumendo valida l’ipotesi nulla. Dichiariamo un risultato significativo se il p-valore è inferiore a un predefinito livello 15 di probabilità di errore di primo tipo. Controlliamo cioè il rischio di dichiarare significativo un risultato quando in realtà non lo è.

Questo modo di procedere coerentemente all’intendimento di Fisher è esplorativo e soggettivo, perché l’interpretazione di quanto forte sia l’evidenza contro l’ipotesi nulla di un determinato valore di p è lasciata al singolo ricercatore. Non è una procedura rigida, diversamente dall’impostazione di Neyman-Pearson che invece è di tipo decisionale e richiede anche la specificazione del livello di errore di secondo tipo, che consiste nel dichiarare non significativo un risultato quando invece è vera l’ipotesi alternativa. Questo approccio è usato nella fase di programmazione di uno studio e nella determinazione della dimensione campionaria. Nella letteratura biomedica tuttavia è dominante la pratica di dividere i risultati in statisticamente “significativi” e “non significativi” senza considerare l’errore di secondo tipo. Questo porta al paradosso di trascurare risultati clinicamente o epidemiologicamente importanti (come grandezza della misura di effetto) ma non significativi per la piccola dimensione campionaria (gli studi piccoli) e sovrastimare l’importanza di risultati in quanto statisticamente significativi. Tale semplificazione ha portato a considerare equivalenti le posizioni di Fisher e di Neyman-Pearson, e ha favorito l’erronea interpretazione del p-valore come la probabilità che sia vera l’ipotesi, e cioè che un risultato significativo voglia dire che è molto improbabile che l’ipotesi nulla sia vera.

Supponiamo di eseguire 1000 test d’ipotesi e che l’ipotesi nulla sia vera per 900 di questi. Supponiamo anche che il nostro studio abbia una numerosità tale per cui la potenza sia dell’80% per l’effetto che vogliamo dimostrare (la potenza è 1 – β, la probabilità di errore di secondo tipo):

Assumendo un livello di significatività del 5%, su 900 risultati 45 saranno dichiarati significativi quando è vera l’ipotesi nulla mentre ben 80 quando è vera l’ipotesi alternativa. Tuttavia, dei 125 risultati dichiarati statisticamente significativi complessivamente, ben il 36% (45) lo sono falsamente! In tal caso la probabilità a posteriori dell’ipotesi nulla è il 36% mentre il p-valore quantifica solo l’evidenza empirica contro l’ipotesi nulla e non tiene conto dell’ipotesi alternativa. Quanto piccolo deve allora essere un p-valore? Ovviamente dipende dall’ipotesi alternativa ma anche dalla probabilità a priori dell’ipotesi nulla (nell’esempio in tabella 900/1000=90%). Goodman ha proposto una semplice trasformazione del p-valore, senza considerare in particolare un’ipotesi alternativa, e di calcolare la probabilità a posteriori sotto diversi scenari. Il vantaggio della sua proposta è duplice: si apprezza direttamente quanto p-valori di 0.05 possono essere ancora troppo grandi per abbandonare l’ipotesi nulla e soprattutto come sia necessario considerare la probabilità a priori dell’ipotesi nulla e di quanto questa si modifichi a seguito dei risultati osservati. Il calcolo è semplice: dato un p-valore si ottiene il minimo fattore di Bayes (MBF) e la probabilità a posteriori come funzione dell’odds a priori (il rapporto tra probabilità a priori dell’ipotesi nulla diviso il suo complemento a uno):

Nella tabella sono riportati a titolo di esempio quattro studi:

Lo studio C fornisce un’evidenza molto forte contro l’ipotesi nulla. Lo studio D potrebbe non essere così stringente se la probabilità a priori in favore dell’ipotesi nulla fosse molto alta (75%) anche se la nostra aspettativa circa l’ipotesi nulla si sposta di molto (per esempio dal 75% si passa al 12%). Lo studio B invece è abbastanza debole perché la probabilità a posteriori dell’ipotesi nulla resta abbastanza alta anche se la nostra aspettativa a priori non fosse già a svantaggio dell’ipotesi nulla stessa. Analoghe considerazioni per lo studio A, dove siamo ancora più dubbiosi circa l’eventuale abbandono dell’ipotesi nulla.

Quanto detto finora si riferisce alla situazione in cui abbiamo un solo risultato, un solo p-valore, e dobbiamo deciderne l’importanza. Il nostro ragionamento ha dovuto far riferimento a una molteplicità di risultati, come se fossimo nella situazione ipotetica della ripetizione infinita del nostro studio e assumendo una certa probabilità (a priori) dell’ipotesi nulla.

Nella pratica ci troviamo molto spesso di fronte a una molteplicità di risultati, di p-valori. Questo come mostreremo è una ricchezza informativa che permette di esplorare nello stesso spirito di Fisher l’evidenza empirica nei confronti dell’ipotesi nulla in modo coerente con le argomentazioni sopra riportate e con la proposta di Goodman. I confronti multipli non sono un problema ma una risorsa. Ma andiamo per gradi.

Il “problema” dei confronti multipli è molto più frequente di quanto si pensi. Infatti anche nel caso dell’analisi di sottogruppi ci troviamo nella necessità di stimare molte misure di effetto, una per ciascun sottogruppo, e di dover valutare se per qualcuno di questi sottogruppi il risultato sia significativo. Nel lavoro di Sterne e Davey-Smith si suggerisce per inciso che l’analisi di sottogruppo deve essere riportata solo se un test “globale” sia risultato significativo.² Non viene chiarito, ma sembra di capire che ci si riferisca a una procedura di controllo della probabilità di errore di primo tipo globale. Nel prossimo paragrafo formalizzeremo questa impostazione.

Che sia importante è chiaramente discusso dagli autori: troppo discredito arriva all’epidemiologia dal riportare sempre risultati significativi, anche solo per alcuni sottogruppi della popolazione. Cosa c’è di sbagliato nel trovare sempre qualcosa di significativo? Nulla, se si adotta una procedura appropriata di controllo della probabilità di falsi positivi tra tutti i confronti dichiarati significativi.

Test Multipli: qualche tecnicismo

La tabella 1 riporta tutti i possibili esiti di m test di ipotesi. H0 indica l’ipotesi nulla, m è fissato a priori e R, il numero di ipotesi nulle rifiutate, è l’unica variabile casuale osservabile.

Quando siamo nella situazione di eseguire molti test d’ipotesi, cerchiamo di controllare che la probabilità di errore non superi il livello desiderato, pur mantenendo più alta possibile la potenza di ciascuno dei singoli test. La quantità generalmente controllata è il FWER (Family Wise Error Rate), cioè la probabilità di avere almeno un risultato falsamente positivo fra tutti i test di ipotesi eseguiti. La procedura più conosciuta è la correzione di Bonferroni:³ fissato ad α il livello di errore di I tipo per m test di ipotesi, la probabilità di errore di primo tipo per ogni singolo test sarà controllata se adottiamo un livello di significatività pari a α/m. Questa procedura garantisce che la probabilità di avere almeno un risultato falsamente positivo fra tutti gli m test di ipotesi eseguiti sia al massimo pari ad α.

Formalmente, in base alle quantità definite nella tabella 1, possiamo definire FWER = Prob(V >= 1).

Benjamini e Hochberg hanno proposto di spostare l’attenzione al valore atteso della proporzione di falsi positivi non tra tutti i test di ipotesi eseguiti ma solo rispetto al numero totale di ipotesi che sono state rifiutate.⁴

Essi hanno usato il termine tasso di false scoperte (FDR) per indicare tale quantità, e in termini formali seguendo la tabella 1 è definito come

dove E[.] indica il valore atteso. Il FDR è quindi il valore atteso del rapporto tra il numero di ipotesi erroneamente rifiutate e il numero totale di ipotesi rifiutate. L’operatore unione U al denominatore garantisce che FDR = 0 se nessuna ipotesi viene rifiutata. Il FDR fattorizza nel prodotto tra il valore atteso di false scoperte, condizionato ad almeno un test dichiarato significativo, e la probabilità di aver almeno un test in cui l’ipotesi nulla sia rifiutata.
Storey definisce il primo fattore nella scomposizione di cui sopra come il tasso di false scoperte “positivo” (pFDR). ⁵ Ci condizioniamo quindi all’occorrenza di almeno un rifiuto di ipotesi (R>0, almeno un risultato “positivo”)

L’interesse nel pFDR sta nel fatto che può essere visto come una probabilità a posteriori Bayesiana che ci permette di definire il q-valore, la controparte Bayesiana del p-valore.
I dettagli sul pFDR si possono trovare nella referenza numero 5. In questa sede il richiamo alla definizione formale del pFDR ci aiuta a definire il q-valore e a compararlo con l’usuale p-valore. Definito con H₀ l’ipotesi nulla, con H₁ l’ipotesi alternativa, con α il livello di probabilità di errore del I tipo, con T una generica statistica test e con Γ la regione di rifiuto, il pFDR può essere scritto come

che è la probabilità a posteriori di H₀. Si noti che il denominatore è la probabilità che la statistica test T cada nella regione di rifiuto

dove π₀ è la probabilità a priori che l’ipotesi nulla sia vera e π₁=1- π₀.

Se invece che considerare la regione di rifiuto Γ definita da un determinato livello di significatività α, si considera quella definita dal  p-valore, cioè dalla probabilità di osservare un valore della statistica test T uguale o più estremo a quello realmente osservato allora

Questa probabilità a posteriori è definita q-valore:
q-valore = Prob (H₀ | T>=T_obs)
e rappresenta l’analogo bayesiano del p-valore
p-valore = Prob (T>=T_obs | H₀)

Il p-valore da indicazione della forza dell’evidenza empirica contro l’ipotesi nulla (è la probabilità minima di errore di I tipo, che si realizza quando rifiutiamo l’ipotesi nulla sulla base del valore osservato della statistica test T). Il q-valore è il minimo pFDR nel quale possiamo incorrere se rifiutiamo l’ipotesi nulla sulla base del valore osservato o più estremo della statistica T; è una misura che tiene conto della molteplicità dei test e rappresenta quanto è probabile l’ipotesi nulla alla luce dei dati osservati.

FWER o FDR?

I due approcci alla correzione per test multipli rispondono a due filosofie diverse.Parafrasando Benjamini e Hochberg il controllo del FWER è importante quando la nostra decisione basata sull’intero insieme dei test potrebbe essere sbagliata quando almeno uno di essi sia stato erroneamente rifiutato. È il caso, per esempio, in cui vogliamo stabilire se la frequenza di tumori cutanei dipende dal colore dei capelli, degli occhi, della pelle o dal fototipo. Ciascuna caratteristica è misurata su scala qualitativa, nominale, le differenti modalità non sono ordinabili. La nostra decisione, che la frequenza di tumori cutanei varia al variare del colore dei capelli, per esempio, potrebbe essere sbagliata quando per almeno una delle modalità di colore dei capelli il test statistico sia stato erroneamente rifiutato. Nel caso della valutazione di variabili categoriche quindi è appropriato un controllo “globale” dell’errore di primo tipo.

Questa non è la situazione più comune. Per esempio, negli studi di genomica migliaia di geni vengono sottoposti allo stesso test per valutare se siano differenzialmente espressi. Ciascun gene però, implica che stiamo valutando un meccanismo diverso e quindi implicitamente una diversa ipotesi nulla. Così come negli studi clinici controllati quando siamo interessati a esiti diversi di un unico trattamento. Ma anche negli studi di epidemiologia ambientale, quando si valutano i rischi relativi per lunghe liste di malattie ma nella stessa popolazione (o viceversa, la stessa malattia in popolazioni diverse). Correggere per la molteplicità dei test seguendo l’approccio FDR è in tutte queste situazioni più appropriato in quanto “le conclusioni generali non sono erronee anche se qualcuna delle ipotesi nulle è falsamente rifiutata”.

La scelta di Barack Obama Leggo oggi su La Repubblica la traduzione dell’articolo di Bob Woodward su The Washington Post che ricostruisce la dinamica dell’azione che ha portato alla morte di Bin Laden. Siamo informati di un percorso in quindici tappe. Le prime sei sono relative alla pista che ha portato alla casa di Abbottabad. La prima evidenza empirica che ci si trovasse di fronte al covo di Bin Laden è fornita da immagini satellitari di un uomo molto alto che esce nel cortile per un’ora o due al giorno. Tale evidenza empirica non è ritenuta sufficiente. La probabilità è troppo bassa per giustificare un’azione. La CIA stabilisce una base ad Abbottabad e sorveglia con grande cautela il sospetto covo. Siamo alla fase finale. L’ultima stima della probabilità che si tratti del covo di Bin Laden è del 60-80%. La tensione è palpabile. Michael Letter, capo del centro nazionale americano anti-terrorismo, ancora in un recente incontro con Obama aveva dato una probabilità del 40%. È bassa? «È vero, ma è il 38% in più rispetto a prima!» Non è quindi il p-valore che viene considerato, cioè la forza dell’evidenza empirica contro l’ipotesi nulla (che Bin Laden NON sia ad Abbottabad), ma la probabilità a posteriori dell’ipotesi nulla (l’ultimo aggiornamento la dava al 20-40%) e il Fattore di Bayes, cioè quanto si è modificata la probabilità dell’ipotesi nulla (il 38% in più) passando dalla probabilità a priori (la penultima stima) a quella a posteriori (l’ultima valutazione). Come si vede le quantità in gioco sono la probabilità a priori, la probabilità a posteriori e il Fattore di Bayes. Nessuna traccia del p-valore e del livello di significatività…

 

Un semplice algoritmo di calcolo del q-valore

Supponiamo di avere m test di ipotesi per ognuno dei quali abbiamo calcolato un p-valore. Ordiniamo questi p-valori dal più piccolo al più grande. I corrispondenti q-valori possono essere ottenuti con un semplice algoritmo:

dove (i) denota l’i-esima ordinata quantità, m è il numero di test e i/m rappresenta Prob(T>=T_obs) valutata al cut off critico dato dall’i-esima statistica test.
Assumendo che la probabilità dell’ipotesi nulla sia pari a 1 (P(H₀)=1 nella formula) definire il q-valore è banale.

Due esempi: calcolo e interpretazione del q-valore

Riportiamo a titolo esplicativo due esempi. Il primo riguarda lo studio su 4 città italiane degli effetti delle ondate di calore.⁶ Il secondo esempio è tratto dal Rapporto sulle aree a rischio della Sardegna.⁷
Nel primo esempio richiamiamo la tabella 4 della referenza 6. Gli autori effettuano 28 test per valutare la presenza di una modificazione d’effetto attraverso l’indice di modificazione di effetto relativa (REM). In tabella 2 riportiamo la lista delle 28 malattie su cui il REM è stato calcolato, i relativi p-valori e la procedura per il calcolo “a mano” dei relativi q-valori.
L’ordine delle malattie della tabella 2 è diverso rispetto a quello riportato dagli autori originariamente⁶ e riflette l’ordinamento dei p-valori (dal più piccolo al più grande) necessario per il calcolo dei q-valori.

Per il calcolo del q-valore abbiamo modificato il p-valore associato a Tumori maligni da 0.000 a 0.001 (il risultato è leggermente più conservativo).
Basandosi sui p-valori 3 malattie presentano una modificazione di effetto: tumori maligni, psicosi e disturbi della conduzione. Fissata una soglia di criticità del 10% secondo il q-valore solo i tumori maligni vengono confermati come modificatori d’effetto. La correzione di Bonferroni da un livello di significatività di α*= 0.00357(0.10/28 il numero di test effettuati): anche in questo caso solo i tumori maligni vengono selezionati.
Per il secondo esempio abbiamo estratto 29 SMR aggiustati per deprivazione per il periodo 1997-2001 che si riferiscono a codici di malattia ICD-IX mutuamente esclusivi per i maschi dell’area di Portoscuso.
Lo scopo è identificare quali malattie si discostano rispetto al riferimento regionale.
In tabella 3 riportiamo il numero osservato e atteso di casi, i p-valori e i q-valori (calcolati con l’algoritmo di Storey) per le 10 malattie con p-valore =<0.10. La correzione di Bonferroni da un livello di significatività di α*= 0.00345 (0.10/29 il numero di cause analizzate): solo 3 p-valori sono sotto α*.

In base ai q-valori 5 malattie vengono identificate come “divergenti”. Di queste 5 malattie dichiarate diverse dall’ipotesi nulla “in media” 0.10×5 è un falso positivo, anche se non possiamo identificare quale dei 5 lo sia.

Conclusioni

In questo contributo abbiamo mostrato come trasformare il p-valore in una probabilità a posteriori, sotto certe assunzioni. Il p-valore è spesso interpretato in modo erroneo e valutare quanto grande è l’evidenza contro l’ipotesi nulla basandosi solo su esso non è semplice. La probabilità a posteriori dell’ipotesi nulla dipende dalla probabilità a priori che assegniamo all’ipotesi nulla stessa. Quando abbiamo un solo test d’ipotesi non resta che definire scenari di probabilità a priori. Si può obiettare che è non solo riduttivo ma anche molto difficile pensare di quantificare la probabilità a priori dell’ipotesi nulla. A noi non pare e abbiamo anche fornito alcuni esempi pratici e realistici. L’idea alla base è la stessa che si applica nel caso dei test diagnostici: il valore predittivo di un test (probabilità a posteriori) dipende dalla prevalenza della malattia (probabilità a priori). Sempre più in epidemiologia ci viene richiesto di supportare decisioni e questa nuova impostazione può aiutarci in questa parte del nostro lavoro. Quando invece abbiamo eseguito tanti test d’ipotesi, come spesso capita negli studi epidemiologici, la distribuzione dei p-valori calcolati contiene essa stessa l’informazione necessaria per calcolare la probabilità a posteriori o q-valore. Non solo. La disponibilità di una molteplicità di test permette anche di sacrificare qualcosa in termini di potenza e partire da una probabilità a priori dell’ipotesi nulla del 100%. Tutto quello che abbiamo scritto è parte di un atteggiamento pragmatico e un Bayesiano rigoroso si potrebbe anche risentire. Ma tant’è, l’importante è cominciare ad abbandonare l’uso acritico del p-valore.

Nota sulle procedure di calcolo del q-valore: L’algoritmo proposto da Storey nel 2002 per il calcolo del q-valore e implementato in R nella libreria “qvalue” produce una stima di P(H0), ovvero della proporzione di ipotesi nulle vere che si basa su un parametro di “tuning”. Se questo parametro viene imposto essere uguale a zero allora il q-valore ottenuto sarà identico a quello dell’uso della formula riportata nel testo quando assumiamo P(H0=1). Nel software STATA 11 il comando “qqvalue” (con l’opzione Simes) riproduce i risultati riportati nella tabella 2. Nei due esempi abbiamo riportato il q-valore calcolato “a mano” e usando la procedura di Storey lasciando che P(H0) venga stimata.

Bibliografia

1. Biggeri A, Catelan D, Barbone F. Reporting uncertainty. Epidemiol Prev 2010;34(5-6):91-5.
2. Stern JAC, Smith DG. Sifting the Evidence. What’s Wrong with Significance Tests? BMJ 2001;322,226-231.
3. Bonferroni C. Teoria Statistica delle Classi e Calcolo delle Probabilità. Pubblicazioni del Regio istituto superiore di scienze economiche e commerciali di Firenze, 1936, 3-62.
4. Benjamini Y, Hochberg, Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society 1995;Series B 57,289-300.
5. Storey JD. The Positive False Discovery Rate: a Bayesian Interpretation and the Q-Value. Annals of Statistics 2003;31,2013-2035.
6. Stafoggia M, Forastiere F, Agostini D, et al. Vulnerability to Heat-Related Mortality A Multicity, Population-Based, Case-Crossover Analysis. Epidemiology 2006; 17:315-323
7. Biggeri A, Lagazio C, Catelan D, Pirastu R, Casson F, Terracini B. Environment and Health in Sardinia. Epidemiol Prev 1996; 30S,1-96.