图模型：从随机变量到社区发现

对于PM2.5数据，我们仅有的信息是各个城市每个小时的PM2.5浓度值，为了研究城市之间PM2.5污染形成的网络结构，一种可能的思路是构造PM2.5相关性矩阵进行研究。但是，众所周知，相关性得出的两两城市之间的PM2.5相互关系并没有排除其余城市的影响，所得出的网络结构缺乏真实性。因此，更合理的做法是排除其他城市影响，研究两两城市之间的条件相关关系。图模型[14]为研究变量之间的条件相关关系提供了有效工具，它被广泛应用于建立各种相互作用的单元之间形成的网络结构，如基因调控网络、蛋白质网络、社交网络等，从而为研究和理解变量之间的相互关系提供了工具和参考。具体地，图模型是一种将图和随机变量的联合概率分布通过马尔科夫性质联系起来的框架，其中节点与随机变量对应，边与随机变量的条件依赖关系对应。具体地，用二元组G=（V，E）表示一个图，其中V=1，…，p表示节点集合，E⊆V×V表示边的集合。本章中，我们将交替使用图和网络来描述变量之间的连接关系。以Gaussian图模型为例，对于Gaussian图模型，即随机变量的联合分布为均值为0、协方差矩阵为Σ的p维多元Gaussian分布，每个变量对应网络中一个节点，变量之间的边即变量之间的条件独立关系由协方差矩阵的逆Θ=（Σ）-1=（θij）的非零元素决定。具体地，节点i和节点j之间有边等价于节点i和节点j不条件独立，同时等价于θij≠0。因此，网络结构与协方差矩阵的逆相对应，我们可以通过估计协方差矩阵的逆来建立城市之间的PM2.5污染网络。

中心点（hub）和社区（community）是真实世界中网络的显著特点，因此对数据进行网络建模时，需要考虑这些特点或先验。中心点和社区缺乏严格定义，我们给出如下描述性定义：中心点指网络中度较大的点，即相比于其他节点，和中心点相连接的边更多，这类节点在网络中起到至关重要的作用，如在社交网络中，中心点往往是具有影响力的公众人物；在生物网络中，中心点往往是对某种功能起到关键作用的基因，甚至可能是治疗某种疾病的潜在靶点。社区指网络中连接紧密的子图，同一社区中的节点和节点连接很紧密，而社区和社区之间的边连接很稀疏。网络社区往往由具有相同特征或属性的节点构成，如社交网络中，社区可能由具有相似兴趣爱好的人构成；生物网络中，社区可能由具有相似功能的基因构成。对于PM2.5污染网络，网络的中心点为对污染影响力大的城市，或者说为与其余城市污染相互影响较大的城市，对于这类城市，需要着重关注和治理；而网络的社区包含了污染具有相似特点的城市，这些城市或许有着相似的经济结构、污染机制，对于不同社区的城市，需要分社区制定治理规则，而对于同一社区的城市污染治理需要相互协调。因此，有必要深入探究PM2.5污染网络的中心点和社区，从而为治理PM2.5污染提供数据支撑和依据。下面我们分别具体介绍实现这两个目标所用到的模型和方法。

1.无标度图模型

Barabasi和Albert[15]研究发现，不同于随机网络，基因网络、神经网络、社会网络等大多数真实网络为无标度网络，即网络中节点的度服从幂律分布。无标度网络的典型特征是网络中大部分节点的度很小，而存在少数度很大的节点，称之为中心点。现有图模型方法大多数聚焦于稀疏网络，即所估计出的网络是稀疏的，但并不具有无标度特点[16-18]。基于网络无标度的特点，Guo，Zhang和Wu[19]提出了无标度图模型，该模型通过构造罚函数充分将网络的无标度先验结合到了模型中，下面我们进行具体介绍。

假设样本矩阵为其中每个样本Xi独立同分布，均服从p维Gaussian分布，即对于i=1，…，n，Xi～N（0，Σ），其中Σ为协方差矩阵，根据前文介绍，在图模型框架下，估计网络对应于估计协方差矩阵的逆Θ=Σ-1=（θij）。对于无标度网络G=（V，E），节点的度服从幂律分布，即P（d）∝d-α，其中α＞0为尺度参数，d为节点的度。在图模型中，节点i的度为，即为θ-i=（θi1，…，θi，i-1，θi，i+1，…，θip）T的l0拟范数，考虑到l0的不连续性及其带来的组合优化问题，[19]中用来近似节点i的度，其中0＜q＜1。在节点度是无标度先验的假设下，[19]提出了如下罚似然方法估计Θ，使得Θ对应网络具有无标度特征：

其中，λ为调控参数，εi＞0保证log函数有定义。我们注意到（1）中罚函数为log型和Lq型罚函数的复合，这相当于引入网络中节点度的分布信息和边的稀疏先验。问题（1）为非凸优化问题，[19]采用重赋权迭代算法求解，重赋权迭代算法是一种高效的求解算法，它通过每次寻找目标函数的Minorization函数将原始非凸问题转化为易于求解的凸优化问题[19]。

2.网络社区发现(www.xing528.com)

网络社区发现方法大致可以分为三类：第一类是基于算法的社区发现方法，这类方法是启发式的，是基于一个事先给定的算法进行节点聚类。第二类是基于优化的社区发现方法，即构造了一个衡量聚类好坏的准则，然后对该准则进行优化。第三类是基于统计模型的社区发现方法，即对网络社区构造了一个统计模型，并基于模型对社区进行统计推断。我们采用基于模块度（Modularity）的优化方法[20]对PM2.5污染网络进行社区发现。模块度是一种衡量社区划分质量的测度，当社区间连接紧密、社区内连接稀疏时，模块度的值高，反之，模块度的值低。具体地，假设网络为由n个节点、m条边构成的无向网络，邻接矩阵为A，节点i的度为di，所在社区为ci，那么模块度Q定义如下：

其中，δ（cv，cw）=I｛cv=cw｝。可以看到，（2）式表示了社区内部的边数超过节点之间以正比于节点度的乘积的概率随机连接的边数，换句话说，表示有社区和无社区（Null Model）情况下边数的差异。因此，边数差异越明显，Q值越大，社区越明显。所以通过最大化Q来寻找社区（ci）。我们注意到（2）式是一个组合优化问题，直接求解是不可行的。有鉴于此，Clauset等[21]提出了快速求解问题（2）的贪婪算法。我们这里采用该方法。具体地，我们采取如下两个步骤进行PM2.5污染网络的社区发现：

步骤一：利用Graphical Lasso[18]构建初始网络。Graphical Lasso对应目标函数如下：

其中，λ为调控参数，我们采用EBIC[22]选择。

步骤二：利用Clauset等提出的快速求解Modularity的算法在步骤一所得网络上进行社区发现。

注意到我们并没有使用无标度图模型所得的网络进行社区发现，其原因归结为无标度图模型聚焦于估计网络的中心，但是无标度图模型不能捕捉网络的社区结构，换句话说，无标度图模型估计得到的网络大部分不具有明显的社区或是基于其的社区发现结果表现不佳。因此，我们采用Graphical Lasso构建初始网络，目的是尽量不破坏原始数据的类结构。