总结

在分析变更点过程生成的数据时,一个关键挑战是确定变更点的数量。由于似然函数的不规则性,经典的贝叶斯信息准则(BIC)统计在这里不能很好地工作。通过贝叶斯因子的渐近逼近,我们导出了变漂移布朗运动模型的修正BIC。修改后的BIC与经典BIC相似,因为第一项由对数似然组成,但不同之处在于对模型维数不利的项。作为一个应用示例,该新统计数据用于分析基于阵列的比较基因组杂交(array-CGH)数据。Array-CGH测量细胞样本每个基因组位置的染色体拷贝数,有助于发现肿瘤细胞中基因组缺失和扩增的区域。与现有方法相比,改进的BIC在准确选择拷贝数变化的区域数量方面表现良好。与现有方法不同,它不依赖于调整参数或密集计算。因此,它是公正的,更容易理解和使用。

1.简介

贝叶斯信息准则(BIC)由施瓦茨(1978)是一种流行的模型选择方法。然而,它的使用在理论上并不适合于变点问题,其中似然函数不满足所需的正则性条件。转换点问题中模型选择的一个方面是确定转换点的数量。在本文中,我们将给出由独立正态分布观测值组成的数据的BIC的修改版本,这些观测值具有常方差和分段常均值。

在变点问题中,没有一种标准的模型选择方法。尽管缺乏理论依据,但BIC之前已经应用于本文中处理的形式的模型(例如,参见Yao,1988年;李,2001)。BIC是从后验模型概率的大样本估计中得出的,后验模型的概率不依赖于模型或参数的先验概率分布。最近的模型选择方法采用了基于计算的贝叶斯方法(例如。,George and McCulloch,1993年;格林,1995年)。他们使用蒙特卡罗方法估计后验模型概率,因此涉及大量计算和先验的主观选择。

最近,Siegmund(2004)提出了一种新的方法来获得贝叶斯因子的大样本近似值,该方法绕过了推导BIC时使用的泰勒展开。使用这种方法,他导出了一种新的类BIC统计,用于在遗传连锁研究中绘制数量性状位点,这是一个与在Ornstein-Uhlenbeck过程的平均值中寻找变化点有关的问题。我们在这里采用的方法与Siegmund建议的方法类似,但我们将处理具有更自然的理论公式的问题,并在连续时间内减少到寻找布朗运动漂移中的变化点。

原始BIC和修改后的BIC都可以解释为惩罚似然方法,尽管应该强调的是,“惩罚”是计算的自然结果,而不是避免过拟合的正则化方法。修改后的BIC与其他惩罚似然方法的简单比较如下第2节.

例如,我们将把新版BIC应用于基于阵列的比较基因组杂交(阵列CGH)数据的分析,该数据定量测量沿基因组线性排列的数千个位置的DNA拷贝数。阵列CGH数据分析的目标是准确检测DNA缺失或扩增区域。阵列CGH研究的更多细节以及与方法的比较Olshen等人(2004)Fridyland等人(2004)在中给出第4节5.

2.高斯变点模型的改进BIC

考虑一系列观察结果= (1,2, … ,T型),其中是独立分布的高斯随机变量:

(1)

这里我们将τ称为该过程的转折点。当然,转换点被限制在集合中

出于显而易见的原因,我们添加了以下限制:j≠μj+ 1。我们假设,对于大样本量,转换点之间的距离足够远

(2)

其中0<第页1< ⋯ <第页< 1. 为了便于说明,我们还将定义τ0=0和τ+ 1=T型.

如果改变点的数量,,则导出了基于最大可能性的程序,以查找和测试变化点(简单的一维情况在James中给出,詹姆斯和西格蒙德,1987年)。我们的目标是.从模型选择的角度来看,对于0≤<∞我们定义论坛为上述高斯模型转换点。因此,论坛是最简单的模型,没有改变论坛与线性增加.模型的拟合论坛观测数据当然会随着。在我们的模型选择过程中,我们将遵循找到最简约模型的原则,该模型可以“很好”地解释数据。

如中所示施瓦兹(1978)Siegmund(2004),我们将通过推导贝叶斯因子的渐近近似值来解决这个问题。这里值得注意的是,尽管以下定理已被证明适用于特定的先验分布,但这些结果适用于非常大的一类先验分布。上一个τ必须采用以下形式如果(τ) =(τ)T型负极,使用C类1<最大值τ(τ) <C类2,其中C类1,C类2是不依赖于均值和方差的先验密度(当方差未知时)必须有一个包含最大似然估计量的开集的支持。然后,作为样本大小T型增加,似然项在贝叶斯因子中占主导地位,先验项凹进到可忽略的余项中。

为了简单起见,我们首先假设方差σ2已知,在不损失通用性的情况下,将其设为1。重新参数化后论坛、模型参数论坛论坛属于子参数空间

(3)

