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

技術要約：置換検定のp値のためのハッシュ増強型適応的マルチレベル分割法

問題提起
コルモゴロフ–スミルノフ検定やマン・ホイットニーのU検定などの非パラメトリックな置換検定は、特定のデータ分布を仮定せずに2つのグループを比較するための、統計的仮説検定の基本となる手法である。厳密なp値の計算が理想的であるが、複雑またはカスタムの検定統計量の場合、可能な置換の数（$\binom{n+m}{n}$）が指数関数的に増大するため、アルゴリズム的に困難（intractable）になることが多い。その結果、研究者はモンテカルロ・サンプリングに頼ることになる。しかし、標準的なモンテカルロ法は、非常に小さなp値（例：$10^{-10}$以下）を推定する場合、相対的な精度が低くなるという問題を抱えており、複数の仮説補正を行った後の統計的検出力を確保するために、非現実的な数のサンプルを必要とする。

さらに、稀な事象の確率を条件付き確率の積へと分解することで事象の確率を推定する手法であるマルチレベル分割（multilevel splitting）を置換検定に適用すると、特有の課題が生じる。置換検定は離散的な空間で動作し、検定統計量の分布はしばしば大きな跳躍（離散性）を示す。このような場合、多くの異なる組み合わせが同一の統計量値を持つことがある。この離散性は、マルチレベル分割に必要な中間レベルの適応的な選択を妨げ、条件付き確率を目標とする定数（例：0.5）付近に維持できない場合に、アルゴリズムの失敗や偏った推定を引き起こす原因となる。

手法
著者らは、置換検定に固有の離散性の問題を解決するために、ハッシュ増強型適応的マルチレベル分割モンテカルロ・アルゴリズムを提案している。コアとなる手法は、以下の3つのコンポーネントで構成される。

1. ハッシュ増強型順序付け: 統計量が同じ値を持つ組み合わせの問題を解決するために、著者らは組み合わせに対する厳密な全順序を導入している。単に統計量 $S(A)$ だけで順序付けるのではなく、統計量値を優先し、統計量が等しい場合にはハッシュ値 $h(A)$ に基づいて $A \succ_{S,h} B$ という順序を定義する。

   • ハッシュ関数 $h(A)$ は、組み合わせを整数にマッピングする。実用的な実装では、各要素に割り当てられた独立したランダム整数の和を使用しており、これによりマルコフ連鎖モンテカルロ（MCMC）ステップ中の更新を $O(1)$ で行うことができる。
   • この増強は、離散分布を事実上「平滑化」し、任意のレベル境界に対して、条件付き確率が推定可能な状態を維持できる有効な次の境界を見つけることを保証する。

2. MCMCを用いた適応的マルチレベル分割: アルゴリズムは、ハッシュ増強型順序に基づく一連の入れ子状のイベントによって定義される、p値 $Pr(S \ge \gamma)$ を条件付き確率の積へと分解することで推定する。

   • 適応的レベル: アルゴリズムは、組み合わせをサンプリングし、それらの分位点（例：中央値）を使用して次のレベル境界を定義する。
   • MCMCサンプリング: 次のレベルのためのサンプルを生成するために、メトロポリス型のMCMC法を採用する。これは、新しい境界条件を満たす既存のサンプルから開始し、それらを摂動（要素の入れ替え）させることで、新しい独立したサンプルを生成する。
   • 再サンプリング戦略: 本論文では、「フル・リサンプリング」（既存の有効なサンプルを複製し、それらすべてに摂動を加える）と「部分的リサンプリング」（サブセットのみに摂動を加える）を評価している。問題の離散性によって導入されるサンプリングバイアスに対しては、フル・リサンプリングの方が堅牢であることが示されている。

3. 対数空間推定: 数値的なアンダーフローを回避して極端な値を扱うために、p値の対数を推定する。アルゴリズムは、一様分布の順序統計量の性質から導出されるディガマ（$\psi$）関数およびトリガマ（$\psi_1$）関数を利用して、対数推定値の平均と分散を計算する。これにより、有効な信頼区間の構築が可能となる。

主な貢献

• アルゴリズムの革新: 離散的な置換空間において適応的マルチレベル分割を正しく機能させるための、ハッシュ増強型順序付けスキームの導入。これは、標準的な手法が失敗する「大きな跳躍」の問題を克服する。
• 堅牢な実装: 再サンプリング戦略（フル vs 部分）およびMCMCの停止基準に関する詳細な分析を行い、スケーリングパラメータ（$\alpha=1$）を用いたフル・リサンプリングが、収束速度とバイアス低減の最適なバランスを提供することを実証した。
• ソフトウェアのリリース: 提案されたアルゴリズムを実装したPythonパッケージ hamstest の開発およびオープンソースとしての公開。これにより、ユーザーはカスタムの検定統計量を定義でき、高パフォーマンスなp値推定のための効率的なインターフェース（C++バックエンドを含む）を利用できる。
• 検証: コルモゴロフ–スミルノフ検定およびマン・ホイットニーのU検定における厳密なアルゴリズムとの詳細な比較を行い、提案手法が有効な信頼区間を伴いつつ、$10^{-243}$ という極めて小さなp値を正確に推定できることを示した。

結果
著者らは、様々なサンプルサイズ（$n+m$ 最大10,000）およびターゲットp値（$10^{-100}$ まで）において、コルモゴロフ–スミルノフ検定とマン・ホイットニーのU検定の両方でアルゴリズムを評価した。

• 精度: 推定されたp値は、真の値（利用可能な場合は厳密な手法により算出）と密接に一致しており、目立ったバイアスは見られなかった。すべてのテストされたパラメータにおいて、推定値の中央値は真のp値と一致した。
• 信頼区間: 構築された95%信頼区間（$\log \hat{p} \pm 2\hat{\sigma}$）は、大多数の独立した実行において真の対数p値を捉えることに成功し、理論的な分散推定の妥当性を裏付けた。
• パフォーマンス: 実行時間は主にp値の大きさとサンプルサイズに依存することが示された。アルゴリズムは、必要なレベル数に対して線形なスケーラビリティを示し、標準的なモンテカルロ法では計算的に不可能なほど極端に小さなp値の推定を可能にした。
• 再サンプリングの比較: 「フル・メディアン」バリアント（中央値境界選択を用いたフル・リサンプリング）は、特に分布の離散性がより顕著なコルモゴロフ–スミルノフ検定において、収束速度とバイアスに対する堅牢性の両面で、部分的リサンプリングのバリアントを一貫して上回った。

意義
本論文は、本研究が置換検定におけるp値推定のための古典的なモンテカルロ・アプローチに対する大幅な改善を提供すると主張している。離散的な空間に特有の落とし穴（すなわち、タイ（同値）が発生する統計量値のために中間レベルを定義できないこと）に対処することで、提案されたアルゴリズムは、任意に小さなp値の正確かつ効率的な推定を可能にする。この能力は、数千のテストが同時に行われ、精緻な偽発見率の制御が極めて微細な裾の確率の信頼できる推定を必要とする、現代のハイスループットなアプリケーション（ゲノミクス、神経科学など）において極めて重要である。著者らは、ハッシュ増強なしの以前のバージョンのアルゴリズムも、遺伝子セット濃縮解析（GSEA）において既に成功を収めていたことを指摘しており、本研究はその有用性をより広範な2標本統計へと拡張し、特定の入力分布で遭遇する不安定性の問題を解決したものである。hamstest パッケージの提供により、これらの高度な統計機能が、実用的な研究ワークフローにおけるユーザー定義の統計量に対して利用可能となる。