• 1. PREPARE
  • 2. WRANGLE
  • 3. EXPLORE
  • 4. MODEL
  • 5. COMMUNICATE
    • 👉 Your Turn ⤵
    • Data Visualization or Table
    • Narrative
    • 🧶 Knit & Check ✅
    • References

1. PREPARE

For the Unit 4 Case Study: Birds of Feather Lead Together, we once again visit the work of Alan Daly and colleagues as we attempt to replicate some of the analyses described in Chapter 9: Network Data and Statistical Models from Social Network Analysis and Education (Carolan 2014). In this case study move beyond the visual and mathematical descriptions of networks explored so far to uncover generative processes, or mechanisms, that attempt to explain how school leaders select peers for collaboration or confidential exchanges.

More specifically, this case study will cover the following topics pertaining to each data-intensive workflow process:

  1. Prepare: Prior to analysis, we’ll take a look at the context from which our data came, formulate some research questions, and get introduced the {statnet} R package for ERGMs.

  2. Wrangle: In section 2 we will learn how to work with matrices, how to dichotomize matrices, and how to convert matrices into more workable formats like our familiar edge-list.

  3. Explore: In section 3, we once again use the {tidygraph} package and the companion {ggraph} package to calculate basic summary stats for our network and will try to replicate and ideally improve upon the sociogram from our course text.

  4. Model: We wrap up our analysis in Section 4 by introducing exponential random graph modules used in Chapter 9 of Carolan (2014), learn to check how well our model “fits our data” and examine diagnostics that might indicate issues with our model.

  5. Communicate: Finally, we’ll select, polish, and narrate a key finding from our analysis that helps to answer or research question(s).

1a. Review the Research

In Social Network Analysis and Education: Theory, Methods & Applications Carolan (2014) makes the following distinctions between mathematical and statistical approaches to social network analysis:

  1. Mathematical approaches focus on what a network of actors “looks like” by describing the network using sociograms and/or network measures such as reciprocity, centrality, and density. However, these approaches “tend to regard the measured relationships and their strengths as accurately reflecting the real, final, or equilibrium status of the network.”

  2. Statistical approaches, or statistical inference, on the other hand, focuses on “assessing the reproducibility or likelihood of an observed pattern” with the goal of explaining and ultimately predicting network structures and outcomes.

For Unit 4, we move beyond previous “mathematical approaches” used so far and explore the use of ERGMs as a “statistical approach” to network analysis. Specifically, we will try to replicate the P* models used in Chapter 9 of Carolan (2014) in order to test hypotheses about individual and network characteristics that may explain who school leaders select for collaboration and confidential exchanges.

Research Questions

For this case study, we will be using school leadership data were collected at two school districts over 3 consecutive years that was used to answer the following research questions below:

  1. Is there a relationship between the frequency of collaboration between school leaders and how often they turn to each other to discuss issues of a confidential nature?
  2. Do school leaders prefer to collaborate with those with whom they have collaborated in the past, or is there some other reason?
  3. Does gender or some other individual attribute predicts confidential exchanges between school leaders, or does some previous relation have a stronger effect?
  4. Does collaboration between leaders explain one’s level of trust in one’s administrative colleagues?
  5. Can we distinguish among different groups of school leaders based on how frequently they collaborate, and if so, are these groupings related to the level at which they work (school versus district)?

Of specific interest for our case study is research question 3 emphasized in bold text above, which asks whether individual attributes predicts tie formation. For example, are males more likely to engage in exchange confidential exchanges with other.

The tendency for individuals to be attracted to and interact with others who share similar characteristics is formally referred to as homophily or assortativity, but more commonly known by the “birds of a feather flock together” proverb.

Homophily has been well documented in education research across a variety of relation types including student friendships, teacher-pupil interactions, and membership in student cohorts (McPherson and Smith-Lovin 1987). Studies of homophily among student friendships, for example, date back to the 1920’s and provide strong evidence that students form ties based on sociodemographic similarities.

Data Collection

For each consecutive year, school district leaders were invited to complete a survey that collected individual:

  • Demographic information (e.g., gender, ethnicity, marital status, age, years of experiences);

  • Network relationships types (e.g., collaboration, confidential exchanges, energy, expertise, leaders approached for support, work-related issues, input, recognition, best practices, and innovation);

  • Frequency of interactions they have with those nominated individuals on a four-point frequency scale ranging from 1 (the least frequent) to 4 (1–2 times a week);

  • Leadership efficacy items were designed based on the Principal Efficacy Scale used in Daly et al. (2011) and Tschannen-Moran and Gareis’s (2004) studies. The efficacy scale includes 18 items rated on a 9-point Likert scale ranging from 1 (None at all) to 9 (A great deal);

  • Trust scale contains eight items rated on a 7-point Likert scale ranging from 1 (Strongly disagree) to 7 (Strongly agree) modified from Tschannen-Moran and Hoy (2003).

Analyses

Two modeling techniques unique to network analysis and demonstrated in Chapter 9 of Carolan (2014) are the quadratic assignment procedure (QAP) and Exponential Random Graph Models (ERGM):

  • QAP/MR: The quadratic assignment procedure (QAP) developed by Hubert (1987) and Krackhardt (1987b) tests the null hypothesis of no correlation between the two networks and adjusts for this dependence between networks by repeatedly permuting the order of rows and columns of one of the networks while keeping the other network intact. The QAP is based on regression models and permutation tests for valued (i.e., continuous) relational variables.

  • P1 and P* (P-Star): Both of these models are some of the first to make use of the ERGM, which provides a basis for comparing whether a network’s observed structural properties occur more frequently than you could expect from chance alone. ERGMs can be used to model ties in complete networks but do so in a manner that explains the presence (or absence) of ties as a function of individual-and/or network-level characteristics. While QAP and MR-QAP procedures control for network structure through permutations, these ERGMs attempt to explain it.

Key Findings

In response to research questions 1 & 2 from above, Carolyn reported that “while collaboration in year 1 does not significantly predict collaboration in year 3, confidential exchanges in year 1 does” and suggested that “collaboration among school leaders provides an important foundation for more sensitive, perhaps even deeper, relations (e.g., confidential exchanges) at a later point in time.”

👉 Your Turn ⤵

For our Unit 4 case study, we focus specifically on generative processes, or micro-level mechanisms, the might help explain who school leaders select to collaborate with or engage in confidential exchanges. Review Chapter 9 of SNA and Education and report out on the findings from the ERGM analysis in response to Research Question 3 below:

  1. Does gender or some other individual attribute predicts confidential exchanges between school leaders, or does some previous relation have a stronger effect?

    • A previous relation has a stronger effect to predict confidential exchanges. Only efficacy attribute showed “that school leaders who believe that they have the capacity to have an effect are more likely to send and receive confidential exchange ties” (p11).

Based on what you know about networks and the context so far, write another research question or hypothesis that can be tested using ERGMs:

  • Do males have more reciprocated ties than females?

Later you’ll have a chance to empirically test it by the end of this case study.

1b. Load Libraries

Recall that packages, or libraries, are shareable collections of R code that can contain functions, data, and/or documentation and extend the functionality of R. You can always check to see which packages have already been installed and loaded into RStudio Cloud by looking at the the Files, Plots, & Packages Pane in the lower right hand corner.

statnet 📦

Similar to the collection of packages contained in the {tidyverse} package, the Statnet Team (Pavel N. Krivitsky et al. 2003-2020) has developed a suite of R packages for the management, exploration, statistical analysis, simulation and vizualization of network data. The statistical modeling framework used in {statnet} relies on Exponential-family Random Graph Models (ERGMs).

As noted in the statnet tutorial by the same authors, the exponential-family random graph models (ERGMs) are a general class of models based in exponential-family theory for specifying the probability distribution for a set of random graphs or networks. Within this framework, one can—among other tasks:

  1. Define a model for a network that includes covariates representing features like homophily, mutuality, triad effects, and a wide range of other structural features of interest;

  2. Obtain maximum-likehood estimates for the parameters of the specified model for a given data set;

  3. Test individual coefficients, assess models for convergence and goodness-of-fit and perform various types of model comparison; and

  4. Simulate new networks from the underlying probability distribution implied by the fitted model.

👉 Your Turn ⤵