下一个定理给出了这种情况下贝叶斯因子的大样本近似值。

定理1:  论坛 是中定义的模型 方程式(1) .让σ= 1并假设(,θ)先穿制服 论坛,然后

(4)

哪里

该定理的证明简图将在支持的Web附录中给出。渐近地,选择具有最大后验概率的模型相当于选择具有最大贝叶斯因子的模型。我们将调用的右侧方程式(4),其中我们忽略了O(运行)第页(1) 余数项,即“修改后的BIC”。我们建议作为一个模型选择程序,选择最大化修改后BIC的模型。

很容易看出,修改后的BIC的第一项是在论坛可以解释右侧的剩余术语(负数)方程(4)作为对模型复杂性的“惩罚”。经典的BIC程序也可以进行类似的解释。与经典BIC和惩罚似然法相比(例如。,Tibshirani,1996年;Birgéand Massart,2001年;顾和王,2003;拉维耶,2005年),修改后的BIC在处罚条款上有所不同。所有惩罚似然方法都选择最大化形式准则的模型

哪里是的对数似然函数论坛是最大似然参数估计值,以及第页是惩罚函数。在经典的BIC中,第页=d日日志T型/2,其中d日是模型的尺寸论坛.在一些其他方法中(例如。,Tibshirani,1996年;Birgéand Massart,2001年;顾和王,2003),第页取决于收缩参数β,该参数必须由用户选择,通常通过交叉验证。拉维耶(2005)提出第页= β日志T型并给出了选择β的数据自适应方法。

与经典BIC一样,修改后的BIC为第页,无需用户选择收缩率参数。然而,虽然在BIC中,每个参数为模型贡献了一个维度,但对于修改后的BIC,惩罚项不再遵循简单的公式。如果我们重新参数化论坛,其中0=第页0<第页1< ⋯ <第页+ 1=1,则从方程(4)可以写为

(5)

第二学期方程式(5)当变化点均匀分布在(0,T型),当它们尽可能靠近时最小化(例如。,论坛)以下为:

(6)
(7)

因为假设(2),方程式(7)在我们的模型下是不可能的,因此修改后的BIC的惩罚项的粗略近似为第页≈ 3日志(T型)/2. 然而,由于原木生长缓慢(T型),即使对于大型数据集方程式(5)仍然可以对处罚做出重大贡献。方程式(5),我们将为每个“跳跃”参数δ赋予一个维度。然后,我们在一维和二维之间为每个转换点位置参数做出贡献,准确的贡献取决于它们之间的相对位置。这与变点模型维度的实际经验一致(例如。,Birgéand Massart,2001年;拉维耶,2005年).

在应用中,对于未知方差的情况,我们需要一个修改的BIC版本,这在下一个定理中给出。

定理2: 在统一的前提下{μ, τ, δ : = 1, … ,},和非信息性先验1/σ对于σ,

(8)

哪里

(9)

 论坛,论坛 定义见定理1

第二个定理中的近似形式类似于定理1中的近似。第一项仍然是最大对数似然比,在方差未知的情况下,密度是Student的t吨而不是高斯。处罚条款,论坛,保持不变。近似值第二行中的其他项都是有序的日志(T型),来自于对扰动参数σ的积分。有关此定理的证明,请参见张(2005).

3.实施细节

这里我们详细描述了修改后的BIC如何应用于数据{1, … , T型}. 假设方差未知。如果每个固定达到某个预冲最大值M(M),最大似然估计论坛计算了假定的变点位置,然后可以直接应用定理2中的公式计算修正的BIC(m)= 1, … , M(M)。然后我们选择*最大化修改后的BIC(m)。我们的方法返回的分段是论坛.最大值M(M)是一种算法上的便利,可以避免不必要的计算。修改后的BIC是一致的,只要M(M)大于真正的更改点数。根据我们的经验,BIC(m)增大到最大值,然后减小。我们试着设置M(M)也许是迭代的,以便我们清楚地看到这个模式。

剩下的问题是如何计算论坛对于给定的虽然原则上希望考虑{1, … , T型}按点t吨1<t吨2⋯ <t吨,除了非常小的值外,这在计算上对所有值都是不可行的.因为递归方法Olshen等人(2004)是非常快速的,并且从他们的示例中可以非常有效地检测到真正的变化点(也许还有错误的变化点),我们使用了他们的循环二进制分割(CBS)算法的以下修改(并且有类似的经验)。支持Web附录中给出了一些定义和技术细节。

  • 1

    论坛.

  • 2

    While期间论坛

    • (a)

      对于每个论坛,请执行以下操作:

      • 如果论坛:让X是方波变化的最大似然统计量(方程式(3)在Web附录中){ 论坛, … ,  论坛},并让,s′是最大可能的变更点位置。然后,

        • 如果在中作为转换点不重要{ 论坛, … ,  论坛},让第页=s′X是单个变化的似然比第页;

        • 否则,如果s′在中作为转换点不重要{ 论坛, … ,  论坛},让第页=X是单个变化的似然比第页;

        • 否则,让第页= {,s′}和X=X′.

      • 否则,当论坛X是单个变化的对数似然比统计(方程式(2)在Web附录中){ 论坛, … ,  论坛},并让第页是最大可能的变点位置。

    • (b)

      选择拥有最大的X,并添加第页到有序列表论坛.

  • 三。

    丢弃0,T型论坛.

