Function to estimate the average micro mediated effect (AMME).

Description

AMME implements parametric and nonparametric estimation routines to estimate the average mediated micro effect. It requires two models. The first is a generative network model (i.e., a model where the dyad, dyad-time period, or dyad-group is the unit of analysis) of the form f(A_{ij}|T_{ij},Z_{ij}), where A is a cross-sectional or longitudinal network or group of longitudinal or cross-sectional networks, T is the possibly endogenous network selection process of interest and Z is a matrix of possibly endogenous confounding selection mechanisms.

The second model is a cross-sectional or longitudinal macro model (i.e., a model where the unit of analysis is a node, subgraph, or network or a combination of nodes, subgraphs, and networks measured collected from multiple settings [such as distinct schools or organizations]) of the form g(Y_i|M_{i},X_{i},T_{i}), where Y_i is the outcome variable, M_i is the mediating macro variable, X_i is a matrix of control variables that possibly vary as a function of selection process T_{ij}, and T_i is the optional unit-level measure of T_{ij}. The AMME is the change in Y_i when T_{ij} allowed to vary versus set to 0 because of an associated change in M_i. The AMME is given by

AMME=\frac{1}{2n} y_i(T_i(t),M_i(T_{ij}),X_i(t))-y_i(T_i(t),M_i(0),X_i(t))

, where n is the number of observations and t=0,T_{ij}. AMME currently accepts the following micro models: glm, glmer, ergm, btergm, sienaFit, rem.dyad, and netlogit objects. The following macro model objects are accepted: lm, glm, lmer, glmer, gam, plm, and lnam objects. Pooled estimation for multiple network models is also implemented for ergm and sienaFit micro models. Both parametric and nonparametric estimation are available.

Usage

AMME(micro_model,
      macro_model,
      micro_process,
      mediator,
      macro_function,
      link_id,
      object_type=NULL,
      controls=NULL,
      control_functions=NULL,
      interval=c(0,1),
      nsim=500,
      algorithm="parametric",
      silent=FALSE,
      full_output=FALSE,
      SAOM_data=NULL,
      SAOM_var=NULL,
      time_interval=NULL,
      covar_list=NULL,
      edgelist=NULL,
      net_logit_y=NULL,
      net_logit_x=NULL,
      group_id=NULL,
      node_numbers=NULL,
      sensitivity_ev=TRUE)

Arguments

Details

Estimates the AMME over the provided intervals. Standard errors and confidence intervals are based on the sampling distribution of simulated values, which are calculated either parametrically or nonparametrically according to algorithm. Parametric estimation is typically faster, but cannot be used for nonparametric network models (e.g., quadratic assignment procedure).

macro_function and control_functions make up the core utilites of AMME. macro_function calculates the mediating variable of interest, while control_functions calculates all control variables that vary as a function of micro_process and potentially confound the effect of mediator. When controls are left NULL, then AMME estimates the AMME without accounting for confounding variables. Specifying controls and control_functions ensures that estimates of the AMME account for alternative pathways from micro_process to the outcome variable in macro_model. In cases where micro_process is included as a predictor variable in macro_model, this can be specified by including the netmediate helper function identity_function into control_functions.

netmediate currently supports functions calculated on igraph and network objects, which should be specified using the object_type argument. These may be functions inherent to the statnet and igraph software package or they may be functions from other packages that accept network/igraph objects. The functions provided to macro_function and control_functions may also be user-defined functions that accept network or igraph objects as inputs and return a numeric value or vector of numeric values as output. It is also possible to over-ride the network and igraph object requirements within a user function. To do so, set the object_type argument (or relevant element within the object_type argument when object_type is a list) to either network or igraph and then define a user-function that accepts a network or igraph object as its input, converts the object to the desired data structure, calculates the statistic of interest, and returns a numeric value or vector of numeric values. See examples below for an illustration.

By default, the AMME is calculated by averaging over the distribution of simulated values. If full_output is set to TRUE, the distribution of simualted statistics is returned. This may be useful when the median or mode of the simulated distribution is required or if the researcher wants to inspect the distributional shape of simulated values.

When the sensitivity_ev argumetn is left at TRUE (the default), AMME will return an E-value and a risk ratio assessing how robust the AMME estimate is to a hypothetical omitted confounding variable. A higher value for either metric means that an omitted variable would have to a large effect on the macro outcome to nullify the effect of the parameterized micro process. See Duxbury and Zhou (2025) for details on interpretation.

AMME also supports pooled estimation for when multiple ergm or sienaFit objects are used as the micro_model. To use pooled estimation, the model parameter should be specified as a list of ergm or sienaFit objects. If using sienaFit, the SAOM_data argument will also need to be specified as an ordered list with elements corresponding to entries in the list of sienaFit objects. Similarly, the SAOM_var parameter will need to be specified as a list of lists, where each entry in the list is, itself, a list containing all varCovar and varDyadCovar objects used to calculate macro statistics of interest. Note that SAOM_var should not be provided if the macro statistic of interest is not a function of the variables contained in varCovar and varDyadCovar.

Value

If full_output=FALSE, then a table is returned with the AMME, its standard error, confidence interval, and p-value.

If full_output=TRUE, then a list is returned with the following three elements.

Author(s)

Duxbury, Scott W. Associate Professor, University of North Carolina–Chapel Hill, Department of Sociology.

Zhao, Xin (Louis), PhD Student, University of North Carolina–Chapel Hill, Department of Sociology