Let’s load the {statnet} package that we’ll be using in the Model section of our case study for bullets 1-3 above, as well as the following packages that will be using for our network analysis:

  • tidyverse
  • readxl
  • igraph
  • tidygraph
  • ggraph
  • skimr
  • janitor
library(knitr)
library(tidyverse)
library(janitor)
library(skimr)
library(ggraph)
library(igraph)
library(tidygraph)
library(readxl)
library(network)
library(ergm)
library(kableExtra)
library(texreg)

2. WRANGLE

For our data wrangling this week, we’ll focus on working with network data stored as an adjacency matrix. Our primary goals for this section are learning how to:

  1. Import Data. In this section, we use the {readxl} package from the tidyverse to read in a matrix and node attributes.

  2. Dichotomize a Matrix. As described in Chapter 9, we’ll recode our edge values to 1s and 0s changing our valued matrix to a binary matrix.

  3. Create Network Graph. Finally, we’ll convert our matrix to an edge-list and store both our edges and node attributes as a network igraph object in preparation for analysis.

2a. Import Data

One of our primary goals is to replicate the ERGM analysis from Chapter 9: Network Data and Statistical Models (Carolan 2014). Specifically, we’ll aim to reproduce the results from Table 9.4. Results of P1 and P*Analyses for the Dichotomized and Directed School Leaders Confidential Exchanges Network Year 3.

To do so, we’ll need to import two Excel files from the Social Network Analysis and Education companion site. The first file contains our edges stored as a square matrix (more on this later) and the second file is a standard rectangular data frame that contains attributes for each node, i.e. school leaders. These files are included in the data folder of your R Studio project. A description of each file from the companion website is linked above and each data file is linked below:

  1. School Leaders Data Chapter 9_d. This adjacency matrix reports on “confidential help” ties among 43 school leaders in year 3 of a three-year study. This is a directed valued (weighted) network measured on five-point scale ranging from 0 to 4, with higher values indicating more frequent collaborations (1–2 times/week). These data are used throughout Chapter 9.

  2. School Leaders Data Chapter 9_e. This rectangular matrix consists of four attribute vectors for 43 school leaders. Following the first ID column, the matrix includes an efficacy score, trust score, and indicators for whether one works at the district-level and is male (1 = yes, 0 = no). These attribute variables will be used as covariates in our ERGMs later in this walkthrough.

Since we are working with Excel files, we’ll need to use the read_excel() function from {readxl} tidyverse package to import our data. Let’s import the School Leaders Data Chapter 9_e.xlsx node file located in the data/ folder first:

leader_nodes <- read_excel("data/School Leaders Data Chapter 9_e.xlsx", 
                           col_types = c("text", "numeric", "numeric", "numeric", "numeric")) |>
  clean_names()

leader_nodes
## # A tibble: 43 × 5
##    id    efficacy trust district_site  male
##    <chr>    <dbl> <dbl>         <dbl> <dbl>
##  1 1         6.06  4                1     0
##  2 2         6.56  5.63             1     0
##  3 3         7.39  4.63             1     0
##  4 4         4.89  4                1     0
##  5 5         6.06  5.75             0     1
##  6 6         7.39  4.38             0     0
##  7 7         5.56  3.63             0     1
##  8 8         7.5   5.63             1     1
##  9 9         7.67  5.25             0     0
## 10 10        6.64  4.78             0     0
## # … with 33 more rows

Note that we specified the ID column as “text” and the remaining columns as “numeric.” By default, the readxl() function would have recognized each column as numeric, but the first column indicates the “names” for each school leader and we’ll be using this column to assign names to our columns and rows for our adjacency matrix in just a bit.

👉 Your Turn ⤵

Recall from above that our relations, or edges, are stored as a valued adjacency matrix in which columns and rows consist of the same actors and each cell contains information about the tie between each pair of actors. In our case, the tie is a directed and valued “arc” where the value indicates the strength of the relationship.

Use the read_excel() function to import the School Leaders Data Chapter 9_d.xlsx file, add an argument setting the column names to FALSE since our file is a simple matrix with no header or column names, and assign the matrix to a variable named leader_matrix:

Hint: Type ?read_excel into the console and check the arguments section to find the name of the argument used to set column headers.

leader_matrix <-  read_excel("data/School Leaders Data Chapter 9_d.xlsx", col_names = FALSE)

Use the code chunk below to inspect the matrix you just imported and answer the questions that follow:

leader_matrix
## # A tibble: 43 × 43
##     ...1  ...2  ...3  ...4  ...5  ...6  ...7  ...8  ...9 ...10 ...11 ...12 ...13
##    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1     0     0     0     0     0     0     0     0     0     0     0     0     0
##  2     0     0     0     0     0     0     0     0     0     0     0     0     0
##  3     3     1     0     1     1     1     1     2     1     1     1     1     1
##  4     1     0     0     0     0     0     0     0     0     0     0     0     0
##  5     0     0     0     0     0     0     1     0     3     0     0     0     0
##  6     0     0     0     0     0     0     0     0     0     0     0     0     0
##  7     2     1     1     1     1     1     4     1     1     3     1     1     1
##  8     1     1     4     1     1     1     1     0     2     3     1     1     1
##  9     1     0     0     0     4     0     0     1     0     0     0     0     0
## 10     1     1     3     1     1     1     4     4     1     0     3     1     1
## # … with 33 more rows, and 30 more variables: ...14 <dbl>, ...15 <dbl>,
## #   ...16 <dbl>, ...17 <dbl>, ...18 <dbl>, ...19 <dbl>, ...20 <dbl>,
## #   ...21 <dbl>, ...22 <dbl>, ...23 <dbl>, ...24 <dbl>, ...25 <dbl>,
## #   ...26 <dbl>, ...27 <dbl>, ...28 <dbl>, ...29 <dbl>, ...30 <dbl>,
## #   ...31 <dbl>, ...32 <dbl>, ...33 <dbl>, ...34 <dbl>, ...35 <dbl>,
## #   ...36 <dbl>, ...37 <dbl>, ...38 <dbl>, ...39 <dbl>, ...40 <dbl>,
## #   ...41 <dbl>, ...42 <dbl>, ...43 <dbl>

  1. Is our network a square matrix? What does this tell you about the number of nodes in our network?

    • Yes, this is a square matrix because it has the same number of columns as rows.

  2. Does our network contain self-loops? How do you know?

    • Yes, because there are nodes that are connected to themselves.

  3. What do the values in each cell of our matrix indicate?

    • The value indicates the strength of the relationship.

2b. Dichotomize Matrix

Since one of our goals is a crude replication of the analysis demonstrated in Chapter 9, and since working with “valued networks” is a little more complex than binary networks that simply indicate the presence or absence of a tie, we will need to “dichotomize” our matrix. In practice this simply entails converting our valued matrix to a binary matrix with cells containing just 1s or 0s.

In their article, Valued Ties Tell Fewer Lies: Why Not To Dichotomize Network Edges With Thresholds, Thomas and Blitzstein (2018) highlight several reasons why the dichotomization procedure is appealing in an investigation aside from convenience and simplicity:

  1. Use of Exclusively Binary Methods. Several classes of models have been designed to incorporate binary information directly, including the exponential random graph model.

  2. Ease of Input and Data Collection. The need to classify continuously-valued quantities into a set of discrete groups is widespread throughout all of science and technology, particularly because of the associated need to make clear decisions based on this information.

  3. Ease of Output in Graphical Representations. The visual appeal of graphs and networks has contributed to much of the field’s attention in the past decade. When plotting a graphical structure, a clever choice of threshold can illuminate which nodes are most central, which connections the most vital.

  4. Sparsity of Structure. In data where there are very few natural zeroes (if any), di- chotomization provides a way to select for a small number of connections which are thought to be of the greatest importance to the system, or to nominate a number of ties for more in-depth study.

  5. Binning To Address Nonlinearity and Reduce Noise. If there is a nonlinear relationship in the data, binning the data into distinct ordinal categories has many advantages, namely, the reduction of total mean-squared error, and a corresponding increase in power for detecting a true non-zero relationship over an improperly specified linear analysis.

As the title of their article suggests, however, Thomas and Blitzstein (2011) argue that the motivations for dichotomization should be revisited as dichotomization produces a range of problematic issues.

With that in mind, and before we can dichotomize our “matrix,” we first need to convert it to a matrix object recognized by R using the as.matrix() function. We can then check to data format of our leader_matrix using the class() function:

