Sunbelt2015: tergm_tutorial.R

File tergm_tutorial.R, 8.6 KB (added by handcock, 4 years ago)

tergm R code

Line 
1
2## ----setup, include=FALSE------------------------------------------------
3library(tergm)
4library(knitr)
5knitr::opts_chunk$set(comment=NA)
6
7
8## ----,eval=FALSE---------------------------------------------------------
9## install.packages('statnet')
10## library(statnet)
11
12
13## ----,eval=FALSE---------------------------------------------------------
14## install.packages('tergm') # will install the network package
15## install.packages('sna')
16
17
18## ----,eval=FALSE---------------------------------------------------------
19## update.packages('name.of.package')
20
21
22## ----,eval=FALSE---------------------------------------------------------
23## install.packages('coda')
24
25
26## ----,results='hide', message=FALSE--------------------------------------
27library(statnet)
28
29
30## ----,results='hide', message=FALSE--------------------------------------
31library(tergm)
32library(sna)
33library(coda)
34
35
36## ----,eval=FALSE---------------------------------------------------------
37## # latest versions:  tergm 3.3, ergm 3.4.0, network 1.12.0, networkDynamic 0.7.1 (as of 6/19/2015)
38## sessionInfo()
39
40
41## ------------------------------------------------------------------------
42set.seed(0)
43
44
45## ----,results='hide'-----------------------------------------------------
46data("florentine")
47ls()
48
49
50## ------------------------------------------------------------------------
51plot(flobusiness)
52
53
54## ------------------------------------------------------------------------
55summary(flobusiness~edges+gwesp(0, fixed=T))
56
57
58## ------------------------------------------------------------------------
59fit1 <- ergm(flobusiness~edges+gwesp(0,fixed=T))
60summary(fit1)
61
62
63## ------------------------------------------------------------------------
64sim1 <- simulate(fit1,nsim=1,
65          control=control.simulate.ergm(MCMC.burnin=1000))
66plot(sim1)
67
68
69## ----,eval=FALSE---------------------------------------------------------
70##         stergm(my.network,
71##             formation= ~edges+kstar(2),
72##             dissolution= ~edges+triangle
73##         )
74
75
76## ------------------------------------------------------------------------
77data(samplk)
78ls()
79
80
81## ------------------------------------------------------------------------
82samp <- list()
83samp[[1]] <- samplk1
84samp[[2]] <- samplk2
85samp[[3]] <- samplk3
86
87
88## ------------------------------------------------------------------------
89plot(samplk1)
90
91
92## ----,results='hide'-----------------------------------------------------
93samp.fit <- stergm(samp,
94        formation= ~edges+mutual+cyclicalties+transitiveties,
95        dissolution = ~edges+mutual+cyclicalties+transitiveties,
96        estimate = "CMLE",
97  times=1:3
98        )
99
100
101## ------------------------------------------------------------------------
102summary(samp.fit)
103
104
105## ----,results='hide'-----------------------------------------------------
106samp.fit.2 <- stergm(samp,
107  formation= ~edges+mutual+cyclicalties+transitiveties,
108        dissolution = ~edges+mutual+cyclicalties+transitiveties,
109        estimate = "CMLE",
110  times=1:2
111        )
112
113
114## ------------------------------------------------------------------------
115theta.diss <- log(9)
116
117
118## ----,results='hide',message=FALSE,warning=FALSE-------------------------
119X11()
120stergm.fit.1 <- stergm(flobusiness,
121        formation= ~edges+gwesp(0,fixed=T),
122        dissolution = ~offset(edges),
123        targets="formation",
124        offset.coef.diss = theta.diss,
125        estimate = "EGMME",
126        control=control.stergm(SA.plot.progress=TRUE)
127        )
128dev.off()
129
130
131## ----,eval=FALSE---------------------------------------------------------
132## mcmc.diagnostics(stergm.fit.1)
133
134
135## ------------------------------------------------------------------------
136stergm.fit.1
137names(stergm.fit.1)
138stergm.fit.1$formation
139stergm.fit.1$formation.fit
140summary(stergm.fit.1)
141
142
143## ------------------------------------------------------------------------
144stergm.sim.1 <- simulate.stergm(stergm.fit.1, nsim=1,
145    time.slices = 1000)
146
147
148## ------------------------------------------------------------------------
149stergm.sim.1
150
151
152## ------------------------------------------------------------------------
153net500 <- network.collapse(stergm.sim.1,at=500)
154net500
155
156
157## ------------------------------------------------------------------------
158plot(net500)
159
160
161## ------------------------------------------------------------------------
162summary(flobusiness~edges+gwesp(0,fixed=T))
163colMeans(attributes(stergm.sim.1)$stats)
164
165
166## ------------------------------------------------------------------------
167plot(attributes(stergm.sim.1)$stats)
168
169
170## ------------------------------------------------------------------------
171plot(as.matrix(attributes(stergm.sim.1)$stats))
172
173
174## ------------------------------------------------------------------------
175stergm.sim.1.df <- as.data.frame(stergm.sim.1)
176names(stergm.sim.1.df)
177stergm.sim.1.df[1,]
178mean(stergm.sim.1.df$duration)
179
180
181## ----, eval=FALSE, show=FALSE--------------------------------------------
182## install.packages('ndtv')
183## library(ndtv)
184## stergm.sim.1a <- simulate.stergm(stergm.fit.1, nsim=1,
185##     time.slices = 100)
186## slice.par=list(start = 0, end = 25, interval=1, aggregate.dur=1, rule="any")
187## compute.animation(stergm.sim.1a, slice.par = slice.par)
188## render.par=list(tween.frames=5,show.time=T,
189##     show.stats="~edges+gwesp(0,fixed=T)")
190## wealthsize <- log(get.vertex.attribute(flobusiness, "wealth")) * 2/3
191## render.animation(stergm.sim.1a,render.par=render.par,
192##     edge.col="darkgray",displaylabels=T,
193##     label.cex=.8,label.col="blue",
194##     vertex.cex=wealthsize)
195## x11()
196## ani.replay()
197
198
199## ------------------------------------------------------------------------
200theta.diss.100 <- log(99)
201
202
203## ------------------------------------------------------------------------
204summary(fit1)
205theta.form <- fit1$coef
206theta.form
207
208
209
210## ------------------------------------------------------------------------
211theta.form[1] <- theta.form[1] - theta.diss.100
212
213
214## ------------------------------------------------------------------------
215stergm.sim.2 <- simulate(flobusiness, formation=~edges+gwesp(0,fixed=T),
216        dissolution=~edges, monitor="all",
217        coef.form=theta.form, coef.diss=theta.diss.100,
218        time.slices=50000)
219
220
221## ------------------------------------------------------------------------
222summary(flobusiness~edges+gwesp(0,fixed=T))
223colMeans(attributes(stergm.sim.2)$stats)
224stergm.sim.dm.2 <- as.data.frame(stergm.sim.2)
225mean(stergm.sim.dm.2$duration)
226plot(attributes(stergm.sim.2)$stats)
227
228
229## ------------------------------------------------------------------------
230msm.net <- network.initialize(500, directed=F) 
231msm.net %v% 'race' <- c(rep(0,250),rep(1,250))
232msm.net
233
234
235## ------------------------------------------------------------------------
236msm.form.formula <- ~edges+nodematch('race')+degree(0)+ concurrent
237msm.target.stats <- c(225,187,180,90)
238
239
240## ------------------------------------------------------------------------
241msm.diss.formula <- ~offset(edges)+offset(nodematch("race"))
242
243
244## ------------------------------------------------------------------------
245msm.theta.diss <- c(2.944, -0.747)
246
247
248## ----,results='hide',warning=FALSE---------------------------------------
249X11()
250msm.fit <- stergm(msm.net,
251        formation= msm.form.formula,
252        dissolution= msm.diss.formula,
253        targets="formation",
254        target.stats= msm.target.stats,
255        offset.coef.diss = msm.theta.diss,
256        estimate = "EGMME",
257        control=control.stergm(SA.plot.progress=TRUE,
258            SA.init.gain=0.005)
259)
260dev.off()
261
262
263## ------------------------------------------------------------------------
264mcmc.diagnostics(msm.fit)
265
266
267## ------------------------------------------------------------------------
268summary(msm.fit)
269
270
271## ----, warning=FALSE, message=FALSE--------------------------------------
272msm.sim <- simulate(msm.fit,time.slices=1000)
273
274
275## ------------------------------------------------------------------------
276colMeans(attributes(msm.sim)$stats)
277msm.target.stats
278
279
280## ------------------------------------------------------------------------
281race <- msm.net %v% 'race'
282msm.sim.dm <- as.data.frame(msm.sim)
283names(msm.sim.dm)
284mean(msm.sim.dm$duration[race[msm.sim.dm$tail] ==  race[msm.sim.dm$head]])
285mean(msm.sim.dm$duration[race[msm.sim.dm$tail] !=  race[msm.sim.dm$head]])
286
287
288## ----, eval=FALSE--------------------------------------------------------
289## slice.par=list(start = 0, end = 100, interval=1, aggregate.dur=1, rule="any")
290## compute.animation(msm.sim, slice.par = slice.par)
291##
292## render.par=list(tween.frames=5,show.time=T)
293## render.animation(msm.sim,render.par=render.par,
294##     edge.col="darkgray", displaylabels=F, vertex.col="race")
295## x11()
296## ani.replay()
297
298