References

Duxbury, Scott W. 2024. "Micro-macro Mediation Analysis in Social Networks." Sociological Methodology.

Duxbury, Scott W., and Xin Zhao. Working paper. "Sensitivity Tests for Micro-Macro Network Analysis."

See Also

MEMS ergm.mma mediate

Examples

##############################
#   Basic AMME specifications
#############################


####create ERGM generative model
library(statnet)
data("faux.mesa.high")
ergm_model<-ergm(faux.mesa.high~edges+
                   nodecov("Grade")+
                   nodefactor("Race")+
                   nodefactor("Sex")+
                   nodematch("Race")+
                   nodematch("Sex")+
                   absdiff("Grade"))


###create node-level data for second stage analysis with
node_level_data<-data.frame(grade=faux.mesa.high%v%"Grade",
                            race=faux.mesa.high%v%"Race",
                            sex=faux.mesa.high%v%"Sex",
                            degree=degree(faux.mesa.high))

node_level_data$senior<-0
node_level_data$senior[node_level_data$grade==max(node_level_data$grade)]<-1
node_level_data$v_id<-1:network.size(faux.mesa.high) #define ID for each observation

probit_model<-glm(senior~race+sex+degree,
                data=node_level_data,
                family=binomial(link="probit"))

###estimate the indirect effect of grade homophily on senior status acting through degree centrality
  #in a model with no network control variables
AMME(micro_model=ergm_model,
     macro_model=probit_model,
     micro_process="absdiff.Grade",
     mediator="degree",
     macro_function=degree,
     link_id=node_level_data$v_id, #specify vertex IDs
     object_type="network",
     interval=c(0,1),
     nsim=50,
     algorithm="parametric",
     silent=FALSE)

#use nonparametric estimation for a generalized additive model
library(gam)

gam_model<-gam(senior~race+sex+s(degree),
               data=node_level_data)

AMME(micro_model=ergm_model,
     macro_model=gam_model,
     micro_process="absdiff.Grade",
     mediator="s(degree)",
     macro_function=degree,
     link_id=node_level_data$v_id,
     object_type="network",
     interval=c(0,1),
     nsim=50,
     algorithm="nonparametric",
     silent=FALSE)



###estimate AMME with linear network autocorrelation model

lnam_model<-lnam(node_level_data$grade,
                 x=as.matrix(node_level_data[,4:5]),
                 W1=as.sociomatrix(faux.mesa.high))


AMME(micro_model=ergm_model,
     macro_model=lnam_model,
     micro_process="absdiff.Grade",
     mediator="degree",
     macro_function=degree,
     link_id=node_level_data$v_id,
     object_type="network",
     interval=c(0,1),
     nsim=50,
     algorithm="parametric",
     silent=FALSE)




############################
#   Including controls
###########################

##single control
node_level_data<-data.frame(grade=faux.mesa.high%v%"Grade",
                            race=faux.mesa.high%v%"Race",
                            sex=faux.mesa.high%v%"Sex",
                            degree=degree(faux.mesa.high),
                            betweenness=betweenness(faux.mesa.high))

node_level_data$senior<-0
node_level_data$senior[node_level_data$grade==max(node_level_data$grade)]<-1
node_level_data$v_id<-1:network.size(faux.mesa.high) #define ID for each observation

probit_model<-glm(senior~race+sex+degree+betweenness,
                  data=node_level_data,
                  family=binomial(link="probit"))


AMME(micro_model=ergm_model,
     macro_model=probit_model,
     micro_process="absdiff.Grade",
     mediator="degree",
     macro_function=degree,
     link_id=node_level_data$v_id, #specify vertex IDs
     controls="betweenness", #should match model output exactly
     control_functions=betweenness,
     object_type="network",
     interval=c(0,1),
     nsim=50,
     algorithm="parametric",
     silent=FALSE)



##multiple controls
##include an AR 1 parameter to make it a nonlinear network autocorrelation model
node_level_data$AR1<-as.sociomatrix(faux.mesa.high)%*%node_level_data$senior
probit_model<-glm(senior~race+sex+degree+betweenness+AR1,
                  data=node_level_data,
                  family=binomial(link="probit"))

#specify user function
ar_function<-function(x){
  return(as.sociomatrix(x)%*%node_level_data$senior)
}


AMME(micro_model=ergm_model,
     macro_model=probit_model,
     micro_process="absdiff.Grade",
     mediator="degree",
     macro_function=degree,
     link_id=node_level_data$v_id,
     controls=c("betweenness","AR1"), #should match model output exactly
     control_functions=list(betweenness,ar_function), #provide functions as a list
     object_type="network",
     interval=c(0,1),
     nsim=50,
     algorithm="parametric",
     silent=FALSE)



##using identity_function when micro_process has a direct effect on y
  #to use identity_function, the control and micro_process need to have the same
  #name and the macro control variable has to be numeric

node_level_data$Sex<-as.numeric(as.factor(node_level_data$sex))
logit_model<-glm(senior~race+Sex+degree+betweenness+AR1,
                  data=node_level_data,
                  family=binomial)



AMME(micro_model=ergm_model,
     macro_model=logit_model,
     micro_process="nodefactor.Sex.M",
     mediator="degree",
     macro_function=degree,
     link_id=node_level_data$v_id,
     controls=c("betweenness","AR1","Sex"), #should match model output exactly
     control_functions=list(betweenness,ar_function,identity_function),
     object_type="network",
     interval=c(0,1),
     nsim=50,
     algorithm="parametric",
     silent=FALSE)





