perceived.ineq <- gapply(sim_net, c(1,2), 1:N, FUN = neighborsMean)
true.ineq <- rep(mean(wage), N)
diff.ineq <- perceived.ineq - true.ineq
diff <- mean(colMeans(diff.ineq, na.rm=TRUE))
diff.low <- mean(colMeans(diff.ineq[class==1,], na.rm=TRUE))
diff.middle <- mean(colMeans(diff.ineq[class==2,], na.rm=TRUE))
diff.high <- mean(colMeans(diff.ineq[class==3,], na.rm=TRUE))
results<- rbind(results, c(low,middle,hom,diff,diff.low, diff.middle,diff.high))
}
}
}
library(statnet)
library(network)
library(ecodist)
library(ergm)
library(ineq)
library(fields)
neighborsGini <- function(arg1){
return(Gini(wage[arg1]))
}
neighborsWageVar <- function(arg1){
return(var(wage[arg1]))
}
neighborsMean <- function(args1){
return(mean(wage[arg1]))
}
N <- 500
reps<-2
results <- NA
homophily <- c(0,1,2,5)
shareMiddle <-seq(0.05,0.5,0.05)
shareLower <- seq(0.05,0.5,0.05)
net<-network(N, directed=FALSE, density=0.1)
for (hom in homophily){
for (middle in shareMiddle) {
for (low in shareLower) {
low.size <- round(low * N)
middle.size <- round(middle * N)
high.size <- (N - low.size - middle.size)
class <- c(rep(1,low.size), rep(2, middle.size), rep(3, high.size))
wage <- class
wage[class == 1]<- 100
wage[class == 2]<- 200
wage[class == 3]<- 300
net %v% "class" <- class
net %v% "wage" <- wage
# Status update
print("Homophily")
print(hom)
print("Low")
print(low)
print("Middle")
print(middle)
flush.console()
sim_net<- simulate(net ~ edges  + nodematch("class", diff=F), constraints=~degreedist, coef=c(0, hom), nsim=reps)
for (one_net in sim_net){
set.network.attribute(one_net, "loops",TRUE)
add.edges(one_net, c(seq(1,N)),c(seq(1,N)))
}
perceived.ineq <- gapply(sim_net, c(1,2), 1:N, FUN = neighborsMean)
true.ineq <- rep(mean(wage), N)
diff.ineq <- perceived.ineq - true.ineq
diff <- mean(colMeans(diff.ineq, na.rm=TRUE))
diff.low <- mean(colMeans(diff.ineq[class==1,], na.rm=TRUE))
diff.middle <- mean(colMeans(diff.ineq[class==2,], na.rm=TRUE))
diff.high <- mean(colMeans(diff.ineq[class==3,], na.rm=TRUE))
results<- rbind(results, c(low,middle,hom,diff,diff.low, diff.middle,diff.high))
}
}
}
wage
mean(wage)
neighborsMean <- function(args1){
return(mean(wage[arg1]))
}
N <- 500
reps<-2
results <- NA
homophily <- c(0,1,2,5)
shareMiddle <-seq(0.05,0.5,0.05)
shareLower <- seq(0.05,0.5,0.05)
net<-network(N, directed=FALSE, density=0.1)
for (hom in homophily){
for (middle in shareMiddle) {
for (low in shareLower) {
low.size <- round(low * N)
middle.size <- round(middle * N)
high.size <- (N - low.size - middle.size)
class <- c(rep(1,low.size), rep(2, middle.size), rep(3, high.size))
wage <- class
wage[class == 1]<- 100
wage[class == 2]<- 200
wage[class == 3]<- 300
net %v% "class" <- class
net %v% "wage" <- wage
# Status update
print("Homophily")
print(hom)
print("Low")
print(low)
print("Middle")
print(middle)
flush.console()
sim_net<- simulate(net ~ edges  + nodematch("class", diff=F), constraints=~degreedist, coef=c(0, hom), nsim=reps)
for (one_net in sim_net){
set.network.attribute(one_net, "loops",TRUE)
add.edges(one_net, c(seq(1,N)),c(seq(1,N)))
}
perceived.ineq <- gapply(sim_net, c(1,2), 1:N, FUN = neighborsMean)
true.ineq <- rep(mean(wage), N)
diff.ineq <- perceived.ineq - true.ineq
diff <- mean(colMeans(diff.ineq, na.rm=TRUE))
diff.low <- mean(colMeans(diff.ineq[class==1,], na.rm=TRUE))
diff.middle <- mean(colMeans(diff.ineq[class==2,], na.rm=TRUE))
diff.high <- mean(colMeans(diff.ineq[class==3,], na.rm=TRUE))
results<- rbind(results, c(low,middle,hom,diff,diff.low, diff.middle,diff.high))
}
}
}
for (hom in homophily){
for (middle in shareMiddle) {
for (low in shareLower) {
low.size <- round(low * N)
middle.size <- round(middle * N)
high.size <- (N - low.size - middle.size)
class <- c(rep(1,low.size), rep(2, middle.size), rep(3, high.size))
wage <- class
wage[class == 1]<- 100
wage[class == 2]<- 200
wage[class == 3]<- 300
net %v% "class" <- class
net %v% "wage" <- wage
# Status update
print("Homophily")
print(hom)
print("Low")
print(low)
print("Middle")
print(middle)
sim_net<- simulate(net ~ edges  + nodematch("class", diff=F), constraints=~degreedist, coef=c(0, hom), nsim=reps)
for (one_net in sim_net){
set.network.attribute(one_net, "loops",TRUE)
add.edges(one_net, c(seq(1,N)),c(seq(1,N)))
}
flush.console()
perceived.ineq <- gapply(sim_net, c(1,2), 1:N, FUN = neighborsMean)
true.ineq <- rep(mean(wage), N)
diff.ineq <- perceived.ineq - true.ineq
diff <- mean(colMeans(diff.ineq, na.rm=TRUE))
diff.low <- mean(colMeans(diff.ineq[class==1,], na.rm=TRUE))
diff.middle <- mean(colMeans(diff.ineq[class==2,], na.rm=TRUE))
diff.high <- mean(colMeans(diff.ineq[class==3,], na.rm=TRUE))
results<- rbind(results, c(low,middle,hom,diff,diff.low, diff.middle,diff.high))
}
}
}
for (hom in homophily){
for (middle in shareMiddle) {
for (low in shareLower) {
low.size <- round(low * N)
middle.size <- round(middle * N)
high.size <- (N - low.size - middle.size)
class <- c(rep(1,low.size), rep(2, middle.size), rep(3, high.size))
wage <- class
wage[class == 1]<- 100
wage[class == 2]<- 200
wage[class == 3]<- 300
net %v% "class" <- class
net %v% "wage" <- wage
# Status update
print("Homophily")
print(hom)
print("Low")
print(low)
print("Middle")
print(middle)
flush.console()
sim_net<- simulate(net ~ edges  + nodematch("class", diff=F), constraints=~degreedist, coef=c(0, hom), nsim=reps)
for (one_net in sim_net){
set.network.attribute(one_net, "loops",TRUE)
add.edges(one_net, c(seq(1,N)),c(seq(1,N)))
}
perceived.ineq <- gapply(sim_net, c(1,2), 1:N, FUN = neighborsWageVar)
true.ineq <- rep(var(wage), N)
diff.ineq <- perceived.ineq - true.ineq
diff <- mean(colMeans(diff.ineq, na.rm=TRUE))
diff.low <- mean(colMeans(diff.ineq[class==1,], na.rm=TRUE))
diff.middle <- mean(colMeans(diff.ineq[class==2,], na.rm=TRUE))
diff.high <- mean(colMeans(diff.ineq[class==3,], na.rm=TRUE))
results<- rbind(results, c(low,middle,hom,diff,diff.low, diff.middle,diff.high))
}
}
}
neighborsMean <- function(arg1){
return(mean(wage[arg1]))
}
N <- 500
reps<-2
results <- NA
homophily <- c(0,1,2,5)
shareMiddle <-seq(0.05,0.5,0.05)
shareLower <- seq(0.05,0.5,0.05)
net<-network(N, directed=FALSE, density=0.1)
for (hom in homophily){
for (middle in shareMiddle) {
for (low in shareLower) {
low.size <- round(low * N)
middle.size <- round(middle * N)
high.size <- (N - low.size - middle.size)
class <- c(rep(1,low.size), rep(2, middle.size), rep(3, high.size))
wage <- class
wage[class == 1]<- 100
wage[class == 2]<- 200
wage[class == 3]<- 300
net %v% "class" <- class
net %v% "wage" <- wage
# Status update
print("Homophily")
print(hom)
print("Low")
print(low)
print("Middle")
print(middle)
flush.console()
sim_net<- simulate(net ~ edges  + nodematch("class", diff=F), constraints=~degreedist, coef=c(0, hom), nsim=reps)
for (one_net in sim_net){
set.network.attribute(one_net, "loops",TRUE)
add.edges(one_net, c(seq(1,N)),c(seq(1,N)))
}
perceived.ineq <- gapply(sim_net, c(1,2), 1:N, FUN = neighborsMean)
true.ineq <- rep(mean(wage), N)
diff.ineq <- perceived.ineq - true.ineq
diff <- mean(colMeans(diff.ineq, na.rm=TRUE))
diff.low <- mean(colMeans(diff.ineq[class==1,], na.rm=TRUE))
diff.middle <- mean(colMeans(diff.ineq[class==2,], na.rm=TRUE))
diff.high <- mean(colMeans(diff.ineq[class==3,], na.rm=TRUE))
results<- rbind(results, c(low,middle,hom,diff,diff.low, diff.middle,diff.high))
}
}
}
N <- 100
reps<-2
results <- NA
homophily <- c(0,1,2,5)
shareMiddle <-seq(0.05,0.5,0.05)
shareLower <- seq(0.05,0.5,0.05)
net<-network(N, directed=FALSE, density=0.1)
for (hom in homophily){
for (middle in shareMiddle) {
for (low in shareLower) {
low.size <- round(low * N)
middle.size <- round(middle * N)
high.size <- (N - low.size - middle.size)
class <- c(rep(1,low.size), rep(2, middle.size), rep(3, high.size))
wage <- class
wage[class == 1]<- 100
wage[class == 2]<- 200
wage[class == 3]<- 300
net %v% "class" <- class
net %v% "wage" <- wage
# Status update
print("Homophily")
print(hom)
print("Low")
print(low)
print("Middle")
print(middle)
flush.console()
sim_net<- simulate(net ~ edges  + nodematch("class", diff=F), constraints=~degreedist, coef=c(0, hom), nsim=reps)
for (one_net in sim_net){
set.network.attribute(one_net, "loops",TRUE)
add.edges(one_net, c(seq(1,N)),c(seq(1,N)))
}
perceived.ineq <- gapply(sim_net, c(1,2), 1:N, FUN = neighborsMean)
true.ineq <- rep(mean(wage), N)
diff.ineq <- perceived.ineq - true.ineq
diff <- mean(colMeans(diff.ineq, na.rm=TRUE))
diff.low <- mean(colMeans(diff.ineq[class==1,], na.rm=TRUE))
diff.middle <- mean(colMeans(diff.ineq[class==2,], na.rm=TRUE))
diff.high <- mean(colMeans(diff.ineq[class==3,], na.rm=TRUE))
results<- rbind(results, c(low,middle,hom,diff,diff.low, diff.middle,diff.high))
}
}
}
results<-results[-1,]
topo<- as.data.frame(results)
names(topo)[names(topo)=="V1"] <- "shareLower"
names(topo)[names(topo)=="V2"] <- "shareMiddle"
names(topo)[names(topo)=="V3"] <- "homophily"
names(topo)[names(topo)=="V4"] <- "diff"
names(topo)[names(topo)=="V5"] <- "diffLow"
names(topo)[names(topo)=="V6"] <- "diffMiddle"
names(topo)[names(topo)=="V7"] <- "diffHigh"
topoBackup <- topo
shareLower <- seq(min(shareLower), max(shareLower), .005)
shareMiddle <- seq (min(shareMiddle), max(shareMiddle), .005)
topo.loess <- loess (diff ~  shareLower * shareMiddle, topo, degree = 2, span = 0.2)
interpolated <- predict(topo.loess, expand.grid(shareLower=shareLower,shareMiddle=shareMiddle))
range(interpolated)
shareLower <- seq(min(shareLower), max(shareLower), .005)
shareMiddle <- seq (min(shareMiddle), max(shareMiddle), .005)
topo.loess <- loess (diff ~  shareLower * shareMiddle, topo, degree = 2, span = 0.2)
interpolated <- predict(topo.loess, expand.grid(shareLower=shareLower,shareMiddle=shareMiddle))
range(interpolated)
par(mfrow=c(2,3))
par(mar=c(3,3,3,3))
for (hom in homophily){
topo<- topoBackup[topoBackup$homophily==hom,]
topo.loess <- loess (diff ~  shareLower * shareMiddle, topo, degree = 2, span = 0.2)
interpolated <- predict(topo.loess, expand.grid(shareLower=shareLower,shareMiddle=shareMiddle))
title <- paste("all, homophily =", hom)
image.plot(x=shareLower, y=shareMiddle, z = interpolated, xlim=range(0.05,.5), ylim=range(0.05,.5), zlim = range(-5,5), main=title, legend.only=FALSE)
topo.loess <- loess (diffLow ~ shareLower * shareMiddle, topo, degree = 2, span = 0.2)
interpolated <- predict(topo.loess, expand.grid(shareLower=shareLower,shareMiddle=shareMiddle))
title <- paste("low, homophily =", hom)
image.plot(x=shareLower, y=shareMiddle, z = interpolated, xlim=range(0.05,.5), ylim=range(0.05,.5), zlim = range(-5,5), main=title, legend.only=FALSE)
plot.new()
topo.loess <- loess (diffMiddle ~ shareLower * shareMiddle, topo, degree = 2, span = 0.2)
interpolated <- predict(topo.loess, expand.grid(shareLower=shareLower,shareMiddle=shareMiddle))
title <- paste("middle, homophily =", hom)
image.plot(x=shareLower, y=shareMiddle, z = interpolated, xlim=range(0.05,.5), ylim=range(0.05,.5), zlim = range(-5,5), main=title, legend.only=FALSE)
topo.loess <- loess (diffHigh ~ shareLower * shareMiddle, topo, degree = 2, span = 0.2)
interpolated <- predict(topo.loess, expand.grid(shareLower=shareLower,shareMiddle=shareMiddle))
title <- paste("high, homophily =", hom)
image.plot(x=shareLower, y=shareMiddle, z = interpolated, xlim=range(0.05,.5), ylim=range(0.05,.5), zlim = range(-5,5), main=title, legend.only=FALSE)
plot.new()
}
topo<- topoBackup[topoBackup$homophily==hom,]
topo.loess <- loess (diff ~  shareLower * shareMiddle, topo, degree = 2, span = 0.2)
interpolated <- predict(topo.loess, expand.grid(shareLower=shareLower,shareMiddle=shareMiddle))
range(interpolated)
topo.loess <- loess (diffLow ~ shareLower * shareMiddle, topo, degree = 2, span = 0.2)
interpolated <- predict(topo.loess, expand.grid(shareLower=shareLower,shareMiddle=shareMiddle))
range(interpolated)
topo.loess <- loess (diffMiddle ~ shareLower * shareMiddle, topo, degree = 2, span = 0.2)
interpolated <- predict(topo.loess, expand.grid(shareLower=shareLower,shareMiddle=shareMiddle))
range(interpolated)
title <- paste("middle, homophily =", hom)
image.plot(x=shareLower, y=shareMiddle, z = interpolated, xlim=range(0.05,.5), ylim=range(0.05,.5), zlim = range(-5,5), main=title, legend.only=FALSE)
topo.loess <- loess (diffHigh ~ shareLower * shareMiddle, topo, degree = 2, span = 0.2)
interpolated <- predict(topo.loess, expand.grid(shareLower=shareLower,shareMiddle=shareMiddle))
range(interpolated)
title <- paste("high, homophily =", hom)
image.plot(x=shareLower, y=shareMiddle, z = interpolated, xlim=range(0.05,.5), ylim=range(0.05,.5), zlim = range(-5,5), main=title, legend.only=FALSE)
plot.new()
par(mfrow=c(2,3))
par(mar=c(3,3,3,3))
for (hom in homophily){
topo<- topoBackup[topoBackup$homophily==hom,]
topo.loess <- loess (diff ~  shareLower * shareMiddle, topo, degree = 2, span = 0.2)
interpolated <- predict(topo.loess, expand.grid(shareLower=shareLower,shareMiddle=shareMiddle))
range(interpolated)
title <- paste("all, homophily =", hom)
image.plot(x=shareLower, y=shareMiddle, z = interpolated, xlim=range(0.05,.5), ylim=range(0.05,.5), zlim = range(-100,100), main=title, legend.only=FALSE)
topo.loess <- loess (diffLow ~ shareLower * shareMiddle, topo, degree = 2, span = 0.2)
interpolated <- predict(topo.loess, expand.grid(shareLower=shareLower,shareMiddle=shareMiddle))
range(interpolated)
title <- paste("low, homophily =", hom)
image.plot(x=shareLower, y=shareMiddle, z = interpolated, xlim=range(0.05,.5), ylim=range(0.05,.5), zlim = range(-100,100), main=title, legend.only=FALSE)
plot.new()
topo.loess <- loess (diffMiddle ~ shareLower * shareMiddle, topo, degree = 2, span = 0.2)
interpolated <- predict(topo.loess, expand.grid(shareLower=shareLower,shareMiddle=shareMiddle))
range(interpolated)
title <- paste("middle, homophily =", hom)
image.plot(x=shareLower, y=shareMiddle, z = interpolated, xlim=range(0.05,.5), ylim=range(0.05,.5), zlim = range(-100,100), main=title, legend.only=FALSE)
topo.loess <- loess (diffHigh ~ shareLower * shareMiddle, topo, degree = 2, span = 0.2)
interpolated <- predict(topo.loess, expand.grid(shareLower=shareLower,shareMiddle=shareMiddle))
range(interpolated)
title <- paste("high, homophily =", hom)
image.plot(x=shareLower, y=shareMiddle, z = interpolated, xlim=range(0.05,.5), ylim=range(0.05,.5), zlim = range(-100,100), main=title, legend.only=FALSE)
plot.new()
}
library(statnet)
library(foreign)
data(florentine)
flobusiness
flomarriage
help (package = sna)
betweenness(flomarriage)
plot(flomarriage)
help plot.network
help network.plot
? plot.network
flomarriage
flomarriage %v% vertex.names
flomarriage %v% wealth
flomarriage %v% "wealth"
flomarriage %v% "vertex.names"
plot(flomarriage, label = flomarriage %v% "vertex.names")
betweeness(flobusiness)
betweenness(flobusiness)
betweenness(flomarriage)
?betweenness
betweenness(flomarriage, gmode="graph")
betweenness(flomarriage)
plot(flomarriage, label = flomarriage %v% "vertex.names")
help (package=network)
delete.vertices(flomarriage, 12)
betweenness(flomarriage)
plot(flomarriage)
install.packages("NetCluster")
edit(hclust)
edit(.Fortran)
edit(hclust)
?hclust
?cor
edit(equiv.dist)
libary(sna)
library(sna)
edit(equiv.dist)
?sna
# install and load necessary R packages
if (!require("foreign",character.only = TRUE)) { install.packages("foreign",  repos = "http://cran.rstudio.com/") }
if (!require("network",character.only = TRUE)) { install.packages("network", repos = "http://cran.rstudio.com/") }
if (!require("ergm",character.only = TRUE)) { install.packages("ergm", repos = "http://cran.rstudio.com/") }
suppressPackageStartupMessages(library("foreign", quietly=TRUE, verbose=FALSE, warn.conflicts=FALSE))
suppressPackageStartupMessages(library("ergm", quietly=TRUE, verbose=FALSE, warn.conflicts=FALSE))
suppressPackageStartupMessages(library("network", quietly=TRUE, verbose=FALSE, warn.conflicts=FALSE))
# install and load necessary R packages
if (!require("foreign",character.only = TRUE)) { install.packages("foreign",  repos = "http://cran.rstudio.com/") }
if (!require("statnet",character.only = TRUE)) { install.packages("statnet", repos = "http://cran.rstudio.com/") }
suppressPackageStartupMessages(library("foreign", quietly=TRUE, verbose=FALSE, warn.conflicts=FALSE))
suppressPackageStartupMessages(library("statnet", quietly=TRUE, verbose=FALSE, warn.conflicts=FALSE))
# set R working directory to Stata working directory
setwd("C:/Users/Thomas/Dropbox/STATA/nwcommands")
# load network and attributes
data <- read.dta("ergdata.dta")
netsize <- dim(data)[1]
netmat <- as.matrix(data[unlist(strsplit("net1 net2 net3 net4 net5 net6 net7 net8 net9 net10 net11 net12 net13 net14 net15 net16 net17 net18 net19 net20"," "))])
netsym <- all(netmat == t(netmat))
net<- network(netmat, directed = !netsym)
attrs <- dim(data)[2] - dim(data)[1]
for (i in 1:attrs){
att <- data[,i]
set.vertex.attribute(net, colnames(data[i]),att)
}
# run ERGM
model <- net~edges
try({summary(model)
ergmresults<-ergm(model, , verbose=TRUE)
summary(ergmresults)})
write.table(geterrmessage(),'ergerror.txt')
# produce output files: ergcoefs.csv, ergcov.csv, ergmodel.csv,  ergstats.csv
write.csv(round(summary(ergmresults)$coefs, 3), 'ergcoefs.csv', na='.')
write.csv(summary(ergmresults)$asycov, 'ergcov.csv', na='.')
write.csv(round(as.data.frame(summary(model)), 3), 'ergmodel.csv', na='.')
capture.output(summary(ergmresults)$control, file='ergcontrol.csv')
fileCon<-file('ergstats.csv', 'w')
writeLines('vertices,edges,directed,numcoeff,coeff,iterations,estimate,aic,bic,samplesize,message', fileCon)
writeLines(paste(network.size(net),',',network.edgecount(net),',',net[[2]]$directed,',',dim(summary(ergmresults)$coef)[1],',',gsub(',','',toString(names(summary(model)))),',',summary(ergmresults)$iterations, ',',summary(ergmresults)$estimate,',',round(summary(ergmresults)$aic,2),',',round(summary(ergmresults)$bic,2),',',summary(ergmresults)$samplesize,',',summary(ergmresults)$message), fileCon)
close(fileCon)
# install and load necessary R packages
if (!require("foreign",character.only = TRUE)) { install.packages("foreign",  repos = "http://cran.rstudio.com/") }
if (!require("statnet",character.only = TRUE)) { install.packages("statnet", repos = "http://cran.rstudio.com/") }
suppressPackageStartupMessages(library("foreign", quietly=TRUE, verbose=FALSE, warn.conflicts=FALSE))
suppressPackageStartupMessages(library("statnet", quietly=TRUE, verbose=FALSE, warn.conflicts=FALSE))
# set R working directory to Stata working directory
setwd("C:/Users/Thomas/Dropbox/STATA/nwcommands")
# load network and attributes
data <- read.dta("ergdata.dta")
netsize <- dim(data)[1]
netmat <- as.matrix(data[unlist(strsplit("net1 net2 net3 net4 net5 net6 net7 net8 net9 net10 net11 net12 net13 net14 net15 net16 net17 net18 net19 net20"," "))])
netsym <- all(netmat == t(netmat))
net<- network(netmat, directed = !netsym)
attrs <- dim(data)[2] - dim(data)[1]
for (i in 1:attrs){
att <- data[,i]
set.vertex.attribute(net, colnames(data[i]),att)
}
# run ERGM
model <- net~edges
try({summary(model)
ergmresults<-ergm(model, , verbose=TRUE)
summary(ergmresults)})
write.table(geterrmessage(),'ergerror.txt')
# produce output files: ergcoefs.csv, ergcov.csv, ergmodel.csv,  ergstats.csv
write.csv(round(summary(ergmresults)$coefs, 3), 'ergcoefs.csv', na='.')
write.csv(summary(ergmresults)$asycov, 'ergcov.csv', na='.')
write.csv(round(as.data.frame(summary(model)), 3), 'ergmodel.csv', na='.')
# install and load necessary R packages
if (!require("foreign",character.only = TRUE)) { install.packages("foreign",  repos = "http://cran.rstudio.com/") }
if (!require("statnet",character.only = TRUE)) { install.packages("statnet", repos = "http://cran.rstudio.com/") }
suppressPackageStartupMessages(library("foreign", quietly=TRUE, verbose=FALSE, warn.conflicts=FALSE))
suppressPackageStartupMessages(library("statnet", quietly=TRUE, verbose=FALSE, warn.conflicts=FALSE))
# set R working directory to Stata working directory
setwd("C:/Users/Thomas/Dropbox/STATA/nwcommands")
# load network and attributes
data <- read.dta("ergdata.dta")
netsize <- dim(data)[1]
netmat <- as.matrix(data[unlist(strsplit("net1 net2 net3 net4 net5 net6 net7 net8 net9 net10 net11 net12 net13 net14 net15 net16 net17 net18 net19 net20"," "))])
netsym <- all(netmat == t(netmat))
net<- network(netmat, directed = !netsym)
attrs <- dim(data)[2] - dim(data)[1]
for (i in 1:attrs){
att <- data[,i]
set.vertex.attribute(net, colnames(data[i]),att)
}
# run ERGM
model <- net~edges
try({summary(model)
ergmresults<-ergm(model, , verbose=TRUE)
summary(ergmresults)})
write.table(geterrmessage(),'ergerror.txt')
# produce output files: ergcoefs.csv, ergcov.csv, ergmodel.csv,  ergstats.csv
write.csv(round(summary(ergmresults)$coefs, 3), 'ergcoefs.csv', na='.')
write.csv(summary(ergmresults)$asycov, 'ergcov.csv', na='.')
write.csv(round(as.data.frame(summary(model)), 3), 'ergmodel.csv', na='.')
capture.output(summary(ergmresults)$control, file='ergcontrol.csv')
fileCon<-file('ergstats.csv', 'w')
writeLines('vertices,edges,directed,numcoeff,coeff,iterations,estimate,aic,bic,samplesize,message', fileCon)
writeLines(paste(network.size(net),',',network.edgecount(net),',',net[[2]]$directed,',',dim(summary(ergmresults)$coef)[1],',',gsub(',','',toString(names(summary(model)))),',',summary(ergmresults)$iterations, ',',summary(ergmresults)$estimate,',',round(summary(ergmresults)$aic,2),',',round(summary(ergmresults)$bic,2),',',summary(ergmresults)$samplesize,',',summary(ergmresults)$message), fileCon)
close(fileCon)