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

5.2 经典随机图模型：Erdos-Renyi 模型R实现

已知节点总数，以及在两个任意顶点之间绘制一条边的概率

> library(sand)
> set.seed(42)
> ger <- sample_gnp(100,0.02)
> plot(ger,layout=layout_in_circle,vertex.label=NA)

特征一：随机生成图一般不是连通图，但往往有巨型组件。对于组件，直接用decompose函数，会列出所有组件，不直观，先用table查看一下

> is.connected(ger)
[1] FALSE
> table(sapply(decompose(ger),vcount))
 1  2  3  4 71 
15  2  2  1  1 
> plot(decompose(ger)[[1]])

特征二：图平均度与期望值近似
有N个节点的图，一个节点最大的边数是N1，则在概率随机图的边的期望值是（N-1）*P，所以100个节点的图的边期望是99*0.02=1.98
而平均度：

> mean(degree(ger))
[1] 1.9

特征三：度分布均匀
hist(degree(ger),col="lightblue")

特征四：节点对之间的路径长度普遍大

> mean_distance(ger)
[1] 5.276511
> data("karate")
> mean_distance(karate)
[1] 5.754011
> diameter(ger)
[1] 14

特征五：内聚性低

> mean_distance(ger)
[1] 5.276511
> data("karate")
> mean_distance(karate)
[1] 5.754011
> transitivity(ger)
[1] 0.01639344
> transitivity(karate)
[1] 0.2556818

5.3 广义随机图模型

定义图的集合，集合中的图，节点数量相同，节点的度分布相同，但节点的连接顺序不同。

dec <- c(2,2,2,2,3,3,3,3)
g1 <- sample_degseq(dec,method = "vl")
g2 <- sample_degseq(dec,method = "vl")
par(mfrow=(c(1,2)))
plot(g1)
plot(g2)
isomorphic(g1,g2)

随机图的意义在于，将随机图作为一个评价标准，从而找出图的显著特征。比如，用已知图的度分布制作一个随机图，然后进行比较：

> fyeast <- sample_degseq(degree(yeast))
> all(degree(fyeast)==degree(yeast))
[1] TRUE
> mean_distance(fyeast)
[1] 3.533058
> mean_distance(yeast)
[1] 5.095629
> diameter(fyeast)
[1] 8
> diameter(yeast)
[1] 15
> table(sapply(decompose(fyeast),vcount))

   2 2591 
  13    1 
> table(sapply(decompose(yeast),vcount))

   2    3    4    5    6    7 2375 
  63   13    5    6    1    3    1 
> transitivity(fyeast)
[1] 0.03826303
> transitivity(yeast)
[1] 0.4686178

5.4.1 小世界模型

用固定概率生成边，是最简单的一个随机模型，若对边的生成过程进行规定，总结出边的生成规律，就是基于机制的网络图
随机图模型正确模拟真实世界网络中的平均路径长度，但是低估了聚类系数。
为解决这个问题，Duncan J. Watts和Steven Strogatz 于1998年提出了小世界模型。
该网络拥有较高的聚类系数，同时有较短的平均路径长度。
为了创建具有这两个属性的网络图，Watts 和 Strogatz 建议从具有格结构lattice structure的图开始，然后随机“重新布线”一小部分边。更具体地说，在这个模型中，我们从一组 $N_v$个顶点开始，以周期性的方式排列，并将每个顶点连接到其每边的邻居的 r 个。然后，对于每条边，以概率 p 独立地移动该边的一端以入射到另一个顶点，在那里统一选择新顶点，但要注意避免循环和多边的构造。
最初的小世界模型将所有顶点放在一个圆上，并将每个顶点连接到平面距离最近的顶点。如果连接的邻居数超过 2，则会出现三元组，因为每个顶点都链接到其邻居及其邻居的邻居 。在这里，隐含的行为假设是行为人倾向于与邻居的邻居有联系。用技术术语来说，他们更喜欢传递闭合。实质上，传递性可能是由于社会参与者在地理上彼此靠近（我认识住在我隔壁的人，他们认识他们的邻居，我也认识他们）或组织上（我认识人和他们在我的组织中认识的人），或者在另一种意义上，例如，因为行为人有共同的兴趣（有些人被联系在一起是因为他们分享网络分析）。换句话说，传递性隐含着社会纽带的形成和保持很重要的语境。
小世界随机图的统计模型固定了每个顶点所链接的环上的顶点数量和附近邻居的数量。这有效地固定了原始小世界模型中的边数。模型的概率部分只涉及边的重新布边，随机选择一条边，以及随机选择一条边重新连接到的顶点。每条边和每个顶点都具有相同的被选中概率，因此该过程可以解释为两个伯努利过程，其中一个参数设置一条边被重新布边的概率。请注意，选择一个顶点来接收重新布边的边的概率由网络的大小固定，因为所有顶点具有相等的概率。
在小世界随机图中，顶点之间的平均路径距离很短，但很难说有多短。已知平均路径距离会随着顶点数量的增加呈对数增加对数增加，但如果我们只有一个网络，这并不能告诉我们太多。此外，聚类系数相对较高，因此如果我们想量化社会网络的小世界特征，我们可以将平均路径长度除以聚类系数平均路径长度除以聚类系数。该值越低，网络包含的平均路径长度越低，聚类越多。

