Accelerating Bayesian Phylogenetic Inference via Delayed Acceptance Sequential Monte Carlo with Random Forest Surrogates

技術的サマリー：ランダムフォレスト代理モデルを用いた遅延受容 Sequential Monte Carlo によるベイズ系統推論の高速化

問題定義
ベイズ系統学における主要な目的は、系統樹（トポロジーと分岐長）および進化パラメータに関する事後分布を推定することである。これには、扱いにくい数の樹トポロジーにわたる総和と、高次元の分岐長にわたる積分を必要とする事後分布の計算が含まれる。マルコフ連鎖モンテカルロ（MCMC）法や Sequential Monte Carlo（SMC）法はこれらの分布を近似するための標準的な手法であるが、大規模な配列データにおいては重大なボトルネックに直面する：尤度関数の反復評価が計算上極めて高コストであることである。適応的 MCMC 移動や変分近似などの既存のアプローチはこれを緩和しようとしたが、正確な尤度計算の計算コストは推論時間を支配する主要因のままである。

手法
著者らは、樹空間の探索を加速するために、**遅延受容（Delayed Acceptance: DA）メカニズムをSequential Monte Carlo（SMC）**に統合し、**ランダムフォレスト（Random Forest: RF）**代理モデルを利用する新たな枠組みを提案する。

1. 遅延受容 MCMC カーネル:
   中核的な革新は、却下されることが予想される提案に対して真の尤度を計算することを回避するように設計された、2段階（実質的に3ステップ）の受容手順である。

   • ステップ 1（代理モデルによるチェック）: 新しい状態 $x^*$ が提案されると、アルゴリズムはまず機械学習モデルによって予測された対数尤度の変化 $\Delta x^*$ を評価する。もし $\Delta x^* < \kappa$ （ここで $\kappa$ はユーザーが指定する負の閾値）であれば、提案は真の尤度、事前分布、または提案比率を計算することなく即座に却下される。
   • ステップ 2（疑似受容）: 提案が最初のチェックを通過した場合、アルゴリズムは事前分布と提案の比率を代理尤度 $f(x^*)$ と組み合わせて計算する。2 番目の受容確率 $\tilde{\alpha}$ が計算される。
   • ステップ 3（真の尤度）: 提案が最初の 2 つの段階を通過した場合にのみ、最終的な受容確率 $\bar{\alpha}$ を決定するために正確な尤度 $p(D|x^*)$ が計算される。
   • 調整パラメータ（$\delta$）: 計算速度と事後分布の精度の間のトレードオフを制御するために、代理対数尤度の変化に調整パラメータ $\delta$ が追加される。$\delta$ が大きいほど、誤った却下（第 2 種の誤り）は減少するが、真の尤度評価の回数は増加する。

2. ランダムフォレスト代理モデル:
   系統樹トポロジーは離散的な状態空間に存在するため、標準的な機械学習は直接尤度をモデル化できない。代わりに、著者らは樹の移動によって引き起こされる対数尤度の変化をモデル化する。

   • 特徴量抽出: この手法は、樹の提案（具体的には eSPR、stNNI、Multiplier、および global Multiplier 移動）を表す 35 の連続特徴量を抽出する。これらの特徴量には、全局的な樹統計量、部分サブツリー情報、再配置経路の詳細、およびヒューリスティックな指標（例：「長い分岐の引き寄せリスク」）が含まれる。
   • モデル学習: ランダムフォレスト回帰モデルは、これらの特徴量に基づいて対数尤度の変化（$\Delta \log L$）を予測するように訓練される。訓練データは、高温および低温の事後領域の両方をカバーすることを保証する、アンニールされた SMC アルゴリズムのパイロット実行を通じて収集される。

3. 遅延受容 SMC（DA-SMC）:
   DA-MCMC カーネルは SMC 枠組みに埋め込まれる。アルゴリズムは、一連の中間ターゲット分布 $\pi_r(x) \propto p(x)[1(\Delta x > \kappa) \min\{p(D|x), f(x)\}]^{\phi_r}$ を通じて進化させる粒子の集団を維持する。

   • 伝播: 粒子は DA-MCMC カーネルを用いて伝播される。
   • 重み付け: 重みは、真の尤度が計算されていない場合に代理尤度を利用した中間ターゲットの比率に基づいて更新される。
   • リサンプリング: 相対有効サンプルサイズ（rESS）が閾値を下回った場合、粒子はリサンプリングされる。
   • 周辺尤度: この枠組みは、代理モデル強化された事後分布の周辺尤度に対する不偏推定量を提供し、モデル選択に使用できる。

