近年来, 图模型^[1]成为人工智能、统计机器学习和生物信息等领域的研究热点。图模型是一类用图来表示随机变量联合概率分布的模型, 它被广泛应用于建立各种相互作用的单元之间形成的网络结构, 如基因调控网络、蛋白质网络、社交网络等。图模型的构成要素为节点和边, 其中节点与随机变量对应, 边与变量之间的条件依赖关系对应。具体的, 用二元组G=(V, E)表示一个图, 其中V代表节点, E代表连接节点的边的集合。根据边是否有方向, 图模型分为无向图和有向图。本文聚焦于从数据中估计网络的结构, 利用n个服从的样本估计精度矩阵。为了构造可解释的图, 众多学者提出将正则化思想应用到图模型的模型选择和参数估计问题中。Meinshausen等^[2]提出了邻域选择方法, 将图模型的正则化问题转化为求解p个L₁正则化^[3]问题, 证明了对于稀疏的高维图模型, 邻域选择的变量选择具有一致性。Yuan等^[4]提出了一种罚似然函数方法, 称为Graphical Lasso, 使得精度矩阵转化为正定矩阵, 可以同时进行模型选择和参数估计。Friedman等^[5]提出了分块坐标下降法来求解Graphical Lasso。

模型加入L₁罚项从贝叶斯角度来讲等价于假设每条边所对应参数分布相同, 且独立于其他所有边, 说明图中大多数节点边数大致相等。这在真实网络中是不现实的。Barabási等^[6]研究发现大多数真实网络为无标度网络, 节点度分布服从幂律分布, 大部分节点度很小, 存在度非常大的少数节点, 这种度非常大的节点称为中心点(Hub)。因此, 研究有中心点的网络结构学习有实际意义。

本文聚焦于具有Hub网络的结构学习问题。在邻域选择框架下, 为了使所得网络具有Hub, 本文采用Elastic net^[7]罚项。Elastic net罚项兼有Lasso回归和岭回归的优点, 既能达到变量选择的目的, 又具有很好的群组效应, 从而使所得网络更容易产生Hub。对于所得正则化模型, 采用坐标下降法求解模型。最后, 模拟数据实验和实际数据实验均说明所提方法有效、实用。

1 基于Elastic net的邻域选择模型 1.1 邻域选择方法(neighborhood selection)

高斯图模型^[8]是基于高斯分布假设的无向图模型。假设X=(X₁, …, X_p)~N_pμ, Σ, 其中μ为均值, Σ为协方差矩阵。Θ=(θ_ij)_{1≤i, j≤p}=Σ^-1表示协方差矩阵的逆矩阵, 称为精度矩阵, 则由马尔科夫性质, θ_ij=0表示变量X_i与X_j条件独立, 同时对应节点i和节点j之间无边相连。

X=(X₁, …, X_p)~N_pμ, Σ, 则由多元正态分布的性质得给定除X_i以外的变量X_-i, X_i的条件分布为

式中, θ_i=Σ_{-i, -i}^-1Σ_{-i, i}∈R^p-1, σ_i²=Σ_{i, i}-Σ_{-i, i}(Σ_{-i, -i})^-1Σ_{-i, i}, 那么自然地可将X_i表示为X_-i的线性函数, 即

(1)

式中, ε_i~N(0, σ_i²)。令Θ=Σ^-1, 则由分块矩阵求逆公式得

因此, (1)式中回归系数向量θ_i中的非零元素与精度矩阵中对应列中的非零元素一致。假设有服从N(0, Σ)的n个样本, 即可估计网络的结构。经典的估计方法如最大似然估计和最小二乘估计所得参数不是稀疏的, 因而往往对应完全连接的图。这样的网络并不能有效研究变量之间的相互关系, 且仅限于样本个数n大于维数p的情形。当p远大于n时, 通过引入稀疏性先验来降低参数的自由度, 使得参数可估计，而且使所得网络是稀疏的，更具有可解释性。关于高维高斯图模型, Meishausen和Bühlmann基于(1)式提出了邻域选择方法, 即用L₁正则化方法估计(1)式中参数, 得到每个节点的邻居节点, 如下

式中, 为了方便, 本文后面X_i表示的都是第i个变量的n个样本组成的向量, θ_i是变量X_-i的系数。但是这样得到的精度矩阵不具有对称性, 为此, 在p个L₁回归之后, 需要将最终结果对称化。具体的, 当θ_ij≠0且θ_ji≠0时, 认为节点i和节点j有边连接, 这称为and-rule; 当θ_ij≠0或θ_ji≠0时, 认为节点i和节点j有边连接, 这称为or-rule。

1.2 Hub网络的估计

如前言所述, Hub网络的特点是有少数节点度很大, 因此, Hub网络对应的精度矩阵具有潜在组结构, 每一行或每一列为一组。而基于L₁的邻域选择方法并不具有组效应。为了使所得网络具有Hub, 在基于L₁的邻域选择框架下, 又引入Ridge罚项, 特别地, 模型如下

(2)

