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.
This is a network with 17 nodes and 41 edges, representing the information retrieval network. It is a directed network, with each row representing one edge and the two columns representing the sender and the receiver. For instance, node 15 goes to node 3 to retrieve information; node 14 goes node 4 to retrieve information. It is a directed network, and the network density is 0.1507353. Our attribute information becomes the vertex attribute in the base network. The other network — allocation network — has the same nodes/actors as our base network, but it has 23 edges instead.
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.
The two visualizations sized by expertise value and in-degree centrality support hypothesis 3a. In the expertise value visualization, nodes 4, 5, 6, 9, 16, 17 are the ones who have high expertise, and most people go to retrieve information from them. Similarly in the in-degree centrality visualization, these nodes are larger, meaning that they have high in-degree centrality. Therefore, the graphs support the hypothesis that individuals tend to retrieve information from other members with high expertise.
Understanding the Covariate Network
Let’s explore the summary statistics of our co-variate network.
In your own words, explain what this network respresents and its relationship to the other network and the attribute information.
This network has 17 nodes and 23 edges, representing the allocation of information. Compared to the previous information retrieval network, information is being passed from nodes with low expertise to nodes with high expertise (we know for example 4, 6, 8, 9, 16, 17 have high expertise from the previous analysis).
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.
The retrieval and allocation networks pretty much support hypothesis 4. We can tell from the allocation network that most of the nodes that individual allocate information to are the same ones that they retrieve information from. However, there is a few exceptions: for example, individual allocates information to node 11 but no one retrieve information from node 11. This might be explained by the attribute since node 11 is with low expertise.
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[\theta'g(y)]/k(\theta)
\]
What term in the equation do the ERGM terms or network statistics correspond to?
g(y).
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.
Y is the network realization, similar to a random variable. y is the observed network. Pr(Y=y) means the probability of the generated network is the same as the observed network. g(y) is a vector of network statistics. θ’ is a vector of coefficients corresponding to g(y). k(θ) is a normalizing factor calculated by summing up the exponential value over all possible network configurations.
Hypothesis 1: Individuals are less likely to retrieve information from those who retrieve information from them.
The coefficient for “mutual” is negative, suggesting that a tie is less likely. However, p=0.14 means that it is not statistically significant. So we do not support hypothesis 1.
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.
The coefficient for “transitive” is positive, suggesting that a tie is more likely. p=0.02 is significant at 0.05, meaning that we can support hypothesis 2a at 95% confident level.
Hypothesis 3a: Individuals tend to retrieve information from other members with high expertise.
The coefficient for “nodeicov” is positive, suggesting that a tie is more likely. Also, p is less then 0.001, meaning that we can support hypothesis 3a with almost 100% confidence.
Hypothesis 3b: Individuals with low expertise tend to retrieve information from many others.
The coefficient for “nodeocov” is positive, suggesting that a tie is more likely. However, p=0.43 means that it is not statistically significant. So we do not support hypothese 3b.
Hypothesis 4: Individuals tend to retrieve information from members to whom they allocate information.
The coefficient for “edgecov” is positive, suggesting that a tie is more likely. Also, p=0.002 means that we can support hypothesis 4 at 99% confident level.
## 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.
The coefficient is positive, suggesting a tie is more likely. However, p=0.09 means that it is not statistically significant. So we do not support hypothesis 2b.
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.”)
Yes. The MCMC process converged to around center for each ERGM terms, which is at a desired state. We can tell from the plots that they all sort of converged at the center.
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 networks
hist(model.tridist, xlim=c(0,140), breaks = 20) # Plots that triangle distribution as a histogram
CRIeq.tri <- summary(CRIeq ~ triangle) # Saves the CRIeq triangle data from the summary to the CRI.eq variable
arrows(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?
The distribution of transitive triangles in the similation is different than the triangles in the observed network.Since the histogram represents the grpahs, we can tell from the histogram that the distribution looks different in the simulation compared to the observed network. It is not like a normal distribution.
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 orientation
plot(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.
Black dots and lines are observed network; grey dots and lines denote simulated networks. If the model is a good fit, black line should fall between or close enough to the grey lines. Therefore, the simulated network shows a good fit between the observed network because we can tell from the four GOF graphs that the black lines all fall in the 95% confidence intervals.
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.