在上述第2(a)步中,假定变化点的重要性s′可以通过固定的p值阈值进行评估,我们将其设置为0.05。该算法对该阈值具有鲁棒性。有关这些似然比统计数据及其重要性计算的描述,请参阅Web附录。上述算法与原始CBS算法的不同之处在于,当一个预定的数字出现时,递归停止已找到个更改点,而不是在无法检测到更多“重要”更改点时。此外,如Olshen等人(2004),我们在步骤2(a)中纠正了以下情况中的边缘效应s′靠近序列的边缘。虽然此算法不能保证是最优的,但根据我们的经验,它是准确且快速的。

4.仿真研究

我们首先使用模拟数据来测试修改后的BIC作为模型选择标准的合理性。为了评估性能,我们将使用经典BIC和Olshen等人(2004年)Fridlyand等人(2004)。与后两种方法进行比较的一个原因是,它们的设计都考虑到了阵列CGH数据的分析,我们将在第5节。我们首先简要描述这两种方法的工作原理。

Olshen等人(2004)开发了一种改进的二进制分割算法,称为CBS,用于估计变化点位置。CBS几乎与中给出的算法完全相同第3节,除了它对递归使用了不同的停止标准:如果没有一个分段包含显著的变化,则停止搜索,这是基于一些预冲p值阈值α。然而,他们发现这种停止标准有高估变化点数量的倾向,因此在CBS算法的末尾添加了一个修剪步骤。修剪步骤取决于阈值参数γ。由于改变γ可以将报告的变化点的数量限制为任何值,并且对于单个研究,没有明确的公正的选择γ的方法,因此我们将只考虑基本的CBS算法和阈值停止规则来进行选择.

Fridlyand等人(2004)使用隐马尔可夫模型(HMM)分割阵列CGH数据。在这种方法中,每个时间点的未知平均值是隐藏状态,而变化点是不同状态之间的转换。找到转换点需要首先对HMM进行数据训练,然后计算最可能的状态序列。该方法中的一个关键参数是HMM中的状态数,它是使用Akaike信息准则(AIC)选择的。对于阵列CGH数据分析,Fridland等人发现最好在AIC之上添加一个合并步骤,以进一步减少状态数。就像修剪一样Olshen等人(2004)该状态合并步骤可以调整为使该方法任意保守,并且没有明确的选择修剪参数的标准。因此,我们将使用HMM方法和AIC进行模型选择。

我们使用的模拟数据集是由Olshen等人(2004).模型如下:

(10)

其中εN个(0, 0.04),T型=497,μ更改六次τ= (138, 225, 242, 299, 308, 332), μ0=−0.18,以及δ= (0.26, 0.99, 1.6, 0.69, 0.85, 0.53). 第二任期方程式(10)是一个正弦趋势分量,旨在使该数据集更具挑战性,最初由添加Olshen等人(2004)模拟阵列CGH数据中存在的周期性趋势。正如他们的论文中所述,我们将对趋势参数进行实验∈{0,0.01,0.025}分别对应于长周期无趋势和局部趋势。使用该模型生成的数据系列的示例以及所使用的每个趋势参数如所示图1.

图1

模拟模型的一个蒙特卡罗样本示例(10)。趋势参数为顶部、中部和底部数字分别为0、0.025和0.01。这条线描绘了加在平均值上的趋势。

表1显示了每种不同方法对模型生成的仿真数据的测试结果(10)。该表列出了每种方法在100次蒙特卡洛复制中找到的时间百分比转换点,带= 5, 6, 7, 8. 从这些结果可以清楚地看出,改进后的BIC比传统BIC更具特异性,同时又不牺牲灵敏度。这是意料之中的,因为经典的BIC具有错误的惩罚条款,对更改点参数的惩罚不够。结果还表明,修正BIC的未知方差版本比已知方差版本更保守。这是由于对数函数在第一项中的阻尼效应方程式(8)因此,与BIC的方差已知版本相比,罚款条款的权重相对较大。比较方程式(4)(8)注意,根据大数定律,不锈钢背景2Tx(发送)对一些人来说x个这取决于δ第页因此(T型−1)日志(1+不锈钢背景/不锈钢重量(wg)) T型日志(1+x个)因此方程式(8)就第一项而言方程式(4)按对数(1)+x个)。对于更具统计意义的变化,这种影响更为显著。

表1

