# \dontshow{
require(ergm)
require(network)
require(sna)
set.seed(21093)
a1<-network::as.network(matrix(c(rbinom(10, 1,.3),
rbinom(10, 1,.3),
rbinom(10, 1,.3),
rbinom(10, 1,.3),
rbinom(10, 1,.3),
rbinom(10, 1,.3),
rbinom(10, 1,.3),
rbinom(10, 1,.3),
rbinom(10, 1,.3),
rbinom(10, 1,.3)),
nrow=10,ncol=10))
network::set.vertex.attribute(a1,"var.1",rbinom(10,1,.3))
a<-ergm(a1~edges+nodeifactor("var.1")+nodeofactor("var.1"))
MEMS(a,
micro_process="nodeifactor.var.1.1",
macro_function=gtrans,
nsim=20,
silent=TRUE)
# }
# \donttest{
########################################
# ERGM examples and basic utilities
#######################################
####start with a simple model
library(statnet)
data("faux.mesa.high")
model1<-ergm(faux.mesa.high~edges+
nodecov("Grade")+
nodefactor("Race")+
nodefactor("Sex")+
nodematch("Race")+
nodematch("Sex")+
absdiff("Grade"))
##calculate the MEMS when the absolute difference in grade is changed from an interval of 0 to 1
#with default specifications for gtrans
MEMS(model1,
micro_process="absdiff.Grade",
macro_function = gtrans,
object_type = "network",
nsim=100,
interval=c(0,1),
silent=FALSE,
algorithm = "parametric")
#call an argument from gtrans by specifying it as a function
#use nonparametric estimation
MEMS(model1,
micro_process="absdiff.Grade",
macro_function = function(x){gtrans(x,measure="strongcensus")},
object_type = "network",
nsim=100,
interval=c(0,1),
silent=FALSE,
algorithm = "nonparametric")
####calculate the MEMS using igraph
MEMS(model1,
micro_process="absdiff.Grade",
macro_function = function(x){igraph::transitivity(x,type="local")},
object_type = "igraph",
nsim=100,
interval=c(0,1),
silent=FALSE,
algorithm = "parametric")
##specify a user function that counts the number of communities
community_counts<-function(x){
walktrap<-igraph::walktrap.community(x) #use walktrap community detection
return(length(unique(walktrap$membership))) #return the number of communities
}
MEMS(model1,
micro_process="absdiff.Grade",
macro_function = community_counts,
object_type = "igraph",
nsim=100,
interval=c(0,1),
silent=FALSE,
algorithm = "parametric")
##calculate a function using exogenous node attributes
assortativity_grade<-function(x){
require(igraph)
return(assortativity_nominal(x,V(x)$Grade))
}
MEMS(model1,
micro_process="absdiff.Grade",
macro_function = assortativity_grade,
object_type = "igraph",
nsim=100,
interval=c(0,1),
silent=FALSE,
algorithm = "parametric")
##specify a user function that does not depend on either igraph or statnet
#assuming a network input object, we have
manual_user_function<-function(x){
x<-as.sociomatrix(x)
return(colSums(x))
}
MEMS(model1,
micro_process="absdiff.Grade",
macro_function = manual_user_function,
object_type = "network",
nsim=100,
interval=c(0,1),
silent=FALSE,
algorithm = "parametric")
####estimation for POOLED ERGM
data("faux.magnolia.high")
model2<-ergm(faux.magnolia.high~edges+
nodecov("Grade")+
nodefactor("Race")+
nodefactor("Sex")+
nodematch("Race")+
nodematch("Sex")+
absdiff("Grade"))
MEMS(list(model1,model2),
micro_process="absdiff.Grade",
macro_function = assortativity_grade,
object_type = "igraph",
nsim=50,
interval=c(0,1),
silent=FALSE,
algorithm = "parametric")
#################################
# Estimation with GLM and GLMER
#################################
library(btergm)
#use models 1 and 2 from examples above
glm_dat<-edgeprob(model1)
glm_dat2<-edgeprob(model2)
glm_dat2<-glm_dat2[,-c(4)]
##create stacked dataset for the purposes of grouped estimation
glm_dat$net_id<-"mesa" #specify ID for each network
glm_dat2$net_id<-"magnolia"
glm_dat<-rbind(glm_dat,glm_dat2)
##estimate as a linear probability model
net_glm<-glm(tie~nodecov.Grade+
nodefactor.Race.Hisp+
nodefactor.Race.NatAm+
nodefactor.Race.Other+
nodefactor.Sex.M+
nodematch.Race+
nodematch.Sex+
absdiff.Grade,
data=glm_dat)
MEMS(net_glm,
micro_process="nodematch.Race", #should be written as in netlogit output
macro_function = function(x){gtrans(x)},
object_type = "network",
nsim=100,
interval=c(0,.5),
silent=FALSE,
full_output = FALSE,
algorithm = "parametric",
group_id=glm_dat$net_id, #provide network ID for estimation
node_numbers =c(network.size(faux.mesa.high), #provide the number of nodes in each network
network.size(faux.magnolia.high)))
##estimate as a multilevel model
library(lme4)
net_glmer<-glmer(tie~nodecov.Grade+
nodefactor.Race.Hisp+
nodefactor.Race.NatAm+
nodefactor.Race.Other+
nodefactor.Sex.M+
nodematch.Race+
nodematch.Sex+
absdiff.Grade+
(1|net_id),
data=glm_dat,
family=gaussian)
MEMS(net_glmer,
micro_process="nodematch.Race", #should be written as in netlogit output
macro_function = function(x){gtrans(x)},
object_type = "network",
nsim=50,
interval=c(0,.5),
silent=FALSE,
full_output = FALSE,
algorithm = "parametric",
group_id=glm_dat$net_id,
node_numbers =c(203,974))
##############################################
##nonparametric estimation for bootstrap TERGM
##############################################
library(btergm)
data(alliances)
ally_data<-list(LSP[[1]],
LSP[[2]],
LSP[[3]])
bt_model<-btergm(ally_data~edges+
gwesp(.7,fixed=T)+
mutual,R=200)
MEMS(bt_model,
micro_process="gwesp.OTP.fixed.0.7",
macro_function = gtrans,
object_type = "network",
nsim=50,
interval=c(0,1),
silent=FALSE,
algorithm = "nonparametric")
################################
# Parametric estimation using SAOM
##################################
library(RSiena)
#specify 3 wave network panel data as DV
network_list<-array(c(s501,s502,s503),dim = c(50,50,3))
Network<-sienaDependent(network_list)
Smoking<-varCovar(s50s)
Alcohol<-varCovar(s50a)
SAOM.Data<-sienaDataCreate(Network=Network,Smoking,Alcohol)
#specify
SAOM.terms<-getEffects(SAOM.Data)
SAOM.terms<-includeEffects(SAOM.terms,egoX,altX,sameX,interaction1="Alcohol")
SAOM.terms<-includeEffects(SAOM.terms,egoX,altX,sameX,interaction1="Smoking")
SAOM.terms<-includeEffects(SAOM.terms,transTies,inPop)
create.model<-sienaAlgorithmCreate(projname="netmediate",
nsub=5,
n3=2000)
##estimate the model using siena07
SAOM_model<-siena07(create.model,
data=SAOM.Data,
effects=SAOM.terms,
verbose=TRUE)
SAOM_model
##basic specification for reciprocity effects on outdegree distribution
MEMS(SAOM_model,
micro_process="reciprocity", #should be written as in SIENA output
macro_function = function(x){igraph::degree(x,mode="out")},
object_type = "igraph",
interval=c(0,.5),
SAOM_data=SAOM.Data,
silent=FALSE,
algorithm = "parametric")
##include user functions on time varying covariates
assortativity_smoking<-function(x){
return(assortativity_nominal(x,V(x)$Smoking))
}
MEMS(SAOM_model,
micro_process="reciprocity",
macro_function = assortativity_smoking,
object_type = "igraph",
interval=c(0,.5),
SAOM_data=SAOM.Data,
SAOM_var=list(Smoking=Smoking,Alcohol=Alcohol), #Smoking and Alcohol are varCovar objects
silent=FALSE,
full_output = FALSE,
algorithm = "parametric")
###Pooled SAOM
MEMS(list(SAOM_model,SAOM_model),
micro_process="reciprocity",
macro_function = gtrans,
object_type = "network",
interval=c(0,.5),
SAOM_data=list(SAOM.Data,SAOM.Data),
silent=FALSE,
full_output = FALSE,
nsim=100,
algorithm = "parametric")
#Pooled SAOM with user functions and time varying attributes
assortativity_smoking<-function(x){
return(assortativity_nominal(x,V(x)$Smoking))
}
MEMS(list(SAOM_model,SAOM_model),
micro_process="reciprocity",
macro_function = assortativity_smoking,
object_type = "igraph",
interval=c(0,.5),
SAOM_data=list(SAOM.Data,SAOM.Data),
SAOM_var=list(list(Smoking=Smoking,Alcohol=Alcohol),
list(Smoking=Smoking,Alcohol=Alcohol)),
silent=FALSE,
full_output = FALSE,
nsim=100,
algorithm = "parametric")
################################################
## Selection and Influence in SAOM when analyzing
## co-evolution of networks and behavior
################################################
##Example Moran decomposition
library(RSiena)
###run the model--taken from RSiena scripts
# prepare first two waves of s50 data for RSiena analysis:
(thedata <- sienaDataCreate(
friendship = sienaDependent(array(
c(s501,s502),dim=c(50,50,2))),
drinking = sienaDependent(s50a[,1:2])
))
# specify a model with (generalised) selection and influence:
themodel <- getEffects(thedata)
themodel <- includeEffects(themodel,name='friendship',gwespFF)
themodel <- includeEffects(themodel,name='friendship',simX,interaction1='drinking')
themodel <- includeEffects(themodel,name='drinking',avSim,interaction1='friendship')
themodel
# estimate this model:
estimation.options <- sienaAlgorithmCreate(projname='results',cond=FALSE,seed=1234567)
(theresults <- siena07(estimation.options,data=thedata,effects=themodel))
##calculate MEMS for selection effect
#Uses Moran_dv--a function internally called by netmediate
#to calculate change in amount of network autocorrelation
#as a function of both endogenous behavior and network dependent
#variables
MEMS(theresults,
micro_process="drinking similarity",
macro_function =Moran_dv,
object_type = "network",
SAOM_data = thedata,
silent=FALSE,
nsim=50)
#just influence
MEMS(theresults,
micro_process="drinking average similarity",
macro_function =Moran_dv,
object_type = "network",
SAOM_data = thedata,
silent=FALSE,
nsim=50)
##joint effect of selection and influence
MEMS(theresults,
micro_process=c("drinking similarity","drinking average similarity"),
macro_function =Moran_dv,
object_type = "network",
SAOM_data = thedata,
silent=FALSE,
nsim=500)
#######################################
# Relational event models using relevent
#######################################
set.seed(21093)
library(relevent)
##generate a network with 15 discrete time periods
#example based on relevent rem.dyad example
library(relevent)
roweff<-rnorm(10) #Build rate matrix
roweff<-roweff-roweff[1] #Adjust for later convenience
coleff<-rnorm(10)
coleff<-coleff-coleff[1]
lambda<-exp(outer(roweff,coleff,"+"))
diag(lambda)<-0
ratesum<-sum(lambda)
esnd<-as.vector(row(lambda)) #List of senders/receivers
erec<-as.vector(col(lambda))
time<-0
edgelist<-vector()
while(time<15){ # Observe the system for 15 time units
drawsr<-sample(1:100,1,prob=as.vector(lambda)) #Draw from model
time<-time+rexp(1,ratesum)
if(time<=15) #Censor at 15
edgelist<-rbind(edgelist,c(time,esnd[drawsr],erec[drawsr]))
else
edgelist<-rbind(edgelist,c(15,NA,NA))
}
effects<-c("CovSnd","FERec")
##estimate model
fit.time<-rem.dyad(edgelist,10,effects=effects,
covar=list(CovSnd=roweff),
ordinal=FALSE,hessian=TRUE)
###aggregate estimation
MEMS(fit.time,
micro_process="CovSnd.1", #should be written as in relevent output
macro_function = function(x){sna::degree(x)},
object_type = "network",
nsim=10,
interval=c(0,.5),
silent=FALSE,
covar_list=list(CovSnd=roweff), #covariate effects
time_interval="aggregate", ##aggregated estimation
edgelist=edgelist,
algorithm = "parametric")
##time interval estimation
##estimation with time intervals
MEMS(fit.time,
micro_process="CovSnd.1",
macro_function = function(x){igraph::degree(x)},
object_type = "igraph",
nsim=10,
interval=c(0,.1),
silent=TRUE,
covar_list=list(CovSnd=roweff),
time_interval=c(0,5,10,15), #specify three time intervals, 0 - 5, 5 - 10, and 10 - 15
algorithm = "parametric")
########################################################
# Network regression with quadratic assignment procedure
########################################################
library(sna)
##generate network data
set.seed(21093)
x<-rgraph(20,4)
y.l<-x[1,,]+4*x[2,,]+2*x[3,,]
y.p<-apply(y.l,c(1,2),function(a){1/(1+exp(-a))})
y<-rgraph(20,tprob=y.p)
nl<-netlogit(y,x,reps=100)
summary(nl)
MEMS(nl,
micro_process="x2", #should be written as in netlogit output
macro_function = function(x){degree(x)},
object_type = "igraph",
nsim=20,
interval=c(0,1),
silent=FALSE,
full_output = FALSE,
net_logit_y=y,
net_logit_x=x,
algorithm = "nonparametric")
# }
Run the code above in your browser using DataLab