Sunbelt2015: ergmEgoTut.R

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