################################
#   More complex data structures
###############################


###############################
# AMME with longitudinal data
##############################

#bootstrap TERGM and panel data model
library(btergm)
library(plm)
data(alliances)

ally_data<-list(LSP[[1]],
                LSP[[2]],
                LSP[[3]])

#fit bootstrap TERGM with 200 replications
bt_model<-btergm(ally_data~edges+
                   gwesp(.7,fixed=T)+
                   mutual,R=200)


#create node data
ally_node_data<-data.frame(outdeg=c(rowSums(LSP[[1]]),rowSums(LSP[[2]]),rowSums(LSP[[3]])),
                           indeg=c(colSums(LSP[[1]]),colSums(LSP[[2]]),colSums(LSP[[3]])))

ally_node_data$v_id<-rep(rownames(LSP[[1]]),3) #create node IDS
ally_node_data$t_id<-c(rep(1, nrow(ally_data[[1]])), #create time IDS
                       rep(2, nrow(ally_data[[1]])),
                       rep(3, nrow(ally_data[[1]])))
ally_node_data$link_id<-paste(ally_node_data$v_id,ally_node_data$t_id)#create node-panel identifiers

ally_node_data$v_id<-as.factor(as.character(ally_node_data$v_id))

#estimate a linear model with node fixed effects
lm_model<- lm(outdeg~indeg +v_id,
          data = ally_node_data)



AMME(micro_model=bt_model,
     macro_model=lm_model,
     micro_process="gwesp.OTP.fixed.0.7",
     mediator="indeg",
     macro_function=function(x){degree(x,cmode="indegree")},
     link_id=ally_node_data$link_id, #provide node-panel identifiers
     object_type="network",
     interval=c(0,1),
     nsim=11,
     algorithm="nonparametric",
     silent=FALSE)


##include controls at different units of analysis
  #include global transitivity statistic at each network panel
transitivity_list<-c(gtrans(as.network(LSP[[1]])),
                     gtrans(as.network(LSP[[2]])),
                     gtrans(as.network(LSP[[3]])))


ally_node_data$transitivity<-c(rep(transitivity_list[1],nrow(LSP[[1]])),
                               rep(transitivity_list[2],nrow(LSP[[2]])),
                               rep(transitivity_list[3],nrow(LSP[[3]])))



lm_model<- lm(outdeg~indeg+transitivity +v_id,
              data = ally_node_data)



AMME(micro_model=bt_model,
     macro_model=lm_model,
     micro_process="gwesp.OTP.fixed.0.7",
     mediator="indeg",
     macro_function=function(x){degree(x,cmode="indegree")},
     link_id=list(ally_node_data$link_id,ally_node_data$t_id),#list of IDs for nodes and time
     controls="transitivity",
     control_functions = gtrans,
     object_type="network",
     interval=c(0,1),
     nsim=11,
     algorithm="nonparametric",
     silent=FALSE)




#SAOM and panel data model with PLM package
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 SAOM
SAOM_model<-siena07(create.model,
                        data=SAOM.Data,
                        effects=SAOM.terms,
                        verbose=TRUE)


##create node-level data
node_level_data<-data.frame(smoking=s50s[,1], #smoking behavior for DV
                            alcohol=s50a[,1],
                            v_id=rownames(s501), #unique node IDS
                            wave="Wave 1",       #unique time IDS
                            outdegree=rowSums(s501),
                            indegree=colSums(s501),
                            AR1=s501%*%s50s[,1],  #assign network autocorrelation
                            gcc=gtrans(as.network(s501)))

node_level_data<-rbind(node_level_data,data.frame(smoking=s50s[,2],
                                                  alcohol=s50a[,2],
                                                  v_id=rownames(s502),
                                                  wave="Wave 2",
                                                  outdegree=rowSums(s502),
                                                  indegree=colSums(s502),
                                                  AR1=s502%*%s50s[,2],
                                                  gcc=gtrans(as.network(s502))))



node_level_data<-rbind(node_level_data,data.frame(smoking=s50s[,3],
                                                  alcohol=s50a[,3],
                                                  v_id=rownames(s503),
                                                  wave="Wave 3",
                                                  outdegree=rowSums(s503),
                                                  indegree=colSums(s503),
                                                  AR1=s503%*%s50s[,3],
                                                  gcc=gtrans(as.network(s503))))


##create unique identifiers for node-panel
node_level_data$unique_ids<-paste(node_level_data$v_id,node_level_data$wave)

##estimate one-way fixed effects model with PLM
library(plm)
FE_model<-plm(smoking~alcohol+outdegree+indegree+AR1+gcc,
               data=node_level_data,
               index=c("v_id","wave"))



##create AR function to provide to AMME
ar_function<-function(x){return(as.sociomatrix(x)%*%(x%v%"Smoking"))}


AMME(micro_model=SAOM_model,
     macro_model=FE_model,
     micro_process="reciprocity",
     mediator="indegree",
     macro_function=function(x){degree(x,cmode="indegree")},
     link_id=list(node_level_data$unique_id,node_level_data$unique_id,
                      node_level_data$unique_id,node_level_data$wave),
     object_type="network",
     controls=c("outdegree","AR1","gcc"),
     control_functions=list(function(x){degree(x,cmode="outdegree")},ar_function,gtrans),
     interval=c(0,.1),
     nsim=500,
     algorithm="parametric",
     silent=FALSE,
     SAOM_data = SAOM.Data,
     SAOM_var=list(Smoking=Smoking,Alcohol=Alcohol)) #provide var_list