sample_smallworld(dim, size, nei, p, loops = FALSE, multiple = FALSE)

Arguments
dim	  Integer constant, the dimension of the starting lattice.
size  Integer constant, the size of the lattice along each dimension.
nei   Integer constant, the neighborhood within which the vertices of the lattice will be connected.
p     Real constant between zero and one, the rewiring probability.
sample_smallworld(1,25,5,0.05) %>% plot()

生成的图， $N_v=25$ 即有 25 个顶点，$r=5$即每个节点的邻居数是5，重连概率 $p=0.05$
再看小世界图的特性，先把重连概率p设为0，此时得到的图是一个内聚性很高的网格图，但节点间距离很大，还不能算小世界网络：

> gws <- sample_smallworld(1,100,5,0)
> la <- layout.kamada.kawai(gws)
> plot(gws,layout=la,vertex.size=1)> diameter(gws)
[1] 10
> mean_distance(gws)
[1] 5.454545
> transitivity(gws)
[1] 0.6666667

以随机方式重连少量的边，明显减少顶点之间的距离，同时仍保持类似的高水平聚类

> gws <- sample_smallworld(1,100,5,0.05)
> la <- layout.kamada.kawai(gws)
> plot(gws,layout=la,vertex.size=1)
> diameter(gws)
[1] 5
> mean_distance(gws)
[1] 2.709091
> transitivity(gws)
[1] 0.502424

同样布局，图的样子是这样的：

用重复实验验证一下小世界模型的特性

