网络数据统计分析笔记|| 网络拓扑结构推断

前情回顾：

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

我们永远不能保证我们收集到的数据是完全的，在网络图中也是一样。我们可能只有网络图中部分潜在边的状态，而不是全部信息。此外，我们甚至无法直接确定是否存在变得信息，而只能对节点和边的属性进行测量，认为他们某种程度上可以预测边的存在状态。这时，从可得数据中构建一个网络图的任务可以自然地视作一种统计推断。

我们的目标是基于数据中的信息和其他先验信息，使用统计建模和推断技术从网络集合中选出一个合适的成员，能够最好地代表系统的基本状态。目前主要有三类网络拓扑推断（Network topology inference）问题：

链路预测：已知所有节点和部分节点，推断某边是否存在
关联网络推断：知道节点的信息，不知道任何边的信息
层析拓扑结构推断：只知道部分节点的信息

链路预测

链路预测(Link prediction)中，我们希望使用观察到的边的存在状态的子集，以及可能有的节点属性的信息，预测网络中节点对之间潜在的边是否存在。如传染病的传播路径，我们可以推断出来病毒是如何以及可能的传播路径。

我们还是以法国政客博客网络fblog为例来说明问题，网络中没对节点共同邻居数量可以按照如下方式计算。

library(sand)
nv <- vcount(fblog)
ncn <- numeric()
A <- as_adjacency_matrix(fblog)

for(i in (1:(nv-1))){
  ni <- ego(fblog, 1, i) # Neighborhood of graph vertices
  nj <- ego(fblog, 1, (i+1):nv)
  nbhd.ij <- mapply(intersect, ni, nj, SIMPLIFY=FALSE)
  temp <- unlist(lapply(nbhd.ij, length)) - 
    2*A[i, (i+1):nv]
  ncn <- c(ncn, temp) 
}

# CHUNK 2
# install.packages('vioplot')
library(vioplot)
Avec <- A[lower.tri(A)]
vioplot(ncn[Avec==0], ncn[Avec==1], 
   names=c("No Edge", "Edge"),
   col="magenta")
title(ylab="Number of Common Neighbors")

library(ROCR)
pred <- prediction(ncn, Avec)
perf <- performance(pred, "auc")
slot(perf, "y.values")
# ---
## [[1]]
## [1] 0.9275179
# ---

关联网络推断（Association networks inference ）

当我们将某个数据集表示为网络时，定义边的规则通常是两个相邻节点的某些属性具有足够程度的关联。这种关联网络可以再很多领域中看到。当存在采样误差的时候，有必要在构建关联网络的过程中纳入统计原理和方法。

举一个例子

我们还是用Ecoli.data数据来做演示，其中有153个基因的表达量，另一个是一个邻接矩阵，概括了我们对于Ecoli中实际调控关系（不完全）的认识。该数据是在手机实验数据同期从RegulonDB数据库收集的。

rm(list=ls())
data(Ecoli.data)
ls()
# ---
## [1] "Ecoli.expr" "regDB.adj"
# ---
Ecoli.expr[1:5,1:5]
              acrR_b0464_at ada_b2213_at adiY_b4116_at alpA_b2624_at appY_b0564_at
yebF__U_N0075      7.439417     8.079050      7.692263      7.569850      4.701593
zipA__U_N0075      7.786913     8.236573      7.611393      7.463183      4.650187
b2618_U_N0075      7.770360     8.040610      8.457703      7.287090      4.701950
bcp___U_N0075      7.648287     8.041520      8.488933      7.463150      4.672690
cpxR__U_N0075      7.713530     8.061450      8.285527      7.626580      4.665150

 regDB.adj[1:5,1:5]
              acrR_b0464_at ada_b2213_at adiY_b4116_at alpA_b2624_at appY_b0564_at
acrR_b0464_at             0            0             0             0             0
ada_b2213_at              0            0             0             0             0
adiY_b4116_at             0            0             0             0             0
alpA_b2624_at             0            0             0             0             0
appY_b0564_at             0            0             0             0             0

heatmap(scale(Ecoli.expr), Rowv=NA)

可以看到基因表达有一些局部的特点，但是不是很清晰。

library(igraph)
g.regDB <- graph_from_adjacency_matrix(regDB.adj,
                                       "undirected")
summary(g.regDB)

IGRAPH 8169821 UN-- 153 209 -- 
+ attr: name (v/c)

plot(g.regDB, vertex.size=3, vertex.label=NA)