leader_matrix <- leader_matrix |>
  as.matrix()

class(leader_matrix)
## [1] "matrix" "array"

Both the collaboration and confidential help network data were dichotomized by recoding values originally coded as a 3 and 4 recoded to 1, indicating the presence of a directed tie for both relations, and zero otherwise.

To dichotomize our our matrix, the following code will “assign” 0’s to all values in our matrix that are less than or equal to 2, and 1’s to all values that are greater the or equal to 3:

leader_matrix[leader_matrix <= 2] <- 0

leader_matrix[leader_matrix >= 3] <- 1

Add Row and Column Names

Before we can convert to an edge-list, we will also need to add the names of our nodes to the columns and rows of our matrix. These are stored in the ID column of our leader_node data frame. We can use the $ operator to select these names and assign to our leader_matrix using the rownames() and colnames() functions respectively:

rownames(leader_matrix) <- leader_nodes$id

colnames(leader_matrix) <- leader_nodes$id

leader_matrix
##    1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
## 1  0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 2  0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 3  1 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1
## 4  0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  1
## 5  0 0 0 0 0 0 0 0 1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 6  0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 7  0 0 0 0 0 0 1 0 0  1  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0
## 8  0 0 1 0 0 0 0 0 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  1  1
## 9  0 0 0 0 1 0 0 0 0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 10 0 0 1 0 0 0 1 1 0  0  1  0  0  0  0  0  0  1  0  1  0  0  0  0  0  0  0  0
## 11 0 0 0 0 0 0 0 0 0  1  0  0  0  0  0  0  0  0  0  0  1  0  1  0  0  0  0  0
## 12 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 13 0 0 0 0 0 0 0 1 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 14 0 0 0 0 0 0 0 0 0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 15 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0
## 16 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0
## 17 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  1  0  0
## 18 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 19 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0
## 20 0 0 0 0 0 0 1 0 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 21 1 0 0 0 0 0 0 0 0  0  1  1  0  0  1  0  0  0  0  0  1  0  0  0  0  0  0  0
## 22 0 0 0 0 0 0 0 1 0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  1
## 23 0 0 0 0 1 0 0 0 0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 24 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0
## 25 0 0 1 0 0 0 0 1 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 26 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0
## 27 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0
## 28 1 1 1 1 0 1 0 1 1  1  0  0  1  0  0  0  0  1  0  0  1  1  0  0  1  0  1  1
## 29 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 30 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 31 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0
## 32 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0
## 33 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0
## 34 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  1  0  0  0  1  0  0  0  0  0  1
## 35 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0
## 36 1 0 1 1 0 0 0 1 0  0  0  0  0  0  0  1  0  1  0  0  0  0  0  0  0  1  1  1
## 37 0 0 1 0 0 0 0 0 0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 38 1 0 1 0 0 0 0 1 0  0  1  0  0  0  0  0  0  1  0  0  0  1  0  0  0  1  1  1
## 39 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 40 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 41 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 42 0 0 0 0 0 0 0 0 1  0  0  0  0  0  1  0  0  0  0  0  0  0  0  1  0  0  1  0
## 43 0 0 0 0 1 0 0 1 0  0  0  0  0  0  1  0  0  1  0  0  1  1  1  0  0  0  0  0
##    29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
## 1   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 2   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 3   0  0  0  0  0  0  0  1  0  1  0  0  0  0  0
## 4   0  0  0  0  0  1  0  1  0  0  0  0  0  0  0
## 5   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 6   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 7   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 8   0  0  0  0  0  0  0  0  1  1  0  0  0  0  0
## 9   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 10  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0
## 11  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0
## 12  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 13  0  0  0  0  0  0  1  0  0  0  0  0  0  1  0
## 14  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0
## 15  0  0  0  0  0  0  0  0  0  0  1  0  0  1  1
## 16  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 17  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 18  1  0  0  0  0  1  0  0  0  0  0  0  0  0  0
## 19  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0
## 20  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 21  0  0  0  0  1  0  1  0  1  1  0  0  0  0  1
## 22  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 23  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 24  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0
## 25  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 26  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0
## 27  0  0  0  1  0  0  0  0  0  0  0  0  0  1  0
## 28  1  0  0  1  0  1  0  1  1  1  0  0  0  0  0
## 29  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 30  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 31  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 32  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 33  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 34  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 35  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 36  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0
## 37  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0
## 38  0  0  0  0  1  1  0  1  1  0  1  0  0  0  1
## 39  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 40  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 41  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0
## 42  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1
## 43  0  0  0  1  1  1  0  0  1  0  0  0  0  1  0

2c. Create Network with Attributes

Recall that edge-lists contain a row for each dyad consisting of at minimum two columns with the name of each actor, and which conveniently can also contain other information or attributes about the relationship such as edge weight, timestamps, or other contextual information as demonstrated in Unit 1.

Edge-lists also have the advantage of being easier to work with when using network packages in R.

Get Edges

The {igraph} package introduced in Unit 1 has a convenient get.data.frame() function for extracting an edge list from a matrix, but first we need to convert our matrix to an igraph network object.

adjacency_matrix <- graph.adjacency(leader_matrix,
                                    diag = FALSE)

class(adjacency_matrix)
## [1] "igraph"
adjacency_matrix
## IGRAPH 09dde33 DN-- 43 143 -- 
## + attr: name (v/c)
## + edges from 09dde33 (vertex names):
##  [1] 3 ->1  3 ->27 3 ->28 3 ->36 3 ->38 4 ->16 4 ->28 4 ->34 4 ->36 5 ->9 
## [11] 7 ->10 7 ->20 8 ->3  8 ->10 8 ->25 8 ->27 8 ->28 8 ->37 8 ->38 9 ->5 
## [21] 9 ->14 10->3  10->7  10->8  10->11 10->18 10->20 10->38 11->10 11->21
## [31] 11->23 11->35 13->8  13->35 13->42 14->37 15->21 15->39 15->42 15->43
## [41] 16->22 17->26 18->29 18->34 19->31 20->7  20->10 21->1  21->11 21->12
## [51] 21->15 21->33 21->35 21->37 21->38 21->43 22->8  22->16 22->28 23->5 
## [61] 23->11 24->27 24->42 25->3  25->8  26->17 26->41 27->24 27->32 27->42
## [71] 28->1  28->2  28->3  28->4  28->6  28->8  28->9  28->10 28->13 28->18
## + ... omitted several edges

Note that we included the diag = FALSE argument which converts all values along the diagonal to 0s, thereby removing self-loops. I believe they may have been included when calculating the network descriptives included in Table 9.1, and our own descriptives may not match exactly, but we’ll need to remove these from our network for the ERGM analysis so this will save us a step.

Now we can use the get.data.frame() function to covert our matrix to a standard edge-list:

leader_edges <- get.data.frame(adjacency_matrix) |>
  mutate(from = as.character(from)) |>
  mutate(to = as.character(to))