# 设置小世界函数中的p参数值，不直接设，用10的指数来实现
> steps <- seq(-4,-0.5,0.1)
# 定义两个向量，对应36个概率参数，向量长度设为36
> (len <- length(steps))
[1] 36
> (cl <- numeric(len))
 [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> (apl <- numeric(len))
 [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# 定义内层循环的次数
> n <- 100
# 外层循环遍历所有p
> for (i in 1:len) {
+   cltmp <- numeric(n)
+   apltmp <- numeric(n)
# 在每个p中，都实验100次，因为 sample_smallworld(）是随机图，每次结果不同
+   for (j in 1:n) {
+     g <- sample_smallworld(1,1000,10,10^steps[i])
+     cltmp[j] <- transitivity(g)
+     apltmp[j] <- mean_distance(g)
+   }
# 固定p实验100次后，计算两个指标的均值
+   cl[i] <- mean(cltmp)
+   apl[i] <- mean(apltmp)
+ }
# 出图
> plot(steps,cl/max(cl),col="red")
> lines(steps,apl/max(apl),col="blue")

图中绘制了平均路径长度和传递性聚类系数的归一化版本的近似期望值，表明在 p 的很大范围内，网络表现出较小的平均距离（蓝线），同时保持高水平的内聚性（红色圆圈）。
在pajek中，Network> Create Random Network 菜单中的 Small World 命令根据 Watts-Strogatz 小世界模型创建一个无向随机图。它比IGRAPH少一个参数在对话框中，首先输入随机图中的顶点数，然后指定顶点在每一侧链接到的邻居数——例如，在此对话框中输入 3 会创建一个顶点链接到的环每边三个最近的邻居，所以平均顶点度数为
6。最后，输入 0 到 1 之间的概率，一条边将被重新布边，即随机接收一个新端点。重新布边/添加概率的用较低值-如0.01 到 0.10 - 足以在随机图中获得
相对较短的路径。

5.4.2 优先连接模型 Preferential Attachment Models

许多网络随时间增长或以其他方式演化的。
伯努利模型和小世界模型都存在一个问题：度数分布不像许多社会网络那样向右倾斜。社会网络通常包含很少的度非常高的顶点以及许多度低的顶点。优
先连接模型通过简单地假设顶点更喜欢链接到具有更高度的顶点来解决这个问题。
优先连接的统计模型通过重复次数和每次重复添加的新边数重复次数和每次重复添加的新边数直接或间接固定随机图的顶点和边数，假设每次重复增加一个新顶点。概率部分涉及为新边选择顶点。
在这个领域，通常基于节点偏好vertex preference、适应度 fitness、复制copying、老化age等概念，为网络在任何给定时间点如何变化指定一个简单的机制。
这类机理中最有名的是优先链接 Preferential Attachment是一类过程，在这类过程中，某些量（通常是某种形式的财富或信贷）是根据一些人或事物已经拥有的量来分配的，从而使那些已经富有的人比那些不富有的人得到更多。”优先链接”是描述该过程的众多名称中最贴近其本质含义的名称。它还被称为“Yule过程”、“优势积累”、“富人越来越富” ，以及说得不那么确切的“马太效应”。它们也与吉布拉定律有关。优先链接之所以受到科学家的关注，主要是因为它能在适当的条件下产生幂律分布 Power Law Distributions。
具体算法是：
以 $N_v^{(0)}$ 个顶点和$N_e^{(0)}$ 条边的初始图 $G^{(0)}$ 开始。然后，在阶段 t = 1, 2, …，对前图 $G^{(t-1)}$ 通过添加度数 $m\ge 1$ 的新顶点来创建新图 $G^{(t)}$，其中 m 条新边会与 $G^{(t-1)}$中的 m 个不同节点连接，并且新顶点与某个给定顶点 v 的连接概率由下式给出
$\frac{d_v}{{\sum }_{v^{\prime }\in V}{d}_{v^{\prime }}}$
即，在每个阶段，新节点，从高度数优先的方式，连接到m 个现有节点。在 t 次迭代之后，生成的图 $G^{(t)}$ 将具有 $N_v^{(t)}=N_v^{(0)}+t$ 顶点和 $N_e^{(t)}=N_e^{(0)}+tm$ 条边。而且，由于倾向于优先连接的趋势，直觉上我们会期望随着 t 的增加，会逐渐出现一些度值较高的节点。

> set.seed(42)
> gpa <- sample_pa(100,directed = FALSE)
> plot(gpa,layout=layout.circle,vertex.label=NA,vertex.size=4)

检查一下度分布：

> hist(degree(gpa),col = "lightblue")
> summary(degree(gpa))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00    1.00    1.00    1.98    2.00    9.00

验证一下优先连接模型的属性

> gba <- sample_pa(100,directed = FALSE)
> mean_distance(gba)
[1] 5.208081
> diameter(gba)
[1] 12
> transitivity(gba)
[1] 0

5.5.1 评估网络中的社团数量

> data("karate")
> cv <- vcount(karate)
> ev <- ecount(karate)
> dg <- degree(karate)
> n <- 1000
> ng <- numeric(n)
> for (i in 1:n) {
+   grg <- sample_gnm(cv,ev)
+   crg <- cluster_fast_greedy(grg)
+   ng[i] <- length(crg)
+ }
> gng <- numeric(n)
> for (i in 1:n) {
+   grg <- sample_degseq(dg,method = "vl")
+   crg <- cluster_fast_greedy(grg)
+   gng[i] <- length(crg)
+ }
> rl <- c(ng,gng)
> indx <- c(rep(0,n),rep(1,n))
> table(indx,rl)/n
    rl
indx     3     4     5     6     7     8
   0 0.010 0.213 0.527 0.214 0.034 0.002
   1 0.010 0.232 0.559 0.184 0.015 0.000
> counts <- table(indx,rl)/n
> barplot(counts,beside = TRUE,col=c("red","blue"))

5.5.2 评估小世界属性

> library(igraphdata)
> data("macaque")
> summary(macaque)
IGRAPH f7130f3 DN-- 45 463 -- 
+ attr: Citation (g/c), Author (g/c), shape (v/c), name (v/c)
> library(igraph)
> ccd <- function(g){
+   a <- as.matrix(get.adjacency(g))
+   s <- a+t(a)
+   dg <- degree(g,mode = c("total"))
+   nm <- diag(s %*% s %*% s)
+   dem <- diag(a %*% a)
+   dem <- 2*(dg * (dg-1) -2*dem)
+   cl <- mean(nm /dem)
+   return(cl)
+ }
> n <- 1000
> nc <- vcount(macaque)
> ne <- ecount(macaque)
> clrg <- numeric(n)
> aplrg <- numeric(n)
> for (i in 1:n) {
+   grg <- sample_gnm(nc,ne,directed = TRUE)
+   clrg <- ccd(grg)
+   aplrg <- mean_distance(grg)
+ }
> summary(clrg)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.2323  0.2323  0.2323  0.2323  0.2323  0.2323 
> summary(aplrg)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.838   1.838   1.838   1.838   1.838   1.838 
> ccd(macaque)
[1] 0.5501073
> mean_distance(macaque)
[1] 2.148485