跨模型MCMC算法的混合效率及其在貝葉斯系統發育學中的應用

Mixing efficiency of trans-model Markov chain Monte Carlo algorithms with applications in Bayesian phylogenetics

跨模型MCMC算法的混合效率及其在貝葉斯系統發育學中的應用

https://arxiv.org/pdf/2607.07188

跨模型馬爾可夫鏈蒙特卡洛（MCMC）算法廣泛應用于貝葉斯推斷，在貝葉斯系統發育學中尤為重要，因為在該領域中系統發育樹代表著不同的統計模型。盡管該算法具有極大的靈活性，但其混合效率差異巨大，且目前對其了解甚少。在此，我們運用數學分析與模擬方法，探討跨模型MCMC提議的混合效率，包括模型提議概率及模型參數的提議核。我們的分析印證了這一直觀認識：應當優先提議后驗概率較高的模型，并盡可能從后驗分布中提議參數值。我們的研究結果為構建高效的跨模型MCMC算法提供了指導原則。這些原則被應用于系統發育重建的MCMC算法中，并使用了靈長類和哺乳類的兩個真實數據集。

引言

自20世紀90年代引入分子系統發育學以來（Rannala and Yang, 1996; Yang and Rannala, 1997; Mau and Newton, 1997; Li et al., 2000），貝葉斯推斷已成為該領域最流行的統計方法之一（Chen et al., 2014; Yang, 2014）。許多復雜的分子序列進化模型已在流行程序中得以實現，例如 MRBAYES（Ronquist et al., 2012）、BEAST（Bouckaert et al., 2014）和 PHYLOBAYES（Lartillot et al., 2009），用于估計參數、比較序列進化模型以及推斷物種系統發育。系統發育重建是貝葉斯模型選擇的一個例子，因為系統發育樹（連同序列進化模型）規定了似然函數并對應不同的似然模型，而每棵樹中的分支長度（以及表征進化模型的參數）則是模型中的參數。計算是通過跨模型馬爾可夫鏈蒙特卡洛（MCMC）算法（Metropolis et al., 1953; Hastings, 1970; Green, 1995）實現的，該算法在樹之間移動，且 MCMC 訪問每棵樹的頻率即是其后驗概率的估計值。

已有研究指出，在現代系統基因組數據集中，推斷出的樹或分支（clades）的后驗概率通常約為 ~100%（例如，Yang and Zhu, 2018），這促使一些研究人員質疑當數據集很大且模型被誤設時（例如，Thomson and Brown, 2022），樹的后驗概率的實用性。這些都是重要的問題，序列進化模型的選擇以及先驗的設定和影響（Yang, 2014）也是如此，但這些超出了本文的范圍。在此，我們專注于當數據、先驗和模型固定時，從后驗分布進行 MCMC 采樣的效率。

針對跨模型 MCMC 移動，已引入了多種構建方法。在積空間（product-space）構建（Carlin and Chib, 1995）中，馬爾可夫鏈的狀態包括所有模型中的所有參數。當鏈處于一個模型中時，其他模型的參數對似然或后驗沒有影響，但被視為偽參數（pseudo-parameters）并使用偽先驗（pseudo-priors）進行更新。Green（Green, 1995）提出了維數匹配（dimension-matching）的概念，以允許在可逆跳躍 MCMC（rjMCMC）算法中在不同大小的模型之間進行移動。復合空間（composite-space）框架（Godsill, 2001）允許參數在模型之間任意重疊，并將 rjMCMC 和積空間構建作為特例包含在內。飽和空間（saturated-space）方法（Brooks et al., 2003）與此類似，它擴充了“小”模型的狀態空間，以達到與“最大”模型相同的維數。不同的跨模型采樣方案在（Dellaportas et al., 2002）中進行了比較。幾種其他的跨模型推斷算法（Grenander and Miller, 1994; Stephens, 2000; Cappé et al., 2003）可被視為 rjMCMC 的特定版本。對于模型內推斷問題，許多作者討論了

利用后驗的局部信息來指導提議（proposals）的好處（Zanella, 2020）。例如，局部梯度可用于引導提議朝向高概率區域，如在 Metropolis 調整的 Langevin 算法（MALA, (Roberts and Rosenthal, 2008)）和哈密頓蒙特卡洛（HMC）（Neal, 2011; Girolami and Calderhead, 2011）中。最近，HMC 正被改編用于提議模型間的移動（Nishimura et al., 2020）。