################################
# AMME with pooled ERGM and SAOM
################################



#pooled ERGM
  #fit two ERGMs to two networks
data("faux.mesa.high")
model1<-ergm(faux.mesa.high~edges+
               nodecov("Grade")+
               nodefactor("Race")+
               nodefactor("Sex")+
               nodematch("Race")+
               nodematch("Sex")+
               absdiff("Grade"))

data("faux.magnolia.high")
model2<-ergm(faux.magnolia.high~edges+
               nodecov("Grade")+
               nodefactor("Race")+
               nodefactor("Sex")+
               nodematch("Race")+
               nodematch("Sex")+
               absdiff("Grade"))


#create node level data
node_level_data<-data.frame(grade=faux.mesa.high%v%"Grade",
                            sex=faux.mesa.high%v%"Sex",
                            degree=degree(faux.mesa.high),
                            betweenness=betweenness(faux.mesa.high),
                            gcc=gtrans(faux.mesa.high),
                            net_id="Mesa")

node_level_data$senior<-0
node_level_data$senior[node_level_data$grade==max(node_level_data$grade)]<-1
node_level_data$v_id<-1:network.size(faux.mesa.high)


node_level_data2<-data.frame(grade=faux.magnolia.high%v%"Grade",
                            sex=faux.magnolia.high%v%"Sex",
                            degree=degree(faux.magnolia.high),
                            betweenness=betweenness(faux.magnolia.high),
                            gcc=gtrans(faux.magnolia.high),
                            net_id="Magnolia")

node_level_data2$senior<-0
node_level_data2$senior[node_level_data$grade==max(node_level_data2$grade)]<-1
node_level_data2$v_id<-206:(network.size(faux.magnolia.high)+205)
node_level_data<-rbind(node_level_data,node_level_data2)


#estimate glm macro model with an AR 1 process
probit_model<-glm(senior~sex+degree+betweenness+gcc,
                data=node_level_data,
                family=binomial(link="probit"))



AMME(micro_model=list(model1,model2),
     macro_model=probit_model,
     micro_process="nodematch.Sex",
     mediator="degree",
     macro_function=degree,
     link_id=list(node_level_data$v_id,node_level_data$v_id,node_level_data$net_id),
     object_type="network",
     controls=c("betweenness","gcc"),
     control_functions=list(betweenness,gtrans),
     interval=c(0,1),
     nsim=50,
     algorithm="parametric",
     silent=FALSE)



##pooled SAOM with control functions using time varying covariates

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 SAOM
SAOM_model<-siena07(create.model,
                    data=SAOM.Data,
                    effects=SAOM.terms,
                    verbose=TRUE)


##create node-level data
node_level_data<-data.frame(smoking=s50s[,1], #smoking behavior for DV
                            alcohol=s50a[,1],
                            v_id=rownames(s501), #unique node IDS
                            wave="Wave 1",       #unique time IDS
                            outdegree=rowSums(s501),
                            indegree=colSums(s501),
                            AR1=s501%*%s50s[,1],  #assign network autocorrelation
                            gcc=gtrans(as.network(s501)))

node_level_data<-rbind(node_level_data,data.frame(smoking=s50s[,2],
                                                  alcohol=s50a[,2],
                                                  v_id=rownames(s502),
                                                  wave="Wave 2",
                                                  outdegree=rowSums(s502),
                                                  indegree=colSums(s502),
                                                  AR1=s502%*%s50s[,2],
                                                  gcc=gtrans(as.network(s502))))



node_level_data<-rbind(node_level_data,data.frame(smoking=s50s[,3],
                                                  alcohol=s50a[,3],
                                                  v_id=rownames(s503),
                                                  wave="Wave 3",
                                                  outdegree=rowSums(s503),
                                                  indegree=colSums(s503),
                                                  AR1=s503%*%s50s[,3],
                                                  gcc=gtrans(as.network(s503))))


#recycle the same model for illustrative purposes
node_level_data$net_ID<-"Model 1"
node_level_data<-rbind(node_level_data,node_level_data)
node_level_data$net_ID[151:300]<-"Model 2"

##create unique identifiers for node-panel
  #ID for node-panel-model
node_level_data$unique_id<-paste(node_level_data$v_id,node_level_data$wave,node_level_data$net_ID)
  #ID for panel-model
node_level_data$unique_waves<-paste(node_level_data$wave,node_level_data$net_ID)

#estimate a linear network autocorrelation model with node fixed effects
FE_model<-lm(smoking~alcohol+outdegree+indegree+AR1+gcc+v_id,
              data=node_level_data)



##create user function calculate AR1 process on time varying node attributes
ar_function<-function(x){return(as.sociomatrix(x)%*%(x%v%"Smoking"))}

##estimate AMME
AMME(micro_model=list(SAOM_model,SAOM_model), #provide list of sienaFit objects
     macro_model=FE_model,
     micro_process="reciprocity",
     mediator="indegree",
     macro_function=function(x){degree(x,cmode="indegree")},
     link_id=list(node_level_data$unique_id,node_level_data$unique_id,
                        node_level_data$unique_id,node_level_data$unique_waves),
     object_type="network",
     controls=c("outdegree","AR1","gcc"),
     control_functions=list(function(x){degree(x,cmode="outdegree")},ar_function,gtrans),
     interval=c(0,.1),
     nsim=100,                  #parametric estimation requires more simulations than coefficients
     algorithm="parametric",
     silent=FALSE,
     SAOM_data = list(SAOM.Data,SAOM.Data), #list of siena objects
     SAOM_var=list(list(Smoking=Smoking,Alcohol=Alcohol),#provide var_list
                  list(Smoking=Smoking,Alcohol=Alcohol)))





