Statistical Analysis of Network Data with R(第二版) 上机实操 4

4.2.3 边的中介中心度（Edge Betweenness centrality）

为每条边分配一个值，该值反映穿过该边的最短路径的数量。

> library(igraphdata)
> data("karate")
> ed <- edge.betweenness(karate) %>% order(decreasing = TRUE)
> E(karate)[ed][1:3]
+ 3/78 edges from 4b458a1 (vertex names):
[1] Actor 20--John A   Mr Hi   --Actor 20 Mr Hi   --Actor 32

作图

> ed <- edge.betweenness(karate)
> cl <- c(rep("blue",52),"red",rep("blue",25))
> plot(karate,edge.color=cl,edge.width=log(ed))

加上节点的中心度指标，似乎边中介中心度最高的边，连接的两个节点，节点中心度的层级并不相同

> plot(karate,edge.color=cl,edge.width=log(ed),
+      vertex.size=10*sqrt(authority.score(karate)$vector))

4.3 表征网络内聚性 Characterizing Network Cohesion

对应网络图中关系定义了边，内聚性反映顶点子集“粘在一起”的程度
内聚性意味着一个社会网络包含许多联系。人与人之间的更多联系会产生更紧密的结构
其指标有：三元组（ triads）、巨型组件 (giant components)、团 (cliques)

4.3.1 子图和普查

定义图内聚性的一种方法是先定义感兴趣的子图的类型。而子图中最典型的例子是团。
团是完全子图，子图中的所有顶点都由边连接，因此是完全内聚的顶点子集。
在社会网络分析中，重叠团的结构被认为代表社会圈而不是单个团，被认为是有内聚性的子群。团或完整的三元组是网络中最密集的部分或“骨骼”，因此
重叠团的结构被认为是其“骨架”。

（1）团 cliques

igraph包的cliques()函数，在没有设置最大值和最小值参数时，会返回所有子图，而三元组必然以二元组为基础，四元组首先必需是三元组，所以叫普查（census）：
书中示例

> table(sapply(cliques(karate),length))1  2  3  4  5
34 78 45 11  2

其实，igraph有一个函数可以实现该功能

> clique_size_counts(karate)
[1] 34 78 45 11  2

上例中，一元组实际是节点的数量，二元组实际是边的数量，在一些网络分析理论中是忽略不计的，直接把三元组作为最小的图。图karate中最大的团是五元组，共有2个：
书中示例，找出最大团

> cliques(karate)[sapply(cliques(karate),length)==5]
[[1]]
+ 5/34 vertices, named, from 4b458a1:
[1] Mr Hi    Actor 2  Actor 3  Actor 4  Actor 14[[2]]
+ 5/34 vertices, named, from 4b458a1:
[1] Mr Hi   Actor 2 Actor 3 Actor 4 Actor 8

同样，igraph有一个函数可以实现该功能

> largest.cliques(karate)
[[1]]
+ 5/34 vertices, named, from 4b458a1:
[1] Actor 2 Mr Hi   Actor 4 Actor 3 Actor 8[[2]]
+ 5/34 vertices, named, from 4b458a1:
[1] Actor 2  Mr Hi    Actor 4  Actor 3  Actor 14

结果完全相同。书中示例也不是无用，比如要找出所有第二大团：

> cliques(karate)[sapply(cliques(karate),length)==4]
[[1]]
+ 4/34 vertices, named, from 4b458a1:
[1] Mr Hi   Actor 2 Actor 3 Actor 4[[2]]
+ 4/34 vertices, named, from 4b458a1:
[1] Actor 24 Actor 30 Actor 33 John A  ···[[11]]
+ 4/34 vertices, named, from 4b458a1:
[1] Mr Hi   Actor 2 Actor 3 Actor 8

如上面所说，如果排除尺寸是1的团，大尺寸的团必定包含小尺寸的团，那些不被任何更多图包含的团，叫极大团（ maximum cliques）
书中示例

> table(sapply(maximal.cliques(karate), length))2  3  4  5
11 21  2  2

igraph提供了几个关于团的计数函数