对不同趋势参数值的模拟数据的测试结果。这些列列出了每种方法找到5、6、7或>8个更改点的100次模拟的次数。正确的数字是6。

趋势方法5678个+
传统BIC077149
Olshen等人(2004年)0916
Fridlyand等人(2004)083134
修正的BIC(定理1)09550
修正的BIC(定理2)0970
传统BIC0583012
Olshen等人(2004)0662311
Fridlyand等人(2004)07027
改进的BIC(定理1)088111
修正的BIC(定理2)09820
传统BIC0383527
Olshen等人(2004)0681121
Fridlyand等人(2004)180136
修正的BIC(定理1)07819
修正的BIC(定理2)29440
趋势方法5678个+
传统BIC077149
Olshen等人(2004)0916
Fridlyand等人(2004)083134
修正的BIC(定理1)09550
修正的BIC(定理2)0970
传统BIC0583012
Olshen等人(2004)0662311
Fridland等人(2004年)07027
修正的BIC(定理1)088111
修正的BIC(定理2)09820
传统BIC0383527
Olshen等人(2004)0681121
Fridlyand等人(2004)180136
修正的BIC(定理1)07819
修正的BIC(定理2)29440
表1

对不同趋势参数值的模拟数据的测试结果。这些列列出了每种方法找到5、6、7或>8个更改点的100次模拟的次数。正确的数字是6。

趋势方法5678个+
传统BIC077149
Olshen等人(2004)0916
Fridlyand等人(2004)083134
修正的BIC(定理1)09550
修正的BIC(定理2)0970
传统BIC0583012
Olshen等人(2004年)0662311
Fridlyand等人(2004)07027
修正的BIC(定理1)088111
修正的BIC(定理2)09820
传统BIC0383527
Olshen等人(2004)0681121
Fridlyand等人(2004)180136
改进的BIC(定理1)07819
修正的BIC(定理2)29440
趋势方法5678个+
传统BIC077149
Olshen等人(2004)0916
Fridlyand等人(2004)083134
修正的BIC(定理1)09550
修正的BIC(定理2)0970
传统BIC0583012
Olshen等人(2004)0662311
Fridlyand等人(2004)07027
修正的BIC(定理1)088111
修正的BIC(定理2)09820
传统BIC0383527
Olshen等人(2004)0681121
Fridland等人(2004年)180136
修正的BIC(定理1)07819
修正的BIC(定理2)29440

这些仿真结果表明,改进的BIC算法在不增加剪枝步骤的情况下,性能优于基本的CBS算法,当加入正弦趋势时,改进效果更加明显。在这个数据集上,修改后的BIC似乎比基本HMM方法更可取Fridlyand等人(2004)然而,我们的模拟研究表明,HMM方法的结果取决于模型参数的初始化(参见支持Web附录)。我们选择不与修剪的CBS或以状态为目标的HMM方法进行比较,因为没有系统或明显的方法来选择修剪(或状态合并)阈值。

5.阵列CGH数据分析

对于均质细胞系,给定基因组位置的DNA拷贝数是每个细胞在该位置的基因组DNA拷贝数。在人类和其他二倍体生物体中,常染色体的DNA拷贝数通常为2。然而,已经观察到肿瘤细胞的基因组显示出拷贝数扩增或减少的区域。准确识别这些区域对研究肿瘤进展很重要,也可能有助于肿瘤分类(Albertson等人,2003年).

最近开发的阵列CGH技术提供了同时定量测量基因组上数千个位置的DNA拷贝数的方法(Pinkel等人,1998年;Pollack等人,1999年)。阵列CGH实验背后的技术与cDNA基因表达实验背后的技术相似。通常,测试基因组DNA池(例如,来自肿瘤细胞样本的基因组DNA)和二倍体参考基因组DNA池用染料(例如,Cy3和Cy5)进行不同标记。这两个染料标记的样本混合并杂交到一个微阵列芯片上,该芯片上发现了数千个基因组靶点,每个靶点都映射到一个已知的基因组位置。然后扫描杂交微阵列芯片,并使用图像分析软件计算每个基因组靶点的测试和参考荧光强度的比率。染料的强度比是用染料标记的DNA样品丰度比的替代品。为了我们当前的数据分析目的,我们可以将阵列CGH实验视为为每个细胞样本生成一个有序的序列(t吨,t吨)对,其中t吨代表基因组中的位置t吨表示日志2来自该位置的基因组靶点的测试与参考点强度的比率。

阵列CGH数据分析的第一步是准确确定拷贝数变化的区域。这个问题有两个主要方面:像差位置的准确估计和像差数量的确定。基于其在Coriel数据集上的性能Snijders等人(2001年)现有方法(例如。,Fridlyand等人,2004年;Olshen等人,2004年;Hsu等人,2005年;Wang等人,2005年)一旦仔细选择了像差的数量,似乎或多或少都同意像差的位置。然而,分段数的估计涉及建模的更微妙方面,是当前方法之间分歧的主要来源。研究者Picard等人(2005年)对阵列CGH数据的几种不同模型选择标准进行了比较,并说明了这一步骤对整体分割结果的重要性。我们提出了改进的BIC作为阵列CGH数据分析的模型选择程序。尽管给出了算法第3节将修改后的BIC与CBS算法配对,它可以作为模型选择过程来补充大多数分割算法。

