StatnetPapers: tutorial.2.r

File tutorial.2.r, 3.4 KB (added by lxwang, 5 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
6library("statnet")
7data("faux.magnolia.high")
8fmh <- faux.magnolia.high
9fmh
10plot(fmh, displayisolates = FALSE)
11table(component.dist(fmh)$csize)
12summary(fmh)
13plot(fmh, displayisolates = FALSE, vertex.col = "Grade")
14fmh.degreedist <- table(degree(fmh, cmode = "indegree"))
15fmh.degreedist
16summary(fmh ~ degree(0:8))
17help(package = "sna")
18summary(fmh ~ degree(0:8, "Sex"))
19summary(fmh ~ triangle)
20summary(fmh ~ edges + triangle)
21mixingmatrix(fmh, "Grade")
22gr <- fmh %v% "Grade"
23table(gr)
24save.image()
25
26model1 <- ergm(fmh ~ edges)
27summary(model1)
28names(model1)
29model1$coef
30model1$mle.lik
31model2 <- ergm(fmh ~ edges + nodematch("Grade") + nodematch("Race") +
32   nodematch("Sex"))
33summary(model2)
34model2$mle.lik
35sim2 <- simulate(model2, burnin = 1e+6, verbose = TRUE, seed = 9)
36mixingmatrix(sim2, "Race")
37mixingmatrix(fmh, "Race")
38plot(summary(fmh ~ degree(0:10)), type = "l", lty = 1, lwd = 2,
39   xlab = "Degree", ylab = "Count")
40lines(summary(sim2 ~ degree(0:10)), lty = 2, lwd = 3)
41legend("topright", legend = c("Observed", "Simulated"), lwd = 3,
42   lty = 1:2)
43c(fmh = summary(fmh ~ triangle), sim2 = summary(sim2 ~ triangle))
44
45model3 <- ergm(fmh ~ edges + triangle, seed = 99)
46pdf("model3diagnostics.pdf")
47mcmc.diagnostics(model3)
48dev.off()
49model3 <- ergm(fmh ~ edges + triangle, verbose = TRUE, seed = 99)
50model3.take2 <- ergm(fmh ~ edges + triangle, MCMCsamplesize = 1e+5,
51   interval = 1000, verbose = TRUE, seed = 88)
52model3.take3 <- ergm(fmh ~ edges + triangle, maxit = 25, seed = 888,
53   control = control.ergm(steplength = 0.25), verbose = TRUE)
54
55model4.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)
59model4.take1$coef
60model4.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)
64model4.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)
68c(model4.take1$mle.lik, model4.take2$mle.lik, model4.take3$mle.lik)
69model4 <- model4.take3
70model4$coef
71
72sim4 <- simulate(model4, burnin = 1e+5, interval = 1e+5,
73     nsim = 100, verbose = TRUE, seed = 321)
74class(sim4)
75names(sim4)
76class(sim4$networks[[1]])
77model4.tridist <- sapply(sim4$networks, function(x) summary(x ~ triangle))
78
79fmh.tri <- summary(fmh ~ triangle)
80r <- range(c(model4.tridist, fmh.tri*c(.98, 1.02)))
81hist(model4.tridist, xlab = "Triangles", xlim=r)
82arrows(fmh.tri, 20, fmh.tri, 5, col = "red", lwd = 3)
83
84sum(fmh.tri <= model4.tridist)
85gof4.deg <- gof(model4 ~ degree, verbose = TRUE, burnin = 1e+5,
86     interval = 1e+5, seed = 246)
87plot(gof4.deg)
88gof4.deg
89gof4.esp.dist <- gof(model4 ~ espartners + distance, verbose = TRUE,
90     burnin = 1e+5, interval = 1e+5, seed = 642)
91getOption("device")(width=8,height=4)
92par(mfrow = c(1, 2))
93plot(gof4.esp.dist, plotlogodds = TRUE)
94
95