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

기술 요약: 순열 검정 p-값에 대한 해시 증강 적응형 다단계 분할법

문제 정의
Kolmogorov–KS 및 Mann–Whitney U 검정과 같은 비모수적 순열 검정은 특정 데이터 분포를 가정하지 않고 두 집단을 비교하는 데 있어 기초적인 통계적 가설 검정 방법이다. 엄밀한 p-값 계산이 이상적이지만, 복잡하거나 사용자 정의된 검정 통계량의 경우 가능한 순열의 수($\binom{n+m}{n}$)가 기하급수적으로 많아져 알고리즘적으로 다루기 어려운 경우가 많다. 결과적으로 연구자들은 몬테카를로 샘플링에 의존하게 된다. 그러나 표준 몬테카를로 방법은 매우 작은 p-값(예: $10^{-10}$ 이하)을 추정할 때 상대적 정확도가 떨어지며, 여러 가설 교정을 거친 후 통계적 검정력을 확보하기 위해 필요한 샘플 수가 비현실적으로 많아지는 문제가 발생한다.

또한, 희귀 사건 확률을 추정하기 위해 설계된 기법인 다단계 분할법(multilevel splitting)을 순열 검정에 적용할 때는 특정한 어려움이 따른다. 순열 검정은 이산적인 공간에서 작동하며, 이 과정에서 검정 통계량의 분포는 큰 도약(discreteness)을 보일 수 있다. 이러한 경우, 동일한 통계치 값을 갖는 서로 다른 조합들이 다수 존재할 수 있다. 이러한 이산성은 다단계 분할법에 필요한 중간 단계(intermediate levels)의 적응적 선택을 방해하여, 목표 상수(예: 0.5) 근처에서 조건부 확률을 유지할 수 없게 만듦으로써 알고리즘을 실패하게 하거나 편향된 추정치를 생성하게 할 수 있다.

방법론
저자들은 순열 검정의 이산성 문제를 해결하기 위해 해시 증강 적응형 다단계 분질법 몬테카를로 알고리즘을 제안한다. 핵심 방법론은 세 가지 주요 구성 요소로 이루어진다:

1. 해시 증강 정렬 (Hash-Augmented Ordering): 여러 조합이 동일한 통계치 값을 공유하는 문제를 해결하기 위해, 저자들은 조합에 대한 엄격한 전순서(total order)를 도입한다. 단순히 통계량 $S(A)$에 의해서만 정렬하는 대신, 통계치 값을 우선하고 통계치가 같을 경우 해시 값 $h(A)$에 따라 $A \succ_{S,h} B$로 정의한다.

   • 해시 함수 $h(A)$는 조합을 정수로 매핑한다. 실질적인 구현에서는 각 요소에 할당된 독립적인 무작위 정수의 합을 사용하여, MCMC(Markov Chain Monte Carlo) 단계 동안 $O(1)$의 업데이트가 가능하도록 한다.
   • 이러한 증강은 이산 분포를 효과적으로 "매끄럽게(smoothing)" 만들어, 어떤 레벨 경계에서도 조건부 확률을 추정 가능한 상태로 유지할 수 있는 유효한 다음 경계를 찾을 수 있도록 보장한다.

2. MCMC를 이용한 적응형 다단계 분할: 알고리즘은 해시 증강 정렬에 의해 정의된 일련의 중첩된 사건들을 통해 $Pr(S \ge \gamma)$를 분해하여 p-값을 추정한다.

   • 적응형 레벨 (Adaptive Levels): 알고리즘은 자명한 레벨에서 시작하여 조합을 샘플링하고, 그 분위수(예: 중앙값)를 사용하여 다음 레벨 경계를 정의한다.
   • MCMC 샘플링: 다음 레벨을 위한 샘플을 생성하기 위해, 알고리즘은 메트로폴리스 스타일의 MCMC 절차를 채택한다. 이는 새로운 경계 조건을 만족하는 기존 샘플에서 시작하여, 요소를 교환(swapping)함으로써 새로운 독립 샘플을 생성한다.
   • 재샘플링 전략: 논문은 "전체 재샘플링(full resampling, 기존의 유효한 샘플을 복제하고 모두 변형함)"과 "부분 재샘플링(partial resampling, 일부만 변형함)"을 평가한다. 전체 재샘플링이 문제의 이산성으로 인해 발생하는 샘플링 편향에 대해 더 견고한 것으로 나타났다.

