# \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"))
compare_MEMS(partial_model=a,
full_model=a,
micro_process="nodeifactor.var.1.1",
macro_function=gtrans,
nsim=20,
silent=TRUE)
# }
# \donttest{
##############
# Not run
###############
library(statnet)
library(igraph)
data("faux.mesa.high")
####################
###mediation analysis
####################
#how much of the effect of racial homophily on transitivity
#is explained by triadic closure effects?
model<-ergm(faux.mesa.high~edges+nodecov("Grade")+nodefactor("Race")+
nodefactor("Sex")+nodematch("Race")+nodematch("Sex")+absdiff("Grade"))
model2<-ergm(faux.mesa.high~edges+nodecov("Grade")+nodefactor("Race")+
nodefactor("Sex")+nodematch("Race")+nodematch("Sex")+absdiff("Grade")+
gwesp(.5,fixed=TRUE))
compare_MEMS(partial_model=model,
full_model=model2,
micro_process="nodematch.Race",
macro_function=transitivity,
object_type = "igraph",
silent=FALSE,
algorithm="parametric")
########################
# Effect size comparison
########################
#Is the effect of racial homophily on transitivity larger or smaller
#than the effect of grade homophily?
compare_MEMS(partial_model=model2,
full_model=model2,
micro_process="nodematch.Race",
micro_process2="absdiff.Grade",
macro_function=transitivity,
object_type = "igraph",
silent=FALSE,
algorithm="parametric")
###################
# Robustness check
##################
#Are MEMS estimates derived from ERGM robust to alternative MPLE estimation strategies?
model_MPLE<-ergmMPLE(faux.mesa.high~edges+nodecov("Grade")+nodefactor("Race")+
nodefactor("Sex")+nodematch("Race")+nodematch("Sex")+absdiff("Grade")+
gwesp(.5,fixed=TRUE),
output="fit")
compare_MEMS(partial_model=model2,
full_model=model_MPLE,
micro_process="gwesp.fixed.0.5",
macro_function=transitivity,
object_type = "igraph",
silent=FALSE,
algorithm="parametric",
sensitivity_ev=FALSE)
###Compare between SAOM and TERGM
#treating behavioral (smoking) autocorrelation
#as outcome
library(RSiena)
##we'll load RSiena since we're using data from here.
###create a list of adjacency matrices
network_array<-list(s501,s502,s503)
smoking<-as.data.frame(s50s) ##load smoking data
##for our analysis, we'll look at binary smoking behavior
for(i in 1:ncol(smoking)){
smoking[,i][smoking[,i]>1]<-2
}
alcohol<-as.data.frame(s50a) ##we'll use alcohol consumption as a covariate as well
########################################
# Co-evolution model
#######################################
##create "sienaDependent" object,
TLSnet<-sienaDependent(array(c(network_array[[1]],
network_array[[2]],
network_array[[3]]),
dim=c(50,50,3)))
TLSbeh<-sienaDependent(as.matrix(smoking),type="behavior")
#set covariates
Alcohol<-varCovar(as.matrix(alcohol))
###create dataset, but specify network AND behavior
SAOM.Data<-sienaDataCreate(Network=TLSnet,
Behavior=TLSbeh,
Alcohol)
###Create the effects object
SAOM.terms<-getEffects(SAOM.Data)
###We'll start by specifying the NETWORK function
SAOM.terms<-includeEffects(SAOM.terms,egoX,altX,absDiffX,interaction1="Alcohol")
SAOM.terms<-includeEffects(SAOM.terms,egoX,altX,sameX,interaction1="Behavior")
SAOM.terms<-includeEffects(SAOM.terms,transTies,inPop)
###Now let's specify the BEHAVIOR function
SAOM.terms<-includeEffects(SAOM.terms,effFrom,name="Behavior",
interaction1="Alcohol")
SAOM.terms<-includeEffects(SAOM.terms,totSim,name="Behavior",
interaction1="Network")
SAOM.terms<-includeEffects(SAOM.terms,isolate,
name="Behavior",interaction1="Network")
#estimate the model
create.model<-sienaAlgorithmCreate(projname="Co-evolution_output",
seed=21093,
nsub=4,
n3=1000)
TLSmodel<-siena07(create.model,
data=SAOM.Data,
effects=SAOM.terms,
verbose=TRUE,
returnDeps=TRUE)
TLSmodel
#now fit the TERGM
library(statnet)
net_list<-list(as.network(s501),as.network(s502),as.network(s503))
net_list[[1]]<-network::set.vertex.attribute(net_list[[1]],"smoking",s50s[,1])
net_list[[2]]<-network::set.vertex.attribute(net_list[[2]],"smoking",s50s[,2])
net_list[[3]]<-network::set.vertex.attribute(net_list[[3]],"smoking",s50s[,3])
net_list[[1]]<-network::set.vertex.attribute(net_list[[1]],"alcohol",s50a[,1])
net_list[[2]]<-network::set.vertex.attribute(net_list[[2]],"alcohol",s50a[,2])
net_list[[3]]<-network::set.vertex.attribute(net_list[[3]],"alcohol",s50a[,3])
TERGM_1<-tergm(net_list~Form(
~edges+
mutual+
gwesp(.5,fixed=TRUE)+
gwidegree(.5,fixed=TRUE)+
nodeicov("smoking")+
nodeocov("smoking")+
nodematch("smoking")+
nodeicov("alcohol")+
nodeocov("alcohol")+
absdiff("alcohol")),
estimate="CMLE"
)
#create network autocorrelation function for TERGM
Moran_tergm<-function(x){
y<-network::get.vertex.attribute(x,"smoking")
return(nacf(x,y,type="moran",lag=1)[2])
}
#test difference in TERGM and SAOM estimates of the MEMS for
#network selection on same smoking behavior for
#smoking similarity at the aggregate level
compare_MEMS(partial_model=TLSmodel,
full_model=TERGM_1,
micro_process="same Behavior",
macro_function =Moran_dv,
micro_process2="Form~nodematch.smoking",
macro_function2=Moran_tergm,
object_type = "network",
SAOM_data = SAOM.Data,
silent=FALSE)
# }
Run the code above in your browser using DataLab