基于组MCP和复合MCP的人脑功能超网络分析及抑郁症分类研究

薛晓倩，李 瑶，梁家瑞， ，孙 超，郭 浩

1(太原理工大学 信息与计算机学院，山西 晋中 030600)2(太原理工大学 软件学院，山西 晋中 030600)

1 引 言

研究发现，大脑的自发神经元活动和低频血氧水平依赖(Blood Oxygen Level-Dependent，BOLD)信号有密切的联系[1]，这表明在静息状态下，BOLD信号构建的功能连接网络能够反映大脑区域之间的交互活动.在过去的几年中，静息态功能磁共振成像技术不断发展，脑功能网络作为大脑交互的简化表示已被广泛地应用在抑郁症[2]、阿尔茨海默症[3]等精神疾病的研究中.

在传统的功能连接方法中，往往是图形化模型[4]和基于相关性分析方法[5]等占据主要位置，然而它们都存在不足之处：图形化模型依赖于图的先验知识[6]；基于相关性分析方法则存在两个主要问题：1)其构建的连接不保证绝对可靠，可根据自定义阈值选择连接，因此存在虚假连接[7]；2)由于其连接网络是两两关系的映射，因此不能解决大脑信息传递的高阶互动过程[6].

相关的研究证明，大脑的信息传递在局部场电位、神经元同位素示踪以及皮层活动中均存在明显的高阶交互[8].因此，超网络被提出[9].脑区对应网络中的节点，如果一个脑区与其他脑区之间有交互，则认为这些脑区中存在超边.

传统的超网络通过LASSO方法进行构建，然而该方法存在两个问题：1)LASSO中的惩罚函数对系数的过强压缩，导致模型中目标变量回归系数并不是无偏估计或近似无偏估计；2)LASSO是单一变量选择模型，缺乏组效应的解释能力，这使得所构建的脑功能超网络将有可能失去一些重要的具有明显神经学意义解释的连接.目前对后者的解决方法主要是通过引入分组的方法模型，如组LASSO[10]、 稀疏组LASSO[11]以及弹性网[12，13]方法.然而组LASSO和弹性网均是基于LASSO方法的延伸，而稀疏组LASSO则是组LASSO方法的延伸，因此上述方法同样存在惩罚函数对系数的过强压缩的问题，导致构建的超网络也许过于严格，从而失去一些重要的连接.

因此，本文提出两种基于极小极大凹惩罚(The Minimax Concave Penalty，MCP)的无偏稀疏模型：组MCP[14]和复合MCP[15，16]，通过解决稀疏线性回归模型进行超网络构建.MCP具有变量选择连续性和无偏性，不会出现惩罚函数对系数的过强压缩的问题.两种方法均能解决组效应问题，不同的是，组MCP方法将MCP惩罚仅应用到组级，解决了组效应问题；而复合MCP方法则把MCP惩罚既用作组级惩罚函数，又作为组内惩罚函数，分别在组级和组间对变量进行选择.两种方法对系数的压缩程度不同，构建出来的超网络也不同.

研究结果表明，两种方法均优于传统方法，而两种方法所构建的脑功能超网络结构有较大差异，考虑到是两种惩罚函数对回归系数的压缩方式导致的，复合MCP方法构建的超网络的超边数相较于组MCP方法分布集中，而组MCP方法的超边分布较为分散；选取超网络组间差异拓扑指标为分类特征，利用支持向量机(Support Vector Machine，SVM)方法构建分类模型.复合MCP具有最好的分类效果和较高的分类权重；最后，进一步分析了方法模型中的各个参数，探究方法参数和分类器参数以及显著性脑区对方法是否产生影响.本文提出的方法一定程度上解决了无偏稀疏脑功能超网络构建过程中的组效应问题，所构建的脑功能超网络可以更好地表达抑郁症患者与正常对照的结构差异，具有重要的理论意义和临床价值.

2 材料和方法

2.1 数据采集和预处理

本文严格按照山西医学伦理委员会的要求(会议号：2012013)，在该项研究实施之前，与每位参与者均达成了书面协议.根据赫尔辛基宣言，书面知情同意书由实验中的每一个被试签署.总共招募了66名被试，其中包括38名首发，无用药重度抑郁症患者(Major Depression Disorder，MDD)(15名男性)和28名健康右利手志愿者(13名男性).静息状态下，应用3T磁共振扫描仪(Siemens Trio 3-Tesla scanner，Siemens，Erlangen，Germany)对他们进行功能磁共振成像(fMRI)扫描.被试的具体情况可参照表1.

