Resources: introToSNAinR_sunbelt_2012.R

File introToSNAinR_sunbelt_2012.R, 22.4 KB (added by morrism, 4 years ago)

tutorial code for workshop 1. Intro

Line 
1
2# Preliminaries
3
4## Download and attaching \statnet and associated packages
5##      First, open \R.
6##      Then install the \statnet installer. At the \R cursor, type:
7       
8install.packages("statnet")
9
10## Now, and in the future, you can install/update \statnet at any point using:
11
12update.packages("statnet")
13
14## Attach \statnet to your \R session by typing:
15
16library(statnet)
17
18## Download and install supplemental packages
19
20## The rgl package sometimes has problems installing on specific platforms. 
21## If you can't install it, don't worry---you need it for 3D network visualization,
22## but not for anything else.  Consult the \R website for more information on these or other contributed packages.
23
24install.packages("rgl")
25
26
27## Set a specific working directory for this tutorial if you wish
28##      If you are using Windows, go to File > Change Dir and choose the directory.
29##      On a Mac, go to Misc > Change Working Directory and choose the directory.
30##      Or, you can just use the setwd() command:
31
32setwd("full.path.for.the.folder")
33
34
35######################################################################################
36# Introduction to basic \R syntax
37
381 + 3 # evaluation
39a <- 3 # assignment
40a # evaluation
41a<-3 # spacing does not matter
42a <- 3 # spacing does not matter
43sqrt(a) # use the square root function
44b <- sqrt(a) # use function and save result
45b
46c # evaluate something that is not there
47a == b # is a equivalent to b?
48a != b # is a not equal to b?
49ls() # list objects in the global environment
50
51help.start() # get help with R generally
52?sqrt # get specific help for a function
53??sqrt # looking for help pertaining to sqrt
54apropos("sq") # it's on the tip of my tongue...
55
56rm(a) # remove a single object
57rm(list=ls()) # remove everything from the environment
58
59
60# Vectors and matrices in \R
61
62## Creating vectors using the ``combine" operator
63
64a <- c(1,3,5) # create a vector by combining values
65a
66a[2] # select the second element
67b <- c("one","three","five") # also works with strings
68b
69b[2]
70a <- c(a,a) # can apply recursively
71a
72a <- c(a,b) # mixing types---what happens?
73a # all converted to the same type
74
75## Sequences and replication
76
77a <- seq(from=1,to=5,by=1) # from 1 to 5 the slow way
78b <- 1:5 # a shortcut!
79a==b # all TRUE
80rep(1,times=5) # a lot of ones
81rep(1:5,times=2) # repeat an entire sequence
82rep(1:5,each=2) # same, but element-wise
83rep(1:5,times=5:1) # can vary the count of each element
84
85## Any, all, and which (with vectors)
86
87a <- 1:5 # create a vector
88a>2 # some TRUE, some FALSE
89any(a>2) # are any elements TRUE?
90all(a>2) # are all elements TRUE?
91which(a>2) # which indices are TRUE?
92
93## From vectors to matrices
94
95a <- matrix(data=1:25, nrow=5, ncol=5) # create a matrix the "formal" way
96a
97a[1,2] # select a matrix element (two dimensions)
98a[1,] # just the first row
99all(a[1,]==a[1,1:5]) # show the equivalence
100a[,2] # can also perform for columns
101a[2:3,3:5] # select submatrices
102a[-1,] # nice trick: negative numbers omit cells!
103a[-2,-2] # get rid of row two, column two
104
105
106b <- cbind(1:5,1:5) # another way to create matrices
107b
108d <- rbind(1:5,1:5) # can perform with rows, too
109d
110cbind(b,d) # no go: must have compatible dimensions!
111dim(b) # what were those dimensions, anyway?
112dim(d)
113NROW(b)
114NCOL(b)
115cbind(b,b) # combining two matrices
116
117
118t(b) # can transpose b
119cbind(t(b),d) # now it works
120rbind(t(b),d) # now it works
121
122
123# Element-wise operations
124
125## Most arithmetic operators are applied element-wise
126
127a <- 1:5
128a + 1 # addition
129a * 2 # multiplication
130a / 3 # division
131a - 4 # subtraction
132a ^ 5 # you get the idea...
133
134
135a + a # also works on pairs of vectors
136a * a
137a + 1:6 # problem: need same length
138
139a <- rbind(1:5,2:6) # same principles apply to matrices
140b <- rbind(3:7,4:8)
141a + b
142a / b
143a %*% t(b) # matrix multiplication
144
145## Logical operators (generally) work like arithmetic ones
146
147a > 0 # each value greater than zero?
148a == b # corresponding values equivalent?
149a != b # corresponding values not equivalent?
150!(a == b) # same as above
151(a>2) | (b>4) # the OR operator
152(a>2) & (b>4) # the AND operator
153(a>2) || (b>4) # beware the "double-pipe"!
154(a>2) && (b>4) # (and the "double-ampersand"!)
155
156## Ditto for many other basic transformations
157
158log(a)
159exp(b)
160sqrt(a+b) # note that we can nest statements!
161log((sqrt(a+b)+a)*b) # as recursive as we wanna be
162
163
164# Data frames
165
166## Build a data frame from scratch, containing three columns of data
167
168d <- data.frame(income=1:5,sane=c(T,T,T,T,F),name=LETTERS[1:5])
169d
170d[1,2] # acts a lot like a matrix!
171d[,1]*5
172d[-1,]
173d$sane # can use dollar sign notation
174d$sane[3]<-FALSE # making changes
175d
176d[2,3] # shows factors for string values
177d$name <- LETTERS[1:5] # eliminate evil factors by overwriting
178d[2,3]
179
180## If you want to do without factors
181
182d <- data.frame(income=1:5,sane=c(T,T,T,T,F),name=LETTERS[1:5],stringsAsFactors=FALSE)
183d
184d <- as.data.frame(cbind(1:5,2:6)) # can create from matrices
185d
186is.data.frame(d) # how can we tell it's not a matrix?
187is.matrix(d) # the truth comes out
188
189
190# Finding built-in data sets
191# Many packages have built-in data for testing and educational purposes
192
193data() # lists them all
194?USArrests # get help on a data set
195data(USArrests) # load the data set
196USArrests # view the object
197
198
199# Elementary visualization
200
201## R's workhorse is the ``plot" command
202
203plot(USArrests$Murder,USArrests$UrbanPop) # using dollar sign notation
204plot(USArrests$Murder,USArrests$UrbanPop,log="xy") # log-log scale
205
206## Adding plot title and axis labels
207
208plot(USArrests$Murder,USArrests$Assault,xlab="Murder",ylab="Assault",main="USArrests")
209
210## Can also add text
211
212plot(USArrests$Murder,USArrests$Assault,xlab="Murder",ylab="Assault", main="USArrests",type="n")
213text(USArrests$Murder,USArrests$Assault,rownames(USArrests))
214
215## Histograms and boxplots are often helpful
216
217hist(USArrests$Murder)
218boxplot(USArrests)
219
220
221
222########################################################################################
223# network Objects: Importing, Exploring, and Manipulating Network Data
224
225# Built-in Network Datasets
226
227library(network) # Make sure that network package is loaded
228data(package="network") # List available datasets in network package
229data(flo) # Load a built-in data set; see ?flo for more
230flo # Examine the flo adjacency matrix
231
232## For more information
233
234?data
235?flo
236
237
238# Importing Relational Data
239
240# Be sure to be in the directory where you stored the data for the workshop.
241# If you've not set the working directory, you must do so now.
242# See above for how to do this.
243
244list.files() # Check what's in the working directory
245
246## Read an adjacency matrix (\R stores it as a data frame by default)
247
248relations <- read.csv("relationalData.csv",header=FALSE,stringsAsFactors=FALSE)
249relations
250
251## Here's a case where matrix format is preferred
252
253relations <- as.matrix(relations) # convert to matrix format
254
255## Read in some vertex attribute data (okay to leave it as a data frame)
256
257nodeInfo <- read.csv("vertexAttributes.csv",header=TRUE,stringsAsFactors=FALSE)
258nodeInfo
259
260## Since our relational data has no row/column names, let's set them now
261
262rownames(relations) <- nodeInfo$name
263colnames(relations) <- nodeInfo$name
264relations
265
266## For more information
267
268?list.files
269?read.csv
270?as.matrix
271?rownames
272?colnames
273
274
275# Creating network Objects
276
277nrelations<-network(relations,directed=FALSE) # Create a network object based on flo
278nrelations # Get a quick description of the data
279nempty <- network.initialize(5) # Create an empty graph with 5 vertices
280nempty # Compare with nrelations
281
282## For more information
283
284?network
285?as.network.matrix
286
287
288# Description and Visualization
289
290summary(nrelations) # Get an overall summary
291network.dyadcount(nrelations) # How many dyads in nflo?
292network.edgecount(nrelations) # How many edges are present?
293network.size(nrelations) # How large is the network?
294as.sociomatrix(nrelations) # Show it as a sociomatrix
295nrelations[,] # Another way to do it
296
297
298plot(nrelations,displaylabels=T) # Plot with names
299plot(nrelations,displaylabels=T,mode="circle") # A less useful layout...
300library(sna) # Load the sna library
301gplot(nrelations) # Requires sna
302
303## For more information
304
305?summary.network
306?network.dyadcount
307?network.edgecount
308?as.sociomatrix
309?as.matrix.network
310?is.directed
311?plot.network
312
313
314# Working With Edges
315
316## Adding edges
317
318g <- network.initialize(5) # Create an empty graph
319g[1,2] <- 1 # Add an edge from 1 to 2
320g[2,] <- 1 # Add edges from 2 to everyone else
321g # Examine the result
322m <- matrix(0, nrow=5, ncol=5) # Create an adjacency matrix
323m[3,4:5] <- 1 # Add entries from 3 to 4 and 5
324g[m>0] <- 1 # Add more entries
325g
326
327## Delete edges
328
329g[3,5] <- 0 # Remove the edge from 3 to 5
330g # It's gone!
331g[,] <- 0 # Remove all edges
332g # Now, an empty graph
333
334## Testing adjacency
335
336nrelations[4,5] # Emma to Sarah?
337nrelations[4,] # Entire Emma row
338nrelations[1:4,5:8] # Subsets are possible
339nrelations[-9,-9] # Leaving Gil out of the picture
340
341## Setting edge values
342
343m<-nrelations[,] # copy over the relational structure
344m[upper.tri(m)>0]<-rep(1:3,times=3) # give different edge values
345m<-symmetrize(m,rule="upper")
346m
347nrelations %e% "strength" <- m # Add the valued ties back in
348
349## Retrieving edge values
350
351list.edge.attributes(nrelations) # See what???s available
352nrelations %e% "strength" # Use the %e% operator
353as.sociomatrix(nrelations,attrname="strength") # Can also do it this way
354nrelations[,] # Original tie structure is preserved
355
356## For more information
357
358?network.extraction
359?add.edge
360?delete.edges
361?delete.vertices
362?get.edges
363?upper.tri
364?symmetrize
365
366
367# Network and Vertex Attributes
368
369## Add some attributes
370
371nrelations %v% "id" <- nodeInfo$id # Add in our vertex attributes
372nrelations %v% "age" <- nodeInfo$age
373nrelations %v% "sex" <- nodeInfo$sex
374nrelations %v% "handed" <- nodeInfo$handed
375nrelations %v% "lastDocVisit" <- nodeInfo$lastDocVisit
376
377## Listing attributes
378
379list.vertex.attributes(nrelations) # List all vertex attributes
380list.network.attributes(nrelations) # List all network attributes
381
382## Retrieving attributes
383
384nrelations %v% "age" # Retrieve vertex ages
385nrelations %v% "id" # Retrieve vertex ids
386
387## For more information
388
389?attribute.methods
390
391
392
393########################################################################################
394# Classical Network Analysis with the sna Package for R
395
396# Getting started
397
398library(sna) # Load the sna library
399library(help="sna") # See a list of help pages for sna package
400load("introToSNAinR_sunbelt_2012.Rdata") # Load supplemental workshop data
401
402## For more information
403
404?help.start
405?library
406?sna
407
408
409# Network data in \sna
410
411## The \sna package can handle network data in many forms.
412## For instance, the function gden() calculates network density; we can use it on a network object,
413## an adjacency matrix, a list of such matrices, etc.
414
415data(flo)
416flo # Adjacency form
417gden(flo)
418nflo<-network(flo,directed=FALSE) # Network form
419gden(nflo)
420
421## The \sna package also supports a special kind of matrix called an ``sna edgelist."
422## These are three-column matrices, each row of which represents an edge
423## (via its sender, recipient, and value, respectively).
424## These ``sna edgelists" have special attributes that indicate their size, vertex names (if any),
425## and bipartite status (if applicable).
426
427eflo<-as.edgelist.sna(flo) # Coerce flo to an sna edgelist
428eflo
429attr(eflo,"n") # How many vertices are there?
430attr(eflo,"vnames") # Are there vertex names?
431as.sociomatrix.sna(eflo) # Can transform back w/ as.sociomatrix.sna
432
433## ``sna edgelists" can be handy with large data sets (as a simple alternative to network objects).
434## To make one, just add an ``$n$" attribute to a valid three-column matrix!
435
436mat<-cbind(rep(2,4),3:6,rep(1,4)) # Create edges from 2 to 3:6
437attr(mat,"n")<-6 # Set total number of vertices to 6
438mat
439gden(mat) # Can now pass to sna routines
440as.sociomatrix.sna(mat) # Can see in adjacency form
441
442## For more information
443
444?as.edgelist.sna
445?as.sociomatrix.sna
446?attr
447?sna
448
449
450# Network visualization with gplot()
451
452## Plotting the data we imported earlier
453
454gplot(nrelations)
455gplot(nrelations,displaylabels=TRUE) # This time with labels
456
457## Let's color the nodes in sex-stereotypic colors
458
459nodeColors<-ifelse(nodeInfo$sex=="F","hotpink","dodgerblue")
460gplot(relations,gmode="graph",displaylabels=TRUE,vertex.col=nodeColors)
461
462## Using data we just loaded in, plot the contiguity among nations in 1993
463## (from the Correlates of War (CoW) project See http://www.correlatesofwar.org/ for more information.)
464
465gplot(contig_1993) # The default visualization
466gplot(contig_1993, usearrows=FALSE) # Turn off arrows manually
467gplot(contig_1993, gmode="graph") # Can also tell gplot the data is undirected
468
469## We can add labels to the vertices
470
471gplot(contig_1993, gmode="graph",displaylabels=TRUE,label.cex=0.5,label.col="blue")
472
473## Other ways to specify the labeling
474
475gplot(contig_1993, gmode="graph",label=colnames(contig_1993[,]),label.cex=0.5,label.col="blue")
476
477## or
478
479gplot(contig_1993, gmode="graph",label=network.vertex.names(contig_1993),label.cex=0.5,label.col="blue")
480
481## Here's an example of directed data---militarized interstate disputes (MIDs) for 1993
482
483gplot(mids_1993,label.cex=0.5,label.col="blue",displaylabels=TRUE)
484
485## All those isolates can get in the way---we can suppress them using displayisolates
486
487gplot(mids_1993,label.cex=0.5,label.col="blue",displaylabels=TRUE,displayisolates=FALSE)
488
489## The default layout algorithm is that of Frutchterman-Reingold (1991), can use others
490
491gplot(mids_1993,label.cex=0.5,label.col="blue",
492displaylabels=TRUE,displayisolates=FALSE,mode="circle") # The infamous circle
493
494## or perhaps
495
496gplot(mids_1993,label.cex=0.5,label.col="blue",
497displaylabels=TRUE,displayisolates=FALSE,mode="mds") # MDS of position similarity
498
499## When a layout is generated, the results can be saved for later reuse:
500
501coords <- gplot(contig_1993) # Capture the magic of the moment
502coords # Show the vertex coordinates
503
504## Saved (or a priori) layouts can be used via the coord argument
505
506gplot(contig_1993, gmode="graph",label=colnames(contig_1993[,]),label.cex=0.5,label.col="blue",coord=coords)
507
508## When the default settings are insufficient, interactive mode allows for tweaking
509
510coords <- gplot(contig_1993, interactive=TRUE) # Modify and save
511gplot(contig_1993,coord=coords,displaylabels=TRUE,
512gmode="graph",label.cex=0.5,label.col="blue") # Should reproduce the modified layout
513
514## For more information
515
516?gplot
517?gplot.layout
518
519
520# Three-dimensional visualization with gplot3d() (requires the rgl package)
521
522gplot3d(contig_1993,displaylabels=TRUE) # Experience the future!
523
524## Other layouts are possible here, too:
525
526gplot3d(contig_1993,displaylabels=TRUE,mode="kamadakawai")
527
528## For more information
529?gplot3d
530?gplot3d.layout
531
532
533
534# Basic centrality indices: degree, betweenness, and closeness
535
536## We begin with the simplest case: degree centrality
537
538degree(mids_1993) # Default: total degree
539ideg <- degree(mids_1993, cmode="indegree") # Indegree for MIDs
540odeg <- degree(mids_1993, cmode="outdegree") # Outdegree for MIDs
541all(degree(mids_1993) == ideg+odeg) # In + out = total?
542
543## Once centrality scores are computed, we can handle them using standard \R methods:
544
545plot(ideg, odeg, type="n", xlab="Incoming MIDs", ylab="Outgoing MIDs") # Plot ideg by odeg
546abline(0, 1, lty=3)
547text(jitter(ideg), jitter(odeg), network.vertex.names(contig_1993), cex=0.75, col=2)
548
549## Plot simple histograms of the degree distribution:
550
551par(mfrow=c(2,2)) # Set up a 2x2 display
552hist(ideg, xlab="Indegree", main="Indegree Distribution", prob=TRUE)
553hist(odeg, xlab="Outdegree", main="Outdegree Distribution", prob=TRUE)
554hist(ideg+odeg, xlab="Total Degree", main="Total Degree Distribution", prob=TRUE)
555par(mfrow=c(1,1)) # Restore display
556
557## Centrality scores can also be used with other \sna routines, e.g., gplot()
558
559gplot(mids_1993, vertex.cex=(ideg+odeg)^0.5/2, vertex.sides=50,label.cex=0.4,
560        vertex.col=rgb(odeg/max(odeg),0,ideg/max(ideg)),displaylabels=TRUE)
561
562## Betweenness and closeness are also popular measures
563
564bet <- betweenness(contig_1993, gmode="graph") # Geographic betweenness
565bet
566gplot(contig_1993, vertex.cex=sqrt(bet)/25, gmode="graph") # Use w/gplot
567clo <- closeness(contig_1993) # Geographic closeness
568clo # A large world after all?
569
570## Can use \sna routines to explore alternatives to the common measures
571
572closeness2 <- function(x){ # Create an alternate closeness function!
573geo <- 1/geodist(x)$gdist # Get the matrix of 1/geodesic distance
574diag(geo) <- 0 # Define self-ties as 0
575apply(geo, 1, sum) # Return sum(1/geodist) for each vertex
576}
577
578clo2 <- closeness2(contig_1993) # Use our new function on contiguity data
579hist(clo2, xlab="Alt. Closeness", prob=TRUE) # Much better behaved!
580cor(clo2, bet) # Correlate with betweenness
581plot(clo2, bet) # Plot the bivariate relationship
582all(clo2/185==closeness(contig_1993,cmode="suminvundir")) # Actually, we support this in sna!
583
584## For more information
585
586?betweenness
587?bonpow
588?closeness
589?degree
590?evcent
591?graphcent
592?infocent
593?prestige
594?stresscent
595
596
597# From centrality to centralization
598
599centralization(mids_1993, degree, cmode="indegree") # Do MIDs concentrate?
600centralization(contig_1993, evcent) # Eigenvector centralization
601
602## For more information
603
604?centralization
605
606
607# Elementary graph-level indices
608
609gden(mids_1993) # Density
610grecip(mids_1993) # Dyadic reciprocity
611grecip(mids_1993, measure="edgewise") # Edgewise reciprocity
612gtrans(mids_1993) # Transitivity
613
614## For more information
615
616?gden
617?grecip
618?gtrans
619
620
621# Subgraph census routines
622
623dyad.census(mids_1993) # M,A,N counts
624dyad.census(contig_1993) # No As in undirected graphs
625triad.census(mids_1993) # Directed triad census
626triad.census(contig_1993, mode="graph") # Undirected triad census
627kpath.census(mids_1993, maxlen=6, tabulate.by.vertex=FALSE) # Count paths of length <=6
628kcycle.census(mids_1993, maxlen=6, tabulate.by.vertex=FALSE) # Count cycles of length <=6
629clique.census(mids_1993, tabulate.by.vertex=FALSE, enumerate=FALSE) # Find maximal cliques
630
631## Can also look at more detailed tabulation/co-membership for paths/cycles/cliques
632
633kpath.census(mids_1993, maxlen=4) # Show tabulation by vertex
634
635indirect <- kpath.census(mids_1993, maxlen=6, dyadic.tabulation="sum")$paths.bydyad
636gplot(indirect,label.cex=0.4,vertex.cex=0.75,
637displaylabels=TRUE,edge.col=rgb(0,0,0,0.25)) # Plot indirect MIDs
638
639## Component information can be obtained in various ways
640
641components(mids_1993) # Strong component count
642components(mids_1993, connected="weak") # Weak component count
643cd <- component.dist(mids_1993, connected="weak") # Get weak components
644cd
645
646## Component sizes
647
648plot(1:length(cd$cdist),cd$cdist,xlab="Size",ylab="Frequency")
649
650## Who's in the largest component?
651
652cl <- component.largest(mids_1993, connected="weak")
653cl
654
655## Plot the largest weak component
656
657gplot(mids_1993[cl,cl], boxed.lab=FALSE, label.cex=0.5, label.col=4,label=network.vertex.names(mids_1993)[cl])
658
659## Likewise, many routines exist for handling isolates
660
661is.isolate(mids_1993, 3) # Is the third vertex (BHM) an isolate?
662is.isolate(mids_1993, "BHM") # The peaceful islands?
663is.isolate(mids_1993, "USA") # Perhaps less so....
664isol <- isolates(mids_1993) # Get the entire list of isolates
665isol
666
667network.vertex.names(mids_1993)[isol] # Which countries are isolates?
668
669## Another way to remove isolates from sociograms
670
671gplot(mids_1993[-isol,-isol], label.cex=0.5, label.col=4,label=network.vertex.names(mids_1993)[-isol])
672
673## For more information
674
675?clique.census
676?components
677?component.dist
678?dyad.census
679?is.isolate
680?isolates
681?kcycle.census
682?kpath.census
683?triad.census
684
685
686# Elementary random graph generation
687
688rgraph(10) # A uniform random digraph of order 10
689rgraph(10, tprob=3/9) # Homogeneous Bernoulli w/mean degree 3
690rgraph(10, tprob=3/9, mode="graph") # As above, but undirected
691rgnm(1, 10, 20) # Uniform conditional on order, edges
692rguman(1, 10, mut=0.5, asym=0.25, null=0.25) # Homogeneous multinomial on dyad census
693rguman(1, 10, mut=0, asym=1, null=0) # An extreme case: random tournament
694gplot3d(rgws(1,50,1,2,0)) # A Watts-Strogatz process - baseline
695gplot3d(rgws(1,50,1,2,0.05)) # ...with rewiring probability 0.05
696gplot3d(rgws(1,50,1,2,0.2)) # ...with rewiring probability 0.2
697
698## For more information
699
700?rgbn
701?rgmn
702?rgraph
703?rguman
704?rgws
705
706
707# Basic connectivity/distance measurement, and cohesion
708
709g <- rgraph(20, tprob=3/19) # Start with a random digraph
710g
711is.connected(g) # Is g strongly connected?
712is.connected(g, connected="weak") # How about weakly connected?
713geodist(g) # Get information on geodesics
714reachability(g) # Return the reachability matrix
715symmetrize(g) # Symmetrize g using the "or" rule
716symmetrize(g, rule="strong") # Symmetrize g using the "and" rule
717
718## Several ways to get relatively cohesive groups
719
720clique.census(g) # Maximal clique census
721kcores(g) # k-cores (by degree)
722cutpoints(g) # find articulation points
723
724## Showing cohesion information can aid visualization. Here, show critical positions
725
726gplot(contig_1993,vertex.col=2+cutpoints(contig_1993,mode="graph",return.indicator=T))
727
728## Show the nesting of cores
729
730kc<-kcores(contig_1993,cmode="indegree")
731gplot(contig_1993,vertex.col=rainbow(max(kc)+1)[kc+1])
732
733## Showing members of the 5-core only
734
735gplot(contig_1993[kc>4,kc>4],vertex.col=rainbow(max(kc)+1)[kc[kc>4]+1])
736
737## For more information
738
739?bicomponent.dist
740?cutpoints
741?geodist
742?kcores
743?is.connected
744?reachability
745?symmetrize
746
747
748# Positional analysis
749
750## Generate a structural equivalence clustering of the CoW alliance data
751
752gplot(alliances_1993,gmode="graph",vertex.cex=0.5) #An initial look....
753ec <- equiv.clust(alliances_1993, mode="graph",plabels=network.vertex.names(alliances_1993))
754ec # The clustering
755plot(ec) # Plot the dendrogram
756rect.hclust(ec$cluster, h=20)
757
758## Use the clustering to form an structural equivalence (SE) blockmodel
759
760bm <- blockmodel(alliances_1993, ec, h=20)
761bm
762
763## We can display the blockmodel in several ways
764
765plot.sociomatrix(alliances_1993[bm$order.vector,bm$order.vector],drawlab=FALSE)
766
767## Extract the block image
768
769bimage <- bm$block.model
770bimage
771bimage[is.nan(bimage)] <- 1
772
773## Visualizing the block image (with self-reflexive ties)
774
775gplot(bimage, diag=TRUE, edge.lwd=bimage*5,
776      vertex.cex=sqrt(table(bm$block.membership))/2,
777      gmode="graph", vertex.sides=50,
778      vertex.col=gray(1-diag(bimage)))
779
780
781
782