蒙特卡罗方法在多粒子群分析中的应用

摘要

我们提出了一种新的非参数贝叶斯方法,用于在观测单元具有来自多个来源的数据的情况下进行聚类分析。我们的方法使用粒子Gibbs采样器进行推理,其中在Gibbs采样器中使用条件粒子滤波器联合更新簇分配,从而改善MCMC链的混合。我们开发了几种提高算法计算性能的方法。这些方法可以在不牺牲精度的前提下,实现一个数量级的性能改进,并且可以更广泛地应用于具有单个数据集的混合模型的贝叶斯推理。我们将我们的算法应用于243例肾肾透明细胞癌患者的风险队列的发现,使用癌症基因组图谱中的样本,这些样本有基因表达、拷贝数变异、DNA甲基化、蛋白质表达和microRNA数据。我们确定了4种不同的一致性亚型,并显示它们是生存率的预后因素(\(p<0.0001\)).

介绍

聚类分析可以广泛地描述为在数据集中推断底层组结构的任务。群的定义使得一个群中的观测值与其他群集中的观测值更相似。聚类分析已经被证明是一种流行的探索性工具,并在市场细分、机器学习、数据压缩和基因组分析等领域找到了应用领域。

在分析基因组数据时,我们可以根据基因组成推断出患特定疾病的患者的风险队列,或者我们可以通过推断基因组来帮助了解其功能。然而,这些应用领域提出了在其他情况下通常不会遇到的问题,特别是由于每个观察单位(即前一示例中的患者)可能具有来自多个数据源的数据,例如基因表达、DNA甲基化或拷贝数变异。这些数据源提供了对底层过程的补充但不同的观点,因此,对这些数据的分析能够在一个单一的综合分析中包含这些数据源是至关重要的。