在貝葉斯系統發育學中，在簡約性和似然樹搜索中開發的分支交換算法，如最近鄰互換（NNI）和子樹修剪與重接（SPR）（Swofford et al., 1996），已被改編為 MCMC 提議算法，以向當前樹引入局部的隨機改變（Lakner et al., 2008; Hohna et al., 2008; Yang, 2014）。人們已致力于設計局部知情（locally-informed）的 MCMC 移動，相比于以均勻概率選擇候選樹的盲目提議，這些移動能更好地反映后驗分布。例如，可以優先選擇較短的內部分支進行樹擾動（Rannala and Yang, 2017），并根據所得樹的簡約性得分（Yang, 2014; Zhang et al., 2020）、分支的條件概率（Hohna and Drummond, 2012）或在預運行（pilot run）期間估計的分割（splits，即由樹上的內部分支定義的物種二分法）的后驗概率（Meyer, 2021），使用權重對重接的目標分支進行采樣。

我們在本文中探討這些問題。在構建關于比較兩個模型（其中一個模型包含參數）這一一般情況的理論之前，我們首先分析了兩個跨模型 MCMC 的簡單示例，其目標分布分別為均勻分布和正態分布。隨后，我們研究了用于隨機樹搜索的 SPR（子樹修剪與重接）算法的幾種變體，以通過這些簡單示例來闡釋所發展的理論。我們重點關注跨模型 MCMC 算法的兩個特征：(i) 模型跳躍概率
以及 (ii) MCMC 樣本用于估計模型后驗概率的效率。

跨模型 MCMC 算法概述

跨模型 MCMC 算法

模型跳躍概率

我們將模型跳躍概率定義為在模型和參數的目標分布上的平均值：

均勻分布示例

高斯分布示例

高斯示例的一般情況

6 MCMC 混合效率

我們利用這一理論來闡釋一個常見的觀察結果：隨著數據量的增加，跨模型算法傾向于變得越來越低效，這是因為新模型中參數的提議并未被調整以緊密匹配后驗分布。當數據量增加時，模型內的后驗分布變得更加集中，而為新模型提議的參數值越來越可能錯過后驗分布的眾數（mode），從而導致提議的模型被拒絕，即使它比當前模型具有更高的后驗概率。

模型提議概率

含參數的兩模型情形下的最優馬爾可夫鏈

系統發育學中的 SPR 提議

基線算法包含一個 SPR 移動（圖 6b&b'）用于更新樹結構，以及一個樹內移動（within-tree move）用于更新分支長度，后者涉及隨機采樣一個分支并通過乘數對其進行改變。我們考慮了五種用于更新分支長度的樹內移動變體，包括：A1，使用一個 (2n-3) 維移動來改變所有分支長度；A2-1，使用一系列一維（1-D）移動來改變分支長度；A2-2，使用一系列一維移動并更新樹長（即所有分支長度之和）；A3-1，使用一系列一維 Bactrian-Laplace 移動來改變分支長度（Yang and Rodriguez, 2013）；以及 A3-2，使用一系列一維移動并利用 Bactrian-Laplace 提議進行樹長更新。詳見 SI 擴展方法。雖然不同的樹內提議預計會影響分支長度的混合（mixing），但它們對跨樹算法（cross-tree algorithm）的混合效率影響甚微（表 2&3）。

我們探索了不同的策略來提議替代的樹拓撲結構并為新樹生成分支長度，包括：B1 (NNI?)，直接轉移分支長度的 NNI；B2 (NNI?)，修改焦點分支（focal branch）長度的 NNI；B3 (NNI?)，使用乘數修改焦點分支周圍所有五個分支長度的 NNI；B4 (NNI_bw)，帶分支權重的 NNI；B5 (NNI_bw&pw)，帶分支權重和簡約權重（parsimony weights）的 NNI；B6 (NNI_local)，Larget 和 Simon (1999) 的局部移動；B7 (SPR_bw)，帶分支權重的 SPR；B8 (SPR_bw&pw)，同時帶分支權重和簡約權重的 SPR。詳見 SI 擴展方法。