leader_edges
##     from to
## 1      3  1
## 2      3 27
## 3      3 28
## 4      3 36
## 5      3 38
## 6      4 16
## 7      4 28
## 8      4 34
## 9      4 36
## 10     5  9
## 11     7 10
## 12     7 20
## 13     8  3
## 14     8 10
## 15     8 25
## 16     8 27
## 17     8 28
## 18     8 37
## 19     8 38
## 20     9  5
## 21     9 14
## 22    10  3
## 23    10  7
## 24    10  8
## 25    10 11
## 26    10 18
## 27    10 20
## 28    10 38
## 29    11 10
## 30    11 21
## 31    11 23
## 32    11 35
## 33    13  8
## 34    13 35
## 35    13 42
## 36    14 37
## 37    15 21
## 38    15 39
## 39    15 42
## 40    15 43
## 41    16 22
## 42    17 26
## 43    18 29
## 44    18 34
## 45    19 31
## 46    20  7
## 47    20 10
## 48    21  1
## 49    21 11
## 50    21 12
## 51    21 15
## 52    21 33
## 53    21 35
## 54    21 37
## 55    21 38
## 56    21 43
## 57    22  8
## 58    22 16
## 59    22 28
## 60    23  5
## 61    23 11
## 62    24 27
## 63    24 42
## 64    25  3
## 65    25  8
## 66    26 17
## 67    26 41
## 68    27 24
## 69    27 32
## 70    27 42
## 71    28  1
## 72    28  2
## 73    28  3
## 74    28  4
## 75    28  6
## 76    28  8
## 77    28  9
## 78    28 10
## 79    28 13
## 80    28 18
## 81    28 21
## 82    28 22
## 83    28 25
## 84    28 27
## 85    28 29
## 86    28 32
## 87    28 34
## 88    28 36
## 89    28 37
## 90    28 38
## 91    31 19
## 92    32 24
## 93    33 21
## 94    34 18
## 95    34 22
## 96    34 28
## 97    34 29
## 98    35 21
## 99    36  1
## 100   36  3
## 101   36  4
## 102   36  8
## 103   36 16
## 104   36 18
## 105   36 26
## 106   36 27
## 107   36 28
## 108   36 34
## 109   37  3
## 110   37 14
## 111   37 38
## 112   38  1
## 113   38  3
## 114   38  8
## 115   38 11
## 116   38 18
## 117   38 22
## 118   38 26
## 119   38 27
## 120   38 28
## 121   38 33
## 122   38 34
## 123   38 36
## 124   38 37
## 125   38 39
## 126   38 43
## 127   42  9
## 128   42 15
## 129   42 24
## 130   42 27
## 131   42 43
## 132   43  5
## 133   43  8
## 134   43 15
## 135   43 18
## 136   43 21
## 137   43 22
## 138   43 23
## 139   43 32
## 140   43 33
## 141   43 34
## 142   43 37
## 143   43 42

👉 Your Turn ⤵

Recall from Unit 2 that we introduced the {tidygraph} package for preparing and summarizing our Twitter network. Tidygraph includes the full functionality of igraph in a tidy API giving you access to almost all of the dplyr verbs plus a few more, developed for use with relational data.

Similar to Unit 2, use the tbl_graph() function to convert our leader_edges and leader_nodes data frames into a network graph object, by including the following arguments and supplying the appropriate code:

  • edges = expects a data frame, in our case leader_edges, containing information about the edges in the graph. The nodes of each edge must either be in a to and from column, or in the two first columns like the data frame we provided.

  • nodes = expects a data frame, in our case leader_nodes, containing information about the nodes in the graph. If to and/or from are characters or names, like in our data frames, then they will be matched to the column named according to node_key in nodes, if it exists, or matched to the first column in the node list.

  • directed = specifies whether the constructed graph be directed.

# YOUR CODE HERE
leader_graph <- tbl_graph(edges = leader_edges,
                          nodes = leader_nodes,
                          directed = TRUE)


leader_graph
## # A tbl_graph: 43 nodes and 143 edges
## #
## # A directed simple graph with 4 components
## #
## # Node Data: 43 × 5 (active)
##   id    efficacy trust district_site  male
##   <chr>    <dbl> <dbl>         <dbl> <dbl>
## 1 1         6.06  4                1     0
## 2 2         6.56  5.63             1     0
## 3 3         7.39  4.63             1     0
## 4 4         4.89  4                1     0
## 5 5         6.06  5.75             0     1
## 6 6         7.39  4.38             0     0
## # … with 37 more rows
## #
## # Edge Data: 143 × 2
##    from    to
##   <int> <int>
## 1     3     1
## 2     3    27
## 3     3    28
## # … with 140 more rows

Congrats! You made it to the end of data wrangling section and are ready to start analysis!

3. EXPLORE

In Section 3, we use the {tidygraph} package for retrieving network descriptives and the {ggraph} package to create a network visualization to help illustrate these metrics. Specifically, in this section we will:

  1. Examine Basic Descriptives. We focus primarily on actors and edges in this walkthrough, including whether or not ties were reciprocated and node degree, an important and fairly intuitive measure of centrality.

  2. Make a Sociogram. Finally, we wrap up the explore phases by learning to plot a network and tweak key elements like the size, shape, and position of nodes and edges to better at communicating key findings.

3a. Examine Descriptives

As noted in SNA and Education (Carolan 2014), many analyses of social networks are primarily descriptive and aim to either represent the network’s underlying social structure through data-reduction techniques or to characterize network properties through network measures.

In the analyses described in Chapter 9, descriptives were limited to the mean and standard deviation for: in-degree, out-degree, trust and efficacy measures. The proportion of male leaders was also reported. For section 3a, let’s see if we can reproduce these descriptives.

Calculate Node Degree

Recall that degree is the number of ties to and from an ego, or in the case of a “simple graph” like ours, degree is simply the number of people to whom someone is connected. In a directed network, in-degree is the number of connections or ties received, whereas out-degree is the number of connections or ties sent.

The activate() function from the {tidygraph} package allows us to treat the nodes in our network object as if they were a standard data frame to which we can then apply tidyverse functions such as select(), filter(), and mutate().

We can use the mutate() functions to create new variables for nodes such as measures of degree, in-degree, and out-degree using the centrality_degree() function in the {tidygraph} package.

Run the following code to add in- and out-degree measures to each of our nodes and examine the output:

leader_measures <- leader_graph |>
  activate(nodes) |>
  mutate(in_degree = centrality_degree(mode = "in")) |>
  mutate(out_degree = centrality_degree(mode = "out"))

leader_measures
## # A tbl_graph: 43 nodes and 143 edges
## #
## # A directed simple graph with 4 components
## #
## # Node Data: 43 × 7 (active)
##   id    efficacy trust district_site  male in_degree out_degree
##   <chr>    <dbl> <dbl>         <dbl> <dbl>     <dbl>      <dbl>
## 1 1         6.06  4                1     0         5          0
## 2 2         6.56  5.63             1     0         1          0
## 3 3         7.39  4.63             1     0         7          5
## 4 4         4.89  4                1     0         2          4
## 5 5         6.06  5.75             0     1         3          1
## 6 6         7.39  4.38             0     0         1          0
## # … with 37 more rows
## #
## # Edge Data: 143 × 2
##    from    to
##   <int> <int>
## 1     3     1
## 2     3    27
## 3     3    28
## # … with 140 more rows

👉 Your Turn ⤵

We now see that, in addition to the previously included attributes trust and efficacy, in-degree and out-degree measures have been added to the nodes in our network. But what if we also want to know to total number of “alters” each “ego” is connected to, i.e. the total number of individuals each school leader are connected to?

Modify the code below to calculate degree for each school leader in our network. Hint: centrality_degree is a wrapper for igraph::degree() and the mode = argument can be found in the corresponding help documentation.

leader_measures <- leader_graph |>
  activate(nodes) |>
  mutate(in_degree = centrality_degree(mode = "in")) |>
  mutate(out_degree = centrality_degree(mode = "out"))

leader_measures
## # A tbl_graph: 43 nodes and 143 edges
## #
## # A directed simple graph with 4 components
## #
## # Node Data: 43 × 7 (active)
##   id    efficacy trust district_site  male in_degree out_degree
##   <chr>    <dbl> <dbl>         <dbl> <dbl>     <dbl>      <dbl>
## 1 1         6.06  4                1     0         5          0
## 2 2         6.56  5.63             1     0         1          0
## 3 3         7.39  4.63             1     0         7          5
## 4 4         4.89  4                1     0         2          4
## 5 5         6.06  5.75             0     1         3          1
## 6 6         7.39  4.38             0     0         1          0
## # … with 37 more rows
## #
## # Edge Data: 143 × 2
##    from    to
##   <int> <int>
## 1     3     1
## 2     3    27
## 3     3    28
## # … with 140 more rows

Summarize Node Measures

We can also use the activate() function combined with the as_tibble() function introduced in our previous unit to extract our new measures to a separate data frame we’ll call node_measures so we can inspect our nodes individually and later use to calculate some summary stats:

node_measures <- leader_measures |> 
  activate(nodes) |>
  as_tibble()

node_measures
## # A tibble: 43 × 7
##    id    efficacy trust district_site  male in_degree out_degree
##    <chr>    <dbl> <dbl>         <dbl> <dbl>     <dbl>      <dbl>
##  1 1         6.06  4                1     0         5          0
##  2 2         6.56  5.63             1     0         1          0
##  3 3         7.39  4.63             1     0         7          5
##  4 4         4.89  4                1     0         2          4
##  5 5         6.06  5.75             0     1         3          1
##  6 6         7.39  4.38             0     0         1          0
##  7 7         5.56  3.63             0     1         2          2
##  8 8         7.5   5.63             1     1         8          7
##  9 9         7.67  5.25             0     0         3          2
## 10 10        6.64  4.78             0     0         5          7
## # … with 33 more rows