表1 被试的具体信息Table 1 Specific information of participants

山西医科大学第一医院完成数据采集任务，并且由精通磁共振技术的放射科医师完成扫描任务.被试在进行扫描时也有相应的要求，即被试需要在放松的状态下闭上眼睛，但不能入睡，需要保持清醒，也不进行特定的思考.每个扫描的结果是248个连续的EPI功能图像，其中的扫描参数有如下设置：33 axial slices，echo time(TE)= 30ms，repetition time(TR)=2000ms，thickness/skip = 4/0 mm，flip angle = 90°，matrix = 64×64 mm，field of view(FOV)=192×192 mm.前10个功能图像的时间序列由于被试对环境的自适应性以及初始磁共振信号的不稳定性而被丢弃.

表1中的a值表示通过双值双尾t检验得到P值；b值表示通过双尾Pearson卡方检验得到P值.用SPM8[注]http：//www.fil.ion.ucl.ac.uk/spm.来完成数据的预处理过程.首先，头动校正和时间片校正必不可少，而在校正过程中，抑郁组和对照组中分别出现2例转动大于3度或者头动大于3毫米的被试，因此弃除这些被试的扫描数据.需要注意的是，最后的66例样本数据中不包含那些被丢弃的数据.接着，图像经过优化仿射变换会被标准化到MNI(Montreal Neurological Institute，MNI)的标准空间中.最后，为降低生物高频噪音和低频漂移的影响，对数据进行带通滤波(0.01-0.10Hz)和线性降维.

脑区分割时使用AAL[17](Anatomical Automatic Labeling，AAL)模板，得到90个感兴趣解剖区域(Region Of Interest，ROI)，它们均等分布在两个脑半球中，这一模板已广泛应用在同类研究中[18].脑网络中的一个节点可以由一个感兴趣 解剖区域表示，因此90个感兴趣解剖区域可以看成90个节点分布在脑网络当中.首先，按照超图的基本定义[19]，提取各脑区的时间序列信号，脑区之间的关联关系可通过计算获得；然后，由于头动校正和白质信号以及平均脑脊髓液会对信号产生影响，因此需要排除它们的影响，进行脑区时间序列的回归分析；最后，依据脑区的时间序列，在两种方法模型下构建脑功能超网络.

2.2 超网络构建

变量选择一直是构建模型的热门方法，自LASSO方法被提出以来，各种变量选择方法层出不穷.在对LASSO模型分析研究后发现，最初人们只是着眼于其仅仅作为一种单变量选择方法，无法解决组效应问题，因而提出各种组级选择方法.事实上，LASSO模型不仅存在以上问题，其惩罚函数对系数还具有过强压缩的特性，从而导致模型中目标变量回归系数的有偏估计.而随后的诸多方法则是在对LASSO方法缺乏组效应解释能力情况下的改进，同样存在系数的过强压缩问题.

因此，本文提出使用MCP惩罚函数，该函数具有LASSO不具有的良好性质，将其应用到组级，既不会对系数过度压缩，又能解决组效应问题.基于MCP惩罚，提出两种变量选择方法构建超网络，分别是组MCP方法和复合MCP方法.

2.2.1 组MCP方法构建超网络

组MCP是将MCP惩罚应用到组级上形成的.MCP是一种单一变量选择方法，具有同步选择一致性和渐近无偏性[20].这种性质意味着该模型渐近等效于最大似然模型的拟合，在该模型中，事先已知真正非零系数的恒等式[21].将MCP扩展到组级，就得到了组MCP，能在组级上对变量进行选择.

组MCP定义如下[14]：

(1)

(2)

构建方法如下：首先，使用协调下降算法实现组MCP的稀疏线性回归；然后，将每一个ROI表示为一个节点，通过变化λ的值来构建超边，根据αk中的非零项构建超边，通过在指定范围内变化λ值可以产生一定数目的超边；最后，将不同λ值构建出来的网络合并为一个，最终的超网络是一个具有90个节点，810条超边的90×810的矩阵.γ采用默认值3.