# 极大团个数
> maximal.cliques.count(karate)
[1] 36
> count_max_cliques(karate)
[1] 36
# 最大团的尺寸
> clique_num(karate)
[1] 5

在pajek中，普查团的工作被分成几个步骤，首先要设置团的形状，比如要普查图的三元组，需要先做一个三元组的图，用.net格式：

可视化展示结果是这样的：

因为在普查的过程中，有向图和无向团差别极大，或者将要普查的元组设置为有向的，比如：

出图是这样的：

或者将原图转换为无向图：

把要普查的元组文件放到network下拉框的第一里，原图放到第二个network下拉框中，注意次序不能换，然后执行：
Networks> Fragment (First in Second)
即可完成普查

（2）核 core

完全子图：满足任意两点都恰有一条边相连的子图，也叫团，就是说，团中的节点不仅彼此间都有边相连，每个节点的度书都大于等于团的尺寸；
把团的条件弱化一下：子图中的节点的度数都大于等于某个数，但这些节点之间不必相互连通，这种子图，叫核。核中节点的最小度数为k，此时，该子图便称为k-核。
从网络中移除最低的 k 核，直到网络分解成相对密集的组件核，直到网络分解成相对密集的组件。然后，每个组件被认为是一个内聚子组，因为它
在组件内至少有 k 个邻居。在（非常）大的网络中，这是寻找有内聚性的子组的有效方法这是寻找有内聚性的子组的有效方法。
之所以说k-核不一定是连通的，是因为一个图中，可以同时存在多个3-核或n-核，他们位于图的不同位置，核与核之间是不连通的，但核内部是连通的，不然成不了核
plot(karate,vertex.color=coreness(karate))

用sna包制作辐射图，注意需要先把igraph包的图文件转换成邻接矩阵，再借助sna包把邻接矩阵转换成sna包能处理的矩阵格式，最后用sna包出图

> a <- as_adjacency_matrix(karate,sparse = FALSE)
> p <- network::as.network.matrix(a)
> sna::gplot.target(p,coreness(karate),usearrows=FALSE,vertex.col=coreness(g))

在 Pajek 中，使用Network> Create Partition> k-Core> Input, Output, All 检测 k-core。该命令产生一个partition ，将每个顶点分配给它出现的最高 k 核。

上图顶点颜色和顶点内的数字显示了 k 核partition 。
使用 k 核partition ，您可以轻松地从网络中删除低 k 核以提取网络中最密集的部分。
比如，只保留k-核大于3的子图，在 Partition 下拉列表中选择前面命令生成的 k-核partition 并执行
Operations> Network + Partition> Extract> SubNetwork Induced by Union of Selected Clusters
在弹出窗口中，输入“3-4”

结果是这样的：

随后，使用 Network> Create Partition> Components> Strong 命令检查选定的 k 核级别是否拆分为两个或多个组件。

（3）二元组 dyad

虽然二元组代表边，但在有向图中，通过对二元组的普查，可以发现单方行为：

> is.simple(aidsblog)
[1] FALSE
> al <- simplify(aidsblog)
> dyad_census(al)
$mut #双向
[1] 3$asym #单向
[1] 177$null #空（无边）
[1] 10405

这表明博客间的引用是单向行为。

（4）三元组 triad

igraph包同时提供了三元组普查函数triad_census( )，有向三元组共有16种状态，所以返回结果分16个元素

> triad_census(al)[1] 484621  20588    429   2195     39     74      1    112      4      0      2      0
[13]     15      0      0      0

具体16种三元组类型：

4.3.2 密度和相对频率

（1）密度 density

密度：实际存在的边与可能的边的比例。
设H是一个图的子图，无环简单无向图的密度公式：
den(H)=∣EH∣∣VH∣(∣VH∣−1)/2den(H)=\frac{|E_H|}{|V_H|(|V_H|-1)/2}den(H)=∣VH∣(∣VH∣−1)/2∣EH∣
有向图的密度公式：
den(H)=∣EH∣∣VH∣(∣VH∣−1)den(H)=\frac{|E_H|}{|V_H|(|V_H|-1)}den(H)=∣VH∣(∣VH∣−1)∣EH∣
团clique是完全子图，所以它的密度是1.
计算图的密度很简单：