3. 로그 공간 추정 (Log-Space Estimation): 알고리즘은 수치적 언더플로우 없이 극단적인 값을 처리하기 위해 p-값의 로그를 추정한다. 또한, 균등 분포의 순서 통계량 특성으로부터 유도된 디감마($\psi$) 및 트리감마($\psi_1$) 함수를 활용하여 로그 추정치의 평균과 분산을 계산한다. 이를 통해 유효한 신뢰 구간을 구축할 수 있다.

주요 기여

• 알고리즘 혁신: 순열의 이산 공간에서 적응형 다단계 분할법이 정상적으로 작동할 수 있도록 해시 증강 정렬 체계를 도입함으로써, 표준 방식이 실패하는 "큰 도약(large jump)" 문제를 극복하였다.
• 견고한 구현: 재샘플링 전략(전체 vs 부분) 및 MCMC 정지 기준에 대한 상세한 분석을 통해, 특정 스케일링 파라미터($\alpha=1$)를 사용하는 전체 재샘플링이 수렴 속도와 편향 감소 사이의 최적의 균형을 제공함을 입증하였다.
• 소프트웨어 공개: 제안된 알고리즘을 구현한 Python 패키지 **hamstest**를 개발하여 오픈 소스로 공개하였다. 이를 통해 사용자는 사용자 정의 검정 통계량을 정의할 수 있으며, 고성능 p-값 추정을 위한 효율적인 인터페이스(C++ 백엔드 포함)를 제공받을 수 있다.
• 검증: Kolmogorov–KS 및 Mann–Whitney U 검정에 대해 엄밀한 알고리즘과 비교하여, 제안된 방법이 유효한 신뢰 구간과 함께 $10^{-243}$만큼 작은 p-값까지 정확하게 추정함을 입증하였다.

결과
저자들은 다양한 표본 크기($n+m$ 최대 10,000)와 목표 p-값(최저 $10^{-100}$)에 대해 Kolmogorov–KS 및 Mann–Whitney U 검정 모두에서 알고리즘을 평가하였다.

• 정확도: 추정된 p-값은 (가능한 경우 엄밀한 방법으로 계산된) 실제 값과 밀접하게 일치하였으며, 눈에 띄는 편향은 발견되지 않았다. 모든 테스트 파라미터에 걸쳐 추정치의 중앙값이 실제 p-값과 일치하였다.
• 신뢰 구간: $\log \hat{p} \pm 2\hat{\sigma}$를 기반으로 구축된 95% 신뢰 구간은 대부분의 독립적인 실행에서 실제 로그 p-값을 포착하였으며, 이는 이론적 분산 추정치의 타당성을 입증한다.
• 성능: 실행 시간은 주로 p-값의 크기와 표본 크기에 의존하는 것으로 나타났다. 알고리즘은 필요한 레벨 수에 대해 선형적 확장성을 보여주었으며, 이는 표준 몬테카를로로는 계산이 불가능한 극도로 작은 p-값을 추정하는 것을 가능하게 한다.
• 재샘플링 비교: "전체-중앙값(full-median)" 변형(중앙값 경계 선택을 사용하는 전체 재샘플링)은 특히 분포의 이산성이 더 뚜렷한 Kolmogorov–KS 검정에서 부분 재샘플링 변형보다 수렴 속도와 편향에 대한 견고성 측면에서 일관되게 우수한 성능을 보였다.

의의
본 논문은 이 작업이 순열 검정의 p-값 추정에 있어 고전적인 몬테카를로 접근법보다 크게 개선된 방식임을 주장한다. 제안된 알고리즘은 이산 공간의 구체적인 난점, 즉 결합된 통계치 값으로 인해 중간 단계를 정의할 수 없는 문제를 해결함으로써, 임의로 작은 p-값을 정확하고 효율적으로 추정할 수 있게 한다. 이러한 능력은 수천 개의 테스트가 동시에 수행되는 현대의 고처리량 응용 분야(예: 유전체학, 신경과학)에서 매우 중요하며, 여기서는 정밀한 허위 발견율(false discovery rate) 제어를 위해 신뢰할 수 있는 꼬리 확률 추정이 필수적이다. 저자들은 해시 증강이 없는 이전 버전의 알고리즘이 이미 유전자 세트 풍부성 분석(GSEA)에서 성공적이었음을 언급하며, 본 연구는 이를 더 넓은 범위의 이표본 통계량으로 확장하고 특정 입력 분포에서 발생하는 불안정성 문제를 해결했다고 밝히고 있다. hamstest 패키지의 가용성은 이러한 고급 통계적 역량이 실제 연구 워크플로에서 사용자 정의 통계량을 위해 접근 가능하도록 보장한다.