无似然马尔可夫链蒙特卡罗耦合的高效近似贝叶斯计算

摘要

近似贝叶斯计算（ABC）技术允许在复杂的人口模型中进行推断，但计算效率低下。提出了马尔可夫链蒙特卡罗（MCMC）方法(马约拉姆 等。2003)但它存在计算问题和较差的混合。我们提出了几种方法学的发展，以克服这种MCMC方法的缺点，从而实现比标准ABC更大的计算进步。其主要思想是放宽MCMC内的容差，以允许良好的混合，但通过对输出进行二次采样和回归调整的组合，保持对后验值的良好近似。我们还建议使用偏最小二乘（PLS）变换来选择信息统计。我们的方法的准确性在有迁移和无迁移的两个种群的发散情况下进行了检验。在这种情况下，我们的ABC–MCMC方法需要相当低的计算时间才能达到与传统ABC相同的精度。然后，我们将我们的方法应用于一个更复杂的案例，估计三个非洲人口之间的分化时间和迁移率。

随着大规模基因分型技术的出现(例如,绿色 等。2006;征收 等。2007)，基因数据可以以前所未有的规模产生(阿尔茨舒勒 等。2005;巴斯塔曼特 等。2005)现在可以在基因组上的数百个位点上对个体和群体的遗传变异性进行常规检测(罗森博格 等。2002;威廉姆森 等。2005;贝克特和普热沃尔斯基2007;弗雷泽 等。2007). 这些大型数据集为更好地理解影响包括人类在内的许多物种多样性的进化力量提供了希望，也为确定过去选择性事件中涉及的基因组区域提供了希望(阿尼西莫娃和利伯莱斯2007;尼尔森 等。2007). 然而，需要考虑人口的人口历史，以将其影响与选择的影响区分开来(哈德迪尔 等。2005;尼尔森 等。2005;比斯瓦斯和阿基2006). 因此，能够从中性遗传数据中正确估计过去的人口统计学数据或同时估计人口统计学和选择似乎很重要(例如,威廉姆森 等。2005). 在过去10年中，突变和人口统计学参数的统计估计有了显著改进，尤其是使用了贝叶斯方法和全似然方法(博蒙特 等。2002;马约拉姆和塔瓦雷2006). 然而，这些方法仍然局限于可以计算其可能性的相对简单的模型，或者局限于可以在合理的时间内分析的小数据集。在实际模型下处理大数据集和估计人口统计参数仍然存在问题，在这些情况下经常使用“良好”方法（参见，例如,马思 等。2004;沙夫纳 等。2005;Plagnol公司和墙壁2006).

