一种脑功能连接分析的稀疏贝叶斯网络与Granger双约束方法与流程

本发明涉及fmri数据的功能网络连接分析，特别是涉及一种脑功能连接分析的稀疏贝叶斯网络与granger双约束方法。

背景技术：

功能磁共振成像(functionalmagneticresonanceimaging，fmri)是一种神经影像学技术，具有高空间分辨率、可重复检测、无创伤性等优点，已成为脑功能研究和脑疾病诊断的重要工具之一。数据驱动的独立成分分析(independentcomponentanalysis,ica)方法能够从fmri数据中提取出被试(subject)的脑空间激活区(spatialmap,sm)成分及其对应的时间过程(timecourse,tc)成分。基于tc成分进行的脑功能连接(functionalconnectivity)分析，可探究不同脑区之间相互作用与协调的模式，发现脑疾病(如精神分裂症、阿尔兹海默病、抑郁症、躁郁症等)患者与健康对照被试在脑功能连接方面的显著异常，进而用于脑疾病研究与辅助诊断。

功能连接分为无向的功能连接及有向的功能连接。其中，无向的功能连接简称功能连接，强调的是不同脑区之间的相关性。有向的功能连接又称为有效连接(effectiveconnectivity)或因果连接，强调的是两个脑区之间的因果关系。格兰杰因果分析(grangercausalityanalysis，简称granger)和稀疏贝叶斯网络(sparsebayesiannetwork,sbn)是两种有效连接分析方法。granger方法采用的是多元自回归模型，根据两个时间序列的过去值对现在值的预测结果，判断两者之间的granger因果关系(见c.w.j.granger,investigatingcausalrelationsbyeconometricmodelsandcross-spectralmethods,econometrica37(3):424-4381969)。假设xt和yt为两个零均值的平稳时间序列，那么因果模型可以表示为

其中，ap，bp，cp和dp为自回归参数，εt和ηt是两个零均值不相关的白噪声序列，e[εt]＝0，e[ηt]＝0，且

xt和yt使用其过去的m个值进行预测，m也可以理解为自回归模型的阶数。如果使用了时间序列yt的过去值yt-j,j＝1,…,m去预测时间序列xt，比仅仅使用xt的过去值xt-j,j＝1,…,m去预测时间序列xt的噪声方差更小，则可以说yt是xt的granger原因；反之亦然。

贝叶斯网络(bayesiannetwork,bn)是一种基于图论和概率论的数学模型，通过一些变量的信息来推断其它变量的概率信息，适于分析多变量间具有不确定性的复杂关系，已广泛应用于遗传学、生态学、社会科学和脑科学等领域。bn结构是一个有向无环图(directedacyclicgraph,dag)，由代表随机变量的节点及连接这些节点的有向边(由父节点指向其子节点)构成，表示了各变量之间的依赖或独立关系。huang等人提出了一种sbn结构学习算法(huangs,lij,yej,etal.,asparsestructurelearningalgorithmforgaussianbayesiannetworkidentificationfromhigh-dimensionaldata,ieeetransactionsonpatternanalysisandmachineintelligence,2013,35(6):1328-1342.)。其基本做法是，增加了一项l1范数惩罚项以施加连接的稀疏性约束；同时增加了另一项惩罚项以保障bn是一个有向无环图。稀疏约束可以有效地减少bn结构学习的计算量，同时sbn将传统bn结构的两步学习(先确定每个节点可能的父节点，再从中识别最终的父节点)简化为一步学习，即直接确定所有变量的父节点，降低了丢失父节点的可能，提高了结构学习的精度。理论分析和实验测试结果均表明，sbn结构学习算法较之现有十种bn结构学习算法在精度和效率方面均有提升。