對基于 NNI 的移動的比較表明，在樹之間移動時最好不要改變分支長度，盡管差異很小。B1 (NNI?)、B2 (NNI?) 和 B3 (NNI?) 都使用分支長度的直接轉移，其表現遠優于 B6 (NNI_local)，后者通過合并和拆分分支來形成新的分支長度。例如，對于 ψ/η 數據集，B1 (NNI?) 的效率是 B6 (NNI_local) 的 20 倍 (= 0.1778/0.0089)；對于 mt 數據集，效率高 16.6 倍 (= 0.0053/0.00032)。此前已指出，直接轉移產生的分支長度更接近最大似然估計 (MLE)，因此優于合并-拆分方法 (Yang, 2014, p.283-4; 另見圖 S1)。NNI 變體 B1-B5 的表現也遠優于 SPR（基線算法和 B7 SPR_bw）。例如，對于 ψ/η 數據集，B1 (NNI?) 的效率是基線算法的 21 倍 (= 0.1778/0.0086)；對于 mt 數據集，效率高 20 倍 (= 0.0053/0.00027)。這些差異歸因于兩個因素。首先，如上所述，NNI 變體 B1-B5 使用直接轉移，而 SPR 使用合并與拆分來為新樹生成分支長度。其次，NNI 是比 SPR 更“小”的移動，因此提議的樹中有更高比例屬于高概率樹，從而導致更高的模型跳躍率。

基線算法包含一個 SPR 移動（圖 6b&b'）用于更新樹結構，以及一個樹內移動（within-tree move）用于更新分支長度，后者涉及隨機采樣一個分支并通過乘數對其進行改變。我們考慮了五種用于更新分支長度的樹內移動變體，包括：A1，使用一個 (2n-3) 維移動來改變所有分支長度；A2-1，使用一系列一維（1-D）移動來改變分支長度；A2-2，使用一系列一維移動并更新樹長（即所有分支長度之和）；A3-1，使用一系列一維 Bactrian-Laplace 移動來改變分支長度（Yang and Rodriguez, 2013）；以及 A3-2，使用一系列一維移動并利用 Bactrian-Laplace 提議進行樹長更新。詳見 SI 擴展方法。雖然不同的樹內提議預計會影響分支長度的混合（mixing），但它們對跨樹算法（cross-tree algorithm）的混合效率影響甚微（表 2&3）。

我們探索了不同的策略來提議替代的樹拓撲結構并為新樹生成分支長度，包括：B1 (NNI?)，直接轉移分支長度的 NNI；B2 (NNI?)，修改焦點分支（focal branch）長度的 NNI；B3 (NNI)，使用乘數修改焦點分支周圍所有五個分支長度的 NNI；B4 (NNI_bw)，帶分支權重的 NNI；B5 (NNI_bw&pw)，帶分支權重和簡約權重（parsimony weights）的 NNI；B6 (NNI_local)，Larget 和 Simon (1999) 的局部移動；B7 (SPR_bw)，帶分支權重的 SPR；B8 (SPR_bw&pw)，同時帶分支權重和簡約權重的 SPR。詳見 SI 擴展方法。

對基于 NNI 的移動的比較表明，在樹之間移動時最好不要改變分支長度，盡管差異很小。B1 (NNI?)、B2 (NNI?) 和 B3 (NNI?) 都使用分支長度的直接轉移，其表現遠優于 B6 (NNI_local)，后者通過合并和拆分分支來形成新的分支長度。例如，對于 ψ/η 數據集，B1 (NNI?) 的效率是 B6 (NNI_local) 的 20 倍 (= 0.1778/0.0089)；對于 mt 數據集，效率高 16.6 倍 (= 0.0053/0.00032)。此前已指出，直接轉移產生的分支長度更接近最大似然估計 (MLE)，因此優于合并-拆分方法 (Yang, 2014, p.283-4; 另見圖 S1)。NNI 變體 B1-B5 的表現也遠優于 SPR（基線算法和 B7 SPR_bw）。例如，對于 ψ/η 數據集，B1 (NNI?) 的效率是基線算法的 21 倍 (= 0.1778/0.0086)；對于 mt 數據集，效率高 20 倍 (= 0.0053/0.00027)。這些差異歸因于兩個因素。首先，如上所述，NNI 變體 B1-B5 使用直接轉移，而 SPR 使用合并與拆分來為新樹生成分支長度。其次，NNI 是比 SPR 更“小”的移動，因此提議的樹中有更高比例屬于高概率樹，從而導致更高的模型跳躍率。

使用分支權重優先改變短內部分支周圍的樹結構，提高了 NNI 和 SPR 的跨樹混合效率。在 ψ/η 數據集中，B4 (NNI_bw) 的效率是 B1 (NNI?) 的 2.2 倍（=0.3969/0.1778），盡管在 mt 數據集中這種影響要小得多（0.0057/0.0053 = 1.08）。同樣，在 ψ/η 數據集中，B7 (SPR_bw) 比 B0（SPR 基線）效率高 2.0 倍（=0.0172/0.0086），而在 mt 數據集中影響較小（0.00028/0.00027 = 1.04）。數據集之間的差異可能是因為 ψ/η 樹（圖 7a）中內部分支長度存在巨大差異，而 mt 數據集中的分支長度則更加均勻（圖 7b）。請注意，我們的分支權重（公式 S90）是任意的，可能更適合一個數據集而不是另一個。