Now let’s view some basic summary statistics for each of the variables using the handy summary() function included in the R {base} package and the skim() functions from the {skimr} package.

summary(node_measures)
##       id               efficacy         trust       district_site   
##  Length:43          Min.   :4.610   Min.   :3.630   Min.   :0.0000  
##  Class :character   1st Qu.:5.670   1st Qu.:4.130   1st Qu.:0.0000  
##  Mode  :character   Median :6.780   Median :4.780   Median :0.0000  
##                     Mean   :6.649   Mean   :4.783   Mean   :0.4186  
##                     3rd Qu.:7.470   3rd Qu.:5.440   3rd Qu.:1.0000  
##                     Max.   :8.500   Max.   :5.880   Max.   :1.0000  
##       male          in_degree       out_degree    
##  Min.   :0.0000   Min.   :0.000   Min.   : 0.000  
##  1st Qu.:0.0000   1st Qu.:2.000   1st Qu.: 1.000  
##  Median :0.0000   Median :3.000   Median : 2.000  
##  Mean   :0.4419   Mean   :3.326   Mean   : 3.326  
##  3rd Qu.:1.0000   3rd Qu.:5.000   3rd Qu.: 4.000  
##  Max.   :1.0000   Max.   :8.000   Max.   :20.000
skim(node_measures)
Data summary
Name node_measures
Number of rows 43
Number of columns 7
_______________________
Column type frequency:
character 1
numeric 6
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
id 0 1 1 2 0 43 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
efficacy 0 1 6.65 1.10 4.61 5.67 6.78 7.47 8.50 ▅▅▃▇▃
trust 0 1 4.78 0.71 3.63 4.13 4.78 5.44 5.88 ▆▆▅▆▇
district_site 0 1 0.42 0.50 0.00 0.00 0.00 1.00 1.00 ▇▁▁▁▆
male 0 1 0.44 0.50 0.00 0.00 0.00 1.00 1.00 ▇▁▁▁▆
in_degree 0 1 3.33 2.12 0.00 2.00 3.00 5.00 8.00 ▅▇▂▅▂
out_degree 0 1 3.33 4.26 0.00 1.00 2.00 4.00 20.00 ▇▁▁▁▁

It looks like our summary stats for Year 3 confidential exchanges network and attribute data are pretty close to those reported in Table 9.1 copied from Carolan (2014). For example, our mean total efficacy and trust scores were 6.65 and 4.78 respectively, compared to 6.64 and 4.77 as reported in Chapter 9.

The average in/out-degree for our school leaders data is 3.26, however, which is a little lower than that reported by Carolyn. Again, this is likely a function of our removal of self-loops. To test this theory, we could simply change the diag = FALSE argument added to the graph.adjacency() function above to TRUE to include self-loops and then try rerunning all the code above.

School/District-Level Stats

Since we are now working with a standard data frame, we can also apply {dyplr} functions like group_by() and summarise() to calculate basic summary stats such as counts, mean, and standard deviation as follows:

node_measures |>
  group_by(district_site) |>
  summarise(n = n(),
            mean = mean(in_degree), 
            sd = sd(in_degree)
            )
## # A tibble: 2 × 4
##   district_site     n  mean    sd
##           <dbl> <int> <dbl> <dbl>
## 1             0    25  2.36  1.58
## 2             1    18  4.67  2.09

We see that our measures are pretty close, but not an exact match. Again, this is likely due to the self-loops we excluded. For example, our average in-degree for district-level (coded “1”) and school-level (coded “0”) leaders is 4.66 and 2.36 respectively, but reported as 4.72 and 2.64 by Carolan (2014).

👉 Your Turn ⤵

Use the code chunk below and additional chunks if needed to try and replicate the school and district level findings for out-degree, trust and efficacy measures.

outdegree <- node_measures|>
  group_by(district_site) |>
  summarise(n = n(),
            mean = mean(out_degree), 
            sd = sd(out_degree)
            )
outdegree
## # A tibble: 2 × 4
##   district_site     n  mean    sd
##           <dbl> <int> <dbl> <dbl>
## 1             0    25  2.44  3.01
## 2             1    18  4.56  5.41
trust2 <- node_measures|>
  group_by(district_site) |>
  summarise(n = n(),
            mean = mean(trust), 
            sd = sd(trust)
            )
trust2
## # A tibble: 2 × 4
##   district_site     n  mean    sd
##           <dbl> <int> <dbl> <dbl>
## 1             0    25  4.90 0.719
## 2             1    18  4.62 0.688
efficacy <- node_measures|>
  group_by(district_site) |>
  summarise(n = n(),
            mean = mean(efficacy), 
            sd = sd(efficacy)
            )
efficacy
## # A tibble: 2 × 4
##   district_site     n  mean    sd
##           <dbl> <int> <dbl> <dbl>
## 1             0    25  7.02 0.953
## 2             1    18  6.13 1.12

How close your result to those reported by Carolyn Table 9.1?

  • Out degree if off at 2.44 & 4.44 compared to Carolan’s 2.72 and 4.61
  • Trust is correct to 4.90 for school level and <.01 off at 4.61 compared to Carolan’s 4.60
  • efficacyis also <.01 for both the school level and district level at 7.01 vs Carolan’s 7.02 and for the district level, I got 6.13 vs.6.12

3b. Visualize Network Properties

In Chapter 9, Carolan (2014) depicts a directed and dichotomous sociogram of the collaboration network for year 3 is shown in .

👉 Your Turn ⤵

Try creating a Year 3 Confidential Exchange Network by modifying the code below and tweaking the included function/arguments or adding new ones for layouts, nodes, and edges to make our plot either more “aesthetically pleasing” or more purposeful in what it’s trying to communicate.

#YOUR CODE HERE
leader_measures |>
ggraph(layout = "fr") + 
  geom_edge_link(aes(alpha = stat(index))) +
  geom_node_point(aes(size = efficacy, alpha = efficacy, color = district_site))

After you are satisfied with your sociogram, answer the following questions?

  1. How is the Year 3 confidential exchange network similar to the collaboration network? How is it different?

    • Both graphs have 43 node and there is a strong tie strength between both graphs. However because we took the loops out the in-degree is lower then the collaboration graph. Additionally, there is a outlier subgroup.

Congrats! You made it to the end of the Explore section and are ready to learn a little ab out modeling network selection processes using ERGMs! Before proceeding further, knit your document and check to see if you encounter any errors.

4. MODEL

Recall from Carolan (2014) that network models and the questions they are designed to address are categorized into three different analytical emphases:

  1. Relationship-level models that focus on the ties between actors in complete networks;

  2. Models that predict individual actors’ attributes; and,

  3. Actor-level models that emphasize the differences within and among groups of actors within a complete network.

Exponential Random Graph Models fall into category 1 and are relationship-level models that focus on how formation of dyadic ties within a complete relational network can be explained by both structural and attribute variables. ERGMs provide a means to compare whether a network’s observed structural properties occur more frequently than you could expect from chance alone. More specifically, ERGMS provide a way to determine whether observed network properties like reciprocity and assortativity occur by chance as a result of other network properties.

While the technical aspects of estimating ERGMs are complex, their interpretation is pretty straightforward. In this section we are interested exploring how these models can be used to make inferences about the social processes at work in the School Leaders data. For example, in this section we will explore the following questions posed in Chapter 9:

  1. Reciprocity. If one school leader turns to another to discuss something confidential, is the latter likely to reciprocate?

  2. Transitivity. Are school leaders more likely to confide in someone if they both confide in the same school leader?

  3. Attribute Effect. Do a leader’s gender and efficacy score, predict a confidential exchange between two leaders?

4a. Loading Network Data

As we’ve discovered in this course, network data can come in many different forms — ties can be stored as edgelists or matrices, saved as .csv or excel files, and converted to a variety of network R objects. Attributes for the nodes, ties and dyads in can also various forms and added in various ways to objects in R. For the {ergm} package, however, data will need to be transformed into a object of the class “network” — the format that statnet uses to store and work with network data.