#################################
# AMME with nested data
################################

####create dyad-level data

library(lme4)
library(btergm)
##use small data to simplify estimation
glm_dat<-edgeprob(model1)
glm_dat$net_id<-"mesa"
glm_dat2<-edgeprob(model2)
glm_dat2$net_id<-"magnolia"
glm_dat<-rbind(glm_dat,glm_dat2[,-c(4)])


##estimate micro model as glm for btoh networks using pooled ERGM data
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)


#create macro data
node_level_data<-data.frame(grade=faux.mesa.high%v%"Grade",
                            sex=faux.mesa.high%v%"Sex",
                            degree=degree(faux.mesa.high),
                            betweenness=betweenness(faux.mesa.high),
                            gcc=gtrans(faux.mesa.high),
                            net_id="Mesa")

node_level_data$senior<-0
node_level_data$senior[node_level_data$grade==max(node_level_data$grade)]<-1
node_level_data$v_id<-1:network.size(faux.mesa.high)


node_level_data2<-data.frame(grade=faux.magnolia.high%v%"Grade",
                             sex=faux.magnolia.high%v%"Sex",
                             degree=degree(faux.magnolia.high),
                             betweenness=betweenness(faux.magnolia.high),
                             gcc=gtrans(faux.magnolia.high),
                             net_id="Magnolia")

node_level_data2$senior<-0
node_level_data2$senior[node_level_data$grade==max(node_level_data2$grade)]<-1
node_level_data2$v_id<-206:(network.size(faux.magnolia.high)+205)
node_level_data<-rbind(node_level_data,node_level_data2)


#estimate glm macro model
probit_model<-glm(senior~sex+degree+betweenness+gcc,
                  data=node_level_data,
                  family=binomial(link="probit"))




AMME(micro_model=net_glm,
     macro_model=probit_model,
     micro_process="nodematch.Sex",
     mediator="degree",
     macro_function=degree,
     link_id=list(node_level_data$v_id,node_level_data$v_id,node_level_data$net_id),
     object_type="network",
     controls=c("betweenness","gcc"),
     control_functions=list(betweenness,gtrans),
     interval=c(0,1),
     nsim=50,
     algorithm="parametric",
     silent=FALSE,
     group_id=glm_dat$net_id,
     node_numbers = c(network.size(faux.mesa.high),
                      network.size(faux.magnolia.high)))




###using glmer for micro model
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)

probit_glmer<-glm(senior~sex+degree+betweenness+gcc,
                data=node_level_data,
                family=binomial(link="probit"))


AMME(micro_model=net_glm,
     macro_model=probit_glmer,
     micro_process="nodematch.Sex",
     mediator="degree",
     macro_function=degree,
     link_id=list(node_level_data$v_id,node_level_data$v_id,node_level_data$net_id),
     object_type="network",
     controls=c("betweenness","gcc"),
     control_functions=list(betweenness,gtrans),
     interval=c(0,1),
     nsim=50,
     algorithm="parametric",
     silent=FALSE,
     group_id=glm_dat$net_id,
     node_numbers = c(network.size(faux.mesa.high),
                      network.size(faux.magnolia.high)))

Function to estimate the micro effect on macro structure (MEMS).

Description

MEMS implements parametric and nonparametric estimation routines to estimate the micro effect on macro structure when using a generative network model (i.e., a model where the dyad, dyad-time period, or dyad-group is the unit of analysis). The MEMS is defined in postestimation as a function of the possibly endogenous micro process X, which is assumed to be a predictor in the micro model of the form A=f(\theta X + \gamma ^TZ), where Z is a matrix of possibly endogenous controls and A is the network of interest. The MEMS is given by

MEMS=\sum_i \frac{M(\theta, X, \gamma, Z)_i-M(\gamma, Z)_i}{n}

, for n observations. Tuning parameters can be assigned to toggle the strength of \theta in model-implied estimates of MEMS. MEMS currently accepts glm, glmer, ergm, btergm, sienaFit, rem.dyad, and netlogit objects and implements both parametric and nonparametric estimation. Pooled estimation for multiple network models is also implemented for ergm and sienaFit objects.

Usage

MEMS(model,
      micro_process,
      macro_function,
      object_type=NULL,
      interval=c(0,1),
      nsim=500,
      algorithm="parametric",
      silent=FALSE,
      full_output=FALSE,
      SAOM_data=NULL,
      SAOM_var=NULL,
      time_interval=NULL,
      covar_list=NULL,
      edgelist=NULL,
      net_logit_y=NULL,
      net_logit_x=NULL,
      group_id=NULL,
      node_numbers=NULL,
      mediator=NULL,
      link_id=NULL,
      controls=NULL,
      control_functions=NULL,
      sensitivity_ev=TRUE)

Arguments

.

Details

Estimates the MEMS over the provided intervals. If the macro statistic is calculated on the node or subgraph levels or on multiple network observations, the aMEMS is provided instead. Standard errors and confidence intervals are based on the sampling distribution of simulated values, which are calculated either parametrically or nonparametrically according to algorithm. Parametric estimation is typically faster, but cannot be used for nonparametric network models (e.g., quadratic assignment procedure).