主要な貢献

• アルゴリズム枠組み: 高価な尤度評価の頻度を特に削減するように調整された、ベイズ系統学のための遅延受容メカニズムと SMC の統合。
• 特徴量エンジニアリング: 機械学習予測のために樹の移動（eSPR、stNNI など）を特徴づけるための包括的な 35 のトポロジーおよび分岐長特徴量のセットの開発。
• 代理モデル化: 事前選択なしに高次元特徴量および非線形関係に対する頑健性を示す、尤度変化の予測へのランダムフォレストの適用。
• パラメータ選択戦略: 訓練段階において誤った却下率（第 2 種の誤り）を制御することにより、調整パラメータ $\delta$ を選択するためのデータ駆動型の手順。
• ソフトウェア実装: 提案された手法を実装する Python パッケージ（DAphyloSMC）の開発。

結果
著者らは、シミュレーションデータおよび実世界のデータセットの両方で手法を検証した：

• シミュレーション研究:
  • モデル選択: DA-SMC は、周辺尤度を推定することにより、候補（JC69、K2P、GTR）の中から真のデータ生成モデル（K2P）を正しく特定し、過小適合と過剰パラメータ化の両方を効果的に罰した。
  • パラメータ感度: 固定された総予算の下では、収束にとって粒子数（$K$）を増加させることよりも、反復回数（$R$）を増加させることがより重要であることが判明した。
  • 効率性: DA-SMC は標準的なアンニール SMC（ASMC）と比較して大幅な時間節約を達成し、速度向上は分類群の数が増加するにつれてより顕著になった。尤度評価ステップの並列化はさらに性能を加速した。
  • 精度: この手法は、尤度評価の回数を大幅に削減しながら、樹トポロジーの頑健な推定（Robinson-Foulds および Partition Metrics によって測定）を維持した。
• 実データ分析:
  • 3 つのデータセット（Primates、M336、M1809）に適用され、同一の計算予算の下で MrBayes と比較された。
  • DA-SMC は、樹距離指標（RF および PM）の面で MrBayes を一貫して上回り、より高い対数尤度および周辺尤度推定量を達成した。
  • M1809 データセットの場合、この手法は誤検出率を制御したまま、提案の 46.64% に対して正確な尤度計算を回避し、精度を犠牲にすることなく顕著な計算効率を示した。

意義と主張
本論文は、提案された DA-SMC アプローチが大規模なベイズ系統推論における計算ボトルネックに対する実用的な解決策を提供すると主張している。機械学習を活用してプロセスの初期段階で不適切な提案をフィルタリングすることにより、この手法は尤度評価に関連する計算負担を大幅に削減する。著者らは、彼らのアプローチが以下の点で優れていると主張している：

1. 探索の高速化: 樹空間を探索するために必要な尤度評価の回数を大幅に削減する。
2. 精度の維持: 事後分布推論およびモデル選択の能力の質を維持し、同一の計算予算内ではしばしば標準的な MCMC（MrBayes）を上回る。
3. 理論的保証の提供: 特定の条件下（完全な代理モデル）において、この手法は正確な事後分布を回復する。著者らは、SMC 推定量の一致性および漸近正規性、ならびに代理モデル強化分布に対する周辺尤度推定量の不偏性に関する理論的証明を提供する。
4. 柔軟性の提供: この枠組みはモジュール式であり、将来的には他の機械学習モデルや深層学習アーキテクチャとの統合が可能であり、SMC に固有の並列化をサポートする。

著者らは限界に関しては謙虚であり、代理モデルの不正確さが誤った受容（第 1 種の誤り）または誤った却下（第 2 種の誤り）につながる可能性を認め、これらは調整パラメータ $\delta$ によって管理されていると述べている。また、現在の実装は手作りの特徴量に依存しているが、将来の研究では樹空間のより洗練された表現を探求できる可能性があると指摘している。