由于探针按基因组定位的线性排序和染色体畸变的不连续性,阵列CGH数据的拷贝数估计问题自然适合于变化点检测框架。我们将为数据假设以下高级模型

(11)

哪里如果是离散拷贝数的单调递增函数c(c)t吨在位置t吨和εt吨具有共同方差σ的独立同分布高斯噪声2.引入原因如果对数荧光强度比的平均值随拷贝数比的对数增加而增加,但不相等,这可能是由于污染或镶嵌所致(Pollack等人,1999年)。此模型类似于Olshen等人(2004)。假设共同方差的动机可以在Picard等人(2005年),虽然我们的方法也可以适应更一般的情况,即方差不同,但任意两个方差的比率是已知的。

与之前的研究一样,我们首先使用了Snijders等人(2001年)测试修改后的BIC。该数据集包括对15个成纤维细胞系(其中我们将使用Olshen等人,2004年)。由于该数据集由染色体拷贝数发生变化的纯二倍体细胞系组成,这些细胞系以前以光谱核型分析为特征,肉眼很容易检测到,因此它主要用作原理证明。根据假阳性和假阴性的数量,我们的方法的结果列在表2并与使用Olshen等人(2004年)Fridlyand等人(2004)请注意修剪步骤Olshen等人(2004)未应用。根据细胞系GM03563和GM01750的数据,绘制出每种列出方法所选模型的拟合平均值图2从这些结果中,我们可以看到修改后的BIC在这个简单的数据集上表现合理,比Olshen等人(2004)同时不影响灵敏度。HMM算法Fridlyand等人(2004)也表现良好,但它涉及用户选择训练参数和阈值,因此可能不是一个公平的比较。仔细检查图2表明Olshen等人(2004年)我们的结果是由数据中出现的弯曲局部趋势引起的。Olshen等人认为,这些局部趋势可能有生物起源。

表2

来自Snijders等人(2001年)对于的CBS算法Olshen等人(2004)(阈值p值=0.01,无修剪),HMM算法Fridlyand等人。(2004;状态合并阈值0.35),定理2的修正BIC,以及定理3的有限层模型的修正BIC

Olshen等人。Fridlyand等人。修改BIC(定理2)修改BIC(定理3)
细胞系FP公司FN公司FP公司FN公司FP公司FN公司FP公司FN公司
GM03563型140005020
GM05296号110404040
GM01750号80202020
GM03134号120406060
GM13330公司170005020
GM01535号60202020
GM07081号61212121
通用汽车1303170207010
GM01524号90002000
Olshen等人。Fridlyand等人。修改BIC(定理2)修改BIC(定理3)
细胞系FP公司FN公司FP公司FN公司FP公司FN公司FP公司FN公司
GM03563号140005020
GM05296号110404040
GM01750号80202020
GM03134号120406060
GM13330公司170005020
克0153560202020
GM07081号61212121
GM13031号70207010
GM01524号90002000
表2

来自Snijders等人(2001年)对于的CBS算法Olshen等人(2004)(阈值p值=0.01,无修剪),HMM算法Fridlyand等人。(2004;状态合并阈值0.35),定理2的修正BIC,以及定理3的有限层模型的修正BIC

Olshen等人。Fridlyand等人。修改BIC(定理2)修改BIC(定理3)
细胞系FP公司FN公司FP公司FN公司FP公司FN公司FP公司FN公司
GM03563号140005020
GM05296号110404040
GM01750号80202020
GM03134号120406060
GM13330公司170005020
GM01535号60202020
克0708161212121
GM13031号70207010
GM01524号90002000
Olshen等人。Fridlyand等人。修改BIC(定理2)修改BIC(定理3)
细胞系FP公司FN公司FP公司FN公司FP公司FN公司FP公司FN公司
GM03563号140005020
GM05296号110404040
GM01750号80202020
GM03134号120406060
通用13330170005020
GM01535号60202020
GM07081号61212121
GM13031号70207010
GM01524号90002000
图2

Coriel数据集中GM03563和GM01750细胞系拟合平均值的比较。第1行:没有修剪的CBS算法。第2行:状态合并阈值设置为0.35的HMM算法。第3行:修改了定理2的BIC。第4行:修改定理3的BIC,选择转换点数量和级别数量。