2.2.2 复合MCP方法构建超网络

组MCP方法只能反映组级选择，未体现组间选择，因此引入复合MCP[15]方法，它也是近似无偏稀疏模型.通过将MCP惩罚同时用作内部惩罚和外部惩罚，可以实现双级选择，既能在组间选择变量，又能选择组内的重要变量.复合MCP方法的这种惩罚方式意味着它既可以在组间产生稀疏度，又能在组内保持稀疏性，因而可以同时选择组间变量和该组内的变量.

复合MCP方法定义如下[16]：

(3)

式中各参数与上述公式(1)中的相应参数有相同的含义.φλ，η(·)为MCP惩罚：

(4)

复合MCP方法首先把MCP惩罚应用于各个组间变量选择上，然后在组级加入MCP惩罚，变量是否进入模型取决于模型自身以及它所属于的组.因此，该变量也许进入模型，也可能不进入模型而被消除.实验过程与组MCP相似，不同之处在于进行稀疏线性回归模型求解时的惩罚函数为复合MCP，通过在一定范围内变化λ的值可以构建一定数目的超边.γ1和γ2均采用默认值3.

2.3 特征提取

构建功能连接超网络之后，选择超网络中每个顶点的属性值作为特征.对已构建的超网络进行指标计算，聚类系数如今已广泛用于衡量网络的局部属性，这3类指标统称为HCC指标，该指标下有3种聚类系数，分别具有不同的定义[23]，分别记作HCC1、HCC2以及HCC3.

(5)

(6)

(7)

公式(5)计算与节点v不存在连接的相邻节点的数目.节点用u，t，v表示，N(v)是集合，指超边中除了节点v以外，还包括其它节点的集合.如果E表示边集，ei表示某一条超边，则I(u，t，v)=1当且仅当∃ei∈E，u，t∈ei，但v∉ei；否则I(u，t，v)=0.HCC1查找不包含u的邻居之间的连接，HCC1的优点是，在这个集合中发现的任何交互都可能表示邻居之间的真实连接.它的缺点在于可能过分关注于那些次要的共享连接，这些连接与u的交互没有什么关系.

公式(6)计算与节点v存在连接关系的相邻节点的数目，式中各相同参数的含义同上.如果∃ei∈E且u，t，v∈ei，则I′(u，t，v)=1.HCC2查找那些包含u的邻居的连接，这种方式找到的边真实地反应了u和邻居之间的聚集.但是需要注意这种连接可能只是与u共享连接的人工数据.

公式(7)计算超边之间的重叠量，具体指节点v的相邻超边，式中各相同参数的含义同上.|e|表示超边中包含的节点数目；S(v)指超边集合，这些超边中都含有节点v，且S(v)={ei∈E：v∈ei}.通过邻域超边的重叠量来衡量邻域的密度.其分子表示与u相关联的超边的顶点个数.分母表示这种重叠可能的数目.HCC3则通过邻域超边的重叠量来衡量邻域的密度.与以上两种定义都不同，它从节点的角度来定义重叠量.

超网络的局部聚类属性可由HCC1、HCC2以及HCC3 3种属性从不同的角度反映得出.在两种方法构建的超网络中，将每一种聚类系数作为特征，分别提取.每个被试都有90个脑区，为了简便，计算出被试的平均聚类系数(对90个脑区进行平均，分别记为平均HCC1，平均HCC2以及平均HCC3).将这3种不同定义下的聚类系数应用于两种方法所构建的超网络，最后分别把3种聚类系数进行融合.

2.4 特征选择和分类

本文利用统计学方法进行特征选择.在样本量较小的情况下，使用Kolmogorov-Smirnov(KS)[24]非参数置换检验来分析两组数据是否具有差异是一种常用的方法.在本文中，对抑郁症患者和正常人的超网络在计算HCC聚类系数后获得的270个节点属性进行组间KS非参数检验，并进行FDR校验(q=0.05)，通过非参数检验分别得出每个被试的270个脑区的p值大小.选取抑郁症患者和正常人的p<0.05的显著性特征，并将这些显著性差异特征作为分类特征，进行分类.