macro_function is the workhorse component of MEMS. The function should calculate the macro statistic of interest. netmediate currently supports functions calculated on igraph and network objects, which should be specified as using the object_type argument. These may be functions inherent to the statnet and igraph software package or they may be functions from other packages that accept network/igraph objects. They may also be user-defined functions that accept network or igraph objects as input and return a numeric value or vector of numeric values as output. It is also possible to over-ride the network and igraph object requirements within a user function. To do so, set the object_type argument to either network or igraph and then define a user-function that accepts a network or igraph object as its input, converts the object to the desired data structure, calculates the statistic of interest, and finally returns a numeric value or vector of numeric values. See examples below for an illustration.

By default, the MEMS is provided by averaging over the distribution of simulated values. If full_output is set to TRUE, the entire distribution of simualted statistics is returned. This may be useful when the median or mode of the simulated distribution is required or if the researcher wants to inspect the distributional shape of simulated values. When the sensitivity_ev argumetn is left at TRUE (the default), MEMS will return an E-value and a risk ratio assessing how robust the MEMS estimate is to a hypothetical omitted confounding variable. A higher value for either metric means that an omitted variable would have to a large effect on the macro outcome to nullify the effect of the parameterized micro process. See Duxbury and Zhou (2025) for details on interpretation.

MEMS also supports pooled estimation for multiple ergm or sienaFit objects. To use pooled estimation, the model parameter should be specified as a list of ergm or sienaFit objects. If using sienaFit, the SAOM_data argument will also need to be specified as an ordered list with elements corresponding to entries in the list of sienaFit objects. Similarly, the SAOM_var parameter will need to be specified as a list of lists, where each entry in the list is, itself, a list containing all varCovar and varDyadCovar objects used to calculate macro statistics of interest. Note that SAOM_var should not be provided if the macro statistic of interest is not a function of the variables contained in varCovar and varDyadCovar.

When estimating a relational event model with a rem.dyad object, time_interval can be specified to provide exact time intervals over which to induce unique networks. This utility is often useful when combining rem.dyad estimation with AMME when the macro_model is panel data with coarse timing information. The same behavior can be obtained when estimating a relational event model using glm or glmer by assigning the desired time intervals in the model matrix and then providing the vector of time intervals to the group_id parameter when calling MEMS.

Value

If full_output=FALSE, then a table is returned with the MEMS, its standard error, confidence interval, and p-value.

If full_output=TRUE, then a list is returned with the following three elements.

Author(s)

Duxbury, Scott W. Associate Professor, University of North Carolina–Chapel Hill, Department of Sociology.

Zhao, Xin (Louis). PhD Candidate, University of North Carolina–Chapel Hill, Department of Sociology

References

Duxbury, Scott W. 2024. "Micro Effects on Macro Structure in Social Networks." Sociological Methodology.

Wertsching, Jenna, and Scott W. Duxbury. Working paper. "Micro Effects on Macro Structure: Identification, Comparison between Nested Models, and Sensitivity Tests for Functional Form."

Duxbury, Scott W., and Xin Zhao. Working paper. "Sensitivity Tests for Micro-Macro Network Analysis."

See Also

AMME ergm.mma mediate

Examples

########################################
# 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")

Function to calculate Moran's first order network autocorrelation in co-evolution SAOM.

Description

A function to calculate Moran's first order network autocorrelation in co-evolution SAOM using both the endogenous dependent variables (behavior and network functions).

Usage

Moran_dv(network)

Arguments

Value

No return value, used internally with MEMS and AMME

Function to compare micro effect on macro structure (MEMS) estimates between models.

Description

compare_MEMS implements parametric and nonparametric routines to compare MEMS estimate between models. When compared between nested models, compare_MEMS results can be interpreted as the portion of a MEMS explained by a mediating or confounding variable. When compared between models with distinct functional forms and the same specification, compare_MEMS results can be interpreted as the sensitivity of MEMS results to decision about model functional form as described by Wertsching and Duxbury (2025).

compare_MEMS can also be used as a test of the difference between two MEMS estimates within the same model. This is useful when researchers watn to formally evaluate whether the MEMS for one explanatory micro mechanism is significantly larger or smaller than a second explanatory micro mecahnism.

The difference in MEMS is the change in MEMS after one or more micro-processes are included into a model or, in the case of sensitivity tests, when the functional form is changed. Let MEMS_p represent the MEMS obtained from a model that omits one or more intervening variables and MEMS_f be the MMES obtained from a model that includes the intervening variable(s). The change in MEMS is given

\Delta MEMS=MEMS_p-MEMS_f

. MEMS_p and MEMS_f may also be have the same specification but use distinct functional forms or other modeling decisions in the case of sensitivity tests. Tuning parameters can be assigned to toggle the strength of \theta in model-implied estimates of MEMS. MEMS currently accepts glm, glmer, ergm, btergm, sienaFit, rem.dyad, and netlogit objects and implements both parametric and nonparametric estimation. Pooled estimation for multiple network models is also implemented for ergm and sienaFit objects.

Usage