为了进一步测试修改后的BIC,我们使用了来自Snijders等人(2003年)该数据集背后的实验研究了错配修复(MMR)基因的缺陷如何影响肿瘤的染色体畸变谱。在我们的测试中,我们将使用BT474和MCF7细胞系,这两种细胞系都具有复杂的畸变特征。本实验的真实染色体拷贝数未知,因此无法直接验证结果。我们在该数据集上测试修改后的BIC的目标有两个:首先,通过肉眼验证修改后的BIC在更复杂的数据集上的表现是否合理;其次,在更具挑战性的场景中,比较修改后的BIC和HMM方法的行为。

这个更复杂的数据集表明,我们的程序是一个比HMM更保守的变化点检测程序。对于BT474细胞系,修改后的BIC在39个变化点处达到峰值,而对于HMM,如果不采取状态合并步骤,则发现103个变化点。状态合并后Fridlyand等人(2004)报告了该数据系列的69个变化点。对于MCF7细胞系,未合并的HMM发现102个变化点,在合并步骤后减少到95,而修改的BIC在55处达到峰值。这两个数据系列及其通过每种方法估计的平均值如所示图34(图中绘制了下一节的结果,就发现的变化点数量而言,这与本节的结果非常相似。)由于真实拷贝数未知,我们无法从这些结果中判断哪种方法更“准确”。然而,通过仔细检查图,很明显,这三种方法产生的分割效果截然不同。这表明,虽然分割方法很容易就Coriel细胞系等简单数据集达成一致,但对于更复杂的数据集,统计方法的选择至关重要。

图3

MMR数据集中BT474细胞系拟合平均值的比较。顶部面板:使用中提出的HMM算法估计的平均值Fridlyand等人(2004),而没有状态合并步骤。中间面板:与顶部面板相同,但合并步骤的阈值参数设置为0.35。底部面板:定理3中的统计信息用于选择转换点和级别的数量;使用CBS算法发现的变化点。

图4

MMR数据集中MC57细胞系拟合平均值的比较。顶部面板:使用中提出的HMM算法估计的平均值Fridlyand等人(2004),没有状态合并步骤。中间面板:与顶部面板相同,但合并步骤的阈值参数设置为0.35。底部面板:定理3中的统计信息用于选择转换点和级别的数量;使用CBS算法发现的变化点。

6.限制变化幅度

在数组CGH数据中,基础染色体拷贝数仅以整数的形式变化。因此,对于许多数据集,基础段意味着{t吨}应采用有限数量的值。具有相同基本拷贝数组成的基因组不同片段也应获得相同的估计平均值,而不是每个片段的估计平均数略有不同。定理1和2中修改后的BIC没有考虑这些限制。按型号(1),在每个更改点,进程平均值可以跳到任何值。因此,如果变化点,有可能,也肯定会有最大的可能性,+1不同水平的估计平均值。

为了为数组CGH数据创建更具体的模型,我们希望限制平均值的跳跃。在这里,我们提出了有限层模型论坛,其中每个变化点平均值必须跳到k个水平:

(12)

哪里j∈{0, … , k个}是段的级别分配在转换点之间t吨−1t吨。我们强制执行0=0和δ0=0,以便在第一个变化点之前,平均值处于基线水平μ。当然,总层数,k个+1,必须不大于部分的总数,+ 1. 什么时候?k个=我们有以前的无限制跳跃模型。为了使这个新问题的可能性最大化,我们不仅要搜索可能的变化点,还要搜索段到级别的可能赋值.给定,k个,可能值的域

(13)

由于模型的尺寸随着k个,我们现在有一个二维模型选择问题。以下定理给出了BIC的修改版本,用于同时选择两个模型k个当方差未知时。

定理三:  论坛 是有m个变化点和k个不同水平的模型。然后,

(14)

哪里

这个定理的证明与定理2的证明非常相似,将被省略。的结构方程式(14),我们称之为修改后的BIC转换点和k个级别,几乎与定理2中修改的BIC的结构完全相同。关键的变化是,在定理2中,两个相邻变化点之间的每个线段在计算不锈钢背景不锈钢重量(wg),此处是分配给同一级别的所有段论坛属于同一组。在剩余条款的某些部分中方程式(14),被替换为k个因为δ现在是k个-维度。特别是处罚条款论坛现在正常k个日志(T型)。因此,惩罚期限随着变化点数量和级别数量的增加而增加。

在Coriel和MMR数据集上测试了有限层模型的修正BIC。与之前一样,使用CBS算法查找变化点,然后将分段平均值聚类为k个组使用k个-意味着。然而,由于跳跃大小的限制,CBS经常会导致以下方面的次优选择论坛,但这里使用它的速度和简单性。这一新模型的性能与其他模型一起列在表2图1对于Coriel数据集,极限水平模型的性能优于定理2的无限制跳跃模型,但不能完全解决局部趋势的问题。对于MMR数据集,有限级别模型并没有显著减少检测到的更改点数量,而是显著减少了级别数量。

7.讨论