近似贝叶斯计算（ABC）框架(塔瓦雷 等。1997;普里查德 等。1999;博蒙特 等。2002)基于简单拒绝算法的，已应用于非模式生物和人类各种进化模型中人口统计学参数的估计(埃斯图 等。2004;塔尔蒙 等。2004;Excoffier公司 等。2005年a;希克森 等。2006;法贡德斯 等。2007;罗森布拉姆 等。2007). 一般原则首先概述于塔瓦雷 等。(1997)就是模拟数据(D类′）与观察结果类似(D类)给定模型下的样本大小和位点数量，带参数(θ)从一些先前的分布中提取。如果D类'与相同D类参数被存储，否则被丢弃，保留的参数用于估计后验分布。因为它不太可能模拟D类'与相同D类对于大型数据集或复杂模型(马约拉姆和塔瓦雷2006)，已提出(普里查德 等。1999)用一组汇总统计数据替换数据秒并且如果模拟的摘要统计信息，则保留特定的模拟S′与观察到的汇总统计数据非常接近秒.说明两者之间的差异秒和S′,博蒙特 等. (2002)最近提出执行局部加权线性回归来计算后验分布。经证明，该回归调整步骤可大大提高估计值。然而，即使有了这种改进，ABC方法的计算效率也不是很高，因为它需要模拟数百万个样本，其中大多数样本（通常≥99%）将被丢弃用于参数估计。

最近，马约拉姆 等。(2003)提出了另一种类似的无障碍方法，其中模拟直接嵌入到马尔可夫链蒙特卡罗（MCMC）框架中D类'如果它们与观测数据相等，则可以接受D类否则将收敛到右后向分布

⁠对于接受率太小的复杂数据集，他们建议用汇总统计代替全部数据，并在与数据足够接近的情况下接受新的参数值，就像传统ABC方法一样。与此方法相关的一个问题是定义接受模拟需要多大程度的接近，这些条件（i）接受率，（ii）链的混合和收敛，以及（iii）老化期，因为如果起点在具有低可能性的区域中，则可能需要非常大量的模拟才能具有第一个可接受的步骤。其他类似的无罩方法，如顺序蒙特卡罗(西森 等。2007)或人口蒙特卡罗（M.A。博蒙特、J.-M。科尔尼埃、J.-M。马林和C.P。罗伯特，未发表的结果），最近提出解决收敛速度慢的问题，并更有效地探索多模态后验。然而，这些无需任何防护罩的方法很少用于复杂的进化模型（但请参阅拉特曼 等。2007用于比较蛋白质网络的进化动力学）。

在本文中，我们提出了一种新的方法，借鉴了基于拒绝算法的传统ABC和无似然MCMC的最佳特性，有效地处理了上述所有问题。我们在不同的人口隔离和迁移模型的情况下测试了我们的方法尼尔森和韦克利(2001)，但不要试图将我们的方法与完全似然方法进行比较。我们的新方法最终应用于三个非洲人口的人口隔离和迁移模型。

方法

我们首先描述马约拉姆 等。的（2003）基于汇总统计（SS）的无相似性MCMC算法，以下称为SS–MCMC，显示了其潜在问题。然后，我们提出了从传统ABC中借用的解决方案，这导致了一种称为ABC–MCMC的新算法。

SS–MCMC算法：

给出一些观察数据D类在给定模型下生成M（M）由一组（未知）参数定义θ事先分配

⁠,马约拉姆 等。(2003)已经表明后验分布

\（f（mathbf{mathrm{{theta}}}{vert}D）

可以从马尔可夫（M）链的样本中获得，而无需通过以下算法生成似然：

• M1.提议脱离当前状态θ到一个新的状态θ′根据转换内核

  \（q（\mathbf{\mathrm{{theta}}}{\rightarrow}\mathbf{\mathbm{{theta}}{^\prime}）
  ⁠.

• M2.模拟数据D类'模型下M（M），使用新参数θ′.

• M3.如果D类′ =D类，转至步骤M4；否则保持状态θ并返回步骤M1。

• M4.接受状态θ′有可能

  \（h（\mathbf{\mathrm{{theta}}}{\rightarrow}\mathbf{\mathm{{theta}}{^\prime}}）{=}\mathrm{min}{\mathrm{{theta}}）/\mathrm{{\pi}}
  ⁠; 否则保持状态θ.转至步骤M1。

因为它不太可能模拟D类'与相同D类对于大型数据集或复杂模型(马约拉姆和塔瓦雷2006),马约拉姆 等。(2003)建议替换数据D类通过一组足够的汇总统计数据秒和条件D类′ =D类根据不太严格的条件

在步骤M3中，其中

\（\mathrm{{\delta}}{\mathrm{{\varepsilon}}}\）

是任意较小的距离秒'和秒请注意，如果摘要统计信息不足，即统计信息无法捕获参数数据中包含的全部信息θ，得到的后验值只是真实后验分布的近似值(马约拉姆和塔瓦雷2006). 使用此新步骤M3，马约拉姆 等。(2003)声称该SS–MCMC链的平稳分布为

\（f（\mathbf{\mathrm{{theta}}}{\vert}\mathrm{{delta}}{leq}\mathr{{delta}}{{\mathr m{{varepsilon}}）\）

⁠，这应该是真实后部的良好近似值

\（f（mathbf{mathrm{{theta}}}{vert}D）

如果

\（\mathrm｛\delta｝｝_｛\mathrm｛\varepsilon｝｝｝｝\）

很小(马约拉姆 等。2003).

SS–MCMC算法的问题：

摘要统计信息和阈值的使用

在修改后的步骤M3中，介绍了我们在下面解决的一些问题。

1. 因为马尔可夫链只能在以下条件下开始

   \（\mathrm{{\delta}}{\leq}\mathrm{{\delta}}_{\mathrm2{{\varepsilon}}}\）
   ⁠，首先满足此条件所需的步骤数未定义，它取决于δ的（未知）分布，
   \（\mathrm{{\pi}}（\mathr{{\delta}}{\vert}\mathbf{\mathrm{S}}，\mathbf{\mathr}{\theta}}）
   如下所示。

2. 阈值的选择

   \（\mathrm{{\delta}}{\mathrm{{\varepsilon}}}\）
   是很重要的，因为过大的公差间隔会导致链由前一个控制
   \（\mathrm{{\pi}}（\mathbf{\mathrm{{\theta}}}）
   ⁠，可能性不大θ将经常被接受。另一方面，值太小会导致接受率很低。此外，链的接受率与θ(西森 等。2007)，使其在低可能性区域非常粘滞，防止使用太小的阈值
   \（\mathrm{{\delta}}{\mathrm{{\varepsilon}}}\）
   ⁠.

3. 最后，估计的准确性取决于汇总统计数据和距离函数的选择(汉密尔顿 等。2005). 在无法计算模型可能性的复杂模型中，很难为所有参数找到足够的统计数据，因此，定义一组最佳汇总统计数据是一个重要且尚未解决的问题。

对SS–MCMC算法的修改：

校准步骤：

我们建议通过执行一系列n个模拟（我们通常使用n个=10000），其中每次参数都是从其获得之前随机抽取的

⁠，经验近似值

\（\mathrm{{\pi}}（\mathr{{\delta}}{\vert}\mathbf{\mathrm{S}}，\mathbf{\mathr}{\theta}}）

⁠通过此校准步骤，我们可以方便地定义公差水平ε和阈值距离

\（\mathrm{{\delta}}{\mathrm{{\varepsilon}}}\）

这样的话

\（P（\mathrm{{delta}}{\leq}\mathrm{{delta}}{{\mathrm2{{varepsilon}}}）{=}\mathr{{varesilon}}

⁠例如，通过设置ε=0.01，我们定义了一个值

\（\mathrm{{delta}}{0.01}\）

低于该条件

\（{\Vert}\mathbf{\mathrm{S}{^\prime}}{-}\mathbf{\mathrm{S}}{\Vert}{\leq}\mathrm{{\delta}}{{\mathm{{\varepsilon}}}）

对于所有随机模拟的数据集的1%将为真。然后我们应该能够使用任何模拟

\（\mathrm{{\delta}}{\leq}\mathrm{{\delta}}_{\mathrm2{{\varepsilon}}}\）

作为链条的起点。这些n个仿真还用于调整建议范围，从而调整过渡核q个将马尔可夫链中新参数值的建议范围设置为统一的宽度范围

\（\mathrm{{\varphi}}\）

以标准偏差单位表示，在n个ε保留模拟。如果提案范围过大，链可能会跳转到以下状态

\（\mathrm{{\delta}}{>}\mathrm{{\delta}}_{\mathrm2{{\varepsilon}}}\）

⁠这导致接受率较低，因此有效样本数量较少。如果建议范围太小，链可能永远不会探索整个参数空间，因此可能导致重复之间的估计参数存在较大差异。注意，旨在增加混合的标准MCMC技术，如模拟回火(博尔托特 等。2007)，可以很容易地包含在这里介绍的框架中。

结合SS–MCMC和ABC：

我们建议通过推出SS–MCMC长链来解决第二个问题（接受率低）秒公差相对较大(即, ε = 0.01). 如所示附录，平稳分布

该SS–MCMC链与通过简单拒绝算法（如传统ABC）获得的相同，使用相同的公差水平ε。因此，我们建议对大小为的子样本进行参数估计t吨由与马尔可夫链生成的最小距离δ相关的模拟组成，正如通常在ABC框架中所做的那样(普里查德 等。1999;博蒙特 等。2002)，并执行局部回归调整，其中样本由其关联的

\（\mathrm{{delta}}{i}\）

-值(博蒙特 等。2002)，大大提高了估计的质量。这种方法允许链具有较大的接受率，同时由于回归调整步骤，使最终估计对先验不太敏感(博蒙特 等。2002). 因此，这种两步方法应能得到与具有容差的简单拒绝算法相同的平稳分布

\（\mathrm{{\varepsilon}}{^\prime}{=}\mathrm{{\varepsilon}}（t/s）\）

⁠因此，与简单的拒绝取样器相比，SS–MCMC的主要增益是需要更少的模拟才能获得t吨容差ε′下的样品：秒+n个SS–MCMC下的仿真与秒/ABC下的ε模拟，意味着计算成本在理论上减少了一个因素

\（s{[}\mathrm{{\varepsilon}}（s{+}n）{]}^{{-}1}\)

⁠例如，校准步骤为n个=10000，一串秒=90000步，公差ε=1%，预期收益是计算成本的90倍。因此，我们预计，在改进的SS–MCMC框架下，总共有100000个模拟将对应于传统ABC框架下的900万个模拟，但这种增益存在实际限制，这是运行小ε值链的固有困难（见下文）。在下文中，我们使用术语ABC–MCMC来指定对SS–MCMC链输出的给定部分进行ABC回归调整的过程。

汇总统计的偏最小二乘变换：

我们建议通过使用偏最小二乘（PLS）回归方法来解决信息汇总统计的选择问题（参见，例如,特尼豪斯 等。1995;布莱斯泰和Strimmer公司2007). 与主成分分析（PCA）一样，PLS从高维数据集中提取正交成分X（X）但除此之外，通过最大化预测变量和响应变量的协方差矩阵，选择这些分量来适当解释响应变量的可变性（参见，例如,布莱斯泰和Strimmer公司2007). 在当前ABC上下文中，预测变量是原始汇总统计数据，响应变量是模型参数。PLS组件数量的选择通常基于离开验证程序(梅维克和韦伦斯2007)，检查回归预测参数的均方根误差（RMSE）。PLS分析的结果是，与大量潜在相关汇总统计的初始集合相比，独立分量的数量应大大减少，其中一些与任何参数几乎无关，因此只会给欧几里德距离增加噪声。PLS变换的另一个优点是，它保证了PLS变换后的汇总统计矩阵是非奇异的，这是在执行最终局部加权线性回归估计参数时所必需的(博蒙特 等。2002). 实际上，我们建议在每个模拟和观测数据集上计算一组相对较大的汇总统计数据（假设为参数的信息）。这个n个在校准步骤中进行的随机模拟被用于计算PLS分量，然后PLS分量被用于变换在马尔可夫链期间生成的模拟数据集上计算的汇总统计。我们使用免费提供的R包“PLS”中的例程“plsr”计算PLS组件，并选择一组最佳的k个部件符合免验证程序(梅维克和韦伦斯2007). 由于PLS变换假设参数和统计数据之间存在线性关系，我们首先应用了多元Box–Cox变换(盒子和考克斯1964)在定义PLS组件之前，分别对每个统计数据进行分析。请注意，可以想象出选择适当的汇总统计数据的其他方法，如根据其包含是否大大提高了推论的质量对其进行评分，正如最近提出的那样(乔伊斯和马约拉姆2008).

ABC–MCMC算法：

我们在这里描述了ABC–MCMC算法（AM），其中包括与无可能性的普通MCMC方法相比提出的改进：

• AM1.执行n个带参数的仿真θ′从他们的先验中随机抽取，每次计算他们相关的汇总统计数据集S′.

• AM2.根据n个 θ′和S′统计的Box-Cox变换后的向量。

• AM3.适用于所有人n个仿真，转换摘要统计信息S′进入之内k个保留的PLS组件，如S′_请.转换观察到的摘要统计信息秒作为

  \（\mathbf{\mathrm{S}}_{\mathlm{PLS}}\）
  并计算
  \（p_{n}（\mathrm{{delta}}{\vert}\mathbf{\mathrm{S}}_{\mathr{PLS}}，\mathbf{\mathrem{{theta}}）
  ⁠.

• AM4.固定ε，估计

  \（\mathrm{{\delta}}{\mathrm{{\varepsilon}}}\）
  从
  \（p_{n}（\mathrm{{delta}}{\vert}\mathbf{\mathrm{S}}_{\mathr{PLS}}，\mathbf{\mathrem{{theta}}）
  ⁠，并设置转换内核参数的建议范围
  \（q（\mathbf{\mathrm{{theta}}}{\rightarrow}\mathbf{\mathbm{{theta}}{^\prime}}）
  基于φ和参数之间的可变性
  \（n\mathrm{{\varepsilon}}\）
  保留仿真。

• AM5.启动总长度的MCMC链秒从一个位置θ从中随机选择

  \（n\mathrm{{\varepsilon}}\）
  最接近的仿真D类.设置我= 0.

• AM6.如果现在在θ，建议动议

  \（\mathbf{\mathrm{{theta}}{^\prime}}\）
  根据转换内核
  \（q（\mathbf{\mathrm{{theta}}}{\rightarrow}\mathbf{\mathbm{{theta}}{^\prime}}）
  ⁠.增量我.

• AM7.模拟D类'基于

  \（\mathbf{\mathrm{{theta}}{^\prime}}\）
  ⁠.计算汇总统计S′并将其转换为S′_请.

• AM8.如果

  \（\mathrm{{delta}}{{i}{=}{\Vert}\mathbf{\mathrm{S}{^\prime}}{\mathrm{PLS}}{-}\mathbf{\mathrm{S}}}{
  ⁠，停留在θ然后转到AM6。

• AM9.接受

  \（\mathbf{\mathrm{{theta}}{^\prime}}\）
  有可能
  \（\mathrm{min}（1，\mathrm{{pi}}（\mathbf{mathrm}{theta}}{^\prime}}）q \mathrm{{theta}}}{\rightarrow}\mathbf{\mathrm{{theta}}{^\prime}}））
  ⁠; 否则留在θ.

• AM10.如果我<秒，转到AM6。

• AM11.来自秒链的样本，保留大小为的子样本t吨由关联最小的仿真组成

  \（\mathrm{{delta}}{i}\）
  -值并丢弃其他值标准时间样品。

• AM12.执行ABC回归调整(博蒙特 等。2002)在上t吨保留样本以估计参数。

为了防止链条在步骤AM5后停留在给定的起始位置，我们选择一个新的初始值θ-如果链未移动到新值，则返回值

在20个提案之后。还要注意，我们在步骤AM6中同时更新所有参数。如上所述，建议范围的宽度针对每个参数进行独立调整，并设置为根据

\（n\mathrm{{\varepsilon}}\）

步骤AM4中保留的仿真。以下内容博蒙特 等。(2002)，我们使用欧几里德范数作为距离函数‖·‖，但没有标准化

\（\mathbf{\mathrm{S}}_{\mathlm{PLS}}\）

⁠.

并行化ABC–MCMC：

在传统ABC方法下执行的模拟可以很容易地分布在多个CPU之间。因此，我们实现了SS–MCMC算法的并行化版本，如下所示：在10000个随机模拟（可以很容易地并行化）的初始校准阶段之后，我们用相同的ε-和φ值。然后将10条链连接起来，并使用局部回归调整来估计5000个“最佳”模拟的后验密度(博蒙特 等。2002)完成ABC–MCMC算法。通过这种方法，100000个仿真中的90%可以分布在10个不同的CPU上，而更多的仿真可以分布在更多的CPU上。

插图和应用：

我们测试了ABC–MCMC方法，并将其应用于人口迁移差异模型(尼尔森和韦克利2001). 在这个模型中，我们假设T型过去的几代人，一个庞大的祖先群体N个_A类分裂成两个大小的种群N个₁和N个₂然后，移民在两个人口之间以一定的速度交换（展望未来）米₁₂和米₂₁。该方法已应用于几个不同的数据集(例如,贝克特和普热沃尔斯基2007;霍尔泽尔 等。2007). 必须指出的是，它是最复杂的种群遗传模型之一，可以完全实现该模型(嘿和尼尔森2007;库纳2009). 然而，我们不打算在这里将我们的方法与完全似然法进行比较，因为本文的目的是改进传统的ABC方法，而不是与全似然法竞争。

我们首先在一个没有迁移的简单案例中比较了传统ABC方法和新的ABC–MCMC方法(米₁₂=米₂₁= 0). 所有参数均来自统一的先验信息：每个群体的基因拷贝数N个∼单位[030000]和发散时间T型∼单位[016000]代。我们用SIMCOAL 2.0程序模拟遗传数据(拉瓦尔和Excoffier公司2004)每个群体25个二倍体个体，每个群体50个非连锁微卫星。在纯逐步突变模型下模拟微卫星数据。而我们使用固定的突变率μ=5×10⁻⁴在这种情况下，允许突变率在应用案例中发生变化，如下所述。共有100个测试数据集，其参数值从先前分布中随机抽取，用于比较ABC和ABC–MCMC。然后，我们将传统ABC方法和新的ABC–MCMC方法在更复杂的情况下与在单位[0, 0.003]. 同样，通过从先前分布中提取参数值，对100个测试数据集进行了随机模拟。

汇总统计：

使用软件包Arlequin 3.1(Excoffier公司 等。2005年b)，我们计算了每个群体中等位基因数量的平均偏差和标准偏差(K)等位基因大小的范围(R（右）)，预期杂合度(H（H）)，Garza–Williamson统计(加尔扎和威廉姆森2001)修改为GW=K/(R（右）+ 1) (Excoffier公司 等。2005年a)以及计算为GW*的GW的另一个修改=K/(R（右）_总计+1），其中R（右）_总计是在所有抽样人群中计算的等位基因大小范围。使用GW*统计数据背后的想法是，它应该反映特定人群的漂移效应，因为分母在所有人群中都是相同的。此外，还计算了两个集合人口的相同统计数据及其标准偏差，但GW*除外，因为在这种情况下GW=GW*。我们还计算了分化指数F类_装货单和遗传距离（δμ）²(戈尔茨坦 等。1995)在两个群体之间。因此，我们计算了所有观测和模拟数据集的总计31个汇总统计数据。根据校准步骤中进行的10000次模拟，从汇总统计数据中提取PLS成分。我们使用了R包PLS(梅维克和韦伦斯2007)找到要使用的PLS组件的适当数量（有迁移和无迁移的情况下都是10个）。

测量方法的准确性：

由于模拟的计算要求比通过回归参数的汇总统计获得的后验分布的估计要高得多(Excoffier公司 等。2005年a)，我们比较了ABC和ABC–MCMC采样器的精确度。我们测量了整个后验分布的均方根积分误差（RMISE），定义为

其中μ_k个是的真实值k个th参数和

\（f（\mathrm{{theta}}{k}{vert}s）

是估计的后部密度。我们还计算了三个点估计值（模式、平均值和中值）的误差，作为点估计值和真实参数值之间的绝对差值。为了评估估计的整体质量，我们测量了所有参数上这些精度度量的几何平均值，与使用100000次模拟进行的传统ABC方法相比。更正式地说，如果

\（x{i{\mathrm{{ABC}}}）

是对我th参数(⁠

\（i｛＝｝1，｛\ldots｝，｛\，｝n\）

⁠)则ABC–MCMC方法的相对精度RA定义为

\（\mathrm{RA}{=}\左（{\sum}_{i{=}1}^{n} x个_{i_{\mathrm{ABC}{-}\mathrm{MCMC}}}/x_{i{\mathr{ABC}}\右）^{1/n}\）

⁠.

不同方法得到的后验分布的覆盖特性也值得检验。通过覆盖率，我们的意思是在给定的可信区间内真实参数值出现的时间比例。例如，80%和95%可信区间应包括概率分别为0.8和0.95的真实参数。换句话说，真实参数值的后分位数应均匀分布在[0，1]中(厨师 等。2006). 为了检查不同方法获得的后验分布的总体覆盖特性，我们因此计算了100个已知真实参数值的测试数据集的参数后验分位数。然后用经典的Kolmogorov–Smirnov检验来评估它们的一致性。

ABC–MCMC在人类进化中的应用：

我们研究了三个非洲人类种群（属于最近分裂的尼日尔-刚果语系的约鲁巴人和曼登卡人种群，以及姆布蒂·俾格米人种群，以下简称俾格米人）之间的遗传关系，分析了约800个微卫星标记(拉马钱德兰 等。2005). 我们假设了一个进化模型，其中曼登卡人和约鲁巴人最近从尼日尔-刚果祖先种群分化而来，而他们自己更早地从俾格米人种群分化而出（参见图1). 假设所有当前和祖先人口对交换移民，但交换速度不同。我们假设两个尼日尔-刚果人口之间以及俾格米人和祖传的尼日尔/刚果人口之间的迁移是对称的，但俾格米斯人和两个尼日-刚果人口的迁移率是不对称的。为了确保与模拟逐步突变模型的良好拟合，我们只考虑331个四微卫星中的一个子集，这些四微卫星的缺失数据小于5%或编码为缺失数据的不完全重复等位基因。

模型参数从以下均匀分布中得出：每个群体的基因拷贝数，N个∼单位[0, 20,000]; 每代迁移基因的数量，牛米∼ 10^单位[0,2]; 尼日尔-刚果两种群的分化时间_天然气∼单位[0,1000]代；和祖先发散时间_A类∼单位[04000]代。我们对发散时间施加了额外的限制，使得TDIV_A类>TDIV公司_天然气导致TDIV之前不均匀_A类所有模拟都是使用SIMCOAL 2.0程序进行的，假设一个纯逐步突变模型具有每代平均突变率

拉入单位[2 × 10⁻⁴, 7 × 10⁻⁴]. 我们实现了以伽玛分布的位点特异性突变率(⁠

\（\mathrm{{\alpha}}，\mathrm{{\alpha}}/\mathrm-{{bar{\mu}}}}\）

⁠)绘制形状参数α单位[8, 16]. 计算与上述定义相同的统计数据集，并转换PLS。我们保留了MCMC链和参数估计的前11个PLS组件。

结果和讨论

各种ABC方法的性能：

我们研究了不同采样器相对于ABC方法在无迁移的种群散度模型中的性能(表1). 正如预期的那样，基于500万次模拟的传统ABC估计（ABC–5M）确实比仅在100000次模拟（ABC–100K）上进行的ABC估计有所改进，并且其推断的后验概率具有很好的覆盖率（Kolmogorov–Smirnov后验分位数一致性检验，第页= 0.469). 正如预期的那样，ABC–MCMC方法的性能在很大程度上取决于用于运行MCMC链的建议范围（φ）和公差（ε）值的选择。总的来说，对于所有φ值，随着公差ε的降低，ABC–MCMC的精确度增加。除了非常小的提案范围（φ≤0.01）外，它比ABC–100K要好。在这种情况下，参数的相对误差可能非常大，并且后验分布没有通过后验分位数均匀性测试，这表明在这些情况下，没有充分探索参数空间。总的来说，ABC–MCMC在φ=1和公差ε=0.01的建议范围内表现最佳，其中精度测量结果甚至优于ABC–5M获得的测量结果。在这种情况下，链条的接受率为～30%，且后部的覆盖范围足够(第页= 0.23). 我们注意到，ABC–MCMC取样器的合格率随着公差水平ε急剧下降。例如，当ε=0.0001时，我们只能对极少数的伪观测数据集运行ABC–MCMC链（φ=0.001时最多六个）。此外，接受率还受到提案范围φ的强烈影响，φ=0.001比φ=5高出两到三倍。如所示图2对于六个随机数据集，ABC–MCMC方法通常会导致后验分布接近于使用ABC–5M获得的后验分布。使用ABC–100K获得的后验结果在大多数情况下更宽，与中报告的精度测量结果一致表1.

然后，我们在人口差异与迁移模型下比较了不同的ABC采样器。结果显示于支持信息,表S1在质量上非常类似于没有迁移的情况。ε=0.01和φ=1的值再次为ABC–MCMC方法带来了最佳的总体精度，它显示了正确的覆盖特性(第页= 0.29). ABC–MCMC的精度再次与ABC–5M方法的精度相当，证实了与传统ABC相比，该方法大大减少了计算时间。

PLS组件的用途：

在图3，我们比较了不同模型和不同统计集的真实参数的后验分位数分布。而基于PLS分量在ABC–5M下获得的后验分位数分布是均匀分布的（Kolmogorov–Smirnov检验，第页=0.469和第页=0.193（对于有迁移和无迁移的散度模型，分别为0.193），在使用原始统计数据时，情况并非如此。通过考虑ABC–5M中的所有统计数据获得的后验分位数确实不是均匀分布的，并且往往太大(第页=0.002和第页< 10⁻¹²对于有迁移和无迁移的散度模型），这意味着当使用所有统计数据时，真实参数被全局低估。这些结果表明，PLS的使用总体上改善了可信区间的覆盖特性。这是意料之中的，因为许多汇总统计数据可能只包含很少的参数信息，这使得很难计算有意义的距离。当然，根据与参数的理论关系或一些评分程序，精心选择了一小组汇总统计数据(例如,乔伊斯和马约拉姆2008)也可能导致无偏的后验分布。然而，目前的PLS方法似乎提供了一种客观的方法来降低汇总统计空间的维数，同时尽可能保留有关参数的信息。

ABC–MCMC在非洲进化中的应用：

我们估计了三个非洲人口之间差异和迁移模型的参数(图1)基于331颗微卫星的数据，使用我们的并行ABC–MCMC方法，1000个独立链，10000个模拟（包括启动），建议范围φ=1，容差水平ε=0.01。参数的后验分布见图4。我们发现了尼日尔-刚果两个族群最近的差异（142代或～3550年前，基于25年的世代时间），这与西非农业扩张的年龄以及尼日尔–刚果语系的多样化非常吻合(木材 等。2005). 俾格米人和尼日尔-刚果人之间的分歧要古老得多T型_模式 =2135代前（～53400年），与之前的研究基本一致(金塔纳-穆尔西 等。2008;贝尔杜 等。2009)但我们注意到，此时95%的最高后验密度可信区间相当大（1075-3712代）。虽然观察到的数据与约鲁巴人和曼登卡人之间的高基因流动率相一致，但俾格米人与尼日尔-刚果两个人口之间交换的移民总体上要少得多。我们发现一些证据表明，从俾格米人到约鲁巴人的基因流水平高于从另一个方向到俾格米族的基因流，这种不对称性已经在先前对俾格米语和邻近人群之间的基因流进行的分析中观察到了(金塔纳-穆尔西 等。2008;贝尔杜 等。2009). 相反，我们发现俾格米人和曼登卡人之间没有迁移不对称，这可能是由于这两个种群之间的地理距离较大。我们的结果进一步表明，俾格米人和祖传的尼日尔-刚果人之间的基因流动水平更低，这表明当时非洲人口的细分程度更高。人口数量的后验分布都很宽，并且指向相对较大的值，即使是俾格米人N个_模式=10876个基因拷贝。有趣的是，我们发现了四个微卫星的平均突变率(⁠

⁠)比之前估计的要低得多(⁠

\（\mathrm{{bar{{mu}}}}{\sim}6.4{times}10^{{-}4}\)

⁠,日沃托夫斯基 等。2003)，其跨位点变异性相对较低（α≅14）。

结论：

我们在这里展示了如何使用无帽MCMC方法来生成近似的后验分布。自

必须选择相对较大的尺寸，以确保充分混合，SS–MCMC链的稳定分布

\（f（\mathbf{\mathrm{{theta}}}{\vert}\mathrm{{delta}}{leq}\mathr{{delta}}{{\mathr m{{varepsilon}}）\）

可能仍然是后验分布的相对粗略的估计

\（f（\mathbf{\mathrm{{\theta}}}{\vert}\mathbf{\mathr m{S}}）\）

⁠，以及旨在获得

\（f（\mathbf{\mathrm{{\theta}}}{\vert}\mathrm{{delta}}{\trightarrow}0）{\approx}f（\mathbf{\tathrm{\theta}}}}{\vert}\mathbf{\mathr{S}}）

需要在SS–MCMC链的输出上执行。链的过渡机制需要微调，以使其适当混合，并导致具有足够覆盖率的后验。公差水平ε=0.01和建议范围φ=1（标准偏差）的组合似乎在两种不同复杂度的模型下提供了总体最佳结果。相似的值φ和ε应确保在其他情况下良好混合，因为建议范围φ和初始公差区间ε的宽度以通用方式表示：φ以参数保留值的标准偏差单位表示，因此，它应该针对不同的参数进行放大，并根据观测数据进行调整，而ε只是与观测值任意接近的模拟的一部分，它既不取决于汇总统计的选择，也不取决于模型的参数化。最后请注意，在回归调整之前，我们没有对ABC–MCMC链的输出进行任何细化，因为拒绝步骤预计会删除大部分自相关，如果链收敛，这不会影响估计。

请注意，我们报告的后验分布通常具有非常相似的模态值，但比在完全似然IMa下获得的值略宽(嘿和尼尔森2007)方法（未显示结果）。虽然我们不建议在存在完全似然法的情况下使用ABC方法，但我们的结果表明，对于没有基于概率的估计可用的复杂场景，可以使用ABC–MCMC进行相对良好的研究，并且其计算成本仅为传统ABC的一小部分。此改进和类似改进(例如，文学硕士。博蒙特、J.-M。科尔尼埃、J.-M。马林和C.P。罗伯特，未发表的结果）应该非常有用，因为ABC技术在示范遗传学研究中的应用越来越多（参见，例如,塔尔蒙 等。2004;Excoffier公司 等。2005年a;汉密尔顿 等。2005;陈 等。2006;希克森 等。2006;施赖纳 等。2006;帕斯夸尔 等。2007;勒格拉 等。2007;罗森布拉姆 等。2007;科尔尼埃 等。2008;Neuenschwander公司 等。2008). 实际上，使用现有的仿真程序，可以在几个小时内为任意复杂的模型建立仿真文件(例如,哈德逊2002;拉瓦尔和Excoffier公司2004;科尔尼埃 等。2008)允许人们专注于现实的进化模型，而不仅仅局限于为其开发特定程序的模型。PLS转换还应该允许人们从大量的汇总统计数据中提取尽可能多的信息，同时保持问题的维度相对较低。虽然ABC框架下的参数估计仍然需要大量的计算时间，但它应该允许进化遗传学家合理地估计他们真正感兴趣的参数，而不是要求他们将兴趣转移到可以获得完全似然解的问题上。

附录

马约拉姆 等. (2003)已经正式证明了参数的后验分布θ给定模型的

可以从MCMC的平稳分布中获得，无似然，其中数据D类'根据当前参数值模拟

\（\mathbf{\mathrm{{theta}}{^\prime}}\）

有可能被接受

\（h（\mathbf{\mathrm{{theta}}}{\rightarrow}\mathbf{\mathm{{theta}}{^\prime}}）{=}\mathrm{min}{\mathrm{{theta}}）/\mathrm{{\pi}}

如果D类'等于观测数据D类并以其他方式拒绝。当数据被一组汇总统计数据替换时秒,马约拉姆 等。(2003)建议接受概率参数

\（h（\mathbf{\mathrm{{theta}}}{\rightarrow}\mathbf{\mathbm{{theta}}{^\prime}}）\）

if摘要统计S′模拟自

\（\mathbf{\mathrm{{theta}}{^\prime}}\）

任意接近观察到的汇总统计数据（如果

\（{\Vert}\mathbf{\mathrm{S}{^\prime}}{-}\mathbf{\mathrm{S}}{\Vert}{=}\mathrm{{\delta}}{leq}\mathr{{\delta}}_{\mathr}{\varepsilon}}}）

⁠). 在这种情况下，他们认为这种SS–MCMC链的平稳分布是

