Sunbelt2016: Valued.R

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