# Replication code for Leifeld and Schneider (2012), AJPS.
# Note that the estimates can only be reproduced approximately
# due to internal changes in the statnet package.
# preparatory steps
library("statnet")
library("xergm")
library("texreg")
seed <- 12345
set.seed(seed)
data("chemnet")
# create confirmed network relation
sci <- scito * t(scifrom) # equation 1 in the AJPS paper
prefsim <- dist(intpos, method = "euclidean") # equation 2
prefsim <- max(prefsim) - prefsim # equation 3
prefsim <- as.matrix(prefsim)
committee <- committeediag(committee) <- 0 # the diagonal has no meaning
types <- types[, 1] # convert to vector
# create network objects and store attributes
nw.pol <- network(pol) # political/stratgic information exchange
set.vertex.attribute(nw.pol, "orgtype", types)
set.vertex.attribute(nw.pol, "betweenness",
betweenness(nw.pol)) # centrality
nw.sci <- network(sci) # technical/scientific information exchange
set.vertex.attribute(nw.sci, "orgtype", types)
set.vertex.attribute(nw.sci, "betweenness",
betweenness(nw.sci)) # centrality
# ERGM: model 1 in the AJPS paper; only preference similarity
model1 <- ergm(nw.pol ~ edges + edgecov(prefsim),
control = control.ergm(seed = seed))
summary(model1)
# ERGM: model 2 in the AJPS paper; complete model
model2 <- ergm(nw.pol ~
edges +
edgecov(prefsim) +
mutual +
nodemix("orgtype", base = -7) +
nodeifactor("orgtype", base = -1) +
nodeofactor("orgtype", base = -5) +
edgecov(committee) +
edgecov(nw.sci) +
edgecov(infrep) +
gwesp(0.1, fixed = TRUE) +
gwdsp(0.1, fixed = TRUE),
control = control.ergm(seed = seed)
)
summary(model2)
# ERGM: model 3 in the AJPS paper; only preference similarity
model3 <- ergm(nw.sci ~ edges + edgecov(prefsim),
control = control.ergm(seed = seed))
summary(model3)
# ERGM: model 4 in the AJPS paper; complete model
model4 <- ergm(nw.sci ~
edges +
edgecov(prefsim) +
mutual +
nodemix("orgtype", base = -7) +
nodeifactor("orgtype", base = -1) +
nodeofactor("orgtype", base = -5) +
edgecov(committee) +
edgecov(nw.pol) +
edgecov(infrep) +
gwesp(0.1, fixed = TRUE) +
gwdsp(0.1, fixed = TRUE),
control = control.ergm(seed = seed)
)
summary(model4)
# regression table using the texreg package
screenreg(list(model1, model2, model3, model4))
# goodness of fit using the xergm package
gof2 <- gof(model2, roc = FALSE, pr = FALSE)
gof2 # print gof output
plot(gof2) # visual inspection of GOF
gof4 <- gof(model4, roc = FALSE, pr = FALSE)
gof4
plot(gof4)
# MCMC diagnostics
pdf("diagnostics2.pdf")
mcmc.diagnostics(model2)
dev.off()
pdf("diagnostics4.pdf")
mcmc.diagnostics(model4)
dev.off()
Run the code above in your browser using DataLab