Sunbelt2016: tergm_tutorial.R

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