Sunbelt2015: ergm_tutorial.R

File ergm_tutorial.R, 11.0 KB (added by kirkli, 4 years ago)

ergm R code

Line 
1## ----setup, include=FALSE------------------------------------------------
2library(ergm)
3library(knitr)
4knitr::opts_chunk$set(cache=T, 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.1.2 and network 1.9.0 (as of 6/17/2014)
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')) # observed statistics for the model
66flomodel.03 <- ergm(flomarriage~edges+nodecov('wealth'))
67summary(flomodel.03)
68
69## ------------------------------------------------------------------------
70data(faux.mesa.high)
71mesa <- faux.mesa.high
72
73## ------------------------------------------------------------------------
74mesa
75par(mfrow=c(1,1)) # Back to 1-panel plots
76plot(mesa, vertex.col='Grade')
77legend('bottomleft',fill=7:12,legend=paste('Grade',7:12),cex=0.75)
78
79## ------------------------------------------------------------------------
80fauxmodel.01 <- ergm(mesa ~edges + nodematch('Grade',diff=T) + nodematch('Race',diff=T))
81summary(fauxmodel.01)
82
83## ------------------------------------------------------------------------
84table(mesa %v% 'Race') # Frequencies of race
85mixingmatrix(mesa, "Race")
86
87## ------------------------------------------------------------------------
88summary(mesa ~edges + nodematch('Grade',diff=T) + nodematch('Race',diff=T))
89
90## ------------------------------------------------------------------------
91data(samplk)
92ls() # directed data: Sampson's Monks
93samplk3
94plot(samplk3)
95summary(samplk3~edges+mutual)
96
97## ------------------------------------------------------------------------
98sampmodel.01 <- ergm(samplk3~edges+mutual)
99summary(sampmodel.01)
100
101## ------------------------------------------------------------------------
102missnet <- network.initialize(10,directed=F)
103missnet[1,2] <- missnet[2,7] <- missnet[3,6] <- 1
104missnet[4,6] <- missnet[4,9] <- missnet[5,6] <- NA
105summary(missnet)
106
107# plot missnet with missing edge colored red.
108tempnet <- missnet
109tempnet[4,6] <- tempnet[4,9] <- tempnet[5,6] <- 1
110missnetmat <- as.matrix(missnet)
111missnetmat[is.na(missnetmat)] <- 2
112plot(tempnet,label = network.vertex.names(tempnet),edge.col = missnetmat)
113
114summary(missnet~edges)
115summary(ergm(missnet~edges))
116
117## ------------------------------------------------------------------------
118missnet_bad <- missnet
119missnet_bad[4,6] <- missnet_bad[4,9] <- missnet_bad[5,6] <- 0
120summary(missnet_bad)
121summary(ergm(missnet_bad~edges))
122
123## ----eval=FALSE----------------------------------------------------------
124## help('ergm-terms')
125
126## ------------------------------------------------------------------------
127flomodel.03.sim <- simulate(flomodel.03,nsim=10)
128class(flomodel.03.sim)
129summary(flomodel.03.sim)
130length(flomodel.03.sim)
131flomodel.03.sim[[1]]
132plot(flomodel.03.sim[[1]], label= flomodel.03.sim[[1]] %v% "vertex.names")
133
134## ------------------------------------------------------------------------
135flomodel.03.gof <- gof(flomodel.03~degree + esp + distance)
136flomodel.03.gof
137plot(flomodel.03.gof)
138
139## ------------------------------------------------------------------------
140mesamodel.02 <- ergm(mesa~edges)
141mesamodel.02.gof <- gof(mesamodel.02~degree + esp + distance, nsim=10)
142plot(mesamodel.02.gof)
143
144
145## ------------------------------------------------------------------------
146summary(flobusiness~edges+degree(1))
147fit <- ergm(flobusiness~edges+degree(1))
148mcmc.diagnostics(fit)
149
150## ---- eval=FALSE---------------------------------------------------------
151## fit <- ergm(flobusiness~edges+degree(1),
152## control=control.ergm(MCMC.interval=1)
153
154## ------------------------------------------------------------------------
155data('faux.magnolia.high')
156magnolia <- faux.magnolia.high
157plot(magnolia, vertex.cex=.5)
158summary(magnolia~edges+triangle)
159
160## ---- eval=F-------------------------------------------------------------
161## fit <- ergm(magnolia~edges+triangle)
162
163## ----eval=T--------------------------------------------------------------
164fit <- ergm(magnolia~edges+triangle, control=control.ergm(MCMLE.maxit=2))
165
166
167## ---- eval=T, results='hide',fig.show='asis'-----------------------------
168mcmc.diagnostics(fit)
169
170## ------------------------------------------------------------------------
171fit <- ergm(magnolia~edges+gwesp(0.25,fixed=T),verbose=T)
172mcmc.diagnostics(fit)
173
174## ------------------------------------------------------------------------
175fit <- ergm(magnolia~edges+gwesp(0.25,fixed=T)+nodematch('Grade')+
176              nodematch('Race')+nodematch('Sex'),
177            control = control.ergm(MCMC.samplesize=50000,MCMC.interval=1000),
178            verbose=T)
179
180## ------------------------------------------------------------------------
181mcmc.diagnostics(fit)
182
183## ------------------------------------------------------------------------
184ego.net <- network.initialize(500, directed=F)
185ego.net %v% 'sex' <- c(rep(0,250),rep(1,250))
186
187## ------------------------------------------------------------------------
188ego.deg <- c(180, 245, 60, 15)                          # node distn
189ego.mixmat <- matrix(c(164,44,26,176)/2, nrow=2, byrow=T)       # adjusted tie distn
190
191## ------------------------------------------------------------------------
192ego.edges <- sum(ego.mixmat)
193ego.sexmatch <- ego.mixmat[1,1]+ego.mixmat[2,2]
194
195## ------------------------------------------------------------------------
196ego.target.stats <- c(ego.edges, ego.sexmatch)
197ego.target.stats
198
199## ------------------------------------------------------------------------
200ego.fit <- ergm(ego.net ~ edges + nodematch('sex'),
201                target.stats = ego.target.stats)
202
203## ------------------------------------------------------------------------
204summary(ego.fit)
205
206## ------------------------------------------------------------------------
207ego.sim1 <- simulate(ego.fit)
208plot(ego.sim1, vertex.cex=.65, vertex.col="sex")
209
210## ------------------------------------------------------------------------
211rbind(sim=summary(ego.sim1 ~ degree(c(0:3))), obs=ego.deg)
212mixingmatrix(ego.sim1, "sex")
213ego.mixmat
214
215## ------------------------------------------------------------------------
216ego.sim100 <- simulate(ego.fit, nsim=100)
217ego.sim100
218
219## ------------------------------------------------------------------------
220sim.stats <- attr(ego.sim100,"stats")
221rbind(sim=colMeans(sim.stats), obs=ego.target.stats)
222
223## ------------------------------------------------------------------------
224matplot(1:nrow(sim.stats), sim.stats,
225        pch=c("e","m","0","+"), cex=.65,
226        main="100 simulations from ego.fit model", sub="(default settings)",
227        xlab="Replicate", ylab="frequency")
228abline(h=ego.target.stats, col=c(1:4))
229
230## ------------------------------------------------------------------------
231sim.fulldeg <- summary(ego.sim100 ~ degree(c(0:10)))
232colnames(sim.fulldeg) <- paste("deg",0:10, sep='')
233sim.fulldeg[1:5,]
234
235## ------------------------------------------------------------------------
236sim.deg <- cbind(sim.fulldeg[,1:3], apply(sim.fulldeg[,4:11],1,sum))
237colnames(sim.deg) <- c(colnames(sim.fulldeg)[1:3],"degree3+")
238rbind(sim=colMeans(sim.deg), obs=ego.deg)
239
240## ------------------------------------------------------------------------
241matplot(1:nrow(sim.deg), sim.deg, pch=as.character(0:3), cex=.5,
242        main="Comparing ego.sims to non-targeted degree frequencies",
243        sub = "(only total edges targeted)",
244        xlab = "Replicate", ylab = "Frequencies")
245abline(h=c(180, 245, 60, 15), col=c(1:4))
246
247## ------------------------------------------------------------------------
248ego.isolates <- ego.deg[1]
249ego.target.stats <- c(ego.edges, ego.sexmatch, ego.isolates)
250ego.fit <- ergm(ego.net ~ edges + nodematch('sex') + degree(0),
251                target.stats = ego.target.stats)
252summary(ego.fit)
253
254## ------------------------------------------------------------------------
255ego.sim100 <- simulate(ego.fit, nsim=100,
256                       control=control.simulate.ergm(MCMC.interval=10000))
257sim.stats <- attr(ego.sim100,"stats")
258rbind(sim=colMeans(sim.stats), obs=ego.target.stats)
259
260## ------------------------------------------------------------------------
261sim.fulldeg <- summary(ego.sim100 ~ degree(c(0:10)))
262sim.deg <- cbind(sim.fulldeg[,1:3], apply(sim.fulldeg[,4:11],1,sum))
263colnames(sim.deg) <- c(colnames(sim.fulldeg)[1:3],"degree3+")
264rbind(sim=colMeans(sim.deg), obs=ego.deg)
265
266## ------------------------------------------------------------------------
267matplot(1:nrow(sim.deg), sim.deg, pch=as.character(0:3), cex=.5,
268        main="Comparing ego.sims to non-targeted degree frequencies",
269        sub = "(only 0, 2+ and total edges targeted)",
270        xlab = "Replicate", ylab = "Frequencies")
271abline(h=c(180, 245, 60, 15), col=c(1:4))
272