对模拟和阵列CGH数据的测试结果表明,在本文讨论的问题类型中,修正的BIC是一个合理的模型选择统计量。我们使用的模拟数据来自Olshen等人(2004),从而对我们的方法进行公正的评估。在Coriel数据集中,可以看出修改后的BIC比来自Olshen等人(2004)更稳健地应对本地趋势,这似乎是阵列CGH数据集中的一个问题。在更复杂的MMR数据集中,修改后的BIC明显比Fridlyand等人的方法更保守,并且对数据的拟合更平滑。与其他确定阵列CGH数据中像差数量的方法相比,改进的BIC的优点之一是其公正性,因为与经典BIC一样,它不需要特定的先验或调谐参数。因此,修改后的BIC是一个自然的模型选择标准,与CBS算法的修改一起使用Olshen等人(2004)用于阵列CGH数据分析。

我们尚未考虑的一类通用模型选择方法是蒙特卡罗方法,用于估计模型的后验概率(例如。,格林,1995年)。进行比较研究是很合理的。当它被发明时,BIC具有计算优势,因为缺乏计算能力使得确定p值和后验概率成为一项艰巨的任务。现在计算不再是一个问题,但这些基于蒙特卡罗的方法涉及选择先验参数和调整参数,以及蒙特卡罗估计的收敛问题。因此,它们的透明度较低,如果没有系统的稳健性研究,人们可能会合理地担心,它们与修改后的BIC之间的差异反映了可能是任意的选择,而不是数据中的信息。

在我们的分析中,有三个关键的建模假设,即数据是正态分布、同方差和不相关的。如上所述,同方差的假设可以稍稍减弱,这对于方差随某些协变量系统性变化的数据可能有好处。序列相关性可能是一个严重的问题,但纠正它似乎相当简单。这对于下一代阵列CGH数据可能会变得很重要。我们分析的自然背景是指数分布族,因此正态性假设可以被其他特定分布族所取代。对于适度偏离常态该怎么办(如果有的话)的问题值得研究。阵列CGH数据的尾部似乎比正常情况稍长,尽管稍微不正确的分割可能会夸大残差的尾部。在似然比涉及观测值和的情况下,人们可以期望中心极限定理使正态性假设合理,前提是常数平均值的分段相当长;但我们不能(特别是对于阵列CGH数据)排除非常短片段的可能性。

第6节,我们给出了修改后的BIC的一个版本,它可以选择级别数和更改点数。在定理3中,改变点之间每个段的级别分配是未知的,但结果很容易适应这种信息已知的更简单的情况。例如,我们可能希望限制更改,以便流程在跳转到更改状态后始终返回到基线平均水平。这种情况下修改后的BIC与方程式(14),但是k个将是的函数“基线”段的级别分配将被固定。这样,修改后的BIC在适应模型对变化幅度的限制方面非常灵活。然而,对于这些受限的变更模型,CBS算法不再适合作为搜索方法来查找最佳变更点,因为它没有考虑级别分配。

修改后的BIC和原始BIC之间的一个不同之处是,修改后BIC的余项是一个有界随机变量,而原始BIC的剩余项是有界常数。因此,修改后的BIC不如原始BIC稳定。如中所示施瓦茨(1978)我们省略了修改后的BIC中贝叶斯因子近似值的所有项,这些项不随样本大小增长。然而,在某些情况下,包括订单条款可能是有益的日志k个包含这些术语涉及任意性,并带来一定程度的偏见。包括哪些术语以及在什么情况下可以证明这些术语是合理的,这是未来调查的一个有趣的主题。

最后,我们想强调的是,改进的BIC是一种通用的模型选择方法,它不仅限于阵列CGH数据的分析。

8.补充资料

中过程的Matlab代码第3节可以通过电子邮件从作者处获得。支持Web附录位于http://www.tibs.org/生物特征.

鸣谢

第一作者的研究得到了国防科学与工程研究生奖学金的支持。第二位作者的研究得到了国家科学基金会和国家卫生研究院的资助。

工具书类

阿尔伯森
,
D.J.博士。
,
柯林斯
,
C、。
,
麦考密克
,
F、。
、和
灰色
,
J·W·。
(
2003
).
实体瘤的染色体畸变
.
自然遗传学
 
34
,
369
376
.

比尔热
,
L。
马萨尔
,
第页。
(
2001
).
高斯模型选择
.
欧洲数学学会杂志
 
,
203
268
.

弗里德兰
,
J。
,
斯尼德斯
,
A。
,
平克尔
,
D。
,
艾伯森
,
D.克。
、和
贾恩
,
答:N。
(
2004
).
隐马尔可夫模型在阵列CGH数据分析中的应用
.特殊基因组问题
多元分析杂志
 
90
,
132
153
.

乔治
,
E.I.公司。
麦卡洛赫
,
R.E.公司。
(
1993
).
通过吉布斯采样选择变量
.
美国统计协会杂志
 
