Sunbelt2016: ergm_tutorial.R

File ergm_tutorial.R, 11.6 KB (added by morrism, 18 months ago)

ERGM tutorial R commands Sunbelt 2016

Line 
1## ----setup, include=FALSE------------------------------------------------
2library(ergm)
3library(knitr)
4knitr::opts_chunk$set(cache=F, comment=NA)
5
6## ----eval=FALSE----------------------------------------------------------
7## install.packages('statnet')
8## library(statnet)
9
10## ----eval=FALSE----------------------------------------------------------
11## install.packages('ergm') # will install the network package
12## install.packages('sna')
13
14## ----eval=FALSE----------------------------------------------------------
15## update.packages('name.of.package')
16
17## ----eval=FALSE----------------------------------------------------------
18## install.packages('coda')
19
20## ----results='hide', message=FALSE---------------------------------------
21library(statnet)
22
23## ----results='hide', message=FALSE---------------------------------------
24library(ergm)
25library(sna)
26library(coda)
27
28## ----eval=FALSE----------------------------------------------------------
29## # latest versions:  ergm 3.6.0 and network 1.13.0 (as of 4/03/2016)
30## sessionInfo()
31
32## ------------------------------------------------------------------------
33set.seed(0)
34
35## ------------------------------------------------------------------------
36data(package='ergm') # tells us the datasets in our packages
37
38## ------------------------------------------------------------------------
39data(florentine) # loads flomarriage and flobusiness data
40flomarriage # Let's look at the flomarriage network properties
41par(mfrow=c(1,2)) # Setup a 2 panel plot (for later)
42plot(flomarriage, main="Florentine Marriage", cex.main=0.8) # Plot the flomarriage network
43summary(flomarriage~edges) # Look at the $g(y)$ statistic for this model
44flomodel.01 <- ergm(flomarriage~edges) # Estimate the model
45summary(flomodel.01) # The fitted model object
46
47## ------------------------------------------------------------------------
48summary(flomarriage~edges+triangle) # Look at the g(y) stats for this model
49flomodel.02 <- ergm(flomarriage~edges+triangle)
50summary(flomodel.02)
51
52## ------------------------------------------------------------------------
53class(flomodel.02) # this has the class ergm
54
55names(flomodel.02) # the ERGM object contains lots of components.
56
57## ------------------------------------------------------------------------
58flomodel.02$coef # you can extract/inspect individual components
59
60## ------------------------------------------------------------------------
61wealth <- flomarriage %v% 'wealth' # %v% references vertex attributes
62wealth
63summary(wealth) # summarize the distribution of wealth
64plot(flomarriage, vertex.cex=wealth/25, main="Florentine marriage by wealth", cex.main=0.8) # network plot with vertex size proportional to wealth
65summary(flomarriage~edges+nodecov('wealth'))
66# observed statistics for the model
67flomodel.03 <- ergm(flomarriage~edges+nodecov('wealth'))
68summary(flomodel.03)
69
70## ------------------------------------------------------------------------
71data(faux.mesa.high)
72mesa <- faux.mesa.high
73
74## ------------------------------------------------------------------------
75mesa
76par(mfrow=c(1,1)) # Back to 1-panel plots
77plot(mesa, vertex.col='Grade')
78legend('bottomleft',fill=7:12,legend=paste('Grade',7:12),cex=0.75)
79
80## ------------------------------------------------------------------------
81fauxmodel.01 <- ergm(mesa ~edges + nodematch('Grade',diff=T) + nodematch('Race',diff=T))
82summary(fauxmodel.01)
83
84## ------------------------------------------------------------------------
85table(mesa %v% 'Race') # Frequencies of race
86mixingmatrix(mesa, "Race")
87
88## ------------------------------------------------------------------------
89summary(mesa ~edges + nodematch('Grade',diff=T) + nodematch('Race',diff=T))
90
91## ------------------------------------------------------------------------
92data(samplk)
93ls() # directed data: Sampson's Monks
94samplk3
95plot(samplk3)
96summary(samplk3~edges+mutual)
97
98## ------------------------------------------------------------------------
99sampmodel.01 <- ergm(samplk3~edges+mutual)
100summary(sampmodel.01)
101
102## ------------------------------------------------------------------------
103missnet <- network.initialize(10,directed=F)
104missnet[1,2] <- missnet[2,7] <- missnet[3,6] <- 1
105missnet[4,6] <- missnet[4,9] <- missnet[5,6] <- NA
106summary(missnet)
107
108# plot missnet with missing edge colored red.
109tempnet <- missnet
110tempnet[4,6] <- tempnet[4,9] <- tempnet[5,6] <- 1
111missnetmat <- as.matrix(missnet)
112missnetmat[is.na(missnetmat)] <- 2
113plot(tempnet,label = network.vertex.names(tempnet),edge.col = missnetmat)
114
115summary(missnet~edges)
116summary(ergm(missnet~edges))
117
118## ------------------------------------------------------------------------
119missnet_bad <- missnet
120missnet_bad[4,6] <- missnet_bad[4,9] <- missnet_bad[5,6] <- 0
121summary(missnet_bad)
122summary(ergm(missnet_bad~edges))
123
124## ----eval=FALSE----------------------------------------------------------
125## help('ergm-terms')
126
127## ------------------------------------------------------------------------
128flomodel.03.sim <- simulate(flomodel.03,nsim=10)
129class(flomodel.03.sim)
130summary(flomodel.03.sim)
131length(flomodel.03.sim)
132flomodel.03.sim[[1]]
133plot(flomodel.03.sim[[1]], label= flomodel.03.sim[[1]] %v% "vertex.names")
134
135## ------------------------------------------------------------------------
136flo.03.gof.model <- gof(flomodel.03 ~ model)
137flo.03.gof.model
138plot(flo.03.gof.model)
139
140## ------------------------------------------------------------------------
141flo.03.gof.global <- gof(flomodel.03 ~ degree + esp + distance)
142flo.03.gof.global
143plot(flo.03.gof.global)
144
145## ------------------------------------------------------------------------
146mesamodel.b <- ergm(mesa~edges)
147plot(gof(mesamodel.b ~ model, nsim=10))
148plot(gof(mesamodel.b ~ degree + esp + distance, nsim=10))
149
150## ------------------------------------------------------------------------
151summary(flobusiness ~ edges+degree(1))
152fit <- ergm(flobusiness ~ edges+degree(1))
153mcmc.diagnostics(fit)
154
155## ---- eval=FALSE---------------------------------------------------------
156## ergm(flobusiness ~ edges+degree(1),
157##      control=control.ergm(MCMC.interval=1))
158
159## ------------------------------------------------------------------------
160data('faux.magnolia.high')
161magnolia <- faux.magnolia.high
162plot(magnolia, vertex.cex=.5)
163summary(magnolia ~ edges+triangle)
164
165## ---- eval=F-------------------------------------------------------------
166## ergm(magnolia ~ edges+triangle)
167
168## ----eval=T--------------------------------------------------------------
169fit.mag.01 <- ergm(magnolia ~ edges+triangle, control=control.ergm(MCMLE.maxit=2))
170
171
172## ---- eval=T, results='hide',fig.show='asis'-----------------------------
173mcmc.diagnostics(fit.mag.01)
174
175## ------------------------------------------------------------------------
176fit.mag.02 <- ergm(magnolia ~ edges+gwesp(0.25,fixed=T))
177mcmc.diagnostics(fit.mag.02)
178
179## ------------------------------------------------------------------------
180gof.mag.02.model <- gof(fit.mag.02, GOF = ~model)
181gof.mag.02.model
182plot(gof.mag.02.model)
183
184
185## ------------------------------------------------------------------------
186fit.mag.03 <- ergm(magnolia ~ edges+gwesp(0.25,fixed=T)
187                   +nodematch('Grade')+nodematch('Race')+nodematch('Sex'),
188               control = control.ergm(MCMC.interval=20000), eval.loglik = F)
189summary(fit.mag.03)
190
191## ------------------------------------------------------------------------
192mcmc.diagnostics(fit.mag.03)
193
194## ------------------------------------------------------------------------
195plot(gof(fit.mag.03, GOF=~model, control=control.gof.ergm(nsim=200)))
196
197## ------------------------------------------------------------------------
198plot(gof(fit.mag.03, GOF = ~ degree + esp + distance))
199
200## ------------------------------------------------------------------------
201ego.net <- network.initialize(500, directed=F)
202ego.net %v% 'sex' <- c(rep(0,250),rep(1,250))
203
204## ------------------------------------------------------------------------
205ego.deg <- c(180, 245, 60, 15)                          # node distn
206ego.mixmat <- matrix(c(164,44,26,176)/2, nrow=2, byrow=T)       # adjusted tie distn
207
208## ------------------------------------------------------------------------
209ego.edges <- sum(ego.mixmat)
210ego.sexmatch <- ego.mixmat[1,1]+ego.mixmat[2,2]
211
212## ------------------------------------------------------------------------
213ego.target.stats <- c(ego.edges, ego.sexmatch)
214ego.target.stats
215
216## ------------------------------------------------------------------------
217ego.fit <- ergm(ego.net ~ edges + nodematch('sex'),
218                target.stats = ego.target.stats)
219
220## ------------------------------------------------------------------------
221summary(ego.fit)
222
223## ------------------------------------------------------------------------
224ego.sim1 <- simulate(ego.fit)
225plot(ego.sim1, vertex.cex=.65, vertex.col="sex")
226
227## ------------------------------------------------------------------------
228rbind(sim=summary(ego.sim1 ~ degree(c(0:3))), obs=ego.deg)
229mixingmatrix(ego.sim1, "sex")
230ego.mixmat
231
232## ------------------------------------------------------------------------
233ego.sim100 <- simulate(ego.fit, nsim=100)
234ego.sim100
235
236## ------------------------------------------------------------------------
237sim.stats <- attr(ego.sim100,"stats")
238rbind(sim=colMeans(sim.stats), obs=ego.target.stats)
239
240## ------------------------------------------------------------------------
241matplot(1:nrow(sim.stats), sim.stats,
242        pch=c("e","m","0","+"), cex=.65,
243        main="100 simulations from ego.fit model", sub="(default settings)",
244        xlab="Replicate", ylab="frequency")
245abline(h=ego.target.stats, col=c(1:4))
246
247## ------------------------------------------------------------------------
248sim.fulldeg <- summary(ego.sim100 ~ degree(c(0:10)))
249colnames(sim.fulldeg) <- paste("deg",0:10, sep='')
250sim.fulldeg[1:5,]
251
252## ------------------------------------------------------------------------
253sim.deg <- cbind(sim.fulldeg[,1:3], apply(sim.fulldeg[,4:11],1,sum))
254colnames(sim.deg) <- c(colnames(sim.fulldeg)[1:3],"degree3+")
255rbind(sim=colMeans(sim.deg), obs=ego.deg)
256
257## ------------------------------------------------------------------------
258matplot(1:nrow(sim.deg), sim.deg, pch=as.character(0:3), cex=.5,
259        main="Comparing ego.sims to non-targeted degree frequencies",
260        sub = "(only total edges targeted)",
261        xlab = "Replicate", ylab = "Frequencies")
262abline(h=c(180, 245, 60, 15), col=c(1:4))
263
264## ------------------------------------------------------------------------
265ego.isolates <- ego.deg[1]
266ego.target.stats <- c(ego.edges, ego.sexmatch, ego.isolates)
267ego.fit <- ergm(ego.net ~ edges + nodematch('sex') + degree(0),
268                target.stats = ego.target.stats)
269summary(ego.fit)
270
271## ------------------------------------------------------------------------
272ego.sim100 <- simulate(ego.fit, nsim=100,
273                       control=control.simulate.ergm(MCMC.interval=10000))
274sim.stats <- attr(ego.sim100,"stats")
275rbind(sim=colMeans(sim.stats), obs=ego.target.stats)
276
277## ------------------------------------------------------------------------
278sim.fulldeg <- summary(ego.sim100 ~ degree(c(0:10)))
279sim.deg <- cbind(sim.fulldeg[,1:3], apply(sim.fulldeg[,4:11],1,sum))
280colnames(sim.deg) <- c(colnames(sim.fulldeg)[1:3],"degree3+")
281rbind(sim=colMeans(sim.deg), obs=ego.deg)
282
283## ------------------------------------------------------------------------
284matplot(1:nrow(sim.deg), sim.deg, pch=as.character(0:3), cex=.5,
285        main="Comparing ego.sims to non-targeted degree frequencies",
286        sub = "(only 0, 2+ and total edges targeted)",
287        xlab = "Replicate", ylab = "Frequencies")
288abline(h=c(180, 245, 60, 15), col=c(1:4))
289