# Resources: introToSNAinR_sunbelt_2012.R

File introToSNAinR_sunbelt_2012.R, 22.4 KB (added by morrism, 6 years ago) |
---|

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 | |

8 | install.packages("statnet") |

9 | |

10 | ## Now, and in the future, you can install/update \statnet at any point using: |

11 | |

12 | update.packages("statnet") |

13 | |

14 | ## Attach \statnet to your \R session by typing: |

15 | |

16 | library(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 | |

24 | install.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 | |

32 | setwd("full.path.for.the.folder") |

33 | |

34 | |

35 | ###################################################################################### |

36 | # Introduction to basic \R syntax |

37 | |

38 | 1 + 3 # evaluation |

39 | a <- 3 # assignment |

40 | a # evaluation |

41 | a<-3 # spacing does not matter |

42 | a <- 3 # spacing does not matter |

43 | sqrt(a) # use the square root function |

44 | b <- sqrt(a) # use function and save result |

45 | b |

46 | c # evaluate something that is not there |

47 | a == b # is a equivalent to b? |

48 | a != b # is a not equal to b? |

49 | ls() # list objects in the global environment |

50 | |

51 | help.start() # get help with R generally |

52 | ?sqrt # get specific help for a function |

53 | ??sqrt # looking for help pertaining to sqrt |

54 | apropos("sq") # it's on the tip of my tongue... |

55 | |

56 | rm(a) # remove a single object |

57 | rm(list=ls()) # remove everything from the environment |

58 | |

59 | |

60 | # Vectors and matrices in \R |

61 | |

62 | ## Creating vectors using the ``combine" operator |

63 | |

64 | a <- c(1,3,5) # create a vector by combining values |

65 | a |

66 | a[2] # select the second element |

67 | b <- c("one","three","five") # also works with strings |

68 | b |

69 | b[2] |

70 | a <- c(a,a) # can apply recursively |

71 | a |

72 | a <- c(a,b) # mixing types---what happens? |

73 | a # all converted to the same type |

74 | |

75 | ## Sequences and replication |

76 | |

77 | a <- seq(from=1,to=5,by=1) # from 1 to 5 the slow way |

78 | b <- 1:5 # a shortcut! |

79 | a==b # all TRUE |

80 | rep(1,times=5) # a lot of ones |

81 | rep(1:5,times=2) # repeat an entire sequence |

82 | rep(1:5,each=2) # same, but element-wise |

83 | rep(1:5,times=5:1) # can vary the count of each element |

84 | |

85 | ## Any, all, and which (with vectors) |

86 | |

87 | a <- 1:5 # create a vector |

88 | a>2 # some TRUE, some FALSE |

89 | any(a>2) # are any elements TRUE? |

90 | all(a>2) # are all elements TRUE? |

91 | which(a>2) # which indices are TRUE? |

92 | |

93 | ## From vectors to matrices |

94 | |

95 | a <- matrix(data=1:25, nrow=5, ncol=5) # create a matrix the "formal" way |

96 | a |

97 | a[1,2] # select a matrix element (two dimensions) |

98 | a[1,] # just the first row |

99 | all(a[1,]==a[1,1:5]) # show the equivalence |

100 | a[,2] # can also perform for columns |

101 | a[2:3,3:5] # select submatrices |

102 | a[-1,] # nice trick: negative numbers omit cells! |

103 | a[-2,-2] # get rid of row two, column two |

104 | |

105 | |

106 | b <- cbind(1:5,1:5) # another way to create matrices |

107 | b |

108 | d <- rbind(1:5,1:5) # can perform with rows, too |

109 | d |

110 | cbind(b,d) # no go: must have compatible dimensions! |

111 | dim(b) # what were those dimensions, anyway? |

112 | dim(d) |

113 | NROW(b) |

114 | NCOL(b) |

115 | cbind(b,b) # combining two matrices |

116 | |

117 | |

118 | t(b) # can transpose b |

119 | cbind(t(b),d) # now it works |

120 | rbind(t(b),d) # now it works |

121 | |

122 | |

123 | # Element-wise operations |

124 | |

125 | ## Most arithmetic operators are applied element-wise |

126 | |

127 | a <- 1:5 |

128 | a + 1 # addition |

129 | a * 2 # multiplication |

130 | a / 3 # division |

131 | a - 4 # subtraction |

132 | a ^ 5 # you get the idea... |

133 | |

134 | |

135 | a + a # also works on pairs of vectors |

136 | a * a |

137 | a + 1:6 # problem: need same length |

138 | |

139 | a <- rbind(1:5,2:6) # same principles apply to matrices |

140 | b <- rbind(3:7,4:8) |

141 | a + b |

142 | a / b |

143 | a %*% t(b) # matrix multiplication |

144 | |

145 | ## Logical operators (generally) work like arithmetic ones |

146 | |

147 | a > 0 # each value greater than zero? |

148 | a == b # corresponding values equivalent? |

149 | a != 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 | |

158 | log(a) |

159 | exp(b) |

160 | sqrt(a+b) # note that we can nest statements! |

161 | log((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 | |

168 | d <- data.frame(income=1:5,sane=c(T,T,T,T,F),name=LETTERS[1:5]) |

169 | d |

170 | d[1,2] # acts a lot like a matrix! |

171 | d[,1]*5 |

172 | d[-1,] |

173 | d$sane # can use dollar sign notation |

174 | d$sane[3]<-FALSE # making changes |

175 | d |

176 | d[2,3] # shows factors for string values |

177 | d$name <- LETTERS[1:5] # eliminate evil factors by overwriting |

178 | d[2,3] |

179 | |

180 | ## If you want to do without factors |

181 | |

182 | d <- data.frame(income=1:5,sane=c(T,T,T,T,F),name=LETTERS[1:5],stringsAsFactors=FALSE) |

183 | d |

184 | d <- as.data.frame(cbind(1:5,2:6)) # can create from matrices |

185 | d |

186 | is.data.frame(d) # how can we tell it's not a matrix? |

187 | is.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 | |

193 | data() # lists them all |

194 | ?USArrests # get help on a data set |

195 | data(USArrests) # load the data set |

196 | USArrests # view the object |

197 | |

198 | |

199 | # Elementary visualization |

200 | |

201 | ## R's workhorse is the ``plot" command |

202 | |

203 | plot(USArrests$Murder,USArrests$UrbanPop) # using dollar sign notation |

204 | plot(USArrests$Murder,USArrests$UrbanPop,log="xy") # log-log scale |

205 | |

206 | ## Adding plot title and axis labels |

207 | |

208 | plot(USArrests$Murder,USArrests$Assault,xlab="Murder",ylab="Assault",main="USArrests") |

209 | |

210 | ## Can also add text |

211 | |

212 | plot(USArrests$Murder,USArrests$Assault,xlab="Murder",ylab="Assault", main="USArrests",type="n") |

213 | text(USArrests$Murder,USArrests$Assault,rownames(USArrests)) |

214 | |

215 | ## Histograms and boxplots are often helpful |

216 | |

217 | hist(USArrests$Murder) |

218 | boxplot(USArrests) |

219 | |

220 | |

221 | |

222 | ######################################################################################## |

223 | # network Objects: Importing, Exploring, and Manipulating Network Data |

224 | |

225 | # Built-in Network Datasets |

226 | |

227 | library(network) # Make sure that network package is loaded |

228 | data(package="network") # List available datasets in network package |

229 | data(flo) # Load a built-in data set; see ?flo for more |

230 | flo # 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 | |

244 | list.files() # Check what's in the working directory |

245 | |

246 | ## Read an adjacency matrix (\R stores it as a data frame by default) |

247 | |

248 | relations <- read.csv("relationalData.csv",header=FALSE,stringsAsFactors=FALSE) |

249 | relations |

250 | |

251 | ## Here's a case where matrix format is preferred |

252 | |

253 | relations <- as.matrix(relations) # convert to matrix format |

254 | |

255 | ## Read in some vertex attribute data (okay to leave it as a data frame) |

256 | |

257 | nodeInfo <- read.csv("vertexAttributes.csv",header=TRUE,stringsAsFactors=FALSE) |

258 | nodeInfo |

259 | |

260 | ## Since our relational data has no row/column names, let's set them now |

261 | |

262 | rownames(relations) <- nodeInfo$name |

263 | colnames(relations) <- nodeInfo$name |

264 | relations |

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 | |

277 | nrelations<-network(relations,directed=FALSE) # Create a network object based on flo |

278 | nrelations # Get a quick description of the data |

279 | nempty <- network.initialize(5) # Create an empty graph with 5 vertices |

280 | nempty # Compare with nrelations |

281 | |

282 | ## For more information |

283 | |

284 | ?network |

285 | ?as.network.matrix |

286 | |

287 | |

288 | # Description and Visualization |

289 | |

290 | summary(nrelations) # Get an overall summary |

291 | network.dyadcount(nrelations) # How many dyads in nflo? |

292 | network.edgecount(nrelations) # How many edges are present? |

293 | network.size(nrelations) # How large is the network? |

294 | as.sociomatrix(nrelations) # Show it as a sociomatrix |

295 | nrelations[,] # Another way to do it |

296 | |

297 | |

298 | plot(nrelations,displaylabels=T) # Plot with names |

299 | plot(nrelations,displaylabels=T,mode="circle") # A less useful layout... |

300 | library(sna) # Load the sna library |

301 | gplot(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 | |

318 | g <- network.initialize(5) # Create an empty graph |

319 | g[1,2] <- 1 # Add an edge from 1 to 2 |

320 | g[2,] <- 1 # Add edges from 2 to everyone else |

321 | g # Examine the result |

322 | m <- matrix(0, nrow=5, ncol=5) # Create an adjacency matrix |

323 | m[3,4:5] <- 1 # Add entries from 3 to 4 and 5 |

324 | g[m>0] <- 1 # Add more entries |

325 | g |

326 | |

327 | ## Delete edges |

328 | |

329 | g[3,5] <- 0 # Remove the edge from 3 to 5 |

330 | g # It's gone! |

331 | g[,] <- 0 # Remove all edges |

332 | g # Now, an empty graph |

333 | |

334 | ## Testing adjacency |

335 | |

336 | nrelations[4,5] # Emma to Sarah? |

337 | nrelations[4,] # Entire Emma row |

338 | nrelations[1:4,5:8] # Subsets are possible |

339 | nrelations[-9,-9] # Leaving Gil out of the picture |

340 | |

341 | ## Setting edge values |

342 | |

343 | m<-nrelations[,] # copy over the relational structure |

344 | m[upper.tri(m)>0]<-rep(1:3,times=3) # give different edge values |

345 | m<-symmetrize(m,rule="upper") |

346 | m |

347 | nrelations %e% "strength" <- m # Add the valued ties back in |

348 | |

349 | ## Retrieving edge values |

350 | |

351 | list.edge.attributes(nrelations) # See what???s available |

352 | nrelations %e% "strength" # Use the %e% operator |

353 | as.sociomatrix(nrelations,attrname="strength") # Can also do it this way |

354 | nrelations[,] # 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 | |

371 | nrelations %v% "id" <- nodeInfo$id # Add in our vertex attributes |

372 | nrelations %v% "age" <- nodeInfo$age |

373 | nrelations %v% "sex" <- nodeInfo$sex |

374 | nrelations %v% "handed" <- nodeInfo$handed |

375 | nrelations %v% "lastDocVisit" <- nodeInfo$lastDocVisit |

376 | |

377 | ## Listing attributes |

378 | |

379 | list.vertex.attributes(nrelations) # List all vertex attributes |

380 | list.network.attributes(nrelations) # List all network attributes |

381 | |

382 | ## Retrieving attributes |

383 | |

384 | nrelations %v% "age" # Retrieve vertex ages |

385 | nrelations %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 | |

398 | library(sna) # Load the sna library |

399 | library(help="sna") # See a list of help pages for sna package |

400 | load("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 | |

415 | data(flo) |

416 | flo # Adjacency form |

417 | gden(flo) |

418 | nflo<-network(flo,directed=FALSE) # Network form |

419 | gden(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 | |

427 | eflo<-as.edgelist.sna(flo) # Coerce flo to an sna edgelist |

428 | eflo |

429 | attr(eflo,"n") # How many vertices are there? |

430 | attr(eflo,"vnames") # Are there vertex names? |

431 | as.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 | |

436 | mat<-cbind(rep(2,4),3:6,rep(1,4)) # Create edges from 2 to 3:6 |

437 | attr(mat,"n")<-6 # Set total number of vertices to 6 |

438 | mat |

439 | gden(mat) # Can now pass to sna routines |

440 | as.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 | |

454 | gplot(nrelations) |

455 | gplot(nrelations,displaylabels=TRUE) # This time with labels |

456 | |

457 | ## Let's color the nodes in sex-stereotypic colors |

458 | |

459 | nodeColors<-ifelse(nodeInfo$sex=="F","hotpink","dodgerblue") |

460 | gplot(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 | |

465 | gplot(contig_1993) # The default visualization |

466 | gplot(contig_1993, usearrows=FALSE) # Turn off arrows manually |

467 | gplot(contig_1993, gmode="graph") # Can also tell gplot the data is undirected |

468 | |

469 | ## We can add labels to the vertices |

470 | |

471 | gplot(contig_1993, gmode="graph",displaylabels=TRUE,label.cex=0.5,label.col="blue") |

472 | |

473 | ## Other ways to specify the labeling |

474 | |

475 | gplot(contig_1993, gmode="graph",label=colnames(contig_1993[,]),label.cex=0.5,label.col="blue") |

476 | |

477 | ## or |

478 | |

479 | gplot(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 | |

483 | gplot(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 | |

487 | gplot(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 | |

491 | gplot(mids_1993,label.cex=0.5,label.col="blue", |

492 | displaylabels=TRUE,displayisolates=FALSE,mode="circle") # The infamous circle |

493 | |

494 | ## or perhaps |

495 | |

496 | gplot(mids_1993,label.cex=0.5,label.col="blue", |

497 | displaylabels=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 | |

501 | coords <- gplot(contig_1993) # Capture the magic of the moment |

502 | coords # Show the vertex coordinates |

503 | |

504 | ## Saved (or a priori) layouts can be used via the coord argument |

505 | |

506 | gplot(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 | |

510 | coords <- gplot(contig_1993, interactive=TRUE) # Modify and save |

511 | gplot(contig_1993,coord=coords,displaylabels=TRUE, |

512 | gmode="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 | |

522 | gplot3d(contig_1993,displaylabels=TRUE) # Experience the future! |

523 | |

524 | ## Other layouts are possible here, too: |

525 | |

526 | gplot3d(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 | |

538 | degree(mids_1993) # Default: total degree |

539 | ideg <- degree(mids_1993, cmode="indegree") # Indegree for MIDs |

540 | odeg <- degree(mids_1993, cmode="outdegree") # Outdegree for MIDs |

541 | all(degree(mids_1993) == ideg+odeg) # In + out = total? |

542 | |

543 | ## Once centrality scores are computed, we can handle them using standard \R methods: |

544 | |

545 | plot(ideg, odeg, type="n", xlab="Incoming MIDs", ylab="Outgoing MIDs") # Plot ideg by odeg |

546 | abline(0, 1, lty=3) |

547 | text(jitter(ideg), jitter(odeg), network.vertex.names(contig_1993), cex=0.75, col=2) |

548 | |

549 | ## Plot simple histograms of the degree distribution: |

550 | |

551 | par(mfrow=c(2,2)) # Set up a 2x2 display |

552 | hist(ideg, xlab="Indegree", main="Indegree Distribution", prob=TRUE) |

553 | hist(odeg, xlab="Outdegree", main="Outdegree Distribution", prob=TRUE) |

554 | hist(ideg+odeg, xlab="Total Degree", main="Total Degree Distribution", prob=TRUE) |

555 | par(mfrow=c(1,1)) # Restore display |

556 | |

557 | ## Centrality scores can also be used with other \sna routines, e.g., gplot() |

558 | |

559 | gplot(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 | |

564 | bet <- betweenness(contig_1993, gmode="graph") # Geographic betweenness |

565 | bet |

566 | gplot(contig_1993, vertex.cex=sqrt(bet)/25, gmode="graph") # Use w/gplot |

567 | clo <- closeness(contig_1993) # Geographic closeness |

568 | clo # A large world after all? |

569 | |

570 | ## Can use \sna routines to explore alternatives to the common measures |

571 | |

572 | closeness2 <- function(x){ # Create an alternate closeness function! |

573 | geo <- 1/geodist(x)$gdist # Get the matrix of 1/geodesic distance |

574 | diag(geo) <- 0 # Define self-ties as 0 |

575 | apply(geo, 1, sum) # Return sum(1/geodist) for each vertex |

576 | } |

577 | |

578 | clo2 <- closeness2(contig_1993) # Use our new function on contiguity data |

579 | hist(clo2, xlab="Alt. Closeness", prob=TRUE) # Much better behaved! |

580 | cor(clo2, bet) # Correlate with betweenness |

581 | plot(clo2, bet) # Plot the bivariate relationship |

582 | all(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 | |

599 | centralization(mids_1993, degree, cmode="indegree") # Do MIDs concentrate? |

600 | centralization(contig_1993, evcent) # Eigenvector centralization |

601 | |

602 | ## For more information |

603 | |

604 | ?centralization |

605 | |

606 | |

607 | # Elementary graph-level indices |

608 | |

609 | gden(mids_1993) # Density |

610 | grecip(mids_1993) # Dyadic reciprocity |

611 | grecip(mids_1993, measure="edgewise") # Edgewise reciprocity |

612 | gtrans(mids_1993) # Transitivity |

613 | |

614 | ## For more information |

615 | |

616 | ?gden |

617 | ?grecip |

618 | ?gtrans |

619 | |

620 | |

621 | # Subgraph census routines |

622 | |

623 | dyad.census(mids_1993) # M,A,N counts |

624 | dyad.census(contig_1993) # No As in undirected graphs |

625 | triad.census(mids_1993) # Directed triad census |

626 | triad.census(contig_1993, mode="graph") # Undirected triad census |

627 | kpath.census(mids_1993, maxlen=6, tabulate.by.vertex=FALSE) # Count paths of length <=6 |

628 | kcycle.census(mids_1993, maxlen=6, tabulate.by.vertex=FALSE) # Count cycles of length <=6 |

629 | clique.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 | |

633 | kpath.census(mids_1993, maxlen=4) # Show tabulation by vertex |

634 | |

635 | indirect <- kpath.census(mids_1993, maxlen=6, dyadic.tabulation="sum")$paths.bydyad |

636 | gplot(indirect,label.cex=0.4,vertex.cex=0.75, |

637 | displaylabels=TRUE,edge.col=rgb(0,0,0,0.25)) # Plot indirect MIDs |

638 | |

639 | ## Component information can be obtained in various ways |

640 | |

641 | components(mids_1993) # Strong component count |

642 | components(mids_1993, connected="weak") # Weak component count |

643 | cd <- component.dist(mids_1993, connected="weak") # Get weak components |

644 | cd |

645 | |

646 | ## Component sizes |

647 | |

648 | plot(1:length(cd$cdist),cd$cdist,xlab="Size",ylab="Frequency") |

649 | |

650 | ## Who's in the largest component? |

651 | |

652 | cl <- component.largest(mids_1993, connected="weak") |

653 | cl |

654 | |

655 | ## Plot the largest weak component |

656 | |

657 | gplot(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 | |

661 | is.isolate(mids_1993, 3) # Is the third vertex (BHM) an isolate? |

662 | is.isolate(mids_1993, "BHM") # The peaceful islands? |

663 | is.isolate(mids_1993, "USA") # Perhaps less so.... |

664 | isol <- isolates(mids_1993) # Get the entire list of isolates |

665 | isol |

666 | |

667 | network.vertex.names(mids_1993)[isol] # Which countries are isolates? |

668 | |

669 | ## Another way to remove isolates from sociograms |

670 | |

671 | gplot(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 | |

688 | rgraph(10) # A uniform random digraph of order 10 |

689 | rgraph(10, tprob=3/9) # Homogeneous Bernoulli w/mean degree 3 |

690 | rgraph(10, tprob=3/9, mode="graph") # As above, but undirected |

691 | rgnm(1, 10, 20) # Uniform conditional on order, edges |

692 | rguman(1, 10, mut=0.5, asym=0.25, null=0.25) # Homogeneous multinomial on dyad census |

693 | rguman(1, 10, mut=0, asym=1, null=0) # An extreme case: random tournament |

694 | gplot3d(rgws(1,50,1,2,0)) # A Watts-Strogatz process - baseline |

695 | gplot3d(rgws(1,50,1,2,0.05)) # ...with rewiring probability 0.05 |

696 | gplot3d(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 | |

709 | g <- rgraph(20, tprob=3/19) # Start with a random digraph |

710 | g |

711 | is.connected(g) # Is g strongly connected? |

712 | is.connected(g, connected="weak") # How about weakly connected? |

713 | geodist(g) # Get information on geodesics |

714 | reachability(g) # Return the reachability matrix |

715 | symmetrize(g) # Symmetrize g using the "or" rule |

716 | symmetrize(g, rule="strong") # Symmetrize g using the "and" rule |

717 | |

718 | ## Several ways to get relatively cohesive groups |

719 | |

720 | clique.census(g) # Maximal clique census |

721 | kcores(g) # k-cores (by degree) |

722 | cutpoints(g) # find articulation points |

723 | |

724 | ## Showing cohesion information can aid visualization. Here, show critical positions |

725 | |

726 | gplot(contig_1993,vertex.col=2+cutpoints(contig_1993,mode="graph",return.indicator=T)) |

727 | |

728 | ## Show the nesting of cores |

729 | |

730 | kc<-kcores(contig_1993,cmode="indegree") |

731 | gplot(contig_1993,vertex.col=rainbow(max(kc)+1)[kc+1]) |

732 | |

733 | ## Showing members of the 5-core only |

734 | |

735 | gplot(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 | |

752 | gplot(alliances_1993,gmode="graph",vertex.cex=0.5) #An initial look.... |

753 | ec <- equiv.clust(alliances_1993, mode="graph",plabels=network.vertex.names(alliances_1993)) |

754 | ec # The clustering |

755 | plot(ec) # Plot the dendrogram |

756 | rect.hclust(ec$cluster, h=20) |

757 | |

758 | ## Use the clustering to form an structural equivalence (SE) blockmodel |

759 | |

760 | bm <- blockmodel(alliances_1993, ec, h=20) |

761 | bm |

762 | |

763 | ## We can display the blockmodel in several ways |

764 | |

765 | plot.sociomatrix(alliances_1993[bm$order.vector,bm$order.vector],drawlab=FALSE) |

766 | |

767 | ## Extract the block image |

768 | |

769 | bimage <- bm$block.model |

770 | bimage |

771 | bimage[is.nan(bimage)] <- 1 |

772 | |

773 | ## Visualizing the block image (with self-reflexive ties) |

774 | |

775 | gplot(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 |