分类性能通过留一交叉验证进行评估.训练集使用K折交叉验证(K-fold Cross Validation，K-CV)对参数c和g进行寻优[25].在训练集验证中，c和g的最优值取分类准确率最高的那组值，并由此建立N个不同的模型.在标准化分类特征的平均值和标准偏差后，最后的分类结果取不同模型的分类准确率的平均值.

3 结果和分析

3.1 超网络结构分析

为确定显著差异是否存在于不同方法构建的超网络中，对两种方法下的超网络结构进行了比较，并进行了如下的实验：

首先分别分析两种方法的正常人和抑郁症患者的超网络的超边，计算它们的超边的边度.结果表明，两种方法下的超网络结构存在差异，具体说来，复合MCP方法构建的超网络的超边度大多数分布在2-13范围内；组MCP方法构建的关于正常人和抑郁症患者的超网络的超边度分布范围较广，分布比较离散，超边度分布情况如图1所示.

在每一个被试的超网络中，计算出它们(对90个脑区进行平均)的平均聚类系数，分别记作平均HCC1，平均HCC2和平均HCC3.对抑郁组和正常对照组进行非参数置换检验，对于平均HCC1，平均HCC2以及平均HCC3，分别比较两种方法构建的超网络的差异.在两种方法所构建的超网络中，平均聚类系数的分布情况如图2所示.

图1 抑郁组和正常对照组的超边度分布图Fig.1 Hyper-edge distribution of depression group and normal control group

图2 抑郁组和正常对照组的平均聚类系数Fig.2 Average clustering coefficient of depression group and normal control group

结果表明，无论是正常对照组还是抑郁组，两种方法下的平均聚类系数关于HCC2和HCC3分别具有显著性差异，而关于HCC1的显著性差异只存在于正常对照组中，抑郁组中则未发现显著性差异(p>0.05，FDR校验，q=0.05).因此，两种方法构建的超网络存在结构上的差异.分析其潜在的原因，组MCP方法是在组级上进行变量选择，而复合MCP是双级选择，既在组级上进行变量选择，又选择组内重要变量，从而造成了超网络结构上的不同.

3.2 差异性脑区

在两种不同的方法下分别进行超网络构建并提取特征，对于每一个特征，进行非参数置换检验对所有被试评估抑郁症患者和正常人之间的差异.表2和表3分别列出了两种不同超网络构建方法下所得到的具有显著差异的大脑区域.

表2 组MCP方法下具有显著差异的大脑区域Table 2 Brain regions with significant differences under the group MCP method

表3 复合MCP方法下具有显著差异的大脑区域Table 3 Brain regions with significant differences under the composite MCP method

表中黑体字表示对应的聚类系数的P值<0.05.两种方法利用统计分析获得的差异性脑区也不相同，其中不重叠的脑区数目较多，主要包括左侧中央前回、右侧中央前回、右侧额中回、右侧补充运动区、左侧内侧和旁扣带脑回、右侧杏仁核、左侧枕上回、左侧顶上回、左侧顶下缘角回、右侧中央旁小叶、左侧豆状壳核和左侧颞极：颞上回.在后面的研究中，复合MCP方法得出了最高的分类准确率，因此本文对该方法下的差异性脑区进行讨论，且验证了复合MCP方法下的很多脑区都与已有文献一致，证明它们和抑郁症病理有着密切联系，具体见表4.

表4 已有文献中关于抑郁症的脑区Table 4 Brain areas of depression in the literature

3.3 分类表现

为评估两种超网络构建方法的分类性能，本文使用了分类准确率、敏感度以及特异度3个特定指标.除此之外，由于不平衡数据集所带来的膨胀性能，将特异度和敏感度的算术平均值定义为BAC.

同时，将文中所提出的方法与广泛应用至脑疾病中的LASSO方法和组LASSO以及稀疏组LASSO方法进行了对比.同样也使用分类准确性、敏感度以及特异度3个指标.结果证明本文所提出的两种方法的分类准确度均优于已有方法.比较结果如表5所示.

表5 不同方法的分类表现Table 5 Classification performance of different methods

