Hash-augmented adaptive multilevel splitting Monte Carlo algorithm for accurate estimation of two-sample permutation test p-values

Sintesi Tecnica: Hash-Augmented Adaptive Multilevel Splitting per i p-Value dei Test di Permutazione

Definizione del Problema
I test di permutazione non parametrici, come i test di Kolmogorov–Smirnov e Mann–Whitney U, sono fondamentali per l'inferenza statistica, in particolare per confrontare due gruppi senza assumere una specifica distribuzione dei dati. Sebbene il calcolo esatto del p-value sia l'ideale, esso è spesso algoritmicamente intrattabile per statistiche di test complesse o personalizzate a causa del numero esponenziale di possibili permutazioni ($\binom{n+m}{n}$). Di conseguenza, i ricercatori si affidano al campionamento Monte Carlo. Tuttavia, i metodi Monte Carlo standard soffrono di una scarsa accuratezza relativa quando stimano p-value molto piccoli (ad esempio, $10^{-10}$ o inferiori), richiedendo un numero impraticabile di campioni per ottenere la potenza statistica dopo molteplici correzioni per test multipli.

Inoltre, l'applicazione del multilevel splitting — una tecnica progettata per stimare le probabilità di eventi rari scomponendole in un prodotto di probabilità condizionate — ai test di permutazione introduce sfide specifiche. I test di permutazione operano in spazi discreti in cui la distribuzione della statistica di test presenta spesso grandi salti (discretizzazione). In tali casi, molte combinazioni distinte di dati possono produrre valori della statistica identici. Questa discretizzazione può impedire la selezione adattiva di livelli intermedi necessari per il multilevel splitting, causando il fallimento dell'algoritmo o la produzione di stime distorte quando la probabilità condizionata di passare al livello successivo non può essere mantenuta vicino a una costante target (ad esempio, 0,5).

Metodologia
Gli autori propongono un algoritmo Monte Carlo di hash-augmented adaptive multilevel splitting per affrontare i problemi di discretizzazione inerenti ai test di permutazione. La metodologia principale coinvolge tre componenti chiave:

1. Ordinamento con Augmentazione tramite Hash (Hash-Augmented Ordering): Per risolvere il problema delle molteplici combinazioni che condividono lo stesso valore della statistica, gli autori introducono un ordine totale rigoroso sulle combinazioni. Invece di ordinare esclusivamente in base alla statistica $S(A)$, definiscono un ordine $A \succ_{S,h} B$ basato prima sul valore della statistica e, in caso di parità, sul valore di un hash $h(A)$.

   • La funzione di hash $h(A)$ mappa le combinazioni in interi. Un'implementazione pratica utilizza la somma di interi casuali indipendenti assegnati a ciascun elemento, consentendo aggiornamenti in $O(1)$ durante le fasi di Markov Chain Monte Carlo (MCMC).
   • Questa augmentazione "ammorbidisce" efficacemente la distribuzione discreta, garantendo che per ogni confine di livello sia possibile trovare un valido confine successivo in modo che la probabilità condizionata rimanga stimabile.

2. Adaptive Multilevel Splitting con MCMC: L'algoritmo stima il p-value $Pr(S \ge \gamma)$ scomponendolo in un prodotto di probabilità condizionate su una sequenza di eventi annidati definiti dall'ordine aumentato tramite hash.

   • Livelli Adattivi (Adaptive Levels): Partendo da un livello banale, l'algoritmo campiona le combinazioni e utilizza i loro quantili (ad esempio, la mediana) per definire il successivo confine di livello.
   • Campionamento MCMC: Per generare campioni per il livello successivo, l'algoritmo impiega una procedura MCMC di tipo Metropolis. Parte da combinazioni esistenti che soddisfano la nuova condizione di confine e le perturba (scambiando elementi) per generare nuovi campioni indipendenti.
   • Strategie di Resampling: Il documento valuta il "full resampling" (duplicazione dei campioni validi esistenti e perturbazione di tutti essi) rispetto al "partial resampling" (perturbazione solo di un sottoinsieme). Il "full resampling" risulta essere più robusto contro i bias di campionamento introdotti dalla natura discreta del problema.