使用簡約得分（公式 S91）來采樣目標分支提高了混合效率。在 NNI 移動中，對于 ψ/η 數據集，B5 NNI_bw&pw 比 B4 NNI_bw 效率高 2.1 倍，但在 mt 數據集中沒有影響。對于 SPR 移動，在 ψ/η 數據集中，B8 SPR_bw&pw 比 B7 SPR_bw 效率高 3.5 倍（= 0.0608/0.0172），在 mt 數據集中高 1.7 倍（= 0.00047/0.00028）。

總體而言，使用這兩個數據集的測試證實了我們的理論分析。在 NNI 移動中，直接轉移分支長度比合并-拆分方法具有更高的效率，這證實了為新模型提議接近模型內后驗眾數（mode）的模型參數（分支長度）的重要性。內部分支的分支權重和目標分支的簡約權重都改變了模型提議概率（q_kk'），傾向于那些接近當前樹且可能具有更高后驗概率的樹。根據我們的理論分析，這些策略預計會提高跨模型算法的混合效率，并且發現在兩個數據集中都對跨樹混合效率有重大影響。

我們使用 B4 (NNI_bw) 作為跨樹移動，來研究樹內移動和跨樹移動之間的最佳分工。對于樹內移動，我們使用了基線方法（Baseline），即隨機改變一個分支長度。對于 ψ/η-globin 數據集，最佳的跨樹提議概率為 0.7–0.8，而對于 mt 數據集，約為 0.5（圖 8a&b）。

樹內移動和跨樹移動之間計算工作量的分配

與其在每次 MCMC 迭代中都進行一次樹內移動和一次跨樹移動，不如這里我們以概率 p? 采樣跨樹移動（以 p? = 1 - p? 采樣樹內移動）。對于 ψ/η 數據集，最佳的 p? 約為 0.75，對于 mt 數據集為 0.55（圖 8）。

模型內（within-model）與跨模型（trans-model）MCMC 算法之間的差異

模型內算法和跨模型算法之間存在許多差異（Yang, 2014; Nascimento et al., 2017）。首先，對于模型內移動，可以將步長（例如，滑動窗口提議的窗口大小）設置得足夠小，使得接受率約為 100%。然而，對于跨模型移動，不存在這種微小步長的概念，此外接受率受到后驗模型概率的約束（公式 7）。如果 MAP（最大后驗概率）模型具有 99% 的后驗概率，P_jump 不能超過 2(1 - 0.99) = 2%，否則鏈將不會足夠頻繁地訪問 MAP 模型以達到正確的后驗概率。其次，對于模型內移動，中間接受率（30-40%）是最佳的（Gelman et al., 1996）。對于跨模型移動，活躍鏈（具有高 P_jump）通常比懶惰鏈（具有低 P_jump）更有效，我們應該努力實現高接受率，盡管我們的理論確定了最大 P_jump 既非最大混合效率的必要條件也非充分條件的情況。雖然對于模型內算法，接受率約為 0 幾乎總是表明存在混合問題（例如，窗口大小可能太大），但這對于跨模型移動不一定意味著混合問題，因為它可能是由于極端的后驗概率造成的。雖然 MCMC 理論經常強調算法對復雜參數空間的普遍適用性，但在我們的目標是開發高效的跨模型算法時，應考慮到模型內和跨模型算法之間的差異。

當前系統發育程序中使用的某些提議是樹內移動和跨樹移動的混合（例如，B6 NNI_local），這使得調整步長以改善性能變得尷尬。例如，Lakner 等人（Lakner et al., 2008）發現總體接受率并不是算法效率的良好指標，而拓撲結構變化的接受率則是。這顯然是因為（Lakner et al., 2008）中的接受率是對樹內移動（其中中間接受率最佳）和跨樹移動（其中高接受率是可取的）的平均，因此總體接受率并不是一個非常有用的性能指標。

系統發育中 MCMC 算法混合效率的度量

除了基于方差比的效率度量（公式 11）外，（Lakner et al., 2008）和（Hohna et al., 2008）使用的一種替代度量是基于分割（splits）的概率。運行一個極長的“參考鏈”以收集分割及其“真實”概率。然后運行測試 MCMC 算法固定的迭代次數（N），并計算“真實分割概率”（p_i 對應分割 i）與測試鏈的估計值（p?_i）之間的距離。例如，距離可以定義為分割概率的最大差值，

原文鏈接：https://arxiv.org/pdf/2607.07188