HW05
JUNG KYU HYEONG
2020 6 21
1. Implement all codes for latent network model.
library(latentnet)## Loading required package: network## network: Classes for Relational Data
## Version 1.16.0 created on 2019-11-30.
## 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.## Loading required package: ergm##
## ergm: version 3.10.4, created on 2019-06-10
## Copyright (c) 2019, 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
## Skye Bender-deMoll, University of Washington
## Chad Klumb
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("ergm").## NOTE: Versions before 3.6.1 had a bug in the implementation of the bd()
## constriant which distorted the sampled distribution somewhat. In
## addition, Sampson's Monks datasets had mislabeled vertices. See the
## NEWS and the documentation for more details.## NOTE: Some common term arguments pertaining to vertex attribute and
## level selection have changed in 3.10.0. See terms help for more
## details. Use 'options(ergm.term=list(version="3.9.4"))' to use old
## behavior.##
## latentnet: version 2.10.5, created on 2020-03-20
## Copyright (c) 2020, Pavel N. Krivitsky
## Mark S. Handcock, University of California -- Los Angeles
## with contributions from
## Susan M. Shortreed
## Jeremy Tantrum
## Peter D. Hoff
## Li Wang
## Kirk Li, University of Washington
## Jake Fisher
## Jordan T. Bates
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("latentnet").
## NOTE: BIC calculation prior to latentnet 2.7.0 had a bug in the calculation of the effective number of parameters. See help(summary.ergmm) for details.
## NOTE: Prior to version 2.8.0, handling of fixed effects for directed networks had a bug. See help("ergmm-terms") for details.data("sampson")
network.vertex.names(samplike)## [1] "John Bosco" "Gregory" "Basil" "Peter" "Bonaventure"
## [6] "Berthold" "Mark" "Victor" "Ambrose" "Romauld"
## [11] "Louis" "Winfrid" "Amand" "Hugh" "Boniface"
## [16] "Albert" "Elias" "Simplicius"samplike %v% "group"## [1] "Turks" "Turks" "Outcasts" "Loyal" "Loyal" "Loyal"
## [7] "Turks" "Loyal" "Loyal" "Loyal" "Loyal" "Turks"
## [13] "Outcasts" "Turks" "Turks" "Turks" "Outcasts" "Outcasts"samplike.fit=ergmm(samplike ~ euclidean(d=2),verbose=TRUE)## Generating initial values for MCMC:
## Computing geodesic distances... Finished.
## Computing MDS locations... Finished.
## Computing other initial values... Finished.
## Finding the conditional posterior mode... Finished.
## Burning in... Finished.
## Starting sampling run... Finished.
## Post-processing the MCMC output:
## Fitting the MKL locations... Finished.
## MKL MBC is not available or non-latent-cluster model.
## Performing Procrustes transformation... Finished.#We fit samplike network data using latent space model with 2-dim Euclidean space#if we use euclidean(d=2),we do not add any group structure in latent space model
#if we use euclidean(d=2),we do not add any group structure in latent position of vertices
#if we use euclidean(d=2,G=3),we add a group structure in latent position of vertices
#do not use d=1
#because we do not have any covariate, we only consider the intercept term.
#logit(pi_{i})=\alpha-|z_i - z_j|
#when we have covariate, we need to change them to the matrix form by using either main
#effect, homophily effect or absolute difference
#logit(pi_{i})=\alpha+\beta x_{ij}-|z_i - z_j|samplike.fit$mcmc.pmode$Z## [,1] [,2]
## [1,] 0.3603455 0.07665688
## [2,] -0.3362347 1.28367030
## [3,] -1.6434457 -1.60823125
## [4,] 2.4307748 -0.59872128
## [5,] 0.9203736 -0.30627673
## [6,] 1.9041286 -0.59138112
## [7,] -0.6730192 1.39097225
## [8,] 1.6677833 -1.32844644
## [9,] 1.0783131 -1.57407550
## [10,] 1.6574102 -1.33609880
## [11,] 2.7711144 0.73853083
## [12,] -0.5425823 1.59002931
## [13,] -0.6589592 -1.15220444
## [14,] 0.4301689 1.73783700
## [15,] 0.1482861 1.80861789
## [16,] -0.9999788 1.50469586
## [17,] -1.9100392 -1.31141963
## [18,] -1.4461493 -1.15082996#Z is latent positiononeL=samplike %v% "group"
oneL[oneL=="Turks"]="T"
oneL[oneL=="Outcasts"]="O"
oneL[oneL=="Loyal"]="L"
oneL[c(1,7,15)]="W"
oneL## [1] "W" "T" "O" "L" "L" "L" "W" "L" "L" "L" "L" "T" "O" "T" "W" "T" "O" "O"oneLcolors=c("red","blue","black","green")[match(oneL, c("T","O","L","W"))]
plot(samplike.fit,label=oneL,vertex.col = oneLcolors,
what="pmean",main="MCMC positions",print.formula = FALSE,labels=TRUE)
#due to the invariance property of the latent space, positions can be rotated
#however, if you rotate your latent space plot, you will have same plot
#that is just the difference in representationModel Fitting: Eigenmodel
library(sand)## Loading required package: igraph##
## Attaching package: 'igraph'## The following objects are masked from 'package:network':
##
## %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## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## union## Loading required package: igraphdata##
## Statistical Analysis of Network Data with R
## Type in C2 (+ENTER) to start with Chapter 2.data(lazega)
summary(lazega)## This graph was created by an old(er) igraph version.
## Call upgrade_graph() on it to use with the current igraph version
## For now we convert it on the fly...## IGRAPH NA 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)library(eigenmodel)
lazega=upgrade_graph(lazega)
A=get.adjacency(lazega,sparse=FALSE)
lazega.leig.fit1=eigenmodel_mcmc(A, R=2, S=10000, burn=10000)## 5 percent done (iteration 1000) Mon Jun 22 00:12:49 2020
## 10 percent done (iteration 2000) Mon Jun 22 00:12:55 2020
## 15 percent done (iteration 3000) Mon Jun 22 00:13:00 2020
## 20 percent done (iteration 4000) Mon Jun 22 00:13:06 2020
## 25 percent done (iteration 5000) Mon Jun 22 00:13:12 2020
## 30 percent done (iteration 6000) Mon Jun 22 00:13:17 2020
## 35 percent done (iteration 7000) Mon Jun 22 00:13:23 2020
## 40 percent done (iteration 8000) Mon Jun 22 00:13:29 2020
## 45 percent done (iteration 9000) Mon Jun 22 00:13:34 2020
## 50 percent done (iteration 10000) Mon Jun 22 00:13:40 2020
## 55 percent done (iteration 11000) Mon Jun 22 00:13:46 2020
## 60 percent done (iteration 12000) Mon Jun 22 00:13:51 2020
## 65 percent done (iteration 13000) Mon Jun 22 00:13:57 2020
## 70 percent done (iteration 14000) Mon Jun 22 00:14:03 2020
## 75 percent done (iteration 15000) Mon Jun 22 00:14:08 2020
## 80 percent done (iteration 16000) Mon Jun 22 00:14:14 2020
## 85 percent done (iteration 17000) Mon Jun 22 00:14:22 2020
## 90 percent done (iteration 18000) Mon Jun 22 00:14:29 2020
## 95 percent done (iteration 19000) Mon Jun 22 00:14:35 2020
## 100 percent done (iteration 20000) Mon Jun 22 00:14:42 2020#R=dimension
#This model is without any covariatesame.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=10000, burn=10000)## 5 percent done (iteration 1000) Mon Jun 22 00:14:50 2020
## 10 percent done (iteration 2000) Mon Jun 22 00:14:56 2020
## 15 percent done (iteration 3000) Mon Jun 22 00:15:03 2020
## 20 percent done (iteration 4000) Mon Jun 22 00:15:09 2020
## 25 percent done (iteration 5000) Mon Jun 22 00:15:16 2020
## 30 percent done (iteration 6000) Mon Jun 22 00:15:22 2020
## 35 percent done (iteration 7000) Mon Jun 22 00:15:28 2020
## 40 percent done (iteration 8000) Mon Jun 22 00:15:34 2020
## 45 percent done (iteration 9000) Mon Jun 22 00:15:40 2020
## 50 percent done (iteration 10000) Mon Jun 22 00:15:47 2020
## 55 percent done (iteration 11000) Mon Jun 22 00:15:53 2020
## 60 percent done (iteration 12000) Mon Jun 22 00:15:59 2020
## 65 percent done (iteration 13000) Mon Jun 22 00:16:05 2020
## 70 percent done (iteration 14000) Mon Jun 22 00:16:11 2020
## 75 percent done (iteration 15000) Mon Jun 22 00:16:17 2020
## 80 percent done (iteration 16000) Mon Jun 22 00:16:23 2020
## 85 percent done (iteration 17000) Mon Jun 22 00:16:30 2020
## 90 percent done (iteration 18000) Mon Jun 22 00:16:36 2020
## 95 percent done (iteration 19000) Mon Jun 22 00:16:42 2020
## 100 percent done (iteration 20000) Mon Jun 22 00:16:48 2020#Model2:with a covariate for common practice
#same.prac covariate informationsame.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=10000, burn=10000)## 5 percent done (iteration 1000) Mon Jun 22 00:16:54 2020
## 10 percent done (iteration 2000) Mon Jun 22 00:17:00 2020
## 15 percent done (iteration 3000) Mon Jun 22 00:17:06 2020
## 20 percent done (iteration 4000) Mon Jun 22 00:17:12 2020
## 25 percent done (iteration 5000) Mon Jun 22 00:17:18 2020
## 30 percent done (iteration 6000) Mon Jun 22 00:17:25 2020
## 35 percent done (iteration 7000) Mon Jun 22 00:17:31 2020
## 40 percent done (iteration 8000) Mon Jun 22 00:17:37 2020
## 45 percent done (iteration 9000) Mon Jun 22 00:17:43 2020
## 50 percent done (iteration 10000) Mon Jun 22 00:17:49 2020
## 55 percent done (iteration 11000) Mon Jun 22 00:17:56 2020
## 60 percent done (iteration 12000) Mon Jun 22 00:18:04 2020
## 65 percent done (iteration 13000) Mon Jun 22 00:18:11 2020
## 70 percent done (iteration 14000) Mon Jun 22 00:18:18 2020
## 75 percent done (iteration 15000) Mon Jun 22 00:18:24 2020
## 80 percent done (iteration 16000) Mon Jun 22 00:18:30 2020
## 85 percent done (iteration 17000) Mon Jun 22 00:18:37 2020
## 90 percent done (iteration 18000) Mon Jun 22 00:18:43 2020
## 95 percent done (iteration 19000) Mon Jun 22 00:18:48 2020
## 100 percent done (iteration 20000) Mon Jun 22 00:18:54 2020#Model 3: with a covariate for shared office locationIn 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.
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]
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
par(mfrow=c(1,3))
plot(lazega,layout=lat.sp.1,vertex.color=v.colors,vertex.shape=v.shapes, vertex.size=v.size,vertex.label=v.label)
plot(lazega,layout=lat.sp.2,vertex.color=v.colors,vertex.shape=v.shapes, vertex.size=v.size,vertex.label=v.label)
plot(lazega,layout=lat.sp.3,vertex.color=v.colors,vertex.shape=v.shapes, vertex.size=v.size,vertex.label=v.label)
#The first 2 visualization appear to be clustered into two groups distinguisher largely common office location(check color)
#The last visualization appear to be only one main cluster
#that means we can explain thie network based on not office location but practice
#If we try to explain the network usin office location, then it is less structured
#Also, we can find how nuch we can explain the network using eigenvalueapply(lazega.leig.fit1$L_postsamp,2,mean)## [1] 0.9993555 0.2324392apply(lazega.leig.fit2$L_postsamp,2,mean)## [1] 0.9097744 -0.1158853apply(lazega.leig.fit3$L_postsamp,2,mean)## [1] -0.1066956 0.3585928#goodness-of-fit
##To evaluate the goodness-of-fit,you can use ROC curve
perm.index = sample(1:630)
nfolds = 5
nmiss = 630/nfolds
Avec = A[lower.tri(A)]
Avec.pred1 =numeric(length(Avec))
for(i in 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
model.fit1 = eigenmodel_mcmc(Y, R=2, S=10000, burn=10000)
model.pred1 = model.fit1$Y_postmean
model.pred.vec1 = model.pred1[lower.tri(model.pred1)]
Avec.pred1[A.miss.index] = model.pred.vec1[A.miss.index]
}## 5 percent done (iteration 1000) Mon Jun 22 00:19:01 2020
## 10 percent done (iteration 2000) Mon Jun 22 00:19:07 2020
## 15 percent done (iteration 3000) Mon Jun 22 00:19:13 2020
## 20 percent done (iteration 4000) Mon Jun 22 00:19:19 2020
## 25 percent done (iteration 5000) Mon Jun 22 00:19:25 2020
## 30 percent done (iteration 6000) Mon Jun 22 00:19:30 2020
## 35 percent done (iteration 7000) Mon Jun 22 00:19:36 2020
## 40 percent done (iteration 8000) Mon Jun 22 00:19:42 2020
## 45 percent done (iteration 9000) Mon Jun 22 00:19:48 2020
## 50 percent done (iteration 10000) Mon Jun 22 00:19:54 2020
## 55 percent done (iteration 11000) Mon Jun 22 00:19:59 2020
## 60 percent done (iteration 12000) Mon Jun 22 00:20:05 2020
## 65 percent done (iteration 13000) Mon Jun 22 00:20:11 2020
## 70 percent done (iteration 14000) Mon Jun 22 00:20:17 2020
## 75 percent done (iteration 15000) Mon Jun 22 00:20:23 2020
## 80 percent done (iteration 16000) Mon Jun 22 00:20:28 2020
## 85 percent done (iteration 17000) Mon Jun 22 00:20:34 2020
## 90 percent done (iteration 18000) Mon Jun 22 00:20:40 2020
## 95 percent done (iteration 19000) Mon Jun 22 00:20:46 2020
## 100 percent done (iteration 20000) Mon Jun 22 00:20:52 2020
## 5 percent done (iteration 1000) Mon Jun 22 00:20:58 2020
## 10 percent done (iteration 2000) Mon Jun 22 00:21:04 2020
## 15 percent done (iteration 3000) Mon Jun 22 00:21:09 2020
## 20 percent done (iteration 4000) Mon Jun 22 00:21:15 2020
## 25 percent done (iteration 5000) Mon Jun 22 00:21:21 2020
## 30 percent done (iteration 6000) Mon Jun 22 00:21:27 2020
## 35 percent done (iteration 7000) Mon Jun 22 00:21:33 2020
## 40 percent done (iteration 8000) Mon Jun 22 00:21:38 2020
## 45 percent done (iteration 9000) Mon Jun 22 00:21:44 2020
## 50 percent done (iteration 10000) Mon Jun 22 00:21:50 2020
## 55 percent done (iteration 11000) Mon Jun 22 00:21:56 2020
## 60 percent done (iteration 12000) Mon Jun 22 00:22:02 2020
## 65 percent done (iteration 13000) Mon Jun 22 00:22:08 2020
## 70 percent done (iteration 14000) Mon Jun 22 00:22:14 2020
## 75 percent done (iteration 15000) Mon Jun 22 00:22:20 2020
## 80 percent done (iteration 16000) Mon Jun 22 00:22:25 2020
## 85 percent done (iteration 17000) Mon Jun 22 00:22:31 2020
## 90 percent done (iteration 18000) Mon Jun 22 00:22:37 2020
## 95 percent done (iteration 19000) Mon Jun 22 00:22:43 2020
## 100 percent done (iteration 20000) Mon Jun 22 00:22:49 2020
## 5 percent done (iteration 1000) Mon Jun 22 00:22:55 2020
## 10 percent done (iteration 2000) Mon Jun 22 00:23:00 2020
## 15 percent done (iteration 3000) Mon Jun 22 00:23:06 2020
## 20 percent done (iteration 4000) Mon Jun 22 00:23:12 2020
## 25 percent done (iteration 5000) Mon Jun 22 00:23:18 2020
## 30 percent done (iteration 6000) Mon Jun 22 00:23:24 2020
## 35 percent done (iteration 7000) Mon Jun 22 00:23:29 2020
## 40 percent done (iteration 8000) Mon Jun 22 00:23:35 2020
## 45 percent done (iteration 9000) Mon Jun 22 00:23:41 2020
## 50 percent done (iteration 10000) Mon Jun 22 00:23:47 2020
## 55 percent done (iteration 11000) Mon Jun 22 00:23:52 2020
## 60 percent done (iteration 12000) Mon Jun 22 00:23:58 2020
## 65 percent done (iteration 13000) Mon Jun 22 00:24:04 2020
## 70 percent done (iteration 14000) Mon Jun 22 00:24:10 2020
## 75 percent done (iteration 15000) Mon Jun 22 00:24:16 2020
## 80 percent done (iteration 16000) Mon Jun 22 00:24:22 2020
## 85 percent done (iteration 17000) Mon Jun 22 00:24:28 2020
## 90 percent done (iteration 18000) Mon Jun 22 00:24:34 2020
## 95 percent done (iteration 19000) Mon Jun 22 00:24:40 2020
## 100 percent done (iteration 20000) Mon Jun 22 00:24:45 2020
## 5 percent done (iteration 1000) Mon Jun 22 00:24:51 2020
## 10 percent done (iteration 2000) Mon Jun 22 00:24:57 2020
## 15 percent done (iteration 3000) Mon Jun 22 00:25:03 2020
## 20 percent done (iteration 4000) Mon Jun 22 00:25:09 2020
## 25 percent done (iteration 5000) Mon Jun 22 00:25:14 2020
## 30 percent done (iteration 6000) Mon Jun 22 00:25:20 2020
## 35 percent done (iteration 7000) Mon Jun 22 00:25:26 2020
## 40 percent done (iteration 8000) Mon Jun 22 00:25:32 2020
## 45 percent done (iteration 9000) Mon Jun 22 00:25:38 2020
## 50 percent done (iteration 10000) Mon Jun 22 00:25:43 2020
## 55 percent done (iteration 11000) Mon Jun 22 00:25:50 2020
## 60 percent done (iteration 12000) Mon Jun 22 00:25:56 2020
## 65 percent done (iteration 13000) Mon Jun 22 00:26:02 2020
## 70 percent done (iteration 14000) Mon Jun 22 00:26:08 2020
## 75 percent done (iteration 15000) Mon Jun 22 00:26:14 2020
## 80 percent done (iteration 16000) Mon Jun 22 00:26:20 2020
## 85 percent done (iteration 17000) Mon Jun 22 00:26:26 2020
## 90 percent done (iteration 18000) Mon Jun 22 00:26:32 2020
## 95 percent done (iteration 19000) Mon Jun 22 00:26:38 2020
## 100 percent done (iteration 20000) Mon Jun 22 00:26:44 2020
## 5 percent done (iteration 1000) Mon Jun 22 00:26:49 2020
## 10 percent done (iteration 2000) Mon Jun 22 00:26:55 2020
## 15 percent done (iteration 3000) Mon Jun 22 00:27:01 2020
## 20 percent done (iteration 4000) Mon Jun 22 00:27:07 2020
## 25 percent done (iteration 5000) Mon Jun 22 00:27:13 2020
## 30 percent done (iteration 6000) Mon Jun 22 00:27:18 2020
## 35 percent done (iteration 7000) Mon Jun 22 00:27:24 2020
## 40 percent done (iteration 8000) Mon Jun 22 00:27:30 2020
## 45 percent done (iteration 9000) Mon Jun 22 00:27:36 2020
## 50 percent done (iteration 10000) Mon Jun 22 00:27:42 2020
## 55 percent done (iteration 11000) Mon Jun 22 00:27:48 2020
## 60 percent done (iteration 12000) Mon Jun 22 00:27:53 2020
## 65 percent done (iteration 13000) Mon Jun 22 00:27:59 2020
## 70 percent done (iteration 14000) Mon Jun 22 00:28:05 2020
## 75 percent done (iteration 15000) Mon Jun 22 00:28:11 2020
## 80 percent done (iteration 16000) Mon Jun 22 00:28:17 2020
## 85 percent done (iteration 17000) Mon Jun 22 00:28:23 2020
## 90 percent done (iteration 18000) Mon Jun 22 00:28:28 2020
## 95 percent done (iteration 19000) Mon Jun 22 00:28:34 2020
## 100 percent done (iteration 20000) Mon Jun 22 00:28:40 2020library(ROCR)
pred1 = prediction(Avec.pred1, Avec)
perf1 = performance(pred1, 'tpr', 'fpr')
par(mfrow=c(1,1))
plot(perf1, col='blue', lwd=3)
perf1.auc = performance(pred1, 'auc')
slot(perf1.auc, 'y.values')## [[1]]
## [1] 0.80459272. Fit Karate network using the latent space model. (only include intercept + distance term, no covariates)
data("karate")
karate=upgrade_graph(karate)
V(karate)## + 34/34 vertices, named, from 4b458a1:
## [1] Mr Hi Actor 2 Actor 3 Actor 4 Actor 5 Actor 6 Actor 7 Actor 8
## [9] Actor 9 Actor 10 Actor 11 Actor 12 Actor 13 Actor 14 Actor 15 Actor 16
## [17] Actor 17 Actor 18 Actor 19 Actor 20 Actor 21 Actor 22 Actor 23 Actor 24
## [25] Actor 25 Actor 26 Actor 27 Actor 28 Actor 29 Actor 30 Actor 31 Actor 32
## [33] Actor 33 John AA=get.adjacency(karate,sparse=F)
v.attrs=get.data.frame(karate,what="vertices")
v.attrs## Faction name label color
## Mr Hi 1 Mr Hi H 1
## Actor 2 1 Actor 2 2 1
## Actor 3 1 Actor 3 3 1
## Actor 4 1 Actor 4 4 1
## Actor 5 1 Actor 5 5 1
## Actor 6 1 Actor 6 6 1
## Actor 7 1 Actor 7 7 1
## Actor 8 1 Actor 8 8 1
## Actor 9 2 Actor 9 9 2
## Actor 10 2 Actor 10 10 2
## Actor 11 1 Actor 11 11 1
## Actor 12 1 Actor 12 12 1
## Actor 13 1 Actor 13 13 1
## Actor 14 1 Actor 14 14 1
## Actor 15 2 Actor 15 15 2
## Actor 16 2 Actor 16 16 2
## Actor 17 1 Actor 17 17 1
## Actor 18 1 Actor 18 18 1
## Actor 19 2 Actor 19 19 2
## Actor 20 1 Actor 20 20 1
## Actor 21 2 Actor 21 21 2
## Actor 22 1 Actor 22 22 1
## Actor 23 2 Actor 23 23 2
## Actor 24 2 Actor 24 24 2
## Actor 25 2 Actor 25 25 2
## Actor 26 2 Actor 26 26 2
## Actor 27 2 Actor 27 27 2
## Actor 28 2 Actor 28 28 2
## Actor 29 2 Actor 29 29 2
## Actor 30 2 Actor 30 30 2
## Actor 31 2 Actor 31 31 2
## Actor 32 2 Actor 32 32 2
## Actor 33 2 Actor 33 33 2
## John A 2 John A A 2karate.s=network::as.network(as.matrix(A),directed=F)
network::set.vertex.attribute(karate.s,'Faction',v.attrs$Faction)
network::set.vertex.attribute(karate.s,'color',v.attrs$color)
karate.s %v% 'Faction'## [1] 1 1 1 1 1 1 1 1 2 2 1 1 1 1 2 2 1 1 2 1 2 1 2 2 2 2 2 2 2 2 2 2 2 2we can see faction value is 1 or 2
karate.fit = ergmm(karate.s ~ euclidean(d=2, G=2), verbose=T)## Generating initial values for MCMC:
## Computing geodesic distances... Finished.
## Computing MDS locations... Finished.
## Computing other initial values... Finished.
## Finding the conditional posterior mode... Finished.
## Burning in... Finished.
## Starting sampling run... Finished.
## Post-processing the MCMC output:
## Performing label-switching... Finished.
## Fitting the MKL locations... Finished.
## Fitting MBC conditional on MKL locations... Finished.
## Performing Procrustes transformation... Finished.We fit karate network data using latent space model with 2-dim Euclidean space
because of number of faction, we add a group structure in latent position of vertices
karate.fit$mcmc.pmode$Z## [,1] [,2]
## [1,] -1.0447232 0.72976993
## [2,] -0.5926773 1.31309443
## [3,] 0.7113302 1.66272851
## [4,] -0.5845566 0.71264067
## [5,] -1.4676213 0.14515291
## [6,] -0.7403098 -0.60881497
## [7,] -0.9822659 0.03743523
## [8,] -0.8332010 1.16203878
## [9,] 0.8154492 2.57840617
## [10,] 1.8417689 3.94521898
## [11,] -1.4504418 0.38880470
## [12,] -1.4309960 2.00660736
## [13,] -1.4692962 0.81176987
## [14,] -0.3149204 1.04292082
## [15,] 3.0879321 2.78874198
## [16,] 0.5145647 3.93651398
## [17,] -1.9486429 -0.40507666
## [18,] -1.8232886 1.55630065
## [19,] 2.0417128 2.18538884
## [20,] -1.2740561 1.71916484
## [21,] 1.3883340 3.29265917
## [22,] -2.2262300 0.20006291
## [23,] 0.1013299 2.69877504
## [24,] 2.2163109 4.23968080
## [25,] 3.3135998 2.61418833
## [26,] 3.3425219 2.40474163
## [27,] 1.8925817 2.87415218
## [28,] 3.3120627 2.92660763
## [29,] 0.9392122 2.63704475
## [30,] 2.6786600 4.11057074
## [31,] 0.9334938 2.89583995
## [32,] 2.4302593 2.80379074
## [33,] 1.6459879 2.97538751
## [34,] 1.6766565 2.91086033oneL=karate.s %v% 'Faction'
oneL[oneL==1]='H'
oneL[oneL==2]='J'
oneL[[1]] = "Hi"
oneL[[34]] = "John"
oneL## [1] "Hi" "H" "H" "H" "H" "H" "H" "H" "J" "J"
## [11] "H" "H" "H" "H" "J" "J" "H" "H" "J" "H"
## [21] "J" "H" "J" "J" "J" "J" "J" "J" "J" "J"
## [31] "J" "J" "J" "John"oneLcolors = c("red", "blue", "yellow", "white")[match(oneL, c("H","J","Hi","John"))]
plot(karate.fit, label=oneL, vertex.col=oneLcolors, what="pmean",
main="MCMC position", print.formula = FALSE, labels=TRUE)
due to the invariance property of the latent space, positions can be rotated
however, if you rotate your latent space plot, you will have same plot
that is just the difference in representation.
we only consider the only include intercept + distance term (no covariates)
1.Hi and John located in each group center
2.we can check that each 2 groups are divided well in the plot
3.we can see that groups are formed according to faction
3. Fit Karate network using model based eigen decomposition model. (intercept + Z)
summary(karate)## IGRAPH 4b458a1 UNW- 34 78 -- Zachary's karate club network
## + attr: name (g/c), Citation (g/c), Author (g/c), Faction (v/n), name
## | (v/c), label (v/c), color (v/n), weight (e/n)library(eigenmodel)
karate.fit1=eigenmodel_mcmc(A,R=2,S=10000,burn=10000)## 5 percent done (iteration 1000) Mon Jun 22 00:29:05 2020
## 10 percent done (iteration 2000) Mon Jun 22 00:29:11 2020
## 15 percent done (iteration 3000) Mon Jun 22 00:29:16 2020
## 20 percent done (iteration 4000) Mon Jun 22 00:29:21 2020
## 25 percent done (iteration 5000) Mon Jun 22 00:29:27 2020
## 30 percent done (iteration 6000) Mon Jun 22 00:29:32 2020
## 35 percent done (iteration 7000) Mon Jun 22 00:29:38 2020
## 40 percent done (iteration 8000) Mon Jun 22 00:29:44 2020
## 45 percent done (iteration 9000) Mon Jun 22 00:29:49 2020
## 50 percent done (iteration 10000) Mon Jun 22 00:29:55 2020
## 55 percent done (iteration 11000) Mon Jun 22 00:30:00 2020
## 60 percent done (iteration 12000) Mon Jun 22 00:30:06 2020
## 65 percent done (iteration 13000) Mon Jun 22 00:30:11 2020
## 70 percent done (iteration 14000) Mon Jun 22 00:30:17 2020
## 75 percent done (iteration 15000) Mon Jun 22 00:30:22 2020
## 80 percent done (iteration 16000) Mon Jun 22 00:30:28 2020
## 85 percent done (iteration 17000) Mon Jun 22 00:30:33 2020
## 90 percent done (iteration 18000) Mon Jun 22 00:30:39 2020
## 95 percent done (iteration 19000) Mon Jun 22 00:30:45 2020
## 100 percent done (iteration 20000) Mon Jun 22 00:30:50 2020basic setting is same as problem2
Using eigen decomposition model, fit karate netwrok (intercept+Z)
lat.sp = eigen(karate.fit1$ULU_postmean)$vec[,1:2]
plot(karate, layout=lat.sp, vertex.color=oneLcolors)
1.we can check that each 2 groups are divided well in the plot
2.we can see that groups are formed according to faction
3.Unlike problem2, Hi and John located in each the corner
apply(karate.fit1$L_postsamp, 2, mean)## [1] -0.7608995 1.2522484there is one eigenvalue that clearly dominates the other, corresponding to the axis on which we obtain a clear separation between the two groups.#in textbook
Unlike in PCA, the eigenvalues are not of decreasing order since they were estimated using MCMC
We see that the second eigenvalue dominates the solution. Thus, the second latent variable is relatively more important than the first one # Modern Psychometrics with R