compare_MEMS(partial_model,
      full_model,
      micro_process,
      micro_process2=NULL,
      macro_function,
      macro_function2=NULL,
      object_type=NULL,
      interval=c(0,1),
      nsim=500,
      algorithm="parametric",
      silent=FALSE,
      full_output=FALSE,
      sensitivity_ev=TRUE,
      SAOM_data=NULL,
      SAOM_var=NULL,
      time_interval=NULL,
      covar_list=NULL,
      edgelist=NULL,
      net_logit_y=NULL,
      net_logit_x=NULL,
      group_id=NULL,
      node_numbers=NULL,
      mediator=NULL,
      link_id=NULL,
      controls=NULL,
      control_functions=NULL)

Arguments

Details

Compares MEMS estimates between two models. If one or more confounding or intervening variables are excluded or included between models, the change in MEMS can be interpreted as the portion of the MEMS explained by one or more confounding or intervening variable. If two models are provided with the same specification but a distinct functional form, the change in MEMS is a sensitivty test of how much the MEMS estimate changes because of a model decision. This can be useful, for example, when comparing TERGM and SAOM estimates as each models make distinct assumptions about sources of network change and the temporal ordering of tie changes.

compare_MEMS functionality inherits directly from the MEMS command. See the MEMS page for more details.

Value

If full_output=FALSE, then a table is returned with the change MEMS, its standard error, confidence interval, and p-value, and the same results for the partial and complete MEMS.

If full_output=TRUE, then a list is returned with the following three elements.

Author(s)

Duxbury, Scott W. Associate Professor, University of North Carolina–Chapel Hill, Department of Sociology.

References

Duxbury, Scott W. 2024. "Micro Effects on Macro Structure in Social Networks." Sociological Methodology.

Wertsching, Jenna, and Scott W. Duxbury. Working paper. "Micro Effects on Macro Structure: Identification, Comparison between Nested Models, and Sensitivity Tests for Functional Form."

Duxbury, Scott W., and Xin Zhao. Working paper. "Sensitivity Tests for Micro-Macro Network Analysis."

See Also

AMME MEMS ergm.mma mediate mediate_MEMS

Examples

##############
# 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)

Function to map micro_process onto macro_model within calls to AMME.

Description

A function to control for a node-level micro_process in AMME estimation.

Usage

identity_function(x)

Arguments

Value

No return value, used internally with AMME

Function to conduct mediation analysis with micro effects on macro structure (MEMS).

Description

mediate_MEMS implements parametric and nonparametric routines to compare MEMS estimate between models. Largely a wrapper for compare_MEMS, the sole novel functionality of mediate_MEMS is provided when a user specifies the model_comparison argument to be TRUE. When model_comparison is set to TRUE, mediate_MEMS compares the direct, indirect, and total MEMS estimates to those obtained from partial_model2 and full_model2. This can be used as a sensitivity test to researchers' choice of model. It can also provide the basis for testing differences in direct, total, and indirect effect sizes within the same model by setting partial_model2 to equal partial_model and full_model2 to equal full_model and provide a second explanatory micro process to the micro_process2 argument. In these cases, the difference between "models" captures the differences effect size of the direct, indirect, and total MEMS estimates for two distinct explanatory micro processes. See Wertsching and Duxbury (2025) for details.

The difference in MEMS is the change in MEMS after one or more micro-processes are included into a model or, in the case of sensitivity tests, when the functional form is changed. Let MEMS_p represent the MEMS obtained from a model that omits one or more intervening variables and MEMS_f be the MMES obtained from a model that includes the intervening variable(s). The change in MEMS is given

\Delta MEMS=MEMS_p-MEMS_f

. MEMS_p and MEMS_f may also be have the same specification but use distinct functional forms or other modeling decisions in the case of sensitivity tests. Tuning parameters can be assigned to toggle the strength of \theta in model-implied estimates of MEMS. MEMS currently accepts glm, glmer, ergm, btergm, sienaFit, rem.dyad, and netlogit objects and implements both parametric and nonparametric estimation. Pooled estimation for multiple network models is also implemented for ergm and sienaFit objects.

Usage

mediate_MEMS(partial_model,
      full_model,
      micro_process,
      macro_function,
      model_comparison=FALSE,
      partial_model2=NULL,
      full_model2=NULL,
      micro_process2=NULL,
      macro_function2=NULL,
      object_type=NULL,
      interval=c(0,1),
      nsim=500,
      algorithm="parametric",
      silent=FALSE,
      full_output=FALSE,
      sensitivity_ev=TRUE,
      SAOM_data=NULL,
      SAOM_var=NULL,
      time_interval=NULL,
      covar_list=NULL,
      edgelist=NULL,
      net_logit_y=NULL,
      net_logit_x=NULL,
      group_id=NULL,
      node_numbers=NULL,
      mediator=NULL,
      link_id=NULL,
      controls=NULL,
      control_functions=NULL)

Arguments

Details

Compares MEMS estimates between two models. If one or more confounding or intervening variables are excluded or included between models, the change in MEMS can be interpreted as the portion of the MEMS explained by one or more confounding or intervening variable. If two models are provided with the same specification but a distinct functional form, the change in MEMS is a sensitivty test of how much the MEMS estimate changes because of a model decision. This can be useful, for example, when comparing TERGM and SAOM estimates as each models make distinct assumptions about sources of network change and the temporal ordering of tie changes.

If a single pair of models are compared, mediate_MEMS results will be identical to compare_MEMS results. However, if model_comparison is set to TRUE and a second pair of models are provided using a different set of control variables or a distinct functional form, mediate_MEMS will also return sensitivity tests evaluating whether the partial, full, and indirect MEMS estimates are significantly different between the two sets of models. This functionality can also be used to formally test differences in partial, full, and indirect effect sizes by setting partial_model2 to be the same as partial model and full_model2 to be the same as full_model and setting micro_process2 to be a different variable than micro_process.

