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

Resumen Técnico: División Multinivel Adaptativa Aumentada por Hash para Valores p de Pruebas de Permutación

Planteamiento del Problema
Las pruebas de permutación no paramétricas, como las de Kolmogorov–Smirnov y Mann–Whitney U, son fundamentales para las pruebas de hipótesis estadística, particularmente al comparar dos grupos sin asumir una distribución de datos específica. Si bien el cálculo exacto del valor p es ideal, suele ser algorítmicamente intratable para estadísticas de prueba complejas o personalizadas debido al número exponencial de permutaciones posibles ($\binom{n+m}{n}$). En consecuencia, los investigadores dependen del muestreo de Monte Carlo. Sin embargo, los métodos estándar de Monte Carlo sufren de una pobre precisión relativa al estimar valores p muy pequeños (por ejemplo, $10^{-10}$ o menores), lo que requiere un número impráctico de muestras para lograr potencia estadística tras múltiples correcciones de hipótesis.

Además, la aplicación de la división multinivel (multilevel splitting) —una técnica diseñada para estimar probabilidades de eventos raros mediante la descomposición de una probabilidad en un producto de probabilidades condicionales— a las pruebas de permutación introduce desafíos específicos. Las pruebas de permutación operan en espacios discretos donde la distribución de la estadística de prueba suele presentar saltos grandes (discretidad). En tales casos, muchas combinaciones distintas de datos pueden generar valores de la estadística idénticos. Esta discretidad puede impedir la selección adaptativa de niveles intermedios requerida para la división multinivel, causando que el algoritmo falle o produzca estimaciones sesgadas cuando la probabilidad condicional de avanzar al siguiente nivel no puede mantenerse cerca de una constante objetivo (por ejemplo, 0.5).

Metodología
Los autores proponen un algoritmo de Monte Carlo de división multinivel adaptativa aumentada por hash para abordar los problemas de discretidad inherentes a las pruebas de permutación. La metodología central involucra tres componentes clave:

1. Ordenamiento Aumentado por Hash: Para resolver el problema de múltiples combinaciones que comparten el mismo valor estadístico, los autores introducen un orden total estricto sobre las combinaciones. En lugar de ordenar únicamente por el valor de la estadística $S(A)$, definen un orden $A \succ_{S,h} B$ basado primero en el valor de la estadística y luego por un valor de hash $h(A)$ si las estadísticas son iguales.

   • La función de hash $h(A)$ mapea las combinaciones a enteros. Una implementación práctica utiliza una suma de enteros independientes asignados a cada elemento, lo que permite actualizaciones de $O(1)$ durante los pasos de Monte Carlo por Cadenas de Markov (MCMC).
   • Este aumento efectivamente "suaviza" la distribución discreta, asegurando que para cualquier límite de nivel, se pueda encontrar un siguiente límite válido de modo que la probabilidad condicional sea estimable.

2. División Multinivel Adaptativa con MCMC: El algoritmo estima el valor p $Pr(S \ge \gamma)$ descomponiéndolo en un producto de probabilidades condicionales sobre una secuencia de eventos anidados definidos por el orden aumentado por hash.

   • Niveles Adaptativos: Partiendo de un nivel trivial, el algoritmo muestrea combinaciones y utiliza sus cuantiles (por ejemplo, la mediana) para definir el siguiente límite de nivel.
   • Muestreo MCMC: Para generar muestras para el siguiente nivel, el algoritmo emplea un procedimiento MCMC de estilo Metropolis; comienza con muestras existentes que satisfacen la nueva condición de límite y las perturba (intercambiando elementos) para generar nuevas muestras independientes.
   • Estrategias de Remuestreo: El artículo evalúa el "remuestreo completo" (duplicar las muestras válidas existentes y perturbar todas ellas) frente al "remuestreo parcial" (perturbar solo un subconjunto). Se encontró que el remuestreo completo es más robusto contra los sesgos de muestreo introducidos por la naturaleza discreta del problema.