对于同一种fmri数据，不同分析方法所检测的脑功能连接应该既有区别又有联系。然而，在现有脑功能连接分析中，各种方法的分析大多彼此独立，导致不同方法的分析结果也相互独立，难于比较；而且，sbn方法尚未用于fmri数据的脑功能连接分析。因此，进行sbn和granger两种方法的联合分析，挖掘二者之间的内在联系，提取鲁棒而稳定的脑功能连接，对深入推进脑功能研究和脑疾病诊断极为重要。

技术实现要素：

本发明的目的在于，提供一种对fmri数据联合进行sbn与granger有效连接分析的方法，揭示两种方法检测结果之间的内在联系，提取稳定的脑网络连接。

本发明的技术方案是，对于从多被试fmri数据中提取的tc成分，同时进行sbn分析与granger因果连接分析。通过对sbn分析方法施加回归系数约束，以及对granger因果连接分析方法施加因果一致性(causalitycoherence)约束，因果一致性也就是因果性强度，寻找两种有效连接方法间的内在一致性，提取两种方法共有的网络连接，作为多被试fmri数据的稳定脑功能连接加以输出。

对于一个有p个节点的bn(图1中的bn有5个节点)，这里用一个p×p矩阵g来表示其结构。gij＝1表示节点xi到xj有一个有向弧，xi是xj的父节点；gij＝0表示节点xi到xj无连接。另外，用一个p×p矩阵p来表示bn的有向路径，如果节点xi到xj有一个有向路径，则pij＝1，反之则为0。如果节点x1是节点x2的父节点，且节点x2是节点x3的父节点，那么p13＝1，而g13＝0。

当每个节点都满足多元正态分布时，节点xi可用下面公式(5)中的回归系数模型来表示：

其中，t(xi)表示节点xi所有父节点的集合；εi为高斯白噪声，均值为0，方差为βi＝(β1i,…,βpi)^t是节点xi的回归系数向量，βij∈(0,1)，j＝1,…,p；b＝[β1,…,βp]表示回归系数矩阵。

sbn方法通过学习回归系数矩阵b，达到确定bn结构的目的。βi非零元素越少，则bn连接越稀疏。huang等人证明，对于一对节点xi和xj，dag的充分必要条件是βij×pji＝0。所以，sbn方法的目标函数如公式(6)所示：

其中，表示节点xj属于集合x/i(此节点集中不包含xi)，即xj是与xi不同的节点；λ1为稀疏约束参数，λ2为dag约束参数。sbn方法的实现流程如图2所示，具体步骤如下：

步骤1，输入包含p个节点数据的数据矩阵x∈r^p×n，停机参数ε，以及约束参数λ1，λ2；λ2>(n-1)²p/λ1-λ1，n为每个节点数据的长度。λ1可根据对稀疏性的要求设定，λ1越大则bn结构越稀疏；ε可根据对结果精度的要求设置，ε越小精度越高；

步骤2，初始化结构矩阵g∈r^p×p，没有先验结构信息时，令g＝0；有先验结构信息时，若节点i到节点j有一个有向弧，则gij＝1，反之gij＝0；初始化回归系数矩阵b＝[β1,…,βp]∈r^p×p，βi＝(β1i,…,βpi)^t，i＝1,…,p，令b＝0；

步骤3，令迭代次数k＝0；

步骤4，对于回归系数矩阵b，设置当前学习列数i＝1；

步骤5，对结构矩阵g进行广度优先搜索，得到有向路径矩阵p：如果节点i到j有一个有向路径，则pij＝1，反之则pij＝0；

步骤6，对回归系数矩阵b第i列βi^(k)的每个元素进行学习，得到βi^(k+1)：设置当前学习元素j＝1；

步骤7，利用公式(7)学习得到βi^(k+1)中的第j个元素

其中x/(i,j)表示不包含第i个节点和第j个节点数据的数据矩阵x，βi/j表示不包含第j个元素的第i列回归系数βi；

步骤8，判断j是否小于p，若是，则j＝j+1，并跳转到步骤7；若否，则跳转到步骤9；

步骤9，应用βi^(k+1)更新g矩阵：如果那么令gji＝1；