mediate_MEMS functionality inherits directly from the compare_MEMS and MEMS commands. See the compare_MEMS and MEMS pages for more details.

Value

If full_output=FALSE and model_comparison=FALSE, then a table is returned with the change MEMS, its standard error, confidence interval, and p-value, and the same results for the partial and complete MEMS.

If full_output=TRUE and model_comparison=FALSE, then a list is returned with the following three elements.

If full_output=FALSE and model_comparison=TRUE, then a list is returned with the following elements.

If full_output=TRUE and model_comparison=TRUE, then a list is returned with the following three elements.

Author(s)

Duxbury, Scott W. Associate Professor, University of North Carolina–Chapel Hill, Department of Sociology.

References

Duxbury, Scott W. 2024. "Micro Effects on Macro Structure in Social Networks." Sociological Methodology.

Wertsching, Jenna, and Scott W. Duxbury. Working paper. "Micro Effects on Macro Structure: Identification, Comparison between Nested Models, and Sensitivity Tests for Functional Form."

Duxbury, Scott W., and Xin Zhao. Working paper. "Sensitivity Tests for Micro-Macro Network Analysis."

See Also

AMME MEMS ergm.mma mediate compare_MEMS

Examples

##############
# Not run
###############
library(statnet)
library(igraph)
data("faux.mesa.high")

#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 results from compare_MEMS and mediate_MEMS

compare_MEMS(partial_model=model,
             full_model=model2,
             micro_process="nodematch.Race",
             macro_function=transitivity,
             object_type = "igraph",
             silent=FALSE,
             algorithm="parametric")


mediate_MEMS(partial_model=model,
             full_model=model2,
             micro_process="nodematch.Race",
             macro_function=transitivity,
             object_type = "igraph",
             silent=FALSE,
             algorithm="parametric")




###test sensitivity to MPLE versus MCMC MLE estimation

modela<-ergmMPLE(faux.mesa.high~edges+nodecov("Grade")+nodefactor("Race")+
              nodefactor("Sex")+nodematch("Race")+nodematch("Sex")+absdiff("Grade"),
              output="fit")

model2a<-ergmMPLE(faux.mesa.high~edges+nodecov("Grade")+nodefactor("Race")+
               nodefactor("Sex")+nodematch("Race")+nodematch("Sex")+absdiff("Grade")+
               gwesp(.5,fixed=TRUE),
               output="fit")




mediate_MEMS(partial_model=model,
             full_model=model2,
             micro_process="nodematch.Race",
             model_comparison  = TRUE,
             partial_model2=modela,
             full_model2=model2a,
             micro_process2="nodematch.Race",
             macro_function=transitivity,
             object_type = "igraph",
             silent=FALSE,
             algorithm="parametric")


##compare direct, total, and indirect effect sizes

mediate_MEMS(partial_model=model,
             full_model=model2,
             micro_process="nodematch.Race",
             model_comparison  = TRUE,
             partial_model2=model,
             full_model2=model2,
             micro_process2="absdiff.Grade",
             macro_function=transitivity,
             object_type = "igraph",
             silent=FALSE,
             algorithm="parametric")





#####################################################
# More complicated sensitivty test using macro function 2
#####################################################

#are the direct, total, and indirect MEMS of
  #network selection on similar smoking behavior on
  #students' shared smoking robust to distinct
  #model choices?

###Compare between SAOM and TERGM treating behavioral
 #(smoking) autocorrelation as outcome, smoking homophily
 #as treatment, and triadic closure as mediator




########################################
#       Co-evolution SAOM
#######################################

library(RSiena)
network_array<-list(s501,s502,s503)
smoking<-as.data.frame(s50s)

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


##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")



###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 WITHOUT transitive ties

create.model<-sienaAlgorithmCreate(projname="Co-evolution_output",
                                   seed=21093,
                                   nsub=4,
                                   n3=1000)
TLSmodel_notrans<-siena07(create.model,
                  data=SAOM.Data,
                  effects=SAOM.terms,
                  verbose=TRUE,
                  returnDeps=TRUE)


#include transTies
SAOM.terms<-includeEffects(SAOM.terms,transTies,inPop)
TLSmodel<-siena07(create.model,
                  data=SAOM.Data,
                  effects=SAOM.terms,
                  verbose=TRUE,
                  returnDeps=TRUE)





#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_nogwesp<-tergm(net_list~Form(
  ~edges+
    mutual+
    gwidegree(.5,fixed=TRUE)+
    nodeicov("smoking")+
    nodeocov("smoking")+
    nodematch("smoking")+
    nodeicov("alcohol")+
    nodeocov("alcohol")+
    absdiff("alcohol")),
  estimate="CMLE"

)


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 direct, total, and indirect
  #MEMS estimates

mediate_MEMS(partial_model=TLSmodel_notrans,
             full_model=TLSmodel,
             micro_process="same Behavior",
             macro_function =Moran_dv,
             model_comparison = TRUE,
             partial_model2=TERGM_1_nogwesp,
             full_model2=TERGM_1,
             micro_process2="Form~nodematch.smoking",
             macro_function2=Moran_tergm,
             object_type = "network",
             SAOM_data = SAOM.Data,
             silent=FALSE,
             nsim=100)