> graph.density(karate)
[1] 0.1390374

但edge_density(graph, loops = FALSE)的输入参数是图，所以要计算图中的一部分，需要先提取成单独的图：

> p <-  induced.subgraph(karate,igraph::neighborhood(karate,1,1)[[1]])
> graph.density(p)
[1] 0.25

（2）相对频率（聚类系数、传递性、互惠性）

所谓相对，指实现数与可能数的比值。
聚类系数描述一个图中的某节点与其相连节点之间聚集成团的程度的一个系数。
聚类系数的目标是比较群组的聚合紧密程度与其能够达到的聚合紧密程度。
一个节点的局部聚类系数体现的是其邻节点也相互连通的可能性。实际上就是算节点u与邻居构成实际组成的三角形数除上最大可能三角形个数，翻译成简单的话就是我的朋友之间相互认识程度。
全局聚类系数以节点三重性 triplets为基础。一个三元组是由两个(开三元组)或三个(闭三元组)无向联系相连接的三个节点。因此，一个三角形图包括三个闭合的三联体，所以下面的公式中，有3倍系数
clT(G)=3TΔ(G)T3(G)cl_T(G)=\frac{3T _\Delta (G)}{T_3(G)}clT(G)=T3(G)3TΔ(G)
其中 T(G) 是图 G 中三角形的数量，而T3(G) 是连接三元组的数量（即，由两条边连接的三个顶点的子图，有时也称为 2 星）。
R语言中，可用：

> transitivity(karate)
[1] 0.2556818
> transitivity(karate,"local",vids = c(1,34))
[1] 0.1500000 0.1102941

如果计算相对聚类的单位不是三元组，而是二元组，指标就被称为互惠性：

# 双向边与非双向边的比例
> reciprocity(aidsblog)
[1] 0.03243243
# 双向边与所以边的比例
> reciprocity(aidsblog,mode = "ratio")
[1] 0.01648352

在 Pajek 中，计算 triad census 非常容易：只需使用Network> Info> Triadic Census 命令。一个对话框询问是否应该报告模型，如果选择此选项，三元组类型（“Type”）、它们的实际频率（“Number of triads [ni]”）以及随机预期的频率（“Expected [ei]”）被报告。此外，显示了实际和预期三元组
数之间的相对差异（“[ni-ei]/ei”），以及检验实际频率等于预期频率的假设的卡方统计量。如果预期频率较低，则此统计数据不可靠。三元组计数和预期
频率也存储在两个向量中。

4.3.3 连通性、割、流

（1）组件

组件（components ）是图中的最大连通子图。
可以说，图是由一个或几个组件构成。反过来，可以利用组件，把图分割：

> is.connected(karate)
[1] TRUE
> decompose(karate)
[[1]]
IGRAPH 6c75550 UNW- 34 78 -- Zachary's karate club network
+ attr: name (g/c), Citation (g/c), Author (g/c), Faction (v/n), name (v/c),
| label (v/c), color (v/n), weight (e/n)
+ edges from 6c75550 (vertex names):[1] Mr Hi  --Actor 2  Mr Hi  --Actor 3  Mr Hi  --Actor 4  Mr Hi  --Actor 5 + ... omitted several edges> g <- karate+vertex(40)
> is.connected(g)
[1] FALSE
> decompose(g)
[[1]]
IGRAPH 70097a1 UNW- 34 78 -- Zachary's karate club network
+ attr: name (g/c), Citation (g/c), Author (g/c), Faction (v/n), name (v/c),
| label (v/c), color (v/n), weight (e/n)
+ edges from 70097a1 (vertex names):[1] Mr Hi  --Actor 2  Mr Hi  --Actor 3  Mr Hi  --Actor 4  Mr Hi  --Actor 5
+ ... omitted several edges[[2]]
IGRAPH 70097bc UNW- 1 0 -- Zachary's karate club network
+ attr: name (g/c), Citation (g/c), Author (g/c), Faction (v/n), name (v/c),
| label (v/c), color (v/n), weight (e/n)
+ edges from 70097bc (vertex names):

