Sunbelt2016: ergm.ego_tutorial.2.R

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