步骤10，判断i是否小于p，若是，则i＝i+1，并跳转到步骤5；若否，则跳转到步骤11；

步骤11，判断||b^k+1-b^k||2是否小于ε，若否，则k＝k+1，并跳转到步骤4；若是，则跳转到步骤12；

步骤12，输出最终的结构矩阵g^new、回归系数矩阵b^new。

本发明双约束分析的具体步骤如下：

第一步，输入从k个被试fmri数据中分别提取的tc成分，记为x′^(k)∈r^l×t,k＝1,…,k，l为tc成分的个数，t为tc成分的时间点数；

第二步，采用截止频率为fl～fh(单位hz)的带通滤波器对x′^(k)进行滤波，得到滤波后的时间过程成分x^(k)∈r^l×t，k＝1,…,k；带通滤波器截止频率fl～fh可根据tc成分的频率特性进行选择，对于静息态fmri数据，可选为0.015～0.15(单位hz)；

第三步，对于sbn分析方法，逐个被试计算：令k＝1；

第四步，设置停机参数ε，约束参数λ1，λ2，利用图2所示的sbn方法实现流程，由滤波数据x^(k)∈r^l×t计算第k个被试共l个节点的sbn回归系数矩阵i,j＝1,…,l，i≠j；

第五步，判断k是否小于k，若是，则k＝k+1，并跳转到第四步；若否，则跳转到第六步；

第六步，剔除非显著连接，保留显著连接：对于k个被试的各个相同连接的i,j＝1,…,l，i≠j，例如分别利用单样本t检验计算t值，保留p<0.05即t>tth的显著性连接，剔除其他连接；tth是在自由度为k-1的情况下，p<0.05所对应的t值，可查表得到。例如，k＝82时，tth＝1.99；

第七步，对于显著连接，计算k个被试的各回归系数的平均值，记为

第八步，约束回归系数：挑选出回归系数的连接，剔除的连接；β′可根据对的约束强度进行设置，约束强度越大，β′越大；

第九步，对于granger因果连接分析方法，逐个被试计算：令k＝1；

第十步，计算第k个被试数据x′^(k)中每两个tc成分间的因果一致性和i,j＝1,…,l，i≠j。根据公式(1)(2)(3)(4)，若tc成分i与成分j的时间序列用xt与yt表示，那么：

式中，

第十一步，判断k是否小于k，若是，则k＝k+1，并跳转到第十步；若否，则跳转到第十二步；

第十二步，剔除非显著连接，保留显著连接：对于k个被试的各个相同连接的因果一致性参数和i,j＝1,…,l，i≠j，例如分别利用单样本t检验计算t值，保留p<0.05即t>tth的显著性连接，剔除其他连接；tth的意义同第六步介绍；

第十三步，对于显著连接，计算k个被试的各因果一致性的平均值，记为和

第十四步，约束因果一致性差值：挑选出的连接，剔除的连接；c′可根据对的约束强度进行设置，约束强度越大，c′越大；

第十五步，对于上述sbn分析与granger因果分析得到的连接结果，参见第八步和第十四步结果，检测两种方法的相同连接；

第十六步，输出两种方法的相同连接，作为稳定连接。

本发明所达到的效果和益处是，通过利用双约束联合分析方法，建立sbn与granger两种有效连接分析方法之间的桥梁，提取稳定而鲁棒的脑功能连接，为基于fmri数据的脑功能研究和脑疾病诊断提供更好的技术支持。例如，以默认网络(defaultmodenetwork,dmn)内部7个子成分为对象(见图4-图8)，以本发明提取的稳定连接(见图8)为目标网络，采用granger有效连接分析方法对40个健康被试对照组与42个精神分裂症患者组进行有效连接分析，结果表明，精神分裂症患者组与健康对照组在目标网络的连接结果上存在显著差异，可用于精神分裂症分类。

附图说明

图1是五节点bn有向无环图示例；

图2是sbn方法的实现流程图；

