Resources: ergm_code.r

File ergm_code.r, 7.6 KB (added by lxwang, 5 years ago)
Line 
1# SECTION 0.  GETTING STARTED
2
3setwd('C:\\Users\\morrism\\Desktop\\Martina\\Presentations\\a_CONFERENCES\\Sunbelt\\2012\\Workshops')
4
5#install.packages('statnet')
6library(statnet)
7
8#install.packages('network')
9#install.packages('ergm')
10#install.packages('sna')
11#library(network)
12#library(ergm)
13#library(sna)
14#update.packages('name.of.package')
15#install.packages('coda')
16#library(coda)
17
18# SECTION 1.  STATISTICAL NETWORK MODELING; THE ERGM COMMAND AND ERGM OBJECT.
19
20
21data(package='ergm')                            # tells us the datasets in our packages
22data(florentine)                                        # loads flomarriage & flobusiness data
23flomarriage                                             # Let's look at the flomarriage data
24plot(flomarriage)                                       # Let's view the flomarriage network
25
26flomodel.01 <- ergm(flomarriage~edges)                  # fit model
27flomodel.01                                                             # look at the model
28
29summary(flomodel.01)                                            # look in more depth
30flomodel.02 <- ergm(flomarriage~edges+triangle)
31summary(flomodel.02)
32class(flomodel.02)                              # this has the class ergm
33names(flomodel.02)                              # let's look straight at the ERGM obj.
34flomodel.02$coef                                # the $ allows you to pull an element out from
35flomodel.02$formula                     # a list
36flomodel.02$mle.lik
37
38wealth <- flomarriage %v% 'wealth'              # the %v% extracts vertex attributes from a wealth                                              # network
39plot(flomarriage, vertex.cex=wealth/25) # network plot with vertex size
40flomodel.03 <- ergm(flomarriage~edges+nodecov('wealth'))
41summary(flomodel.03)
42
43data(samplk)                                    # Let's try a model or two on
44ls()                                                    # directed data: Sampson's Monks
45samplk3
46plot(samplk3)
47sampmodel.01 <- ergm(samplk3~edges+mutual)# Is there a statistically significant
48summary(sampmodel.01)                           # tendency for ties to be reciprocated
49# ("mutuality")?
50
51data(faux.mesa.high)                            # Let's try a larger network
52mesa <- faux.mesa.high
53plot(mesa)
54mesa
55plot(mesa, vertex.col='Grade')
56legend('bottomleft',fill=7:12,legend=paste('Grade',7:12),cex=0.75)
57
58fauxmodel.01 <- ergm(mesa ~edges + nodematch('Grade',diff=T) +
59        nodematch('Race',diff=T))
60summary(fauxmodel.01)
61table(mesa %v% "Race")                          # Frequencies of race
62mixingmatrix(mesa, "Race")
63
64missnet <- network.initialize(10,directed=F)
65missnet[1,2] <- missnet[2,7] <- missnet[3,6] <- 1
66missnet[4,6] <- missnet[4,9] <- NA
67missnet
68plot(missnet)
69ergm(missnet~edges)
70ergm(missnet~edges+degree(2))
71missnet[4,6] <- missnet[4,9] <- 0
72ergm(missnet~edges+degree(2))
73
74# SECTION 2.  MODEL TERMS AVAILABLE FOR ergm ESTIMATION and SIMULATION
75
76help('ergm-terms')
77
78# SECTION 3.  NETWORK SIMULATION: THE SIMULATE COMMAND AND NETWORK.LIST OBJECTS.
79
80flomodel.03.sim <- simulate(flomodel.03,nsim=10)
81class(flomodel.03.sim)
82length(flomodel.03.sim)
83flomodel.03.sim[[1]]                    # double brackets pull an element
84plot(flomodel.03.sim[[1]], label= flomodel.03.sim[[1]] %v% "vertex.names")
85
86# SECTION 4.  EXAMINING THE QUALITY OF MODEL FIT - GOF.
87
88flomodel.03.gof <- gof(flomodel.03~degree)
89flomodel.03.gof
90plot(flomodel.03.gof)
91mesamodel.02 <- ergm(mesa~edges)
92mesamodel.02.gof <- gof(mesamodel.02~distance,nsim=10)
93plot(mesamodel.02.gof)
94
95# SECTION 5.  DIAGNOSTICS: TROUBLESHOOTING AND CHECKING FOR  MODEL DEGENERACY
96
97fit <- ergm(flobusiness~edges+degree(1),
98        control=control.ergm(MCMC.interval=1, MCMC.burnin=1000, seed=1))
99mcmc.diagnostics(fit, center=F)
100fit <- ergm(flobusiness~edges+degree(1))
101mcmc.diagnostics(fit, center=F)
102data('faux.magnolia.high')
103magnolia <- faux.magnolia.high
104plot(magnolia, vertex.cex=.5)
105fit <- ergm(magnolia~edges+triangle,
106 control=control.ergm(seed=1))
107mcmc.diagnostics(fit, center=F)
108fit <- ergm(magnolia~edges+triangle,
109 control=control.ergm(seed=1),
110 verbose=T)
111fit <- ergm(magnolia~edges+triangle,seed=1,
112 control = control.ergm(seed=1, MCMC.samplesize=20000),
113 verbose=T)
114mcmc.diagnostics(fit, center=F)
115
116fit <- ergm(magnolia~edges+gwesp(0.5,fixed=T),
117 control =  control.ergm(seed=1))
118mcmc.diagnostics(fit)
119fit <- ergm(magnolia~edges+gwesp(0.5,fixed=T)+nodematch('Grade')+nodematch('Race')+
120        nodematch('Sex'),
121 control = control.ergm(seed=1),
122 verbose=T)
123pdf('diagnostics1.pdf') #Use the recording function if possible, or send to pdf
124mcmc.diagnostics(fit)
125dev.off()
126
127fit <- ergm(magnolia~edges+gwesp(0.25,fixed=T)+nodematch('Grade')+nodematch('Race')+
128        nodematch('Sex'),
129 control = control.ergm(seed=1))
130pdf('diagnostics2.pdf') #Ditto
131mcmc.diagnostics(fit)
132dev.off()
133
134args(ergm)
135fit <- ergm(magnolia~edges+gwesp(0.25,fixed=T)+nodematch('Grade')+nodematch('Race')+
136        nodematch('Sex'),
137 control = control.ergm(seed=1,MCMC.samplesize=50000,MCMC.interval=1000),
138 verbose=T)
139pdf('diagnostics3.pdf')
140mcmc.diagnostics(fit)
141dev.off()
142save.image()
143
144# SECTION 6.  WORKING WITH EGOCENTRICALLY SAMPLED NETWORK DATA
145
146ego.net <- network.initialize(500, directed=F)
147ego.net %v% 'sex' <- c(rep(0,250),rep(1,250))
148ego.deg <- c(180, 245, 60, 15)
149ego.mixmat <- matrix(c(164,     44, 26, 176)/2, nrow=2, byrow=T)
150ego.edges <- sum(ego.mixmat)
151ego.sexmatch <- ego.mixmat[1,1]+ego.mixmat[2,2]
152ego.isolates <- ego.deg[1]
153ego.conc <- ego.deg[3] + ego.deg[4]
154
155
156ego.target.stats <- c(ego.edges, ego.sexmatch)
157ego.deg
158ego.mixmat
159ego.target.stats
160
161ego.fit <- ergm(ego.net ~ edges + nodematch('sex'),
162 target.stats = ego.target.stats)
163summary(ego.fit)
164ego.sim1 <- simulate(ego.fit)
165summary(ego.sim1)
166plot(ego.sim1, vertex.cex=.65, vertex.col="sex")
167
168rbind(summary(ego.sim1 ~ degree(c(0:3))), ego.deg)
169mixingmatrix(ego.sim1, "sex")
170ego.mixmat
171
172# simulate 100 nets
173
174ego.sim100 <- simulate(ego.fit, nsim=100)
175ego.sim100
176summary(ego.sim100)
177sim.stats <- attr(ego.sim100,"stats")
178rbind(colMeans(sim.stats), ego.target.stats)
179matplot(1:nrow(sim.stats), sim.stats,
180  pch=c("e","m"), cex=.65,
181  main="100 simulations from ego.fit model", sub="(default settings)",
182  xlab="Replicate", ylab="frequency")
183abline(h=ego.target.stats, col=c(1:2))
184
185# increase MCMC interval
186
187ego.sim100 <- simulate(ego.fit, nsim=100,
188  control=control.simulate.ergm(MCMC.interval=10000))
189sim.stats <- attr(ego.sim100,"stats")
190matplot(1:nrow(sim.stats), sim.stats,
191  pch=c("e","m"), cex=.65,
192  main="100 simulations from ego.fit model", sub="(MCMC.interval=10000)",
193  xlab="Replicate", ylab="frequency")
194abline(h=ego.target.stats, col=c(1:2))
195
196# compare to non-targeted observed stats
197
198sim.fulldeg <- summary(ego.sim100 ~ degree(c(0:10)))
199sim.fulldeg
200sim.deg <- cbind(sim.fulldeg[,1:3], apply(sim.fulldeg[,4:11],1,sum))
201colnames(sim.deg) <- c(colnames(sim.fulldeg)[1:3],"degree3+")
202rbind(colMeans(sim.deg),ego.deg)
203matplot(1:nrow(sim.deg), sim.deg, pch=as.character(0:3), cex=.5,
204   main="Comparing ego.sims to non-targeted degree frequencies",
205   sub = "(edges+nodematch only)",
206   xlab = "Replicate", ylab = "Frequencies")
207abline(h=ego.deg, col=c(1:4))
208
209# add isolates to the model
210
211ego.isolates <- ego.deg[1]
212ego.target.stats <- c(ego.edges, ego.sexmatch, ego.isolates)
213ego.target.stats
214ego.fit <- ergm(ego.net ~ edges + nodematch('sex') + degree(0),
215 target.stats = ego.target.stats)
216summary(ego.fit)
217
218ego.sim100 <- simulate(ego.fit, nsim=100,
219  control=control.simulate.ergm(MCMC.interval=10000))
220sim.stats <- attr(ego.sim100,"stats")
221rbind(colMeans(sim.stats), ego.target.stats)
222
223
224# compare to non-targeted observed stats
225
226sim.fulldeg <- summary(ego.sim100 ~ degree(c(0:10)))
227sim.fulldeg
228sim.deg <- cbind(sim.fulldeg[,1:3], apply(sim.fulldeg[,4:11],1,sum))
229colnames(sim.deg) <- c(colnames(sim.fulldeg)[1:3],"degree3+")
230rbind(colMeans(sim.deg),ego.deg)
231
232matplot(1:nrow(sim.deg), sim.deg, pch=as.character(0:3), cex=.5,
233   main="Comparing ego.sims to non-targeted degree frequencies",
234   sub = "(edges + nodematch + isolates)",
235   xlab = "Replicate", ylab = "Frequencies")
236abline(h=ego.deg, col=c(1:4))
237
238save.image()