前情回顧：
Gephi網(wǎng)絡(luò)圖極簡教程
Network在單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析中的應(yīng)用
Gephi網(wǎng)絡(luò)圖極簡教程
Network在單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析中的應(yīng)用
網(wǎng)絡(luò)數(shù)據(jù)統(tǒng)計(jì)分析筆記|| 為什么研究網(wǎng)絡(luò)
網(wǎng)絡(luò)數(shù)據(jù)統(tǒng)計(jì)分析筆記|| 操作網(wǎng)絡(luò)數(shù)據(jù)
網(wǎng)絡(luò)數(shù)據(jù)統(tǒng)計(jì)分析筆記|| 網(wǎng)絡(luò)數(shù)據(jù)可視化
網(wǎng)絡(luò)數(shù)據(jù)統(tǒng)計(jì)分析筆記|| 網(wǎng)絡(luò)數(shù)據(jù)的描述性分析
網(wǎng)絡(luò)數(shù)據(jù)統(tǒng)計(jì)分析筆記||網(wǎng)絡(luò)圖的數(shù)學(xué)模型
學(xué)過統(tǒng)計(jì)的我們知道，要執(zhí)行推斷只靠幾個(gè)分布模型是不夠的，還要有模型來做推斷，這就是統(tǒng)計(jì)建模。那么，上一章的經(jīng)典隨機(jī)圖模型，伯努利隨機(jī)圖模型，廣義隨機(jī)圖模型等是我們網(wǎng)絡(luò)圖世界的正態(tài)分布，卡方分布等等等。本章介紹的內(nèi)容，作類比的話，就像之前我們學(xué)過的線性模型，非線性模型，廣義線性模型之流。
為什么要學(xué)習(xí)模型，一個(gè)模型就像一個(gè)思考框架，在我們需要描述問題的時(shí)候，有個(gè)思考的方向。模型，質(zhì)言之，有模具的形狀。之所有模型，是因?yàn)檫@世界是可歸納的。
當(dāng)前主要的三類模型與用于非網(wǎng)絡(luò)數(shù)據(jù)的統(tǒng)計(jì)模型很相似。指數(shù)隨機(jī)圖模型類似標(biāo)準(zhǔn)的回歸模型，尤其是廣義線性模型。類似地，隨機(jī)塊模型受到混合模型的啟發(fā)，其基本形式本質(zhì)上是一些經(jīng)典隨機(jī)圖模型的混合。最后，潛變量網(wǎng)絡(luò)模型是一種基于網(wǎng)絡(luò)的變體模型，同時(shí)使用觀測變量和未觀測變量（即，潛變量）對(duì)結(jié)果（此處為網(wǎng)絡(luò)中邊的缺失與否）進(jìn)行建模。
- 指數(shù)隨機(jī)圖模型
- 一般形式
- 模型界定
- 模型擬合
- 擬合優(yōu)度
- 網(wǎng)絡(luò)塊模型
- 模型界定
- 模型擬合
- 潛變量模型
- 一般形式
- 模型界定
- 模型擬合
- 擬合優(yōu)度
指數(shù)隨機(jī)圖模型
指數(shù)隨機(jī)圖模型（exponential random graph models，ERGM）直接類比經(jīng)典的廣義線性模型（GLM）而建立。這當(dāng)然是為了在模型建立、擬合與比較方面利用和擴(kuò)展業(yè)已成熟的統(tǒng)計(jì)原理與方法。然而，指數(shù)隨機(jī)模型的界定和擬合顯然比標(biāo)準(zhǔn)的廣義線性模型更為微妙。在我們使用模型的時(shí)候，一定要知道自己使用的是什么。應(yīng)用模型或者算法，就像吃飯，要知道自己吃的是什么。應(yīng)用模型不可以離開具體的應(yīng)用場景，不然代碼會(huì)顯得非常簡單乏味枯燥。
一般形式
考慮一個(gè)隨機(jī)圖
其中：
- 每個(gè)都是一個(gè)構(gòu)型，其定義為的一個(gè)節(jié)點(diǎn)子集節(jié)點(diǎn)之間可能的變得集合；
-
,故，若構(gòu)型存在于則為1 ，否則為0.
- 非零值表示在給定剩余部分圖的條件下，與中的所有節(jié)點(diǎn)對(duì)相依。
-
是歸一化常數(shù)。
我們可以從一個(gè)
在R中我們使用擴(kuò)展包是ergm，它是statnet的一部分。由于ergm使用network的網(wǎng)絡(luò)對(duì)象，為了將igraph對(duì)象轉(zhuǎn)化為statnet的格式，首先將網(wǎng)絡(luò)轉(zhuǎn)化為鄰接矩陣和相關(guān)屬性。
library(sand)
data(lazega)
A <- as_adjacency_matrix(lazega)
v.attrs <- as_data_frame(lazega, what="vertices")
之后創(chuàng)建ergm對(duì)象
library(ergm) # Will load package 'network' as well.
lazega.s <- network::as.network(as.matrix(A),
directed=FALSE)
network::set.vertex.attribute(lazega.s, "Office",
v.attrs$Office)
network::set.vertex.attribute(lazega.s, "Practice",
v.attrs$Practice)
network::set.vertex.attribute(lazega.s, "Gender",
v.attrs$Gender)
network::set.vertex.attribute(lazega.s, "Seniority",
v.attrs$Seniority)
network::plot.network(lazega.s,vertex.col = "Seniority")
模型界定
模型界定是在給定模型框架下，進(jìn)一步具體化相關(guān)函數(shù)，建立更加具體的模型。
my.ergm.bern <- formula(lazega.s ~ edges)
my.ergm.bern
lazega.s ~ edges
summary(my.ergm.bern)
edges
115
my.ergm <- formula(lazega.s ~ edges + kstar(2)
+ kstar(3) + triangle)
summary(my.ergm)
edges kstar2 kstar3 triangle
115 926 2681 120
my.ergm <- formula(lazega.s ~ edges
+ gwesp(1, fixed=TRUE))
summary(my.ergm)
edges gwesp.fixed.1
115.0000 213.1753
lazega.ergm <- formula(lazega.s ~ edges
+ gwesp(log(3), fixed=TRUE)
+ nodemain("Seniority")
+ nodemain("Practice")
+ match("Practice")
+ match("Gender")
+ match("Office"))
這個(gè)公式像不像我們線性模型里面的R軟件提供的擬合計(jì)算廣義線性模型的函數(shù)glm()
glm(formula, family = gaussian, data, weights, subset,
na.action, start = NULL, etastart, mustart, offset,
control = list(...), model = TRUE, method = "glm.fit",
x = FALSE, y = TRUE, singular.ok = TRUE, contrasts = NULL, ...)
glm.fit(x, y, weights = rep.int(1, nobs),
start = NULL, etastart = NULL, mustart = NULL,
offset = rep.int(0, nobs), family = gaussian(),
control = list(), intercept = TRUE, singular.ok = TRUE)
模型擬合
set.seed(42)
lazega.ergm.fit <- ergm(lazega.ergm)
Starting maximum pseudolikelihood estimation (MPLE):
Evaluating the predictor and response matrix.
Maximizing the pseudolikelihood.
Finished MPLE.
Starting Monte Carlo maximum likelihood estimation (MCMLE):
Iteration 1 of at most 20:
Optimizing with step length 1.
The log-likelihood improved by 0.5253.
Step length converged once. Increasing MCMC sample size.
Iteration 2 of at most 20:
Optimizing with step length 1.
The log-likelihood improved by 0.07175.
Step length converged twice. Stopping.
Finished MCMLE.
Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
This model was fit using MCMC. To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.
Warning message:
`set_attrs()` is deprecated as of rlang 0.3.0
This warning is displayed once per session.
方差分析
anova(lazega.ergm.fit)
Analysis of Variance Table
Model 1: lazega.s ~ edges + gwesp(log(3), fixed = TRUE) + nodemain("Seniority") +
nodemain("Practice") + match("Practice") + match("Gender") +
match("Office")
Df Deviance Resid. Df Resid. Dev Pr(>|Chisq|)
NULL 630 873.37
Model 1: 7 413.74 623 459.63 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
在圖模型的框架里面也要考慮哪些變量納入到模型中，這時(shí)候就需要看看：
變量貢獻(xiàn)度
summary(lazega.ergm.fit)
==========================
Summary of model fit
==========================
Formula: lazega.s ~ edges + gwesp(log(3), fixed = TRUE) + nodemain("Seniority") +
nodemain("Practice") + match("Practice") + match("Gender") +
match("Office")
Iterations: 2 out of 20
Monte Carlo MLE Results:
Estimate Std. Error MCMC % z value Pr(>|z|)
edges -7.00655 0.67114 0 -10.440 < 1e-04 ***
gwesp.fixed.1.09861228866811 0.59166 0.08554 0 6.917 < 1e-04 ***
nodecov.Seniority 0.02456 0.00620 0 3.962 < 1e-04 ***
nodecov.Practice 0.39455 0.10218 0 3.861 0.000113 ***
nodematch.Practice 0.76966 0.19060 0 4.038 < 1e-04 ***
nodematch.Gender 0.73767 0.24362 0 3.028 0.002463 **
nodematch.Office 1.16439 0.18753 0 6.209 < 1e-04 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Null Deviance: 873.4 on 630 degrees of freedom
Residual Deviance: 459.6 on 623 degrees of freedom
AIC: 473.6 BIC: 504.7 (Smaller is better.)
擬合優(yōu)度
> gof.lazega.ergm <- gof(lazega.ergm.fit)
> # CHUNK 12
> plot(gof.lazega.ergm)
> gof.lazega.ergm
Goodness-of-fit for degree
obs min mean max MC p-value
0 2 0 2.61 8 1.00
1 3 0 2.81 7 1.00
2 2 0 2.67 8 1.00
3 4 0 2.94 9 0.72
4 2 0 2.95 8 0.92
5 4 0 2.73 7 0.62
6 4 0 2.93 8 0.70
7 1 0 3.11 7 0.36
8 1 0 2.65 7 0.46
9 5 0 2.35 7 0.14
10 1 0 1.97 10 0.90
11 1 0 1.80 8 0.94
12 2 0 1.39 7 0.68
13 3 0 1.11 5 0.26
14 0 0 0.69 4 0.92
15 1 0 0.60 4 0.84
16 0 0 0.30 2 1.00
17 0 0 0.24 2 1.00
18 0 0 0.06 1 1.00
19 0 0 0.06 1 1.00
20 0 0 0.02 1 1.00
23 0 0 0.01 1 1.00
Goodness-of-fit for edgewise shared partner
obs min mean max MC p-value
esp0 5 2 8.52 21 0.46
esp1 16 6 15.29 27 0.98
esp2 29 6 20.86 38 0.20
esp3 17 3 21.81 44 0.50
esp4 23 2 18.66 35 0.52
esp5 11 0 12.72 32 0.94
esp6 10 0 7.46 26 0.54
esp7 4 0 4.04 18 0.96
esp8 0 0 1.95 8 0.62
esp9 0 0 0.96 12 1.00
esp10 0 0 0.19 3 1.00
esp11 0 0 0.10 2 1.00
Goodness-of-fit for minimum geodesic distance
obs min mean max MC p-value
1 115 60 112.56 163 1.00
2 275 100 259.30 354 0.96
3 148 45 128.59 198 0.56
4 21 0 26.73 120 1.00
5 2 0 5.06 63 0.76
6 0 0 1.00 27 1.00
7 0 0 0.29 21 1.00
8 0 0 0.06 6 1.00
Inf 69 0 96.41 304 0.96
Goodness-of-fit for model statistics
obs min mean max MC p-value
edges 115.0000 60.00000 112.5600 163.0000 1.00
gwesp.fixed.1.09861228866811 222.9108 80.38272 215.2664 383.0266 0.96
nodecov.Seniority 4687.0000 2494.00000 4552.9400 6505.0000 0.82
nodecov.Practice 359.0000 181.00000 351.0100 506.0000 0.98
nodematch.Practice 72.0000 40.00000 70.7500 107.0000 1.00
nodematch.Gender 99.0000 48.00000 97.4500 140.0000 0.94
nodematch.Office 85.0000 45.00000 83.5000 133.0000 0.98
網(wǎng)絡(luò)塊模型
塊模型是一類帶有種植簇的隨機(jī)圖?！澳改Ｐ汀笔请S機(jī)塊模型(random block model, SBM)，它被廣泛用作社區(qū)檢測(community detection)的典型模型。它可以說是帶有社區(qū)的圖的最簡單模型。由于SBM是一種生成模型，它從社區(qū)的基礎(chǔ)事實(shí)中獲益，這允許在正式的環(huán)境中考慮背景問題。
SBM的歷史很長。Asmentioned早些時(shí)候,模型出現(xiàn)在多個(gè)獨(dú)立科學(xué)同一片藍(lán)天下：座的術(shù)語。近年來似乎已經(jīng)占主導(dǎo)地位的,來自于機(jī)器學(xué)習(xí)和統(tǒng)計(jì)文獻(xiàn)，雖然這些模型通常被稱為種植分區(qū)模型（planted partition model ）在理論電腦科學(xué)[BCLS87、DF89 Bop87],和非均勻隨機(jī)圖模型在數(shù)學(xué)文獻(xiàn)[BJR07]
在某種意義上，SBM的作用與信息論中的離散無記憶通道(DMC)相似。雖然建模外部噪音的任務(wù)可能更容易簡化真實(shí)數(shù)據(jù)集，SBM捕獲了社區(qū)檢測的一些關(guān)鍵瓶頸現(xiàn)象，并承認(rèn)許多可能的改進(jìn)，以提高對(duì)真實(shí)數(shù)據(jù)的適合度。在這里，我們的重點(diǎn)將放在對(duì)核心SBM的基本理解上，而不是深入到改進(jìn)后的擴(kuò)展上。
#install.packages('blockmodels')
library(blockmodels)
set.seed(42)
fblog
IGRAPH 3e87bca UN-- 192 1431 --
+ attr: name (v/c), PolParty (v/c)
+ edges from 3e87bca (vertex names):
[1] jeunesverts.org/bordeaux-- bix.enix.org/ jeunesverts.org/bordeaux-- dominiquevoynet.net/blog
[3] bix.enix.org/ -- www.arnaudcaron.net/ bix.enix.org/ -- dominiquevoynet.net/blog
[5] bix.enix.org/ -- blogs.lesverts.fr/ bix.enix.org/ -- emilien.net/
[7] bix.enix.org/ -- lipietz.net/blog.php3?id_breve=63 bix.enix.org/ -- democratiesansfrontiere.org/
[9] bix.enix.org/ -- www.dsk2007.net bix.enix.org/ -- www.parti-socialiste.fr
[11] bix.enix.org/ -- puteaux.typepad.com bix.enix.org/ -- www.monputeaux.com/
[13] bix.enix.org/ -- www.marielauremeyer.com/ bix.enix.org/ -- remibazillier.blogspot.com/
[15] bix.enix.org/ -- www.temps-reels.net/blogs/paris/ bix.enix.org/ -- www.monneuilly.com/
+ ... omitted several edges
A.fblog <- as.matrix(as_adjacency_matrix(fblog))
# 對(duì)塊模型進(jìn)行伯努利概率分布估計(jì)描述
fblog.sbm <- BM_bernoulli("SBM_sym", A.fblog,
verbosity=0, plotting='')
blockmodels object
model: bernoulli
membership: SBM_sym
network: 192 x 192 scalar network
maximum of ICL: 10 groups
Most usefull fields and methods:
The following fields are indexed by the number of groups:
$ICL : vector of ICL
$PL : vector of pseudo log liklihood
$memberships : list of memberships founds by estimation
each membership is represented object
$model_parameters : models parameters founds by estimation
Estimation methods:
$estimate(reinitalization_effort=1) : to run again estimation with a
higher reinitalization effort
Plotting methods:
$plot_obs_pred(Q) : to plot the obeserved and predicted network for Q groups
$plot_parameters(Q) : to plot the model_parameters for Q groups
Please note that each membership object have a plotting pethod
fblog.sbm$estimate()
這里用積分似然準(zhǔn)則（Integrated ClassificationLikelihood，ICL）選擇網(wǎng)絡(luò)擬合的類別數(shù)量。
塊聚類(或共聚類)旨在同時(shí)對(duì)數(shù)據(jù)表的行和列進(jìn)行分區(qū)，以揭示同構(gòu)塊結(jié)構(gòu)。這種結(jié)構(gòu)可以源于潛在塊模型，該模型提供了數(shù)據(jù)表的概率建模，其塊模式由行和列類定義。對(duì)于連續(xù)數(shù)據(jù)，每個(gè)表?xiàng)l目通常假定遵循高斯分布。對(duì)于給定的數(shù)據(jù)表，通常會(huì)檢查幾個(gè)候選模型:它們可能在集群數(shù)量或offree參數(shù)數(shù)量上有所不同。然后，模型選擇就變成了一個(gè)關(guān)鍵的問題，為基于模型的單向聚類而派生的工具需要進(jìn)行調(diào)整。在單向聚類中，大多數(shù)選擇標(biāo)準(zhǔn)是基于漸近的考慮，由于行和列的雙重性質(zhì)，這些考慮很難在塊聚類中呈現(xiàn)。我們通過基于綜合分類可能性開發(fā)非漸進(jìn)標(biāo)準(zhǔn)（non-asymptotic criterion）來規(guī)避這一問題。 在參數(shù)上定義適當(dāng)?shù)南惹胺植己?，可以以封閉形式計(jì)算此標(biāo)準(zhǔn)。 實(shí)驗(yàn)結(jié)果表明，具有分離良好和中分聚類的大中型數(shù)據(jù)表性能穩(wěn)定。
# CHUNK 14
ICLs <- fblog.sbm$ICL
Q <- which.max(ICLs)
Q
# ---
## [1] 10
# ---
(Z <- fblog.sbm$memberships[[Q]]$Z)
(cl.labs <- apply(Z,1,which.max)) # 可以看做對(duì)節(jié)點(diǎn)的劃分（聚類）
[1] 1 1 1 1 1 1 1 3 3 6 6 6 9 6 6 6 6 6 2 6 6 6 3 6 3 6 6 1 6 3 3 3 6 6 6 6 3 3 6 3 3 6 6 6 6 9 6 2 8 2 8 2 2 2 2 2 2 8 8 9 2
[62] 2 2 8 2 2 2 2 2 2 2 2 2 2 8 2 2 2 2 10 1 4 10 10 4 7 4 4 10 10 1 4 4 4 4 7 10 10 10 10 10 10 9 10 4 4 4 7 7 4 7 10 10 10 10 4 4 4 4 1 1 9
[123] 4 10 7 10 10 4 1 1 10 10 4 10 4 4 1 1 1 1 1 1 10 9 1 1 1 1 1 1 1 1 1 5 5 1 5 5 5 5 5 5 5 1 1 9 2 9 8 9 9 9 1 9 9 1 1 5 5 5 5 5 5
[184] 5 5 5 5 5 5 5 5 5
nv <- vcount(fblog)
summary(Z[cbind(1:nv,cl.labs)])
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.8586 0.9953 0.9953 0.9938 0.9953 0.9953
# CHUNK 18
cl.cnts <- as.vector(table(cl.labs))
alpha <- cl.cnts/nv
alpha
# ---
## [1] 0.18229167 0.14062500 0.05729167 0.10937500
## [5] 0.12500000 0.13020833 0.03125000 0.03645833
## [9] 0.06770833 0.11979167
# ---
# CHUNK 19
Pi.mat <- fblog.sbm$model_parameters[[Q]]$pi
Pi.mat[3,]
# ---
## [1] 0.0030340287 0.0073657690 0.9102251927 0.0009221811
## [5] 0.0009170384 0.0364593875 0.0177621832 0.0024976022
## [9] 0.0431732528 0.0012495852
# ---
# CHUNK 20
ntrials <- 1000
Pi.mat <- (t(Pi.mat)+Pi.mat)/2
deg.summ <- list(ntrials)
for(i in (1:ntrials)){
blk.sz <- rmultinom(1,nv,alpha)
g.sbm <- sample_sbm(nv,pref.matrix=Pi.mat,
block.sizes=blk.sz,
directed=FALSE)
deg.summ[[i]] <- summary(degree(g.sbm))
}
Reduce('+',deg.summ)/ntrials
# ---
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.931 9.165 13.127 15.183 18.896 49.484
# ---
summary(degree(fblog))
# ---
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.00 8.00 13.00 14.91 18.00 56.00
# ---
擬合優(yōu)度
評(píng)估隨機(jī)塊模型的擬合優(yōu)度
plot(fblog.sbm$ICL,xlab="Q",ylab="ICL",type="b")
lines(c(Q,Q),c(min(ICLs),max(ICLs)),col="red",lty=2)
edges <- as_edgelist(fblog,names=FALSE)
neworder<-order(cl.labs)
m<-t(matrix(order(neworder)[as.numeric(edges)],2))
plot(1, 1, xlim = c(0, nv + 1), ylim = c(nv + 1, 0),
type = "n", axes= FALSE, xlab="Classes",
ylab="Classes",
main="Reorganized Adjacency matrix")
rect(m[,2]-0.5,m[,1]-0.5,m[,2]+0.5,m[,1]+0.5,col=1)
rect(m[,1]-0.5,m[,2]-0.5,m[,1]+0.5,m[,2]+0.5,col=1)
cl.lim <- cl.cnts
cl.lim <- cumsum(cl.lim)[1:(length(cl.lim)-1)]+0.5
clip(0,nv+1,nv+1,0)
abline(v=c(0.5,cl.lim,nv+0.5),
h=c(0.5,cl.lim,nv+0.5),col="red")
對(duì)這10個(gè)分群（劃分）的節(jié)點(diǎn)屬性可視化。
g.cl <- graph_from_adjacency_matrix(Pi.mat,
mode="undirected",
weighted=TRUE)
# Set necessary parameters
vsize <- 100*sqrt(alpha)
ewidth <- 10*E(g.cl)$weight
PolP <- V(fblog)$PolParty
class.by.PolP <- as.matrix(table(cl.labs,PolP))
pie.vals <- lapply(1:Q, function(i)
as.vector(class.by.PolP[i,]))
my.cols <- topo.colors(length(unique(PolP)))
# Plot
plot(g.cl, edge.width=ewidth,
vertex.shape="pie", vertex.pie=pie.vals,
vertex.pie.color=list(my.cols),
vertex.size=vsize, vertex.label.dist=0.1*vsize,
vertex.label.degree=pi)
# Add a legend
my.names <- names(table(PolP))
my.names[2] <- "Comm. Anal."
my.names[5] <- "PR de G"
legend(x="topleft", my.names,
fill=my.cols, bty="n")
其中餅圖是博客對(duì)應(yīng)政黨所占的比例，邊的粗細(xì)正比于兩組博客間相互連接的估計(jì)概率。
潛變量模型
隨機(jī)塊模型及其關(guān)鍵創(chuàng)新點(diǎn)在于以節(jié)點(diǎn)類別的形式引入了潛變量（latent variable）。換言之，模型使用了未被觀測到的變量，而這些變量在決定節(jié)點(diǎn)對(duì)之間的鏈接起作用。已知了我們的數(shù)據(jù)有n個(gè)t，那么這個(gè)時(shí)候我們需要一組神奇的變量x,，他看不見摸不到，但是卻實(shí)際的決定了每個(gè)t的狀態(tài)，所以，我們把這個(gè)變量稱為潛變量，潛在的變量，看不見的變量，潛水的變量（_）。這樣，我們就可以得到一個(gè)潛變量和原始變量的聯(lián)合分布。
模型界定
如何描述潛在變量是一個(gè)概率函數(shù)，如何定義這個(gè)函數(shù)反映了應(yīng)用者希望看到哪些潛在變量，目前主要有三個(gè)函數(shù)：
- 潛類別模型（latent class）
- 潛在距離模型（latent distance）
- 特征模型( eigenmodel)
模型擬合
構(gòu)建得到的潛變量模型為分層形式，所以和自然地采用貝葉斯推斷方法。
summary(lazega)
IGRAPH 3e8b2bf UN-- 36 115 --
+ attr: name (v/c), Seniority (v/n), Status (v/n), Gender (v/n), Office (v/n), Years (v/n), Age (v/n), Practice (v/n),
| School (v/n)
print_all(lazega)
IGRAPH 3e8b2bf UN-- 36 115 --
+ attr: name (v/c), Seniority (v/n), Status (v/n), Gender (v/n), Office (v/n), Years (v/n), Age (v/n), Practice (v/n),
| School (v/n)
+ edges (vertex names):
V1 -- V17
V2 -- V7, V16, V17, V22, V26, V29
V3 -- V18, V25, V28
V4 -- V12, V17, V19, V20, V22, V26, V28, V29, V31
V5 -- V18, V24, V28, V31, V32, V33
V6 -- V24, V28, V30, V31, V32
V7 -- V2, V18
V8 --
V9 -- V12, V16, V29
V10 -- V24, V26, V29, V31, V34
V11 -- V17
V12 -- V4, V9, V15, V16, V17, V19, V26, V29, V34
V13 -- V31, V33
V14 -- V16, V17, V25, V28, V30, V32
V15 -- V12, V16, V19, V20, V22, V24, V26, V29, V32, V35, V36
V16 -- V2, V9, V12, V14, V15, V17, V22, V26, V27, V29, V32, V34, V36
V17 -- V1, V2, V4, V11, V12, V14, V16, V19, V22, V24, V25, V26, V28, V29, V34
V18 -- V3, V5, V7, V28, V31, V32, V33, V35
V19 -- V4, V12, V15, V17, V22, V24, V26, V28, V34, V35
V20 -- V4, V15, V22, V26
V21 -- V27
V22 -- V2, V4, V15, V16, V17, V19, V20, V31, V32
V23 --
V24 -- V5, V6, V10, V15, V17, V19, V26, V31, V36
V25 -- V3, V14, V17, V28, V35
V26 -- V2, V4, V10, V12, V15, V16, V17, V19, V20, V24, V27, V32
V27 -- V16, V21, V26
V28 -- V3, V4, V5, V6, V14, V17, V18, V19, V25, V30, V31, V32, V35
V29 -- V2, V4, V9, V10, V12, V15, V16, V17, V34
V30 -- V6, V14, V28, V31
V31 -- V4, V5, V6, V10, V13, V18, V22, V24, V28, V30, V32, V33, V35
V32 -- V5, V6, V14, V15, V16, V18, V22, V26, V28, V31, V33, V35
V33 -- V5, V13, V18, V31, V32
V34 -- V10, V12, V16, V17, V19, V29
V35 -- V15, V18, V19, V25, V28, V31, V32
V36 -- V15, V16, V24
由于潛變量網(wǎng)絡(luò)模型的特征模型形式能夠同時(shí)反映距離和同質(zhì)性因素，對(duì)于我們的數(shù)據(jù)，我們可以分別擬合并比較三種不同的特征模型：
- 沒有成對(duì)協(xié)變量
- 共同執(zhí)業(yè)協(xié)變量
- 共同辦公地點(diǎn)協(xié)變量
不包含成對(duì)的協(xié)變量，潛在特征空間維度為2
library(eigenmodel)
set.seed(42)
A <- as_adjacency_matrix(lazega, sparse=FALSE)
lazega.leig.fit1 <- eigenmodel_mcmc(A, R=2, S=11000,
burn=10000) # Approximate the posterior distribution of parameters in an eigenmodel
5 percent done (iteration 1050) Thu Oct 01 08:59:26 2020
10 percent done (iteration 2100) Thu Oct 01 08:59:31 2020
15 percent done (iteration 3150) Thu Oct 01 08:59:37 2020
20 percent done (iteration 4200) Thu Oct 01 08:59:41 2020
25 percent done (iteration 5250) Thu Oct 01 08:59:45 2020
30 percent done (iteration 6300) Thu Oct 01 08:59:49 2020
35 percent done (iteration 7350) Thu Oct 01 08:59:53 2020
40 percent done (iteration 8400) Thu Oct 01 08:59:57 2020
45 percent done (iteration 9450) Thu Oct 01 09:00:00 2020
50 percent done (iteration 10500) Thu Oct 01 09:00:05 2020
55 percent done (iteration 11550) Thu Oct 01 09:00:09 2020
60 percent done (iteration 12600) Thu Oct 01 09:00:14 2020
65 percent done (iteration 13650) Thu Oct 01 09:00:18 2020
70 percent done (iteration 14700) Thu Oct 01 09:00:23 2020
75 percent done (iteration 15750) Thu Oct 01 09:00:27 2020
80 percent done (iteration 16800) Thu Oct 01 09:00:31 2020
85 percent done (iteration 17850) Thu Oct 01 09:00:35 2020
90 percent done (iteration 18900) Thu Oct 01 09:00:39 2020
95 percent done (iteration 19950) Thu Oct 01 09:00:43 2020
100 percent done (iteration 21000) Thu Oct 01 09:00:47 2020
定義分組
same.prac.op <- v.attr.lazega$Practice %o%
v.attr.lazega$Practice
same.prac <- matrix(as.numeric(same.prac.op
%in% c(1, 4, 9)), 36, 36)
same.prac <- array(same.prac,dim=c(36, 36, 1))
使用分組信息擬合
lazega.leig.fit2 <- eigenmodel_mcmc(A, same.prac, R=2,
+ S=11000,burn=10000)
5 percent done (iteration 1050) Thu Oct 01 09:08:47 2020
10 percent done (iteration 2100) Thu Oct 01 09:08:53 2020
15 percent done (iteration 3150) Thu Oct 01 09:08:58 2020
20 percent done (iteration 4200) Thu Oct 01 09:09:05 2020
25 percent done (iteration 5250) Thu Oct 01 09:09:10 2020
30 percent done (iteration 6300) Thu Oct 01 09:09:15 2020
35 percent done (iteration 7350) Thu Oct 01 09:09:20 2020
40 percent done (iteration 8400) Thu Oct 01 09:09:26 2020
45 percent done (iteration 9450) Thu Oct 01 09:09:31 2020
50 percent done (iteration 10500) Thu Oct 01 09:09:36 2020
55 percent done (iteration 11550) Thu Oct 01 09:09:41 2020
60 percent done (iteration 12600) Thu Oct 01 09:09:47 2020
65 percent done (iteration 13650) Thu Oct 01 09:09:52 2020
70 percent done (iteration 14700) Thu Oct 01 09:09:58 2020
75 percent done (iteration 15750) Thu Oct 01 09:10:03 2020
80 percent done (iteration 16800) Thu Oct 01 09:10:09 2020
85 percent done (iteration 17850) Thu Oct 01 09:10:16 2020
90 percent done (iteration 18900) Thu Oct 01 09:10:22 2020
95 percent done (iteration 19950) Thu Oct 01 09:10:28 2020
100 percent done (iteration 21000) Thu Oct 01 09:10:35 2020
最后，我們將共同辦公地點(diǎn)加到模型里面
same.off.op <- v.attr.lazega$Office %o%
v.attr.lazega$Office
same.off <- matrix(as.numeric(same.off.op %in%
c(1, 4, 9)), 36, 36)
same.off <- array(same.off,dim=c(36, 36, 1))
lazega.leig.fit3 <- eigenmodel_mcmc(A, same.off,
R=2, S=11000, burn=10000)
為了比較網(wǎng)絡(luò)我們?cè)诿總€(gè)模型推斷得到的二維潛在變量空間表示有何不同，我們提取了每個(gè)擬合模型的特征向量。
lat.sp.1 <-
eigen(lazega.leig.fit1$ULU_postmean)$vec[, 1:2]
lat.sp.2 <-
eigen(lazega.leig.fit2$ULU_postmean)$vec[, 1:2]
lat.sp.3 <-
eigen(lazega.leig.fit3$ULU_postmean)$vec[, 1:2]
并用這些坐標(biāo)作為布局在igraph中繪制這一網(wǎng)絡(luò)。
colbar <- c("red", "dodgerblue", "goldenrod")
v.colors <- colbar[V(lazega)$Office]
v.shapes <- c("circle", "square")[V(lazega)$Practice]
v.size <- 3.5*sqrt(V(lazega)$Years)
v.label <- V(lazega)$Seniority
plot(lazega, layout=lat.sp.1, vertex.color=v.colors,
vertex.shape=v.shapes, vertex.size=v.size,
vertex.label=(1:36))
擬合優(yōu)度
使用交叉驗(yàn)證
perm.index <- sample(1:630)
nfolds <- 5
nmiss <- 630/nfolds
Avec <- A[lower.tri(A)]
Avec.pred1 <- numeric(length(Avec))
# 交叉驗(yàn)證過程
for(i in seq(1,nfolds)){
# Index of missing values.
miss.index <- seq(((i-1) * nmiss + 1),
(i*nmiss), 1)
A.miss.index <- perm.index[miss.index]
# Fill a new Atemp appropriately with NA's.
Avec.temp <- Avec
Avec.temp[A.miss.index] <-
rep("NA", length(A.miss.index))
Avec.temp <- as.numeric(Avec.temp)
Atemp <- matrix(0, 36, 36)
Atemp[lower.tri(Atemp)] <- Avec.temp
Atemp <- Atemp + t(Atemp)
# Now fit model and predict.
Y <- Atemp
model1.fit <- eigenmodel_mcmc(Y, R=2,
S=11000, burn=10000)
model1.pred <- model1.fit$Y_postmean
model1.pred.vec <-
model1.pred[lower.tri(model1.pred)]
Avec.pred1[A.miss.index] <-
model1.pred.vec[A.miss.index]
}
使用ROC曲線評(píng)估模型
ROC曲線：接收者操作特征曲線（receiver operating characteristic curve），是反映敏感性和特異性連續(xù)變量的綜合指標(biāo)，roc曲線上每個(gè)點(diǎn)反映著對(duì)同一信號(hào)刺激的感受性。
對(duì)于分類器，或者說分類算法，評(píng)價(jià)指標(biāo)主要有precision，recall，F(xiàn)-score等，以及這里要討論的ROC和AUC
library(ROCR)
pred1 <- prediction(Avec.pred1, Avec)
perf1 <- performance(pred1, "tpr", "fpr")
plot(perf1, col="blue", lwd=3)
# CHUNK 34
perf1.auc <- performance(pred1, "auc")
slot(perf1.auc, "y.values")
[[1]]
[1] 0.819181
http://www.bu.edu/cs/files/2014/05/Kolaczyk.pdf
https://www.r-bloggers.com/2016/05/ergm-tutorial/
https://baike.baidu.com/item/%E7%A4%BE%E4%BC%9A%E7%BD%91%E7%BB%9C%E6%8C%87%E6%95%B0%E9%9A%8F%E6%9C%BA%E5%9B%BE%E6%A8%A1%E5%9E%8B%EF%BC%9A%E7%90%86%E8%AE%BA%E3%80%81%E6%96%B9%E6%B3%95%E4%B8%8E%E5%BA%94%E7%94%A8/20351978
https://blog.csdn.net/weixin_40895857/article/details/105184024
https://scholar.princeton.edu/sites/default/files/bstewart/files/soc-stats-ergms.pdf
https://www.cs.du.edu/~meiyin/Talk.pdf
http://www.princeton.edu/~eabbe/publications/abbe_FNT_2.pdf
https://hal.archives-ouvertes.fr/hal-00730829/file/model_selection_in_block_clustering_by_the_icl.pdf
https://cloud.tencent.com/developer/article/1595157
https://blog.csdn.net/g8015108/article/details/69367602
https://psych-networks.com/meaning-model-equivalence-network-models-latent-variables-theoretical-space/#:~:text=The%20paper%20constructs%20elegant%20representations%20of%20the%20Ising,the%20latent%20variable%20acts%20as%20a%20common%20effect%29.
http://www.stat.cmu.edu/~brian/905-2009/all-papers/hoff-raftery-handcock-2002-jasa.pdf
https://www.ijcai.org/Proceedings/11/Papers/286.pdf
https://www.r-bloggers.com/2011/01/example-8-21-latent-class-analysis/
