网络数据统计分析笔记||网络图的数学模型

前情回顾：
Gephi网络图极简教程
 Network在单细胞转录组数据分析中的应用
 Gephi网络图极简教程
 Network在单细胞转录组数据分析中的应用
 网络数据统计分析笔记|| 为什么研究网络
 网络数据统计分析笔记|| 操作网络数据
 网络数据统计分析笔记|| 网络数据可视化
 网络数据统计分析笔记|| 网络数据的描述性分析

在前面的章节中我们了解到网络图的构建，可视化，以及网络结构的特征化描述。从本章开始，我们将进入网络图建模的主题，在网络数据分析中构建与使用模型。本章主要介绍几种常见的数学模型，就像我们在学统计建模的时候，先要学习几个常见的分布模型一样。关于统计建模的一般性描述见环境与生态统计：R语言应用。

所谓的网络图模型是指：
$\{P_\theta(G),G \in \scr G : \theta \in \Theta \rbrace$

其中 $\scr G$ 是所有可能的图的集合， $P_\theta$ 是 $\scr G$ 上的一个概率分布， $\theta$ 是参数构成的向量，该向量的所有可能取值为 $\Theta$ 。

经典随机图模型
伯努利随机图模型
广义随机图模型
基于机制的网络图模型
- 小世界模型
- 优先连接模型
- BA随机图(igraph::barabasi.game)
评估网络图特征的显著性
- 评估网络社团数量
- 评估小世界性

经典随机图模型

在随机图模型（Random graphs）中，我们模仿这样的一个环境，假如一个团体中有很多的个体，之后两个人随机的认识并且成为朋友，那么随着时间的推移，这个团体会变成什么样子呢？或者说这个以人为节点，边代表好友关系的网络会是什么样子的呢？

随机图模型是后续的基础，也是重要的参考模型
有助于我们通过比较更深入地认识真实网络数据
随机图模型能帮助我们理解随机过程对网络结构能够产生多大程度的影响。

正式地讲，随机图模型通常是指一个给定了集合 $\scr G$ 及其上的均匀概率分布 $P$ 的模型。其重要作用和完备性就像统计建模中的均匀分布一样。

比较常见的随机网路模型是Erdos-Renyi model，可以通过sample_gnp来构建。

library(sand)
set.seed(42)
?sample_gnp
g.er <- sample_gnp(100, 0.02)
layouts <- grep("^layout_", ls("package:igraph"), value=TRUE)[-1] 
# Remove layouts that do not apply to our graph.
layouts <- layouts[!grepl("bipartite|merge|norm|sugiyama|tree", layouts)]

length(layouts)
par(mfrow=c(3,5), mar=c(1,1,1,1))
for (layout in layouts) {
   print(layout)
   l <- do.call(layout, list(g.er)) 
   plot(g.er, edge.arrow.mode=0, layout=l, main=layout) }
# CHUNK 2
is_connected(g.er)
# ---
## [1] FALSE
# ---

查看图中组件和团的情况

table(sapply(decompose(g.er), vcount)) # 组件的普查

 1  2  3  4 71 
15  2  2  1  1 

table(sapply(cliques(g.er), length))  # 团的普查

  1   2   3 
100  95   1

par(mfrow=c(3,7), mar=c(1,1,1,1))
for(i in 1: length(decompose(g.er))){
   l <- do.call(layout_with_lgl, list(decompose(g.er)[[i]]) )
   plot(decompose(g.er)[[i]], edge.arrow.mode=0, layout=l, main=i)
   box(which = "figure", col = "blue", lwd = 1)
   
}
box(which = "outer", col = "red",  lwd = 10)

可以看到我们生成的随机图不是连通的，有一个巨型组件。

经典随机网络的性质包括：平均度与期望值比较接近，度分布均匀，节点对之间最短路径上的节点相对较少等。

# CHUNK 4
mean(degree(g.er))
# ---
## [1] 1.9
# ---

# CHUNK 5
hist(degree(g.er), col="lightblue",
   xlab="Degree", ylab="Frequency", main="")

# CHUNK 6
mean_distance(g.er)
# ---
## [1] 5.276511
# ---
diameter(g.er)
# ---
## [1] 14
# ---

# CHUNK 7
transitivity(g.er)
# ---
## [1] 0.01639344
# ---

广义随机图模型

广义随机图模型是经典随机图模型的一般化，具体地：

定义图集合 $\scr G$ ,包含阶数为 $N_v$ ，且具有给定特征的所有图；
为每个图 $G\in \scr G$ 分配相同的概率。

在Erdos-Renyi模型之外，最常选择的特征是固定度序列。假设对于节点数为8，一半节点的度为2，另4个节点的度为3，从满足条件的图集合中均匀抽取两个。

degs <- c(2,2,2,2,3,3,3,3)
#?sample_degseq
g1 <- sample_degseq(degs, method="vl")
g2 <- sample_degseq(degs, method="vl")
par(mfrow=c(1,2), mar=c(1,1,1,1))
plot(g1, vertex.label=NA)
plot(g2, vertex.label=NA)

isomorphic(g1, g2)
# ---
## [1] FALSE
# ---

c(ecount(g1), ecount(g2))
# ---
## [1] 10 10
# ---

可见两个图并非同构。

我们可以从构建一个与已知图序列相同的图：

data(yeast)
degs <- degree(yeast)
fake.yeast <- sample_degseq(degs, method=c("vl"))
all(degree(yeast) == degree(fake.yeast))
[1] TRUE

plot(yeast,vertex.label=NA)
plot(fake.yeast,vertex.label=NA)

diameter(yeast)
# ---
## [1] 15
# ---
diameter(fake.yeast)
# ---
## [1] 8
# ---

# CHUNK 13
transitivity(yeast)
# ---
## [1] 0.4686178
# ---
transitivity(fake.yeast)
# ---
## [1] 0.04026804
# ---

模拟图直径减少一半，之前的聚类也减少了。

基于机制的网络图模型

随机图模型为我们描述了在不受任何条件控制的条件下的图，可理解为数学模型的背景模型，但是现实世界里的图往往是由特定结构的。基于机制的网络图模型把我们带入了现实世界。其中最著名的需要所小世界模型了。

Small-World Model

小世界模型最经典的特征是既具有规则网络的高聚集性，又具有类似随机网络的小直径。相较随机图模型，小世界模型能够更好地反映真实网络的情况。就像我们人类社会一样，人以群分，六度分隔。

例如在写本笔记的时候：

媒体经常提到COVID-19呼吸道疾病的病例和死亡人数呈“指数”增长，但这些数字暗示了其他东西，一个可能具有幂律属性的“小世界”网络。这将大大不同于疾病的指数增长路径。

在介绍随机网络时提到，随机网络无法解释真实网络中存在的一些情况：局部集聚（较高的集聚系数）和三元闭合（朋友的朋友是朋友）。从网络结构来看，随机网络与真实网络的一大差异便是过低的集聚系数，所以在随机网络模型基础上进行改进时，需要要着重考虑的便是——如何在保留小网络直径这一特点的同时提高集聚系数，使得构建的模型能够对网络局部结构进行更好的刻画。

g.ws <- sample_smallworld(1, 25, 5, 0.05)
par(mfrow=c(3,5), mar=c(1,1,1,1))
for (layout in layouts) {
   print(layout)
   l <- do.call(layout, list(g.ws)) 
   plot(g.ws, edge.arrow.mode=0, layout=l, main=layout) }

dev.off()

小世界的性质：

# CHUNK 15
g.lat100 <- sample_smallworld(1, 100, 5, 0)
transitivity(g.lat100)
# ---
## [1] 0.6666667
# ---

# CHUNK 16
diameter(g.lat100)
# ---
## [1] 10
# ---
mean_distance(g.lat100)
# ---
## [1] 5.454545
# ---

# CHUNK 17
g.ws100 <- sample_smallworld(1, 100, 5, 0.05)
diameter(g.ws100)
# ---
## [1] 5
# ---
mean_distance(g.ws100)
# ---
## [1] 2.748687
# ---
transitivity(g.ws100)
# ---
## [1] 0.5166263
# ---

steps <- seq(-4, -0.5, 0.1)
len <- length(steps)
cl <- numeric(len)
apl <- numeric(len)
ntrials <- 100
function(x) {
  for (i in (1:len)) {
   cltemp <- numeric(ntrials)
   apltemp <- numeric(ntrials)
   for (j in (1:ntrials)) {
     g <- sample_smallworld(1, 1000, 10, 10^steps[i])
     cltemp[j] <- transitivity(g)
     apltemp[j] <- mean_distance(g)
   }
   cl[i] <- mean(cltemp)
   apl[i] <- mean(apltemp)
 }
}
cl <- c(0.710063379997766, 0.709978692214238, 0.709907143256545, 0.709724130686251,
0.709438119171845, 0.709084388626035, 0.708846266062516, 0.70839051192321,
0.707759691875033, 0.707113107172047, 0.706190905933217, 0.705111695935303,
0.703784841816035, 0.702347962546443, 0.699998029666335, 0.696966876979092,
0.693486200400489, 0.689434391992611, 0.683909800354255, 0.676998368887877,
0.669280399418907, 0.657931797843006, 0.645296561052957, 0.628819148376097,
0.609573848258676, 0.585356133848734, 0.55633728515175, 0.521088308467222,
0.480754321558662, 0.433680652553029, 0.378558209487185, 0.318761294080951,
0.25616767320489, 0.193136725458156, 0.134572797222469, 0.0849655822333312)
apl <- c(19.2309712512513, 18.1696153953954, 17.7627795195195, 16.6679303103103,
14.5979843643644, 12.8854347747748, 12.1038251251251, 11.2341607007007,
9.82806862862863, 9.11716512512512, 8.24078900900901, 7.61965115115115,
7.0006803003003, 6.52964078078078, 6.03792086086086, 5.60314024024024,
5.26338822822823, 4.98922616616617, 4.70584088088088, 4.45406394394394,
4.26384696696697, 4.05971747747748, 3.89097137137137, 3.72829705705706,
3.58416728728729, 3.45041687687688, 3.33307095095095, 3.22760666666667,
3.12912454454454, 3.0362582982983, 2.94736488488488, 2.87122124124124,
2.80854338338338, 2.75799567567568, 2.71760948948949, 2.68516954954955)

# CHUNK 19
plot(steps, cl/max(cl), ylim=c(0, 1), lwd=3, type="l",
   col="blue", xlab=expression(log[10](p)),
   ylab="Clustering and Average Path Length")
lines(steps, apl/max(apl), lwd=3, col="red")

优先连接模型

优先连接”（preferential attachment）指的是进入一个网络的新节点倾向于与节点度高的节点相连接。反过来说，一个节点如果已经接受了很多连接，那么它就越容易被新来的节点所连接。

优先连接现象最早是在1925年，由英国统计学家George Udny Yule研究的。后来科学计量之父Derek J. de Solla Price在1976年也研究了这一现象，并把它叫做积累优势（cumulative advantage）。不过，描述优先连接最著名的模型是Albert-Laszlo Barabasi和Reka Albert提出的，所以也被叫做Barabási–Albert模型或BA模型。它的基本形式非常简明：一个新的节点i连接到网络里某个已有节点j的概率，就是节点j的度占全部已有节点的度之和的比重。

BA模型的节点度符合幂律分布，生成的是一个无标度网络（scale-free network）。
网络无标度性的形成有两个基本的要素：一是网络生长，也就是新的节点加入网络的过程；二是网络生长过程当中的优先连接。

set.seed(42)
?sample_pa
g.ba <- sample_pa(100, directed=FALSE)

# CHUNK 21
par(mfrow=c(3,5), mar=c(1,1,1,1))
for (layout in layouts) {
   print(layout)
   l <- do.call(layout, list(g.ba)) 
   plot(g.ws, edge.arrow.mode=0, layout=l, main=layout) }

hist(degree(g.ba), col="lightblue",
   xlab="Degree", ylab="Frequency", main="")

ba网络的性质

# CHUNK 23
summary(degree(g.ba))
# ---
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    1.00    1.00    1.98    2.00    9.00
# ---

# CHUNK 24
mean_distance(g.ba)
# ---
## [1] 5.815556
# ---
diameter(g.ba)
# ---
## [1] 12
# ---

# CHUNK 25
transitivity(g.ba)
# ---
## [1] 0
# ---

评估网络图特征的显著性

如开头所言，随机网络作为网络的背景，它经常用来评估网络特征的显著性：即，待观测的网络与随机网络有多大程度的不一样？

假设我们有一个来自某种观测的图，此处称为 $G^{obs}$ ，而我们对某些结构特征感兴趣，不妨称为 $\eta(·)$ 。在很多情况下，自然会考虑 $\eta(G^{obs})$ 是否是显著的，即在某种意义上是不寻常的和超预期的。这一过程很像我们的统计推断过程统计推断概述。

评估网络社团数量

生成参考分布

与待观测图相同规模和阶数的图
与待观测图有相同的度分布


# CHUNK 26
data(karate)
nv <- vcount(karate)
ne <- ecount(karate)
degs <- degree(karate)

# CHUNK 27
ntrials <- 1000 # 模拟1000次

# 经典随机图
num.comm.rg <- numeric(ntrials)
for(i in (1:ntrials)){
   g.rg <- sample_gnm(nv, ne)
   c.rg <- cluster_fast_greedy(g.rg)
   num.comm.rg[i] <- length(c.rg)
}

# 广义随机图
num.comm.grg <- numeric(ntrials)
for(i in (1:ntrials)){
   g.grg <- sample_degseq(degs, method="vl")
   c.grg <- cluster_fast_greedy(g.grg)
   num.comm.grg[i] <- length(c.grg)
}

rslts <- c(num.comm.rg,num.comm.grg)
indx <- c(rep(0, ntrials), rep(1, ntrials))
counts <- table(indx, rslts)/ntrials
barplot(counts, beside=TRUE, col=c("blue", "red"),
   xlab="Number of Communities",
   ylab="Relative Frequency",
   legend=c("Fixed Size", "Fixed Degree Sequence"))

而真实的我们数据的社团数是：

length(cluster_fast_greedy(karate))
[1] 3

可以说是很显著的了。这时，你要问为什么？

评估小世界性

评估小世界性的一种经典方法是：针对待观测网络以及可能观测到的/经过适当修饰的经典随机图，比较两者聚类系数和平均（最短）路径的长度。如果出现小世界性：

聚类系数大于随机图
平均路径基本相同

评估有向图的小世界性：

library(igraphdata)
data(macaque)
summary(macaque)
# ---
## IGRAPH f7130f3 DN-- 45 463 -- 
## + attr: Citation (g/c), Author (g/c), shape (v/c), 
##   name (v/c)
# ---

# CHUNK 32  有向图聚类系数 
clust_coef_dir <- function(graph) {
   A <- as.matrix(as_adjacency_matrix(graph))
   S <- A + t(A)
   deg <- degree(graph, mode=c("total"))
   num <- diag(S %*% S %*% S)
   denom <- diag(A %*% A)
   denom <- 2 * (deg * (deg - 1) - 2 * denom)
   cl <- mean(num/denom)
   return(cl)
}

# CHUNK 33 # 模拟评估 
ntrials <- 1000
nv <- vcount(macaque)
ne <- ecount(macaque)
cl.rg <- numeric(ntrials)
apl.rg <- numeric(ntrials)
for (i in (1:ntrials)) {
   g.rg <- sample_gnm(nv, ne, directed=TRUE)
   cl.rg[i] <- clust_coef_dir(g.rg)
   apl.rg[i] <- mean_distance(g.rg)
}

# CHUNK 34
summary(cl.rg)
# ---
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##  0.2159  0.2302  0.2340  0.2340  0.2377  0.2548
# ---

# CHUNK 35
summary(apl.rg)
# ---
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   1.810   1.827   1.833   1.833   1.838   1.858
# ---

# CHUNK 36
clust_coef_dir(macaque)
# ---
## [1] 0.5501073
# ---
mean_distance(macaque)
# ---
## [1] 2.148485
# ---

0.5501073 > 0.2548 ; 2.148485 > 1.858 具有一定程度的小世界性质。

https://zhuanlan.zhihu.com/p/146499763
https://zhuanlan.zhihu.com/p/205012648
https://blog.csdn.net/limiyudianzi/article/details/81632139
http://economics.mit.edu/files/4623#:~:text=Generalized%20random%20graph%20models%20%28such%20as%20the%20con,combines%20high%20clustering%20with%20short%20path%20lengths%20is
https://ocw.mit.edu/courses/economics/14-15j-networks-spring-2018/lecture-and-recitation-notes/MIT14_15JS18_lec12.pdf
https://zhuanlan.zhihu.com/p/37121528
https://www.zdnet.com/article/graph-theory-suggests-covid-19-might-be-a-small-world-after-all/
https://www.sohu.com/a/402313767_169228