Statistical Models for Network Graphs
Introduction
Robins andMorris write
“A good model(Statistical network graph) needs to be both estimable from data and a reasonable representation of that data, to be theoretically plausible about the type of effects that might have produced the network, and to be amenable to examining which competing effects might be the best explanation of the data.”
the three main such classes of models developed to date closely parallel more familiar statistical models for non-network datasets.
Type of Network Models.
we can think of two classes of network models:
1. Mathematical/ Probabilistic Models
-> chapter5. Mathematical Models for Network Graphs
2. Statistical Models
-> chapter6. Statistical Models for Network Graphs
Three Canonical classes of Statistical Models
there are network-based version of three canonical classes of statistical models :
- regression model(i.e ERGMs)
- mixture models (i.e Network Blocks Models)
- mixed effects models (i.e Latent Variable Models)
table of contents
- 6.2 Exponential Random Graph Models
General Formulation
Specifying a Model
Model Fitting
Goodness-of-Fit
Model Specification
Model Fitting
Goodness-of-Fit
- 6.4 Latent Network Models
General Formulation
Specifying the Latent Effects
Model Fitting
Goodness-of-Fit
6.2 Exponential Random Graph Models
\[P_{\theta}(Y=y) = (1/k)*\exp{\sum_{H}{\theta}_{H}g_{H}(y)}\]
We will illustrate the construction, fitting, and assessment of ERGMs using the lazega data set on collaboration among lawyers, introduced in Chap. 1.
Within R, easily using the package for ERGMs is the ergm.
first separating the network into adjacency matrix and attributes
lazega is an igraph graph object, undirected. It has the following vertex attributes:
‘name’, ‘Seniority’, ‘Status’ (all 1, meaning partner),
‘Gender’ (1 is man, 2 is woman), ‘Office’ (1 is Boston, 2 is Hartford, 3 is Providence),
‘Years’ (years with the firm), ‘Age’, ‘Practice’ (1 is litigation, 2 is corporate),
and ‘School’ (1 is Harvard or Yale, 2 is University of Connecticut, 3 is other).
See the reference below for more
library(sand)
## Loading required package: igraph
## Loading required package: igraphdata
##
## Statistical Analysis of Network Data with R
## Type in C2 (+ENTER) to start with Chapter 2.
data(lazega)
summary(lazega)
## IGRAPH 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)
A <- get.adjacency(lazega) #Convert a graph to an adjacency matrix(인접행렬만들기)
A
## 36 x 36 sparse Matrix of class "dgCMatrix"
## [[ suppressing 36 column names 'V1', 'V2', 'V3' ... ]]
##
## V1 . . . . . . . . . . . . . . . . 1 . . . . . . . . . . . . . . . . . .
## V2 . . . . . . 1 . . . . . . . . 1 1 . . . . 1 . . . 1 . . 1 . . . . . .
## V3 . . . . . . . . . . . . . . . . . 1 . . . . . . 1 . . 1 . . . . . . .
## V4 . . . . . . . . . . . 1 . . . . 1 . 1 1 . 1 . . . 1 . 1 1 . 1 . . . .
## V5 . . . . . . . . . . . . . . . . . 1 . . . . . 1 . . . 1 . . 1 1 1 . .
## V6 . . . . . . . . . . . . . . . . . . . . . . . 1 . . . 1 . 1 1 1 . . .
## V7 . 1 . . . . . . . . . . . . . . . 1 . . . . . . . . . . . . . . . . .
## V8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## V9 . . . . . . . . . . . 1 . . . 1 . . . . . . . . . . . . 1 . . . . . .
## V10 . . . . . . . . . . . . . . . . . . . . . . . 1 . 1 . . 1 . 1 . . 1 .
## V11 . . . . . . . . . . . . . . . . 1 . . . . . . . . . . . . . . . . . .
## V12 . . . 1 . . . . 1 . . . . . 1 1 1 . 1 . . . . . . 1 . . 1 . . . . 1 .
## V13 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . 1 . .
## V14 . . . . . . . . . . . . . . . 1 1 . . . . . . . 1 . . 1 . 1 . 1 . . .
## V15 . . . . . . . . . . . 1 . . . 1 . . 1 1 . 1 . 1 . 1 . . 1 . . 1 . . 1
## V16 . 1 . . . . . . 1 . . 1 . 1 1 . 1 . . . . 1 . . . 1 1 . 1 . . 1 . 1 .
## V17 1 1 . 1 . . . . . . 1 1 . 1 . 1 . . 1 . . 1 . 1 1 1 . 1 1 . . . . 1 .
## V18 . . 1 . 1 . 1 . . . . . . . . . . . . . . . . . . . . 1 . . 1 1 1 . 1
## V19 . . . 1 . . . . . . . 1 . . 1 . 1 . . . . 1 . 1 . 1 . 1 . . . . . 1 1
## V20 . . . 1 . . . . . . . . . . 1 . . . . . . 1 . . . 1 . . . . . . . . .
## V21 . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . .
## V22 . 1 . 1 . . . . . . . . . . 1 1 1 . 1 1 . . . . . . . . . . 1 1 . . .
## V23 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## V24 . . . . 1 1 . . . 1 . . . . 1 . 1 . 1 . . . . . . 1 . . . . 1 . . . .
## V25 . . 1 . . . . . . . . . . 1 . . 1 . . . . . . . . . . 1 . . . . . . 1
## V26 . 1 . 1 . . . . . 1 . 1 . . 1 1 1 . 1 1 . . . 1 . . 1 . . . . 1 . . .
## V27 . . . . . . . . . . . . . . . 1 . . . . 1 . . . . 1 . . . . . . . . .
## V28 . . 1 1 1 1 . . . . . . . 1 . . 1 1 1 . . . . . 1 . . . . 1 1 1 . . 1
## V29 . 1 . 1 . . . . 1 1 . 1 . . 1 1 1 . . . . . . . . . . . . . . . . 1 .
## V30 . . . . . 1 . . . . . . . 1 . . . . . . . . . . . . . 1 . . 1 . . . .
## V31 . . . 1 1 1 . . . 1 . . 1 . . . . 1 . . . 1 . 1 . . . 1 . 1 . 1 1 . 1
## V32 . . . . 1 1 . . . . . . . 1 1 1 . 1 . . . 1 . . . 1 . 1 . . 1 . 1 . 1
## V33 . . . . 1 . . . . . . . 1 . . . . 1 . . . . . . . . . . . . 1 1 . . .
## V34 . . . . . . . . . 1 . 1 . . . 1 1 . 1 . . . . . . . . . 1 . . . . . .
## V35 . . . . . . . . . . . . . . 1 . . 1 1 . . . . . 1 . . 1 . . 1 1 . . .
## V36 . . . . . . . . . . . . . . 1 1 . . . . . . . 1 . . . . . . . . . . .
##
## V1 .
## V2 .
## V3 .
## V4 .
## V5 .
## V6 .
## V7 .
## V8 .
## V9 .
## V10 .
## V11 .
## V12 .
## V13 .
## V14 .
## V15 1
## V16 1
## V17 .
## V18 .
## V19 .
## V20 .
## V21 .
## V22 .
## V23 .
## V24 1
## V25 .
## V26 .
## V27 .
## V28 .
## V29 .
## V30 .
## V31 .
## V32 .
## V33 .
## V34 .
## V35 .
## V36 .
what : Character constant, whether to return info about vertices, edges, or both.The default is ‘edges’.
get.adjacency() 함수를 이용하여 인접행렬을 만들면, 그래프에는 들어있는 노드에 대한 속성값을 상실한다. 따라서 igraph에서 ergm으로 자료를 이동하면서 속성값도 같이 이동하기 위해 get.data.frame(, what=“vertices”)를 이용하여 노드의 속성을 추출한다. 여기서 what 매개변수를 지정하지 않을 경우 edge list가 출력된다
v.attrs <- get.data.frame(lazega, what="vertices")#Creating igraph graphs from data frames
head(v.attrs)
## name Seniority Status Gender Office Years Age Practice School
## V1 V1 1 1 1 1 31 64 1 1
## V2 V2 2 1 1 1 32 62 2 1
## V3 V3 3 1 1 2 13 67 1 1
## V4 V4 4 1 1 1 31 59 2 3
## V5 V5 5 1 1 2 31 59 1 2
## V6 V6 6 1 1 2 29 55 1 1
str(v.attrs)
## 'data.frame': 36 obs. of 9 variables:
## $ name : chr "V1" "V2" "V3" "V4" ...
## $ Seniority: int 1 2 3 4 5 6 7 8 9 10 ...
## $ Status : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Gender : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Office : int 1 1 2 1 2 2 2 1 1 1 ...
## $ Years : int 31 32 13 31 31 29 29 28 25 25 ...
## $ Age : int 64 62 67 59 59 55 63 53 53 53 ...
## $ Practice : int 1 2 1 2 1 1 2 1 2 2 ...
## $ School : int 1 1 1 3 2 1 3 3 1 3 ...
creating the analogous network object for ergm
## Loading required package: statnet.common
## Loading required package: network
## network: Classes for Relational Data
## Version 1.11.3 created on 2014-12-05.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
## Mark S. Handcock, University of California -- Los Angeles
## David R. Hunter, Penn State University
## Martina Morris, University of Washington
## Skye Bender-deMoll, University of Washington
## For citation information, type citation("network").
## Type help("network-package") to get started.
##
##
## Attaching package: 'network'
##
## The following objects are masked from 'package:igraph':
##
## %c%, %s%, add.edges, add.vertices, delete.edges,
## delete.vertices, get.edge.attribute, get.edges,
## get.vertex.attribute, is.bipartite, is.directed,
## list.edge.attributes, list.vertex.attributes,
## set.edge.attribute, set.vertex.attribute
##
##
## ergm: version 3.2.4, created on 2014-12-13
## Copyright (c) 2014, Mark S. Handcock, University of California -- Los Angeles
## David R. Hunter, Penn State University
## Carter T. Butts, University of California -- Irvine
## Steven M. Goodreau, University of Washington
## Pavel N. Krivitsky, University of Wollongong
## Martina Morris, University of Washington
## with contributions from
## Li Wang
## Kirk Li, University of Washington
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("ergm").
##
## NOTE: If you use custom ERGM terms based on 'ergm.userterms'
## version prior to 3.1, you will need to perform a one-time update
## of the package boilerplate files (the files that you did not write
## or modify) from 'ergm.userterms' 3.1 or later. See
## help('eut-upgrade') for instructions.
앞서 추출한 인접행렬을 이용하여 network 패키지의 네트워크 객체를 만들고, 해당 네트워크 객체에 get.data.frame()을 이용하여 추출한 노드 속성을 set.vertex.attribute()를 이용하여 노드의 속성을 지정하였다.
lazega.s <- network::as.network(as.matrix(A),directed=FALSE)
lazega.s
## Network attributes:
## vertices = 36
## directed = FALSE
## hyper = FALSE
## loops = FALSE
## multiple = FALSE
## bipartite = FALSE
## total edges= 115
## missing edges= 0
## non-missing edges= 115
##
## Vertex attribute names:
## vertex.names
##
## No edge attributes
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::’는 set.vertex.attribute() 함수가 network 패키지에 있다는 의미이다. set.vertex.attribute()함수가 igraph와 network에 존재하기 때문에 패키지를 명시하지 않을 경우 잘못된 결과를 얻을수 있다.
6.2.2 Specifying a Model
We have already seen an example of what is arguably the simplest ERGM, in the form of the Bernoulli random graph model chap.5.2.
(우리는 이미 5.2절에서 ERGM의 가장 간단한 베르누이 랜덤그래프모델의 형태를 보았다)
To see this suppose, we specify that, for a given pair of vertices, the presence or absence of an edge between that pair is independent of the status of possible edges between any other pairs of vertices.
(모든 주어진 점들의 쌍에 대해 모든 점의 쌍들 사이의 가능한 선들의 존재유무에 대해 독립인 상태라고 가정한다. 즉 주어진 하나의 선의 상태는 다른 선과 는 관계없다. )
-> then \({\theta}_H = 0\) for all configurations H involving three or more vertices.
(그래서 모든 점들의 H배열에 대해 영이다.\({\theta}_H =0\) )
As, a result, the ERGM in (6.2) reduces to
(이결과를 통해 (6.2)의 식ERGM을 줄여볼 수 있다.)
\(P_{\theta}(Y=y) = (1/k)\exp{\sum_{i,j}{\theta}_{i,j}y_{i,j}}\) (6.4)
Furthermore, if we assume that the coefficients \({\theta_{ij}}\) are equal to some common value \({\theta}\) are equal to some common value \({\theta}\), then further simplifies to
\(P_{\theta}(Y=y) = (1/k)\exp{{\theta}L(y)}\) (6.5)
where \(L(y) = \sum_{i,j}y_{ij}= N_{e}\) is the number of edges in the graph.
(더욱이, 만일 계수 \({\theta_{i,j}}\)가 몇몇 공통의 값인 \({/theta}\)와 동등하다면, 6.4를 6.5로 단순화시킬수 있다. 이때 L(y)는 그래프에서 선의 수이다. )
To, specify models in ergm, we use the function formula and standard R syntax.
For example, model (6.5) may be specified for the network lazega.s as
m <- matrix(0,4,4)
m[1,2]<-1;m[1,3]<-1;m[2,3]<-1;m[3,4]<-1
m[2,1]<-1;m[3,1]<-1;m[3,2]<-1;m[4,2]<-1
g <- network(m)
g.m <- formula(g~edges)
summary(g.m)
## edges
## 8
str(g.m)
## Class 'formula' length 3 g ~ edges
## ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
summary.statistics(g.m)
## edges
## 8
my.ergm.bern <- formula(lazega.s ~edges)
my.ergm.bern
## lazega.s ~ edges
in which case the statistics L takes the value
summary.statistics(my.ergm.bern)
## edges
## 115
formula() 그래프이 연결선을 반응변수로 하는 모형을 작성한다. 위에서 언급한 edges는 연결선 수를 연결선을 설명하는 변수로 하겠다는 의미이다.
6.3 Network Block Model
Suppose that each vertex i∈ V of a graph G=(V,E) can belong to one of Q classes, say \(g_{1},....g_{Q}\)
(class label q and r of vertices i and j)
A block model for G specifies that each element \(Y_{ij}\) of the adjacency matrix Y is, an independent Bernoulli random variable with probability \({\pi}_{qr}\).
Undirected graph, \({\pi}_{qr} = {\pi}_{rq}\)
The block model is hence a variant of the Bernoulli random graph model, where the probabilities of an edge are restricted to be one of Q^2 possible values \({\pi}_{qr}\).
\[P_{\theta}(Y=y) = (1/k)\exp{\sum_{q,r}{\theta}L_{qr}(y)}\] \({L_{qr}(y)}\) is the number of edges in the observed graph y connecting pairs of vertices of classes q and r.
Nevertheless, the assumption that the class memebership of vertices is known , but in general, we are considered untenable in practice. So we use the SBM.
6.3 Stochastic Block Models(SBM)
This model specifies only that there are Q classes, for some Q, but does not specify the nature of those classes nor the class membership of the individual vertices.
(모델의 Q개의 클래스를 지정할 뿐, 클래의 성질도 각점의 클래스 구성원을 지정하지 않습니다.)
Stochastic block models explicitly parameterize the notion of
- groups/modules, with
- different rates of connetions between/within.
They are a hierarchical version of an ERGM, in which one level assigns group memebership and the second level is block model.
Often used in community detection, where the focus is solely on the inference of group membership, resulting in a graph partitioning algorithm.
More generally, when modeling we may be interested in the parameters, the groups, or both!
SBM: Model Specification
This model specifies only that there are Q classes, for some Q, but does not specify the nature of those classes nor the class membership of the individual vertices.
It dicates simply that the class memebership of each vertex i be determined independently, according to a common distribution on the set{1,…,Q}.
Formally, let \(Z_{iq}\) = 1 if vertex i is of class q, and zero otherwise. Under a stochastic block model, the vectors \(Z_i = (Z_{i1}, . . . ,Z_{iQ})\) are determined independently, where \[P(Z_{iq} = 1) = α_q\] and \[{\sum_{q=1}^{Q}}{α_q}= 1.\]
A stochastic block model is thus, effectively, a mixture of classical random graph models.
this class of models as a ‘mixture model for random graphs’.
SBM: Model Fitting
In the case of stochastic block models, both the (now conditional) edge probabilities \(π_{qr}\) and the class membership probabilities \(α_q\) must estimate.
참조 6.18
the complete-data log-likelihood, is of the form where \(b(y;π) =π^y{1-π}^{1-y}\).
While this may not seem like much of a change over the ordinary block model, the task of model fitting becomes decidedly more complex in this setting.
Traditinally, inference uses the expectation-maximization (EM)algorithm ; more recent work includes the use of variational methods.
To illustrate, we use the network fblog of French political blogs introduced in Chap.3.5.
Note that, we have specifed only that the total number of classes Q be between 2 and 15.
mixer to select the number of classes fit to the network.
library(sand)
library(mixer)
setSeed(42)
data(fblog)
summary(fblog)
## IGRAPH UN-- 192 1431 --
## attr: name (v/c), PolParty (v/c)
fblog.sbm <- mixer(as.matrix(get.adjacency(fblog)),qmin=2, qmax=15)
## Mixer: the adjacency matrix has been transformed in a undirected edge list
#class 2~15를 관찰네트워크그래프와 적합시킨다.
str(fblog.sbm)
## List of 9
## $ method : chr "variational"
## $ nnames :List of 1
## ..$ : chr [1:192] " jeunesverts.org/bordeaux" " bix.enix.org/" " www.arnaudcaron.net/" " dominiquevoynet.net/blog" ...
## $ nnodes : int 192
## $ map : int [1:192] 1 2 3 4 5 6 7 11 12 13 ...
## $ edges : int [1:2, 1:1431] 1 2 2 3 1 4 2 4 3 4 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "row" "col"
## .. ..$ : chr [1:1431] " jeunesverts.org/bordeaux" " bix.enix.org/" " jeunesverts.org/bordeaux" " bix.enix.org/" ...
## $ qmin : num 2
## $ qmax : num 15
## $ output :List of 14
## ..$ :List of 4
## .. ..$ criterion: num -4632
## .. ..$ alphas : num [1:2] 0.778 0.222
## .. ..$ Pis : num [1:2, 1:2] 0.0758 0.0332 0.0332 0.4279
## .. ..$ Taus : num [1:2, 1:192] 1.00 1.22e-10 9.97e-01 2.93e-03 1.00 ...
## ..$ :List of 4
## .. ..$ criterion: num -4408
## .. ..$ alphas : num [1:3] 0.644 0.172 0.184
## .. ..$ Pis : num [1:3, 1:3] 0.0803 0.0198 0.0479 0.0198 0.374 ...
## .. ..$ Taus : num [1:3, 1:192] 1.00 5.17e-05 1.51e-10 9.92e-01 1.00e-10 ...
## ..$ :List of 4
## .. ..$ criterion: num -4158
## .. ..$ alphas : num [1:4] 0.536 0.171 0.163 0.13
## .. ..$ Pis : num [1:4, 1:4] 0.082 0.0238 0.0678 0.0161 0.0238 ...
## .. ..$ Taus : num [1:4, 1:192] 1.00 4.85e-05 1.00e-10 1.12e-06 9.99e-01 ...
## ..$ :List of 4
## .. ..$ criterion: num -4035
## .. ..$ alphas : num [1:5] 0.492 0.0573 0.1711 0.1542 0.1253
## .. ..$ Pis : num [1:5, 1:5] 0.0887 0.0164 0.0255 0.0776 0.0163 ...
## .. ..$ Taus : num [1:5, 1:192] 1.00 1.00e-10 7.36e-05 1.00e-10 1.84e-06 ...
## ..$ :List of 4
## .. ..$ criterion: num -3932
## .. ..$ alphas : num [1:6] 0.3735 0.1387 0.0573 0.1718 0.1334 ...
## .. ..$ Pis : num [1:6, 1:6] 0.10787 0.02466 0.00931 0.02939 0.11626 ...
## .. ..$ Taus : num [1:6, 1:192] 9.75e-01 2.49e-02 1.00e-10 2.27e-04 1.00e-10 ...
## ..$ :List of 4
## .. ..$ criterion: num -3786
## .. ..$ alphas : num [1:7] 0.3538 0.1411 0.0573 0.1457 0.0417 ...
## .. ..$ Pis : num [1:7, 1:7] 0.09008 0.01659 0.00804 0.00526 0.16431 ...
## .. ..$ Taus : num [1:7, 1:192] 9.97e-01 2.76e-03 1.00e-10 1.73e-04 1.00e-10 ...
## ..$ :List of 4
## .. ..$ criterion: num -3735
## .. ..$ alphas : num [1:8] 0.2489 0.1228 0.0573 0.1421 0.0417 ...
## .. ..$ Pis : num [1:8, 1:8] 0.08182 0.02879 0.015 0.00771 0.19897 ...
## .. ..$ Taus : num [1:8, 1:192] 1.00 4.23e-04 1.00e-10 3.84e-05 1.00e-10 ...
## ..$ :List of 4
## .. ..$ criterion: num -3712
## .. ..$ alphas : num [1:9] 0.2496 0.1229 0.0573 0.136 0.0312 ...
## .. ..$ Pis : num [1:9, 1:9] 0.08151 0.02867 0.01491 0.00721 0.09361 ...
## .. ..$ Taus : num [1:9, 1:192] 9.99e-01 4.60e-04 1.00e-10 3.05e-05 1.00e-10 ...
## ..$ :List of 4
## .. ..$ criterion: num -3700
## .. ..$ alphas : num [1:10] 0.2367 0.1217 0.0573 0.1358 0.0312 ...
## .. ..$ Pis : num [1:10, 1:10] 0.08322 0.02885 0.01196 0.00736 0.09861 ...
## .. ..$ Taus : num [1:10, 1:192] 1.00 3.63e-04 1.00e-10 3.23e-05 1.00e-10 ...
## ..$ :List of 4
## .. ..$ criterion: num -3706
## .. ..$ alphas : num [1:11] 0.2367 0.1219 0.0573 0.1358 0.0312 ...
## .. ..$ Pis : num [1:11, 1:11] 0.08319 0.0289 0.01188 0.00758 0.09898 ...
## .. ..$ Taus : num [1:11, 1:192] 1.00 3.48e-04 1.00e-10 2.33e-05 1.00e-10 ...
## ..$ :List of 4
## .. ..$ criterion: num -3693
## .. ..$ alphas : num [1:12] 0.1511 0.1273 0.1081 0.0573 0.1358 ...
## .. ..$ Pis : num [1:12, 1:12] 0.090964 0.005527 0.049616 0.000001 0.002688 ...
## .. ..$ Taus : num [1:12, 1:192] 1.00 1.40e-05 9.64e-06 1.00e-10 1.34e-06 ...
## ..$ :List of 4
## .. ..$ criterion: num -3752
## .. ..$ alphas : num [1:13] 0.1439 0.1302 0.1173 0.0573 0.1357 ...
## .. ..$ Pis : num [1:13, 1:13] 0.098337 0.005721 0.048055 0.000001 0.002781 ...
## .. ..$ Taus : num [1:13, 1:192] 1.00 1.27e-05 4.63e-06 1.00e-10 1.26e-06 ...
## ..$ :List of 4
## .. ..$ criterion: num -3777
## .. ..$ alphas : num [1:14] 0.1555 0.1257 0.0803 0.0573 0.136 ...
## .. ..$ Pis : num [1:14, 1:14] 0.093541 0.006388 0.033716 0.000001 0.002558 ...
## .. ..$ Taus : num [1:14, 1:192] 1.00 2.47e-06 2.83e-06 1.00e-10 1.15e-10 ...
## ..$ :List of 4
## .. ..$ criterion: num -3854
## .. ..$ alphas : num [1:15] 0.1513 0.1285 0.1151 0.0573 0.1356 ...
## .. ..$ Pis : num [1:15, 1:15] 0.090478 0.005411 0.039594 0.000001 0.002601 ...
## .. ..$ Taus : num [1:15, 1:192] 1.00 1.50e-05 7.66e-06 1.00e-10 1.38e-06 ...
## $ directed: logi FALSE
## - attr(*, "class")= chr "mixer"
Examining the model output
fblog.sbm.output <- getModel(fblog.sbm) #get the parameters of a model
names(fblog.sbm.output)
## [1] "q" "criterion" "alphas" "Pis" "Taus"
we see that the network of French blogs has been fit with classes, in estimated propotions(측정된 비율)
fblog.sbm.output$q
## [1] 12
fblog.sbm.output$alphas
## [1] 0.15106842 0.12727445 0.10813837 0.05729167 0.13584314 0.03123703
## [7] 0.12386173 0.09333667 0.01041667 0.02089779 0.12500907 0.01562500
The output from a fitted stochastic block model also allows for the assignment of vertices to classes.
graph partitioning
examining the estimates for the first three vertices in the French blog network (1:3까지의 점이 어떤 class에 속하는지 파악가능)
fblog.sbm.output$Taus[,1:3]
## [,1] [,2] [,3]
## [1,] 9.999747e-01 0.0049141866 9.999814e-01
## [2,] 1.395463e-05 0.0000000001 7.569293e-07
## [3,] 9.639159e-06 0.9942175617 1.772634e-05
## [4,] 1.000000e-10 0.0000000001 1.000000e-10
## [5,] 1.338119e-06 0.0000000001 8.194924e-09
## [6,] 1.000000e-10 0.0000000001 1.000000e-10
## [7,] 3.550935e-07 0.0008682507 6.117352e-08
## [8,] 1.000000e-10 0.0000000001 1.000000e-10
## [9,] 1.000000e-10 0.0000000001 1.000000e-10
## [10,] 1.000000e-10 0.0000000001 1.000000e-10
## [11,] 4.432111e-09 0.0000000001 1.000000e-10
## [12,] 1.000000e-10 0.0000000001 1.000000e-10
we see that an assignment rule based on the maximum a posteriori criterion(사후판정법) would place the first and third vertices in class1, and the second, in class3.
6.3.3. Goodness-of-Fit
In assessing the goodness-of-fit of a stochastic block model we could, of course, use the same types of simulation-based methods we employed in the anaysis of ERGMs.
A selection of summaries produced by the mixer package are displayed in Fig.6.2.
plot(fblog.sbm, classes=as.factor(V(fblog)$PolParty))
ICL : seems to suggest there is some latitude in this choice, with anywhere from 8 to 12 classes being reasonable.
adjacency matrix reorganized according to the estimated partition for a given number of classes 12
(12개의 클래스에 따른 측정값을 이용한 인접행렬)
degree distribution (histogram) and theoretical degree distribution (blue curve) computed from the q-class model parameters (alphas, Pis).
Providing external classes each node displays a pie chart pointing out the classification relevance(분류관련성).
6.4.Latent Network Models
From the perspective of statistical modeling, one key innovation underlying stochastic block models and their extensions is the incorporation of latent variables, in the form of vertex classes.
That is, the use of variabels that a unobserved but which play a role in determining the probability that vertex airs are incident to each other.
Latent Eigenspace Model
Consider the modeling of Y=\([Y_{ij}]\), indicating presence/absence of edges in a graph G, as a classification problem.
A standard logistic regression classifier in this context might be based on models of the form
\[\log[\frac{P_{\beta} (Y_{ij} = 1 | Z_{ij}=z)}{P_{\beta} (Y_{ij} = 0 | Z_{ij}=z)}] = {\beta}^T z\]
where
- \(Z_{ij}\) is a vector of explanatory variabels indexed in the unordered pairs {i,j}, and
- \({\beta}\) is a vector of regression coefficients, assumed common to all pairs.
Problem : This type of model assumes \(Y_{ij}'s\) conditionally independent given \(Z_{ij}'s\). Unlikely to be ture here!
Solulation : Introduce latent vairabels M=[\(M_{ij}\)], and assume conditional independence of \(Y_{ij}'s\) given both \(Z_{ij}'s\) and \(M_{ij}\).
Let \[M=U^{T}ㅅU +E\]be an(unknown) random, symmetric \(N_v * N_v\) matrix. Then replace the standard classification model with
\[\log[\frac{P_{\beta} (Y_{ij} = 1 | Z_{ij}=z)}{P_{\beta} (Y_{ij} = 0 | Z_{ij}=z)}] = {\beta}^T z + m\]
Fitting Latent Eigenspace Models
The latent variable matrix M is intended to capture effects of network structural characteristics or process not already described by the observed explanatory variables \(Z_{ij}\)
Additional distributional assumptions needed for \({\beta}\) and components of M
Hoff Proposes
- U uniform on the space of all \({N_v * N_v}\) orthonormal matrices, and
- diag(ㅅ), E, and \({\beta}\) multivariate Gaussian.
MCMC used to simulate for posterior-based inference.
practice
The packages eigenmodel implements the eigenmodel formulation. The function eigenmodel_mcmc in eigenmodel uses Monte Carlo Markov Chain(MCMC) techniques to simulate from the relevant posterior distribution, using largely conjugate priors to complete the model specification.
The network lazega of collaborations among lawyers allows for demonstration of a number of the concepts we have discussed so far.
Recall that this network involved 36 lawyers, at three different office location, involved in two types of practice(i.e , corporate and litigation).
summary(lazega)
## IGRAPH 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)
Latent network models is able to capture aspects of both distance and homophily, it is interesting to compare the fitted models that we obtain for three different eigenmodels,
(i) no pair-specific covariates,
(ii)a covariate for common practice ,
(iii)a covariate for shared office location.
To fit the model with no pair-specific covariates and a latent eigen-space of Q=2 dismensions is accomplished as follows.
(잠재적 네트워크 모델은 거리와 homophily(유사한거끼리 연결하는 성질) 두 측면을 포착 할 수있다. 이 특징은 3개의 다른 고유모델에 대하여 적합모델 비교에 이용할수 있다.)
(i) no pair-specific covariates
library(eigenmodel)
set.seed(42)
A <- get.adjacency(lazega, sparse=FALSE)
lazega.leig.fit1 <- eigenmodel_mcmc(A,R=2,S=11000,burn=10000)
## 5 percent done (iteration 550) Mon Apr 13 17:47:22 2015
## 10 percent done (iteration 1100) Mon Apr 13 17:47:26 2015
## 15 percent done (iteration 1650) Mon Apr 13 17:47:29 2015
## 20 percent done (iteration 2200) Mon Apr 13 17:47:33 2015
## 25 percent done (iteration 2750) Mon Apr 13 17:47:37 2015
## 30 percent done (iteration 3300) Mon Apr 13 17:47:40 2015
## 35 percent done (iteration 3850) Mon Apr 13 17:47:44 2015
## 40 percent done (iteration 4400) Mon Apr 13 17:47:47 2015
## 45 percent done (iteration 4950) Mon Apr 13 17:47:51 2015
## 50 percent done (iteration 5500) Mon Apr 13 17:47:54 2015
## 55 percent done (iteration 6050) Mon Apr 13 17:47:58 2015
## 60 percent done (iteration 6600) Mon Apr 13 17:48:02 2015
## 65 percent done (iteration 7150) Mon Apr 13 17:48:05 2015
## 70 percent done (iteration 7700) Mon Apr 13 17:48:09 2015
## 75 percent done (iteration 8250) Mon Apr 13 17:48:12 2015
## 80 percent done (iteration 8800) Mon Apr 13 17:48:16 2015
## 85 percent done (iteration 9350) Mon Apr 13 17:48:19 2015
## 90 percent done (iteration 9900) Mon Apr 13 17:48:23 2015
## 95 percent done (iteration 10450) Mon Apr 13 17:48:27 2015
## 100 percent done (iteration 11000) Mon Apr 13 17:48:31 2015
#R : the rank of the approximating factor matrix
#S : number of samples from the Markov chain
#burn : number of initial scans of the Markov chain to be dropped
(ii)a covariate for common practice
In order to include the effects of common practice, we create an array with that information
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)
#벡 터 1,4,9 가 있으면 1 없으면 0으로 인접행렬만들기
same.prac <- array(same.prac,dim=c(36, 36, 1))
#and fit the model with this addtional argument
lazega.leig.fit2 <- eigenmodel_mcmc(A, same.prac, R=2,S=11000,burn=10000)
## 5 percent done (iteration 550) Mon Apr 13 17:48:35 2015
## 10 percent done (iteration 1100) Mon Apr 13 17:48:38 2015
## 15 percent done (iteration 1650) Mon Apr 13 17:48:42 2015
## 20 percent done (iteration 2200) Mon Apr 13 17:48:46 2015
## 25 percent done (iteration 2750) Mon Apr 13 17:48:50 2015
## 30 percent done (iteration 3300) Mon Apr 13 17:48:53 2015
## 35 percent done (iteration 3850) Mon Apr 13 17:48:57 2015
## 40 percent done (iteration 4400) Mon Apr 13 17:49:01 2015
## 45 percent done (iteration 4950) Mon Apr 13 17:49:04 2015
## 50 percent done (iteration 5500) Mon Apr 13 17:49:08 2015
## 55 percent done (iteration 6050) Mon Apr 13 17:49:12 2015
## 60 percent done (iteration 6600) Mon Apr 13 17:49:16 2015
## 65 percent done (iteration 7150) Mon Apr 13 17:49:19 2015
## 70 percent done (iteration 7700) Mon Apr 13 17:49:23 2015
## 75 percent done (iteration 8250) Mon Apr 13 17:49:27 2015
## 80 percent done (iteration 8800) Mon Apr 13 17:49:30 2015
## 85 percent done (iteration 9350) Mon Apr 13 17:49:34 2015
## 90 percent done (iteration 9900) Mon Apr 13 17:49:38 2015
## 95 percent done (iteration 10450) Mon Apr 13 17:49:42 2015
## 100 percent done (iteration 11000) Mon Apr 13 17:49:46 2015
(iii)a covariate for shared office location
Finally, we do similarly for the model that includes information on shared office locations.
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)
## 5 percent done (iteration 550) Mon Apr 13 17:49:50 2015
## 10 percent done (iteration 1100) Mon Apr 13 17:49:53 2015
## 15 percent done (iteration 1650) Mon Apr 13 17:49:57 2015
## 20 percent done (iteration 2200) Mon Apr 13 17:50:00 2015
## 25 percent done (iteration 2750) Mon Apr 13 17:50:04 2015
## 30 percent done (iteration 3300) Mon Apr 13 17:50:08 2015
## 35 percent done (iteration 3850) Mon Apr 13 17:50:11 2015
## 40 percent done (iteration 4400) Mon Apr 13 17:50:15 2015
## 45 percent done (iteration 4950) Mon Apr 13 17:50:19 2015
## 50 percent done (iteration 5500) Mon Apr 13 17:50:22 2015
## 55 percent done (iteration 6050) Mon Apr 13 17:50:26 2015
## 60 percent done (iteration 6600) Mon Apr 13 17:50:29 2015
## 65 percent done (iteration 7150) Mon Apr 13 17:50:33 2015
## 70 percent done (iteration 7700) Mon Apr 13 17:50:37 2015
## 75 percent done (iteration 8250) Mon Apr 13 17:50:40 2015
## 80 percent done (iteration 8800) Mon Apr 13 17:50:44 2015
## 85 percent done (iteration 9350) Mon Apr 13 17:50:48 2015
## 90 percent done (iteration 9900) Mon Apr 13 17:50:51 2015
## 95 percent done (iteration 10450) Mon Apr 13 17:50:55 2015
## 100 percent done (iteration 11000) Mon Apr 13 17:50:59 2015
?eigenmodel_mcmc
## starting httpd help server ... done
In order to compare the representation of the network lazega in each of the underlying two-dimensional latent spaces inferred for these models, we extract the eigenvectors for each fitted model.
(이 모델에서 추정되는 2차원latent spaces의 lazega의 표현을 위해 ,우리는 각각의 맞는 모델의 고유 벡터를 추출합니다.)
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]
and plot the network in igraph using these coordinates as the layout.
(위 레이아웃의 좌표를 이용해 네트워크를 그려라)
Visualizations of network of Lazega’s lawyers, with layout determined according to the inferred latent eigenvectors in model with no pair-specific covariates(letf), acovariate for common practice(center), and a covariate for shared office location(right)
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)
par(mfrow=c(1,3))
plot(lazega, layout=lat.sp.1, vertex.color=v.colors,vertex.shape=v.shapes, vertex.size=v.size)
plot(lazega, layout=lat.sp.2, vertex.color=v.colors,vertex.shape=v.shapes, vertex.size=v.size)
plot(lazega, layout=lat.sp.3, vertex.color=v.colors,vertex.shape=v.shapes, vertex.size=v.size)
the first two are somewhat similar, the third is distinct. first two visualizations appear to be clustered into two main groups distinguished largely by common office lacation(i.e. color), in the third there appears to be only one main cluster.
(첫번째 두 그래프는 비슷하나 세번째는 다소 유사하다.)
(첫번째 두 시각화는 오피스지역에 따라 크게 두개로 분류되나 3번째로 오직한개의 클러스터로 나타난다. )
first two models there is one eigenvalue that clearly dominates the other, corresponding to the axis on which we obtain a clear separation between the two groups.
(처음 두 모델을 명확하게 한개의 고유값이 다른거에 압도적으로 우세하여 축에 대응하는 명확한 부리는 두그룹사이에서 얻을 수 있다.)
whereas for the third model the eigenvalues are comparable in their magnitude.
(반면에 3번째모델의 고유값은 크기가 비슷하다.)
apply(lazega.leig.fit1$L_postsamp, 2, mean)
## [1] 0.2603655 1.0384032
apply(lazega.leig.fit2$L_postsamp, 2, mean)
## [1] 0.9083401 -0.1385321
apply(lazega.leig.fit3$L_postsamp, 2, mean)
## [1] 0.5970403 0.3112896
6.4.4. Goodness-of-Fit
Use the same types of methods we employed in the analysis of ERGMs.
For example, consider the model fit above to the data Lazega with no pair-specific covariates.
After initiating a permutation of the 36*35/2=630 unique off-diagonal elements of the symmetric adjacency matrix We can do the same for the models fit above with pair-specific covariates for common practice and shared office location, respectively, yielding,say,
We can do the same for the models fit above with pair-specific covariates for common practice and shared office location , yielding , say, Avec.pred2 and Avec.pred3.
The result of the predictions generated under each of these three models can be assessed by examination of the corresponding receiver operating characteristic (ROC) curves.
(이러한 세 가지 모델의 각각의 아래에 생성 된 예측 결과는 특성 (ROC) 곡선을 조작하고 해당 수신기의 시험에 의해 평가 될 수있다.)