构建相关网络。

mycorr <- cor(Ecoli.expr)

# CHUNK 9
z <- 0.5 * log((1 + mycorr) / (1 - mycorr))

# CHUNK 10
z.vec <- z[upper.tri(z)]
n <- dim(Ecoli.expr)[1]
corr.pvals <- 2 * pnorm(abs(z.vec), 0,
   sqrt(1 / (n-3)), lower.tail=FALSE)

# CHUNK 11
length(corr.pvals) # 进行了多少次检验
[1] 11628

多重检验检验中p值的校正。

corr.pvals.adj <- p.adjust(corr.pvals, "BH")

# CHUNK 13
length(corr.pvals.adj[corr.pvals.adj < 0.05])
 5227

有 5227对基因是显著的，这数值远远大于真实结果，这个结果是可疑的。幸运的是，之前有1万多个独立的元素，在这种明显是数据丰富的情况下，我们可以使用基于数据的方法来获得零分布。

library(fdrtool)

# CHUNK 15
mycorr.vec <- mycorr[upper.tri(mycorr)]
fdr <- fdrtool(mycorr.vec, statistic="correlation")
Step 1... determine cutoff point
Step 2... estimate parameters of null distribution and eta0
Step 3... compute p-values and estimate empirical PDF/CDF
Step 4... compute q-values and local fdr
Step 5... preppare for plotting

偏相关网络

# CHUNK 16
pcorr.pvals <- matrix(0, dim(mycorr)[1],
    dim(mycorr)[2])
for(i in seq(1, 153)){
   for(j in seq(1, 153)){
     rowi <- mycorr[i, -c(i, j)]
     rowj <- mycorr[j, -c(i, j)]
     tmp <- (mycorr[i, j] -
       rowi*rowj)/sqrt((1-rowi^2) * (1-rowj^2))
     tmp.zvals <- (0.5) * log((1+tmp) / (1-tmp))
     tmp.s.zvals <- sqrt(n-4) * tmp.zvals
     tmp.pvals <- 2 * pnorm(abs(tmp.s.zvals),
       0, 1, lower.tail=FALSE)
     pcorr.pvals[i, j] <- max(tmp.pvals)
   }
}

# CHUNK 17
pcorr.pvals.vec <- pcorr.pvals[lower.tri(pcorr.pvals)]
pcorr.pvals.adj <- p.adjust(pcorr.pvals.vec, "BH")

# CHUNK 18
pcorr.edges <- (pcorr.pvals.adj < 0.05)
length(pcorr.pvals.adj[pcorr.edges])

25

> pcorr.A <- matrix(0, 153, 153)
> pcorr.A[lower.tri(pcorr.A)] <- as.numeric(pcorr.edges)
> g.pcorr <- graph_from_adjacency_matrix(pcorr.A,
+                                         "undirected")
> # CHUNK 20
> intersection(g.regDB, g.pcorr, byname=FALSE)
IGRAPH 51ace56 UN-- 153 4 -- 
+ attr: name (v/c)
+ edges from 51ace56 (vertex names):
[1] yhiW_b3515_at--yhiX_b3516_at rhaR_b3906_at--rhaS_b3905_at marA_b1531_at--marR_b1530_at gutM_b2706_at--srlR_b2707_at
> # CHUNK 21
> fdr <- fdrtool(pcorr.pvals.vec, statistic="pvalue",
+    plot=FALSE)
Step 1... determine cutoff point
Step 2... estimate parameters of null distribution and eta0
Step 3... compute p-values and estimate empirical PDF/CDF
Step 4... compute q-values and local fdr

> pcorr.edges.2 <- (fdr$qval < 0.05)
> length(fdr$qval[pcorr.edges.2])
[1] 25

高斯图网络

library(huge)
set.seed(42)
huge.out <- huge(Ecoli.expr)

# CHUNK 23
huge.opt <- huge.select(huge.out, criterion="ric")
g.huge <- graph_from_adjacency_matrix(huge.opt$refit, 
                                       "undirected")
ecount(g.huge)
# ---
## [1] 0
# ---

# CHUNK 24
huge.opt <- huge.select(huge.out, criterion="stars")
g.huge <- graph_from_adjacency_matrix(huge.opt$refit, "undirected")
summary(g.huge)
# ---
## IGRAPH 0f13a7e U--- 153 759 --
# ---