综合聚类算法为一组观测单元推断出一个聚类结构,这些观测单元的数据来自多个来源。虽然我们希望将这些观测单位分配给集群,以解释每个数据源中的不同变化,但不同数据集的集群结构是相同的并不重要。使这项任务复杂化的是,数据源不必是共享类型的,例如,我们可能需要将连续数据与离散数据集成,因此我们不能简单地将各种数据源连接到单个数据矩阵中。有许多综合聚类分析的方法:在潜在的潜在变量水平上进行聚类分析(参见,例如,Shen等人。2009;Gabasova等人。2017;McParland等人。2014,2017或对多个数据集的独立聚类解决方案求平均值的集成方法(参见Monti等人。2003;洛克和邓森2013, ). 在这些和其他类似的方法中,术语“相关聚类”、“一致聚类”或“多视图聚类”可以代替“综合聚类”。在本文中,我们提出了一种新的基于多数据集集成(MDI)框架的集成聚类算法(Kirk等人。2012)一种灵活的基于模型的综合聚类算法,有利于数据集之间的信息共享。

本文的其余部分结构如下:第一节2介绍MDI;节介绍ParticleMDI,我们在MDI上的开发,它使用条件粒子过滤器更新集群分配;Sects3.13.2提出了几种提高ParticleMDI计算性能的技术,这些技术在单数据集上下文中更为普遍适用4概述了该方法在合成和真实数据集中的应用实例,展示了推断有临床意义的亚组的能力;第三节5总结道。

多数据集集成

MDI(Kirk等人。2012)是一个对潜在不同数据类型的多个数据集进行综合聚类分析的框架。它将标准混合模型推广到多个数据集的上下文中,允许成对数据集的簇结构相互通知。模型推断出数据集之间的信息共享程度,允许根据数据集的结构相似性为不同的数据集对共享不同级别的信息。这些跨数据集比较是在推断的集群分配级别上进行的,这为MDI提供了自然地建模不同数据类型的灵活性。萨维奇等人(2013)演示MDI如何在基因组数据集中发现有临床意义的群体。

我们首先描述一个综合混合模型的生成模型,在这个模型中,没有信息在Eq中的数据集之间共享1(柯克等人。2012). 我们假设K数据集也被观测到n观测单位。数据集k\(púk\)其中收集的特征\((n倍pаk)-维度数据矩阵\(X)具有\(i,j)-第进入\(x{i,j,k}\). 观测值的顺序在数据集中是一致的,例如属于\(X)\(X)对应于同一观察单位,如果我们的目的是推断人群,这些数量是指对同一个人进行的不同测量。在本文中,我们将使用符号\(x{a:b}\)表示序列\(x{a},x{a+1},圆点,x{b-1},x{b}\). 模型是

$$$$$\begin{Align}\ begin{Align}\ begin{Align}x{i,1:p UK,k}| c{i,k},theta&\sim F(\theta{c{c{i i,k},k})\\\c\c{i i,k}|\pi k&\sim多项多项式(\pi{1,k},k},点,\点,pi{N p{1,k},点,pi{N,k{N,k})\\\\pi \\\pi{k{k}=(\pi{1,k}、\dots、\pi{N,k})和\sim Dirichlet(\alpha/N、\dots、\alpha/N)\\\theta{c{i,k},k}&\sim G^0}结束{aligned}\end{aligned}$$
(一)

一个潜在的集群分配变量\(c{i,k}\在{1,2,dots,N}\中)叙述一个观察\(x{i,k}\)给其中一个N混合物模型中的组件。每个组件都有一个特定的概率模型,F,具有特定于群集的参数\(\theta{c{i,k},k}\)具有共轭先验\(G^0\U k\). 集群标签,\(c{i,k}\),具有关联权重\(\pi{1:N,K}\),柯克等人(2012)给Dirichlet过程指定一个有限逼近的先验(Ishwaran和Zarepour2002)它定义了一个Dirichlet多项式分配混合模型(Green和Richardson2001). 该算法不需要预先指定簇的数目,而是要求最大可能的簇数,N,以适应数据。这通常被称为过度拟合混合模型,因为模型指定的簇数可能大于实际数。卢梭和门格森(2011)证明了这种模型的后验分布往往集中在稀疏表示上,其中任何多余的簇不被任何数据点占据。N在所有数据集中共享,则每个数据集中推断的簇数(占用的组件数)不必相同。虽然在实践中N可能被指定为大量的观察,选择往往会由计算时间和可用资源决定;柯克等人(2012)建议用观察数的一半作为折衷。

注意等式1独立地为每个数据集的簇结构建模。为了促进数据集之间的信息共享,Kirk等人(2012)介绍\({K\atopwithdelims()2}\)-尺寸协调参数,\(\varPhi\),反映了成对数据集之间簇结构的依赖程度,如图所示1. 联合分配\(c{i,1},c{i,2},\点,c{i,K}\)在等式中1替换为以下选项

$$$$$\开始{对齐}p(c{i,1},1},c{i,2},点、点、c{i,K}{\varPhi\varPhi)\Proppto{\displaystyle\ prod}{K=1}{K=1}{K}}pi{c{i,K},K}{显示样式\ prod{K=1}{{1}}{{K=1}{K-1}{{1}{{1}}{{1}}{\显示样式{{{{{l=K+1}^K}(1+\phi{kl}\mathbb{1}(c{i,K}=c}{i,l}))\end{aligned}$$
(二)

的每个元素\(\varPhi\)反映了在聚类标签级别的数据集的聚类结构中的成对一致性,其中\(\phi{k,l}\)表示数据集中簇结构之间关系的强度k. 为了捕捉数据集之间的依赖关系\(\varPhi\)值用于增重在两个或多个数据集中属于同一个簇的观测值的可能性。例如,如果\(\phi{1,2}\)如果数据量很大,MDI将更倾向于将数据集1和2中的观测值放在同一个群集中。不同的\(\varPhi\)值允许数据集之间存在不同级别的部分一致性,而不需要严格的一致性,例如\(\phi{k,l}=0\),然后是数据集中的群集分配k不会互相影响。

图1
图1

MDI和ParticleMDI在三数据集情况下的图形模型表示\(K=3).\(x{i,k}\)表示观察在数据集中k集群产生\(c{i,k}\)带参数\{\thetuk,塔卡},它们被赋予了优先权\(克^{(0)}\). 数据集中的簇k有优先分配权,\(\pi\u k\)它们自己被赋予\(Dirichlet(\alpha/N\)优先。这个\(\varPhi{i,j}\)值允许在数据集中进行分配j互相告知。[图来自(Kirk等人。2012)]

算法1展示了Kirk等人的方法(2012)在MDI模型中,采用Gibbs抽样方法在更新聚类分配和更新超参数之间交替进行。集群分配的更新方式与标准中的类似k-均值法:在初始的观测值分配到簇之后,观测值在簇之间迭代移动,条件是所有其他的簇分配保持不变。具体地说,给定当前的观测值分配到簇,MDI迭代每个观测值,计算将其移动到每个簇的边际可能性,并基于这些边际可能性概率为其分配一个聚类标签。这种一次更新聚类标签的方法可能会抑制Gibbs采样链的混合,这意味着算法的后续样本可以高度相关。因此,该算法在探索样本空间时可能比较慢,一旦发现局部最优,就很难摆脱局部最优,因为如果与当前的建议不相似,就很难提出更好的观测值分配给簇。这种高度自相关的MCMC链在短期内可能不代表后验分布,因为它们过度代表了链中的局部勘探区域。从长远来看,这些问题可以通过产生更多的样本和随后的稀释来缓解。然而,这需要产生额外的计算成本来生成这些样本,并且,根据混合被抑制的程度,这些额外样本的生成成本可能过高。

图A

分词

在本文中,我们提出了ParticleMDI,一种在MDI框架内构建的算法,它保持了MDI在建模灵活性方面的优势,但将推断方法从Gibbs采样改为粒子Gibbs采样(Andrieu et al。2010)通过使用条件粒子过滤器更新簇标签,该过滤器允许我们联合更新群集分配组。

粒子滤波器或序贯蒙特卡罗(SMC)通常用于状态空间模型中的推理,我们希望了解潜在变量,\(x\n\),给出一个观察结果,\(是的),当时. 这些方法的目的是按顺序近似密度序列\(p(x{1:n}| y{1:n})\)对于\(第1页)以及边缘可能性序列

$$\begin{aligned}p(y{1:n})=\int\dots\int p(y_n\vert x_n)p(x_n\vert x{1:(n-1)})\,dx_1\dots dx_n\end{aligned}$$

对于\(第1页). 他们通过将问题分解为从一系列中间密度中取样来实现这一点\(p(x{1:i}| y{1:i})\)对于\(i=1,\点,n\). 在每个,该方法首先生成建议分布中的“粒子”,使用\(i-1\). 粒子是随机样本,代表潜在变量的不同潜在实现,\(x{1:n}\). 每个粒子根据重要性权重重新加权\({\xi^{(m)}}\)对于\(m=1,\点,m\)这取决于\(x)以及观测值\(你我\). 括号中上标索引的符号是指所讨论的特定粒子。在重新加权步骤之后,某些粒子的权重可以忽略不计,因此重采样步骤通过移除权重较低的粒子并替换为权重相对较高的粒子,将粒子的加权样本替换为未加权集。存在各种重采样技术(见Hol等人。2006)其中,我们实施了系统的重采样,因为它比其他方法提供了显著改进的混合,并降低了路径退化的水平(肖邦和辛格)2015;格里芬2014). 这样就可以通过集中计算资源来利用这些更重粒子获得的信息。通常不建议在粒子过滤器的每一步进行重采样,因此通常只有当由粒子权重的方差测量的有效样本大小(ESS)降到某个阈值以下时才执行重采样(Liu和Chen1995). 这通常被选为粒子数的一半(Doucet和Johansen2009). 在混合模型的推理中,格里芬(2014)Bouchard-Côe等人(2017)发现自适应重采样在每一步都比重采样带来更高的性能。

然而,聚类分析的任务并不典型地被视为是随着时间的推移而发展的,因此似乎不适合于顺序方法。然而,粒子滤波方法已经成功地应用于这项任务(见,例如,肖邦)2002;费恩黑德2004;格里芬2014Bouchard-Côe等人。2017, ). 该方法包括处理观察指数,1:n,作为人工时间索引。在标准混合模型中,一个潜在变量,\(c{i}\),反映了分配给观测的聚类标记\(x{i}\). 然后,可以将问题表述为\(p(c{1:n}| x{1:n})\)对于\(第1页),通过中间密度\(p(c{1:i}| x{1:i})\)对于\(i=1,\点,n\). 费恩黑德(2004)演示与指定观测值相关联的增量粒子权重\(x)聚集可以计算为

$$$$\开始{对齐}\展示样式{\西安{对齐}\展示样式{\xi{(m)}=\xi{(m)的}\倍p(c{i}{(m)的}p{c{1(i-1)}{(m)的1:(i-1)}{(m)的x{(m)i})=\xi{(m)的}金金金金{a{1{1{1{1{N{N{f f f(x U i i | x{1:1:1:1:1)f(1:i(i i(i i i-1)},c{1:(i-1)}^{(m)}=a)}、\N编号\\\结束{aligned}$$
(三)

哪里f(.)是观察的后验预测密度\(x)考虑到\(x{1:(i-1)}\)为了什么\(c{1:(i-1)}^{(m)}=a\). 这种方法有助于对不同类型的数据进行建模,前提是我们可以分析地得出后验预测分布。

混合模型的SMC方法可以扩展到集成MDI模型,将每个粒子视为所有粒子的簇分配的实现K数据集。也就是说,我们想推断\(p(c{1:n,1:K}| X_1,\dots,X_K)\)将观测值分配给群集时,只根据该数据集中的分配情况。我们更新情商综合上下文如下

$$$$$$$$$$$$$$$$$$$$${对齐}\Si{(m)}\times\displaystyle\ prod{k=1}{{k}{k}p p(c〈c{i i,k,k}{(m)}}(m)},x{1:i,i,1:p UK,k}{1:i,1:p UK,k}{1{k{1{1{1{1{1{1{1{1{1}{k{1{1{1{{l=k+1}^{k}(1+\phi{k,l}\mathbb{1}(c{i,k}^{(m)}=c{i,l}^{(m)}))\n编号\\\end{aligned}$$
(四)

我们推断\(\varPhi\)通过粒子Gibbs采样器(Andrieu等人。2010),一个粒子MCMC扩展Gibbs取样器,其中一个粒子\(c{1:n,1:K}^{(1)}\)用于更新超参数,随后用作粒子过滤器的下一个过程的输入。超参数的条件分布在Kirk et al(2012),尽管我们为\(\varPhi\)在补充材料中。粒子是从粒子滤波器的初始通道采样的,并且在随后的条件粒子滤波器扫描过程中,确保其在重采样步骤中的生存。由于这个粒子是用整个观测序列推断出来的,所以它有助于引导其余粒子朝向好的采样空间的区域。我们在算法2中概述了这种条件粒子滤波器。

该算法的规范是灵活的,因为它允许对任何统计分布进行建模,我们可以对其进行分析,得出后验预测分布。在本文中,我们将连续值数据建模为高斯分布,具有正态伽马先验(见墨菲)2007)我们将离散值数据建模为具有Dirichlet先验的范畴分布(参见Bernardo和Smith)2001). 其他数据类型也可以使用共轭对,如泊松分布和伽马分布。在计算速度上,我们将特征视为独立的,但也可以合并从属特征。将特征视为独立的也允许我们结合(Savage et al。2013). 在这种情况下,特定特征被选择的概率与Bayes因子有关,Bayes因子是指在推断的聚类模型下,该特征中数据的边际似然与所有观测值都属于一个共享聚类的空模型的边际似然之比。

然而,MDI的粒子Gibbs方法计算量更大。虽然一个算法的计算时间和它的混合时间(接近后验分布所需的迭代次数)之间可以进行权衡,但算法并没有在两者之间取得适当的平衡。在本节的剩余部分中,我们将探讨一些提高算法计算效率的技术。

图B

块更新吉布斯采样器

考虑算法2中描述的粒子Gibbs采样器中的传播步骤。每次观察,\(x{i,1:p},k}\),一个群集标签,,与后验预测分布成比例抽样\(f(x{i,1:p_k,k}| c{1:i,k},x{1:(i-1),1:p_k,k})\). 那就是观察在数据集中k分配了一个簇标签,,基于第一个\(一-1\)观察。什么时候?\(我同意)这将对集群标签进行有意义的估计\(c{i,1},圆点,c{i,K}\). 但是,例如,\(i=3\)估计就没那么有意义了。因此,结果中存在一个依赖于数据被合并到粒子过滤器中的顺序。减少这种依赖性的一种方法(如算法2和Griffin中所做的那样)(2014))正在更新\(c{1:(i-1),1},圆点,c{1:(i-1),K}\)在重采样步骤中。虽然这意味着所有的聚类标签都可以使用完整的数据集进行更新,但这大大增加了算法的计算复杂度。在重采样不自适应的方案中,或者在自适应方案中的最坏情况下,这将导致每次迭代所需的操作,,从\({\mathcal{O}}(n)\)\({\mathcal{O}}(n^2)\),假设每个单独步骤都可以在恒定时间内执行。这种解决方案还使浪费计算的计算问题复杂化,即根据很少的先前观测值将观测值分配给簇。

此外,还有一个特定于整合语境的问题。如公式所示4我们根据共享同一个聚类标签的观察结果,为数据集间的一致性膨胀粒子权重。虽然我们在粒子过滤器的过程之间进行标签匹配,但是在一个过程中不能保证簇标签引用了数据集中相同的观测分区。

我们通过增加粒子吉布斯采样器来解决这两个问题,建议在粒子过滤器的每个过程中只更新一部分簇分配。剩余的分配从上一次传递中保持不变;我们将此部分称为\(\rho\),其中\(0<\rho<1\). 我们随机选择指数1的子集:n大小\(\lfloor n\rho\rfloor\)并将这些观测值与标准粒子Gibbs采样器中的参考粒子相同。虽然算法基于相继的蒙特卡罗采样器,数据的顺序,事实上,并不重要。因此,随机改变观察顺序并假设第一个是无关紧要的\(\lfloor n\rho\rfloor\)从以前的粒子群中通过粒子过滤器的分配。这种排列确保了我们不会在每次扫描时对相同的观察部分进行条件调整,并且将固定部分放在开始位置,从而避免了定义在粒子的不同固定部分之间建立桥梁的潜在挑战性任务。我们称之为置换观察指标的函数\(\sigma(\cdot)\),确保对每个数据集执行相同的置换。然后,算法继续进行,未来的分配将根据先前观察到的数据进行更新。也就是说,粒子吉布斯采样器通常交替从\(p(\theta | x{1:n},c{1:n})\)\(p{\theta}(c{1:n}| x{1:n})\)(其中\(\theta\)表示粒子过滤器每次运行时固定的超参数),ParticleMDI从\(p(\theta | x{1:n},c{1:n})\)\{\r{r\r\n西格玛}(RhoIl}\r\n西格玛)(图2).

图2
图2

标准粒子吉布斯扫描()在三粒子、单数据集的情况下。单个参考粒子(\(m=1\))以(如阴影圆圈所示)为条件,保证它在重采样步骤中生存。圆圈中显示的值为\(c{i,1}^{(m)}\). 在我们的增援中(b)簇分配的一个子集是以所有粒子共享一大块观测值为条件的(\(\rho\)这里,3/n)与参考粒子相同。为了使相同的观测值不在每次迭代中保持不变,观测值的顺序根据排列函数进行洗牌\(\sigma(\cdot)\)

在最初的粒子吉布斯论文(Andrieu et al。2010). 考虑到我们可以对数据进行简单的重新排序,这是粒子Gibbs采样器的这种扩充的一个特别合适的应用;粒子Gibbs分裂合并(PGSM)算法(Bouchard-Côet al。2017),其中粒子Gibbs扫描仅限于共享一个簇标签的观测值和两个随机选择的“锚”点。然后,如果两个锚点属于同一个簇,则PGSM算法将继续“拆分”该簇,否则将“合并”它们。

提供\(\rho\)如果选择适当,这种增加可以最小化对数据观察顺序的依赖性。因此,在第2步中,我们提出了贪婪算法,其中在第3步中,我们提出了一个更新的算法。假设每一步都可以在固定的时间内执行,该算法的计算复杂度现在是\({\mathcal{O}}(\lceil n\times(1-\rho\rceil)\),虽然只更新了一部分聚类标签,但与算法2的比较并不完全准确。

图

找到值\(\rho\)然而,“适当选择”可能并不简单。选择\(\rho\)在计算时间和混合时间之间进行权衡,使用较大的值来减少计算时间,但代价是更新较少的集群分配。然而,计算时间和混合时间之间的权衡不是线性的,因为非常低的值可能导致对观测值的合并顺序的依赖。该算法的贪婪本质是基于数据块作为一个整体的充分摘要的条件。首先,在膨胀粒子权重的情况下,这一点很重要,因为粒子将观测值分配给同一个簇标签。重要的是要确保\(\rho\)足够大,以便在星团标签上嵌入一些意义,确保在可行的范围内,星团标签指的是观测单元的相同分区。这应该是引用粒子的情况,因为我们在这个粒子中跨数据集对齐簇标签。因此,\(\rho\)应详细说明,以便在余下的观察中保持这一含义。例如,如果观测值的条件部分很小,足以完全排除许多簇,那么就不能保证分配给这些簇的标签在遇到时会在数据集中一致。同样重要的是\(\rho\)指定一个足够大的值,以克服对数据合并到条件粒子过滤器中的顺序的依赖性。例如,如果我们仅以一次观察为条件,这些数据将不具备足够的权重来指导随后的观察。然而,通过重采样使参考粒子的扩散有助于缓解这两个问题,但可能无法完全消除它们。我们探讨了选择\(\rho\)在门派4.1.

利用粒子滤波中的冗余

粒子滤波方法的精度要求对,粒子数。在混合模型的背景下,格里芬(2014)发现即使使用相对较少的粒子也可以获得良好的性能。在ParticleMDI中,对于每个粒子,我们需要存储指定的簇标签以及汇总统计信息,以便于粒子的传播。这就需要大量的内存访问,这会阻碍算法的计算性能。

然而,并不是所有这些存储都是必需的,因为粒子过滤器的性质会导致粒子中的冗余。重采样步骤会生成某些粒子的副本,同时移除其他粒子。在重新采样之后,粒子将随着附加观测值的到达而传播。传播需要评估每个簇的新观测的后验预测密度,如等式所示4. 这些计算对于粒子的每个副本都是相同的,并且会显著增加算法的计算负担。更具体地说,考虑到我们处理的是一个离散的状态空间,几乎没有什么不同的可能后代对于每个粒子,结果是,即使不重新采样,我们也可能生成许多相同的粒子。

这些冗余已经被费恩黑德讨论过了(2004)他提出了一种对重采样算法的扩充,以尽量减少重复粒子的数量。其基本原理是,例如,一组粒子指数\(\左[1,1,2,2\右]\)各有重量\(\左[0.1,0.2,0.3,0.4\右]\)相当于只有两个具有权重的粒子\(\左[0.3,0.7\右]\)因为可以考虑所有可能的粒子后代。费恩黑德(2004)建议保留标准化权重高于阈值的所有粒子,对剩余粒子执行重新采样,将粒子重复的数量最小化。

我们不是尽量减少冗余,而是试图利用它。为了避免执行冗余计算,我们需要能够识别哪些粒子与其他粒子相同。我们可以从重新采样的粒子索引中识别出那些由于重采样步骤而相同的粒子。对于其他算法,我们建议为每个粒子分配一个运行ID,如算法4所述,并且只为每个唯一的粒子计算一次增量重要性权重并存储在哈希表中。但是,此方法不会捕获所有重复粒子。考虑两个粒子,它们在簇标签排列之前是相同的。这些粒子不会共享粒子ID,因为ID是基于指定的标签构造的。这可以从图中看出,其中粒子1和粒子2可以通过切换分配标签变得相同。同时考虑两个粒子基本相同的情况,例如在10个推断出的簇中,这两个粒子为其中的8个指定了相同的标签,而其余两个簇则具有不同的标签。同样,这些粒子不会被识别为重复的,尽管增量重要性权重的评估在两个粒子中是相同的。为了克服这一点,我们将算法的一般设置从包含不同环境的粒子转移到存在一个全局环境的场景中,每个粒子都索引到其中,如图所示. 因此,当计算增量重要性权重时,我们现在只需要评估属于数据中推断出的每个独特簇的观测值的后验预测密度。

有花纹的

为了检验使用这些方法可以获得的好处,我们研究了随着粒子数的增长计算时间的增长。由于避免冗余计算的好处会随着所分析的特定数据集的复杂性而变化,因此我们考虑两种情况:一个简单的高斯数据集包含两个定义良好的具有均值的聚类\(+5\)\(-5\)分别在100个观察值中,以及Fisher的虹膜数据集(Fisher1936)它包含三个星团,其中两个在150个观测中没有很好地分开。第一个数据集代表了一个最佳情况,说明了我们实现ParticleMDI的方式可能带来的计算时间上的改进,而Iris数据集可能代表在实践中可以预期的更温和的改进。4显示了这些数据集中粒子数对计算时间的经验影响。所有的分析都是针对1000个MCMC样本\(\rho=0.25\),并且给出的结果是100次这样的运行的平均值。在传统的方法中,计算时间会随着粒子数的增加而线性增长。然而,在ParticleMDI下,在两种情况下都明显存在次线性增长,其中\(M=1{,}024\)粒子大约\(15倍)只要完成\(M=2\)对于合成数据,以及\(161倍)和虹膜数据集一样长。这代表大约\(3\%\)\(31.5\%\)分别表示在标准实现下预期的时间,表示计算时间的潜在改进超过一个数量级。即使是在时间尺度的改进中,也要注意到粒子群的识别问题\(M=2\)在计算时间上有所改进。必须指出的是,这些改进在推断的集群分配的准确性方面是没有代价的,因为我们只是避免了对冗余计算的评估。

图3
图3

同一粒子滤波系统的两种表示。显示了一个标准实现,其中每个粒子都包含一个孤立的簇环境。突变步骤包括评估后验预测密度,以便将来观察到每个粒子中的每个簇。b表示ParticleMDI中的实现,其中每个粒子索引到群集的全局环境中。在这种情况下,后验预测密度只需要对全局聚类环境中的每个唯一聚类进行计算

图4
图4

与粒子数增长相关的计算时间的增长。是ParticleMDI在包含两个不同簇的合成高斯数据集上的应用。b是ParticleMDI在Fisher虹膜数据中的应用

提取共识簇

ParticleMDI的输出是每个数据集中每个观测值的一系列簇标签。由于这些簇标签在迭代过程中是可交换的,我们不能简单地计算属于某个特定簇的观测值的后验概率。相反,我们计算后验相似矩阵(PSM)(Monti等人。2003)对于每个K数据集。PSM是一个\(n \乘以n\)谁的矩阵(j)-第个条目的值\(\左[0,1\右]\)表明观察频率j分配给同一个群集。然后,数据集的总体一致性可作为每个PSM的元素平均值(Savage等人,2008年)。2013). PSM可被视为热图,如图所示67. 通过应用分层聚类方法,可以从PSM中推断出最终的一致性群集分配\(1-平方米)作为距离矩阵(Medvedovic et al。2004;Fritsch等人。2009拉斯穆森等人。2009). 从这一分析得到的树状图将产生一个观察的顺序,这样最经常分配给同一个星团的观测对可以在热图中相邻地定位,从而可以启发式地识别潜在的星团结构。

应用

在本节中,我们将展示ParticleMDI在合成和真实数据集中推断簇结构的能力,并探讨第节中描述的粒子Gibbs采样器中块更新对聚类精度的影响3.1. 补充材料中包括更多的例子。

选择的影响\(\rho\)

探索如何选择\(\rho\)影响混合时间和计算时间之间的平衡,我们在单个数据集上下文中运行ParticleMDI,计算预算为1秒。虽然我们承认这种限制在实践中不太可能,但只要在综合数据环境中了解这些相互竞争的目标之间的平衡就足够了。我们生成的合成数据集如下:

  1. 1

    生成簇权重

    $$\begin{aligned}\pi\sim Dirichlet(0.2,0.2,0.2,0.2)\end{aligned}$$
  2. 2

    生成群集分配

    $$\begin{aligned}c{1:100,1}\sim多项式(1,\pi)\end{aligned}$$
  3. 三。

    生成数据集

    $$$\begin{{对齐}x{i,1:32,1}\sim{左{左{开始{阵列{{将{数学{数学{N}}(0,1)&{{}文{为}为}为}c}{i,2}i,2}=1}1}{N}}(2,1)&{{1}文{为}为{c{i i,2,2}2}{数学{N}{{N}}}{}}{}{}(4,1)&{}\text{for}c{i,2}=3\\{\mathcal{N}(-4,1)&{\text{for}c{i,2}=4\\{\mathcal{N}(-2,1)&{\text{for}c{i,2}=5。\右结束{array}\右。}\\x{i,17:32,1}\sim{\mathcal{N}}(0,1)\end{aligned}$$

我们用\(n=100\)我们从ParticleMDI生成样本,计算预算为1s. 我们将[0,1]秒的范围分成250个相等的箱子,对于在每个时间仓结束时生成的最后一个样本,我们计算调整后的Rand索引Rand(1971)作为当前群集分配准确性的度量。5显示这些时间点的中值ARI值范围\(\rho\). 可以看出\(\rho\)-0.75和0.95——表现明显不佳。在这些情况下,ParticleMDI在每次迭代时只更新非常小的一部分数据,因此,从初始分配中转移需要很长时间。相反,选择值\(\rho\)这个值太低-0.01也与性能差有关,因为生成的样本相对较少,收敛速度较慢。价值观\(\rho\)在0.25和0.35之间似乎提供了计算时间和混合时间之间的适当平衡;我们使用的值为\(\rho=0.25\)在剩下的例子中。

图5
图5

变量选择对ParticleMDI输出精度的影响\(\rho\)对于固定的计算时间1s

影响\(\varPhi\)

为了证明参数的影响\(\varPhi\)在允许ParticleMDI跨数据集共享信息时,我们给出了ParticleMDI在Fisher的Iris数据集(Fisher)上的示例应用1936),这是clsuter分析中的一个典型示例应用程序。该数据集包含对三种不同鸢尾属植物的150次观察:刚毛鸢尾属、维吉尼亚鸢尾属和花色鸢尾属,每一种都代表了数据集的相等部分。将ParticleMDI应用于图中这些数据的结果6一个突出的特点是,虽然一个物种的刚毛是明确的可辨认的,其他两个显示出很大程度的重叠。为了解决这个问题,我们在一个综合的环境中进行分析,通过用一个包含真正的聚类标签的额外的分类数据集来补充虹膜数据。虽然我们承认不太可能有一个数据集具有如此完美的地面真相簇结构信号,图6b说明ParticleMDI如何使用来自一个数据集的信息来通知另一个数据集的集群结构。不仅真实的集群结构在整体一致性PSM中很明显,而且在原始数据中推断出的簇结构现在显示了三个物种之间的分离。

图6
图6

ParticleMDI应用于iris数据集的PSM热图()在综合上下文中,虹膜数据由包含真实聚类标签的分类数据补充(b)

癌症患者的亚型发现

MDI的主要动机是综合分析多个生物学数据集。在这一节中,我们将探索ParticleMDI应用于推断一组特定癌症的患者群,并记录每个患者的基因组学和其他“组学”测量值。我们使用这些推断组来回顾性地预测患者的风险状况。

数据

我们考虑了Yuan等人之前分析的癌症基因组图谱泛癌数据集(2014)他发现,当独立考虑时,临床协变量在预测个体生存概率方面的表现大多优于“组学”特征,而综合方法则赋予了额外的预后能力。具体地说,我们关注肾肾透明细胞癌(KIRC)的患者,这是一个由243名患者组成的小组,他们的每种“omic”数据类型以及生存时间数据都是可用的。我们使用了袁等人准备的数据(2014). 在补充资料中,我们提供了ParticleMDI对卵巢浆液性囊腺癌患者数据集的附加应用。虽然ParticleMDI能够同时执行特征选择和聚类分析,但是为了计算目的,我们通过初始减少特征的数量来补充这一点。最初的减少是基于单个集群的可变性,因为这与“集群性”有关(Steinley和Brusco2008)已经在其他类似的分析中进行过(Shen等人。2009;劳勒等人。2016). 对所选特征的事后分析和推断的聚类分布可以潜在地揭示驱动患者预后的特定基因。可用数据如下:

  • 信使RNA(mRNA)-数据平台为Illumina HiSeq 2000 RNA测序V2。在现有的20,203个mRNA特征中,我们选择了其中100个最变异的。我们将这些数据视为连续数据,并将其建模为高斯分布。

  • microRNA(miRNA)-数据平台为Illumina Genome Analyzer/HiSeq 2000。在795个miRNA特征中,我们选择了100个最变异的。我们将这些数据视为连续数据,并将其建模为高斯分布。

  • 反相蛋白阵列(RPPA)-数据平台为MD Anderson反相蛋白质阵列核心平台。在166个特性中,我们保留了100个最具变体的特性。我们将这些数据视为连续数据,并将其建模为高斯分布。

  • DNA甲基化(甲基化)-数据平台是IlluminaInfinium人类DNA甲基化450K。在16484个DNA甲基化特征中,我们选择了100个最变异的。我们在0.85的阈值下对这些特征进行离散化,并去除了“命中率”低于10的所有特征。其余的特征被建模为来自分类分布。

  • 体细胞拷贝数改变(SCNA)已知SCNA在癌症中极为常见,而获得SCNA是癌症的一个已知驱动因素。Trix 0.0是全人类基因组数据平台。共有69个SCNA特征,我们将其离散化,并给出了值\(>0.1\)表示放大倍数和值\(<-0.1\)表示删除。这些值被建模为分类数据。

  • 临床数据临床和管理数据可用于患者,包括年龄、性别、肿瘤分级和存活时间。为了这个分析的目的,我们考虑生存数据和肿瘤分级,作为评估推断的聚类有多有意义的手段。

结果

图中所示的PSMs7显示在五个数据集中推断的群集结构。总体一致的PSM表明存在一个大的簇(左下角),它似乎强烈地反映了SCNA的簇结构和甲基化数据。利用蛋白质数据推断出的结构可以进一步细分这个簇。通过对一致性后验相似性矩阵应用层次聚类分析得到的树状图被切割成四个聚类。后均值\(\varPhi\)数值如图所示8f、 结果表明,在聚类结构上,甲基化和SCNA数据的一致性最强。

尽管在mRNA数据中推断的簇中存在明显的粒度,但是可以在图中看到8a该聚类3通常对应于一组在所选基因中具有特别高水平表达的患者,相对于其他簇,类似地,簇2包含典型低表达的患者。簇3在甲基化数据中也是不同的,图8d、 显示出明显高水平的甲基化8可以看出,属于第3组和第1组的患者比第2组和第4组的患者表现出更大程度的拷贝数改变。肿瘤分级评分如图所示8表明第1组和第3组中的个体通常具有较高的肿瘤分级分数,这表明相对于其他组中的个体,这些个体的预后较差。这是由卡普兰-迈耶生存曲线所证实的9,显示在推断组的生存率有显著差异(对数秩p-价值=\(1.197133乘以10^{-7}\)). Tarone软件测试(Tarone和Ware1977)生存曲线的差异对早期生存曲线的差异也很敏感(\(p<0.0001\)). 我们将这两个结果作为生存分析的显著性检验,通常会受到生存曲线中出现的交叉的限制,正如Li等人所讨论的那样(2015). 值得注意的是,在多个数据集中被认为具有独特特征的第3组,其生存预后明显比其他组差。属于第3组的24名患者的中位生存时间仅为866天,而整个患者队列的中位生存时间为1979天。第1组患者(\(n=79\))以体细胞拷贝数变化高为特征,中位生存时间为1625天。第2组患者(\(n=47\))以mRNA表达水平较低而著称的患者预后较好,中位生存时间为2830天,接近研究期结束。补充材料中进一步分析了推测的团簇结构的生物学意义。

图7
图7

将ParticleMDI应用于KIRC数据集得到的PSMs的热图表示。整体PSM是各数据集特定PSM的元素平均值。集群分配由白色虚线表示

图8
图8

e每个KIRC数据集的热图,以及从ParticleMDI推断出的集群分配,用白色虚线和右侧边缘的灰色条线突出显示。最右边的条形图反映了肿瘤分级分数,较亮的颜色表示更高(更差)的分级分数。后验特征选择概率由下边缘的条表示。f推断的\(\varPhi\)数据元对中的di

图9
图9

癌症数据集的生存曲线分析4.3,显示了每一个推断的集群的不同风险状况。生存曲线差异的对数秩检验给出\(p<0.0001\). 使用survminer创建(Kassambara和Kosinski2018)

运行规范ParticleMDI产生了10万个样本,稀释到每10个样本一个。第一个\(50\%\)的样品因老化而被丢弃。推断的\(\varPhi\)检查数值是否收敛(见补充资料)。总共使用了1024个粒子\(\rho\)设置为0.25。

讨论

我们提出了ParticleMDI,一种利用粒子蒙特卡罗方法对多个数据集进行综合聚类分析的新方法。虽然我们关注的是高斯和分类数据的应用,但我们的方法很容易推广到其他数据类型,在这些数据类型中可以分析地导出后验预测分布,包括具有完全协方差结构的多元高斯模型。我们提出了新的方法来提高我们的算法的计算效率,实验结果表明改进了一个数量级或更高的精度没有损失。我们的方法更广泛地应用于离散状态空间的粒子滤波。我们已经证明了我们的方法在综合数据中发现基本事实的有效性。我们已经将我们的算法应用于来自癌症基因组图谱的真实生物学数据,证明了我们的方法能够推断出临床上有意义的亚组,如其中包含的患者的显著不同的生存情况所示。这些组的成员不受单个数据集中变量的控制,这突出了使用综合方法的重要性。ParticleMDI可以作为统计编程语言Julia(Cunningham et al。2019)

工具书类

  1. Andrieu C,Doucet A,Holenstein R(2010)粒子马尔可夫链蒙特卡罗方法。J R Stat Soc Ser B Stat方法72(3):269–342

    数学网 数学 文章 谷歌学者 

  2. Bernardo JM,Smith AF(2001)贝叶斯理论

  3. Bouchard-CôA,Doucet A,Roth A(2017)《混合模型中贝叶斯推断的粒子吉布斯分裂合并抽样》。机械学习研究杂志18(28):1–39

    数学网 数学 谷歌学者 

  4. 肖邦N(2002)静态模型的序贯粒子滤波方法。生物计量学89(3):539–552

    数学网 数学 文章 谷歌学者 

  5. 肖邦N,Singh SS(2015)关于颗粒吉布斯取样。伯努利21(3):1855-1883

    数学网 数学 文章 谷歌学者 

  6. Cunningham N,Griffin JE,Wild DL,Lee A(2019年),《贝叶斯统计:新挑战和新世代》,2018卷,Springer

  7. Doucet A,Johansen AM(2009)《粒子过滤和平滑教程:十五年后》。Handb非线性滤波器12(656–704):3

    数学 谷歌学者 

  8. Fearnhead P(2004)具有未知数量组分的混合物模型的粒子过滤器。统计计算14(1):11–21

    数学网 文章 谷歌学者 

  9. Fisher-RA(1936)在分类学问题中多重测量的使用。安尤根7(2):179–188

    文章 谷歌学者 

  10. Fritsch A,Ickstadt K等人(2009)改进了基于后验相似矩阵的聚类标准。贝叶斯分析4(2):367–391

    数学网 数学 文章 谷歌学者 

  11. Gabasova E,Reid J,Wernisch L(2017)Clusternomics:异构数据集的集成上下文相关聚类。公共科学图书馆计算机生物学13(10):e1005781

    文章 谷歌学者 

  12. Green PJ,Richardson S(2001)具有和不具有Dirichlet过程的异质性建模。Scand J Stat 28(2):第355–375页

    数学网 数学 文章 谷歌学者 

  13. Griffin J(2014)具有独立增量先验的归一化随机测度混合物的序贯蒙特卡罗方法。统计计算27(1):131–145

    数学网 数学 文章 谷歌学者 

  14. Hol JD,Schon-TB,Gustafsson F(2006),粒子滤波器的重采样算法。非线性统计信号处理研讨会,2006年IEEE,IEEE,第79-82页

  15. Ishwaran H,Zarepour M(2002)Dirichlet过程的精确和近似和表示。Can J Stat 30(2):269–283

    数学网 数学 文章 谷歌学者 

  16. Kassambara A,Kosinski M(2018)survminer:使用“ggplot2”绘制生存曲线。R包版本(4):2

  17. Kirk P,Griffin JE,Savage RS,Ghahramani Z,Wild DL(2012)集成多个数据集的贝叶斯相关聚类。生物信息学28(24):3290–3297

    文章 谷歌学者 

  18. Lawlor N,Fabbri A,Guan P,George J,Karuturi RKM(2016)multiclust:r-package for identification biological related clusters in cancer transcriptome profiles.多库:识别癌症转录组中生物相关簇的r-package。癌症Inf 15:CIN-S38000

    文章 谷歌学者 

  19. 李浩,韩德,侯勇,陈浩,陈志(2015)两种交叉生存曲线的统计推断方法:方法比较。公共科学图书馆一号10(1):e0116774

    文章 谷歌学者 

  20. Liu JS,Chen R(1995)序列插补盲反褶积。美国统计协会杂志90(430):567–576

    数学网 数学 文章 谷歌学者 

  21. Lock EF,Dunson DB(2013)贝叶斯共识聚类。生物信息学29(20):2610–2616

    文章 谷歌学者 

  22. McParland D、Gormley IC、McCormick TH、Clark SJ、Kabudula CW、Collinson MA(2014)使用潜在变量模型根据南非家庭的资产状况对其进行聚类分析。应用统计8(2):747

    数学网 数学 文章 谷歌学者 

  23. McParland D、Phillips CM、Brennan L、Roche HM、Gormley IC(2017)聚类高维混合数据以揭示亚表型:表型和基因型数据的联合分析。统计医学36(28):4548–4569

    数学网 文章 谷歌学者 

  24. Medvedovic M,Yaung K,Bumganner(2004)基于贝叶斯混合模型的复制微阵列数据聚类。生物信息学20(8):1222–1232

    文章 谷歌学者 

  25. Monti S,Tamayo P,Mesirov J,Golub T(2003)共识聚类:基于重采样的基因表达微阵列数据的类发现和可视化方法。马赫学习52(1-2):91-118

    数学 文章 谷歌学者 

  26. Murphy KP(2007)高斯分布的共轭贝叶斯分析。技术代表

  27. Rand-WM(1971)评价聚类方法的客观标准。美国统计协会杂志66(336):846–850

    文章 谷歌学者 

  28. Rasmussen C,de la Cruz B,Ghahramani Z,Wild D(2009)使用Dirichlet过程混合物对基因表达簇中的不确定性进行建模和可视化。IEEE/ACM Trans Comput Biol生物信息6(4):615-628

    文章 谷歌学者 

  29. Rousseau J,Mengersen K(2011)过拟合混合模型中后验分布的渐近行为。J R Stat Soc Ser B Stat方法73(5):689–710

    数学网 数学 文章 谷歌学者 

  30. Savage RS,Ghahramani Z,Griffin JE,Kirk P,Wild DL(2013),通过结合基因组学、转录组学和表观遗传学数据识别胶质母细胞瘤中的癌症亚型。arXiv预印本arXiv:1304.3577

  31. Shen R,Olshen AB,Ladanyi M(2009)使用联合潜在变量模型对多个基因组数据类型进行整合聚类,并应用于乳腺癌和肺癌亚型分析。生物信息学25(22):2906–2912

    文章 谷歌学者 

  32. Steinley D,Brusco MJ(2008),《聚类分析中变量的选择:八种程序的实证比较》。心理测量学73(1):125

    数学网 数学 文章 谷歌学者 

  33. Tarone RE,Ware J(1977)关于生存分布相等的无分布检验。生物计量学64(1):156–160

    数学网 数学 文章 谷歌学者 

  34. Yuan Y,Van Allen EM,Omberg L,Wagle N,Amin Mansour A,Sokolov A,Byers LA,Xu Y,Hess KR,Diao L等人(2014年),评估癌症基因组和蛋白质组数据在不同肿瘤类型中的临床效用。天然生物技术32(7):644

    文章 谷歌学者 

下载参考资料

作者信息

隶属关系

作者

通讯作者

通信对象内森·坎宁安.

附加信息

出版说明

对于出版地图和机构隶属关系中的管辖权主张,Springer Nature保持中立。

电子辅助材料

以下是电子补充材料的链接。

补充材料1(pdf 13710 KB)

权利和权限

开放存取本文根据Creative Commons Attribution 4.0国际许可证获得许可,该许可证允许以任何媒介或格式使用、共享、改编、分发和复制,只要您对原始作者和来源给予适当的信任,提供到Creative Commons许可证的链接,并说明是否进行了更改。本文中的图像或其他第三方材料包含在文章的知识共享许可证中,除非材料的信用额度中另有说明。如果文章的“知识共享”许可证中没有包含材料,并且您的预期用途不受法律法规的允许或超出了允许的使用范围,则您需要直接从版权所有者处获得许可。要查看此许可证的副本,请访问http://creativecommons.org/licenses/by/4.0/.

转载和许可

关于这篇文章

通过十字标记验证货币和真实性

引用这篇文章

Cunningham,N.,Griffin,J.E.和Wild,D.L.ParticleMDI:多数据集聚类分析的粒子蒙特卡罗方法及其在癌症亚型识别中的应用。高级数据分析分类 14日,463–484(2020年)。https://doi.org/10.1007/s11634-020-00401-y

下载引文

关键词

  • 聚类分析
  • 混合模型
  • 贝叶斯推理
  • 粒子蒙特卡罗

数学学科分类

  • 62小时30分
  • 第62页
  • 65摄氏度