式中，θ_i是变量X_-i的p-1维系数向量, 其中每个分量是否为0决定其余节点与当前节点是否相连。λ₁和λ₂是调控参数。值得注意的是, (2)式中的正则项即为Elastic net^[7]正则项。Ridge罚项的引入使得每一个节点对应的边自然地成为一组, L₁罚项的引入使得组内具有稀疏性。且并没有提前假设哪个节点是Hub, 由(2)式可以自动估计出网络的结构, 从而度较大的几个节点即为网络的Hub。对于(2)式的求解, 本文采用坐标下降法^[9]。坐标下降算法的基本思想是依照某种固定顺序, 在每步迭代中, 保持其余系数不变, 只更新一个坐标系数, 更新完所有坐标后进入下一步迭代, 直至收敛。对于(2)式, 每步更新均有显式解, 故算法是高效的。

2 实验

本节分别用真实网络数据和生成的2种网络比较Elastic net(and-rule)、Elastic net(or-rule)、Graphical Lasso、Neighborhood Selection(and-rule)和Neighborhood Selection(or-rule)模型在参数估计、模型选择方面的效果, 并说明调控参数对模型的影响。

2.1 模拟实验

为了评估模型的效果，考虑图模型的邻接矩阵A的2种生成机制。第一种生成带Hub点的网络。先随机选取Hub节点, 令邻接矩阵中相关的行和列的元素等于1的概率是0.8, 即每个节点与Hub点连接的概率是0.8。接下来令A_ij=A_ji=1(i < j)的概率是0.02, 即2个节点有边的概率是0.02, A_ij与A_ji为邻接矩阵A的元素, A_ij=A_ji=1表示节点i和节点j有边相连。A_ij=A_ji=0表示节点i和节点j无边相连。第二种生成无标度网络。一个节点的度是k的概率服从幂律分布P(k)~k^-α, 根据Barabási等^[10]提出的机制:增长性和优先连接性, 直接生成邻接矩阵。对于生成的每个邻接矩阵A, 为了生成精度矩阵Θ, 先构建一个对称矩阵C, 当A_ij=0时, C_ij=0, 即节点i和节点j不相连, 当A_ij=1时, C_ij的取值来自均匀分布U^{[1, 2]}, 也可以取其他值。最后, 令精度矩阵Θ=C+[0.1-λ_min(C)]I_p, λ_min(C)是C的最小特征值, I_p是p×p的单位矩阵以确保Θ的特征值是正数。

按照第一种机制生成3组带Hub点的网络, 样本个数n都是200, 变量p的个数分别是100, 250, 500, 对应Hub点的个数分别是3, 5, 10。按照第二种机制生成3组无标度网络。样本个数n都是200, 变量p的个数分别是100, 250, 500, 对应Hub点的个数分别是3, 5, 10。

比较生成网络和用模型估计得到网络的一些指标衡量模型的效果。第一类带Hub点的网络实验结果如图 1所示。由图 1的结果可以看出Elastic net正则化模型对于带Hub点的网络的估计效果比其他3种模型都要好。第二类无标度网络实验结果如图 2所示, 由结果可以看出Elastic net正则化模型对Hub点的估计效果比其他模型都好。

2.2 真实数据实验

本节将本文方法应用到肺癌的基因数据集, 以检验模型的效果。这些原始数据可在NCBI(national center for biotechnology information)中下载。该数据集包含了7 129个基因表达, 这些基因是由肺癌患者和普通样本组成的96个样本，本文随机选取79个基因。研究表明Hub基因在肺癌基因调控网络中起着重要作用, 特别是它们可能是肺癌进展的潜在生物标志物。因此, 分析每个方法所发现的Hub是有意义的。让调控参数λ在一定范围内变动并拟合模型, 记录在该范围内每个方法对应的每个基因度的总和。将提出的Elastic net正则化模型应用到肺癌的基因数据。选择度排序为前15个基因作为每种方法选择的潜在的Hub基因, 结果如表 1所示。G Lasso表示Graphical Lasso, NeiSel(and)表示Neighborhood Selection(and-rule), NeiSel(or)表示Neighborhood Selection(or-rule), El net(and)表示Elastic net (and-rule), El net(or)表示Elastic net (or-rule)。

通过查阅文献发现, 其中大部分基因已被证实与肺癌相关, 用黑色实心圆圈表示这些基因。从表 1可以看出El net(and)和El net(or)选出的15个Hub基因中有12个已经在相关文献中被证实, 多于其余所有方法。从而说明在该数据集上, Elastic net正则化模型表现更好。

2.3 调控参数实验

Elastic net的罚函数为岭回归罚函数和Lasso罚函数的凸线性组合, 即α|β|₁+(1-α)|β|², 0≤α≤1。α=0时, Elastic net即为岭回归; 当α=1时, Elastic net即为Lasso回归。α表示的是L₁范数惩罚项所占比例, 实验中通过改变α的值控制调控参数, 模拟实验中α的取值为0.001。为了说明调控参数对模型的影响, 改变模拟实验中α的值, 生成n=200, p=100, Hub点的个数为3的带Hub点的网络, α分别取0.1, 0.01和0.000 5。结果如图 3所示。可以看出α=0.1时, 模型的效果差异并不很明显, α=0.01时, 只有估计Hub的边数这一组的效果有明显差异, α=0.000 5时, 模型的效果有明显差异, Elastic net正则化模型对Hub点的估计效果比其他模型都好。

3 结论

本文针对带有中心点的高斯图模型, 提出新的正则化模型, 将邻域选择模型与Elastic net罚项相结合, 构造了新的正则化模型。本文使用坐标下降法求解模型。最后由模拟实验以及真实数据表明所提出的新模型对于带Hub点的高斯图模型结构的估计有较好的效果。