对大型网络，组件分割图似乎不像上面这么容易，需要先筛选一下：

data("yeast")
plot(yeast,vertex.size=3,vertex.label=NA)

> table(sapply(decompose(yeast),vcount))2    3    4    5    6    7 2375 63   13    5    6    1    3    1

R按包含节点数量从多到少的次序存放结果，所以第一个是包含2375个节点的最大组件

> decompose(yeast)[[1]]
IGRAPH d6a950b UN-- 2375 11693 -- Yeast protein interactions, von Mering et al.
+ attr: name (g/c), Citation (g/c), Author (g/c), URL (g/c), Classes (g/x),
| name (v/c), Class (v/c), Description (v/c), Confidence (e/c)
+ edges from d6a950b (vertex names):[1] YLR197W--YDL014W YOR039W--YOR061W YDR473C--YPR178W YOR332W--YLR447C

如果图本身是连通的，decompose无法再分割，如果不连通，每个组件都会分割为一个图。
在Pajek中，执行 Network> Create Partition> Components，用于查找三种类型的组件：强、弱和强

（2）小世界性

小世界性（small world property），简单说，就是图中顶点越多，网络中顶点之间的平均距离越短

> data("yeast")
> mean_distance(yeast)
[1] 5.095629
> ygc <- decompose(yeast)[[1]]
> mean_distance(ygc)
[1] 5.09597
> diameter(ygc)
[1] 15
> diameter(yeast)
[1] 15
> transitivity(ygc)
[1] 0.4686663
> transitivity(yeast)
[1] 0.4686178

从上面的代码可以看出：

巨型组件直接影响图的属性，必要是可用巨型组件的属性值近似图的属性值
图的最长路径是15，但平均路径只有5
transitivity值0.468意味着将近一半的连通三元组闭合成三角形

（3）割点

组件是一个静态的概念，只有子图具有彼此连通的结构，就构成组件。而割则是一个动态的概念：如果删除k个节点（或k条边），图是不是仍然连通？
如果删除最大不超过k个节点后，剩余的子图仍然是连通的，则称“k节点连通（ k-vertex-connected）”
同样，如果删除最大不超过k条边后，剩余的子图仍然是连通的，则称“k边连通（ k-edge-connected）”
此时，最大数k称为图的节点（边）连通度（vertex (edge) connectivity）

> vertex.connectivity(ygc)
[1] 1
> vertex.connectivity(yeast)
[1] 0
> data("karate")
> vertex.connectivity(karate)
[1] 1

连通度是图的整体属性，降维到具体的节点，如果删除某个节点就会破坏图的连通性，则这个节点就叫割点，又叫关节点

> V(ygc)$color <- "blue"
> V(ygc)[articulation.points(ygc)]$color <- "red"
> plot(ygc,vertex.size=3,vertex.label=NA)

上图中红色点就是R计算的割点

（4）门杰定理

门杰定理(Menger’s theorem)又称“Menger定理”，是关于图的连通性的一个定理，门杰定理断言：若X和Y是图G的两个不交的节点子集，k是一个正整数，则G上存在k条分别以X和Y中的节点为端点而且两两除端点外互不交的路，当且仅当每一个XY分离点集包含至少k个节点，上述的XY分离点集指G的这样一个节点子集，若从G上去掉该子集的所有节点，所得到的图不存在连X中节点与Y中节点的路，这个定理是门杰(K.Menger)于1927年首先发现的，于是，由此而得名。
1927年，Menger证明了图的连通性与图中连接两个顶点的不相交路的数目有关。

shortest.paths(karate,3,8)
graph.maxflow(yeast,1,78)
graph.mincut(yeast,1,78)

4.4 图划分Graph Partitioning

