Sunbelt2015: Valued.R

File Valued.R, 13.1 KB (added by krivitsky, 4 years ago)
Line 
1## ----setup, cache=FALSE, echo=FALSE, results='hide'----------------------
2library(knitr)
3opts_chunk$set(echo=TRUE,tidy=TRUE,cache=TRUE,autodep=TRUE)
4options(width=60, useFancyQuotes = FALSE, continue="  ")
5
6## ----install packages, eval=FALSE----------------------------------------
7## install.packages("ergm.count")
8## install.packages("latentnet")
9
10## ----load libraries, results='hide', cache=FALSE-------------------------
11library(ergm.count)
12library(latentnet)
13
14## ----constructing from a sociomatrix-------------------------------------
15data(samplk)
16ls()
17as.matrix(samplk1)[1:5,1:5]
18# Create a sociomatrix totaling the nominations.
19samplk.tot.m<-as.matrix(samplk1)+as.matrix(samplk2)+as.matrix(samplk3)
20samplk.tot.m[1:5,1:5]
21
22# Create a network where the number of nominations becomes an attribute of an edge.
23samplk.tot <- as.network(samplk.tot.m, directed=TRUE, matrix.type="a",
24                           ignore.eval=FALSE, names.eval="nominations" # Important!
25                           )
26# Add vertex attributes.  (Note that names were already imported!)
27samplk.tot %v% "group" <- samplk1 %v% "group" # Groups identified by Sampson
28samplk.tot %v% "group"
29
30# We can view the attribute as a sociomatrix.
31as.matrix(samplk.tot,attrname="nominations")[1:5,1:5]
32
33# Also, note that samplk.tot now has an edge if i nominated j *at least once*.
34as.matrix(samplk.tot)[1:5,1:5]
35
36## ----constructing from an edgelist---------------------------------------
37samplk.tot.el <- as.matrix(samplk.tot, attrname="nominations",
38                           matrix.type="edgelist")
39samplk.tot.el[1:5,]
40# and an initial empty network.
41samplk.tot2 <- samplk1 # Copy samplk1
42delete.edges(samplk.tot2, seq_along(samplk.tot2$mel)) # Empty it out
43samplk.tot2  #We could also have used network.initialize(18)
44
45samplk.tot2[samplk.tot.el[,1:2], names.eval="nominations", add.edges=TRUE] <-
46  samplk.tot.el[,3]
47as.matrix(samplk.tot2,attrname="nominations")[1:5,1:5]
48
49## ----plotting via matrix-------------------------------------------------
50par(mar=rep(0,4))
51samplk.ecol <-
52  matrix(gray(1 - (as.matrix(samplk.tot, attrname="nominations")/3)),
53         nrow=network.size(samplk.tot))
54plot(samplk.tot, edge.col=samplk.ecol, usecurve=TRUE, edge.curve=0.0001,
55     displaylabels=TRUE, vertex.col=as.factor(samplk.tot%v%"group"))
56
57## ----plotting via edge values--------------------------------------------
58par(mar=rep(0,4))
59valmat<-as.matrix(samplk.tot,attrname="nominations") #Pull the edge values
60samplk.ecol <-
61  matrix(rgb(0,0,0,valmat/3),
62         nrow=network.size(samplk.tot))
63plot(samplk.tot, edge.col=samplk.ecol, usecurve=TRUE, edge.curve=0.0001,
64     displaylabels=TRUE, vertex.col=as.factor(samplk.tot%v%"group"),
65     edge.lwd=valmat^2)
66
67## ----zach plot-----------------------------------------------------------
68data(zach)
69zach.ecol <- gray(1 - (zach %e% "contexts")/8)
70zach.vcol <- rainbow(5)[zach %v% "faction.id"+3]
71par(mar=rep(0,4))
72plot(zach, edge.col=zach.ecol, vertex.col=zach.vcol, displaylabels=TRUE)
73
74## ----valued ERGM summary error, eval=FALSE-------------------------------
75## summary(samplk.tot~sum)
76
77## ----valued ERGM summary-------------------------------------------------
78summary(samplk.tot~sum, response="nominations")
79
80## ----help ergm-references,eval=FALSE-------------------------------------
81## help("ergm-references")
82
83## ----DiscUnif reference,fig.height=3-------------------------------------
84y <- network.initialize(2,directed=FALSE) # A network with one dyad!
85## Discrete Uniform reference
86# 0 coefficient: discrete uniform
87sim.du3<-simulate(y~sum, coef=0, reference=~DiscUnif(0,3),
88                  response="w",statsonly=TRUE,nsim=1000)
89# Negative coefficient: truncated geometric skewed to the right
90sim.trgeo.m1<-simulate(y~sum, coef=-1, reference=~DiscUnif(0,3),
91                       response="w",statsonly=TRUE,nsim=1000)
92# Positive coefficient: truncated geometric skewed to the left
93sim.trgeo.p1<-simulate(y~sum, coef=+1, reference=~DiscUnif(0,3),
94                      response="w",statsonly=TRUE,nsim=1000)
95# Plot them:
96par(mfrow=c(1,3))
97hist(sim.du3,breaks=diff(range(sim.du3))*4)
98hist(sim.trgeo.m1,breaks=diff(range(sim.trgeo.m1))*4)
99hist(sim.trgeo.p1,breaks=diff(range(sim.trgeo.p1))*4)
100
101## ----Binomial reference,fig.height=3-------------------------------------
102## Binomial reference
103# 0 coefficient: Binomial(3,1/2)
104sim.binom3<-simulate(y~sum, coef=0, reference=~Binomial(3),
105                     response="w",statsonly=TRUE,nsim=1000)
106# -1 coefficient: Binomial(3, exp(-1)/(1+exp(-1)))
107sim.binom3.m1<-simulate(y~sum, coef=-1, reference=~Binomial(3),
108                        response="w",statsonly=TRUE,nsim=1000)
109# +1 coefficient: Binomial(3, exp(1)/(1+exp(1)))
110sim.binom3.p1<-simulate(y~sum, coef=+1, reference=~Binomial(3),
111                        response="w",statsonly=TRUE,nsim=1000)
112# Plot them:
113par(mfrow=c(1,3))
114hist(sim.binom3,breaks=diff(range(sim.binom3))*4)
115hist(sim.binom3.m1,breaks=diff(range(sim.binom3.m1))*4)
116hist(sim.binom3.p1,breaks=diff(range(sim.binom3.p1))*4)
117
118## ----unbounded references------------------------------------------------
119sim.geom<-simulate(y~sum, coef=log(2/3), reference=~Geometric,
120                   response="w",statsonly=TRUE,nsim=1000)
121mean(sim.geom)
122sim.pois<-simulate(y~sum, coef=log(2), reference=~Poisson,
123                   response="w",statsonly=TRUE,nsim=1000)
124mean(sim.pois)
125
126## ----unbounded histograms, fig.height=3----------------------------------
127par(mfrow=c(1,2))
128hist(sim.geom,breaks=diff(range(sim.geom))*4)
129hist(sim.pois,breaks=diff(range(sim.pois))*4)
130
131## ----outside parameter space, fig.height=3-------------------------------
132par(mfrow=c(1,1))
133sim.geo0<-simulate(y~sum, coef=0, reference=~Geometric,
134                    response="w",statsonly=TRUE,nsim=100,
135                    control=control.simulate(MCMC.burnin=0,MCMC.interval=1))
136mean(sim.geo0)
137plot(c(sim.geo0),xlab="MCMC iteration",ylab="Value of the tie")
138
139## ----help ergm-terms, eval=FALSE-----------------------------------------
140## help("ergm-terms")
141
142## ----monks total fit, results='hide'-------------------------------------
143samplk.tot.nm <-
144  ergm(samplk.tot~sum + nodematch("group",diff=TRUE,form="sum"),
145       response="nominations", reference=~Binomial(3))
146mcmc.diagnostics(samplk.tot.nm)
147
148## ----monks total fit summ------------------------------------------------
149summary(samplk.tot.nm)
150
151## ----create leader attribute---------------------------------------------
152unique(zach %v% "role")
153# Vertex attr. "leader" is TRUE for Hi and John, FALSE for others.
154zach %v% "leader" <- zach %v% "role" != "Member"
155
156## ----zach sum leader fit, results='hide'---------------------------------
157zach.lead <-
158  ergm(zach~sum + nodefactor("leader"),
159       response="contexts", reference=~Poisson)
160mcmc.diagnostics(zach.lead)
161
162## ----zach sum leader fit summ--------------------------------------------
163summary(zach.lead)
164
165## ----zach sum nonzero group fit, results='hide'--------------------------
166samplk.tot.nm.nz <-
167  ergm(samplk.tot~sum + nonzero + nodematch("group",diff=TRUE,form="sum"),
168       response="nominations", reference=~Binomial(3))
169mcmc.diagnostics(samplk.tot.nm.nz)
170
171## ----zach sum nonzero group fit summ-------------------------------------
172summary(samplk.tot.nm.nz)
173
174## ----monks big fit, results='hide'---------------------------------------
175samplk.tot.ergm <-
176  ergm(samplk.tot ~ sum + nonzero + mutual("min") +
177       transitiveweights("min","max","min") +
178       cyclicalweights("min","max","min"),
179       reference=~Binomial(3), response="nominations")
180mcmc.diagnostics(samplk.tot.ergm)
181
182## ----monks big fit summ--------------------------------------------------
183summary(samplk.tot.ergm)
184
185## ----zach big model------------------------------------------------------
186summary(zach~sum+nonzero+nodefactor("leader")+absdiffcat("faction.id")+
187        nodesqrtcovar(TRUE), response="contexts")
188
189## ----zach big fit, results='hide'----------------------------------------
190zach.pois <-
191  ergm(zach~sum+nonzero+nodefactor("leader")+absdiffcat("faction.id")+
192       nodesqrtcovar(TRUE),
193       response="contexts", reference=~Poisson,
194       control=control.ergm(MCMLE.trustregion=100, MCMLE.maxit=50), verbose=TRUE)
195mcmc.diagnostics(zach.pois)
196
197## ----zach big fit summ---------------------------------------------------
198summary(zach.pois)
199
200## ----zach big fit sim, results='hide'------------------------------------
201# Simulate from model fit:
202zach.sim <-
203  simulate(zach.pois, monitor=~transitiveweights("geomean","sum","geomean"),
204           nsim = 1000, statsonly=TRUE)
205
206## ----zach big fit sim summ-----------------------------------------------
207# What have we simulated?
208colnames(zach.sim)
209
210# How high is the transitiveweights statistic in the observed network?
211zach.obs <- summary(zach ~ transitiveweights("geomean","sum","geomean"),
212                    response="contexts")
213zach.obs
214
215## ----zach big fit sim plot,fig.height=3----------------------------------
216par(mar=c(5, 4, 4, 2) + 0.1)
217# 9th col. = transitiveweights
218plot(density(zach.sim[,9]))
219abline(v=zach.obs)
220
221## ----zach big fit sim test-----------------------------------------------
222# Where does the observed value lie in the simulated?
223# This is a p-value for the Monte-Carlo test:
224min(mean(zach.sim[,9]>zach.obs), mean(zach.sim[,9]<zach.obs))*2
225
226## ----eval=FALSE----------------------------------------------------------
227## install.packages('ergm.rank',repos='http://statnet.csde.washington.edu')
228
229## ----setup ergm.rank, cache=FALSE----------------------------------------
230library(ergm.rank)
231
232## ----references ergm.rank, eval=FALSE------------------------------------
233## help("ergm-references", "ergm.rank")
234
235## ----newcomb data--------------------------------------------------------
236data(newcomb)
237as.matrix(newcomb[[1]], attrname="rank")
238as.matrix(newcomb[[1]], attrname="descrank")
239
240## ----newcomb fit1, results='hide'----------------------------------------
241newc.fit1<- ergm(newcomb[[1]]~rank.nonconformity+rank.nonconformity("localAND")+rank.deference,response="descrank",reference=~CompleteOrder,control=control.ergm(MCMLE.trustregion=1000, MCMC.burnin=4096, MCMC.interval=32, CD.conv.min.pval=0.05),eval.loglik=FALSE)
242
243## ----newcomb fit1 summ---------------------------------------------------
244summary(newc.fit1)
245
246## ----newcomb fit1 mcmc, results='hide'-----------------------------------
247mcmc.diagnostics(newc.fit1)
248
249## ----newcomb fit15, results='hide'---------------------------------------
250newc.fit15 <- ergm(newcomb[[15]]~rank.nonconformity+rank.nonconformity("localAND")+rank.deference,response="descrank",reference=~CompleteOrder,control=control.ergm(MCMLE.trustregion=1000, MCMC.burnin=4096, MCMC.interval=32, CD.conv.min.pval=0.05),eval.loglik=FALSE)
251
252## ----newcomb fit15 summ--------------------------------------------------
253summary(newc.fit15)
254
255## ----newcomb fit15 mcmc, results='hide'----------------------------------
256mcmc.diagnostics(newc.fit15)
257
258## ----monks nodematch mle fit, results='hide'-----------------------------
259samplk.nm.l <- ergmm(samplk.tot~nodematch("group",diff=TRUE),tofit="mle", verbose=TRUE)
260
261## ----monks nodematch ergm fit, results='hide'----------------------------
262samplk.nm.e <- ergm(samplk.tot~edges+nodematch("group",diff=TRUE))
263
264## ----monks nodematch compare---------------------------------------------
265summary(samplk.nm.l, point.est="mle", se=TRUE)
266summary(samplk.nm.e)
267
268## ----monks latent fit, results='hide'------------------------------------
269samplk.d2G3<-ergmm(samplk.tot~euclidean(d=2,G=3), verbose=TRUE)
270
271## ----monks latent re fit, results='hide'---------------------------------
272samplk.d2G3r<-ergmm(samplk.tot~euclidean(d=2,G=3)+rreceiver, verbose=TRUE)
273mcmc.diagnostics(samplk.d2G3r)
274
275## ----monks latent re fit plot, fig.height=3------------------------------
276par(mfrow=c(1,2))
277# Extract a clustering
278Z.K.ref <- summary(samplk.d2G3,point.est="pmean")$pmean$Z.K
279# Plot one model, saving positions, using Z.K.ref to set reference clustering.
280Z.ref <- plot(samplk.d2G3, pie=TRUE, Z.K.ref=Z.K.ref)
281# Plot the other model, using Z.ref and Z.K.ref to ensure similar
282# orientation and coloring.
283plot(samplk.d2G3r, rand.eff="receiver", pie=TRUE, Z.ref=Z.ref, Z.K.ref=Z.K.ref)
284
285## ----eval=FALSE----------------------------------------------------------
286## ? families.ergmm
287
288## ----monks count latent fit,results='hide'-------------------------------
289# Bernoulli logit fit (recall)
290# samplk.d2G3 <- ergmm(samplk.tot~euclidean(d=2,G=3))
291# Binomial(trials=3) logit fit
292samplk.ct.d2G3 <-
293  ergmm(samplk.tot~euclidean(d=2,G=3),response="nominations",
294        family="binomial",fam.par=list(trials=3), verbose=TRUE)
295
296## ----monks count latent fit plot, fig.height=3---------------------------
297# Plot them side-by-side, using edge.col argument:
298par(mfrow=c(1,2))
299plot(samplk.d2G3,pie=TRUE, Z.ref=Z.ref, Z.K.ref=Z.K.ref)
300plot(samplk.ct.d2G3,pie=TRUE, Z.ref=Z.ref, Z.K.ref=Z.K.ref, edge.col=samplk.ecol)
301
302## ----zach count latent fit, results='hide'-------------------------------
303zach.d2G2S <- ergmm(zach~nodefactor("leader")+euclidean(d=2,G=2)+rsociality,
304                    response="contexts",family="Poisson",
305                    verbose=TRUE)
306mcmc.diagnostics(zach.d2G2S)
307
308## ----zach count latent fit summ plot-------------------------------------
309summary(zach.d2G2S)
310par(mar=c(5, 4, 4, 2) + 0.1)
311plot(zach.d2G2S, rand.eff="sociality", edge.col=zach.ecol, labels=TRUE)
312