Sunbelt2016: ergm.ego_tutorial.R

File ergm.ego_tutorial.R, 7.2 KB (added by krivitsky, 18 months ago)
Line 
1install.packages('ergm.ego')
2
3
4library(help='ergm.ego')
5
6library(ergm.ego)
7
8help('ergm.ego-terms')
9
10
11
12
13set.seed(1)
14
15
16data("faux.mesa.high")
17mesa <- faux.mesa.high
18
19
20plot(mesa, vertex.col="Grade")
21
22
23mesa.ego <- as.egodata(mesa) # Warning because there were no vertex names in the first place.
24
25
26str(mesa.ego)
27summary(mesa.ego)
28head(mesa.ego$egos)
29head(mesa.ego$alters)
30
31
32write(t(mesa.ego$egos), file="mesa.ego.table.csv", ncol=dim(mesa.ego$egos)[2], sep=",")
33write(t(mesa.ego$alters), file="mesa.alter.table.csv", ncol=dim(mesa.ego$alters)[2], sep=",")
34
35
36mesa.egos <- read.csv("mesa.ego.table.csv")
37head(mesa.egos)
38mesa.alters <- read.csv("mesa.alter.table.csv")
39head(mesa.alters)
40
41
42test <- as.egodata(mesa.egos,alters=mesa.alters,egoIDcol="egoID")
43head(test$egos)
44head(test$alters)
45
46
47table(mesa.ego$egos$Sex, exclude=NULL)
48table(mesa.ego$egos$Race, exclude=NULL)
49plot(table(mesa.ego$egos$Grade), ylab="frequency")
50
51# compare egos and alters...
52
53par(mfrow=c(1,2))
54plot(table(mesa.ego$egos$Race, exclude=NULL)/nrow(mesa.ego$egos),
55      ylab="percent", main="Egos")
56plot(table(mesa.ego$alters$Race, exclude=NULL)/nrow(mesa.ego$alters),
57      ylab="percent", main="Alters")
58
59
60# to get the crosstabulated counts of ties:
61
62mixingmatrix.egodata(mesa.ego,"Grade")
63
64# contrast with the network's mixmatrix:
65
66mixingmatrix(mesa, "Grade")
67
68# to get the row conditional probabilities:
69
70mixingmatrix.egodata(mesa.ego, "Grade", rowprob=T)
71mixingmatrix.egodata(mesa.ego, "Race", rowprob=T)
72
73
74# ties
75nrow(mesa.ego$alters)
76
77# mean degree
78nrow(mesa.ego$egos)/nrow(mesa.ego$alters)
79
80
81## overall degree distribution
82summary(mesa.ego ~ degree(0:20))
83
84## by sex, and race
85summary(mesa.ego ~ degree(0:13, by="Sex"))
86summary(mesa.ego ~ degree(0:13, by="Race"))
87
88
89summary(mesa.ego ~ degree(0:10), scaleto=100000)
90summary(mesa.ego ~ degree(0:10), scaleto=nrow(mesa.ego$egos)*100)
91
92
93
94# to get the frequency counts
95
96degreedist.egodata(mesa.ego)
97
98degreedist.egodata(mesa.ego, by="Sex")
99
100# to get the proportion at each degree level
101
102degreedist.egodata(mesa.ego, by="Sex", prob=T)
103
104
105degreedist.egodata(mesa.ego, brg=T)
106
107degreedist.egodata(mesa.ego, by="Sex", prob=T, brg=T)
108
109
110?control.ergm.ego
111
112
113fit.edges <- ergm.ego(mesa.ego ~ edges)
114summary(fit.edges)
115
116
117names(fit.edges)
118fit.edges$ppopsize
119fit.edges$popsize
120
121
122summary(ergm.ego(mesa.ego ~ edges,
123                 control = control.ergm.ego(ppopsize=1000)))
124
125
126mcmc.diagnostics(fit.edges)
127
128
129plot(gof(fit.edges, GOF="model"))
130
131
132plot(gof(fit.edges, GOF="degree"))
133
134
135
136set.seed(0)
137fit.deg0 <- ergm.ego(mesa.ego ~ edges + degree(0), control=control.ergm.ego(ppopsize=1000))
138summary(fit.deg0)
139
140
141mcmc.diagnostics(fit.deg0)
142
143
144plot(gof(fit.deg0, GOF="model"))
145plot(gof(fit.deg0, GOF="degree"))
146
147
148fit.full <- ergm.ego(mesa.ego ~ edges + degree(0:1)
149                      + nodefactor("Sex")
150                      + nodefactor("Race", base=2)
151                      + nodefactor("Grade")
152                      + nodematch("Sex") + nodematch("Race") + nodematch("Grade"))
153summary(fit.full)
154
155
156mcmc.diagnostics(fit.full)
157
158
159plot(gof(fit.full, GOF="model"))
160plot(gof(fit.full, GOF="degree"))
161
162
163sim.full <- simulate(fit.full)
164summary(mesa.ego ~ edges + degree(0:1)
165                      + nodefactor("Sex")
166                      + nodefactor("Race", base=2)
167                      + nodefactor("Grade")
168                      + nodematch("Sex") + nodematch("Race") + nodematch("Grade"))
169
170summary(sim.full ~ edges + degree(0:1)
171                      + nodefactor("Sex")
172                      + nodefactor("Race", base=2)
173                      + nodefactor("Grade")
174                      + nodematch("Sex") + nodematch("Race") + nodematch("Grade"))
175plot(sim.full, vertex.col="Grade")
176
177
178sim.full2 <- simulate(fit.full, popsize=network.size(mesa)*2)
179summary(mesa~edges + degree(0:1)
180                      + nodefactor("Sex")
181                      + nodefactor("Race", base=2)
182                      + nodefactor("Grade")
183                      + nodematch("Sex") + nodematch("Race") + nodematch("Grade"))*2
184
185summary(sim.full2~edges + degree(0:1)
186                      + nodefactor("Sex")
187                      + nodefactor("Race", base=2)
188                      + nodefactor("Grade")
189                      + nodematch("Sex") + nodematch("Race") + nodematch("Grade"))
190
191
192
193data(faux.magnolia.high)
194faux.magnolia.high -> fmh
195(N <- network.size(fmh))
196
197
198fit.ergm <- ergm(fmh~edges+degree(0:3)+nodefactor("Race")+nodematch("Race")
199                    +nodefactor("Sex")+nodematch("Sex")+absdiff("Grade"))
200
201
202fmh.ego <- as.egodata(fmh)
203head(fmh.ego)
204
205egofit <- ergm.ego(fmh.ego~edges+degree(0:3)+nodefactor("Race")+nodematch("Race")
206                  +nodefactor("Sex")+nodematch("Sex")+absdiff("Grade"), popsize=N,
207                  ppopsize=N)
208
209modse <- function(fit) sqrt(diag(vcov(fit, sources="model"))) # A convendience function.
210
211# Parameters recovered:
212param.compare <- cbind(coef(fit.ergm), coef(egofit)[-1], (coef(fit.ergm)-coef(egofit)[-1])/modse(egofit)[-1])
213colnames(param.compare) <- c("Full nw est", "Ego census est", "diff Z")
214round(param.compare, 3)
215
216
217# MCMC diagnostics.
218mcmc.diagnostics(egofit)
219
220
221# Check whether the model converged to the right statistics:
222plot(gof(egofit, GOF="model"))
223
224
225plot(gof(egofit, GOF="degree"))
226
227
228set.seed(1)
229fmh.egosampN <- sample(fmh.ego, N, replace=TRUE)
230egofitN <- ergm.ego(fmh.egosampN ~ edges+degree(0:3)+nodefactor("Race")
231                                 +nodematch("Race")+nodefactor("Sex")+nodematch("Sex")
232                                 +absdiff("Grade"),popsize=N)
233
234# compare the coef
235param.compare <- cbind(coef(fit.ergm), coef(egofit)[-1], coef(egofitN)[-1])
236colnames(param.compare) <- c("Full nw est", "Ego census est", "Sample N est")
237round(param.compare, 3)
238                 
239# compare the s.e.'s
240se.compare <- cbind(modse(fit.ergm), modse(egofit)[-1], modse(egofitN)[-1])
241colnames(se.compare) <- c("Full nw SE", "Ego census SE", "Sample N SE")
242round(se.compare, 3)
243
244
245fmh.egosampN4 <- sample(fmh.ego, round(N/4), replace=TRUE)
246
247egofitN4 <- ergm.ego(fmh.egosampN4~edges+degree(0:3)+nodefactor("Race")
248                     +nodematch("Race")+nodefactor("Sex")+nodematch("Sex")
249                     +absdiff("Grade"), popsize=N)
250# compare the coef
251param.compare <- cbind(coef(fit.ergm), coef(egofit)[-1], coef(egofitN)[-1], coef(egofitN4)[-1])
252colnames(param.compare) <- c("Full nw est", "Ego census est", "Sample N est", "Sample N/4 est")
253round(param.compare, 3)
254                 
255# compare the s.e.'s
256se.compare <- cbind(modse(fit.ergm), modse(egofit)[-1], modse(egofitN)[-1], modse(egofitN4)[-1])
257colnames(se.compare) <- c("Full nw SE", "Ego census SE", "Sample N SE", "Sample N/4 SE")
258round(se.compare, 3)
259
260
261w <- 1+3*((fmh %v% "Race")!="White")
262fmh.egosampN4w <- sample(fmh.ego, round(N/4), replace=TRUE, prob=w)
263
264head(fmh.egosampN4w)
265egofitN4w<-ergm.ego(fmh.egosampN4w~edges+degree(0:3)+nodefactor("Race")+nodematch("Race")+nodefactor("Sex")+nodematch("Sex")+absdiff("Grade"), popsize=N)
266
267
268# compare the coef
269
270param.compare <- cbind(coef(fit.ergm), coef(egofitN4)[-1], coef(egofitN4w)[-1])
271colnames(param.compare) <- c("Full nw est", "Sample N/4 est", "Oversampled")
272round(param.compare, 3)
273                 
274# compare the s.e.'s
275
276se.compare <- cbind(modse(fit.ergm),  modse(egofitN4)[-1], modse(egofitN4w)[-1])
277colnames(se.compare) <- c("Full nw SE", "Sample N/4 SE", "Oversampled SE")
278round(se.compare, 3)
279
280