Fortunately, we have already prepared our data for quick conversion to a network object and can supply just a couple arguments to the as.network() function from the {network} package used by statnet:

  • x = a matrix giving the network structure in adjacency, incidence, or edgelist form;

  • vertices = an optionaldata.frame containing the vertex attributes. The first column is assigned to the "vertex.names" and additional columns are used to set vertex attributes using their column names.

Similar to the tbl_graph function from {tidygraph} that we used above, let’s add our leader_edges edge-list and our leader_nodes data frame as the x = and vertices = arguments in the as.network() function and assign to a new object called leader_network:

leader_network <- as.network(leader_edges,
                             vertices = leader_nodes)

leader_network
##  Network attributes:
##   vertices = 43 
##   directed = TRUE 
##   hyper = FALSE 
##   loops = FALSE 
##   multiple = FALSE 
##   bipartite = FALSE 
##   total edges= 143 
##     missing edges= 0 
##     non-missing edges= 143 
## 
##  Vertex attribute names: 
##     district_site efficacy male trust vertex.names 
## 
## No edge attributes

Note that we could also have supplied our leader_matrix as the first argument, but using and edgelist vastly simplifies the process of adding node attributes.

Let’s also check the be sure that our leader_network is indeed a network object by using the class() function to identify object type:

class(leader_network)
## [1] "network"

So far so good! Let’s start building our first ERGM model!

4b. Estimate the ERGM Model

We’ll begin by running a simple model both to demonstrate the most commonly used functions for ERGMs and also to replicate the approach in Chapter 9 of Carolan (2014).

The syntax for specifying a model in the ergm package follows R’s standard formula convention: my_network ~ ergm_term_1 + ergm_term_2 + ergm_term_3 and so forth.

This syntax is used for both the summary and ergm functions.

  • The summary function simply returns the numerical values of the network statistics in the model.

  • The ergm function estimates the model with those statistics.

As noted by the Statnet Development Team:

It is good practice to run a summmary command on any model before fitting it with ergm. This is the ERGM equivalent of performing some descriptive analysis on your covariates. This can help you make sure you understand what the term represents, and it can help to flag potential problems that will lead to poor modeling results.

Network Structure Parameters

Let’s start with with a simple model similar to the P1 model in Table 9.4. This model contains the required ergm-term edges that represents the total number of edges in the network, and the ergm-term mutual that examines the tendency for ties to be reciprocated, i.e. “mutuality”.

summary(leader_network ~ edges + mutual)
##  edges mutual 
##    143     36

We see from our summary that our leader_network consists of 143 edges and 36 reciprocated dyads.

Since the ergm() function uses a stochastic MCMC-based estimation algorithm, use the set.seed() function and set the value to 589 (note this could be any number besides our course number) so we produce the same results each time.

Now let’s estimate our model, save the results as ergm_mod_1, and use the summary() function again to take a look at our estimates. Also,

set.seed(589)

ergm_mod_1 <-ergm(leader_network ~ edges + mutual)

summary(ergm_mod_1)
## Call:
## ergm(formula = leader_network ~ edges + mutual)
## 
## Monte Carlo Maximum Likelihood Results:
## 
##        Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges   -3.1101     0.1283      0  -24.23   <1e-04 ***
## mutual   3.1246     0.3003      0   10.40   <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 2503.6  on 1806  degrees of freedom
##  Residual Deviance:  891.8  on 1804  degrees of freedom
##  
## AIC: 895.8  BIC: 906.8  (Smaller is better. MC Std. Err. = 0.7447)

Since ERGMs predict the presence of a network tie, with estimates indicating the importance of each to the presence of a tie, estimated coefficients can be explained in terms similar to logistic regression. That is, positive significant coefficients indicate that the corresponding parameters in the observed network (e.g. reciprocated ties between school leaders), controlling for all other parameters in the model, occur more than would be expected by chance, thus increasing the likelihood that a tie will occur, and vice-versa for negative coefficients.

After several iterations, we see that our model finally converges and suggest that the estimates for our edges and mutual terms are statistically significant. The negative estimate for the edge paramter, as noted by Carolan (2014), implies that the probability of a confidential exchange tie in year 3 is relatively low. The reciprocity parameter, on the other hand, is 3.11 in our model, and indicates a strong tendency for confidential-exchange ties to be reciprocated.

Now let’s add an ergm-term for transitivity. Modeling transitivity, or the “friend of a friend” phenomenon in social networks, is both computationally intensive because of the all the possible “triangles” in a directed network, and also very prone to model degeneration, i.e. when models that fail to converge by reaching a value expected by the parameters included in the model.

Ergm-terms for transitivity such as triangles and transitive in particular are prone to model degeneration. Fortunately, the {ergm} package includes a “more robust way of modeling triangles: the geometrically-weighed edgewise shared partner term (GWESP).”

Let’s add the gwesp term to our model with the suggested defaults from the Statnet Tutorial section on What it looks like when a model fails and take a look at the summary first:

summary(leader_network ~ edges + 
          mutual +
          transitive +
          gwesp(0.25, fixed=T))
##            edges           mutual       transitive gwesp.fixed.0.25 
##         143.0000          36.0000         202.0000         116.2971

As you can see, we have 202 transitive triad types.

Now lets run our model including transitivity and take a look at our estimates. Note that this make take a couple minutes to run. If you’d like to watch the number of iterations of the model set message = TRUE.

ergm_mod_2 <-ergm(leader_network ~ edges + 
                    mutual +
                    gwesp(0.25, fixed=T))

summary(ergm_mod_2)
## Call:
## ergm(formula = leader_network ~ edges + mutual + gwesp(0.25, 
##     fixed = T))
## 
## Monte Carlo Maximum Likelihood Results:
## 
##                  Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges             -3.9745     0.1522      0 -26.117   <1e-04 ***
## mutual             2.2989     0.3171      0   7.250   <1e-04 ***
## gwesp.fixed.0.25   1.0713     0.1305      0   8.207   <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 2503.6  on 1806  degrees of freedom
##  Residual Deviance:  805.7  on 1803  degrees of freedom
##  
## AIC: 811.7  BIC: 828.2  (Smaller is better. MC Std. Err. = 0.487)

Contrary to the analysis in Chapter 9, our model so far suggest that there is a tendency toward transitivity, that is a confidential exchange is likely to occur between two people who both have a confidential exchange with the same individual.

At this point, however, we have only focused on structural mechanisms inherent to networks themselves, or “global features” of the graph independent of actors, and have not looked at individual attributes among actors that may be shaping out network.

Actor Attribute Parameters

As noted by Carolan (2014), ERGMs have advanced to the point where actor attributes, or actor-level covariates, can now be incorporated into model estimation, such as individuals’ demographic (e.g., gender) or behavioral (e.g., efficacy) characteristics. These two attributes examined in Chapter 9 with leadership efficacy measured by survey items based on the Principal Efficacy Scale.

The attributes were included as both sender and receiver effects, which again can be very computationally intensive and thus timely to run, while sometimes prone to model degeneracy without some model fine tuning.

To simplify this approach, I’ve include these attributes using the nodefactor() and nocov() ergm-terms respectively which allow us to test whether those who are Male or have higher efficacy scores are more likely to either send or receive a confidential exchange.

ergm_3 <- ergm(leader_network ~ edges +
                 mutual +
                 gwesp(0.25, fixed=T) +
                 nodefactor('male') +
                 nodecov('efficacy')
               )
  
  
summary(ergm_3)
## Call:
## ergm(formula = leader_network ~ edges + mutual + gwesp(0.25, 
##     fixed = T) + nodefactor("male") + nodecov("efficacy"))
## 
## Monte Carlo Maximum Likelihood Results:
## 
##                   Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges             -4.35415    0.42693      0 -10.199   <1e-04 ***
## mutual             2.29853    0.31875      0   7.211   <1e-04 ***
## gwesp.fixed.0.25   1.05200    0.13281      0   7.921   <1e-04 ***
## nodefactor.male.1  0.10622    0.06547      0   1.622    0.105    
## nodecov.efficacy   0.02230    0.02967      0   0.752    0.452    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 2503.6  on 1806  degrees of freedom
##  Residual Deviance:  802.5  on 1801  degrees of freedom
##  
## AIC: 812.5  BIC: 840  (Smaller is better. MC Std. Err. = 0.4733)

Although there appears to be a very slight, though not statistically significant tendency towards male leaders being more likely to send or receive ties, the primary drivers of network formation appear to be structural features of the network including reciprocity and transitivity.