图3是本发明的流程图；

图4是对82个被试静息态fmri数据7个dmn子网络进行sbn分析所得的显著连接结果图；

图5是对sbn的显著连接结果施加回归系数约束后得到的功能连接图；

图6是对82个被试静息态fmri数据7个dmn子网络进行granger分析所得的显著连接结果图；

图7是对granger的显著连接分析结果施加因果一致性差值约束后得到的有效连接图；

图8是sbn和granger联合分析输出的稳定连接。

具体实施方式

下面结合技术方案和附图，详细叙述本发明的一个具体实施例。

现有82个被试(k＝82)的静息态fmri数据，利用gift工具箱(http://mialab.mrn.org/software/gift/index.html)提供的groupica方法(v.calhoun,t.adali,g.pearlson,andj.pekar,amethodformakinggroupinferencesfromfunctionalmridatausingindependentcomponentanalysis,humanbrainmapping14:140-151,2001)，采用infomax算法分离得到120个脑空间激活区成分及其对应的tc成分。根据smith脑空间激活区参考模板(s.m.smith,p.t.fox,k.l.miller,d.c.glahn,p.m.fox,c.e.mackay,etal.,correspondenceofthebrain’sfunctionalarchitectureduringactivationandrest,proceedingsofthenationalacademyofsciences106(31):13040-13045,2009)，按照相关最大化原则挑选出默认网络(defaultmodenetwork,dmn)的l＝7个感兴趣子网络，提取其对应于每个被试的7个tc成分，每个tc成分的时间点数t＝150，没有先验结构信息。实施本发明的步骤如图3所示，具体如下：

第一步，输入从82个被试fmri数据中分别提取的7个tc成分，记为x′^(k)∈r^7×150,k＝1,…,82；

第二步，采用截止频率为0.015～0.15(单位hz)的带通滤波器对x′^(k)进行滤波，得到滤波后的时间过程成分x^(k)∈r^7×150，k＝1,…,82；

第三步，对于sbn分析方法，逐个被试计算：令k＝1；

第四步，设置停机参数ε＝0.01，约束参数λ1＝50，λ2＝10[(t-1)²l/λ1-λ1]＝30581.4，利用图2所示的sbn方法实现流程，由滤波数据x^(k)∈r^7×150计算第k个被试共7个节点的sbn回归系数矩阵i,j＝1,…,7，i≠j；

第五步，判断k是否小于82，若是，则k＝k+1，并跳转到第四步；若否，则跳转到第六步；

第六步，剔除非显著连接，保留显著连接：对于82个被试的各个相同连接的i,j＝1,…,7，i≠j，分别利用单样本t检验计算t值，保留p<0.05即t>1.99的显著性连接，剔除其他连接，结果如图4所示；

第七步，对于显著连接，计算82个被试的各回归系数的平均值，记为

第八步，约束回归系数：挑选出回归系数的连接，剔除的连接，结果如图5所示；

第九步，对于granger因果连接分析方法，逐个被试计算：令k＝1；

第十步，利用公式(8)(9)，计算第k个被试数据x′^(k)中每两个tc成分间的因果一致性和i,j＝1,…,7，i≠j；

第十一步，判断k是否小于82，若是，则k＝k+1，并跳转到第十步；若否，则跳转到第十二步；

第十二步，剔除非显著连接，保留显著连接：对于82个被试的各个相同连接的因果一致性参数和i,j＝1,…,7，i≠j，分别利用单样本t检验计算t值，保留p<0.05即t>1.99的显著性连接，剔除其他连接，结果如图6所示；

第十三步，对于显著连接，计算82个被试的各因果一致性的平均值，记为和

第十四步，约束因果一致性差值：挑选出的连接，剔除的连接，结果如图7所示；

第十五步，对于上述sbn分析与granger因果分析得到的连接结果，参见图5和图7的结果，检测两种方法的相同连接；

第十六步，输出两种方法的相同连接，见图8，作为稳定连接。