# StatnetPapers: tutorial.2.r

File tutorial.2.r, 3.4 KB (added by lxwang, 6 years ago) |
---|

Line | |
---|---|

1 | # This code is designed to run under R2.8.1 and ergm2.2 |

2 | |

3 | #install.packages("statnet") #uncomment this line if statnet not yet installed |

4 | #setwd("Your Directory Here") #modify, uncomment and use if desired |

5 | |

6 | library("statnet") |

7 | data("faux.magnolia.high") |

8 | fmh <- faux.magnolia.high |

9 | fmh |

10 | plot(fmh, displayisolates = FALSE) |

11 | table(component.dist(fmh)$csize) |

12 | summary(fmh) |

13 | plot(fmh, displayisolates = FALSE, vertex.col = "Grade") |

14 | fmh.degreedist <- table(degree(fmh, cmode = "indegree")) |

15 | fmh.degreedist |

16 | summary(fmh ~ degree(0:8)) |

17 | help(package = "sna") |

18 | summary(fmh ~ degree(0:8, "Sex")) |

19 | summary(fmh ~ triangle) |

20 | summary(fmh ~ edges + triangle) |

21 | mixingmatrix(fmh, "Grade") |

22 | gr <- fmh %v% "Grade" |

23 | table(gr) |

24 | save.image() |

25 | |

26 | model1 <- ergm(fmh ~ edges) |

27 | summary(model1) |

28 | names(model1) |

29 | model1$coef |

30 | model1$mle.lik |

31 | model2 <- ergm(fmh ~ edges + nodematch("Grade") + nodematch("Race") + |

32 | nodematch("Sex")) |

33 | summary(model2) |

34 | model2$mle.lik |

35 | sim2 <- simulate(model2, burnin = 1e+6, verbose = TRUE, seed = 9) |

36 | mixingmatrix(sim2, "Race") |

37 | mixingmatrix(fmh, "Race") |

38 | plot(summary(fmh ~ degree(0:10)), type = "l", lty = 1, lwd = 2, |

39 | xlab = "Degree", ylab = "Count") |

40 | lines(summary(sim2 ~ degree(0:10)), lty = 2, lwd = 3) |

41 | legend("topright", legend = c("Observed", "Simulated"), lwd = 3, |

42 | lty = 1:2) |

43 | c(fmh = summary(fmh ~ triangle), sim2 = summary(sim2 ~ triangle)) |

44 | |

45 | model3 <- ergm(fmh ~ edges + triangle, seed = 99) |

46 | pdf("model3diagnostics.pdf") |

47 | mcmc.diagnostics(model3) |

48 | dev.off() |

49 | model3 <- ergm(fmh ~ edges + triangle, verbose = TRUE, seed = 99) |

50 | model3.take2 <- ergm(fmh ~ edges + triangle, MCMCsamplesize = 1e+5, |

51 | interval = 1000, verbose = TRUE, seed = 88) |

52 | model3.take3 <- ergm(fmh ~ edges + triangle, maxit = 25, seed = 888, |

53 | control = control.ergm(steplength = 0.25), verbose = TRUE) |

54 | |

55 | model4.take1 <- ergm(fmh ~ edges + nodematch("Grade") + |

56 | nodematch("Race") + nodematch("Sex") + gwesp(0, fixed = TRUE), |

57 | MCMCsamplesize = 1e+5, maxit = 15, verbose = TRUE, |

58 | control = control.ergm(steplength = 0.25), seed = 123) |

59 | model4.take1$coef |

60 | model4.take2 <- ergm(fmh ~ edges + nodematch("Grade") + |

61 | nodematch("Race") + nodematch("Sex") + gwesp(0.1, fixed = TRUE), |

62 | MCMCsamplesize = 1e+5, maxit = 15, verbose = TRUE, |

63 | control = control.ergm(steplength = 0.25), seed = 123) |

64 | model4.take3 <- ergm(fmh ~ edges + nodematch("Grade") + |

65 | nodematch("Race") + nodematch("Sex") + gwesp(0.2, fixed = TRUE), |

66 | MCMCsamplesize = 1e+5, maxit = 15, verbose = TRUE, |

67 | control = control.ergm(steplength = 0.25), seed = 123) |

68 | c(model4.take1$mle.lik, model4.take2$mle.lik, model4.take3$mle.lik) |

69 | model4 <- model4.take3 |

70 | model4$coef |

71 | |

72 | sim4 <- simulate(model4, burnin = 1e+5, interval = 1e+5, |

73 | nsim = 100, verbose = TRUE, seed = 321) |

74 | class(sim4) |

75 | names(sim4) |

76 | class(sim4$networks[[1]]) |

77 | model4.tridist <- sapply(sim4$networks, function(x) summary(x ~ triangle)) |

78 | |

79 | fmh.tri <- summary(fmh ~ triangle) |

80 | r <- range(c(model4.tridist, fmh.tri*c(.98, 1.02))) |

81 | hist(model4.tridist, xlab = "Triangles", xlim=r) |

82 | arrows(fmh.tri, 20, fmh.tri, 5, col = "red", lwd = 3) |

83 | |

84 | sum(fmh.tri <= model4.tridist) |

85 | gof4.deg <- gof(model4 ~ degree, verbose = TRUE, burnin = 1e+5, |

86 | interval = 1e+5, seed = 246) |

87 | plot(gof4.deg) |

88 | gof4.deg |

89 | gof4.esp.dist <- gof(model4 ~ espartners + distance, verbose = TRUE, |

90 | burnin = 1e+5, interval = 1e+5, seed = 642) |

91 | getOption("device")(width=8,height=4) |

92 | par(mfrow = c(1, 2)) |

93 | plot(gof4.esp.dist, plotlogodds = TRUE) |

94 | |

95 |