# CHUNK 25
intersection(g.pcorr, g.huge)
# ---
## IGRAPH 6948e99 U--- 153 25 -- 
## + edges from 6948e99:
## [1] 145--146 144--146 112--125 112--113 109--138 
## [6] 108--135  97--111  96--119  92--107  87--104  
## [11]  86-- 87  84--129  81--137  72--141  70-- 75  
## [16]  60--127  46-- 77  42-- 43  38--153  37-- 98  
## [21]  27--123  21-- 33  12--135   9-- 90   3-- 60
# ---

# CHUNK 26
intersection(g.regDB, g.huge, byname=FALSE)
# ---
## IGRAPH 662d795 UN-- 153 21 -- 
## + attr: name (v/c)
## + edges from 662d795 (vertex names):
## [1] yhiW_b3515_at--yhiX_b3516_at
## [2] tdcA_b3118_at--tdcR_b3119_at
## [3] rpoE_b2573_at--rpoH_b3461_at
## [4] rpoD_b3067_at--tyrR_b1323_at
## [5] rhaR_b3906_at--rhaS_b3905_at
## [6] nac_b1988_at --putA_b1014_at
## [7] marA_b1531_at--marR_b1530_at
## [8] malT_b3418_at--rpoD_b3067_at
## + ... omitted several edges
# ---

层析拓扑结构推断(Tomographic network topology inference)

网络层析技术是一种基于网络端到端(End-to-End)的性能测量统计，利用信号处理及统计分析的方法推测网络拓扑结构及内部链路级的丢包率、时延、瓶颈带宽等性能参数，对无法直接观测的网络内部性能进行推断的方法。它具有不依赖路由器协作、开销少的特点，保证了用户隐私和信息安全，减少了测量信息的传输，对于网络测量具有重要意义。传统网络拓扑推断方法的缺陷是要求大量的中间节点配合以及相应的协议支持，这在实际的大型网络中往往是难以办到的。与传统方法相比，基于网络层析的拓扑推断方法只通过对网络端到端的测量来推测网络拓扑结构，其最大的优势在于可以不需要中间节点的协作就能完成网络拓扑结构推断，具有很高的研究价值。

![](https://upload-images.jianshu.io/upload_images/7600498-b22a1360c6766fe5.png?imageMogr2
/auto-orient/strip%7CimageView2/2/w/1240)

data(sandwichprobe)
delaydata[1:5, ]
# ---
##   DelayDiff SmallPktDest BigPktDest
## 1       757            3         10
## 2       608            6          2
## 3       242            8          9
## 4        84            1          8
## 5      1000            7          3
# ---
host.locs
# ---
##  [1] "IST"    "IT"     "UCBkly" "MSU1"   "MSU2"
##  [6] "UIUC"   "UW1"    "UW2"    "Rice1"  "Rice2"
# ---

# CHUNK 28
meanmat <- with(delaydata, by(DelayDiff,
   list(SmallPktDest, BigPktDest), mean))
image(log(meanmat + t(meanmat)), xaxt="n", yaxt="n",
   col=heat.colors(16))
mtext(side=1, text=host.locs,
   at=seq(0.0,1.0,0.11), las=3)
mtext(side=2, text=host.locs,
   at=seq(0.0,1.0,0.11), las=1)

# CHUNK 29
SSDelayDiff <- with(delaydata, by(DelayDiff^2,
   list(SmallPktDest, BigPktDest), sum))

# CHUNK 30
x <- as.dist(1 / sqrt(SSDelayDiff))
myclust <- hclust(x, method="average")

# CHUNK 31
plot(myclust, labels=host.locs, axes=FALSE,
   ylab=NULL, ann=FALSE)

https://cloud.tencent.com/developer/article/1667098
https://www.bbsmax.com/A/gVdnD9rXJW/
https://blog.csdn.net/qq_36384657/article/details/108095179
https://ieeexplore.ieee.org/document/7990048
http://www2.ece.rochester.edu/~gmateosb/ECE442/Slides/block_4_sampling_modeling_inference_part_d.pdf
https://dspace.xmu.edu.cn/bitstream/handle/2288/51403/%e5%9f%ba%e4%ba%8e%e7%bd%91%e7%bb%9c%e5%b1%82%e6%9e%90%e7%9a%84%e6%8b%93%e6%89%91%e6%8e%a8%e6%96%ad%e7%ae%97%e6%b3%95%e7%a0%94%e7%a9%b6.pdf?sequence=1&isAllowed=y