Statistical Learning: Modeling Social Networks
H. K. Tseng
5/20/2025
Rendering Social Networks with igraph
Visualizing graphs and analyzing networks using \(\textsf{igraph}\).
#install.packages(c("igraph", "igraphdata"))
# load packages
library(igraph)##
## Attaching package: 'igraph'## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## unionlibrary(igraphdata)
# Let's first create an undirected network consists of three nodes and three edges using igraph
g1 <- make_graph(c(1,2, 2,3, 3,1), n = 3, directed = FALSE) # Note how edge list c() are defined
class(g1)## [1] "igraph"plot(g1) 
# n (the number of nodes) can be greater than the number of edges in the edge list, meaning that many nodes are unconnected
g2 <- make_graph( edges=c(1,2, 2,3, 3,1), n = 10) # now with 10 vertices, and directed by default
plot(g2) 
# Name the nodes
g3 <- make_graph(c("A","B", "B","C", "C","A")) # `n` will be ignored for edge list with vertex names
# When the edge list has vertex names, the number of nodes is not needed
plot(g3)
# Named nodes without edges --> 'isolates'
g4 <- make_graph( c("A", "B", "B", "C", "B", "C", "A", "A"), isolates=c("D", "E", "F", "G"))
set.seed(123)
plot(g4, edge.arrow.size = .5,
vertex.color = "yellow",
vertex.frame.color = "gray",
vertex.label.color = "black",
vertex.size = 8,
vertex.label.cex = 1.5,
vertex.label.dist = 2, # The distance of the label from the center of the node
edge.curved = 0.2) 
# You can also use graph_from_literal() as a handy way to quickly create small graphs
# Undirected graph
plot(graph_from_literal(A-B, B-C)) # Two edges: node B to A, node B to C
# Directed graph
plot(graph_from_literal(A+-B, B-+C)) 
# Note: the number of dashes doesn't matter
# Use - for undirected edges, +- or -+ for directed edges, + means the direction pointing left or right. Use ++ for a symmetric relationships, and use : for sets of nodes, that is, set ABC to set DEF
plot(graph_from_literal(A-+B:C:D))
# Get edge, node and network attributes
# Plot a network
g5 <- make_graph( c("A", "B", "B", "C", "H", "C", "A", "E", "D", "C", "B", "E"),
isolates=c("F", "G"), directed = TRUE)
plot(g5)
# View the igraph object as an adjacency matrix
as_adjacency_matrix(g5)## 8 x 8 sparse Matrix of class "dgCMatrix"
## A B C H E D F G
## A . 1 . . 1 . . .
## B . . 1 . 1 . . .
## C . . . . . . . .
## H . . 1 . . . . .
## E . . . . . . . .
## D . . 1 . . . . .
## F . . . . . . . .
## G . . . . . . . .# View as a sparse matrix
as_adjacency_matrix(g5, sparse = FALSE) ## A B C H E D F G
## A 0 1 0 0 1 0 0 0
## B 0 0 1 0 1 0 0 0
## C 0 0 0 0 0 0 0 0
## H 0 0 1 0 0 0 0 0
## E 0 0 0 0 0 0 0 0
## D 0 0 1 0 0 0 0 0
## F 0 0 0 0 0 0 0 0
## G 0 0 0 0 0 0 0 0# g5[4, ] # Find Node H's ties with other nodes
# add edge g5['A','H'] = 1
# Get edge list
E(g5) ## + 6/6 edges from 4e6fa6a (vertex names):
## [1] A->B B->C H->C A->E D->C B->E# Count the number of edges
ecount(g5)## [1] 6# Get node (vertex) list
V(g5)## + 8/8 vertices, named, from 4e6fa6a:
## [1] A B C H E D F G# Count the number of nodes (vertices)
vcount(g5)## [1] 8# Find the neighbors of A
V(g5)[.nei("A")]## + 2/8 vertices, named, from 4e6fa6a:
## [1] B E# Find the neighbors of A and B
g5[[c("A","B")]]## $A
## + 2/8 vertices, named, from 4e6fa6a:
## [1] B E
##
## $B
## + 2/8 vertices, named, from 4e6fa6a:
## [1] C E# Alternatively, you can use
neighbors(g5,"A")## + 2/8 vertices, named, from 4e6fa6a:
## [1] B E# Find nodes that D is directed to
E(g5)[.from("D")]## + 1/6 edge from 4e6fa6a (vertex names):
## [1] D->C# Are F and H tied?
E(g5)["F" %--% "H"]## + 0/6 edges from 4e6fa6a (vertex names):# Find the tie that goes from F to H
E(g5)["F" %->% "H"]## + 0/6 edges from 4e6fa6a (vertex names):# Find nodes having ties directed to B
E(g5)[.to("B")]## + 1/6 edge from 4e6fa6a (vertex names):
## [1] A->B# Get node attributes
V(g5)$name## [1] "A" "B" "C" "H" "E" "D" "F" "G"# Add node attributes, for example. country name
V(g5)$country <- c("USA", "UKG", "AUS", "JPN", "FRN", "DEU", "NOR", "RUS")
# Add edge attributes, let's add numerical values to denote the importance of an edge, such as trade volume. We have a total of 6 edges: length(E(g5))
E(g5)$type <- "Million USD" # Applied to all edges
E(g5)$value <- c(25, 10, 7, 8, 4, 1)
# Set these new attributes
g5 <- set_graph_attr(g5, "type", "value")
# Assign value as edge weights
E(g5)$weight <- E(g5)$value
is_weighted(g5)## [1] TRUE# See node and edge attributes
vertex_attr(g5)## $name
## [1] "A" "B" "C" "H" "E" "D" "F" "G"
##
## $country
## [1] "USA" "UKG" "AUS" "JPN" "FRN" "DEU" "NOR" "RUS"edge_attr(g5)## $type
## [1] "Million USD" "Million USD" "Million USD" "Million USD" "Million USD"
## [6] "Million USD"
##
## $value
## [1] 25 10 7 8 4 1
##
## $weight
## [1] 25 10 7 8 4 1# You can also remove (delete) certain attributes: g5 <- delete_graph_attr(g5, "value")
# Quickly check basic properties
g5## IGRAPH 4e6fa6a DNW- 8 6 --
## + attr: type (g/c), name (v/c), country (v/c), type (e/c), value (e/n),
## | weight (e/n)
## + edges from 4e6fa6a (vertex names):
## [1] A->B B->C H->C A->E D->C B->E# Serial number and letters refer to descriptive details of this graph object g5: U means this is an undirected graph; D means it is directed graph; N indicates it is a named graph.
# '--' refers to attributes not applicable to g:
# W would refer to a weighted graph --> edges have a weight attribute
# B would refer to a bipartite graph --> vertices have a type attribute
# 8 refers to the number of nodes (vertices) in g, 6 refers to the number of edges in g
# (g/c) means graph/character, (v/c) vertex/character, (e/n) edge/numeric, (e/c) edge/character
## Network Measures
# list degree centrality
degree(g5, v = V(g5), # The degree of a node records the number of its adjacent edges.
mode = c("all", "out", "in", "total"),
loops = TRUE, normalized = FALSE)## A B C H E D F G
## 2 3 3 1 2 1 0 0# in-degree
degree(g5, mode="in") # How many edges are pointing to a particular node## A B C H E D F G
## 0 1 3 0 2 0 0 0# out-degree
degree(g5, mode="out") # How many edges are from a particular node## A B C H E D F G
## 2 2 0 1 0 1 0 0# Eigenvector Centrality
# Eigenvector centrality is a measure of the influence of a node in a network. Scores are assigned to other nodes in the network based on the idea that connections to high-scoring nodes contribute more to the score of the node than to low-scoring nodes. A high eigenvector score means that a node is connected to many nodes who themselves have high scores.
evcent(g5)$vector## Warning: `evcent()` was deprecated in igraph 2.0.0.
## ℹ Please use `eigen_centrality()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## A B C H E D
## 9.685752e-01 1.000000e+00 3.836254e-01 9.469184e-02 3.084933e-01 5.410962e-02
## F G
## 8.026701e-20 8.026701e-20# Hubs and Authorities: similar measures to Eigenvector centrality
hub.score(g5)$vector## Warning: `hub.score()` was deprecated in igraph 2.0.0.
## ℹ Please use `hits_scores()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## A B C H E D
## 1.0000000000 0.0138481147 0.0000000000 0.0015531986 0.0000000000 0.0008875421
## F G
## 0.0000000000 0.0000000000authority.score(g5)$vector## Warning: `authority.score()` was deprecated in igraph 2.0.0.
## ℹ Please use `hits_scores()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## A B C H E D
## 0.000000000 1.000000000 0.006116148 0.000000000 0.320553925 0.000000000
## F G
## 0.000000000 0.000000000# Closeness
# Calculated as the reciprocal of the sum of the length of the shortest paths between the node and all other nodes in the graph. The more central a node is, the closer it is to all other nodes.
closeness(g5)## A B C H E D F
## 0.01470588 0.09090909 NaN 0.14285714 NaN 0.25000000 NaN
## G
## NaN# Betweeness
# A measure of how often a node is on the shortest path between other nodes in the network. The higher the measure, the more important a node is.
betweenness(g5)## A B C H E D F G
## 0 1 0 0 0 0 0 0# Generate empty graph
empty_g <- make_empty_graph(20)
plot(empty_g, vertex.size = 5, vertex.label=NA)
# Generate full graph
full_g <- make_full_graph(20)
plot(full_g, vertex.size = 5, vertex.label = NA)
# Star graph
star_g <- make_star(20)
plot(star_g, vertex.size = 5, vertex.label = NA) 
# Generate scale-free Barabasi networks
bara_g <- sample_pa(200) # You can also use barabasi.game(200)
plot(degree_distribution(bara_g), xlab = "Node", ylab = "Degree") # Does degree distribution follow power law?
# Erdos-Renyi random graph
erdo_g <- sample_gnm(n = 20, m = 20) # You can also use erdos.renyi.game. n is the number of nodes, m the number of edges
## Options include directed = and loops =
plot(erdo_g, vertex.size = 5, vertex.label = NA) 
# Generate a graph using graph_from_adjacency_matrix()
set.seed(123)
# Generate a sparse matrix by sampling from Bernoulli distribution of size 64.
adj_m <- matrix(sample(0:1, 64, replace = TRUE, prob = c(0.9, 0.1)),
nc = 8) # Number of columns
adj_m## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0
## [3,] 0 1 0 0 0 0 0 0
## [4,] 0 0 1 0 0 0 0 0
## [5,] 1 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0
## [7,] 0 0 0 1 0 0 0 0
## [8,] 0 0 1 1 0 0 0 0g_adj <- graph_from_adjacency_matrix(adj_m,
diag = FALSE, # If this is FALSE then the diagonal is zeroed out first
mode = "max") # Default is directed, the latest version requires one to set mode = "max" in lieu of "undirected"
plot(g_adj)
## Generate graph from uploaded data (usually stored in edgelist format) using graph_from_edgelist()
# install.packages("remotes")
# remotes::install_github("RandiLGarcia/dyadr")
library(dyadr) # We'll be using the dyadic_trade data from Randi Garcia's dyadr package https://rdrr.io/github/RandiLGarcia/dyadr/man/dyadic_trade.html
data(dyadic_trade)
# Take a look at the data
head(dyadic_trade)## ccode1 ccode2 year importer1 importer2 flow1 flow2 source1
## 1 2 20 1920 United States of America Canada 611.86 735.48 1
## 2 2 20 1921 United States of America Canada 335.44 442.99 1
## 3 2 20 1922 United States of America Canada 364.02 502.84 1
## 4 2 20 1923 United States of America Canada 416.00 598.14 1
## 5 2 20 1924 United States of America Canada 399.14 496.32 1
## 6 2 20 1925 United States of America Canada 473.72 608.35 1
## source2
## 1 1
## 2 1
## 3 1
## 4 1
## 5 1
## 6 1# Subset (retrieve data from year 2009)
trade_df <- dyadic_trade[dyadic_trade$year == 2009, ]
# Data wrangling
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:igraph':
##
## as_data_frame, groups, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)## Warning: package 'tidyr' was built under R version 4.4.3##
## Attaching package: 'tidyr'## The following object is masked from 'package:igraph':
##
## crossingd_trade = trade_df %>% dplyr::select(importer1, importer2, flow1) %>%
group_by(importer1, importer2) %>% # Organize the data into dyad format with 3 columns
summarise(USD_m = sum(flow1))## `summarise()` has grouped output by 'importer1'. You can override using the
## `.groups` argument.# Convert to dataframe and remove NAs
d_trade <- as.data.frame(d_trade)
d_trade <- d_trade[d_trade$importer1 != "", ]
d_trade <- d_trade[d_trade$importer2 != "", ]
d_trade <- d_trade[d_trade$USD_m != "", ]
head(d_trade)## importer1 importer2 USD_m
## 2 Afghanistan Australia 25.27
## 3 Afghanistan Bangladesh 4.56
## 4 Afghanistan Brunei 0.03
## 5 Afghanistan Cambodia 1.98
## 6 Afghanistan China 239.26
## 7 Afghanistan Fiji 0.00# Randomly sample 100 obs
set.seed(123)
d_trade <- sample_n(d_trade, 100)
# Convert to matrix
trade_m <- as.matrix(d_trade[, 1:2])
# Rescale trade volume (million USD), it will be used to represent edge width
d_trade$USD_m <- d_trade$USD_m / 80
trade_g <- graph_from_edgelist(trade_m, directed = FALSE)
plot(trade_g,
vertex.size = 6,
edge.width = d_trade[,3], # Set edge.width proportional to bilateral trade volume
vertex.color = "coral")
## Alternatively, you can also generate a network via graph_from_data_frame() that loads nodes (vertices) and edges separately. Here's my example: International Armament Collaboration data
# Load required packages
library(foreign)
# Load the data
vertices <- read.csv(url("https://www.dropbox.com/s/j1o1wh5w2ucvkje/vertices.csv?dl=1"), header = TRUE)
edges <- read.csv(url("https://www.dropbox.com/s/p78ndqv4olktifu/IAC_edges.csv?dl=1"), header = TRUE)
# Recode and generate color labels
is.na(vertices$alliance) <- vertices$alliance == ""
vertices$c_color <- c("white","orange","red","olivedrab1")[vertices$alliance] # Partners' alliance status
vertices$m_color <- ifelse(vertices$mutual_defense > 0,'black','grey80') # If a partner holds mutual defense treaty with the U.S.
color_alliance <- vertices %>%
group_by(alliance, c_color) %>%
summarise(n = n()) %>%
filter(!is.na(c_color))## `summarise()` has grouped output by 'alliance'. You can override using the
## `.groups` argument.color_defense <- vertices %>%
group_by(mutual_defense, m_color) %>%
summarise(n = n()) %>%
filter(!is.na(m_color))## `summarise()` has grouped output by 'mutual_defense'. You can override using
## the `.groups` argument.# Convert into graph format
g <- graph_from_data_frame(d = edges, vertices = vertices, directed = FALSE)
# Mapping
V(g)$color[V(g)$alliance == "NATO"] = "red"
V(g)$color[V(g)$alliance == "MNNA"] = "orange"
V(g)$color[V(g)$alliance == "USA"] = "olivedrab1"
w1 <- E(g)$IAC # Number of IAC cases with the U.S.
w1 <- log(w1) # Log transform
col = ifelse(V(g)$mutual_defense > 0,'black','grey80')
line_type <- ifelse(!is.na(E(g)$mutual_defense), 1, 2)
m1 <- layout_nicely(g)
V(g)$vertex_degree <- igraph::degree(g)
# Plotting
par(mar=c(1,1,1,1))
plot(g,
vertex.size = 7,
vertex.label.color = "black",
vertex.frame.color = col,
vertex.label.cex = 0.4,
edge.color = 'grey30',
edge.width = w1+1,
edge.lty = line_type,
edge.arrow.size = 0.6,
layout = m1)
legend("topleft", legend = c("Alliance type", as.character(color_alliance$alliance), "Mutual defense", as.character(color_defense$mutual_defense)), pch = c(1, 19, 19, 19, 19, 1), col = c(NA, color_alliance$c_color, NA, color_defense$m_color), pt.cex = 1, cex = 0.7, bty = "n", ncol = 1, title = "")
legend("topleft", legend = "", cex = 0.9, bty = "n", ncol = 1,
title = "Nature of Security Ties")
Re-analyze FTA formation using static and temporal Exponential Random Graphs Models (ERGM)
# Load required packages
library(ergm)
library(btergm) # For binary temporal ERGM
library(speedglm) # For the use of glm inside btergm()
library(statnet)
library(network)
# Load FTA data again
FTA <- readRDS(url("https://www.dropbox.com/scl/fi/6luyfnrur7ictisun8zxh/FTA19982019.rds?rlkey=h6af5x79e318vm6q7rj8b8loc&st=bpi38ulh&dl=1"))
# Load cid again for country indexing
cid <- read.csv(url("https://www.dropbox.com/scl/fi/7isucibr6woxr5gnlthu8/countryid.csv?rlkey=cezf0zg1vbkcvxrrvazzu4pjf&st=uxw2yxdw&dl=1"))
# Pretty much the same, but this time, we need to tweak the data a bit by converting the data to network class via network()
for(i in 1:length(FTA)){
colnames(FTA[[i]]) <- cid$comtrade_id
rownames(FTA[[i]]) <- cid$comtrade_id
FTA[[i]] <- network(as.matrix(FTA[[i]]),
matrix.type = 'adjacency', # Convert to adjacency matrix
ignore.eval = TRUE,
directed = TRUE) # Directed network
}
FTA <- FTA[18:22] # Subset the last five years (2015-2019)
# Load dyadic trade data again
trade2015 <- readRDS(url("https://www.dropbox.com/scl/fi/17i8e42xzhgo9dqdmzohn/trade2015.rds?rlkey=5ldvbtn4550vwrbobsu5lfuhx&st=ipp306c7&dl=1"))
trade2016 <- readRDS(url("https://www.dropbox.com/scl/fi/bevq047cg7q6p96qk8rk3/trade2016.rds?rlkey=x19pybshmqf45aqclc8audnqj&st=qn68qrog&dl=1"))
trade2017 <- readRDS(url("https://www.dropbox.com/scl/fi/1a2wvu3o28ur5ce462blr/trade2017.rds?rlkey=izxqr82e8iirr9h73b48i3yfx&st=gavsra6g&dl=1"))
trade2018 <- readRDS(url("https://www.dropbox.com/scl/fi/tvqyzcgqvkfo7kxpkn9j8/trade2018.rds?rlkey=sibcdornmgwzznb5myjfaeqju&st=bnssd0n7&dl=1"))
trade2019 <- readRDS(url("https://www.dropbox.com/scl/fi/gxustrrcx6kvw9137uttd/trade2019.rds?rlkey=jjvgc0ju17465bv6e9jjnebu3&st=cpgy2cry&dl=1"))
## Combine annual dyadic trade into a list object
dyadic_trade <- list(trade2015, trade2016, trade2017, trade2018, trade2019)
## Load country attributes data (again)
country_attr <- readRDS(url("https://www.dropbox.com/scl/fi/ddhf2z5tl9bo8l9jo29vb/country_attr.rds?rlkey=5upfqxcljtuphtuvingnptptm&st=9bc9dddy&dl=1"))
# Setting nodal attributes within network objects using set.vertex.attribute() function, so that these variables will be attached to the array of network objects, FTA, as nodal attributes
for(i in unique(country_attr$year)){
population <- country_attr[country_attr$year == i, c("population")]
trade_open <- country_attr[country_attr$year == i, c("trade_open")]
network::set.vertex.attribute(FTA[[(i-2014)]], value = population, attrname = "population") # Use i - 2014 in order to sort the index in FTA list
network::set.vertex.attribute(FTA[[(i-2014)]], value = trade_open, attrname = "trade_open")
}
set.seed(123)
# Baseline model. select a single year (2015) and estimate tie probability of that year
FTA.ergm <- ergm(FTA[[1]] ~ edgecov(dyadic_trade[[1]]) + nodecov("population") + nodecov("trade_open")) # Here we use edgecov() and nodecov() to create edge and node covariates inside ergm() function## Starting maximum pseudolikelihood estimation (MPLE):## Obtaining the responsible dyads.## Evaluating the predictor and response matrix.## Maximizing the pseudolikelihood.## Finished MPLE.## Evaluating log-likelihood at the estimate.summary(FTA.ergm)## Call:
## ergm(formula = FTA[[1]] ~ edgecov(dyadic_trade[[1]]) + nodecov("population") +
## nodecov("trade_open"))
##
## Maximum Likelihood Results:
##
## Estimate Std. Error MCMC % z value Pr(>|z|)
## edgecov.dyadic_trade[[1]] 0.023837 0.002330 0 10.230 <1e-04 ***
## nodecov.population -0.067472 0.001540 0 -43.819 <1e-04 ***
## nodecov.trade_open 0.001270 0.000184 0 6.898 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 35711 on 25760 degrees of freedom
## Residual Deviance: 23158 on 25757 degrees of freedom
##
## AIC: 23164 BIC: 23189 (Smaller is better. MC Std. Err. = 0)# Wow, previously insignificant variables now become significant! ERGM does appear better at predicting ties.
# Adding network structural characteristics
# density: equivalent to "edges" or "istar(1)" or "ostar(1)" divided by n*(n-1).
# mutual: reciprocity term
# twopath: this term adds one statistic to the model, equal to the number of 2-paths in the network i->j->k
# Let's add the "twopath" term as an exploratory probe. Please do take a look at the verbose log of the Monte Carlo maximum likelihood (MCMLE) process
FTA_structural.ergm <- ergm(FTA[[1]] ~ edgecov(dyadic_trade[[1]]) + nodecov("population") + nodecov("trade_open") + twopath,
control = control.ergm(parallel = 4, # Use 4 CPU cores
parallel.type = "PSOCK",
MCMC.samplesize = 5000)) # Number of network statistics, randomly drawn from a given distribution on the set of all networks, returned by the Metropolis-Hastings algorithm.## Starting maximum pseudolikelihood estimation (MPLE):
## Obtaining the responsible dyads.
## Evaluating the predictor and response matrix.
## Maximizing the pseudolikelihood.
## Finished MPLE.
## Starting Monte Carlo maximum likelihood estimation (MCMLE):
## Iteration 1 of at most 60:## Warning: 'glpk' selected as the solver, but package 'Rglpk' is not available;
## falling back to 'lpSolveAPI'. This should be fine unless the sample size and/or
## the number of parameters is very big.## 1 Optimizing with step length 1.0000.
## The log-likelihood improved by 3.3274.
## Estimating equations are not within tolerance region.
## Iteration 2 of at most 60:
## 1 Optimizing with step length 1.0000.
## The log-likelihood improved by 0.4189.
## Estimating equations are not within tolerance region.
## Iteration 3 of at most 60:
## 1 Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0406.
## Convergence test p-value: 0.0721. Not converged with 99% confidence; increasing sample size.
## Iteration 4 of at most 60:
## 1 Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0133.
## Convergence test p-value: 0.5063. Not converged with 99% confidence; increasing sample size.
## Iteration 5 of at most 60:
## 1 Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0228.
## Convergence test p-value: 0.1542. Not converged with 99% confidence; increasing sample size.
## Iteration 6 of at most 60:
## 1 Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0219.
## Convergence test p-value: 0.0091. Converged with 99% confidence.
## Finished MCMLE.
## Evaluating log-likelihood at the estimate. Fitting the dyad-independent submodel...
## Bridging between the dyad-independent submodel and the full model...
## Setting up bridge sampling...
## Using 16 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 .
## Bridging finished.
##
## This model was fit using MCMC. To examine model diagnostics and check
## for degeneracy, use the mcmc.diagnostics() function.summary(FTA_structural.ergm)## Call:
## ergm(formula = FTA[[1]] ~ edgecov(dyadic_trade[[1]]) + nodecov("population") +
## nodecov("trade_open") + twopath, control = control.ergm(parallel = 4,
## parallel.type = "PSOCK", MCMC.samplesize = 5000))
##
## Monte Carlo Maximum Likelihood Results:
##
## Estimate Std. Error MCMC % z value Pr(>|z|)
## edgecov.dyadic_trade[[1]] 1.813e-02 2.066e-03 0 8.777 <1e-04 ***
## nodecov.population -8.873e-02 1.760e-03 0 -50.427 <1e-04 ***
## nodecov.trade_open 1.896e-05 1.564e-04 0 0.121 0.903
## twopath 1.743e-02 9.250e-04 0 18.837 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 35711 on 25760 degrees of freedom
## Residual Deviance: 22951 on 25756 degrees of freedom
##
## AIC: 22959 BIC: 22991 (Smaller is better. MC Std. Err. = 3.146)# After adding twopaths, trade_open lost its significance.Longitudinal Analysis of FTA formation: BTERGM
# Fit a plain vanilla base btergm
FTA.btergm <- btergm(FTA ~ edgecov(dyadic_trade) + nodecov("population") + nodecov("trade_open"))
# Adding memory term, note that multiple memory term options cannot be used together (you can only choose one)
# NOT RUN
# Dyadic stability (type = "stability"): checks whether both edges and non-edges are stable between previous and current network
FTA_stability.btergm <- btergm(FTA ~ edgecov(dyadic_trade) + nodecov("population") + nodecov("trade_open") + memory(type = "stability"))
# You may also consider positive autoregression (type = "autoregression"): checks whether previous ties are carried over to the current network
# Lagged by 1 period
FTA_auto.btergm <- btergm(FTA ~ edgecov(dyadic_trade) + nodecov("population") + nodecov("trade_open") + memory(type = "autoregression", lag = 1))Results.
# btergm FTA estimation results
FTA.btergm <- readRDS(url("https://www.dropbox.com/scl/fi/x4c2c56o21ryxarlhfjbw/FTA.btergm.rds?rlkey=z7hu0v1z7u2qf21vne1t2r1ix&st=475rvyzi&dl=1"))
FTA_stability.btergm <- readRDS(url("https://www.dropbox.com/scl/fi/0yf039ydu9ohuxzq5ba7r/FTA_stability.btergm.rds?rlkey=gvp889g0r4v4igdcxsxkoe3sw&st=r7estp39&dl=1"))
# Results at a glance
summary(FTA.btergm)## ==========================## Summary of model fit## ==========================##
## Formula: FTA ~ edgecov(dyadic_trade) + nodecov("population") + nodecov("trade_open")## Time steps: 5## Bootstrapping sample size: 500## Estimates and 95% confidence intervals:## Estimate Boot mean 2.5% 97.5%
## edgecov.dyadic_trade[[i]] 5.6619e-02 0.05712394 0.0518 0.0635
## nodecov.population -6.9177e-02 -0.06946778 -0.0735 -0.0654
## nodecov.trade_open -6.2946e-05 -0.00007736 -0.0002 0.0000summary(FTA_stability.btergm)## ==========================## Summary of model fit## ==========================##
## Formula: FTA ~ edgecov(dyadic_trade) + nodecov("population") + nodecov("trade_open") + memory(type = "stability")## Time steps: 4## Bootstrapping sample size: 500## Estimates and 95% confidence intervals:## Estimate Boot mean 2.5% 97.5%
## edgecov.dyadic_trade[[i]] 0.0309516 0.0243176 -0.0014 0.0372
## nodecov.population 0.0144705 0.0183102 0.0117 0.0291
## nodecov.trade_open 0.0018269 0.0015195 -0.0007 0.0029
## edgecov.memory[[i]] 5.8766994 6.0160961 5.3119 7.3824# Only population and structural variables remain significant.Assessing prediction accuracy of BTERGM
# NOT RUN
# Use simulated networks (via gof() function) from the estimation of [t, t+1] data and compare them to the actual network at t+2
# Model 1 (baseline model)
pred1_auc <- list()
for (i in 1:3){ # A total of 3 predictions will be made
train_t <- c(i, i+1) # Set time frame (use two years' data as training data)
dy_trade <- dyadic_trade[train_t] # Subset training data time frame
model1_fit <- btergm(FTA[train_t] ~
edgecov(dy_trade) + nodecov("population") + nodecov("trade_open"), # edge and nodal covariates
R = 1000, # 1000 bootstrap replications
parallel = "multicore", # Parallel and multi-core processing
ncpus = 4) # Number of CPU cores used
model1_test <- i+2 # Number of predictions: 3 predictions
g <- gof(model1_fit, statistics = c(rocpr), nsim = 100, target = FTA[model1_test]) # Retrieve ROC-Precision metric
pred1_auc[[i]] <- g$`Tie prediction` # Store tie prediction result in pred1_auc list
}
# Model 2 (baseline model w/ stability term)
pred2_auc <- list()
for (i in 1:3){
train_t <- c(i+1, i+2)
dy_trade <- dyadic_trade[train_t]
model2_fit <- btergm(FTA[train_t] ~
edgecov(dy_trade) + nodecov("population") + nodecov("trade_open") + + memory(type = "stability"),
R = 1000, # 1000 bootstrap replications
parallel = "multicore", # Parallel and multi-core processing
ncpus = 4) # Number of CPU cores used
model2_test <- i+2
g <- gof(model2_fit, statistics = c(rocpr), nsim = 100, target = FTA[model2_test])
pred2_auc[[i]] <- g$`Tie prediction`
}Create dataframe for plotting
# Load previously generated prediction result list
pred1_auc <- readRDS(url("https://www.dropbox.com/scl/fi/hzws9swlsytyybextjn9r/pred1_auc.rds?rlkey=c2hsmntk3w18o9n4ci3lcvcpk&st=1gmvilhu&dl=1"))
pred2_auc <- readRDS(url("https://www.dropbox.com/scl/fi/ru05hs3xymsqzlfdq8byv/pred2_auc.rds?rlkey=fbr3pmjgzuek7h4x4e9tl46yd&st=l1da7non&dl=1"))
# Take a look at the estimated model AUC
pred1_auc[[1]]## ## ROC model ROC random PR model PR random
## 1 0.5444093 0.5103156 0.2170958 0.1958238# As expected, PRs are generally lower
pred2_auc[[1]]## ## ROC model ROC random PR model PR random
## 1 0.9610709 0.5091163 0.9487324 0.1925754# After accounting for the ties of previous networks, both ROC AUC and PR AUC rose to near .99
# Create dataframe and fill two new data columns with NA
mod1_auc <- data.frame(Model = "Baseline", year = 2017:2019)
mod1_auc$auc.roc <- NA
mod1_auc$auc.pr <- NA
for (i in 1:3){
mod1_auc$auc.roc[i] <- pred1_auc[[i]]$auc.roc # Subset the AUC of ROC
mod1_auc$auc.pr[i] <- pred1_auc[[i]]$auc.pr # Subset the AUC of PR
}
mod2_auc <- data.frame(Model = "Dyadic Stability", year = 2017:2019)
mod2_auc$auc.roc <- NA
mod2_auc$auc.pr <- NA
for (i in 1:3){
mod2_auc$auc.roc[i] <- pred2_auc[[i]]$auc.roc
mod2_auc$auc.pr[i] <- pred2_auc[[i]]$auc.pr
}
# Combine data
df <- rbind(mod1_auc, mod2_auc)
df$year <- as.factor(df$year)
# ROC & PR AUC
library(ggplot2)
roc_auc <- ggplot(df, aes(x = year, y = auc.roc, group = Model)) +
geom_line(aes(linetype = Model, color = Model), lwd = 0.75, alpha = 0.75) +
scale_color_manual(values = c("#F8766D", "#00BFC4")) +
scale_linetype_manual(values = c("dashed", "solid")) +
theme_bw() +
labs(title = "ROC", x ="Year", y = "AUC") +
theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
# PR AUC
pr_auc <- ggplot(df, aes(x = year, y = auc.pr, group = Model)) +
geom_line(aes(linetype = Model, color = Model), lwd = 0.75, alpha = 0.75) +
scale_color_manual(values = c("#F8766D", "#00BFC4")) +
scale_linetype_manual(values = c("dashed", "solid")) +
theme_bw() +
labs(title="Precision-Recall", x = "Year", y = "AUC") +
theme(plot.title = element_text(hjust = 0.5), legend.box.background = element_rect(colour = "black"))
# Load required packages
library(gridExtra)
library(grid)
library(cowplot)
# Render two subplots side-by-side
plot_grid(roc_auc, pr_auc, nrow = 1, rel_widths = c(0.4, 0.6))
Social Network Analysis using Stochastic Actor-Oriented Models (SAOM)