Sadly, we were not able to replicate the findings from Carolan (2014), but we also did not replicate their model exactly, nor did we fine tune our model parameters. Regardless, this provided a practical exercise for demonstrating ERGMs.

👉 Your Turn ⤵

Before moving on to checking the goodness of fit for our final model, try testing to see if shared characteristics among school leaders might help explain tie formation. To test for homophily among gender or location, the ergm() function includes a nodematch() argument similar to nodefactor() used above that could be used to answer the following question:

  1. Controlling for edges, reciprocity, and transitivity, are school leaders more likely to confide in colleagues of their own gender or who work at the same district-site level?
#YOUR CODE HERE
ergm_4 <- ergm(leader_network ~ edges +
                 mutual +
                 gwesp(0.25, fixed=T) +
                 nodematch('male') +
                 nodecov('efficacy')
)

  
  
summary(ergm_4)
## Call:
## ergm(formula = leader_network ~ edges + mutual + gwesp(0.25, 
##     fixed = T) + nodematch("male") + nodecov("efficacy"))
## 
## Monte Carlo Maximum Likelihood Results:
## 
##                  Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges            -4.21148    0.42044      0 -10.017   <1e-04 ***
## mutual            2.29888    0.32391      0   7.097   <1e-04 ***
## gwesp.fixed.0.25  1.06715    0.12779      0   8.351   <1e-04 ***
## nodematch.male   -0.14499    0.16827      0  -0.862    0.389    
## nodecov.efficacy  0.02326    0.02966      0   0.784    0.433    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 2504  on 1806  degrees of freedom
##  Residual Deviance:  805  on 1801  degrees of freedom
##  
## AIC: 815  BIC: 842.5  (Smaller is better. MC Std. Err. = 0.4946)

Once you’ve checked the summary of your model estimates, write your interpretation of the results below:

  • YOUR RESPONSE HERE

We’ve only scratched the surface for the number of network and node attributes properties that can be tested using ERGMs. Take a look at a table of common ergm-terms or type “ergm-terms” into the help menu of the files pane for some additional examples.

4c. Check Model Fit

One test of whether an ergm model is a “good fit” for the data is “how well it reproduces the observed global network properties that are not in the model.” This can be accomplished by choosing a network statistic that is not in the model, and comparing the value of this statistic observed in the original network to the distribution of values we get in simulated networks from our model, using the gof() function to test the goodness-of-fit.

The GOF Function

The gof() function is a bit different than the summary() and ergm() functions, in that it only takes 3 ergm-terms as arguments: degree, esp (edgwise share partners), and distance (geodesic distances). Each of these terms captures an aggregate network distribution, at either the node level (degree), the edge level (esp), or the dyad level (distance).

Let’s go ahead and run the gof() function on our ergm_3 model and plot the results:

ergm_3_gof <- gof(ergm_3)

plot(ergm_3_gof)

Overall, the model appears to fit reasonably well in that black line in our charts (the actual observed network measures) closely follows the aggregate measures generated from the simulations, or permutations of our network.

Check MCNC Diagnostics

One final check on our model is to examine the Monte-Carlo Markov Chain (MCMC) diagnostics to make sure our model is not heading of in the wrong direction and likely to never produce a network similar to what was observed, which is referred to as “model degeneracy.” As the statnet authors note:

When a model is not a good representation of the observed network, the simulated networks produced in the MCMC chains may be far enough away from the observed network that the estimation process is affected. In the worst case scenario, the simulated networks will be so different that the algorithm fails altogether.

Fortunately, our models did not fail to converge, the AIC and BIC indicators seemed to improve with each term added to the model, and the goodness-of-fit seemed to somewhat mirror the global features of our network.

Let run the mcmc.diagnostics() function on our final model anyways and check the results:

mcmc.diagnostics(ergm_3)
## Sample statistics summary:
## 
## Iterations = 7393280:31387648
## Thinning interval = 4096 
## Number of chains = 1 
## Sample size per chain = 5859 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                      Mean     SD Naive SE Time-series SE
## edges              -4.591  51.45   0.6721         2.6523
## mutual             -1.464  16.83   0.2199         0.8491
## gwesp.fixed.0.25   -4.969  63.58   0.8306         3.2645
## nodefactor.male.1  -4.895  59.53   0.7778         3.1039
## nodecov.efficacy  -62.654 696.91   9.1047        36.0250
## 
## 2. Quantiles for each variable:
## 
##                      2.5%     25%     50%    75%   97.5%
## edges               -93.0  -45.00  -4.000  34.00   93.00
## mutual              -30.0  -15.00  -2.000  11.00   31.55
## gwesp.fixed.0.25   -108.3  -56.99  -7.051  41.82  119.66
## nodefactor.male.1  -107.5  -52.00  -5.000  39.00  107.55
## nodecov.efficacy  -1259.6 -616.21 -55.990 454.95 1263.27
## 
## 
## Are sample statistics significantly different from observed?
##                  edges      mutual gwesp.fixed.0.25 nodefactor.male.1
## diff.      -4.59122717 -1.46373101       -4.9694724        -4.8953746
## test stat. -1.73105929 -1.72390925       -1.5222868        -1.5771922
## P-val.      0.08344119  0.08472423        0.1279372         0.1147513
##            nodecov.efficacy Overall (Chi^2)
## diff.          -62.65407919              NA
## test stat.      -1.73918395    4.130851e+01
## P-val.           0.08200241    1.006831e-07
## 
## Sample statistics cross-correlations:
##                       edges    mutual gwesp.fixed.0.25 nodefactor.male.1
## edges             1.0000000 0.9790145        0.9913763         0.9654745
## mutual            0.9790145 1.0000000        0.9813615         0.9488636
## gwesp.fixed.0.25  0.9913763 0.9813615        1.0000000         0.9623458
## nodefactor.male.1 0.9654745 0.9488636        0.9623458         1.0000000
## nodecov.efficacy  0.9988166 0.9781341        0.9906310         0.9634763
##                   nodecov.efficacy
## edges                    0.9988166
## mutual                   0.9781341
## gwesp.fixed.0.25         0.9906310
## nodefactor.male.1        0.9634763
## nodecov.efficacy         1.0000000
## 
## Sample statistics auto-correlation:
## Chain 1 
##               edges    mutual gwesp.fixed.0.25 nodefactor.male.1
## Lag 0     1.0000000 1.0000000        1.0000000         1.0000000
## Lag 4096  0.8563656 0.8611438        0.8607252         0.8547248
## Lag 8192  0.7507988 0.7552559        0.7531766         0.7492259
## Lag 12288 0.6641857 0.6662636        0.6647575         0.6661954
## Lag 16384 0.5869818 0.5872707        0.5853704         0.5909378
## Lag 20480 0.5156311 0.5189978        0.5140322         0.5225789
##           nodecov.efficacy
## Lag 0            1.0000000
## Lag 4096         0.8572204
## Lag 8192         0.7515372
## Lag 12288        0.6653919
## Lag 16384        0.5886474
## Lag 20480        0.5177033
## 
## Sample statistics burn-in diagnostic (Geweke):
## Chain 1 
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##             edges            mutual  gwesp.fixed.0.25 nodefactor.male.1 
##          0.020048         -0.045774         -0.026767          0.049916 
##  nodecov.efficacy 
##         -0.008712 
## 
## Individual P-values (lower = worse):
##             edges            mutual  gwesp.fixed.0.25 nodefactor.male.1 
##         0.9840052         0.9634904         0.9786455         0.9601895 
##  nodecov.efficacy 
##         0.9930490 
## Joint P-value (lower = worse):  0.4086241 .

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

One simple heuristic for interpreting mcmc.diagnostics() results is too look for “balance.” The charts with all the squiggles, to use a non-technical term, should have a scale that ideally is centered on 0 and have roughly straight line running through it, while the chart that looks like a simple line graph should also be ideally centered at 0 and roughly normally distributed on both ends. Overall, the diagnostics indicate there is room for improvement but nothing is wildly off and hence our model did not fail to converge.

👉 Your Turn ⤵

Use the code chunk below to check the model fit and diagnostics for the model the you created in 4b. Estimate the ERGM Model:

ergm_4_gof <- gof(ergm_4)

plot(ergm_4_gof)