3. Stima nello Spazio Logaritmico (Log-Space Estimation): L'algoritmo stima il logaritmo del p-value per gestire valori estremi senza problemi di underflow numerico. Utilizza le funzioni digamma ($\psi$) e trigamma ($\psi_1$) per calcolare la media e la varianza delle stime logaritmiche, derivate dalle proprietà delle statistiche d'ordine delle distribuzioni uniformi. Ciò consente la costruzione di intervalli di confidenza validi.

Contributi Chiave

• Innovazione Algoritmica: L'introduzione di uno schema di ordinamento aumentato tramite hash che permette al multilevel splitting adattivo di funzionare correttamente negli spazi di permutazione discreti, superando il problema dei "grandi salti" dove i metodi standard falliscono.
• Implementazione Robusta: Un'analisi dettagliata delle strategie di resampling (full vs. partial) e dei criteri di arresto MCMC, dimostrando che il full resampling con un parametro di scala specifico ($\alpha=1$) fornisce il miglior equilibrio tra velocità di convergenza e riduzione del bias.
• Rilascio del Software: Lo sviluppo e il rilascio open-source del pacchetto Python hamstest, che implementa l'algoritmo proposto. Esso permette agli utenti di definire statistiche di test personalizzate e fornisce interfacce efficienti (inclusi backend in C++) per la stima ad alte prestazioni dei p-value.
• Validazione: Un confronto rigoroso con algoritmi esatti per i test di Kolmogorov–KS e Mann–Whitney U, dimostrando che il metodo proposto stima accuratamente p-value anche di entità $10^{-243}$ con intervalli di confidenza validi.

Risultati
Gli autori hanno valutato l'algoritmo sia sui test di Kolmogorov–KS che su quelli di Mann–Whitney U attraverso varie dimensioni campionarie ($n+m$ fino a 10.000) e p-value target (fino a $10^{-100}$).

• Accuratezza: I p-value stimati corrispondono strettamente ai valori reali (calcolati tramite metodi esatti ove disponibili), senza alcun bias visibile. La mediana delle stime si allinea al vero p-value in tutti i parametri testati.
• Intervalli di Confidenza: Gli intervoli di confidenza al 95% costruiti (basati su $\log \hat{p} \pm 2\hat{\sigma}$) hanno catturato con successo il vero log p-value nella stragrande maggioranza delle esecuzioni indipendenti, validando le stime teoriche della varianza.
• Performance: I tempi di esecuzione sono risultati dipendere principalmente dall'entità del p-value e dalla dimensione del campione. L'algoritmo ha mostrato una scalabilità lineare rispetto al numero di livelli richiesti, rendendo fecondo stimare p-value estremamente piccoli che sarebbero computazionalmente proibitivi per il Monte Carlo standard.
• Confronto Resampling: La variante "full-median" (full resampling con selezione del confine mediano) ha costantemente superato le varianti di partial resampling in termini di velocità di convergenza e robustezza al bias, particolarmente nel test di Kolmogorov–KS dove la discretizzazione della distribuzione è più pronunciata.

Significatività
Il documento sostiene che questo lavoro rappresenti un miglioramento fondamentale rispetto agli approcci Monte Carlo classici per la stima dei p-value nei test di permutazione. Affrontando le criticità specifiche degli spazi discreti — ovvero l'impossibilità di definire livelli intermedi a causa dei valori di statistica uguali (ties) — l'algoritzione proposto consente la stima accurata ed efficiente di p-value arbitrariamente piccoli. Questa capacità è cruciale per le moderne applicazioni ad alto throughput (ad esempio, genomica, neuroscienze) dove vengono eseguiti simultaneamente miglia di test e dove il controllo preciso dei tassi di falsa scoperta richiede una stima affidabile delle probabilità di coda. Gli autori osservano che una versione precedente di questo algoritmo (senza l'aumento tramite hash) era già stata utilizzata con successo nella Gene Set Enrichment Analysis (GSEA); il presente lavoro estende questa utilità a una classe più ampia di statistiche a due campioni e risolve i problemi di instabilità riscontrati con specifiche distribuzioni di input. La disponibilità del pacchetto hamstest assicura che queste avanzate capacità statistiche siano accessibili per statistiche definite dall'utente nei flussi di lavoro della ricerca pratica.