组MCP和复合MCP方法的分类准确率都能达到85%以上，其中复合MCP方法的分类准确率最高，能达到89.39%.两种方法对变量是否进入模型采取了不同的解决办 法，组MCP方法将惩罚函数运用在组级层面，解决了组效应问题；而复合MCP方法则在组级和组内都应用了惩罚函数，既选择了组间重要变量，又选择出组内重要的变量，因而是比组MCP更为严格的方法模型.以上结果表明在脑功能超网络构建中，组级选择方法能更好地体现人脑中固有的组结构，且双级选择能较好地体现人脑中复杂的层次关系.

图3 两种方法的分类权重Fig.3 Classification weight of the two methods

此外，使用Relief算法[37]对特征的有效性进行验证.该方法会根据每个类别和功能相关性分配不同的权重，特征的分类能力越强，分配的权重越大(反之亦然).本文对两种方法以及传统LASSO方法对应的分类权重进行了比较，结果表明，两种方法的分类权重均高于LASSO方法，复合MCP方法具有最高的分类权重，结果显示在图3中.结果同样表明基于双级选择的复合MCP方法能构建合理有效的超网络.

4 参数讨论

4.1 分组参数k的讨论

两种方法都在预定义分组下进行超网络构建，实验中使用了k-means[38]聚类算法，k表示分组大小，可以将数据预先分成k个组，是预定义分组的重要基础.本实验以3为步长，进行组数不大于最大分组90的分组过程，分别在每一个k值下进行50次实验，选取正确率的算术平均值作为最后的分类结果.结果显示当组数k=48时，复合MCP方法表现出最高的正确率89.394%.当k=24时，组MCP方法表现出最高正确率达到86.364%.图4为两种方法在不同的k值下的正确率.

图4 两种方法的不同k值下的准确率Fig.4 Accuracy of the two methods under different k values

4.2 SVM分类参数c和g的影响

SVM分类器在各领域中的应用都比较广泛，而分类时则经常涉及到核函数的选取问题.由于RBF核函数应用广泛，无论是小样本还是大样本均适用，因此，在分类中选取RBF核函数.SVM模型中有两个参数对分类影响较大，即惩罚因子c和核参数g.c表示调节优化方向中两个指标(间隔大小，分类准确度)偏好的权重，即对误差的宽容度，c过大或过小，其泛化能力都会变差；g隐含地决定了数据映射到新的特征空间后的分布，g越大，支持向量的个数越少，反之则越多，而支持向量的个数影响训练与预测的速度.最优的c和g能使SVM的分类性能达到最佳.对于如何找出最佳的c和g，首先选取某一组给定的c和g的值，利用K-CV方法得到在该组c和g下的分类准确率，然后不断更换c，g的值，最终取分类准 确率最高的那组c和g作为最佳参数.c，g两参数可在[2-8，28]内变化，并将1设置为步长大小.结果显示，分类准确率最高时的c，g参数的值分别为2和0.5，最高准确率是92.308%，如图5所示.

图5 c和g的参数寻优结果Fig.5 Parameter optimization results of c and g

5 结 论

如今静息态功能磁共振成像(R-fMRI)技术已被广泛应用于预测抑郁症等精神疾病中.超网络因能充分反应各个脑区之间的互动而被广泛应用.现有的超网络构建多是基于稀疏表示方法，传统的基于LASSO的稀疏表示方法缺乏解释脑区之间的组效应的能力，而且惩罚函数对系数的过强压缩，导致模型中目标变量回归系数的有偏估计.考虑到这些问题，本文提出两种超网络构建方法：组MCP方法和复合MCP方法，前者只能在组级上进行变量选择；后者既能在组级上进行变量选择，又能选择组间的重要变量，称为双级选择.

实验结果表明，本文所提出的两种方法均优于LASSO方法.在超网络的结构上，两种方法之间存在差异，复合MCP构建的超网络的超边分布更为集中，能更好地体现人脑中复杂的层次关系；在分类性能上，本文所提出的两种方法均优于传统LASSO方法以及组LASSO和稀疏组LASSO方法，并能够获得较高的分类权重.

尽管在目前的研究中，实验已取得较好的结果，但仍存在局限性.首先，在实验中假设组间没有重叠的信息交互，但是在很多情况下，重叠组是存在的；其次，两种方法由于聚类初始种子点的随机选取以及聚类数k的不同会造成网络结构以及分类结果的不唯一，建立更加稳定的超边被期望以进一步改善超网络.