In this lab, we will be testing a number of hypotheses about a
network’s structure using exponential random graph modeling (ERGM)
techniques using the statnet package in R.1 Mark S. Handcock, David R. Hunter, Carter T. Butts, Steven M. Goodreau, and Martina Morris (2003). statnet: Software tools for the Statistical Modeling of Network Data. statnetproject.org; see also?? statnet. For more information about ERGMs, see generally D. Lusher, J. Koskinen, & G. Robins (2012) Exponential Random Graph Models for Social Networks.statnet
provides a comprehensive framework for ERGM-based network modeling,
including tools for model estimation, model evaluation, model-based
network simulation, and network visualization. This functionality is
powered by a central Markov chain Monte Carlo Maximum Likelihood
Estimation (MCMCMLE) algorithm.2 For a great introduction to MCMC grounded in graph theory, see Jeremy Kun, Markov Chain Monte Carlo Without All the Bullshit.
We will analyze the communication behaviors within a team of
seventeen members who were involved in designing military installations.
CRIeq.txt: each team member’s communication to
retrieve information from other team members on the topic of
environmental quality (eq). This is a directed, binary relation.
CAIeq.txt: each team member’s communication to
allocate information to other team members on the topic of environmental
quality (eq). This is a directed, binary relation.
EXeq_cons.txt: each team member’s expertise on the
topic of environmental quality (eq) as perceived on average by all team
members. This is a continuous attribute.
Hypotheses
We will test various hypotheses based on the Theory of Transactive Memory.3See Monge & Contractor (2003) Theories of Communication Networks, 198—203.
Hypothesis 1: Individuals are less likely to retrieve information from those who retrieve information from them.
Hypothesis 2a: Information retrieval tends to be transitive. That is, if individual i retrieves information from individual k, and individual k retrieves information from individual j, then individual i is more likely to retrieve information from individual j.
Hypothesis 2b: Transitivity increases at a sub-linear rate as a function of the number of ties.
Hypothesis 3a: Individuals tend to retrieve information from other members with high expertise.
Hypothesis 3b: Individuals with low expertise tend to retrieve information from many others.
Hypothesis 4: Individuals tend to retrieve information from members to whom they allocate information.
Part I. Building & Visualizing the Networks (30 pts)
The analysis will use three files: the CRIeq.txt as the network file, EXeq_cons.txt as the attribute file, and CAIeq.txt
as the co-variate network file. To begin, we must convert the data
files into matrices, transform those matrices into networks, and attach
the attribute file to our base network.
Understanding the Base Network
Let’s begin by looking at the summary of our base network.
Analysis
In your own words, explain what this network respresents and
its relationship to our attribute information and the other network. ANSWER:
This network represents 17 nodes, or workers, who go to each other in
search of information. In total, there are 41 links between the nodes.
These relationships are directed, meaning they are not necessarily
reciprocal. Some nodes are more sought after than others. The other
network shows the relationships between nodes when allocating
information to others.
Visualization
Before we conduct further analysis, let’s visualize our base network.
Similar to our approach in Lab 2, we will begin by establishing set
coordinates for our nodes in order to simplify visual comparisons.
Base Graph Structure
Next we can visualize our base network.
Base Network: Retrieval of Environmental Quality
Next, we will size the nodes by the their expertise value.
Base Network: Retrieval of Environmental Quality, Nodes Sized by Expertise Score
Let’s compare this visualization to sizing by in-degree centrality.
Base Network: Nodes by Indegree Centrality
Analysis
Consider hypothesis 3a from Part I. Do these visualizations
prove or disprove the hypothesis? In your own words, interpret the
graphs and explain how they support or reject the hypothesis. ANSWER:
It supports the hypothesis. Individual with larger expertise value tend
to have higher in-degree centrality (9, 6, 16, 5, 4)
Understanding the Covariate Network
Let’s explore the summary statistics of our co-variate network.
In your own words, explain what this network represents and its relationship to the other network and the attribute information. ANSWER:
This network represents the same 17 nodes, or workers, who allocate
information to other nodes. There are fewer links between nodes – 23.
These relationships are also directed, meaning they are not necessarily
reciprocal. Again, some nodes are more sought after than others.
Visualization
We will repeat the visualization process for our co-variate network.
Observe the location and distribution of edges in the following
visualization.
Covariate Network: Allocation of Environmental Quality
Next, we will size the nodes by their expertise scores.
Base Network: Allocation of Environmental Quality, Nodes Sized by Expertise Score
Allocation
Analysis
Consider hypothesis 4 from Part I. Do the visualizations of
the retrieval and allocation networks sized by expertise support or
disprove the hypothesis based on visual inspection? In your own words,
interpret the graphs and explain how they support or reject the
hypothesis. ANSWER: Partially, but not as strongly as
hypothesis 3a Nodes with higher expertise tend to have higher in-degree
centrality, but it’s not the case for every high-expertise node. There
are also examples of low expertise nodes that have high in-degree
centrality.
Part II: Constructing & Analyzing the ERGM Model (70 pts)
Next, we’re going to construct an Exponential Random Graph Model. Note that model construction
is integral to this process; ERGM is not a single method of analysis,
but a type of modelling that requires a theoretical grounding specific
to the network and the hypotheses posed by the researcher. While the
base ERGM simulation is set to take up to twenty iterations of
simulations to fit the model by estimating parameters, it will stop at
fewer if the generated networks converge on estimates of our coefficient
values or paramters. ERGMs are used to predict ties as a
function of individual covariates (i.e. attribute data, like “EX”
in our example) or network structure.
We will start with a very basic model, looking only at the
probability of edge formation, otherwise known as density, to
demonstrate how co-efficients can be translated into odds ratios.4 The term for tie density (edges) is often used similarly to an intercept term in a linear regression or other linear model such as R’s glm.
Keep in mind that the base network is about information retrieval. Our
model will primarily allow us to ask what the probability is that an
information retrieval relationship will form between two nodes.
model0<-ergm(CRIeq~edges)
## Evaluating log-likelihood at the estimate.
summary(model0)
##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: CRIeq ~ edges
##
## Iterations: 5 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % p-value
## edges -1.7288 0.1695 0 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 377.1 on 272 degrees of freedom
## Residual Deviance: 230.6 on 271 degrees of freedom
##
## AIC: 232.6 BIC: 236.3 (Smaller is better.)
What do coefficients mean? Coefficients are the change in the
(log-odds) likelihood of a tie for a unit change in a predictor. In our
basic model above, our only predictor is information on the number of
edges in the network. We can see that the coefficient estimate is
negative, suggesting that a tie is more likely not to form than form
(i.e. density is less than .5.) To get a better sense of how less
likely it is for a tie to form, we can translate our log-odds into a
probability.
To translate our estimated coefficient for the edges
parameter into a probability, we can take the inverse-logit. We are
thereby finding the probability that a tie will form in the network,
looking only at the number of ties in the base network to build our
model. We can see below that this probability is equal to our network
density.
plogis(coef(model0)[[1]])
## [1] 0.1507353
network.density(CRIeq)
## [1] 0.1507353
Conceptually, this should be fairly easy to follow: if all we know
about a network is the number of ties we have and we’re attempting to
predict the probability that an edge will exist solely based on that
information, then the probability of an edge existing at any point in
the network equals the density of the network as a whole.5 As we move on from this toy model, keep in mind that the estimate for our edges
parameter will change as we add additional predictors or network
statistics as additional terms will partially explain tie formation.
Base Network ERGM
As we build a network, we can evaluate whether individual network
statistics or node attributes prove our hypotheses and whether they do
so in a way that is significantly different from random chance. Later,
we will evaluate whether the model does a good job of explaining our
observed network.
model1<-ergm(CRIeq~edges# Set the base term based on density/edge formation. +mutual# H1+transitive# H2a: Transitive triads ( type 120D, 030T, 120U, or 300) # What if the tie is part of a transitive triad?+nodeicov("EX")# H3a+nodeocov("EX")# H3b+edgecov(CAIeq)# H4)
## Starting maximum likelihood estimation via MCMLE:
## Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
## This model was fit using MCMC. To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.
Take a look at the ERGM equation discussed in the week 5 slides, reproduced below:
Pr(Y=y)=exp[θ′g(y)]/k(θ)
What term in the equation do the ERGM terms or network statistics correspond to? ANSWER: g(y) is a vector of network statistics.
Explain each ERGM term and its relationship to your
hypotheses. Test your hypotheses (think about whether the sign of your
coefficient suggests a tie is more or less likely) and report whether
your results are significant.6 A parameter is significant if its absolute value is more than twice its Standard Error.
Hypothesis 1: Individuals are less likely to retrieve information from those who retrieve information from them. ANSWER:
According to the documentaion, “mutual” is “equal to the number of
pairs of actors i and j for which (i,j) and (j,i) both exist.” That
number is negative (-1.4515), so the estimate alone would confirm the
hypothesis. However, the std. error is relatively high (0.9915), which
means we can’t say that conclusively.
Hypothesis 2a: Information retrieval tends to be transitive. That is, if individual i retrieves information from individual k, and individual k retrieves information from individual j, then individual i is more likely to retrieve information from individual j. ANSWER:
“Transitive”, which according to the documentation is “the number of
triads in the network that are transitive” is positive (0.2817) and the
standard error is just small enough so that value is significant.
Hypothesis 3a: Individuals tend to retrieve information from other members with high expertise. ANSWER:
According to the documentation, nodeicov calculates the “main effect of
a covariate for in-edges”, in this case, for EX, which is the team
member’s expertise. This estimate is very positive and much high er than
the standard error, meaning the data supports the hypothesis with quite
a large certainty.
Hypothesis 3b: Individuals with low expertise tend to retrieve information from many others. ANSWER:
According to the documentation, nodeocov calculates the “main effect of
a covariate for out-edges”, also in this case, for EX, which is the
team member’s expertise. The estimate value for this is positive, but
the standard error is even larger than the estimate. Meaning the data
does not support this hypothesis.
Hypothesis 4: Individuals tend to retrieve information from members to whom they allocate information. ANSWER:
The funcion “edgecov”, which in this case uses data from CAIeq, has a
positive value much larger than its standard error. So the hypothesis is
strongly supported.
## Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
## This model was fit using MCMC. To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.
Evaluate the remaining hypothesis with your model.
Hypothesis 2b: Transitivity increases at a sub-linear rate as a function of the number of ties. ANSWER:
Again, another case in which the estimate points to a direction of the
hypothesis, but the standard error does not allow for support of the
hypothesis. gwesp has a estimate of 0.6104, but a std. error of 0.3613.
Model Diagnostics
Next, judge convergence of the MCMC processes of the first model, using the mcmc.diagnostics() function. The function will plot the change of model statistics during the last iteration of the MCMC estimation procedure.7
Note that although the edge graphs appear to be periodic, the dips
between whole numbers are due to the fact that edges are always whole
numbers. For each model statistic, the left hand side plot gives
the change of the statistic with iterations, and the right hand side
plot is a histogram of the statistic values. Both are normalized, so the
observed data locate at 0.
mcmc.diagnostics(model1)# Performs the markov chain monte carlo diagnostics
## Warning in formals(fun): argument is not a function
MCMC Diagnostics, Model 1.
MCMC Diagnostics, Model 1.
##
## MCMC diagnostics shown here are from the last round of simulation, prior to computation of final parameter estimates. Because the final estimates are refinements of those used for this simulation run, these diagnostics may understate model performance. To directly assess the performance of the final model on in-model statistics, please use the GOF command: gof(ergmFitObject, GOF=~model).
Repeat the process for the second model.
mcmc.diagnostics(model2)# Performs the markov chain monte carlo diagnostics
## Warning in formals(fun): argument is not a function
MCMC Diagnostics, Model 2.
MCMC Diagnostics, Model 2.
##
## MCMC diagnostics shown here are from the last round of simulation, prior to computation of final parameter estimates. Because the final estimates are refinements of those used for this simulation run, these diagnostics may understate model performance. To directly assess the performance of the final model on in-model statistics, please use the GOF command: gof(ergmFitObject, GOF=~model).
Analysis
Has the MCMC process converged to a desired state for each ERGM term? Explain how you interpret the plots.8 D. Lusher et al, ERGM,
12.3 (“[T]he inferential goal is to center the distribution of
statistics over those of the observed network, thus fitting a model that
we say gives maximal support to the data.”) ANSWER: Yes, it converged to the desired state. The charts show a distribution around the terms’ values.
Model Evaluation
To evaluate the goodness-of-fit for our model, we need to simulate many variations of the model.9See?? simulate for more information.
Let’s visually inspect two of our random networks based on our first model.
Next we’re going to extract the number of triangles from each of the
100 samples, create a histogram of that model, and place a red arrow at
the value of the observed network.
model.tridist<-sapply(1:100, function(x)summary(sim[[x]]~triangle))# Extracts the tiangle data from the simulated networkshist(model.tridist, xlim=c(0,140), breaks=20)# Plots that triangle distribution as a histogramCRIeq.tri<-summary(CRIeq~triangle)# Saves the CRIeq triangle data from the summary to the CRI.eq variablearrows(CRIeq.tri,20, CRIeq.tri, 5, col="red", lwd=3)# Adds an arrow to the plotted histogram
Triangle Distribution
Analysis
Is the distribution of triangles in your simulation a good match with the distribution of triangles in your observed network? ANSWER:
The arrow points to a place in the histogram which is close, but not
exactly, the most frequent distribution. So it seems close, but not
exactly accurate.
Goodness of Fit
Next, we will calculate the Goodness of Fit for several of the
parameters in our model. A p-value closer to one is better; this
represents the difference between the observed networks and simulations.
# -------------------------------------------------------------------------------------------------# Test the goodness of fit of the model# Compiles statistics for these simulations as well as the observed network, and calculates p-values # -------------------------------------------------------------------------------------------------par(mfrow=c(2,2))# Separate the plot window into a 2 by 2 orientationplot(gof)# Plot the goodness of fit
Goodness of Fit
Analysis
Evaluate the plots and summary statitistics of the Goodness
of Fit measures for Model 1. Are the four terms evaluated show a good
fit between the simulated networks and the observed network?10 In general, for configurations in the model, the fit is considered good if │t│≤ 0.1. For configurations not included in the model, the fit is considered good if 0.1<│t│≤ 1, and not extreme if 1< │t│≤
2. For your plot, the dark black line represents the data for the
observed network. The boxplots represent the distribution of
corresponding degrees across the simulated networks, and the soft lines
are the 95% confidence intervals. ANSWER: Out of the
four models, two of them show a good fit: in degree (with the exception
of 5 and 11) and edgewise shared partner. For the other two, the results
are not that good. For out-degree and minimum geodesic distance, almost
half do not fit.
Submitting the Lab (5 pts)
After knitting your file to RPubs, copy the URL and paste it into the
comment field of the Lab 2 Assignment on Canvas. Save this .Rmd file
and submit it in the file portion of your Canvas assignment. Make sure
to review your file and its formatting. Run spell check (built into RStudio) and proofread your answers before submitting. If you can’t publish to RPubs, save your HTML file as a PDF and submit that instead.11 There are many different ways to do this with different browsers. Google it.