\（f（\mathbf{\mathrm{{theta}}}{\vert}\mathrm{{delta}}{leq}\mathr{{delta}}{{\mathr m{{varepsilon}}）\）

⁠我们在下文中说明了这一点。

SS–MCMC的后向分布：

请注意

根据定义，也是在简单ABC算法下保留参数的分布，如果随机生成的参数导致汇总统计，则可以接受这些参数

\（\mathrm{{\delta}}{\leq}\mathrm{{\delta}}_{\mathrm2{{\varepsilon}}}\）

并以其他方式拒绝。这意味着SS–MCMC方法与具有相似容差水平ε的ABC算法具有相同的平稳分布。在图S1，我们通过报告距离的分布，实证证明了这一点

\（\mathrm{Pr}（\mathr m{{delta}}{<}\mathrm{{delta}}_{\mathrm2{{varepsilon}}}）

根据传统ABC拒绝算法以及SS–MCMC生成。

脚注

脚注

致谢

我们感谢Samuel Neuenschwander、Nicolas Ray、Olivier François和Gerald Heckel的有益讨论。我们还要感谢马蒂厄·福尔（Matthieu Foll）、大卫·巴尔丁（David Balding）、朱迪·嘿（Jody Hey）和一位匿名评论员对手稿的评论。我们还感谢Nelson Fagundes建议使用GW*统计数据。这项工作得到了瑞士国家基金会（编号3100A0-11072）向L.E。

工具书类