editor_path <- if (requireNamespace("rstudioapi", quietly = TRUE) &&
rstudioapi::isAvailable()) {
tryCatch(rstudioapi::getSourceEditorContext()$path, error = function(e) "")
} else ""
day8_candidates <- unique(c(
if (nzchar(editor_path)) dirname(editor_path) else character(),
"day8_netify", "1_netify", "."
))
day8_candidates <- day8_candidates[
vapply(file.path(day8_candidates, "_setup.R"), file.exists, logical(1))
]
if (!length(day8_candidates)) {
stop("Cannot find the Day 8 folder. Open the Day 8 QMD from its downloaded folder, ",
"or set the working directory to the folder containing the Day 8 materials.", call. = FALSE)
}
day8_dir <- normalizePath(day8_candidates[[1]], mustWork = TRUE)
source(file.path(day8_dir, "_setup.R"))
require_pkgs(c("netify", "igraph", "ggplot2", "network", "sandwich",
"lmtest", "dplyr", "generics", "digest"))
suppressPackageStartupMessages({
library(netify)
library(igraph)
library(ggplot2)
})
theme_set(theme_bw())
set.seed(6886)Day 8: Netify, Building the Object and Mapping the Whole Course
Advanced Network Analysis · ICPSR
ImportantRun This Setup First
Run the setup chunk below before beginning the walkthrough. It loads the packages, finds the Day 8 folder whether you are using the downloaded day8_netify folder or the internal course source folder, and defines the local data loader used throughout the document.
NoteHow to use this document
Everything in here runs. There is no separate slide deck hiding the code. This is the script, with the reasoning left between the chunks. You can read it from top to bottom afterwards without losing the thread of the analysis.
A few conventions you will see all week:
- Inline numbers are live. When a sentence says a network has some density, that number was computed from the object on the line above it, not typed in by me. So if you re-run with different data, the prose updates too.
- Expensive fits are cached. Later days fit models that take minutes. Those are wrapped in
cache_fit(), which stores the result incache/so the first person pays for it once and everyone after loads it in milliseconds. Today has almost nothing expensive in it, so almost everything just runs. - No absolute paths. Every file this document touches lives in its own folder. Unzip
1_netify/, open this, and it runs.
Today you need three packages, all on CRAN:
install.packages(c("netify", "igraph", "ggplot2", "statnet", "sandwich",
"lmtest", "dplyr", "digest"))1 Where we are: Day 8
Welcome back. The first seven days gave you the foundations of network analysis. Today begins the advanced modeling sequence: getting a coefficient, an interval, a prediction, and a test while taking network dependence seriously. We will get to that problem in the second half of today, and it is the thread that ties Days 8 through 14 together.
Here is where we are going, so the plumbing has a point. Before the first break I will show you a perfectly ordinary-looking regression, the kind in every methods textbook, that screams “significant” about 50% of the time when the true effect is exactly zero. Everything that follows is about how you build the object so that you are not the person who published it.
But before you can model a network you have to have one, as an object, in the right shape, with the right things attached. That sounds like plumbing, and half of it is. The other half is a series of quiet decisions: Is a tie directed? Is it weighted? Who counts as an actor? What does a missing dyad mean? Every one of those decisions changes the model you get at the end. Get the object wrong and no amount of clever modeling later will save you.
So today has two parts:
- Part one: the object. We build networks from raw data with
netify, justify every argument as we go, attach node- and dyad-level covariates, read what the object is telling us, plot it, and hand it off in the exact shape each modeling package wants. This is the substrate the rest of the course stands on. - Part two: the map. Once you have the object, what do you actually do with it? There is a small map of methods, one cell per day for the rest of the week, and it all comes from a single problem. We will draw that map and point at Day 9.
1.1 Let’s Start With Your Research Question
Before netify() gets an argument, the research question has already made most of the important decisions. Consider three familiar political science projects:
- Alliance formation: states are the actors, a defense commitment is the tie, direction depends on how the treaty is coded, and the risk set excludes state-pairs that did not coexist.
- Legislative collaboration: legislators are the actors, cosponsorship or joint voting is the tie, and a Congress or legislative session defines the time window.
- Conflict and cooperation: states or armed groups are the actors, events are the raw records, and aggregation determines whether the outcome is any interaction, the number of interactions, or their intensity.
Those are not cosmetic differences. They produce different outcomes and different populations. A claim about why states initiate conflict needs a directed network. A claim about whether alliances cluster needs a binary or weighted alliance definition. A claim about change needs comparable actor rosters and time periods.
The construction sequence for every example today is:
- Name the substantive event. What happened between which political actors?
- Define the tie. Decide its direction, scale, and time window.
- Define eligibility. Determine which actor-pairs could have produced an observed tie.
- Build the network outcome. Convert the raw records into an adjacency matrix or longitudinal array.
- Attach covariates by actor identity. Do not rely on row order.
- Audit the result. Check actors, ties, missingness, direction, and dimensions before fitting anything.
The code in this session follows those six steps. The point is not simply to create an R object. It is to make the estimand visible before a model quietly inherits it.
1.2 The Applied Compass for Today
Even today, before the main modeling sequence begins, every tool has to earn its place in an applied argument.
| Tool or approach | Question we want to understand | What it contributes | What makes it distinct | What the result sounds like in actual units |
|---|---|---|---|---|
| Network construction | Who could relate to whom, what counts as a tie, and over what period? | Turns political records into a defensible outcome and risk set | It defines the estimand before any estimator sees the data | “One cell records how many cooperative statements country A directed toward country B during one calendar year. A value of 12 means twelve coded statements. A zero means both countries existed and no statement was recorded. A missing value means we cannot treat that country-pair-year as an observed opportunity with no cooperation.” |
| Descriptive network summaries | Which countries initiated the most actions, which countries were the targets of the most actions, and which pairs acted toward one another? | Makes the relational pattern visible and suggests puzzles worth modeling | It describes the observed graph without claiming a data-generating process | “An out-strength of 240 means a country initiated 240 coded cooperative statements across all partner countries that year. An in-strength of 240 means other countries directed 240 such statements toward it. Those are different political actions, not two labels for the same kind of activity.” |
| Naive dyadic regression | How are measured country and country-pair characteristics associated with the number of actions from one country toward another if country-pairs are treated as conditionally independent? | Provides a familiar benchmark and exposes what the independence assumption buys | It uses ordered country-pairs as rows and does not represent shared-country or broader network dependence | “If the outcome is a count of cooperative statements, a coefficient is the expected change in that count for one ordered country-pair when the predictor rises by one unit, holding the listed covariates constant. We should not yet trust its interval because one country contributes a row for every potential partner.” |
| Country-clustered uncertainty preview | How much does recognizing that the same country initiates actions toward many partners change uncertainty in that benchmark? | Recalculates uncertainty under one limited dependence pattern | It changes the covariance estimate, not the fitted mean or the substantive dependence model | “The estimated change in the number of cooperative statements stays exactly the same. Its interval changes after we tell the calculation that all actions initiated by the same country may share unmeasured causes. This still does not account for repeated targets, reciprocal action, or dependence involving third countries.” |
The question at the end of every code block is therefore not only “Did it run?” It is “What social fact did this object or estimator let us describe, and what did it leave unresolved?”
Let’s start abstract, because not every project in this course studies countries.
One more convention: the collapsible boxes titled “the full toolbox” are self-study. They contain the rest of the netify surface in a form you can run at home, and nothing in the main path depends on them.
2 A tie is a decision, not a cell in a spreadsheet
Strip away the substance and every network in this course is the same object: a set of actors, and a set of ties running between them. A tie has a sender (\(i\)) and a receiver (\(j\)). That is the whole vocabulary. It does not matter whether \(i\) and \(j\) are students who named each other as friends, firms that traded, brain regions that co-activated, or countries that signed a treaty: the machinery is identical, and we will say “sender” and “receiver” all week so nobody has to translate out of political science.
The smallest useful example is a classroom. netify ships one: thirty students, and a list of who named whom as a friend. Friendship nomination is directed: you can name someone who does not name you back, so this is a directed network, and it is a good place to meet the mechanics before we touch anything with a foreign ministry in it.
classroom_edges <- load_day8_data("classroom_edges", day8_dir)
classroom_nodes <- load_day8_data("classroom_nodes", day8_dir)
head(classroom_edges)| from | to |
|---|---|
| s01 | s03 |
| s03 | s07 |
| s04 | s07 |
| s06 | s07 |
| s01 | s08 |
| s06 | s08 |
This is an edgelist: one row per tie, a from and a to. This is the rawest form network data comes in, and it is what netify() most wants to eat. The nodes come with attributes too, including gender, grade, and GPA, which we will attach in a moment.
class_net <- netify(
classroom_edges,
actor1 = "from", # The sender column
actor2 = "to", # The receiver column
symmetric = FALSE # Nomination is directed: i -> j is not j -> i
)
class_netThat print-out is worth reading slowly, because it is netify telling you what it thinks you just built. It is unipartite (one kind of actor, not students-and-clubs), asymmetric (directed, because we said so), binary (a tie either exists or it does not; there are no weights yet), and cross-sectional (one snapshot, no time). Thirty actors. Those four facts are the identity of the object, and every one of them was a decision you made or accepted a default on. A small Rosetta note, because the labels on screen are not always the words we say: netify prints Asymmetric (= directed), Binary Weights (= unweighted), and Cross-Sectional (= one time point).
WarningThe message about
missing_to_zero is a modeling decision in disguise
netify mentions that missing dyads were filled with zeros. For describing a classroom that is usually fine: if two students are not listed, they are not friends. But “we never observed this pair” and “this pair is definitely not tied” are different claims about the world, and three arguments control which one you are making: missing_to_zero, diag_to_NA, and nodelist. They are cheap to see on a toy, so let’s not take them on faith:
toy <- data.frame(
from = c("ana", "ana", "ben", "cal"),
to = c("ben", "cal", "ana", "ana"),
mins = c(20, 5, 12, 3) # Minutes talked, say
)
# The defaults: unobserved pairs become 0, the diagonal becomes NA
toy_zero <- netify(toy, actor1 = "from", actor2 = "to",
symmetric = FALSE, weight = "mins")
peek(toy_zero, from = 1:3, to = 1:3)#> ana ben cal
#> ana NA 20 5
#> ben 12 NA 0
#> cal 3 0 NA# The preserve-uncertainty version: unobserved pairs stay NA, and
# `nodelist` adds dee, an actor we know exists but who has no rows
toy_na <- netify(toy, actor1 = "from", actor2 = "to",
symmetric = FALSE, weight = "mins",
missing_to_zero = FALSE,
nodelist = c("ana", "ben", "cal", "dee"))
peek(toy_na, from = 1:4, to = 1:4)#> ana ben cal dee
#> ana NA 20 5 NA
#> ben 12 NA NA NA
#> cal 3 NA NA NA
#> dee NA NA NA NA(peek() gets a proper introduction in the next section: it just shows you a corner of the matrix.) Three things moved. Ben→cal flipped from 0 to NA: we never saw that pair, and the default zero was a claim, not an observation. The diagonal is NA in both because diag_to_NA = TRUE by default: self-ties are usually meaningless, and you flip it for the networks where they are not. And dee, an isolate with no rows in the edgelist, exists at all only because nodelist put her there: without it her non-ties would silently vanish from every density you compute. On a later day the model will believe whichever claim you encoded, so make these three deliberately.
3 Reading the object
You now have a network. Before the summary tools, one flat truth to keep the object from ever feeling like a black box: a netify object is the adjacency matrix, with bookkeeping attached. peek() shows you any corner of it without drowning the console in a thirty-by-thirty print:
peek(class_net, from = 1:5, to = 1:5)#> s01 s02 s03 s04 s05
#> s01 NA 0 1 0 0
#> s02 0 NA 0 0 0
#> s03 0 0 NA 0 0
#> s04 0 0 0 NA 0
#> s05 0 0 0 0 NARows are senders, columns are receivers, a 1 is a nomination, and the diagonal is NA because you cannot name yourself. Any time an object confuses you this week, peek() first: it is the fastest way to see what you actually built.
Now the two functions that tell you almost everything else, both of which you will use every single day of this course.
summary() describes the network as a whole:
summary(class_net)| net | num_actors | density | num_edges | prop_edges_missing | competition_row | competition_col | sd_of_row_means | sd_of_col_means | covar_of_row_col_means | reciprocity | mutual | transitivity |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 30 | 0.0586207 | 51 | 0 | 0.0595925 | 0.0549789 | 0.0529192 | 0.048046 | -0.4628256 | -0.0622711 | 0 | 0.0576923 |
Read a few of those off. The density is 0.059. Of all the directed pairs that could be friends, that fraction actually are, so this is a sparse network, as classrooms are. Reciprocity here is -0.062, and the sign matters: netify reports reciprocity as the correlation between the \(i\to j\) and \(j\to i\) tie values, so it runs from \(-1\) to \(+1\), not \(0\) to \(1\). At -0.06 it is essentially zero: if anything a hair below chance. In fact not a single one of these nominations is returned: summary(class_net)$mutual is exactly 0. That these students name friends who do not name them back is already a substantive feature of this network.
summary() prints a dozen more columns; ignore them for now: density and reciprocity are all you need today, and the rest arrive when the models that use them do.
summary_actor() drops one level down and describes each actor:
sa_class <- summary_actor(class_net)
head(sa_class[, c("actor", "degree_in", "degree_out", "degree_total")])| actor | degree_in | degree_out | degree_total |
|---|---|---|---|
| s01 | 0 | 4 | 4 |
| s02 | 0 | 1 | 1 |
| s03 | 1 | 3 | 4 |
| s04 | 0 | 3 | 3 |
| s05 | 0 | 4 | 4 |
| s06 | 0 | 5 | 5 |
degree_out is how many people a student named; degree_in is how many named them. Those are different questions: being friendly and being popular are not the same, and the fact that a directed network lets you separate them is exactly the kind of structure we will spend the week modeling. The most “popular” student here receives 5 nominations while the busiest nominator sends out 5.
Now attach the node attributes, so those degrees can be explained rather than just ranked:
class_net <- add_node_vars(
class_net,
classroom_nodes,
actor = "student", # The id column in classroom_nodes
node_vars = c("gender", "grade", "gpa")
)
class_netThe object now reports nodal features: gender, grade, and GPA are riding along with each actor. When we get to modeling, “do students befriend others of the same gender?” becomes a coefficient, and it will need those attributes to be attached in exactly this way.
Finally, look at it. A picture of a thirty-node network is not a rigorous object, but it is how you catch the gross mistakes: an isolated component you did not expect, one actor tied to everyone, a direction reversed. And because the attributes are attached, the picture can carry them: any nodal feature can drive any aesthetic through the node_*_by arguments.
plot(class_net, node_color_by = "gender", node_size_by = "gpa")
TipCheckpoint: you can already do the core loop
If you can read the four lines of that print-out (Unipartite / Asymmetric / Binary Weights / Cross-Sectional: netify’s own words for one-mode / directed / unweighted / one time point), run summary() and summary_actor(), and attach node attributes, you have the whole descriptive loop. Everything from here is the same loop on messier, more interesting data.
NoteThe full toolbox: know your object (self-study)
The debugging vocabulary for when an object you built is not what you thought it was. None of this is analysis; all of it saves you at 11pm.
# What am I holding? Seven predicates ...
c(netify = is_netify(class_net),
binary = is_binary(class_net),
longitudinal = is_longitudinal(class_net),
multilayer = is_multilayer(class_net),
bipartite = is_bipartite_netify(class_net),
directed = is_directed_netify(class_net),
symmetric = is_symmetric_netify(class_net))#> netify binary longitudinal multilayer bipartite directed
#> TRUE TRUE FALSE FALSE FALSE TRUE
#> symmetric
#> FALSE# ... and three dimension counters (these work on ANY netify object,
# and they are the documented way to ask, not length())
c(actors = n_actors(class_net), periods = n_periods(class_net),
layers = n_layers(class_net))#> actors periods layers
#> 30 1 1# The attached node attributes, back as a data.frame
head(nodal_data(class_net), 3)| actor | gender | grade | gpa |
|---|---|---|---|
| s01 | F | 12 | 2.66 |
| s02 | F | 9 | 2.89 |
| s03 | M | 12 | 2.81 |
# A full coherence check, run it when something feels off; the usual
# culprit it catches is an actor-name mismatch between your tables
all(unlist(validate_netify(class_net, verbose = FALSE)))#> [1] TRUE# And the escape hatch: it really is just a matrix
m <- as.matrix(class_net)
m[1:3, 1:3]#> s01 s02 s03
#> s01 NA 0 1
#> s02 0 NA 0
#> s03 0 0 NAget_raw() is the same escape hatch for any shape: on a longitudinal object it hands you the full actor × actor × time array.
4 A two-mode interlude: bipartite networks
Everything so far has one kind of actor. Plenty of you have data with two kinds, such as students and clubs, people and events, firms and boards, or words and documents, where ties run only between the modes and never within them. That is a bipartite network, and netify builds it the moment you say so:
affil <- data.frame(
student = c("sam", "sam", "lena", "lena", "omar", "omar", "noor"),
club = c("chess", "debate", "debate", "robotics", "chess",
"robotics", "debate")
)
bip_net <- netify(affil, actor1 = "student", actor2 = "club",
mode = "bipartite")
bip_netThe print-out now says bipartite, and rows and columns are different kinds of thing:
peek(bip_net, from = 1:4, to = 1:3)#> chess debate robotics
#> lena 0 1 1
#> noor 0 1 0
#> omar 1 0 1
#> sam 1 1 0The matrix is a rectangle, not a square: rows are students, and columns are clubs. The familiar “who shares a club with whom” projection is a choice you make later, not something the object does behind your back. A rectangle also reads naturally as tiles:
plot(bip_net, style = "heatmap")
We stay unipartite for the rest of the week, but if your data are two-mode you are not locked out of the modeling half of this course: Day 10’s AME machinery has a bipartite variant, and this object is exactly what it eats.
5 From events to a network
Classrooms come pre-cleaned. Observed data almost never does. A common raw form in the social sciences is event data: a long log of things that happened between pairs of actors, one row per event, and it is not a network until you decide how to turn events into ties.
netify ships an example from ICEWS, a dataset of interactions between countries coded from news reports, with one row per directed country-pair per year. “Country” is just the substantive name for “actor” here: if you study firms or people, read i and j as your senders and receivers.
icews <- load_day8_data("icews", day8_dir)
dim(icews)#> [1] 298376 24head(icews[, c("i", "j", "year", "verbCoop", "matlConf", "i_polity2")])| i | j | year | verbCoop | matlConf | i_polity2 | |
|---|---|---|---|---|---|---|
| 2 | Afghanistan | Albania | 2002 | 6 | 0 | NA |
| 3 | Afghanistan | Albania | 2003 | 1 | 0 | NA |
| 4 | Afghanistan | Albania | 2004 | 10 | 1 | NA |
| 5 | Afghanistan | Albania | 2005 | 0 | 0 | NA |
| 6 | Afghanistan | Albania | 2006 | 6 | 21 | NA |
| 7 | Afghanistan | Albania | 2007 | 3 | 0 | NA |
Nearly three hundred thousand rows. Each has a sender i, a receiver j, a year, and several counts: verbCoop is how many cooperative statements \(i\) made toward \(j\) that year, matlConf how many materially conflictual acts, and so on. There are also attributes measured on the sender: i_polity2 (a regime-type score), i_log_gdp, i_log_pop.
We are not going to model a hundred and fifty countries today: big networks plot badly and distract from the mechanics. So pick a coherent set of major players. This is the first substantive actor-handling decision, and it is yours to make explicitly:
keep <- c(
"United States", "China", "Russian Federation",
"United Kingdom", "France", "Germany", "Japan", "India",
"Brazil", "Israel", "Iran, Islamic Republic Of", "Iraq",
"Egypt", "Saudi Arabia", "Pakistan", "Afghanistan",
"Korea, Republic Of", "Syrian Arab Republic"
)
ic <- icews[icews$i %in% keep & icews$j %in% keep, ]
c(rows = nrow(ic), actors = length(unique(ic$i)),
years = length(unique(ic$year)))#> rows actors years
#> 3978 18 13
NoteNote the country names are exact strings, and that bites everyone once
"Iran, Islamic Republic Of", not "Iran". "Korea, Republic Of", not "South Korea". If you type the friendly name, netify will not error: it will simply build a network that silently omits Iran, and you may not notice until much later. Actor identity is string identity, and reconciling names across sources (ISO codes, COW codes, UN names) is an important and under-taught part of this work. When in doubt, sort(unique(c(ICEWS$I, ICEWS$j))) and read the list.
Now the two decisions that define the tie. Is it directed, and is it weighted?
ic14 <- ic[ic$year == 2014, ]
net14 <- netify(
ic14,
actor1 = "i",
actor2 = "j",
symmetric = FALSE, # Cooperation i->j need not equal j->i
weight = "verbCoop" # The tie carries a count, not just 0/1
)
net14Two arguments carried all the meaning there:
symmetric = FALSE. The United States can direct far more cooperative signals at a small state than it gets back. If you setsymmetric = TRUEon directed data,netifymirrors one triangle of the matrix onto the other and throws half your information away: quietly. So the rule is: symmetry is a claim about the world, not a convenience. Set it to match how the data were generated.weight = "verbCoop". With no weight you get a binary tie: did \(i\) ever cooperate with \(j\)? With a weight, the tie carries how much. Which you want depends on your question, and it changes the model family downstream (a count is not a coin flip). Here we keep the count.
The density is now 0.987: among these major powers essentially everyone signals everyone, which is exactly why a binary network of great-power cooperation would be boring and a weighted one is not.
6 sum_dyads: when a pair shows up more than once
Here is the argument people skip and then get wrong. We just built a single year. But netify needs one tie per pair, so what happens when the raw data has the same pair on multiple rows? That is not a corner case; it is what all event data looks like before you aggregate it.
Watch what happens if we hand netify all thirteen years at once but do not tell it about time. Each country-pair now appears up to thirteen times:
pair <- paste(ic$i, ic$j)
table(table(pair)) # How many rows does each directed pair get?#>
#> 13
#> 306Every pair shows up thirteen times. netify cannot put thirteen numbers in one cell of the matrix, so it makes you decide with sum_dyads:
net_pooled <- netify(
ic,
actor1 = "i",
actor2 = "j",
symmetric = FALSE,
weight = "verbCoop",
sum_dyads = TRUE # Collapse the repeated rows by summing the weight
)
net_pooledsum_dyads = TRUE says “these repeated rows are pieces of one tie: add them up.” Now the weight on the tie from the United States to China is total verbal cooperation from 2002 to 2014. That is a legitimate object: a pooled, cross-sectional summary of a decade. The mean tie weight jumps to 7106, roughly thirteen times the single-year figure, exactly because we stacked thirteen years into it.
Important
sum_dyads is the difference between “a tie” and “an average tie”
Whenever more than one row maps to the same pair, you must decide what the tie means: the sum of the events (sum_dyads = TRUE) or, if you aggregate by hand first, their average, maximum, or most recent value. netify will not guess for you, and the default (FALSE) assumes each pair is already unique. In netify 1.5.3, repeated weighted dyads trigger a warning and an automatic switch to sum_dyads = TRUE, so the rows are summed rather than dropped. Do not rely on that fallback: make the choice explicit so your code says which world you are in.
NoteThe full toolbox: event data before it becomes a network (self-study)
sum_dyads = TRUE aggregates inside netify(). Sometimes you want to aggregate first and inspect the result before you commit to a network. aggregate_dyad() is the standalone version of the same move. netify ships UCDP records of organized violence in Mexico for practicing it: one row per recorded incident, with a best estimate of deaths.
mexico <- load_day8_data("mexico", day8_dir)
mex_agg <- aggregate_dyad(mexico, actor1 = "side_a", actor2 = "side_b",
time = "year", weight = "best",
symmetric = FALSE)
head(mex_agg, 3)| side_a | side_b | year | best |
|---|---|---|---|
| Government of Mexico | Civilians | 1989 | 3 |
| Government of Mexico | Civilians | 1990 | 4 |
| Government of Mexico | Civilians | 1992 | 1 |
mex_net <- netify(mex_agg, actor1 = "side_a", actor2 = "side_b",
time = "year", symmetric = FALSE, weight = "best")
mex_net18,957 raw incident rows collapse to 313 dyad-years before the network exists, and you got to read the aggregation on the way through, which is exactly what you want when the aggregation rule itself is a modeling decision. Should you use the sum, the maximum, or the most recent value? One companion utility: when the aggregation is undirected (symmetric = TRUE treats i→j and j→i as the same pair), gen_symm_id() builds the order-free pair id that makes that work: both orderings of a pair get the same id:
ex <- data.frame(i = c("usa", "china"), j = c("china", "usa"))
gen_symm_id(ex, actor1 = "i", actor2 = "j") # Same id, both rows#> [1] "5:china|3:usa" "5:china|3:usa"(myanmar, a second UCDP event dataset, also ships with the package if you want a homework-sized playground.)
But pooling a decade into one number also throws away the time. Usually the time is the interesting part. So instead of collapsing the years, let’s keep them.
Break. You have built cross-sectional networks, learned what symmetric, weight, and sum_dyads each claim about the world, and met the caching pattern. After the break: time, covariates, and handing the object to a model.
7 Going longitudinal
Give netify a time column and it stops making one network and starts making one network per period, lined up so the same actor is the same row in every year. This is the object that days 11 and 12 need.
net_long <- netify(
ic,
actor1 = "i",
actor2 = "j",
time = "year", # <- The only new argument
symmetric = FALSE,
weight = "verbCoop"
)
net_longNo sum_dyads needed now: once time is a dimension, each pair is unique within a year. The object reports 13 periods (n_periods() is how you ask; n_actors() and n_layers() are its siblings), and summary() on a longitudinal object gives you one row per period. But great-power cooperation is nearly saturated every single year: its density sits near one and barely budges, so it is the wrong place to watch a network move. Read the same countries on material conflict instead:
net_long_conf <- netify(
ic, actor1 = "i", actor2 = "j", time = "year",
symmetric = FALSE, weight = "matlConf" # Conflict, not cooperation
)
sl <- summary(net_long_conf)
sl[, c("net", "num_actors", "density", "mean_edge_weight")]| net | num_actors | density | mean_edge_weight |
|---|---|---|---|
| 2002 | 18 | 0.6830065 | 59.21053 |
| 2003 | 18 | 0.7450980 | 99.46053 |
| 2004 | 18 | 0.7679739 | 78.05106 |
| 2005 | 18 | 0.7483660 | 46.76856 |
| 2006 | 18 | 0.7712418 | 47.05508 |
| 2007 | 18 | 0.7777778 | 62.21849 |
| 2008 | 18 | 0.7549020 | 52.90909 |
| 2009 | 18 | 0.7647059 | 51.22650 |
| 2010 | 18 | 0.7385621 | 53.02212 |
| 2011 | 18 | 0.7777778 | 52.45378 |
| 2012 | 18 | 0.7679739 | 55.08085 |
| 2013 | 18 | 0.7483660 | 62.26638 |
| 2014 | 18 | 0.7777778 | 70.02101 |
Now the panel remembers the news. Density climbs from 0.68 in 2002 to 0.78 by 2014, and the mean edge weight, the average amount of conflict carried on an active tie, jumps from 59 in 2002 to 99 in 2003. Watch 2003: that spike is the Iraq invasion landing in the object. A single cross-section would have hidden it; the longitudinal object lets a historical event show up as a number.
And because the object knows it is a panel, the tools you already have turn into figures for free. plot() facets by period, filtered here to four years so each panel remains legible, and plot() on the summary table turns its columns into time series:
plot(net_long_conf,
time_filter = as.character(c(2002, 2003, 2009, 2014)),
mutate_weight = log1p) # Event counts are skewed; log for readability
plot(sl)
One more question a panel invites: how much does the network change year to year? compare_networks() can return several pairwise comparisons at once. Here we request method = "all", so it computes the ordinary edge-pattern correlations we display below and also runs a quadratic assignment procedure (QAP) test that permutes actor labels jointly across rows and columns:
cn <- compare_networks(net_long_conf, method = "all",
n_permutations = 500, seed = 6886,
return_details = TRUE)
cn#> metric mean sd min max
#> 1 correlation 0.564 0.243 0.079 0.969
#> 2 jaccard 0.746 0.036 0.647 0.819
#> 3 hamming 0.220 0.033 0.154 0.307
#> 4 qap_correlation 0.564 0.243 0.079 0.969
#> 5 spectral 3444.400 3395.200 270.400 1.11e+04#> 2002 2003 2004 2005 2006 2007 2008
#> 2002 1.0000000 0.93150772 0.89333704 0.8499138 0.8816440 0.6879525 0.6183963
#> 2003 0.9315077 1.00000000 0.96853173 0.8321313 0.8453982 0.6874325 0.5243189
#> 2004 0.8933370 0.96853173 1.00000000 0.8841220 0.8695235 0.6930095 0.5668536
#> 2005 0.8499138 0.83213131 0.88412198 1.0000000 0.9171521 0.7052095 0.7019919
#> 2006 0.8816440 0.84539819 0.86952354 0.9171521 1.0000000 0.8007490 0.8221283
#> 2007 0.6879525 0.68743252 0.69300949 0.7052095 0.8007490 1.0000000 0.6656050
#> 2008 0.6183963 0.52431885 0.56685364 0.7019919 0.8221283 0.6656050 1.0000000
#> 2009 0.3558954 0.19798673 0.24403809 0.4523788 0.5972208 0.5042312 0.8645049
#> 2010 0.3606311 0.19545823 0.24226101 0.5243762 0.6205780 0.5439233 0.8202507
#> 2011 0.2599259 0.14431782 0.17969735 0.3511321 0.5129090 0.4280943 0.7964822
#> 2012 0.2345353 0.09436527 0.13460750 0.4018625 0.5076933 0.4522892 0.7825113
#> 2013 0.1783792 0.09181299 0.07931633 0.2544302 0.3188623 0.2717077 0.6029515
#> 2014 0.5986396 0.59029948 0.57465955 0.6224615 0.6689281 0.5714430 0.7200738
#> 2009 2010 2011 2012 2013 2014
#> 2002 0.3558954 0.3606311 0.2599259 0.23453534 0.17837919 0.5986396
#> 2003 0.1979867 0.1954582 0.1443178 0.09436527 0.09181299 0.5902995
#> 2004 0.2440381 0.2422610 0.1796974 0.13460750 0.07931633 0.5746595
#> 2005 0.4523788 0.5243762 0.3511321 0.40186247 0.25443022 0.6224615
#> 2006 0.5972208 0.6205780 0.5129090 0.50769331 0.31886232 0.6689281
#> 2007 0.5042312 0.5439233 0.4280943 0.45228919 0.27170771 0.5714430
#> 2008 0.8645049 0.8202507 0.7964822 0.78251129 0.60295154 0.7200738
#> 2009 1.0000000 0.8457467 0.8013030 0.84964092 0.48784461 0.5211909
#> 2010 0.8457467 1.0000000 0.8395894 0.88364354 0.53624371 0.5015491
#> 2011 0.8013030 0.8395894 1.0000000 0.81431285 0.60281087 0.4476392
#> 2012 0.8496409 0.8836435 0.8143128 1.00000000 0.56893696 0.4548003
#> 2013 0.4878446 0.5362437 0.6028109 0.56893696 1.00000000 0.5425039
#> 2014 0.5211909 0.5015491 0.4476392 0.45480032 0.54250391 1.0000000#> 2002 2003 2004 2005 2006 2007
#> 2002 NA 0.001996008 0.001996008 0.001996008 0.001996008 0.001996008
#> 2003 0.001996008 NA 0.001996008 0.001996008 0.001996008 0.001996008
#> 2004 0.001996008 0.001996008 NA 0.001996008 0.001996008 0.001996008
#> 2005 0.001996008 0.001996008 0.001996008 NA 0.001996008 0.001996008
#> 2006 0.001996008 0.001996008 0.001996008 0.001996008 NA 0.001996008
#> 2007 0.001996008 0.001996008 0.001996008 0.001996008 0.001996008 NA
#> 2008 0.001996008 0.001996008 0.001996008 0.001996008 0.001996008 0.001996008
#> 2009 0.019960080 0.025948104 0.019960080 0.013972056 0.001996008 0.001996008
#> 2010 0.011976048 0.021956088 0.033932136 0.003992016 0.003992016 0.001996008
#> 2011 0.013972056 0.033932136 0.029940120 0.013972056 0.003992016 0.005988024
#> 2012 0.023952096 0.041916168 0.037924152 0.003992016 0.001996008 0.001996008
#> 2013 0.025948104 0.031936128 0.057884232 0.013972056 0.013972056 0.013972056
#> 2014 0.001996008 0.001996008 0.001996008 0.001996008 0.001996008 0.001996008
#> 2008 2009 2010 2011 2012 2013
#> 2002 0.001996008 0.019960080 0.011976048 0.013972056 0.023952096 0.025948104
#> 2003 0.001996008 0.025948104 0.021956088 0.033932136 0.041916168 0.031936128
#> 2004 0.001996008 0.019960080 0.033932136 0.029940120 0.037924152 0.057884232
#> 2005 0.001996008 0.013972056 0.003992016 0.013972056 0.003992016 0.013972056
#> 2006 0.001996008 0.001996008 0.003992016 0.003992016 0.001996008 0.013972056
#> 2007 0.001996008 0.001996008 0.001996008 0.005988024 0.001996008 0.013972056
#> 2008 NA 0.001996008 0.001996008 0.001996008 0.001996008 0.001996008
#> 2009 0.001996008 NA 0.001996008 0.001996008 0.001996008 0.001996008
#> 2010 0.001996008 0.001996008 NA 0.001996008 0.001996008 0.001996008
#> 2011 0.001996008 0.001996008 0.001996008 NA 0.001996008 0.001996008
#> 2012 0.001996008 0.001996008 0.001996008 0.001996008 NA 0.001996008
#> 2013 0.001996008 0.001996008 0.001996008 0.001996008 0.001996008 NA
#> 2014 0.001996008 0.001996008 0.001996008 0.001996008 0.001996008 0.001996008
#> 2014
#> 2002 0.001996008
#> 2003 0.001996008
#> 2004 0.001996008
#> 2005 0.001996008
#> 2006 0.001996008
#> 2007 0.001996008
#> 2008 0.001996008
#> 2009 0.001996008
#> 2010 0.001996008
#> 2011 0.001996008
#> 2012 0.001996008
#> 2013 0.001996008
#> 2014 NA
#> attr(,"n_perm")
#> [1] 500#> 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014
#> 2002 NA 500 500 500 500 500 500 500 500 500 500 500 500
#> 2003 500 NA 500 500 500 500 500 500 500 500 500 500 500
#> 2004 500 500 NA 500 500 500 500 500 500 500 500 500 500
#> 2005 500 500 500 NA 500 500 500 500 500 500 500 500 500
#> 2006 500 500 500 500 NA 500 500 500 500 500 500 500 500
#> 2007 500 500 500 500 500 NA 500 500 500 500 500 500 500
#> 2008 500 500 500 500 500 500 NA 500 500 500 500 500 500
#> 2009 500 500 500 500 500 500 500 NA 500 500 500 500 500
#> 2010 500 500 500 500 500 500 500 500 NA 500 500 500 500
#> 2011 500 500 500 500 500 500 500 500 500 NA 500 500 500
#> 2012 500 500 500 500 500 500 500 500 500 500 NA 500 500
#> 2013 500 500 500 500 500 500 500 500 500 500 500 NA 500
#> 2014 500 500 500 500 500 500 500 500 500 500 500 500 NA
#> attr(,"n_perm")
#> [1] 500cor_mat <- cn$details$correlation_matrix
c(`2002 vs 2003` = round(cor_mat["2002", "2003"], 2),
`2002 vs 2014` = round(cor_mat["2002", "2014"], 2))#> 2002 vs 2003 2002 vs 2014
#> 0.93 0.60The 2002-2003 pair correlates at 0.93; across the full twelve-year gap that drops to 0.6. The conflict network has memory, although one adjacent pair is not enough to establish a smooth rate of decay, which is exactly the structure Day 12 (SAOM) turns into a model instead of a description.
NoteA reference for comparing networks
compare_networks() can compare four different things, and the choice should follow the question. what = "edges" compares the dyad cells themselves; what = "structure" compares graph-level summaries such as density and reciprocity; what = "nodes" checks who entered or left; and what = "attributes" compares the distributions of attached nodal attributes. For a longitudinal object it compares all pairs of waves, not only adjacent waves.
The edge metrics are not interchangeable:
- Correlation asks whether the two sets of edge weights rise and fall together. Shared zeros count, so a sparse pair of networks can correlate highly even when the positive ties overlap only modestly.
- Jaccard thresholds the weights and asks what fraction of ties present in at least one network appear in both. For a weighted network the default
edge_threshold = 0means any positive weight is a tie. - Hamming reports the fraction of eligible dyads whose thresholded tie status differs.
- QAP compares the observed Pearson edge correlation with a reference distribution produced by relabeling actors jointly across rows and columns. It is a test against that relabeling null, not a correction for omitted variables and not a causal design.
The same object can answer the structure and composition questions without being rebuilt:
structure_over_time <- compare_networks(net_long_conf, what = "structure")
actors_over_time <- compare_networks(net_long_conf, what = "nodes")
structure_over_time$summary| network | num_actors | density | num_edges | prop_edges_missing | mean_edge_weight | sd_edge_weight | median_edge_weight | min_edge_weight | max_edge_weight | competition_row | competition_col | sd_of_row_means | sd_of_col_means | covar_of_row_col_means | reciprocity | mutual | transitivity | mean_degree |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2002 | 18 | 0.6830065 | 209 | 0 | 59.21053 | 324.1937 | 7.0 | 1 | 4258 | 0.3070161 | 0.2029327 | 88.53331 | 67.77771 | 0.1973468 | 0.1627361 | 0.7416667 | 0.8287938 | 11.61111 |
| 2003 | 18 | 0.7450980 | 228 | 0 | 99.46053 | 753.6213 | 8.0 | 1 | 11173 | 0.3742165 | 0.3679582 | 182.63191 | 180.82962 | 0.2722463 | 0.2835956 | 0.7272727 | 0.8827362 | 12.66667 |
| 2004 | 18 | 0.7679739 | 235 | 0 | 78.05106 | 439.1223 | 10.0 | 1 | 6222 | 0.2640068 | 0.2329370 | 119.47459 | 110.21162 | 0.7360668 | 0.6533335 | 0.7803030 | 0.8820652 | 13.05556 |
| 2005 | 18 | 0.7483660 | 229 | 0 | 46.76856 | 151.1109 | 9.0 | 1 | 1755 | 0.1692152 | 0.1126077 | 51.51328 | 36.49658 | 0.8106576 | 0.7226663 | 0.8174603 | 0.8492582 | 12.72222 |
| 2006 | 18 | 0.7712418 | 236 | 0 | 47.05508 | 134.0385 | 11.0 | 1 | 1614 | 0.1648921 | 0.1117831 | 52.38751 | 37.56816 | 0.7268191 | 0.6537191 | 0.7611940 | 0.8985812 | 13.11111 |
| 2007 | 18 | 0.7777778 | 238 | 0 | 62.21849 | 202.3437 | 10.5 | 1 | 2031 | 0.1862538 | 0.1041586 | 76.37616 | 46.57518 | 0.5572080 | 0.3917168 | 0.8167939 | 0.8823529 | 13.22222 |
| 2008 | 18 | 0.7549020 | 231 | 0 | 52.90909 | 150.2178 | 10.0 | 1 | 1031 | 0.1822549 | 0.1164830 | 62.06632 | 43.04032 | 0.7283532 | 0.7093538 | 0.7906977 | 0.8707483 | 12.83333 |
| 2009 | 18 | 0.7647059 | 234 | 0 | 51.22650 | 175.9671 | 9.0 | 1 | 1553 | 0.1574567 | 0.1599011 | 54.59171 | 55.24261 | 0.7930298 | 0.7690522 | 0.7727273 | 0.8834056 | 13.00000 |
| 2010 | 18 | 0.7385621 | 226 | 0 | 53.02212 | 135.5804 | 9.0 | 1 | 1018 | 0.1547442 | 0.1077054 | 53.84223 | 39.04081 | 0.7684170 | 0.7558204 | 0.8080000 | 0.8426763 | 12.55556 |
| 2011 | 18 | 0.7777778 | 238 | 0 | 52.45378 | 141.0171 | 10.0 | 1 | 1392 | 0.1463712 | 0.1093430 | 53.67359 | 41.30675 | 0.8105601 | 0.7108185 | 0.7500000 | 0.9035403 | 13.22222 |
| 2012 | 18 | 0.7679739 | 235 | 0 | 55.08085 | 142.6921 | 13.0 | 1 | 948 | 0.1230319 | 0.1090987 | 47.97012 | 42.73139 | 0.7444379 | 0.7665349 | 0.8076923 | 0.8724234 | 13.05556 |
| 2013 | 18 | 0.7483660 | 229 | 0 | 62.26638 | 194.2963 | 9.0 | 1 | 1703 | 0.1627615 | 0.1478441 | 66.60781 | 61.80017 | 0.2931914 | 0.3009831 | 0.7890625 | 0.8681948 | 12.72222 |
| 2014 | 18 | 0.7777778 | 238 | 0 | 70.02101 | 214.6127 | 10.0 | 1 | 1870 | 0.1614356 | 0.1187480 | 77.36401 | 59.76741 | 0.4480364 | 0.5134052 | 0.8167939 | 0.8821586 | 13.22222 |
actors_over_time$summary| network | n_nodes | mean_overlap | mean_jaccard | |
|---|---|---|---|---|
| 2002 | 2002 | 18 | 18 | 1 |
| 2003 | 2003 | 18 | 18 | 1 |
| 2004 | 2004 | 18 | 18 | 1 |
| 2005 | 2005 | 18 | 18 | 1 |
| 2006 | 2006 | 18 | 18 | 1 |
| 2007 | 2007 | 18 | 18 | 1 |
| 2008 | 2008 | 18 | 18 | 1 |
| 2009 | 2009 | 18 | 18 | 1 |
| 2010 | 2010 | 18 | 18 | 1 |
| 2011 | 2011 | 18 | 18 | 1 |
| 2012 | 2012 | 18 | 18 | 1 |
| 2013 | 2013 | 18 | 18 | 1 |
| 2014 | 2014 | 18 | 18 | 1 |
Before interpreting any comparison, check that the networks use the same actor identities, direction, observation window, weight scale, and missing-data rule. An actor who did not exist is not the same as an actor with zero ties, and a one-month event count is not directly comparable to a one-year count.
Note
actor_time_uniform: who is allowed to appear and when
We got a clean panel because all eighteen of these countries appear in all thirteen years. Empirical panels are often ragged: states enter, firms merge, and people drop out of a survey wave. netify handles this with actor_time_uniform. Left at its default it assumes your actor set is the same every period; set it FALSE and composition is inferred from when each actor appears in the data, treating an absent actor as absent rather than as present-with-no-ties. Those are different claims, and a toy makes the difference visible: dee only shows up in year 2:
toy_l <- data.frame(
i = c("ana", "ben", "ana", "ben", "dee", "ana"),
j = c("ben", "ana", "cal", "cal", "ana", "dee"),
yr = c(1, 1, 1, 2, 2, 2))
uniform <- netify(toy_l, actor1 = "i", actor2 = "j", time = "yr",
symmetric = FALSE) # The default
ragged <- netify(toy_l, actor1 = "i", actor2 = "j", time = "yr",
symmetric = FALSE, actor_time_uniform = FALSE)
summary(uniform)[, c("net", "num_actors", "density")]| net | num_actors | density |
|---|---|---|
| 1 | 4 | 0.25 |
| 2 | 4 | 0.25 |
summary(ragged)[, c("net", "num_actors", "density")]| net | num_actors | density |
|---|---|---|
| 1 | 3 | 0.50 |
| 2 | 4 | 0.25 |
peek(uniform, from = 1:4, to = 1:4, time = 1) # dee: a row of zeros#> $`1`
#> ana ben cal dee
#> ana NA 1 1 0
#> ben 1 NA 0 0
#> cal 0 0 NA 0
#> dee 0 0 0 NApeek(ragged, from = 1:4, to = 1:4, time = 1) # dee: not there at all#> $`1`
#> ana ben cal
#> ana NA 1 1
#> ben 1 NA 0
#> cal 0 0 NASame data, different year-1 density, because “Dee sent no ties” and “Dee did not exist yet” are different denominators. If the automatic inference is not the composition you want, actor_pds lets you hand netify an explicit actor-by-period record instead. The model downstream will believe whichever claim you encoded; Day 12 (SAOM) leans on this directly when actors join and leave between waves.
Now the caching pattern, once, on the one object that takes more than a blink. Building the network for all one-hundred-fifty-odd countries across all thirteen years requires appreciable computation:
net_world <- cache_fit("net_world", netify(
icews,
actor1 = "i",
actor2 = "j",
time = "year",
symmetric = FALSE,
weight = "verbCoop",
nodal_vars = c("i_polity2", "i_log_gdp", "i_log_pop")
), dir = file.path(day8_dir, "cache"))
c(actors = n_actors(net_world), periods = n_periods(net_world))#> actors periods
#> 152 13The call stays right there on the screen: there is no hidden second version that ran. cache_fit() hashes that code and saves the result to the Day 8 cache/net_world.rda file, then reloads it in milliseconds until you change the code. At that point it notices and rebuilds, which is how every slow model fit later in the week ships without stalling a live class.
Once the validated world-sized object exists, you rarely rebuild from raw data again: you carve. subset_netify() takes any combination of actors, periods, and layers and returns a new netify object with the covariates still attached:
net_recent <- subset_netify(net_world,
actors = keep, # Our 18 majors
time = as.character(2010:2014)) # Last 5 years
net_recentThis is the pattern every later day’s data prep starts with: build wide once, then subset down to the estimation window you can defend. (Actor names are string-exact here too: misspell "Iran, Islamic Republic Of" and it is silently not in your subset.)
8 Attaching covariates
A network you only ever describe can stay bare. A network you want to model needs its covariates attached, and there are two straightforward ways to do it.
Inline, at build time: name the columns and let netify pull them:
net_cov <- netify(
ic14,
actor1 = "i",
actor2 = "j",
symmetric = FALSE,
weight = "verbCoop",
nodal_vars = c("i_polity2", "i_log_gdp"), # Measured on the actor
dyad_vars = "matlConf", # Measured on the pair
dyad_vars_symmetric = FALSE # Conflict i->j != j->i
)
net_covTwo kinds of covariate, and the distinction is the whole ballgame later today:
- Nodal covariates attach to an actor: regime type, GDP, population. One value per country. These get copied into every dyad that country is in, and that copying is precisely what will cause the trouble in part two.
- Dyadic covariates attach to a pair: material conflict from \(i\) to \(j\), distance between them, whether they share an alliance. One value per ordered pair, and
dyad_vars_symmetrictellsnetifywhether the \(i\to j\) value equals the \(j\to i\) value (distance: yes; conflict: no).
Or after the fact, when your covariates live in a separate table, which they usually do:
node_tab <- unique(ic14[, c("i", "i_polity2", "i_log_gdp")])
names(node_tab)[1] <- "country"
net_bare <- netify(ic14, actor1 = "i", actor2 = "j",
symmetric = FALSE, weight = "verbCoop")
net_bare <- add_node_vars(
net_bare, node_tab, actor = "country",
node_vars = c("i_polity2", "i_log_gdp")
)
net_bare <- add_dyad_vars(
net_bare, ic14, actor1 = "i", actor2 = "j",
dyad_vars = "matlConf", dyad_vars_symmetric = FALSE
)
net_bareSame object, either road. add_node_vars() and add_dyad_vars() matter because they do the alignment for you: they match on actor identity, not on row order, so a covariate can never silently end up attached to the wrong actor. If you have ever built a design matrix by cbind()-ing two tables that turned out to be sorted differently, you know why that guarantee is worth having.
NoteThe full toolbox: reshaping and rescuing objects (self-study)
Five utilities you will not need today but will want the week you meet messy data, plus the doors in from other formats:
# 1. Drop actors with missing covariates (returns the cleaned object)
drop_na_actors(net_cov)
# 2. Combine networks built separately, `align_actors` decides how
# the two actor sets get reconciled ("union" keeps everyone)
eu <- subset_netify(net14, actors = c("United Kingdom", "France", "Germany"))
asia <- subset_netify(net14, actors = c("China", "Japan", "India"))
bind_netifies(eu, asia, align_actors = "union")
# 3. Swap a dyadic covariate in as THE network, demoting the old
# network to a covariate, here conflict becomes the outcome
pivot_dyad_to_network(net_cov, dyad_var = "matlConf")
# 4. Decompose back into plain data.frames (edges + nodes)
dec <- decompose_netify(net_cov)
head(dec$edge_data, 3)| from | to | time | verbCoop | matlConf |
|---|---|---|---|---|
| Afghanistan | China | 1 | 749 | 18 |
| Afghanistan | Egypt | 1 | 4 | 0 |
| Afghanistan | France | 1 | 26 | 9 |
# 5. Arrived with an igraph or network object instead of an edgelist?
# to_netify() is the door in (round-tripped here for the demo)
to_netify(to_igraph(net_cov), weight = "weight")(In the bind_netifies() output, the cross-bloc cells are missing, not zero: the two source networks never observed each other’s dyads, which is missing_to_zero logic doing the right thing.) And if your collaborator hands you a Gephi or Pajek file, there are readers for those too:
read_graphml("your_file.graphml")
read_gml("your_file.gml")
read_pajek("your_file.net")9 Actor Statistics and Plotting with ICEWS
Back to summary_actor(), now on something with more going on than a classroom. On a directed, weighted network it gives you the full menu: in- and out-degree, strength (degree weighted by tie value), betweenness, closeness, hub and authority scores:
sa14 <- summary_actor(net14)
head(sa14[order(-sa14$strength_sum_out),
c("actor", "degree_total", "strength_sum_out",
"strength_sum_in")])| actor | degree_total | strength_sum_out | strength_sum_in | |
|---|---|---|---|---|
| 18 | United States | 34 | 33252 | 29330 |
| 3 | China | 34 | 18755 | 18610 |
| 8 | Iran, Islamic Republic Of | 34 | 17451 | 14027 |
| 14 | Russian Federation | 34 | 16441 | 16717 |
| 11 | Japan | 34 | 10681 | 11203 |
| 7 | India | 34 | 8996 | 9522 |
Sort by out-strength and you get the states that initiated the largest total number of cooperative statements toward other states. Sort by in-strength and you get the states that were the targets of the largest total number of cooperative statements. United States initiated the most in this network. These rankings can name different states because initiating diplomatic cooperation and being the target of it are different political actions. That distinction is the seed of the entire modeling week. Remember it for part two.
Now the picture. Two small moves come first because raw event counts are so right-skewed that a direct plot is one blinding U.S.-China edge and seventeen ghosts. mutate_weights() transforms the tie weights (the original column survives as a dyad feature, so nothing is lost), and add_node_vars() will attach anything, including the actor statistics you just computed, so the figure can be driven by the statistics:
net14_log <- mutate_weights(net14, transform_fn = log1p,
new_name = "log_verbCoop")
net14_log <- add_node_vars(net14_log, sa14, actor = "actor",
node_vars = "strength_sum_out")plot(net14_log, node_size_by = "strength_sum_out")
That is the whole aesthetic grammar: node_size_by, node_color_by, edge_alpha_by, remove_isolates, and a layout argument for the algorithm. This is the same mapping you used to color the classroom by gender, now fed by computed statistics instead of survey columns. For a weighted matrix there is also a more complete view than a node-link diagram: the heatmap shows every cell, including the empty ones, with no layout algorithm in between.
plot(net14_log, style = "heatmap")
The actor statistics themselves deserve a figure too: ranking tables do not survive a slide:
plot_actor_stats(sa14, across_actor = TRUE,
specific_stats = c("strength_sum_out", "strength_sum_in"))
And summary_actor() is extensible: hand any function of the adjacency matrix to other_stats and it becomes a column beside the built-ins. Here is the pattern:
sa_custom <- summary_actor(net14, other_stats = list(
n_partners_out = function(mat) rowSums(mat > 0, na.rm = TRUE)))
head(sa_custom[, c("actor", "degree_total", "n_partners_out")], 3)| actor | degree_total | n_partners_out |
|---|---|---|
| Afghanistan | 31 | 16 |
| Brazil | 31 | 16 |
| China | 34 | 17 |
Remember it when the centrality you need is not on the menu.
NoteDepth: plotting is its own craft, and
netify gives you the seams
plot() returns a ggplot object, so everything you know about ggplot2 applies: you can add layers, change scales, facet a longitudinal network by year. Under the hood netify computes a layout and hands it back through get_node_layout() / get_edge_layout() if you want to freeze coordinates across time so a longitudinal network does not jump around frame to frame. For a static teaching plot the defaults are fine; for a paper figure you will want that control.
When the high-level plot method almost does what you want, net_plot_data() exposes the two data frames behind it. The edge frame contains x1, y1, x2, and y2; the node frame contains x and y. From there it is ordinary ggplot2:
plot_parts <- net_plot_data(net14_log)$net_dfs
ggplot() +
geom_segment(
data = plot_parts$edge_data,
aes(x = x1, y = y1, xend = x2, yend = y2,
alpha = log_verbCoop),
colour = "grey45", linewidth = 0.35
) +
geom_point(
data = plot_parts$nodal_data,
aes(x = x, y = y, size = strength_sum_out),
colour = "#18453B"
) +
scale_alpha_continuous(range = c(0.05, 0.65), guide = "none") +
scale_size_continuous(range = c(2, 9), name = "out-strength") +
coord_equal() +
theme_void()
For a longitudinal figure, compute one layout and reuse it across waves, or use static_actor_positions = TRUE in the high-level plot method. Allowing a new layout in every facet can make stable actors look as though they moved. The manual_plotting vignette in the netify documentation develops this seam layer further.
NoteThe full toolbox: styles, themes, ego networks, layers (self-study)
Styles. Every plot above used the defaults. netify ships preset styles, each supplied as a function that you pass to style =, as well as named palettes:
head(list_network_styles(), 4)| style | description | background |
|---|---|---|
| style_rose | Soft rose, purple, and gold | light |
| style_red_blue | Red and blue accents on a clean background | white |
| style_orange_teal | Orange and teal on gray | gray |
| style_solarized | Solarized color scheme for reduced eye strain | dark |
head(list_palettes(), 4)#> [1] "rose" "red_blue" "cyberpunk" "scientific_blue"plot(class_net, style = style_tufte)
plot_with_style() and the theme_netify() / theme_publication_netify() family cover the rest of the aesthetics, and animate_netify() renders a longitudinal object as an animation. It needs gganimate and additional render time, so it is shown rather than run:
animate_netify(net_long_conf, static_actor_positions = TRUE)Ego networks. ego_netify() carves out one actor’s neighborhood, and the plot method has layouts built for exactly that shape:
usa_ego <- ego_netify(net14, ego = "United States")
plot(usa_ego, layout = "radial")
An ego network contains the focal actor, its directly connected alters, and the alter-to-alter ties observed inside that boundary. Every statistic is therefore local to the boundary. A high betweenness score here means central inside this extracted neighborhood, not globally influential in the parent network.
The boundary is part of the analysis, not a plotting preference. In a weighted network it changes when the tie threshold or neighborhood direction changes:
w14 <- as.matrix(net14)
strong_cut <- quantile(w14[w14 > 0], 0.90, na.rm = TRUE)
usa_out <- ego_netify(net14, ego = "United States",
ngbd_direction = "out")
usa_out_strong <- ego_netify(net14, ego = "United States",
ngbd_direction = "out",
threshold = strong_cut)
c(any_positive_outgoing = n_actors(usa_out),
strongest_outgoing = n_actors(usa_out_strong))#> any_positive_outgoing strongest_outgoing
#> 17 10Report the direction, threshold, time window, and parent network whenever you compare egos. remove_ego_edges() is useful when you want a picture that reveals alter-to-alter ties, but removing those edges changes the graph and should not silently change the object used for summaries. One last distinction matters: carving an ego network out of a complete network is not the same data design as collecting an egocentric survey. The latter may not observe all alter-to-alter ties, so full-network measures do not become valid just because the result can be plotted as a network.
Multilayer. Two relations on the same actors can live in one object, with cooperation and conflict side by side and n_layers() now returning 2:
net14_conf <- netify(ic14, actor1 = "i", actor2 = "j",
symmetric = FALSE, weight = "matlConf")
ml <- layer_netify(list(net14, net14_conf),
layer_labels = c("coop", "conflict"))
mlThe ego_networks and multilayer_networks vignettes take both to full depth; no later day fits models to these shapes, which is why they live in this self-study section.
10 Handing the object to a model
The benefit of doing this through one object rather than ad hoc reshape calls is portability. Every modeling package in this course wants the network in its own shape: statnet and ergm want a network object, igraph wants an igraph, and the latent-variable models in lame want a stack of matrices. You built the network once with the covariates aligned. When you export it, that alignment comes along.
One prerequisite move first. Two of the destinations, ergm on Day 11 and the binary family of models on Day 10, want a 0/1 network, and ours carries counts. binarize() makes the weighted-to-binary decision explicit instead of letting it hide in an export default:
w <- as.matrix(net14)
net14_bin <- binarize(net14, threshold = median(w[w > 0], na.rm = TRUE))
c(density_weighted = summary(net14)$density,
density_binary = summary(net14_bin)$density)#> density_weighted density_binary
#> 0.9869281 0.4934641Thresholding at the median positive weight turns “essentially everyone signals everyone” (density 0.99) into a network of above-typical cooperation (density 0.49) that a binary model can actually learn from. The threshold is your claim: zero, the median, a substantive cutoff. Saying it out loud is the point.
Now the doors. Each export is one line, and the covariates ride along in whatever form the destination expects:
net_stat <- to_statnet(net_cov)
net_stat#> Network attributes:
#> vertices = 18
#> directed = TRUE
#> hyper = FALSE
#> loops = FALSE
#> multiple = FALSE
#> bipartite = FALSE
#> verbCoop: 18x18 matrix
#> matlConf: 18x18 matrix
#> total edges= 302
#> missing edges= 0
#> non-missing edges= 302
#>
#> Vertex attribute names:
#> i_log_gdp i_polity2 vertex.names
#>
#> Edge attribute names:
#> matlConf_e verbCoopnetwork::list.vertex.attributes(net_stat) # The nodal covariates rode along#> [1] "i_log_gdp" "i_polity2" "na" "vertex.names"net_ig <- to_igraph(net_cov)
# ... and igraph's entire toolkit is immediately usable. The three most
# eigen-central states, straight off the exported object:
tail(sort(round(igraph::eigen_centrality(net_ig)$vector, 2)), 3)#> Russian Federation China United States
#> 0.67 0.74 1.00to_statnet() attaches nodal covariates as vertex attributes and dyadic ones as network-level matrices ready for ergm’s edgecov(); to_igraph() hands you every centrality and community routine igraph owns, no re-wrangling. The third door is to_lame(), and Day 10 lives behind it, so today we only read the sign on it: the export splits the covariates into the sender / receiver / dyad arrays ame() needs and even writes the model call you will run:
to_lame(net_cov)$ame_call#> [1] "lame::ame(Y = nl$Y, Xdyad = nl$Xdyad, Xrow = nl$Xrow, Xcol = nl$Xcol, family = \"normal\", nscan = 1000, burn = 500)"And nothing is trapped in the object. unnetify() closes the loop back to the dyadic data.frame we started from, covariates and all: the shape you would want for a merge, a regression table, or a colleague who has never heard of any of this:
head(unnetify(net_cov), 3)| from | to | time | verbCoop | matlConf | from_id | to_id | i_polity2_from | i_polity2_to | i_log_gdp_from | i_log_gdp_to |
|---|---|---|---|---|---|---|---|---|---|---|
| Afghanistan | Brazil | 1 | 0 | 0 | Afghanistan_1 | Brazil_1 | -1 | 8 | 23.7045 | 28.25614 |
| Afghanistan | China | 1 | 749 | 18 | Afghanistan_1 | China_1 | -1 | -7 | 23.7045 | 29.96645 |
| Afghanistan | Egypt | 1 | 4 | 0 | Afghanistan_1 | Egypt_1 | -1 | -4 | 23.7045 | 26.47765 |
The network object is the single source of truth, and the modeling packages are just different views of it. This is the entire reason netify exists, and it is why Day 8 is about the object: get it right once, and days 9 through 14 are downstream of it.
NoteThe full toolbox: three lines for the tidyverse (self-study)
The broom trio works on netify objects too (this needs dplyr installed): as_tibble() and tidy() give the tidy edgelist, glance() the one-row network summary: handy inside a purrr pipeline over many networks.
dplyr::as_tibble(net14)[1:3, ]| from | to | verbCoop |
|---|---|---|
| Afghanistan | Brazil | 0 |
| Afghanistan | China | 749 |
| Afghanistan | Egypt | 4 |
generics::glance(net14)[, 1:6]| net | num_actors | density | num_edges | prop_edges_missing | mean_edge_weight |
|---|---|---|---|---|---|
| 1 | 18 | 0.9869281 | 302 | 0 | 551.0298 |
Break. The object is built. After the break we ask the only question that matters: now that you have a network, what do you do with it? That question is the rest of the course.
11 Now what do you DO with it?
You have spent the descriptive week that precedes this course, and most of today, describing networks: degree, closeness, betweenness, communities, and just now summary_actor() on the trade data we are about to load. Description is useful and it is where everyone starts. But it has a ceiling, and the rest of this course lives above that ceiling. The next example shows where that ceiling is, followed by a map of the models that take us beyond it.
We switch to a cross-sectional directed trade network of thirty countries because it gives us a clean export example with meaningful source-target direction and dyadic covariates. Day 9 uses different applications, but the object-shaping work is the same.
# The .rda ships in this day's data/ folder alongside several day-9
# objects; load it into its own environment and take only the one
# data.frame we need, so nothing else lands in your workspace unannounced
e <- new.env()
load(file.path(day8_dir, "data", "trade_netlet.rda"), envir = e)
trade <- e$trade
trade_net <- netify(
trade, actor1 = "Var1", actor2 = "Var2",
symmetric = FALSE, weight = "trade",
nodal_vars = c("pop1", "gdp1", "polity1"),
dyad_vars = c("conflicts", "distance", "shared_igos"),
dyad_vars_symmetric = c(FALSE, TRUE, TRUE)
)
sa_trade <- summary_actor(trade_net)
head(sa_trade[order(-sa_trade$strength_sum_out),
c("actor", "strength_sum_out", "strength_sum_in")])| actor | strength_sum_out | strength_sum_in | |
|---|---|---|---|
| 30 | USA | 84.38194 | 87.03976 |
| 7 | CHN | 55.54867 | 52.03601 |
| 15 | JPN | 50.45011 | 44.77305 |
| 26 | TAW | 45.96767 | 45.42541 |
| 2 | AUL | 39.70973 | 38.99810 |
| 29 | UKG | 38.33812 | 39.78094 |
11.1 The ceiling of description
Description tells you who trades the most, who sits at the centre, and who clusters with whom. It is necessary, and you should always do it. But the moment you want to say something causal or explanatory, such as “Do democracies trade more with each other?” or “Does sharing an alliance raise trade?”, three questions arrive that description simply cannot answer:
- Why? Is the pattern homophily, reciprocity, transitivity, or a covariate? A centrality ranking cannot tell you.
- How sure? Description hands you a number with no uncertainty around it. No interval, no test.
- What next? Description cannot predict a tie it has not seen.
netify itself hands you the first rung above that ceiling, and it is worth climbing it right now so you can feel exactly where it ends. The question we just posed, whether democracies trade more with one another, is a homophily question, and homophily() answers its descriptive half with a simple dyad-permutation benchmark. It holds the observed ties fixed, computes attribute similarity for every eligible dyad, shuffles those similarity values across dyads, and asks whether connected pairs look more similar than that shuffled reference would suggest. This treats dyads as exchangeable, so it is an exploratory benchmark, not a node-label permutation, a correction for network dependence, or a causal test.
set.seed(6886) # homophily()'s permutations draw from the RNG
h <- homophily(trade_net, attribute = "polity1")
h[, c("homophily_correlation", "similarity_difference", "p_value")]| homophily_correlation | similarity_difference | p_value |
|---|---|---|
| 0.0377607 | 0.5311257 | 0.2527473 |
The observed regime-similarity correlation is a modest 0.038, with a permutation p-value of 0.25: under this crude dyad-exchangeability reference, we do not detect regime homophily. That is a descriptive result under one simple reference distribution, not a covariate-adjusted answer.
The same permutation logic also answers “is this network’s structure surprising at all?”: compare_to_null() simulates networks from a null model at the observed density and asks where your statistic falls among them:
null_recip <- compare_to_null(trade_net,
fn = function(net) c(reciprocity = summary(net)$reciprocity),
n_sim = 200, model = "erdos_renyi", seed = 6886, verbose = FALSE)
null_recip[, c("metric", "observed", "null_mean", "p_value", "extreme")]| metric | observed | null_mean | p_value | extreme |
|---|---|---|---|---|
| reciprocity | 0.9640582 | -0.0047148 | 0.0099502 | TRUE |
Observed reciprocity of 0.96 against a null distribution centred at 0: p = 0.01, about as extreme as two hundred draws can call anything. Trade flows are highly reciprocal under this reference distribution. Day 9’s social relations model gives reciprocity a parameter on a different application, while Days 11 and 13 reuse the broader logic of comparing what we observed with an explicit reference world.
But look hard at what both of these tests have in common: one statistic at a time, nothing held constant, no prediction. A permutation p on homophily cannot tell you whether democracies trade because they are democracies or because they are rich and near each other: it has no way to hold anything else fixed. That is the ceiling. The instant you want covariates adjusted for one another, you reach for a regression. And that is where the one problem of the whole course walks in.
NoteThe full toolbox: attribute analysis and null-model primitives (self-study)
Attribute analysis. homophily() has siblings. mixing_matrix() is the categorical version: who ties with whom, across groups (back to the classroom, where gender is categorical):
mm <- mixing_matrix(class_net, attribute = "gender")
mm$summary_stats[, c("attribute", "assortativity", "diagonal_proportion")]| attribute | assortativity | diagonal_proportion |
|---|---|---|
| gender | 0.1215805 | 0.6666667 |
plot_mixing_matrix(mm)
dyad_correlation() screens dyadic covariates against the outcome network, previewing Day 10’s Xdyad logic. plot_homophily() draws the distribution behind the permutation test you ran above:
dyad_correlation(trade_net, dyad_vars = "distance")[
, c("dyad_var", "correlation", "p_value")]| dyad_var | correlation | p_value |
|---|---|---|
| distance | -0.1645388 | 1.1e-06 |
plot_homophily(h, netlet = trade_net, type = "distribution",
attribute = "polity1")
attribute_report() runs the whole battery, including homophily, mixing, dyadic and centrality correlations, and attribute summaries, in one call. Its output is long, so here is just its one-line verdict; remember this function when Homework 3 asks for attribute EDA:
ar <- attribute_report(trade_net, node_vars = "polity1",
dyad_vars = "distance")
ar$overall_summary$summary_points#> [1] "No strong attribute-network relationships detected with current thresholds"How to read this battery without over-reading it. These functions are useful screens, but none of them is a jointly adjusted network model:
- For a continuous attribute,
homophily(method = "correlation")uses the negative absolute difference between the two actors as its similarity score; the reported homophily correlation relates that score to tie presence. For a categorical attribute, similarity is simply same category versus different category. - A mixing matrix shows where observed ties fall across groups. Raw counts are partly driven by how many dyads each group pairing makes possible, while row normalization answers the directional question “where do ties from this group go?” A large diagonal is not automatically a mechanism or a causal effect.
dyad_correlation()is a bivariate screen. Its usual p-value comes from an ordinary correlation test and does not adjust for dyads sharing actors. It is especially misleading if the proposed predictor is the outcome itself, or a component used to construct the outcome.attribute_report()organizes these diagnostics in one place. It does not make them conditional on one another, and it does not turn centrality into influence.
For longitudinal or multilayer objects, the functions repeat the descriptive calculation by wave or layer. That is useful for finding patterns to investigate, but it is not a model of change over time or dependence across layers.
Null-model primitives. compare_to_null() is a convenience wrapper, and the parts are yours to use directly. simulate() draws whole networks from a null family: "erdos_renyi" (density only), "configuration" (preserves the degree sequence), or "dyad_permutation" (relabels actors and preserves everything else). bootstrap_netlet() puts an interval on any statistic of the observed network:
sims <- simulate(net14, nsim = 3, model = "configuration", seed = 6886)
summary(sims[[1]])$density # Same degrees as net14, shuffled structure#> [1] 0.9869281bootstrap_netlet(trade_net,
fn = function(net) c(reciprocity = summary(net)$reciprocity),
n_boot = 200, seed = 6886, verbose = FALSE)| metric | point | n_valid | mean | sd | lower | upper |
|---|---|---|---|---|---|---|
| reciprocity | 0.9640582 | 200 | 0.9594435 | 0.0151486 | 0.9250803 | 0.9801347 |
Choosing which null to simulate from is itself the modeling decision. “Surprising given density” and “surprising given the degree distribution” are different questions, and Day 11 makes that choice into a career.
11.2 The one problem
The moment you want to explain a tie with a covariate, you write a regression with one row per dyad. And a regression assumes the rows are independent. In a network they are not, because the same actor sits in a whole band of rows at once. Let me make that concrete before I prove it hurts.
n <- 30 # Actors in the trade network
dyad_rows <- n * (n - 1) # Rows in a directed dyadic regression
rows_per_actor <- 2 * (n - 1) # Rows that contain any one given actor
c(dyad_rows = dyad_rows,
rows_per_actor = rows_per_actor,
share = round(rows_per_actor / dyad_rows, 3))#> dyad_rows rows_per_actor share
#> 870.000 58.000 0.067870 rows look like 870 independent facts. But nudge one country, perhaps because Japan trades more than its size predicts, and 58 of those rows, about 7% of the entire regression, move together, because Japan is the sender or the receiver in every one of them. The rows are not independent draws. They are coupled in bands, one band per actor.
Does that coupling actually break anything? Let’s not argue about it: let’s simulate a world where we know the truth and watch the regression lie.
set.seed(6886)
n <- 25
# One experiment: a predictor x with NO true effect on y.
# The common actor-level case, x is an ACTOR attribute (gdp, regime,
# ideal point) copied across every dyad that actor is in.
sim_once <- function(coupled_x = TRUE) {
a <- rnorm(n) # Actor "sociability" -> structures the ERRORS
z <- rnorm(n) # An actor attribute with ZERO true effect
d <- expand.grid(i = 1:n, j = 1:n)
d <- d[d$i != d$j, ]
d$x <- if (coupled_x) z[d$i] else rnorm(nrow(d)) # z[d$i]: actor i's single value STAMPED onto every dyad i sends, this is the coupling
d$y <- a[d$i] + a[d$j] + rnorm(nrow(d)) # y depends on WHO is in it
coef(summary(lm(y ~ x, data = d)))["x", "Pr(>|t|)"]
}
p_actor <- replicate(500, sim_once(coupled_x = TRUE))
p_dyad <- replicate(500, sim_once(coupled_x = FALSE))
c(actor_attribute = mean(p_actor < 0.05),
dyad_specific = mean(p_dyad < 0.05))#> actor_attribute dyad_specific
#> 0.502 0.058Read that carefully, because it is the single most important number today. We built a world where x has no effect on y whatsoever. A well-calibrated test should call it “significant” about 5% of the time. When the predictor is an actor attribute, which GDP, population, democracy, and ideal point all are, naive OLS calls it significant 50% of the time. That is a coin flip wearing a lab coat. When the predictor varies dyad by dyad, the false positive rate is a well-behaved 6%. The damage comes from actor-level coupling, the same coupling the arithmetic above described. One line of the simulation creates it: d$x <- z[d$i] takes an actor-length vector z and indexes it by the dyad’s sender column, stamping each actor’s single value across every dyad it appears in. That is the “copied into every dyad that country is in” move you made by hand back in Section 8, now shown to be the source of the false result.
ImportantThe Checkpoint to Carry Forward
Before you model anything of your own, ask: Is your predictor an actor attribute or a dyadic attribute? If it is an actor attribute, as most of the interesting ones are, your effective sample size is closer to your number of actors than your number of dyads (for the 30-country trade network you built today, that is about 30, not its 870 directed dyads or the 25-actor toy the simulation above used), and a naive p-value is not to be trusted. This is re-runnable and it stands alone; it is the one thing to remember if you remember nothing else from today.
NoteA first taste of one fix, previewing Day 13
You do not have to change the model to make the p-value stop lying: one option is to change only how the standard errors are computed. Cluster them on the sender:
suppressMessages({ library(sandwich); library(lmtest) })
set.seed(6886); n <- 25
fp_clustered <- mean(replicate(500, {
a <- rnorm(n); z <- rnorm(n)
d <- expand.grid(i = 1:n, j = 1:n); d <- d[d$i != d$j, ]
d$x <- z[d$i]; d$y <- a[d$i] + a[d$j] + rnorm(nrow(d))
m <- lm(y ~ x, data = d)
coeftest(m, vcov = vcovCL(m, cluster = ~ i))["x", "Pr(>|t|)"]
}) < 0.05)
fp_clustered#> [1] 0.084Clustering on the sender pulls the false-positive rate from ~50% down to about 8%: without changing a single coefficient, just by telling the standard errors that rows sharing a sender are not independent. It still overshoots 5% because with only ~25 actors there are too few clusters for the robust standard errors to be exactly calibrated. Day 13 treats that small-network limitation directly and develops the full dyadic cluster-robust calculation.
11.3 The map
The dependence will not go away because the dependence is the network: it is the thing you are trying to study, not noise to be scrubbed out. You do not get to ignore it. You get to choose where to put it. There are only a handful of places, and each one is a method, and each method is a day.
The cleanest way to see it is one equation. Write \(\eta_{ij}\) for the conditional linear predictor of a tie, so \(\mathbb E(Y_{ij}\mid\text{covariates, latent terms})=g^{-1}(\eta_{ij})\). Every day for the rest of the week adds, replaces, or refuses one piece:
\[ \eta_{ij} \;=\; \underbrace{\beta^{\top} x_{ij}}_{\text{observed covariates}} \;+\; \underbrace{\,a_i + b_j\,}_{\substack{\text{additive}\\ \text{SRM · Day 9}}} \;+\; \underbrace{\;\alpha_{z_i z_j}\;}_{\substack{\text{discrete roles}\\ \text{blockmodel · Day 9}}} \;+\; \underbrace{-\lVert z_i - z_j\rVert}_{\substack{\text{shared position}\\ \text{distance · Day 10}}} \;+\; \underbrace{\,u_i^{\top}\Lambda\, v_j\,}_{\substack{\text{sender/receiver factors}\\ \text{factor/AME · Day 10}}} \]
This is a map, not a recommendation to put every latent term in one kitchen-sink model. The block, distance, and factor terms are usually alternative ways to represent the same residual structure. For binary ties, \(g^{-1}\) is commonly the inverse logit or Probit; for Gaussian outcomes, the residual error is added on the outcome scale.
The two geometry terms are not the same move, and Day 10 is built around exactly that contrast. The distance idea gives every actor a single latent position \(z_i\) and a tie forms when two actors sit close: a symmetric distance, \(-\lVert z_i - z_j\rVert\), so \(i\to j\) and \(j\to i\) get the same pull. The factor model instead gives each actor a separate sender vector \(u_i\) and receiver vector \(v_j\), so the bilinear term \(u_i^{\top}\Lambda v_j\) can carry directed, asymmetric structure that one shared position cannot. Day 10 opens with the distance intuition, the cleanest picture of latent geometry there is, and then shows, on data where we know the truth, exactly where it breaks and why the factor model contains it and goes further.
Those four added terms are one philosophy: represent dependence with latent actor structure in the conditional mean or link. Conditional on those latent effects, distinct dyads may be independent; after we average over them, dyads that share actors are dependent. This can reduce distortion in the covariate estimates, provided the latent effects are not themselves correlated with the covariates. When an actor trait like GDP is tangled up with how much that actor sends, the random-effect route can bias the very coefficient you care about, and that is where the robust-inference camp pushes back. Three other methods take a different route:
- ERGM (Day 11) does not decompose the variance. It writes the probability of the whole graph directly, \(p(Y\mid\theta)\propto\exp\{\theta^{\top} g(Y)\}\), where \(g(Y)\) counts edges, triangles, shared partners. Dependence goes into the model of the structure itself. Day 11 ends by carrying that same machinery across panel waves: the temporal ERGM.
- SAOM (Day 12) treats the graph as the moving part and models \(P(Y^{(t+1)}\mid Y^{(t)})\) as actors add and drop ties. Dependence goes across time.
- Dyadic cluster-robust inference (Day 13) keeps the naive \(y_{ij}=\beta^{\top}x_{ij}+\varepsilon_{ij}\), with no extra term at all, and changes only how \(\mathrm{Var}(\hat\beta)\) is computed. Dependence goes into the standard errors (you just saw the four-line preview). Day 13 builds this out into a full session on causal inference with network data: what randomization buys you, what spillovers do to your estimand, and which claims an observational design can support.
And then Day 14 refuses to pick a cell at all. It puts the models head to head on one dataset, offers a decision framework for which question each one actually answers, and tours what is out there beyond this course (relational event models, multilayer networks).
The whole map fits in one table. One column, the measurement role, gets explained immediately afterward:
| where the dependence goes | method | the one sentence | the measurement role | day |
|---|---|---|---|---|
| describe it, don’t model | centrality, communities | who is central, who clusters, but not why, no uncertainty, no prediction | not applicable | Days 1-8 |
| into LATENT EFFECTS (additive) | social relations model | some actors just send/receive more; reciprocity, soak it into random effects | who over-sends, who over-receives | 9 |
| into LATENT EFFECTS (discrete) | stochastic blockmodel | actors fall into blocks; ties depend only on block membership | recovered roles / blocs | 9 |
| into LATENT EFFECTS (distance) | latent distance model | actors near each other in a latent space tie, homophily as distance | the map itself: latent positions | 10 |
| into LATENT EFFECTS (factor) | latent factor / AME | actors with similar latent roles tie, captures triangles description missed | latent relations: who leans toward whom | 10 |
| into the MODEL (structure) | ERGM | a coefficient on the triangle asks whether there is closure beyond nodal traits | not applicable | 11 |
| across TIME | SAOM (RSiena) | model how ties are added and dropped wave to wave | not applicable | 12 |
| into the STANDARD ERRORS | dyadic cluster-robust SEs + causal design | keep the naive model, fix the inference, and ask what causal claim survives | not applicable | 13 |
| all of the above, head-to-head | which model when? + the frontier | one dataset, every model, and a framework for choosing among them | not applicable | 14 |
Every later day fills in one or two rows of this table and points back at it. Day 9 takes the additive and discrete latent-effect rows, Day 10 takes the two continuous latent-effect rows, Days 11 through 13 take one row each, and Day 14 runs the whole table on a single dataset. This sequence continues the foundation you built in Days 1 through 7 rather than starting a separate course.
11.4 The latent terms: nuisance, or the thing itself
One column needs more explanation: the measurement role. Everything above treats the latent terms as a nuisance: the \(a_i\)’s, the blocks, the \(z_i\)’s and \(u_i\)’s give the modeled dependence somewhere to live while we estimate the covariate associations. Days 9 and 10 often use them that way. There is also a second use that I find especially interesting: sometimes the latent structure is the estimand. You did not come for \(\beta\) at all; you came for the map.
A useful analogy comes from test theory. Psychometricians call it item response theory, or IRT, and the acronym will come back in a moment. No single exam question measures ability, but forty noisy questions together help locate a student on an unobserved trait. Ideal-point models use the same logic: nobody directly observes ideology, yet hundreds of roll-call votes help place each legislator on a recovered dimension. Networks have a close analogue. Every tie is one noisy observation about a pair, and the positions the model recovers can summarize relational features we do not directly measure, such as affinity or alignment. Who has moved closer to whom in trade, votes, or conflict can be the finding itself, not merely a step toward a coefficient.
So why is this not just IRT with more rows? The inference is harder for three reasons:
- The map is not identified on its own. Rotate, reflect, or slide every position at once and the likelihood does not move. In Day 10’s factor variant, even the overall scale can be traded against the loadings. Only relative geometry means anything, so positions are interpretable only after Procrustes alignment, and the number of dimensions is chosen a priori, not estimated: Day 10 teaches both fixes.
- One graph is one realization. IRT gets many students crossed with many items; you get one network. \(n(n-1)\) cells, but coupled in the bands from Section 11.2, with no repeated items per pair: information about any one position accrues slowly, and “large-\(n\)” theory on a single graph remains nonstandard. Day 10’s uncertainty statements will therefore be more careful than you may be used to.
- The space absorbs whatever the mean model omits. Leave a relevant covariate out and it does not disappear: it can reappear as latent structure, so your “map of affinity” partly reflects the specification you chose. The same problem runs in the other direction. If the latent terms correlate with covariates you did include, the exogeneity argument for \(\beta\) fails too, with bias on the order of \(\mathrm{Cov}(x,\text{latent})/ \mathrm{Var}(x)\). Day 13 takes this up as a design problem.
The first reason fits in four lines, so let’s not take it on faith:
set.seed(6886)
z <- matrix(rnorm(10 * 2), 10, 2) # Ten actors in a 2-d latent space
th <- pi / 3 # Rotate the whole map by 60 degrees
Rm <- matrix(c(cos(th), -sin(th), sin(th), cos(th)), 2, 2)
rot_gap <- max(abs(dist(z) - dist(z %*% Rm)))
rot_gap # Distances are all the model sees#> [1] 4.440892e-164.4e-16: machine zero. Two maps that look completely different on paper are identical to the likelihood, and a reflection does the same trick exactly. Any substantive story you tell about “the x-axis” of a latent space is a story about an arbitrary rotation until you have pinned it down.
This map also has a rival cartography that often appears in econometrics. Econometricians hit the same dependence and made the opposite bet: put it in fixed effects, not random ones. Condition the actor-level heterogeneity out entirely (Graham 2017, Econometrica: the tetrad logit for network formation, which treats each actor’s degree heterogeneity as incidental parameters and conditions on a sufficient statistic for them), bias-correct the incidental-parameter contamination instead (Fernández-Val & Weidner 2018, Annual Review of Economics), or absorb it in two-way exporter/importer terms and fit the gravity model by PPML: Santos Silva and Tenreyro (2006) showed the logged version breaks under heteroskedasticity, and Jochmans (2017) showed that fixed-effect bias in these dyadic models is the same order as the standard error, so naive intervals under-cover. Then use dyadic-robust variance estimates (Cameron and Miller 2014; Tabord-Meehan 2019 supplies the published asymptotics), which is the estimator introduced on Day 13.
The sharpest point of contact is Bai (2009, Econometrica): his interactive fixed effects term \(\lambda_i^{\top} f_t\) is the same bilinear object as AME’s \(u_i^{\top}\Lambda v_j\), estimated as fixed parameters instead of given a prior. The fixed-effects approach never has to assume the heterogeneity is uncorrelated with your covariates, a practical advantage that we will take seriously, at the price that the latent relations are conditioned away rather than measured. So the two camps split cleanly: fixed effects plus robust standard errors is robust to correlated heterogeneity and silent about the relations themselves; random or latent structure is efficient and produces measures of the latent relations, but owes an exogeneity assumption in return. Neither dominates: the choice tracks exactly the question this section just posed: is the latent structure your nuisance, or your estimand? Day 14’s reading list opens the doorway into this literature (de Paula 2017 is the survey to start with).
Keep the one-line version: when the latent structure is your nuisance you mostly need it to stay uncorrelated with your covariates; when it is your estimand you additionally need it identified and its uncertainty quantified. Those are different demands, and Day 10 is where we start paying both bills.
11.5 The Model-Choice Tension
Two cells of that map, “dependence in the model” (ERGM) and “dependence in the variance” (latent space / AME), are not just different tools; they answer different questions, and the field debates how to choose between them. I wrote the AME paper, so I want to state the distinction directly rather than sell you my side: ERGM answers a question about mechanism that AME cannot; AME answers a question about covariate inference and prediction after modeling latent dependence. Pick the one that matches your question, not the one that “wins.” The full version, including ERGM’s old degeneracy reputation and what modern ergm 4.x does about it, is in the box. Day 14 stages the comparison on one dataset.
NoteERGM vs. Latent Space: The Substantive Disagreement
Two cells of that map, “dependence in the model” (ERGM) and “dependence in the variance” (latent space / AME), are not just different tools. They answer different questions, and the field debates how to choose between them. I wrote the AME paper, so I want to state the tradeoff clearly rather than sell you my side:
- ERGM makes dependence the substance. If your theory is about structure, such as “conflict spreads by closing triangles,” “friendship reciprocates,” “popularity begets popularity”: ERGM hands you a coefficient on exactly that mechanism and a hypothesis test for it. A latent-variable model cannot give you a “transitivity effect” to interpret; it absorbs the triangle-generating structure into geometry, and you can fit it beautifully while never learning the mechanism was closure. That is a genuine loss, and it is ERGM’s home turf.
- Latent space / AME makes dependence a nuisance to represent, so the covariate effects are estimated after accounting for the modeled latent structure and the model can predict held-out ties. It is rarely degenerate and scales comfortably to large graphs. But it will not name the structural mechanism for you.
The short version is: ERGM answers a question about mechanism that AME cannot; AME answers a question about covariate inference and prediction after modeling latent dependence. They optimize different things. Pick the one that matches your question, not the one that “wins.” Naive triangle-based ERGMs can be prone to degeneracy, while geometrically weighted terms such as gwesp and gwdsp and modern ergm estimation make many specifications much more stable, though they still require care. ERGMs also remain harder to fit at scale than latent-variable models. Days 10 and 11 let you work through that tradeoff on the same data, and Day 14 puts the two side by side.
11.6 Tomorrow: the social relations model, then blocks
The dependence is part of the network, not noise, and you have to decide where to put it. The cheapest place to start, and the one whose bookkeeping every later method reuses, is the variance. We say some actors just send more and some receive more, some pairs reciprocate, and we let random effects carry all of that so the covariates can return to their conditional role, as long as those random effects are not themselves entangled with the covariates, the caveat we just met. That is the social relations model, and it is the first half of tomorrow. The second half takes the same variance philosophy in a discrete direction: actors sorted into blocks, with ties depending only on which block they occupy. That is the stochastic blockmodel. One day, two ways of putting the dependence in latent actor effects: continuous actor effects before the break, discrete roles after it.
And here is the object we start from: the exact trade network you built two sections ago, handed to the model with one line:
lame_data <- to_lame(trade_net) # Day 10 uses the same matrix-list contract
sapply(lame_data[c("Y", "Xrow", "Xcol", "Xdyad")],
function(x) paste(dim(x), collapse = " x "))#> Y Xrow Xcol Xdyad
#> "30 x 30" "30 x 3" "30 x 3" "30 x 30 x 3"Tomorrow’s very first question is the one this half of today just made unavoidable: What exactly goes wrong if I just run the regression? By the end of Day 9 one coefficient will have flipped sign, with the standard errors left untouched. That sign flip is the p-value lie you saw today, caught and corrected by putting the dependence where it belongs.
11.7 After the Code Runs: Say What Network You Built
Today does not end with a model coefficient. It ends with a defensible measurement statement. Before interpreting density, centrality, reciprocity, or any other network summary, write down four things in ordinary political science language:
- Who can appear? Name the actors and the eligibility rule. A state-year network, for example, should not quietly treat states that do not yet exist as eligible isolates.
- What makes a tie? State the event, threshold, direction, weight, and time window. “Conflict” is too vague if the code actually records at least one material-conflict event initiated during a calendar year.
- What do zero and missing mean? A zero is an observed eligible non-tie. Missing means the tie was not observed. A structural zero means the tie was impossible by design.
- What descriptive pattern deserves explanation? Name the politically interesting concentration, asymmetry, clustering, or change over time, but do not turn that pattern into a mechanism before fitting a model that represents the relevant alternatives.
TipA human version of the result
“We construct an annual directed network of material conflict among states. A tie records at least one material-conflict event initiated by one state against another during the year. State pairs that do not coexist are excluded rather than coded as peaceful dyads. The resulting network is sparse, and a small number of states initiate a large share of the recorded conflict. That concentration motivates separate model terms for how often each state initiates conflict and how often each state is targeted.”
That paragraph does useful substantive work. It identifies the political object, tells the reader what an absence means, reports the main descriptive pattern, and explains why the next model needs actor-level heterogeneity. It does not claim that a high-centrality state is powerful, that reciprocity proves retaliation, or that a descriptive association is causal. Those interpretations require additional theory, comparison, and evidence.
12 Your turn
13 Reading
- Minhas, Dorff, Gallop, Foster, Liu, Tellez & Ward (2022), “Taking Dyads Seriously,” Political Science Research and Methods,
minhas_etal_2022.pdf. The paper behind the whole spine of this course. Sections 1 and 2 are exactly the “why can’t I just run the regression?” argument you saw in part two; the rest is where we are going. - Ward, Stovel & Sacks (2011), “Network Analysis and Political Science,” Annual Review of Political Science:
ward_etal_2011.pdf. A gentle, discipline-spanning tour of why networks need their own methods. Read this one first if today was your first serious contact with network modeling.
NoteAppendix: what we did not show, and where it lives
Today’s main path plus the toolbox asides cover the netify surface a working modeler needs. The deliberate omissions, named so you know they are omissions and not oversights:
- Manual plot assembly:
net_plot_data(),assemble_netify_plot(), thenetify_node()/netify_edge()/netify_label()geom layers, andget_node_layout()/get_edge_layout()for freezing coordinates across panels; themanual_plottingvignette in the netify docs is the walkthrough. - Package internals:
new_netify()and theget_adjacency*()family are the low-level constructors thatnetify()wraps; theinternalsvignette if you are curious. to_lame()/to_amen()in anger: Day 10 exercises the whole pipeline; today you only read$ame_calloff the export.to_dbn()andfrom_lame_fit(): the dynamic-bilinear-network and fitted-model bridges; days 10 and 14 own those.netify_measurements(),get_actor_time_info(),melt(), and the long-formnetify_to_*aliases: superseded for our purposes byn_actors()/n_periods()/n_layers(),subset_netify(), and the shortto_*names you learned today. Same story foradd_vertex_attributes()/add_edge_attributes()(straight aliases of theadd_node_vars()/add_dyad_vars()you already know) and thedecompose_igraph()/decompose_statnet()variants ofdecompose_netify()for foreign objects.- Modeling ego or multilayer objects: the objects themselves appeared in the plotting toolbox aside; the
ego_networksandmultilayer_networksvignettes go deeper, and no later day fits models to them.
Two habits to leave with. First: whenever you set symmetric, weight, sum_dyads, or missing_to_zero, say out loud what claim about the world you are making: the defaults are conveniences, not truths. Second: open ?netify_workflows once tonight. It is the cheat sheet for this entire course: every pipeline recipe (ERGM, lame, temporal, large networks) in one help page, and it will still be there in August when you cannot remember which day we did what.
Session info
sessionInfo()#> R version 4.3.3 (2024-02-29)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 24.04.3 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
#>
#> locale:
#> [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
#> [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
#> [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
#> [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] lmtest_0.9-40 zoo_1.8-15 sandwich_3.1-1 ggplot2_4.0.3 igraph_2.2.2
#> [6] netify_1.5.3
#>
#> loaded via a namespace (and not attached):
#> [1] tensorA_0.36.2.1 tidyr_1.3.2 generics_0.1.4
#> [4] lattice_0.22-5 digest_0.6.39 magrittr_2.0.5
#> [7] statnet.common_4.13.0 evaluate_1.0.5 grid_4.3.3
#> [10] RColorBrewer_1.1-3 fastmap_1.2.0 lame_1.3.5
#> [13] jsonlite_2.0.0 Matrix_1.6-5 ggrepel_0.9.6
#> [16] ggnewscale_0.5.2 network_1.20.0 backports_1.5.0
#> [19] purrr_1.2.2 viridisLite_0.4.3 scales_1.4.0
#> [22] abind_1.4-8 cli_3.6.6 rlang_1.2.0
#> [25] withr_3.0.2 yaml_2.3.12 otel_0.2.0
#> [28] parallel_4.3.3 tools_4.3.3 checkmate_2.3.4
#> [31] coda_0.19-4.1 dplyr_1.2.1 broom_1.0.11
#> [34] posterior_1.6.1 vctrs_0.7.3 R6_2.6.1
#> [37] matrixStats_1.5.0 lifecycle_1.0.5 htmlwidgets_1.6.4
#> [40] pkgconfig_2.0.3 pillar_1.11.1 gtable_0.3.6
#> [43] loo_2.9.0 glue_1.8.1 Rcpp_1.1.1-1.1
#> [46] xfun_0.55 tibble_3.3.1 tidyselect_1.2.1
#> [49] rstudioapi_0.17.1 knitr_1.51 farver_2.1.2
#> [52] patchwork_1.3.2 htmltools_0.5.9 rmarkdown_2.30
#> [55] labeling_0.4.3 compiler_4.3.3 S7_0.2.2
#> [58] distributional_0.6.0