内聚性节点子集是这样一个节点子集：子集内部节点联系紧密，子集之间相对分离。图划分问题，也被称为社团发现 community detection
社团划分的算法比较多，大致可以分为两大类：拓扑分析和流分析。前者一般适用于无向无权网络，思路是社团内部的连边密度要高于社团间。后者适用于有向有权网络，思路是发现在网络的某种流动（物质、能量、信息）中形成的社团结构。这两种分析各有特点，具体应用取决于网络数据本身描述的对象和研究者想要获得的信息。

4.4.1 聚合——层次聚类Hierarchical Clustering

层次聚类分聚合 agglomerative和分裂divisive
层次聚类中常用的模块度 modularity，先将图的节点进行划分partition，当然这个划分不是唯一的。如果原节点被划分成k组，就等于得到k个模块，我们现在只关注模块之间的边，而忽略模块内部的边。注意，只是忽略模块内部的边、不是删除，因为划分只是对节点视角的升维，边的确定最终还要根据具体的节点。
当我们只考虑模块之间的边时，便出现参照系的确定问题，常用的方法是将一个随机数作为参照。由于此时模块数是k，所以确定随机值的节点数是k:
mod(p)=∑k=1k[fkk(p)−fkk∗]2mod(p)=\sum_{k=1}^k[f_{kk}(p)-f_{kk}^*]^2mod(p)=∑k=1k[fkk(p)−fkk∗]2
因为划分不是唯一的，所以上面的公式要得到最优值，需要大量的计算。igraph采用一种快速贪婪算法，并设计了专门函数：

> cluster_fast_greedy(karate)
IGRAPH clustering fast greedy, groups: 3, mod: 0.43
+ groups:$`1`[1] "Actor 9"  "Actor 10" "Actor 15" "Actor 16" "Actor 19" "Actor 21" "Actor 23"[8] "Actor 24" "Actor 25" "Actor 26" "Actor 27" "Actor 28" "Actor 29" "Actor 30"[15] "Actor 31" "Actor 32" "Actor 33" "John A"  $`2`[1] "Mr Hi"    "Actor 2"  "Actor 3"  "Actor 4"  "Actor 8"  "Actor 12" "Actor 13"[8] "Actor 14" "Actor 18" "Actor 20" "Actor 22"$`3`+ ... omitted several groups/vertices

这个函数直接把划分结果告诉我们，他返回类型是igraph的community类，可用membership()函数查看具体成员


> kc <- cluster_fast_greedy(karate)
> length(kc)
[1] 3
> sizes(kc)
Community sizes1  2  3
18 11  5
> membership(kc)Mr Hi  Actor 2  Actor 3  Actor 4  Actor 5  Actor 6  Actor 7  Actor 8  Actor 9 2        2        2        2        3        3        3        2        1
Actor 10 Actor 11 Actor 12 Actor 13 Actor 14 Actor 15 Actor 16 Actor 17 Actor 18 1        3        2        2        2        1        1        3        2
Actor 19 Actor 20 Actor 21 Actor 22 Actor 23 Actor 24 Actor 25 Actor 26 Actor 27 1        2        1        2        1        1        1        1        1
Actor 28 Actor 29 Actor 30 Actor 31 Actor 32 Actor 33   John A 1        1        1        1        1        1        1
> plot(kc,karate)

library(ape)
dendPlot(kc,mode = "phylo")

本部分，plot()自动分色、树形图背后的原理，未找到合适的解释，照搬原书代码。

4.4.2 分裂——谱划分Spectral Partitioning

拉普拉斯特征值法
谱图理论，将图 G 的连通性与某些矩阵的特征分析联系起来。此处介绍拉普拉斯矩阵。
简单地说，拉普拉斯矩阵就是L=D−AL = D − AL=D−A。其中A是邻接矩阵，D是度矩阵。

可见，拉普拉斯矩阵主对角线上的值等于相应节点的度，其他值是邻接矩阵的相反数。
经过推导可以发现拉普拉斯矩阵的第行其实是第个节点在产生扰动时对其他节点产生的收益累积。

> klap <- laplacian_matrix(karate)
> eig <- eigen(klap)
> plot(eig$values,col="blue")