mcmc.diagnostics(ergm_4)
## Sample statistics summary:
## 
## Iterations = 1753088:50331648
## Thinning interval = 16384 
## Number of chains = 1 
## Sample size per chain = 2966 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                     Mean     SD Naive SE Time-series SE
## edges             -5.047  55.74   1.0235         2.2935
## mutual            -1.744  18.15   0.3332         0.7609
## gwesp.fixed.0.25  -6.052  68.56   1.2588         2.8348
## nodematch.male    -1.883  26.22   0.4814         1.0466
## nodecov.efficacy -68.128 754.50  13.8540        31.0979
## 
## 2. Quantiles for each variable:
## 
##                   2.5%     25%     50%    75%  97.5%
## edges              -94  -52.00   -8.00  36.00  105.9
## mutual             -31  -17.00   -3.00  11.00   34.0
## gwesp.fixed.0.25  -109  -64.42  -11.24  43.19  133.7
## nodematch.male     -43  -24.00   -3.00  17.00   52.0
## nodecov.efficacy -1275 -704.80 -102.95 489.79 1427.3
## 
## 
## Are sample statistics significantly different from observed?
##                  edges      mutual gwesp.fixed.0.25 nodematch.male
## diff.      -5.04720162 -1.74376264       -6.0520442    -1.88267026
## test stat. -2.20065937 -2.29176493       -2.1349193    -1.79879464
## P-val.      0.02776015  0.02191921        0.0327676     0.07205117
##            nodecov.efficacy Overall (Chi^2)
## diff.           -68.1277680              NA
## test stat.       -2.1907498    2.462633e+01
## P-val.            0.0284699    1.781574e-04
## 
## Sample statistics cross-correlations:
##                      edges    mutual gwesp.fixed.0.25 nodematch.male
## edges            1.0000000 0.9821387        0.9919954      0.9738180
## mutual           0.9821387 1.0000000        0.9843173      0.9544524
## gwesp.fixed.0.25 0.9919954 0.9843173        1.0000000      0.9664271
## nodematch.male   0.9738180 0.9544524        0.9664271      1.0000000
## nodecov.efficacy 0.9989921 0.9816024        0.9914252      0.9732447
##                  nodecov.efficacy
## edges                   0.9989921
## mutual                  0.9816024
## gwesp.fixed.0.25        0.9914252
## nodematch.male          0.9732447
## nodecov.efficacy        1.0000000
## 
## Sample statistics auto-correlation:
## Chain 1 
##               edges    mutual gwesp.fixed.0.25 nodematch.male nodecov.efficacy
## Lag 0     1.0000000 1.0000000        1.0000000      1.0000000        1.0000000
## Lag 16384 0.6459246 0.6473194        0.6439058      0.6184790        0.6470462
## Lag 32768 0.4394995 0.4507851        0.4418863      0.4156574        0.4407662
## Lag 49152 0.2933342 0.3023061        0.2906277      0.2791341        0.2936363
## Lag 65536 0.1900030 0.1975078        0.1918624      0.1806380        0.1905881
## Lag 81920 0.1198984 0.1204661        0.1183082      0.1140627        0.1202699
## 
## Sample statistics burn-in diagnostic (Geweke):
## Chain 1 
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##            edges           mutual gwesp.fixed.0.25   nodematch.male 
##           0.9557           0.8491           0.9365           1.1736 
## nodecov.efficacy 
##           0.9592 
## 
## Individual P-values (lower = worse):
##            edges           mutual gwesp.fixed.0.25   nodematch.male 
##        0.3392118        0.3958112        0.3489923        0.2405570 
## nodecov.efficacy 
##        0.3374384 
## Joint P-value (lower = worse):  0.3277513 .

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

5. COMMUNICATE

For your final Your Turn, your goal is to distill our analysis from above into a simple “data product” designed to illustrate key findings about changes in the collaboration network over time. For the purposes of this task, imagine that your audience consists of the developers and facilitators of the DLT MOOC-Eds who have limited background in SNA and adapt the following steps accordingly:

  1. Select. Select your ERGM analysis from above, or a new analysis if so motivated, that you think would be interesting or relevant for the target audience and that helps answer our research questions.

  2. Polish. Create and polish a data visualization and/or data table to communicate your selected findings.

  3. Narrate. Write a brief narrative (2-3 paragraphs) to accompany your visualization and/or table that includes the following:

    • The question or questions guiding the analysis;

    • The conclusions you’ve reached based on our findings;

    • How your audience might use this information;

    • How you might revisit or improve upon this analysis in the future.

👉 Your Turn ⤵

Use the code chunk below create a polished table and/or visualization(s) and write a brief narrative in the space that follows.

Data Visualization or Table

m1b <- as.formula( leader_network ~ # The network
                     mutual('male', diff = TRUE) + # How many ties are reciprocated?
                     edges) # Edges term
fit1b <- ergm(m1b)

#summarize model
summary(fit1b)
## Call:
## ergm(formula = m1b)
## 
## Monte Carlo Maximum Likelihood Results:
## 
##                    Estimate Std. Error MCMC % z value Pr(>|z|)    
## mutual.same.male.0   1.5008     0.4711      0   3.186  0.00144 ** 
## mutual.same.male.1   2.9323     0.3695      0   7.935  < 1e-04 ***
## edges               -2.6898     0.1065      0 -25.257  < 1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 2503.6  on 1806  degrees of freedom
##  Residual Deviance:  952.1  on 1803  degrees of freedom
##  
## AIC: 958.1  BIC: 974.6  (Smaller is better. MC Std. Err. = 0.4556)
#print list and interpret P value
screenreg(list(fit1b))
## 
## ===============================
##                     Model 1    
## -------------------------------
## mutual.same.male.0     1.50 ** 
##                       (0.47)   
## mutual.same.male.1     2.93 ***
##                       (0.37)   
## edges                 -2.69 ***
##                       (0.11)   
## -------------------------------
## AIC                  958.08    
## BIC                  974.58    
## Log Likelihood      -476.04    
## ===============================
## *** p < 0.001; ** p < 0.01; * p < 0.05
#plot covariate of Male at 95% Confidence 
plotreg(list(fit1b),
        theme = theme_dark(),
        custom.model.names = "Covariates",
        custom.coef.names = c("Male", "Edge", "Female"),
        ggtitle = "Do males have more reciprocated ties than females?")

Narrative

I decided to explore the ERGM package and learn how to plot confidence intervals. I am still getting the hang of the interpretations and with what fuunctions to use. I will show more of that in the independent study.

My question wondered whether males have more reciprocated ties than females. This would mean our hypothesis is that gender has no effect on mutually reciprocating ties. In the summary output, we can see that males have a higher effect that is positive and significant. Moreover, males have a slightly higher propensity to engage in mutual behavior. Overall, the effects are quite similar for the two groups. In the list view P < .001 and we have strong evidence to reject the null hypothesis.

This information can be used to increase mutual reciprocation between all genders within a district. Being able to receive and send help when needed is important for a positive work environment. Additionally, understanding where the holes are central office can offer support and relationship building workshops. Reciprocity further allows important tasks to get done in a timely manner if you have help getting beyond your Zone of Proximinal development.

🧶 Knit & Check ✅

Congratulations - you’ve completed the Unit 3 case study! To share your work, click the drop down arrow next to the ball of yarn that says “Knit” at the top of this markdown file, then select “Knit top HTML”. Assuming your code contains no errors, this will create a web page in your Files pane that serves as a record of your work.

Once your file has been knitted, you can publish this file online using RPubs, or share the HTML file through another means.

References

Carolan, Brian. 2014. “Social Network Analysis and Education: Theory, Methods & Applications.” https://doi.org/10.4135/9781452270104.
McPherson, J Miller, and Lynn Smith-Lovin. 1987. “Homophily in Voluntary Organizations: Status Distance and the Composition of Face-to-Face Groups.” American Sociological Review, 370–79.
Pavel N. Krivitsky, Mark S. Handcock, David R. Hunter, Carter T. Butts, Chad Klumb, Steven M. Goodreau, and Martina Morris. 2003-2020. Statnet: Software Tools for the Statistical Modeling of Network Data. Statnet Development Team. http://statnet.org.
Thomas, Andrew C, and Joseph K Blitzstein. 2011. “Valued Ties Tell Fewer Lies: Why Not to Dichotomize Network Edges with Thresholds.” arXiv Preprint arXiv:1101.0788.