88
,
881
889
.

绿色
,
P.J.公司。
(
1995
).
可逆跳跃马尔可夫链蒙特卡罗计算与贝叶斯模型确定
.
生物特征
 
82
,
711
732
.

,
C、。
,
J。
(
2003
).
惩罚似然密度估计:直接交叉验证和可扩展近似
.
中国统计局
 
13
,
811
826
.

,
L。
,
自我
,
S.G.公司。
,
树林
,
D。
,
伦道夫
,
T。
,
,
K。
,
德尔罗
,
J·J。
,
,
L。
、和
搬运工
,
第页。
(
2005
).
基于阵列的比较基因组杂交数据的小波去噪
.
生物统计学
 
6
,
211
226
.

詹姆斯
,
B。
,
詹姆斯
,
K.L.公司。
、和
西格蒙德
,
D。
(
1987
).
变更点测试
.
生物特征
 
74
,
71
84
.

拉维埃尔
,
M。
(
2005
).
对变化点问题使用惩罚对比度
.
信号处理
 
85
,
1501
1510
.

,
西。
(
2001
).
DNA分割作为模型选择过程
.英寸
第五届国际计算生物学会议记录
,
T。
 
伦高(Lengauer)
,
D。
 
桑苛夫
,
美国。
 
伊斯特拉尔
,
第页。
 
佩夫兹纳
、和
M。
 
沃特曼
(编辑),
204
210
.
纽约
:
计算机协会
.

奥尔申
,
答:B。
,
文卡特拉曼
,
E.S.公司。
,
卢西托
,
R。
、和
威格勒
,
M。
(
2004
).
基于阵列的DNA拷贝数数据分析的循环二值分割
.
生物统计学
 
5
,
557
572
.

皮卡德
,
F、。
,
罗宾
,
美国。
,
拉维埃尔
,
M。
,
瓦伊斯
,
C、。
、和
多丹
,
J。
(
2005
).
阵列CGH数据分析的统计方法
.
BMC生物信息学
 
6
,
27
.

平克尔
,
D。
,
Seagraves公司
,
R。
,
苏达尔
,
D。
,. (
1998
).
利用比较基因组杂交技术对微阵列进行DNA拷贝数变异的高分辨率分析
.
自然遗传学
 
20
,
207
211
.

波拉克
,
J.R.公司。
,
佩罗
,
C.M.公司。
,
阿利扎德
,
答:A。
,
艾森
,
医学学士。
,
佩加门斯基科夫
,
A。
,
威廉姆斯
,
成本加运费。
,
杰弗里
,
第S条。
,
博茨坦
,
D。
、和
棕色
,
采购订单。
(
1999
).
使用cDNA微阵列对DNA拷贝数变化进行全基因组分析
.
自然遗传学
 
23
,
41
46
.

施瓦兹
,
G.公司。
(
1978
).
估算模型的维数
.
统计年刊
 
6
,
461
464
.

西格蒙德
,
D。
(
1992
).
随机场极大值的尾部近似
.英寸
概率论:1989年新加坡概率会议论文集
,
L·H·Y。
 
,
K.P.公司。
 
,
K。
 
、和
J.-H.公司。
 
(编辑),
147
158
.
柏林
:
德格鲁特
.

西格蒙德
,
D。
(
2004
).
不规则问题中的模型选择:在绘制数量性状基因座中的应用
.
生物特征
 
92
,
785
800
.

斯尼德斯
,
上午。
,
诺瓦克
,
N。
,
Segraves公司
,
R。
,. (
2001
).
用于全基因组DNA拷贝数测量的微阵列组装
.
自然遗传学
 
29
,
263
264
.

斯尼德斯
,
上午。
,
弗里德兰
,
J。
,
曼斯
,
D.A。
,
Segraves公司
,
R。
,
贾恩
,
答:N。
,
平克尔
,
D。
、和
阿尔伯森
,
D.克。
(
2003
).
通过不稳定性和选择形成肿瘤和耐药基因组
.
癌基因
 
22
,
4370
4379
.

提比什拉尼
,
R。
(
1996
).
通过套索回归收缩和选择
.
英国皇家统计学会杂志B辑
 
58
,
267
288
.

,
第页。
,
基姆
,
年。
,
波拉克
,
J。
,
多复变函数
,
B。
、和
提比什拉尼
,
R。
(
2005
).
一种在阵列CGH数据中调用增益和损耗的方法
.
生物统计学
 
6
,
45
58
.

姚明
,
Y.C.(纽约)。
(
1988
).
利用Schwarz准则估计转换点的数量
.
统计与概率信件
 
6
,
181
189
.

,
N.R.(不适用)。
(
2005
).
变点检测和序列比对:基因组学的统计问题
.
加州斯坦福大学统计系博士论文
.

本文根据牛津大学出版社标准期刊出版模式的条款出版和发行(https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)