将该拉普拉斯矩阵的特征值绘图：

# 数据karate的属性中，已经设置了Faction属性的值，所以，本例好像监督学习，可以直接研判划分结果是否准确
> faction <- igraph::get.vertex.attribute(karate, "Faction")
# 设置颜色向量
> fcolor <- 1:(length(faction))
> fcolor[1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
[28] 28 29 30 31 32 33 34
# 依据元数据的属性，将节点颜色分成红蓝两色
> fcolor[faction==1] <- "blue"
> fcolor[faction==2] <- "red"
# 先将图转换成拉普拉斯矩阵
> klap <- laplacian_matrix(karate)
# 对矩阵求特征值，结果是34*34矩阵
> eig <- eigen(klap)
# 本例选33列有自己的原因
> fv <- eig$vectors[,33]
# 出图后发现划分完全准确
> plot(fv,pch=16,col=fcolor)
> abline(0,0,lwd=2,col="lightgray")

# 但选29列，图变成这个样子
> fv <- eig$vectors[,29]
> plot(fv,pch=16,col=fcolor)

n阶矩阵有n个特征值，因为特征值是特征多项式的根,n阶方阵的特征多项式是个n次多项式，根据代数基本定理,n次多项式有且只有n个根(重根按重数计算),这些根可能是实数,也可能是复数。n阶实对称矩阵一定有n个特征值，但这n个特征值包括重根。
如果n阶方阵有n个不同的特征值，则该方阵有n个线性无关的特征向量，可以相似对角化。如果n阶实对称矩阵有n个不同的特征值，则该矩阵有n个正交的特征向量，且实对称矩阵一定可以相似对角化。

社团矩阵特征值 leading eigenvector of the community matrix

> g <- make_full_graph(5)%du% make_full_graph(5)%du% make_full_graph(5)
> cluster_leading_eigen(g)
IGRAPH clustering leading eigenvector, groups: 3, mod: 0.67
+ groups:$`1`[1] 1 2 3 4 5$`2`[1]  6  7  8  9 10$`3`[1] 11 12 13 14 15> cluster_leading_eigen(g) %>% membership()[1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3

从上面的示例中可以看到，cluster_leading_eigen( )函数直接用图作为数据参数，而拉普拉斯法则分成两步：先将图转成拉普拉斯矩阵，再对矩阵求特征值。

4.5 同配性和混合Assortativity and Mixing

同配性（ Assortativity）, 用作考察度值相近的节点是否倾向于互相连接。如果总体上度大的节点倾向于连接度大的节点，那么就称网络的度正相关的，或者成网络是同配的；如果总体上度大的节点倾向于连接度小的节点，那么就称网络的度负相关的，或者成网络是异配的。
同配系数（Assortativity coefficient）是一种基于“度”的皮尔森相关系数，用来度量相连节点对的关系。结果r为正值代表具有相同度的点之间有某种协同关系，负值表示具有不同度数的节点间有某种联系。通常来说，r的值在-1到+1之间，+1表示网络具有很好的同配模式，0表示网络是非同配的，-1表示这个网络完全不同配。
可见，同配性指的是以节点的度为标准、度数值的高低是否与边多少存在相关性。所以同配性是“相关系数概念的一种变体”
如果想单纯根据节点度来分析图的同配性，可直接调用函数：

> assortativity.degree(yeast)
[1] 0.4610798

如果想根据某些节点的属性来计算同配性，因为对同配性来说，节点的属性本身是没用意义的，所以，首先需要根据这个属性，对所有节点做一个划分，然后根据这个划分来计算同配性：

> tv[is.na(tv)] <- 1
> assortativity.nominal(yeast,tv)
[1] 0.5232879
> tv <- (V(yeast)$Class=="P")+1
> tv[is.na(tv)] <- 2
> assortativity.nominal(yeast,tv)
[1] 0.5063002

pajek中，度同配相关系数是网络的一个特征，因此命令位于Network> Info> Degree Assortativity中。其他属性的同配系数，需要提供一个向量参数，所以命令位于Operations> Network + Vector> Info> Assortativity中