前情回顧：
Gephi網(wǎng)絡(luò)圖極簡教程
Network在單細胞轉(zhuǎn)錄組數(shù)據(jù)分析中的應(yīng)用
Gephi網(wǎng)絡(luò)圖極簡教程
Network在單細胞轉(zhuǎn)錄組數(shù)據(jù)分析中的應(yīng)用
網(wǎng)絡(luò)數(shù)據(jù)統(tǒng)計分析筆記|| 為什么研究網(wǎng)絡(luò)
網(wǎng)絡(luò)數(shù)據(jù)統(tǒng)計分析筆記|| 操作網(wǎng)絡(luò)數(shù)據(jù)
網(wǎng)絡(luò)數(shù)據(jù)統(tǒng)計分析筆記|| 網(wǎng)絡(luò)數(shù)據(jù)可視化
網(wǎng)絡(luò)數(shù)據(jù)統(tǒng)計分析筆記|| 網(wǎng)絡(luò)數(shù)據(jù)的描述性分析

在前面的章節(jié)中我們了解到網(wǎng)絡(luò)圖的構(gòu)建，可視化，以及網(wǎng)絡(luò)結(jié)構(gòu)的特征化描述。從本章開始，我們將進入網(wǎng)絡(luò)圖建模的主題，在網(wǎng)絡(luò)數(shù)據(jù)分析中構(gòu)建與使用模型。本章主要介紹幾種常見的數(shù)學模型，就像我們在學統(tǒng)計建模的時候，先要學習幾個常見的分布模型一樣。關(guān)于統(tǒng)計建模的一般性描述見環(huán)境與生態(tài)統(tǒng)計：R語言應(yīng)用。

所謂的網(wǎng)絡(luò)圖模型是指：

其中是所有可能的圖的集合，是上的一個概率分布，是參數(shù)構(gòu)成的向量，該向量的所有可能取值為。

• 經(jīng)典隨機圖模型
• 伯努利隨機圖模型
• 廣義隨機圖模型
• 基于機制的網(wǎng)絡(luò)圖模型
  • 小世界模型
  • 優(yōu)先連接模型
  • BA隨機圖(igraph::barabasi.game)
• 評估網(wǎng)絡(luò)圖特征的顯著性
  • 評估網(wǎng)絡(luò)社團數(shù)量
  • 評估小世界性

經(jīng)典隨機圖模型

在隨機圖模型（Random graphs）中，我們模仿這樣的一個環(huán)境，假如一個團體中有很多的個體，之后兩個人隨機的認識并且成為朋友，那么隨著時間的推移，這個團體會變成什么樣子呢？或者說這個以人為節(jié)點，邊代表好友關(guān)系的網(wǎng)絡(luò)會是什么樣子的呢？

• 隨機圖模型是后續(xù)的基礎(chǔ)，也是重要的參考模型
• 有助于我們通過比較更深入地認識真實網(wǎng)絡(luò)數(shù)據(jù)
• 隨機圖模型能幫助我們理解隨機過程對網(wǎng)絡(luò)結(jié)構(gòu)能夠產(chǎn)生多大程度的影響。

正式地講，隨機圖模型通常是指一個給定了集合及其上的均勻概率分布的模型。其重要作用和完備性就像統(tǒng)計建模中的均勻分布一樣。

比較常見的隨機網(wǎng)路模型是Erdos-Renyi model，可以通過sample_gnp來構(gòu)建。

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)

可以看到我們生成的隨機圖不是連通的，有一個 巨型組件。

經(jīng)典隨機網(wǎng)絡(luò)的性質(zhì)包括：平均度與期望值比較接近，度分布均勻，節(jié)點對之間最短路徑上的節(jié)點相對較少等。

# 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
# ---

廣義隨機圖模型

廣義隨機圖模型是經(jīng)典隨機圖模型的一般化，具體地：

• 定義圖集合,包含階數(shù)為，且具有給定特征的所有圖；
• 為每個圖分配相同的概率。

在Erdos-Renyi模型之外，最常選擇的特征是固定度序列。假設(shè)對于節(jié)點數(shù)為8，一半節(jié)點的度為2，另4個節(jié)點的度為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
# ---

可見兩個圖并非同構(gòu)。

我們可以從構(gòu)建一個與已知圖序列相同的圖：

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
# ---

模擬圖直徑減少一半，之前的聚類也減少了。

基于機制的網(wǎng)絡(luò)圖模型

隨機圖模型為我們描述了在不受任何條件控制的條件下的圖，可理解為數(shù)學模型的背景模型，但是現(xiàn)實世界里的圖往往是由特定結(jié)構(gòu)的?；跈C制的網(wǎng)絡(luò)圖模型 把我們帶入了現(xiàn)實世界。其中最著名的需要所小世界模型了。

• Small-World Model