3. Estimación en el Espacio Logarítmico: El algoritmo estima el logaritmo del valor p para manejar valores extremos sin subdesbordamiento numérico (underflow). Utiliza las funciones digamma ($\psi$) y trigamma ($\psi_1$) para calcular la media y la varianza de las estimaciones logarítmicas, derivadas de las propiedades de las estadísticas de orden de distribuciones uniformes. Esto permite la construcción de intervalos de confianza válidos.

Contribuciones Clave

• Innovación Algorítmica: La introducción de un esquema de ordenamiento aumentado por hash que permite que la división multinivel adaptativa funcione correctamente en espacios de permutación discretos, superando el problema de los "saltos grandes" donde los métodos estándar fallan.
• Implementación Robusta: Un análisis detallado de las estrategias de remuestreo (completo vs. parcial) y los criterios de parada de MCMC, demostrando que el remuestreo completo con un parámetro de escala específico ($\alpha=1$) proporciona el mejor equilibrio entre velocidad de convergencia y reducción de sesgo.
• Lanzamiento de Software: El desarrollo y lanzamiento de código abierto del paquete de Python hamstest, que implementa el algoritmo propuesto. Permite a los usuarios definir estadísticas de prueba personalizadas y proporciona interfaces eficientes (incluyendo backends en C++) para la estimación de valores p de alto rendimiento.
• Validación: Comparación rigurosa contra algoritmos exactos para las pruebas de Kolmogorov–Smirnov y Mann–Whitney U, demostrando que el método propuesto estima con precisión valores p tan pequeños como $10^{-243}$ con intervalos de confianza válidos.

Resultados
Los autores evaluaron el algoritmo en pruebas tanto de Kolmogorov–KS como de Mann–Whitney U a través de varios tamaños de muestra ($n+m$ hasta 10,000) y valores p objetivo (hasta $10^{-100}$).

• Exactitud: Los valores p estimados coincidieron estrechamente con los valores reales (calculados mediante métodos exactos cuando estaban disponibles), sin sesgo visible. La mediana de las estimaciones se alineó con el valor p real en todos los parámetros probados.
• Intervalos de Confianza: Los intervalos de confianza del 95% construidos (basados en $\log \hat{p} \pm 2\hat{\sigma}$) capturaron con éxito el logaritmo del valor p real en la gran mayoría de las ejecuciones independientes, validando las estimaciones teóricas de la varianza.
• Rendimiento: Se demostró que los tiempos de ejecución dependen principalmente de la magnitud del valor p y del tamaño de la muestra. El algoritmo demostró una escalabilidad lineal respecto al número de niveles requeridos, lo que lo hace factible para estimar valores p extremadamente pequeños que son computacionalmente prohibitivos para el Monte Carlo estándar.
• Comparación de Remuestreo: La variante "full-median" (remuestreo completo con selección de límite por mediana) superó consistentemente a las variantes de remuestreo parcial en términos de velocidad de convergencia y robustez al sesgo, particularmente para la prueba de Kolmogorov–$\text{KS}$, donde la discretidad de la distribución es más pronunciada.

Significancia
El artículo afirma que este trabajo representa una mejora mayor sobre los enfoques clásicos de Monte Carlo para la estimación de valores p en pruebas de permutación. Al abordar los obstáculos específicos de los espacios discretos —a saber, la incapacidad de definir niveles intermedios debido a los valores de la estadística empatados— el algoritmo propuesto permite la estimación precisa y eficiente de valores p arbitrariamente pequeños. Esta capacidad es crítica para aplicaciones modernas de alto rendimiento (por ejemplo, genómica, neurociencia) donde se realizan miles de pruebas simultáneamente, y el control preciso de las tasas de falso descubrimiento requiere una estimación fiable de las probabilidades de cola. Los autores señalan que una versión anterior de este algoritmo (sin el aumento por hash) ya había tenido éxito en el Análisis de Enriquecimiento de Conjuntos de Genes (GSEA); el presente trabajo extiende esta utilidad a una clase más amplia de estadísticas de dos muestras y resuelve los problemas de inestabilidad encontrados con distribuciones de entrada específicas. La disponibilidad del paquete hamstest asegura que estas capacidades estadísticas avanzadas sean accesibles para estadísticas definidas por el usuario en flujos de trabajo de investigación prácticos.