跨模型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