小世界模型最經(jīng)典的特征是既具有規(guī)則網(wǎng)絡(luò)的高聚集性，又具有類似隨機網(wǎng)絡(luò)的小直徑。相較隨機圖模型，小世界模型能夠更好地反映真實網(wǎng)絡(luò)的情況。就像我們?nèi)祟惿鐣粯?，人以群分，六度分隔?/p>

例如在寫本筆記的時候：

媒體經(jīng)常提到COVID-19呼吸道疾病的病例和死亡人數(shù)呈“指數(shù)”增長，但這些數(shù)字暗示了其他東西，一個可能具有冪律屬性的“小世界”網(wǎng)絡(luò)。這將大大不同于疾病的指數(shù)增長路徑。

在介紹隨機網(wǎng)絡(luò)時提到，隨機網(wǎng)絡(luò)無法解釋真實網(wǎng)絡(luò)中存在的一些情況：局部集聚（較高的集聚系數(shù)）和三元閉合（朋友的朋友是朋友）。從網(wǎng)絡(luò)結(jié)構(gòu)來看，隨機網(wǎng)絡(luò)與真實網(wǎng)絡(luò)的一大差異便是過低的集聚系數(shù)，所以在隨機網(wǎng)絡(luò)模型基礎(chǔ)上進行改進時，需要要著重考慮的便是——如何在保留小網(wǎng)絡(luò)直徑這一特點的同時提高集聚系數(shù)，使得構(gòu)建的模型能夠?qū)W(wǎng)絡(luò)局部結(jié)構(gòu)進行更好的刻畫。

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()

小世界的性質(zhì)：

# 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")

• 優(yōu)先連接模型

優(yōu)先連接”（preferential attachment）指的是進入一個網(wǎng)絡(luò)的新節(jié)點傾向于與節(jié)點度高的節(jié)點相連接。反過來說，一個節(jié)點如果已經(jīng)接受了很多連接，那么它就越容易被新來的節(jié)點所連接。

優(yōu)先連接現(xiàn)象最早是在1925年，由英國統(tǒng)計學家George Udny Yule研究的。后來科學計量之父Derek J. de Solla Price在1976年也研究了這一現(xiàn)象，并把它叫做積累優(yōu)勢（cumulative advantage）。不過，描述優(yōu)先連接最著名的模型是Albert-Laszlo Barabasi和Reka Albert提出的，所以也被叫做Barabási–Albert模型或BA模型。它的基本形式非常簡明：一個新的節(jié)點i連接到網(wǎng)絡(luò)里某個已有節(jié)點j的概率，就是節(jié)點j的度占全部已有節(jié)點的度之和的比重。

BA模型的節(jié)點度符合冪律分布，生成的是一個無標度網(wǎng)絡(luò)（scale-free network）。
網(wǎng)絡(luò)無標度性的形成有兩個基本的要素：一是網(wǎng)絡(luò)生長，也就是新的節(jié)點加入網(wǎng)絡(luò)的過程；二是網(wǎng)絡(luò)生長過程當中的優(yōu)先連接。

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網(wǎng)絡(luò)的性質(zhì)

# 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
# ---

評估網(wǎng)絡(luò)圖特征的顯著性

如開頭所言，隨機網(wǎng)絡(luò)作為網(wǎng)絡(luò)的背景，它經(jīng)常用來評估網(wǎng)絡(luò)特征的顯著性：即，待觀測的網(wǎng)絡(luò)與隨機網(wǎng)絡(luò)有多大程度的不一樣？

假設(shè)我們有一個來自某種觀測的圖，此處稱為，而我們對某些結(jié)構(gòu)特征感興趣，不妨稱為。在很多情況下，自然會考慮是否是顯著的，即在某種意義上是不尋常的和超預(yù)期的。這一過程很像我們的統(tǒng)計推斷過程統(tǒng)計推斷概述。

• 評估網(wǎng)絡(luò)社團數(shù)量

生成參考分布

• 與待觀測圖相同規(guī)模和階數(shù)的圖
• 與待觀測圖有相同的度分布

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

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

# 經(jīng)典隨機圖
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"))

而真實的我們數(shù)據(jù)的社團數(shù)是：

length(cluster_fast_greedy(karate))
[1] 3

可以說是很顯著的了。這時，你要問為什么？

• 評估小世界性

評估小世界性的一種經(jīng)典方法是：針對待觀測網(wǎng)絡(luò)以及可能觀測到的/經(jīng)過適當修飾的經(jīng)典隨機圖，比較兩者聚類系數(shù)和平均（最短）路徑的長度。如果出現(xiàn)小世界性：

• 聚類系數(shù)大于隨機圖
• 平均路徑基本相同

評估有向圖的小世界性：

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  有向圖聚類系數(shù) 
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 具有一定程度的小世界性質(zhì)。

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