Day 14: Which Model, When? The Frontier

Advanced Network Analysis · ICPSR

Author

Shahryar Minhas

Published

July 30, 2026

How to Use This Document

This QMD is the main walkthrough and analysis script. The download also includes three short helper files for setup, model comparison, and the walk-embedding demonstration. The walkthrough calls them where needed, and the rendered document already contains every result, so nobody needs to run a model during class. If you choose to rerun a chunk later, the included cache and the seed 6886 reproduce the displayed analysis.

The expensive fits are wrapped in cache_fit(), which stores them in cache/. The downloadable materials include those fitted objects, so the walkthrough loads them without refitting. Delete cache/ only when you deliberately want to estimate everything again from scratch.

Two packages appear today for the first time: dcr for the dyadic uncertainty correction and sir for the influence segment. lame returns from Days 9 and 10:

install.packages(c("netify", "blockmodels", "ergm", "network",
                   "ggplot2", "ggrepel", "ggraph",
                   "coda", "digest", "remotes",
                   "multiwayvcov"))   # dcr uses this dependency
remotes::install_github("s7minhas/lame")
remotes::install_github("jscarlson/dcr")    # Dyadic cluster-robust SEs
remotes::install_github("netify-dev/sir")   # Social influence regression

Where this day sits on the map. It doesn’t pick a cell; it runs the whole map. The syllabus promised, as its first objective, that you would “develop an ability to compare the relative merits of the various models covered for a given empirical application.” Today is where that promise is kept: one dataset through every compatible cross-sectional lens, a decision framework you can carry out of the room, and a practical tour of what exists beyond the course.

Today’s Path

Segment
The question that outlasts the course
One dataset: the rebel-cooperation network, third pass
Lens 0: the GLM, and shared-actor robust standard errors
Lens 1: the SRM: where does the actor dependence live?
Lens 2: the blockmodel: do actors come in kinds?
Lens 3: AME: a relational surface on top of the covariates
Lens 4: the ERGM: does a shared-partner term add information?
The model comparison table: one network, every cross-sectional lens
The decision framework
The frontier: new directions beyond the models covered in class
The frontier: influence: who moves whom, with a cached SIR example
The frontier: multilayer networks
The frontier: scale, and the adjacent field
The frontier: the econometrics of networks
Close: the map, filled in
Your turn
Reading

The course rundown, model comparison table, decision framework, and frontier sections on prediction and multiple relation layers form the main path. The first three organize the course, while the last two show how familiar ideas extend beyond the models fit in class.

1 The Question That Outlasts the Course

On Day 8, we began this part of the course with a claim: dyadic observations are not independent, and every method we teach makes a different decision about where to represent that dependence. Since then you have represented it with additive effects, discrete roles, continuous relational profiles, graph statistics, temporal models, uncertainty corrections, and causal designs. Each day added another part of the map.

Today there is no next cell. Today is the question you will actually face the week you get home:

I have network data and a question. What do I fit?

Notice the order. The data do not pick the model, and neither does your discipline’s habit. The question picks the model. A common mistake in applied network work is to answer a measurement question with an uncertainty correction, or a mechanism question with a descriptive tool, without noticing that the target changed.

So today has three movements:

One dataset, every compatible cross-sectional lens. We take the Syrian rebel cooperation network and fit the GLM, dyadic-robust inference, SRM, blockmodel, AME, and ERGM. The underlying binary network stays fixed, the covariates align where each family permits them, and the fitted network models face one shared goodness-of-fit yardstick. The temporal and causal tools return in the course map rather than being forced onto data that cannot answer their questions. The point is not to declare that one model won. It is to see the different jobs these tools perform on the same network.
The decision framework. A flowchart and a table keyed to the course map. This is the artifact to screenshot.
The frontier. A separate tour of methods beyond the course, each with a named entry point so “go read about X” becomes an actionable sentence.

1.1 What We Built Across Days 8 Through 13

The course has not been a parade of packages. Each day added a tool because a different feature of relational data became part of the question. Read these tables from left to right: the applied question comes first, the model target comes second, and the limitation keeps the result in its proper lane.

1.1.1 Days 8 Through 10: Build the Data, Then Represent Unobserved Structure

Day and tool	The question it answers	What it learns	What to remember
Day 8: `netify` and dyadic data construction	What exactly are the nodes, ties, missing dyads, directions, layers, and time points?	A coherent relational object and the descriptive patterns that motivate a model	Converting events to ties, combining relation types, binarizing counts, or turning missing dyads into zeros changes the outcome before estimation begins
Day 9: SRM with `lame::ame()` or `lame()`	Do some actors form or receive ties broadly more often than others, beyond measured predictors?	Covariate associations, additive actor effects, dyadic reciprocity, and their uncertainty	Additive effects absorb broad activity and popularity, but they do not by themselves represent third-order affinities
Day 9: `blockmodels` and NetMix	Can a small number of discrete relational roles summarize who connects to whom at one observation or across a panel?	Cross-sectional memberships or actor-year role mixtures, plus role-to-role tie probabilities	A block is a recurring partner-list type, not automatically a community, ideology, or institutional category; the Day 9 NetMix fit did not impose a transition model between years
Day 10: latent distance, latent factor models, and AME	Does a measured association survive broad actor heterogeneity and continuous residual relationship structure?	Latent-distance positions or a signed low-rank relational surface, together with coefficients, actor effects, and model-based uncertainty	Distance and factor models encode different geometry; their axes are not observed traits, and good fit alone does not establish causality

1.1.2 Days 11 Through 13: Put Structure in the Graph, Across Time, or in the Design

Day and tool	The question it answers	What it learns	What to remember
Day 11: ERGM	Which named graph configurations receive greater conditional weight after the other modeled features are held fixed?	Coefficients on change statistics through a graph likelihood, plus simulated-network checks	A positive closure term is a conditional structural association, not proof of a friend-of-a-friend behavioral mechanism
Day 11: TERGM	How do current ties depend on prior networks and current covariates across observed waves?	Temporal-network coefficients under a tie-oriented transition model	Stability combines persistence of ties and non-ties; the specification and time spacing define the temporal story
Day 12: SAOM	Which actor-level micro-changes receive more relative choice weight as a network and behavior co-evolve?	Rate and evaluation parameters by simulated method of moments	A coefficient changes the relative attractiveness of a candidate micro-step, not the unconditional probability that a tie exists
Day 13: Causal inference on networks	What comparison could identify an effect when one unit’s treatment may affect other units?	Design-based exposure contrasts and observational strategies when their assumptions are defended	Causal work begins with an intervention, exposure mapping, estimand, support, and a design or explicit identification strategy

The Software Handles, in One Place

netify() builds and checks the relational object and converts it to the structures expected by other packages. lame::ame() fits cross-sectional SRM and AME models; lame() handles longitudinal versions, while latentnet fits the latent-distance model introduced on Day 10. blockmodels fits cross-sectional discrete roles. NetMix estimates actor-year mixed memberships in longitudinal data; the Day 9 fit allowed those mixtures to differ by year but did not add a transition model between them. ergm fits cross-sectional graph models, and btergm is the main course implementation for TERGMs. RSiena fits SAOMs. Day 14 introduces dcr, which leaves a fitted coefficient unchanged and recalculates its covariance for repeated actors. The Day 13 design-based examples use known assignment probabilities, exposure mappings, and Horvitz-Thompson estimators; the observational negative-control example uses GMM to estimate a confounding bridge. These package names are handles for the methods, not substitutes for stating the question and estimand.

Four Jobs That Should Not Be Blurred Together

Data construction defines the outcome and who was at risk of a tie. A dependence model changes what the conditional mean or graph distribution can represent. An uncertainty correction such as DCR changes the covariance estimate without changing the fitted coefficient. A causal design or identification strategy says which comparison can be interpreted as an intervention effect. Better fit at one job does not quietly solve the other three.

The One Sentence of the Day

A model is a lens, not a truth: each one answers the question it was built for, so state your question first and pick the lens second.

1.2 How Each Estimator Learns From the Same Network

The rebel-cooperation data will not change across the cross-sectional comparison, but each procedure treats the network as evidence about a different pattern. The GLM looks for measured dyadic associations. The SRM looks for broad actor differences. The blockmodel looks for a small set of recurring relational types. AME looks for continuous pair-specific compatibility. The ERGM asks whether named graph configurations are unusually common. DCR is different because it does not change the fitted relationship model at all.

Lens	Statistical target	What moves during estimation
GLM	Maximize the sum of dyad-level Bernoulli log-likelihood contributions	Regression coefficients
SRM	In this capstone, sample the posterior for covariate coefficients, additive actor effects, and variance components, plus dyadic reciprocity when the outcome is directed	Conditional mean and uncertainty
Blockmodel	Maximize a variational lower bound on the observed-network likelihood	Actor membership probabilities and block-to-block tie probabilities
AME	In this capstone, sample the posterior for coefficients, additive effects, the low-rank surface, and variance components	Conditional mean, latent surface, and uncertainty
ERGM	Maximize the graph likelihood; MCMC approximates the expected graph statistics in its score	Coefficients on covariates and graph configurations
TERGM with `btergm`	Estimate a conditional tie-transition model by maximum pseudolikelihood and use a temporal bootstrap for uncertainty	Transition coefficients and bootstrap uncertainty
SAOM	Use simulated method of moments to make expected simulated change statistics match the observed targets	Rate and evaluation parameters

The TERGM and SAOM rows complete the course map, although they cannot be fitted to today’s single snapshot. DCR is not a model-fitting procedure because it has no new likelihood, posterior, moment condition, or coefficient target. The DCR section below introduces it as a separate covariance calculation around the GLM estimate.

That table explains why the coefficient on the same covariate can move. The models are not merely applying different standard errors to one fixed regression. Most of them rebuild the conditional mean around a different representation of network dependence.

The distinction between a fitting target and a search algorithm is worth keeping. The blockmodel’s target is an approximate likelihood bound, while variational EM is the update routine. The ERGM’s target is a graph likelihood, while MCMC provides the network expectations needed to optimize it. The capstone SRM and AME fits use MCMC to target posterior distributions. The fast longitudinal lame examples from Days 9 and 10 instead used ALS to obtain point estimates of a declared loss. The SAOM estimator repeatedly simulates the unobserved micro-step path and adjusts parameters until the expected simulated statistics match the observed change statistics closely enough.

For every lens, we will follow the same sequence:

Name the substantive question that motivates the lens.
State the estimand or conditional association.
Identify the unit of statistical information.
Walk through what the estimator updates and what it holds fixed.
Read the coefficient in the model’s own conditional language.
Check the feature that would make us abandon or expand the model.

1.2.1 What Each Lens Should Let Us Say

Lens	The applied question	The result we should be able to say out loud
GLM with DCR	Which measured associations remain after allowing uncertainty to reflect repeated organizations?	“Cooperation is more common for pairs in the same location, pairs closer in ideology, and pairs containing larger organizations. Their dyadic-robust intervals remain separated from zero.”
SRM	How much of the tie pattern reflects persistent differences in organizations’ overall propensity to cooperate?	“ANF and ASIM cooperate with unusually many organizations even after the measured predictors enter, and the actor effects account for much of the network’s hub structure.”
Blockmodel	Do organizations occupy a small number of relational types with distinct mixing patterns?	“The fit separates ANF and ASIM into a two-organization core and places the other 28 organizations together. That is a coarse hub-versus-periphery split, not evidence of two broad coalition communities.”
AME	Does a covariate association survive additive actor heterogeneity and latent pairwise compatibility?	“Shared location remains positive and power difference remains negative, while the ideological-distance interval crosses zero after latent pair-specific structure enters.”
ERGM	Is an explicit shared-partner configuration still needed after actor activity enters?	“In the diagnostic ERGM that includes a sociality term for every organization, the remaining GWESP association is too imprecise to separate from zero, and simulations cover the two selected scalar yardsticks. Because this coefficient conditions on a different specification, it should not be treated as directly comparable with the GWESP coefficient from the simpler model.”
Influence model	Whose past activity predicts whose later activity?	“After the focal pair’s own history and measured characteristics enter, states’ later conflict targets resemble the earlier targets of some allies and diplomatic partners.”

If we cannot complete the right-hand sentence without hiding behind the package vocabulary, we have not yet interpreted the model.

2 One Dataset: The Rebel-Cooperation Network, Third Pass

The data come from Gade, Gabbay, Hafez, and Kelly (2019) and cover cooperation among 30 armed organizations in the Syrian civil war from July 2012 through June 2015. A symmetric binary tie means that the pair claimed at least one tactical joint operation. The published AME analysis uses the square root of the number of claimed joint operations; this capstone simplifies the count to any versus none. The organization-level covariates are ideology score and group size. The pair-level covariates are ideological distance, power difference, shared location, and shared sponsor.

You have met this substantive setting before, but not always with the same outcome or specification. Day 9 used the Syrian organizations to introduce actor effects and blockmodels. Day 10 added latent factors. Day 14 adds a dyadic uncertainty correction and fits the cross-sectional lenses to the same 30-organization binary cooperation network so that they can be judged side by side. SAOM and TERGM remain on the course map, but this one-wave application cannot estimate change between waves.

load("data/gadeData.rda")          # A dyadic data.frame named `data`
data$coop <- as.numeric(data$coopActions > 0)

net <- netify(data, actor1 = "actor1", actor2 = "actor2",
              symmetric = TRUE, weight = "coop",
              nodal_vars = c("averageId_actor1", "size_actor1"),
              dyad_vars  = c("ideol_diff", "powerdiff", "loc", "spons"))
A <- as.matrix(net); diag(A) <- NA

Day 8 established the descriptive routine we use before fitting any model:

plot(net, edge_color = "lightgrey", add_points = FALSE,
     add_label = TRUE, label_color_var = "averageId_actor1",
     label_size_var = "size_actor1") +
  scale_color_gradient(low = "#A7A8AA", high = "#18453B") +
  labs(color = "avg ideology", size = "group size")

A labeled network of 30 Syrian armed organizations. ANF and ASIM are the most visibly connected organizations, and several organizations sit at the sparse edge of the graph. — The rebel cooperation network. Label size is group size, color is average ideology. The visual story is the modeling story: a few very active hubs, a quiet periphery, and ties bunching among ideological kin.

dens <- mean(A, na.rm = TRUE)
deg  <- rowSums(A, na.rm = TRUE)

The density is 0.195: 85 of the 435 possible pairs cooperate. The number of partners runs from 1 to 25, and that spread is the network’s most visible feature. ASIM cooperates with 25 other organizations, while 6 organizations have two or fewer partners.

The shared yardstick. To hold every model to the same standard, we build one common function around the actor-heterogeneity and transitivity-dependence checks used on Day 9. On this symmetric network it yields two scalar checks: actor heterogeneity, measured by the standard deviation of row means, and a centered, standardized third-order transitivity-dependence moment. The second quantity is not the ordinary triangle count, global clustering coefficient, or full shared-partner distribution.

calc_gof <- function(mat) {
  mat <- as.matrix(mat)
  actor_sd <- sd(rowMeans(mat, na.rm = TRUE), na.rm = TRUE)
  E <- mat - mean(mat, na.rm = TRUE)
  D <- 1 * (!is.na(E)); E[is.na(E)] <- 0
  s3 <- sd(c(mat), na.rm = TRUE)^3
  trans <- sum(diag(E %*% t(E) %*% E)) /
           (sum(diag(D %*% t(D) %*% D)) * s3)
  out <- c(actor_sd, trans)
  out[is.na(out)] <- 0
  setNames(out, c("actor_sd", "trans"))
}
gof_obs <- calc_gof(A)
round(gof_obs, 3)

#> actor_sd    trans 
#>    0.184    0.007

checkpoint(density = round(dens, 3),
           actor_sd = round(gof_obs["actor_sd"], 3),
           trans = round(gof_obs["trans"], 4))

#> ------------------------------------------------------------------
#> CHECKPOINT: density = 0.195   |   actor_sd = 0.184   |   trans = 0.007
#> ------------------------------------------------------------------

The centered transitivity score can be negative because it averages products of tie values after the overall density has been subtracted. A negative value means less third-order dependence than the density-based reference, not a negative number of triangles. For a symmetric network, the directed cycle and transitivity formulas are algebraically identical because \(E=E^\top\), so we report only trans rather than duplicate the same statistic under two labels.

Hold one number in your head all morning: actor_sd is 0.184. That is the hub-versus-periphery spread in the plot above, and it will be the statistic that separates the lenses. Notice too that trans is only 0.0067. Later, actor-heterogeneity specifications reproduce that statistic and the residual GWESP association becomes imprecise. That pattern is compatible with hub structure generating many observed triangles, but it does not establish the generating mechanism or prove that closure is absent.

3 Lens 0: Start With the Dyadic GLM

What question does it answer? “Conditional on my covariates, which pair-level traits go with cooperation?” It is not a network model. It is the familiar baseline that the network models will extend.

Because the outcome is undirected, the GLM must also be undirected. Keep one row for each unordered organization pair and enter an actor-level covariate as the sum of its values at the two endpoints. That matches the \(\beta(x_i+x_j)\) contribution used by the symmetric SRM, AME, and ERGM node-covariate terms.

pair_key <- ifelse(
  as.character(data$actor1) < as.character(data$actor2),
  paste(data$actor1, data$actor2, sep = "\r"),
  paste(data$actor2, data$actor1, sep = "\r")
)
glm_data <- data[!duplicated(pair_key), ]
rownames(glm_data) <- NULL
glm_data$averageId_sum <-
  glm_data$averageId_actor1 + glm_data$averageId_actor2
glm_data$size_sum <- glm_data$size_actor1 + glm_data$size_actor2

glm_fit <- glm(coop ~ ideol_diff + powerdiff + loc + spons +
                 averageId_sum + size_sum,
               data = glm_data, family = binomial)
round(summary(glm_fit)$coef, 3)

#>               Estimate Std. Error z value Pr(>|z|)
#> (Intercept)     -6.153      1.087  -5.659    0.000
#> ideol_diff      -0.253      0.128  -1.970    0.049
#> powerdiff       -0.119      0.033  -3.600    0.000
#> loc              3.235      1.026   3.154    0.002
#> spons           -0.300      0.424  -0.707    0.480
#> averageId_sum    0.277      0.081   3.418    0.001
#> size_sum         0.126      0.023   5.395    0.000

The coefficient estimates describe conditional associations in the usual GLM language. The ordinary standard errors, however, treat all 435 pair rows as independent. That is difficult to defend because each organization appears repeatedly. For example, the ANF-ASIM row and every other row containing ANF may share unmeasured ANF-specific shocks.

3.1 Why Shared Organizations Change the Uncertainty Calculation

Dyadic cluster-robust standard errors, usually shortened to DCRSEs or DCR, allow the regression errors for two dyads to be correlated whenever those dyads share either organization. DCR does not refit the regression and does not add actor effects. It keeps every GLM coefficient exactly as estimated above, then recalculates the covariance matrix used for standard errors, confidence intervals, and hypothesis tests.

In practical terms, DCR says that ANF-ASIM and ANF-Ahrar al-Sham should not count as two completely independent pieces of evidence. It collects the score contributions associated with each organization and uses their shared movement in a sandwich covariance calculation. Students do not need to reproduce that matrix algebra to understand the job it performs: same fitted associations, an uncertainty calculation that allows rows sharing an actor to move together.

Now keep the GLM point estimates and recompute only their uncertainty:

invisible(capture.output(          # dcr prints progress; keep it quiet
  dcr_out <- dcr::dcr(model = glm_fit, dyad_mem1 = "actor1",
                      dyad_mem2 = "actor2", data = glm_data)))
naive_se <- summary(glm_fit)$coef[, 2]
tab_dcr <- data.frame(
  estimate  = round(coef(glm_fit), 3),
  naive_se  = round(naive_se, 3),
  dcr_se    = round(dcr_out$dcrse, 3),
  inflation = round(dcr_out$dcrse / naive_se, 2))
tab_dcr

	estimate	naive_se	dcr_se	inflation
(Intercept)	-6.153	1.087	1.117	1.03
ideol_diff	-0.253	0.128	0.077	0.60
powerdiff	-0.119	0.033	0.084	2.54
loc	3.235	1.026	0.849	0.83
spons	-0.300	0.424	0.437	1.03
averageId_sum	0.277	0.081	0.122	1.51
size_sum	0.126	0.023	0.038	1.60

3.2 What DCR Does Not Fix

DCR is useful, but its job is narrow. It does not repair an omitted actor trait, a wrong functional form, or a missing network process. It still assumes that dyads with no shared organization are independent, so broad shocks or more distant network dependence remain outside the correction. Its large-sample justification concerns the number of actors, and 30 organizations is a small setting for that approximation. It also does not turn a conditional association into a causal effect.

A corrected standard error can be larger or smaller than the ordinary standard error because the estimated shared-actor covariances can have either sign. The point is not that DCR always makes inference more conservative. The point is that it uses an uncertainty calculation suited to rows that overlap through actors.

The one-liner behind the next paragraph

infl_med <- median(tab_dcr$inflation[-1])   # Drop the intercept

The median ratio of dyadic-robust to ordinary standard errors across the substantive terms is 1.27, and the change is not uniform: the power-difference standard error grows 2.5 times. Same location, ideological similarity, and combined group size remain clearly associated with cooperation, while the power-difference interval widens across zero. Carlson, Incerti, and Aronow (2024) explain the logic and limitations of this correction in dyadic research.

What this lens cannot do: everything else on today’s list. It has no account of why the rows are correlated; no actor effects, no groups, no geometry, no closure. And that is a checkable claim, not a slogan: simulate networks from the fitted GLM and score them against the yardstick, using one fitted probability for each unordered pair and then flipping independent coins:

set.seed(6886)
phat <- predict(glm_fit, type = "response")
P <- matrix(0, 30, 30, dimnames = dimnames(A))
ii <- match(glm_data$actor1, rownames(P))
jj <- match(glm_data$actor2, colnames(P))
P[cbind(ii, jj)] <- phat
P[cbind(jj, ii)] <- phat                   # Same undirected pair probability
gof_glm <- t(replicate(200, {
  S <- matrix(0, 30, 30)
  S[upper.tri(S)] <- rbinom(choose(30, 2), 1, P[upper.tri(P)])
  S <- S + t(S); diag(S) <- NA; calc_gof(S)
}))
rbind(observed = round(gof_obs, 3),
      round(apply(gof_glm, 2, quantile, c(0.025, 0.975)), 3))

#>          actor_sd  trans
#> observed    0.184  0.007
#> 2.5%        0.096 -0.021
#> 97.5%       0.140  0.064

The first clear miss is actor heterogeneity. The simulated networks put actor_sd in [0.096, 0.14] against an observed 0.184. Independent dyadic draws do not reproduce the observed hub structure, even with the covariates included. Hold that miss in mind because each model in the next set represents it differently. The GLM does reproduce this particular transitivity statistic, but that does not show that it reproduces every triangle-related feature or that its independence assumption is correct. This lens estimates conditional dyadic associations and adjusts their uncertainty for shared actors. It does not model the remaining network structure.

4 Lens 1: The SRM: Where Does the Actor Dependence Live?

What question does it answer? “How much of the network reflects persistent differences in organizations’ overall propensity to cooperate, and what happens to the covariate results once every organization receives an additive effect?” Because this network is undirected, there is one sociality effect per organization rather than separate sender and receiver effects.

The Day 9 additive model structure returns here: the same covariates, plus a random effect \(a_i\) for every group. The fitted object is included in the cache:

dl <- to_lame(net)
gade_srm <- cache_fit("gade_srm_lame135", ame(
  Y = dl$Y, Xdyad = dl$Xdyad, Xrow = dl$Xrow,
  family = "binary", R = 0, symmetric = TRUE,
  seed = 6886, burn = 10000, nscan = 10000, odens = 10,
  plot = FALSE, verbose = FALSE))

ess(gade_srm)     # Effective sample size per coefficient; the Day 9 habit

#>             intercept averageId_actor1_node      size_actor1_node 
#>                   250                   831                   371 
#>       ideol_diff_dyad        powerdiff_dyad              loc_dyad 
#>                   355                   366                   212 
#>            spons_dyad 
#>                   762

srm_ci <- t(apply(gade_srm$BETA, 2, quantile, c(0.025, 0.5, 0.975)))
round(srm_ci, 3)

#>                         2.5%    50%  97.5%
#> intercept             -5.694 -4.100 -2.753
#> averageId_actor1_node -0.009  0.175  0.361
#> size_actor1_node       0.064  0.115  0.181
#> ideol_diff_dyad       -0.527 -0.246 -0.041
#> powerdiff_dyad        -0.193 -0.110 -0.049
#> loc_dyad               0.889  1.777  2.976
#> spons_dyad            -1.141 -0.473  0.161

Read it the Day-9 way. Ideological distance has a clearly negative interval (\([-0.53, -0.04]\)), matching the negative association in the correctly symmetrized GLM, and shared location is strongly positive. The variance components say where the dependence lives:

vc <- apply(gade_srm$VC, 2, median)
round(vc, 3)

#>    va    ve 
#> 0.506 1.000

apm_top <- sort(gade_srm$APM, decreasing = TRUE)[1:3]
round(apm_top, 2)

#> ASIM  ANF  ATB 
#> 2.01 1.39 0.78

The actor variance va is 0.51 against a probit residual fixed at 1, confirming that organizations differ substantially in their overall tendency to cooperate. The additive effects also identify which organizations cooperate unusually widely after the measured predictors enter: ASIM and ANF have the two largest fitted effects.

G_srm <- gade_srm$GOF
gof_band <- function(G) apply(G[-1, , drop = FALSE], 2,
                              quantile, c(0.025, 0.975))
round(gof_band(G_srm)[, c("sd.rowmean", "trans.dep")], 3)

#>       sd.rowmean trans.dep
#> 2.5%       0.147    -0.094
#> 97.5%      0.207     0.066

Both selected scalar yardsticks land inside the SRM’s posterior-predictive bands. The \(a_i\) term is built to represent actor heterogeneity, so it closes the actor_sd miss posted by the GLM. The observed centered transitivity-dependence moment also falls inside the band, which means this one check does not detect residual third-order dependence beyond what the fitted covariates and actor effects reproduce. It does not prove that hubs generated the triangles or that closure is absent. Day 9 used square-root counts of claimed joint operations and ended by showing what an additive model still missed in that specification. Here the outcome is a binary indicator of any claimed joint operation and the fitted SRM covers both chosen checks. Goodness-of-fit belongs to a particular outcome, specification, and dataset, so the contrast is exactly why we simulate from the model we actually fitted.

5 Lens 2: The Blockmodel: Do Actors Come in Kinds?

What question does it answer? “Can a small number of relational roles summarize which organizations cooperate with one another?” This is the question from Day 9’s second half:

set.seed(6886)                    # VEM restarts draw random inits
A0 <- A; diag(A0) <- 0
bm <- BM_bernoulli("SBM_sym", A0, verbosity = 0, plotting = "")
bm$estimate()
K_star <- which.max(bm$ICL)

K_star

#> [1] 2

Z <- apply(bm$memberships[[K_star]]$Z, 1, which.max)
table(Z)

#> Z
#>  1  2 
#> 28  2

# VEM label order is not guaranteed across machines; index by size,
# never by raw label, or a label swap silently inverts the story:
tiny <- as.integer(names(which.min(table(Z))))
big  <- setdiff(seq_len(K_star), tiny)
round(bm$model_parameters[[K_star]]$pi, 2)

#>      [,1] [,2]
#> [1,] 0.11 0.74
#> [2,] 0.74 0.96

rownames(A0)[Z == tiny]

#> [1] "ANF"  "ASIM"

old_par <- par(mfrow = c(1, 2), mar = c(4.5, 4.5, 3, 1))
set.seed(6886)
g_block <- igraph::graph_from_adjacency_matrix(A0, mode = "undirected", diag = FALSE)
hub_label_angle <- rep(0, nrow(A0))
hub_label_angle[rownames(A0) == "ANF"] <- pi
plot(g_block, vertex.color = ifelse(Z == tiny, "#18453B", "#A7A8AA"),
     vertex.size = ifelse(Z == tiny, 15, 7),
     vertex.label = ifelse(Z == tiny, rownames(A0), NA),
     vertex.label.color = "#18453B", vertex.label.cex = 0.9,
     vertex.label.dist = ifelse(Z == tiny, 2.7, 0),
     vertex.label.degree = hub_label_angle,
     edge.color = "grey75", main = "Observed network by fitted role")
role_order <- c(big, tiny)
role_probability <-
  bm$model_parameters[[K_star]]$pi[role_order, role_order, drop = FALSE]
image(1:2, 1:2, role_probability,
      axes = FALSE, col = rev(hcl.colors(20, "Greens 3")), zlim = c(0, 1),
      main = "Fitted cooperation probability")
axis(1, at = 1:2, labels = c("Periphery", "Core"))
axis(2, at = 1:2, labels = c("Periphery", "Core"), las = 2)
heatmap_cells <- expand.grid(x = 1:2, y = 1:2)
heatmap_values <- as.vector(role_probability)
text(heatmap_cells$x, heatmap_cells$y,
     labels = sprintf("%.2f", heatmap_values),
     col = ifelse(heatmap_values > 0.55, "white", "black"),
     cex = 1.2)

Two panels. The network panel highlights ANF and ASIM as a small fitted core. The heatmap reports cooperation probabilities of 0.11 within the periphery, 0.74 between core and periphery, and 0.96 within the core. — The fitted two-block summary. Left: the observed cooperation network, with ANF and ASIM in the two-organization core. Right: the fitted probability of cooperation for each role pairing. The model finds a hub-versus-periphery split, not two cohesive communities.

par(old_par)

ICL selects 2 blocks. One contains ANF and ASIM, the two organizations with the most cooperation partners. Their fitted probability of cooperating with each other is 0.96, and their fitted probability of cooperating with a member of the other block is 0.74. The remaining 28 organizations have a fitted within-block probability of only 0.11. The model has found a core-periphery summary: two highly connected organizations and a large, sparse remainder. It has not uncovered two broad communities with dense cooperation inside each. The hard partition leaves only one possible pair within the two-organization core. The variational estimator also uses the organizations’ soft membership probabilities, so the 0.96 estimate is not literally calculated from that one dyad alone. Even so, a role represented this clearly by only two organizations is fragile and should be treated as a coarse summary rather than a well-estimated population rate.

This result develops the blockmodel lesson from Day 9 without claiming it is the same fit. Day 9 used blockmodels on a 31-organization binary joint-operation network and introduced NetMix separately on a longitudinal alliance panel. Here we fit a plain two-block model to the 30-organization binary cooperation outcome used by every capstone lens. A plain SBM has one cooperation rate for every role pairing and no organization-specific activity term. It can isolate the two largest hubs, but it cannot represent the smoother differences among the other 28 organizations. The useful conclusion is not “there are two rebel types.” It is that a two-role model compresses this network into a tiny core and a heterogeneous periphery, and that compression is too coarse for the remaining actor variation. A degree-corrected blockmodel could preserve discrete roles while adding actor-specific activity, but we would need to fit and check that model before claiming it solves this miss.

6 Lens 3: AME: A Relational Surface on Top of the Covariates

What question does it answer? “After the covariates and the additive effects, is there residual third-order structure, and which pairs have unexpectedly high or low compatibility?” This is Day 10’s workhorse, and the longer fitted chain is included in the cache:

gade_ame <- cache_fit("gade_ame2_lame135", ame(
  Y = dl$Y, Xdyad = dl$Xdyad, Xrow = dl$Xrow,
  family = "binary", R = 2, symmetric = TRUE,
  seed = 6886, burn = 10000, nscan = 50000, odens = 10,
  plot = FALSE, verbose = FALSE))

ame_beta_ess <- ess(gade_ame)
ame_beta_ess

#>             intercept averageId_actor1_node      size_actor1_node 
#>                   281                   197                   294 
#>       ideol_diff_dyad        powerdiff_dyad              loc_dyad 
#>                   679                   248                   243 
#>            spons_dyad 
#>                  1765

ame_ci <- t(apply(gade_ame$BETA, 2, quantile, c(0.025, 0.5, 0.975)))
round(ame_ci, 3)

#>                         2.5%    50%  97.5%
#> intercept             -9.185 -6.435 -4.014
#> averageId_actor1_node -0.366  0.122  0.480
#> size_actor1_node       0.068  0.167  0.263
#> ideol_diff_dyad       -0.746 -0.274  0.114
#> powerdiff_dyad        -0.458 -0.251 -0.121
#> loc_dyad               1.601  3.136  5.539
#> spons_dyad            -0.811 -0.029  0.747

The chain saves 5000 coefficient draws after thinning. Their effective sample sizes range from 197 to 1765 because adjacent MCMC draws remain correlated. These values describe mixing within this one coefficient chain; they do not establish agreement across starting values. The displayed summaries are adequate for an exploratory teaching comparison, but a research analysis should run independently initialized chains, compare posterior summaries, and extend any chain whose tail intervals remain unstable.

Put the SRM and AME intervals side by side. Shared location remains positive and power difference remains negative. Ideological distance is the result that changes: its interval is clearly negative under the additive SRM but widens to \([-0.75, 0.11]\) once the multiplicative term competes to explain the same pairs. The data do not cleanly separate “organizations closer in ideology cooperate more” from “cooperating organizations share an unmeasured relational profile.” Day 10 made the same identification point in the longitudinal ICEWS application: adding a latent surface can change the conditional association, but it does not reveal a unique causal explanation.

Day 9’s compare_fits() helper, pointed at today’s two fits: The table behind that paragraph in one call

compare_fits(gade_srm, gade_ame, names = c("srm", "ame"))

term	est_srm	est_ame	width_ratio	moved
intercept	-4.127	-6.472	1.76
averageId_actor1_node	0.176	0.103	2.29
size_actor1_node	0.117	0.166	1.67
ideol_diff_dyad	-0.256	-0.288	1.77	interval now includes zero
powerdiff_dyad	-0.113	-0.261	2.34
loc_dyad	1.828	3.238	1.89
spons_dyad	-0.478	-0.033	1.20

(The compare_fits() helper works on two lame fits with the same terms. The width_ratio column shows how the intervals change, and the moved column identifies terms whose interval now includes zero.)

The practical coefficient result is straightforward. Shared location remains associated with cooperation, but the ideology result becomes uncertain after the model is allowed to represent recurring partner patterns. The plotted factor surface is an exploratory posterior-mean illustration. It highlights the two highly connected organizations and contains a negative-eigenvalue direction that cannot be summarized by a simple map where nearby organizations are always more likely to cooperate.

lam <- diag(gade_ame$L)
round(lam, 1)

#> [1] 133.8 -42.0

Mu <- gade_ame$ULUPM; diag(Mu) <- NA    # Posterior-mean signed latent surface
top_ix <- which(abs(Mu) == max(abs(Mu), na.rm = TRUE), arr.ind = TRUE)
top_pair <- rownames(Mu)[top_ix[1, ]]
c(median_abs = round(median(abs(Mu), na.rm = TRUE), 2),
  max_abs = round(max(abs(Mu), na.rm = TRUE), 1))

#> median_abs    max_abs 
#>       1.36      40.70

top_pair

#> [1] "ASIM" "ANF"

U_display <- sweep(gade_ame$U, 2, sqrt(abs(lam)), `*`)
udf <- data.frame(u1 = U_display[, 1], u2 = U_display[, 2],
                  name = rownames(gade_ame$U),
                  ideol = dl$Xrow[, "averageId_actor1"],
                  size = dl$Xrow[, "size_actor1"])
ggplot(udf, aes(u1, u2, color = ideol, size = size)) +
  geom_point(alpha = 0.75) +
  ggrepel::geom_text_repel(aes(label = name), size = 3,
                           max.overlaps = 12, show.legend = FALSE) +
  scale_color_gradient(low = "#A7A8AA", high = "#18453B") +
  labs(x = "assortative factor score", y = "disassortative factor score",
       color = "avg ideology", size = "group size")

A two-dimensional factor-score plot of the Syrian organizations. ASIM lies far to the positive side of the assortative dimension, while ANF also has a relatively high assortative score. The vertical dimension represents a disassortative contribution. — Posterior-mean factor profiles from the rank-2 model, weighted by the square root of each eigenvalue’s absolute magnitude. The horizontal dimension contributes assortatively and the vertical dimension disassortatively. Read the signed fitted latent surface, not Euclidean distance or compass direction, as the relational result.

The largest entry in this posterior-mean latent surface belongs to ASIM–ANF, the same two organizations the blockmodel placed in its small core. This point estimate suggests that much of the fitted residual structure is concentrated around ANF and ASIM. The saved object contains coefficient draws but not latent-factor draws, so it does not supply intervals or effective sample sizes for this named-actor geometry. Treat the plot as exploratory rather than as evidence that every displayed position is stable.

Why a Mixed-Sign Factor Surface Can Represent More Than Distance

6.0.1 The Positive Axis

A positive eigenvalue means that organizations with scores in the same direction receive a higher latent compatibility score. This is an assortative pattern. A distance model can represent something similar by placing organizations with compatible profiles near one another.

6.0.2 The Negative Axis

A negative eigenvalue reverses that logic. Organizations at opposite ends of the axis can receive a higher compatibility score, while organizations on the same end can receive a lower score. This can represent organizations that relate similarly to a set of prominent partners without cooperating with one another. That is a form of stochastic equivalence without cohesion.

6.0.3 Why Distance Cannot Reproduce It

After the intercept and additive actor effects are removed, a Euclidean distance component is conditionally positive semidefinite. It cannot represent a disassortative negative-eigenvalue direction (Hoff 2007). The AME result is a signed inner product, not straight-line distance. The displayed axes are weighted by \(|\Lambda|^{1/2}\) so that they reflect their contribution to \(U\Lambda U^\mathsf{T}\), but their rotation, sign, and raw scale are not uniquely identified. In this small binary network, extreme probit-scale entries also warn us not to interpret their numerical magnitude literally. For research, save latent draws, compare invariant fitted surfaces across independently initialized chains, and base claims on prediction and goodness-of-fit rather than on the compass direction of a plotted axis.

round(gof_band(gade_ame$GOF)[, c("sd.rowmean", "trans.dep")], 3)

#>       sd.rowmean trans.dep
#> 2.5%        0.17    -0.054
#> 97.5%       0.20     0.060

The two selected scalar yardsticks fall inside the AME posterior-predictive bands, as they did for the SRM. The extra machinery nevertheless changes one coefficient conclusion: the ideological-distance interval includes zero once the latent surface competes to explain the same ties. This is useful specification sensitivity, not proof that the plotted latent positions are uniquely recovered.

7 Lens 4: The ERGM: Does a Shared-Partner Term Add Information?

What question does it answer? “Holding the rest of the observed network and the other included statistics fixed, does a proposed tie receive greater conditional odds when it adds shared partners or another specified configuration?” Day 11’s family asks a different question from the preceding latent-effects models. The latent models represent broad residual structure, while the ERGM assigns coefficients to named configurations. Neither move by itself identifies the behavioral process that produced those configurations.

A first specification includes the same covariates plus geometrically weighted shared partners:

snet <- to_statnet(net)
erg_base <- cache_fit("gade_ergm_base",
  ergm(snet ~ edges + edgecov("ideol_diff") + edgecov("powerdiff") +
         edgecov("loc") + edgecov("spons") +
         nodecov("averageId_actor1") + nodecov("size_actor1") +
         gwesp(0.5, fixed = TRUE),
       control = control.ergm(seed = 6886)))

eb <- summary(erg_base)$coefficients
round(eb[, c(1, 2, 5)], 3)

#>                          Estimate Std. Error Pr(>|z|)
#> edges                      -7.202      1.131    0.000
#> edgecov.ideol_diff         -0.252      0.121    0.037
#> edgecov.powerdiff          -0.087      0.029    0.002
#> edgecov.loc                 2.757      1.026    0.007
#> edgecov.spons              -0.378      0.384    0.326
#> nodecov.averageId_actor1    0.177      0.061    0.004
#> nodecov.size_actor1         0.077      0.020    0.000
#> gwesp.fixed.0.5             1.378      0.295    0.000

The question for this ERGM sequence is simple: does the fitted distribution still show a shared-partner association after we account for organizations that cooperate unusually widely? These are two patterns that can look similar. A process of choosing partners through partners can produce shared-partner clustering, but hubs also create many shared partners simply because they connect to many organizations. The coefficient describes the conditional network association. It does not show that organizations chose partners because they had partners in common.

In this first specification, gwesp is 1.38 with p 3e-06. Within this model, proposed ties that add shared partners receive more conditional weight. The model does not yet give every organization its own baseline propensity to cooperate, so this coefficient cannot tell us whether the shared-partner pattern remains after the two hubs are represented adequately.

But we built a yardstick precisely so that “taken at face value” is not the last step. Simulate from the fitted ERGM and score it:

sims_base <- simulate(erg_base, nsim = 200, seed = 6886)
gof_ergb <- t(sapply(sims_base, function(s) {
  m <- as.matrix(s); diag(m) <- NA; calc_gof(m)
}))
gof_summary <- function(G, observed = gof_obs) {
  rbind(
    observed = observed,
    simulated_median = apply(G, 2, median),
    `2.5%` = apply(G, 2, quantile, 0.025),
    `97.5%` = apply(G, 2, quantile, 0.975)
  )
}
round(gof_summary(gof_ergb), 3)

#>                  actor_sd trans
#> observed            0.184 0.007
#> simulated_median    0.139 0.088
#> 2.5%                0.117 0.058
#> 97.5%               0.158 0.144

The medians show the direction of the mismatch, and the 95% simulation intervals show its size. The observed actor-heterogeneity score is 0.184, above the simulated median of 0.139 and above the interval’s upper bound of 0.158. The observed transitivity score is 0.007, below the simulated median of 0.088 and below the lower bound of 0.058. This specification produces too little variation in organizations’ partner counts and too much transitivity on this measure.

The specification has no organization-specific activity terms, even though ANF and ASIM cooperate far more widely than the other organizations. Because hubs mechanically generate many shared partners, the gwesp term is being asked to help represent a network feature that the model has not described directly. The coefficient is the fitted conditional association implied by this specification, but its substantive interpretation is limited because the specification misses both focal network features.

7.1 First Try a Compact Degree Adjustment

A common next step is to add gwdegree(). This adds one curved statistic summarizing the network’s degree distribution. It is a compact way to represent degree heterogeneity, but it is not the same as giving each organization its own activity parameter.

erg_deg <- cache_fit("gade_ergm_deg",
  ergm(snet ~ edges + edgecov("ideol_diff") + edgecov("powerdiff") +
         edgecov("loc") + edgecov("spons") +
         nodecov("averageId_actor1") + nodecov("size_actor1") +
         gwdegree(0.5, fixed = TRUE) + gwesp(0.5, fixed = TRUE),
       control = control.ergm(seed = 6886)))

ed <- summary(erg_deg)$coefficients
round(ed[, c(1, 2, 5)], 3)

#>                          Estimate Std. Error Pr(>|z|)
#> edges                      -9.232      1.330    0.000
#> edgecov.ideol_diff         -0.228      0.111    0.039
#> edgecov.powerdiff          -0.097      0.030    0.001
#> edgecov.loc                 2.801      1.028    0.006
#> edgecov.spons              -0.305      0.393    0.437
#> nodecov.averageId_actor1    0.175      0.065    0.007
#> nodecov.size_actor1         0.083      0.022    0.000
#> gwdeg.fixed.0.5             2.501      0.886    0.005
#> gwesp.fixed.0.5             2.260      0.461    0.000

c(AIC_base = round(AIC(erg_base), 1), AIC_deg = round(AIC(erg_deg), 1))

#> AIC_base  AIC_deg 
#>    331.2    324.1

sims_deg <- simulate(erg_deg, nsim = 200, seed = 6886)
gof_ergd <- t(sapply(sims_deg, function(s) {
  m <- as.matrix(s); diag(m) <- NA; calc_gof(m)
}))
round(gof_summary(gof_ergd), 3)

#>                  actor_sd trans
#> observed            0.184 0.007
#> simulated_median    0.134 0.091
#> 2.5%                0.108 0.037
#> 97.5%               0.158 0.139

AIC is 7.1 points lower, but the observed actor-heterogeneity and transitivity scores remain outside their respective 95% simulation intervals. A lower AIC says that this degree term improves relative in-sample fit among these two specifications, which use the same outcome and covariate set. It does not show that the curved degree term has a substantive interpretation or that the model reproduces the network features we care about.

Do not interpret the change in gwesp from 1.38 to 2.26 as closure becoming stronger. Adding gwdegree() changes the model and therefore changes what the GWESP coefficient conditions on. More importantly, the new model still fails the two simulation checks. The useful conclusion is simply that one compact degree term did not solve the mismatch.

7.2 Give Every Organization Its Own Baseline

The next specification uses sociality() to give each organization its own baseline propensity to cooperate, analogous to the actor effect in the SRM. The network is undirected, so each organization receives one activity parameter rather than separate sender and receiver parameters.

The nodal main-effect covariates must come out of this specification. With one constrained sociality parameter for nearly every organization in this undirected network, the nodal main effects are linear combinations of the activity indicators and cannot be estimated separately. The dyadic covariates remain because they vary across pairs rather than only across organizations. This means that the change in GWESP cannot be attributed solely to replacing gwdegree() with sociality(): the conditioning set also no longer contains group size and average ideology.

erg_soc <- cache_fit("gade_ergm_soc",
  ergm(snet ~ edges + edgecov("ideol_diff") + edgecov("powerdiff") +
         edgecov("loc") + edgecov("spons") +
         sociality(nodes = -1) + gwesp(0.5, fixed = TRUE),
       control = control.ergm(seed = 6886)))

es <- summary(erg_soc)$coefficients
round(es["gwesp.fixed.0.5", c(1, 2, 5)], 3)

#>   Estimate Std. Error   Pr(>|z|) 
#>      0.184      0.465      0.692

sims_soc <- simulate(erg_soc, nsim = 200, seed = 6886)
gof_ergs <- t(sapply(sims_soc, function(s) {
  m <- as.matrix(s); diag(m) <- NA; calc_gof(m)
}))
round(gof_summary(gof_ergs), 3)

#>                  actor_sd  trans
#> observed            0.184  0.007
#> simulated_median    0.190 -0.018
#> 2.5%                0.172 -0.075
#> 97.5%               0.207  0.061

Read this third model on its own. Its conditional GWESP estimate is 0.18, with a 95% interval of [-0.73, 1.1]. The coefficient interval includes zero, so this specification does not detect a clear conditional shared-partner association after organization-specific activity is included.

For the simulation check, the answer to “is this fine?” is yes for these two statistics, but not as a general verdict on the model. The observed actor-heterogeneity score is 0.184, compared with a simulated median of 0.19 and a 95% interval of [0.172, 0.207]. The observed transitivity score is 0.007, compared with a median of -0.018 and an interval of [-0.075, 0.061]. Neither check detects a mismatch. The median is useful because it shows where a typical simulation lands; interval coverage alone only says that the observed statistic is not unusual in these 200 simulations.

This is still a deliberately saturated teaching diagnostic with 29 organization-specific sociality parameters and only 85 observed ties. Because those activity terms closely track the degree sequence, good actor-heterogeneity fit is partly built into the specification rather than an independent validation of it.

That gives us a focused substantive reading: within this highly flexible specification, the remaining conditional shared-partner association is imprecise, and the two selected simulation checks do not reveal a mismatch. This does not prove that closure never occurs and does not identify why ANF and ASIM became hubs. It also does not establish broad goodness-of-fit beyond the two statistics we checked.

There is also a tradeoff. The model uses 29 sociality parameters for 30 organizations, and it can no longer estimate nodal main effects such as group size. If group size is the main quantity of interest, this is not the appropriate final specification.

Diagnostics, Decay Choice, and What We Have Not Checked

All three cached ERGMs completed without an estimation failure. Rechecking the retained MCMC samples found no Geweke burn-in warning at the joint level, and fresh gof(..., GOF = ~model) simulations reproduced the statistics included in each fitted model. Those checks support using the fitted objects, but they do not make the final specification parsimonious or substantively correct.

We fix the gwesp() and gwdegree() decay at 0.5 so that the teaching sequence changes the representation of actor heterogeneity while keeping the geometric weighting convention constant. The value is a teaching convention, not a substantively established constant. A research analysis should examine plausible fixed values or estimate the decay when the model is stable enough to support that.

The custom check uses 200 simulations and only two summaries. Tail estimates from 200 draws are somewhat noisy, and passing these checks does not establish fit to the full degree distribution, edgewise shared partners, geodesic distances, components, isolates, triad census, or substantively important mixing. Use more simulations for a research result and add the relevant native ergm::gof() panels. The purpose here is narrower: follow the same two yardsticks across model families.

What This ERGM Sequence Does and Does Not Show

Four points matter. First, each coefficient belongs to its own specification. Do not subtract or rank the GWESP coefficients as though they estimated the same quantity under the same conditions. Second, the covariate set also changes when sociality() enters because the nodal main effects are no longer separately identifiable. Third, the final model handles the focal simulation mismatch using 29 organization-specific activity parameters, so matching degree heterogeneity is partly mechanical. Other degree specifications and decay values are possible. Fourth, a clear conditional shared-partner association would still describe a network pattern, not prove that organizations formed ties through a behavioral friend-of-a-friend process.

8 The Model Comparison Table: One Network, Every Compatible Cross-Sectional Lens

First, the yardstick scoreboard. The GLM’s simulations are already in hand from lens 0 (gof_glm, the day’s first NO); the one missing simulation is the blockmodel, scored from its fitted block-tie rates:

set.seed(6886)
Pi <- bm$model_parameters[[K_star]]$pi
Pb <- Pi[Z, Z]                            # Each pair's block-implied rate
gof_blk <- t(replicate(200, {
  S <- matrix(0, 30, 30)
  S[upper.tri(S)] <- rbinom(choose(30, 2), 1, Pb[upper.tri(Pb)])
  S <- S + t(S); diag(S) <- NA; calc_gof(S)
}))

band_txt <- function(x, d = 3) {
  q <- quantile(x, c(0.025, 0.5, 0.975))
  paste0(round(q[2], d), " [", round(q[1], d), ", ",
         round(q[3], d), "]")
}
covered <- function(x, obs) {
  q <- quantile(x, c(0.025, 0.975))
  ifelse(obs >= q[1] & obs <= q[2], "yes", "NO")
}
score <- data.frame(
  model = c("glm (independence)", "srm (additive effects)",
            "SBM (discrete roles)", "AME (additive + low-rank)",
            "ergm (gwesp)", "ergm (gwdeg + gwesp)",
            "ergm (sociality + gwesp)"),
  actor_sd = c(band_txt(gof_glm[, "actor_sd"]),
               band_txt(gade_srm$GOF[-1, "sd.rowmean"]),
               band_txt(gof_blk[, "actor_sd"]),
               band_txt(gade_ame$GOF[-1, "sd.rowmean"]),
               band_txt(gof_ergb[, "actor_sd"]),
               band_txt(gof_ergd[, "actor_sd"]),
               band_txt(gof_ergs[, "actor_sd"])),
  hit_a = c(covered(gof_glm[, "actor_sd"], gof_obs["actor_sd"]),
            covered(gade_srm$GOF[-1, "sd.rowmean"], gof_obs["actor_sd"]),
            covered(gof_blk[, "actor_sd"], gof_obs["actor_sd"]),
            covered(gade_ame$GOF[-1, "sd.rowmean"], gof_obs["actor_sd"]),
            covered(gof_ergb[, "actor_sd"], gof_obs["actor_sd"]),
            covered(gof_ergd[, "actor_sd"], gof_obs["actor_sd"]),
            covered(gof_ergs[, "actor_sd"], gof_obs["actor_sd"])),
  trans = c(band_txt(gof_glm[, "trans"]),
            band_txt(gade_srm$GOF[-1, "trans.dep"]),
            band_txt(gof_blk[, "trans"]),
            band_txt(gade_ame$GOF[-1, "trans.dep"]),
            band_txt(gof_ergb[, "trans"]),
            band_txt(gof_ergd[, "trans"]),
            band_txt(gof_ergs[, "trans"])),
  hit_t = c(covered(gof_glm[, "trans"], gof_obs["trans"]),
            covered(gade_srm$GOF[-1, "trans.dep"], gof_obs["trans"]),
            covered(gof_blk[, "trans"], gof_obs["trans"]),
            covered(gade_ame$GOF[-1, "trans.dep"], gof_obs["trans"]),
            covered(gof_ergb[, "trans"], gof_obs["trans"]),
            covered(gof_ergd[, "trans"], gof_obs["trans"]),
            covered(gof_ergs[, "trans"], gof_obs["trans"])))
score

model	actor_sd	hit_a	trans	hit_t
glm (independence)	0.119 [0.096, 0.14]	NO	0.022 [-0.021, 0.064]	yes
srm (additive effects)	0.179 [0.147, 0.207]	yes	-0.013 [-0.094, 0.066]	yes
SBM (discrete roles)	0.164 [0.14, 0.186]	yes	-0.055 [-0.107, -0.016]	NO
AME (additive + low-rank)	0.187 [0.17, 0.2]	yes	0.002 [-0.054, 0.06]	yes
ergm (gwesp)	0.139 [0.117, 0.158]	NO	0.088 [0.058, 0.144]	NO
ergm (gwdeg + gwesp)	0.134 [0.108, 0.158]	NO	0.091 [0.037, 0.139]	NO
ergm (sociality + gwesp)	0.19 [0.172, 0.207]	yes	-0.018 [-0.075, 0.061]	yes

Each simulation cell reports median [2.5th percentile, 97.5th percentile]. The median shows the center of the simulated distribution, while the interval shows whether the observed statistic is unusual under this particular check. These are separate percentile intervals for each statistic, not a simultaneous band that guarantees joint coverage across the whole table.

checkpoint(obs_actor_sd = round(gof_obs["actor_sd"], 3),
           obs_trans = round(gof_obs["trans"], 4),
           models_covering_both = sum(score$hit_a == "yes" &
                                      score$hit_t == "yes"))

#> ------------------------------------------------------------------
#> CHECKPOINT: obs_actor_sd = 0.184   |   obs_trans = 0.007   |   models_covering_both = 3
#> ------------------------------------------------------------------

Read the scoreboard carefully. The GLM covers this particular transitivity statistic. That is useful, but it does not show that the model reproduces every triangle-related feature or that its independence assumption is correct. The SBM also places the observed actor_sd at the edge of its simulation band, so two roles can approximate the broad hub-versus-periphery spread. Its simulations nevertheless miss the transitivity statistic in the opposite direction. A simulation check can reveal a clear miss. Coverage only means that this particular statistic was not unusual under the fitted simulations.

One comparability caveat matters: the SRM and AME rows are posterior-predictive bands that include parameter uncertainty, while the GLM, SBM, and ERGM rows simulate at fixed point estimates. The Bayesian bands are therefore wider by construction. The ERGM base and gwdegree misses remain clear because they miss in direction, not only by a small amount, but the SBM’s edge-of-band coverage should be read with this difference in mind. The ERGM family gets three rows because we watched its specification change. Moving from the GLM to the SRM adds one actor effect per organization, and adding sociality() gives the ERGM an analogous actor-specific adjustment. The plain SBM compresses actor heterogeneity into two coarse roles; we did not fit a degree-corrected SBM here.

The scoreboard aligns the fit checks. The next display aligns coefficient summaries without pretending that their numerical scales are the same. Every family prints its coefficient table differently: summary() matrices for the GLM and ERGM, posterior BETA draws for lame, and a robust-variance vector for DCR. _compare_helpers.R has one row_*() function per format, and each returns the same five columns: model, term, estimate, lower bound, and upper bound.

source("_compare_helpers.R")   # row_glm, row_dcr, row_ame, row_ergm, row_saom, forest_plot
rows <- rbind(row_glm(glm_fit),
              row_dcr(glm_fit, dcr_out),
              row_ame(gade_srm, "srm (additive)"),
              row_ame(gade_ame, "ame (R = 2)"),
              row_ergm(erg_base, "ergm (initial gwesp)"),
              row_ergm(erg_soc, "ergm (sociality + gwesp)"))
head(rows, 3)

model	term	estimate	lo	hi
glm (naive)	(baseline)	-6.1533602	-8.2844276	-4.0222928
glm (naive)	ideol_diff	-0.2529838	-0.5046566	-0.0013110
glm (naive)	powerdiff	-0.1186774	-0.1832987	-0.0540562

Warning

The horizontal lengths are not comparable across rows because the models use different scales and conditioning sets. Use this display only to track the direction of a result and whether each model’s interval includes zero, and only when the underlying covariate contrast is aligned.

forest_plot(rows)

A faceted coefficient-interval plot comparing aligned covariate results across the GLM, DCR, SRM, AME, initial ERGM, and sociality ERGM. Shared location is positive throughout, while the ideology interval includes zero in the AME and final ERGM. — Five families, one forest. Each facet is a covariate, each row a model, each segment a 95% interval on that family’s own scale. Read sign and interval-vs-zero across rows; never read size across families.

This is a specification-sensitivity display, not a method-performance comparison. The GLM and ERGM rows use logit coefficients, the lame rows use probit coefficients, and the ERGM conditions on the rest of the graph. A shorter segment is therefore not a smaller substantive effect. Read only the aligned qualitative results: shared location is positive in every row; ideological distance has an interval away from zero in the GLM, DCR, SRM, and initial ERGM, but its interval includes zero in both the AME and final sociality ERGM; and the DCR row keeps the GLM estimate while changing its uncertainty. The baseline, sociality, and gwesp terms are omitted because they have no cross-family analogue.

And the take-home part: each row_*() needs only the fitted object (plus, for row_dcr, the DCR output); point them at your fits and forest_plot() works unchanged. The file also has a row_saom() for your HW 4 sienaFit objects, with a loud comment you should take seriously: SAOM rate and evaluation parameters answer “how do ties change?”, not “which ties exist?”, so they get their own plot, never a facet beside the tie models.

Now the table that matters more than the scoreboard: what question each lens asked, and what it answered here:

Lens	Its question	Its answer on this network
Descriptives	Who is central, who clusters?	Two dominant hubs, ideological bunching, low density
GLM + DCR SEs	Which pair traits go with ties, with actor-shared uncertainty?	Same location, ideological similarity, overall ideological position, and combined group size have intervals separated from zero; the power-difference SE changes most, by 2.5x
SRM	How much do organizations differ in their overall tendency to cooperate?	A great deal: ANF and ASIM have unusually broad cooperation patterns, and ideological similarity remains associated with cooperation after actor effects enter
SBM	Do actors come in discrete kinds?	A two-organization core and a large periphery are useful as a coarse summary, but two roles cannot represent the remaining differences among organizations
AME	Is there residual pair-specific structure?	Yes, but it is concentrated: the latent surface again identifies the small hub core, and the ideology interval includes zero once this surface competes to explain the same ties
ERGM	Is a shared-partner term still informative beyond covariates and actor activity?	Not detected in the sociality specification; gwesp becomes imprecise and both yardstick statistics are covered, which is compatibility rather than proof of no closure
(SAOM / TERGM)	How do ties change over waves?	Not askable; this is a single snapshot; the wave models of Days 11–12 need panels

The fitted lenses produce several apparent tensions, and most become understandable once each answer is attached to its own question and specification. The apparent contradiction, “the ERGM found closure while the SRM says it is actor activity,” weakens when actor activity is represented in both models. That is a substantive reconciliation, not a guarantee that different defensible specifications will always agree. Where the answers converge, as several lenses do when they point to the same two hubs, that convergence is the strongest empirical claim the day produces. The decision framework turns that comparison into a repeatable process.

The Model Comparison in One Line

Many apparent disagreements shrink once we attach each result to the question it answers and the quantities it conditions on. The remaining differences, such as the SRM and AME interval verdicts for ideological distance, are specification sensitivity to report, not a contest with one automatic winner.

9 The Decision Framework

Here is the artifact to screenshot. Start from your question, not your data format, and walk down:

%%{init: {"theme": "base", "themeVariables": {"primaryColor": "#E8F0EC", "primaryTextColor": "#18453B", "primaryBorderColor": "#18453B", "secondaryColor": "#EFF6D8", "tertiaryColor": "#FFFFFF", "lineColor": "#18453B"}}}%%
flowchart LR
  Q["what is your question?"]
  Q --> D["describe /<br/>visualize"]
  Q --> M["measure a latent<br/>quantity"]
  Q --> B["a coefficient<br/>on a covariate"]
  Q --> S["is local structure<br/>a mechanism?"]
  Q --> T["how do ties<br/>change?"]
  Q --> C["causal effect<br/>of a treatment"]

  D --> D1["centrality, communities,<br/>netify summaries: Day 8"]

  M --> M1{"kinds or<br/>gradients?"}
  M1 -->|"discrete roles"| M2["blockmodel positions:<br/>Day 9"]
  M1 -->|"continuous profiles"| M3["latent distance or AME factors,<br/>plus additive effects: Day 10"]

  B --> B1{"model the dependence<br/>or just fix the SEs?"}
  B1 -->|"absorb it"| B2{"binary or valued?<br/>one net or panel?"}
  B2 --> B3["SRM / AME family:<br/>Days 9 and 10; lame handles<br/>valued + longitudinal"]
  B1 -->|"inference only"| B4["dyadic cluster-robust SEs:<br/>Day 14"]

  S --> S1["ERGM + GOF:<br/>Day 11"]

  T --> T1["SAOM: Day 12;<br/>TERGM: Day 11"]

  C --> C1{"can you defend<br/>interference + exposure<br/>assumptions?"}
  C1 -->|"yes (design)"| C2["design-based estimators:<br/>Day 13 literacy list"]
  C1 -->|"not really"| C3["observational claim:<br/>association + sensitivity;<br/>DCR handles shared-actor uncertainty"]

The decision framework. Every terminal box is a course day or a named frontier method with an entry point in the per-family reading table near the close.

And the same framework as a table, with the follow-up questions that pick the variant:

Your question	First tool	The follow-up questions that matter
Describe / visualize	Centrality, community detection, `netify` summaries	One network or many? Weighted?; Purely a description, no uncertainty
Measure a latent quantity (positions, roles, affinities)	Blockmodel (kinds), latent distance, or AME factors (gradients)	Discrete or continuous? Directed? Do you need uncertainty on the fitted latent quantity? Use Procrustes alignment only for coordinate models where rotation and reflection leave the fit unchanged
Estimate \(\beta\) while representing dependence	SRM/AME	Binary, ordinal, count? (`family=` in `lame`) · one snapshot or panel? (`lame` is longitudinal) · concerned that the conclusion depends on how dependence is represented? Compare a dyadic model with DCR and a structural model, while keeping their distinct conditioning sets clear
Adjust uncertainty without a structural dependence model	Dyadic cluster-robust SEs (`dcr`)	Is your estimand really just the conditional association?
A named network configuration	ERGM	Which local configurations encode the proposed dependence, and does GOF survive? Valued ties -> valued ERGM (ancillary)
Tie dynamics across observed waves	SAOM or TERGM	Actor-oriented micro-changes (SAOM, Day 12) or a tie-oriented network transition (TERGM, Day 11)? Co-evolution of behavior? Composition change?
Which observable features make actors influential, and on whom	SIR (influence regression); Section 11	Directed longitudinal counts? Can you name the observed features (\(W_r\)) that could carry the influence?
Causal effect	Day 13’s ladder	Randomized? Then exposure mappings and design. Observational? Say “we cannot defend the assumptions” when you cannot; and mean it

Use the following three examples to practice the framework. For each one, identify the primary branch, the strongest runner-up model, and the result that would change the choice.

The education researcher. Three waves of friendship nominations in two classrooms, plus each student’s reading score; the question is whether kids befriend similar-scoring peers or converge after befriending. Walk: the question is dynamic, selection versus influence, with panel waves, so the framework lands on the SAOM, and specifically on Day 12’s co-evolution machinery. Separate ERGMs for each wave would describe the networks but could not use the sequence of friendship and reading-score changes to distinguish the two stories.

The public-health researcher. One snapshot of patient-transfer counts among 120 hospitals, and the question is whether shared ownership predicts flows, controlling for capacity. Walk: estimate the ownership association while representing dependence; valued (counts, so a Poisson family in lame); directed. Use an SRM/AME with a count family, plus a DCR pass on the dyadic model to see whether the substantive conclusion survives both routes.

The political economist. Twenty years of sanctions and trade, and the draft’s abstract says “the effect of sanctions on trade.” Walk: the word effect sends you down the causal branch. The identification difficulty is specific: sanctions and trade can respond to the same unmeasured dyadic shocks, both networks evolve over time, and one state’s sanction policy can alter outcomes for other state pairs. Without a design or assumptions that address treatment assignment, interference, and time-varying confounding, the fitted network model estimates a conditional association rather than a causal sanctions effect. The defensible landing spot is to model that association with the dependence represented, run a sensitivity analysis tied to the main threats, and write “associated with” rather than “caused.”

The Framework’s One Rule

Before fitting, say whether the question concerns description, measurement, a covariate association, a network configuration, change over time, or a causal effect. That decision determines what the model must represent and what the result can mean.

9.1 Four Causal Lessons to Carry Forward

The four applications from Day 13 were different, but each gives us a question that travels to other research settings. Nickerson’s household voting study asks whether an intervention given to one person can affect someone connected to them. Paluck, Shepherd, and Aronow’s school study distinguishes a student’s own assignment, exposure through assigned peers, and the wider context of attending a program school. Ichino and Schündeln’s Ghana study asks whether an intervention suppresses behavior or simply moves it to nearby places. Egami and Tchetgen Tchetgen’s peer-GPA study asks whether connected people looked alike before any influence could occur.

These become four practical habits. First, distinguish a unit’s own assignment, indirect exposure through assigned neighbors, and broader setting-level exposure. Second, map every route through which treatment could reach an outcome. Third, look outside the treated unit or location for displacement. Fourth, test whether friendship, alliance, trade, or geographic proximity connected cases that already resembled one another. These habits do not require researchers to reproduce every specialized estimator from the original studies. They require researchers to define the comparison more carefully.

9.2 A Practical Checklist for Observational IR and Comparative Research

Most researchers cannot randomize sanctions, alliances, election monitoring, protest diffusion, or institutional reform. The useful lesson from Day 13 is therefore not “run an expensive experiment.” It is to make the comparison and its weaknesses explicit before choosing a model.

Start with six questions. What changes? Define the treatment or exposure precisely enough that another researcher could code it. Whose outcome may change, and when? Put the treatment before the outcome. Compared with what? Name the untreated units, earlier periods, or differently exposed cases supplying the comparison. Why might those cases already differ? List the causes that could affect both treatment and outcome. Can treatment reach a unit through its network neighbors? If it can, define the exposure rather than pretending every unit is simply treated or untreated. What can the study still teach us if the causal comparison is not credible? A careful conditional association, prediction, or descriptive network pattern can still matter.

9.2.1 Turn Those Questions Into a Modeling Plan

Consider the question, Do sanctions reduce later trade? Define treatment as a sanction imposed before the trade outcome, not as a vague indicator that the pair was ever sanctioned. Restrict the comparison to state-pair periods with meaningful overlap. Measure pre-sanction trade, security relations, regime ties, and economic conditions because they may influence both sanctions and later trade. Ask whether sanctions imposed by allies also affect the focal pair; if so, that is a network exposure that belongs in the design. Examine pre-treatment outcomes or credible negative controls for evidence that hidden selection remains.

Only then choose the outcome model. Repeated states and dyads call for an uncertainty calculation or dependence model suited to those repeated units. Temporal network dependence may call for a panel model. Latent actor positions may help proxy hidden selection when omitted traits leave a recoverable footprint in the observed network. None of those steps substitutes for the comparison. DCR can improve uncertainty for shared actors, and SRM or AME can represent unobserved relational structure, but neither makes sanctions as-if randomly assigned.

The same workflow carries to comparative research. For election monitoring, ask why observers were placed in particular locations before comparing monitored and unmonitored returns. For policy diffusion, define whether exposure means adoption by any neighbor, a proportion of neighbors, or a weighted set of influential neighbors. For protest or conflict spillovers, specify the geographic, organizational, or alliance network through which the spillover is supposed to travel. In every case, the network should enter because it defines exposure, selection, dependence, or the outcome, not simply because a network measure is available.

9.2.2 What the Headache Example Teaches About Negative Controls

Egami and Tchetgen Tchetgen (2024) study whether students do better academically when their friends have higher GPAs. The outcome they care about is the student’s later GPA. Headaches are not a second outcome of substantive interest. They provide one clue about unmeasured background factors shared through friendship or school context.

The problem is that students choose their friends. Shared motivation, family resources, stress, school conditions, or other unmeasured background factors could affect both which friends a student has and the student’s later GPA. Friends’ headaches should not directly change the student’s later GPA, but they may carry traces of that same background. If friends’ headaches predict the student’s later GPA, the analysis is finding a relationship that the peer-GPA story cannot explain.

The student’s earlier GPA provides a second clue. It reflects background influences on academic performance, but friends’ GPA cannot travel backward in time and change a GPA that was already recorded. Friends’ headaches therefore give a clue from the peer environment, while the student’s earlier GPA gives a clue from the student’s outcome history.

The basic lesson is diagnostic: choose a variable that may share hidden causes with the focal exposure or outcome but should not itself produce the focal effect. If friends’ headaches still predict later GPA, or friends’ GPA appears to predict a GPA recorded earlier, either unmeasured network confounding remains or the proposed negative control does not satisfy its exclusion assumptions. Egami and Tchetgen Tchetgen show that, under additional assumptions, the two clues can also help adjust the peer-GPA estimate. That procedure is more demanding than simply adding headaches and earlier GPA to an ordinary regression.

If the comparison remains weak after those checks, scale the claim to what the analysis can support. Report the modeled association, show how it changes across defensible specifications and negative controls, and identify the unmeasured process that could still explain it. That is more useful than attaching causal language to a model that solved only the dependence problem.

10 The Frontier: Node2vec Turns Network Position Into Predictors

Suppose we are interested in prediction rather than a coefficient or a complete model of how a network formed. We might want to predict a hidden cooperation tie, classify an organization, or find actors with similar network positions. Node2vec gives a prediction model a compact summary of where each actor sits in the network.

We can make that concrete with the Syrian armed-organization network used throughout this walkthrough. The network contains 30 organizations and 85 recorded cooperation ties. Al-Nusra Front, or ANF, has 20 cooperation partners, Ahrar al-Sham Islamic Movement, or ASIM, has 25, and the two organizations share 18 partners. Their local neighborhoods overlap substantially, although their complete network positions are not identical.

We now run the prediction exercise. In each of 30 splits, we hide 17 recorded cooperation ties and 17 non-ties. The hidden ties are removed before any network features are constructed. We compare three logistic prediction models: simple degree and shared-neighbor summaries, low-dimensional walk coordinates, and both sets of predictors together.

This compact demonstration compares network-derived features with one another. It does not include ideology, power, location, or sponsorship in the fitted baseline. A research application should add those observed pair characteristics and ask whether the walk coordinates improve prediction beyond both the simple network summaries and the measured predictors.

The coordinates use node2vec’s biased-walk idea with \(p=1\) and \(q=0.5\). For transparency, the code turns the resulting walk co-occurrences into a positive pointwise mutual information table and compresses that table with SVD. This is a matrix-factorization implementation of the walk-embedding idea, not the original skip-gram software implementation. The distinction matters for exact replication, although both procedures use sampled network contexts to produce actor coordinates.

Open the repeated held-out prediction exercise

source("_node2vec_demo.R")
node2vec_cv <- cache_fit("gade_node2vec_cv", {
  node2vec_holdout(
    as.matrix(net),
    repeats = 30,
    holdout_ties = 17,
    rank = 4,
    p = 1,
    q = 0.5,
    seed = 6886
  )
})
node2vec_summary <- data.frame(
  model = c(
    "Degree and shared-neighbor summaries",
    "Walk coordinates",
    "Both sets of predictors"
  ),
  mean_auc = round(
    colMeans(node2vec_cv[c("baseline", "walk_coordinates", "combined")]),
    2
  ),
  split_sd = round(
    vapply(
      node2vec_cv[c("baseline", "walk_coordinates", "combined")],
      sd,
      numeric(1)
    ),
    2
  )
)
node2vec_summary

	model	mean_auc	split_sd
baseline	Degree and shared-neighbor summaries	0.82	0.06
walk_coordinates	Walk coordinates	0.60	0.12
combined	Both sets of predictors	0.79	0.07

AUC equals 0.50 for chance ranking and 1.00 for perfect ranking. The simple summaries average 0.82; the walk coordinates average 0.60; and combining them averages 0.79. On this small, hub-dominated network, the walk coordinates do not improve on degree and shared-neighbor information. That negative result is useful: the embedding does not earn its extra complexity in this application.

10.1 What Node2vec Does, One Step at a Time

Build the test set first. Reserve some ties and non-ties for evaluation, and remove the reserved ties before learning any network features.
Take many short computer-generated trips through the remaining network. A trip starts at one organization and follows a few observed cooperation ties.
Notice which organizations repeatedly occur in similar network neighborhoods. If trips starting from ANF and ASIM often pass through many of the same organizations, their summaries will tend to look similar.
Place each organization in a small learned space. The short row of numbers is the organization’s coordinates, much like factor scores. Organizations encountered in similar network neighborhoods tend to receive nearby coordinates.
Turn actor summaries into a prediction. Combine the two organizations’ numerical summaries into pair-level predictors, fit a model on the visible pairs, and test it on the pairs hidden in step 1.

The computer-generated trips are a way to sample network neighborhoods. They are not organizations moving through the conflict, information spreading, or one group influencing another. The learned coordinates describe relative network position under node2vec’s definition of neighborhood similarity. Distances, dot products, or pairwise combinations of the coordinates can become predictors. The individual axes can rotate or change signs and are not automatically ideology, power, or alliance commitment, so their value must be judged by the task they help perform.

The useful sentence after the model runs

If walk-based features improve prediction on ties that were held out before the features were created, we can say that the larger network contains useful predictive information beyond the baseline variables. Here they do not, so the appropriate conclusion is that the simple network summaries are the better predictive choice for this task.

10.2 How This Fits Into the Machine-Learning Side of Network Analysis

Each method produces a short coordinate row for every observed actor. What changes is the network information used to decide which actors should receive similar rows.

Method	The practical question it uses to define similarity
DeepWalk, Perozzi et al. (2014)	Which actors repeatedly appear near one another during ordinary short walks?
LINE, Tang et al. (2015)	Which actors are directly connected, or have similar immediate neighbor lists?
Node2vec, Grover and Leskovec (2016)	Should a walk stay local and emphasize the same neighborhood, or travel outward and emphasize similar roles in different neighborhoods?

For the Syrian network, DeepWalk would place organizations close when ordinary walks repeatedly encounter them in similar contexts. LINE would emphasize direct cooperation and immediate partner lists. Node2vec can be tuned toward organizations in the same local coalition or toward organizations that occupy similar hub or brokerage roles in different portions of the network. These are alternative definitions of network similarity, not three estimates of the same underlying quantity.

10.2.1 NetMF and GraphSAGE Solve Two Different Practical Problems

NetMF makes the walk calculation visible. Build a large table in which row \(i\), column \(j\) records how often walks from actor \(i\) reach actor \(j\), relative to how often \(j\) would appear by chance. Actors with similar rows have similar sampled neighborhoods. SVD then compresses each long row into a few coordinates. NetMF is useful for understanding what the walk-based methods are approximately compressing.

GraphSAGE handles a new actor. Ordinary node2vec stores one fitted coordinate row for every actor used in training. It has no stored row for an organization that enters later. GraphSAGE instead learns a reusable recipe: combine the actor’s observed traits with summaries of its observed neighbors, then turn that information into coordinates. If a new armed organization appears and we observe its traits and initial partners, the same recipe can construct its features.

10.3 From Actor Features to a Research Result

The numerical summaries are not the result. They become useful only after we connect them to a clearly defined outcome and test them on cases the algorithm did not see.

Research task	A concrete example	What a credible test requires
Predict a hidden relationship	In the Syrian organization network, do two organizations’ network summaries help identify held-back cooperation ties?	Reserve ties and non-ties, remove test ties before learning the summaries, and compare against degree, shared partners, and measured dyadic variables
Predict an actor characteristic	Brown et al. (2021) use networks of links among Twitter users and internet domains to predict ideological labels	Keep the test actors in the graph so ordinary node2vec can embed them, hide their labels from the classifier, and compare against useful non-network predictors
Generalize to a new actor	GraphSAGE can summarize an actor not used in training by combining its measured traits with information from its observed neighbors	Hold out entire actors or later periods, rather than randomly hiding entries from actors already seen

For link prediction, we need one row per possible pair. A simple row can contain the dot product of the two actors’ node2vec summaries, their coordinate-by-coordinate differences, or both. We then fit an ordinary prediction model on the visible pairs and ask how well it identifies the pairs that were hidden before node2vec was trained. The exact pair construction matters less than preserving that order.

These are three different jobs. Success at predicting hidden ties does not establish that the same features predict ideology, and success on actors already represented in training does not establish that the method will work for a new organization or a later period.

10.3.1 NetMF Makes the Day 10 Connection Visible

Imagine a very large table. Each row is an organization. Each column represents a possible network context, and each cell records how strongly short trips from that organization reach that context. Two organizations with similar rows repeatedly lead us into similar parts of the network.

Keeping that entire table would be cumbersome. The same compression idea used in the Day 10 movie example lets us approximate it with a few numbers per organization. Qiu et al. (2018) showed that DeepWalk, LINE, and node2vec can be understood through this kind of matrix-factorization lens. NetMF makes the matrix view explicit for important members of this family.

10.4 Node2vec and the Day 10 Latent Factor Model Are Cousins

Both methods compress complicated relational patterns into a few numbers per actor. The important difference is what they compress and why.

Question	Node2vec and related graph ML	Latent factor model or AME
Main purpose	Build actor features for a prediction task	Model a declared dyadic outcome
What supplies the information	Network neighborhoods created by the walk procedure	Unexplained patterns in observed dyadic outcomes
Where measured predictors usually enter	In a later prediction model	Directly in the dyadic outcome model
What a good result looks like	Better performance on properly held-out cases	Useful coefficient estimates, uncertainty, prediction, and network goodness-of-fit

10.4.1 The Technical Connection, and Where It Stops

Node2vec first converts the observed graph into a training dataset of node-context pairs. Biased random walks generate sequences, and nodes that occur within a chosen window of one another become positive context pairs. The algorithm then learns a source vector \(z_i\) and a context vector \(c_j\) so observed context pairs receive large dot products and sampled alternatives receive small ones. With negative sampling, a common objective is

\[ \max_{Z,C} \left[ \sum_{(i,j)\in\mathcal D^+}\log \sigma(z_i^\mathsf T c_j) + \sum_{(i,k)\in\mathcal D^-}\log \sigma(-z_i^\mathsf T c_k) \right]. \]

Here, \(\mathcal D^+\) is created by co-occurrence in the walks, \(\mathcal D^-\) contains sampled non-context pairs, and \(\sigma(s)=1/(1+\exp(-s))\). This objective is a classifier for sampled walk contexts. It is not a Bernoulli likelihood for the original network tie \(Y_{ij}\), even though both use dot products and logistic functions.

Qiu et al. (2018) show why this still has a matrix-factorization interpretation. For DeepWalk and related negative-sampling methods, the dot products approximate entries in a transformed node-context matrix built from powers of a random-walk transition matrix, averaged across the context window and adjusted for node frequencies and negative sampling. Schematically,

\[ z_i^\mathsf T c_j \approx M_{ij}^{\text{walk}}, \]

where \(M^{\text{walk}}\) summarizes how often a walk from \(i\) reaches the context of \(j\) relative to its baseline frequency. It is not simply the adjacency matrix \(A\). Node2vec changes the transition process through its return and exploration settings, so its implicit matrix depends on a second-order walk and is not identical to the closed-form DeepWalk matrix used by NetMF.

The Day 10 AME model starts from a different statistical object:

\[ g\{\mathbb E(Y_{ij}\mid X,a,b,U,V)\} = x_{ij}^{\mathsf T}\beta+a_i+b_j+u_i^{\mathsf T}v_j. \]

For a binary network, \(g\) is a logit or probit link and the likelihood is written for the observed dyadic outcomes \(Y_{ij}\). The additive effects represent broad sending and receiving propensities, while \(u_i^\mathsf T v_j\) represents residual pair-specific compatibility after the measured predictors and additive effects enter. Bayesian AME estimates a joint posterior for these quantities; point-estimation versions optimize the corresponding outcome loss.

The genuine family resemblance is the low-rank dot product. Both approaches replace a large relational table with a small number of coordinates per actor, both are unchanged by compatible rotations of their coordinate systems, and both can use dot products as pair-level scores. The distinction is what the dot product reconstructs. Node2vec reconstructs sampled network-context proximity. AME reconstructs the conditional mean of a declared dyadic outcome after covariates and additive actor effects.

That distinction changes interpretation and evaluation. Node2vec coordinates are normally treated as generated predictors and judged on a held-out task. Standard errors for a later regression do not automatically account for the fact that the coordinates were estimated from the same network. AME supplies model-based uncertainty under its declared likelihood or posterior, but that uncertainty is only as credible as the outcome family, dependence specification, prior choices, and convergence. Neither method makes an unidentified covariate effect causal.

The comparison also changes for directed networks. Node2vec’s source and context vectors can be different, much like the \(U\) and \(V\) sides of a directed factor model, but their roles are defined by walk prediction rather than by sender and receiver effects in an outcome model. Standard node2vec is transductive: it stores coordinates for actors observed during training. GraphSAGE is inductive because it learns an aggregation rule that can construct coordinates for a new actor from available attributes and neighbors.

The Day 10 model writes the expected dyadic outcome as a function of measured predictors, broad actor activity, and a low-rank relational term. That setup is useful when we want to interpret coefficients, separate actor-level heterogeneity from pair-specific compatibility, and examine model-based uncertainty. Node2vec is useful when the immediate goal is scalable feature construction for a clearly defined prediction task.

GraphSAGE and many later graph neural networks push the machine-learning route farther by learning a rule for combining an actor’s measured traits with information from nearby actors. Unlike standard node2vec, that rule can be applied to a new actor if its traits and neighbors are available. The added flexibility can help prediction, but it also makes careful held-out evaluation and strong simple baselines even more important.

The warranted applied claim is narrow. If the features improve held-out cooperation prediction, say that the network contains predictive information not captured by the baseline variables. Do not name an embedding axis “ideology,” infer influence from proximity, or treat predictive performance as an explanation of why ties formed.

11 The Frontier: Social Influence Regression

Most of the fitted models in Days 9 through 12 focused on who connects to whom or how those ties change. Social influence regression, or SIR, asks a different longitudinal question: whose prior activity helps predict whose later activity?

Minhas and Hoff (2026) apply SIR to monthly material-conflict events among countries recorded by ICEWS. Material conflict includes physical attacks, destruction of property, and other coercive actions. The outcome is a directed count: how many such events state \(i\) directs toward state \(j\) during a month.

11.1 Separate the Pair’s Own History From the Wider Network

The direct part of the model describes the focal state pair. It includes conflict from \(i\) to \(j\) in the previous month, conflict in the reverse direction as a measure of possible retaliation, geographic distance, joint democracy, alliance status, trade, and verbal cooperation. These terms ask familiar questions such as whether recent conflict persists, whether states retaliate, and whether nearby or allied states have different expected conflict levels.

SIR then adds a wider network question. Suppose state \(k\) directed conflict toward state \(j\) last month. Does knowing that state \(i\) is allied with \(k\), trades with \(k\), or exchanges cooperative statements with \(k\) help predict whether \(i\) directs conflict toward \(j\) this month? That is sender-side influence. It asks whether related senders develop similar target patterns over time.

Receiver-side influence reverses the perspective. Suppose a source directed conflict toward state \(k\) last month. If state \(j\) is geographically close to \(k\), does that source also tend to direct conflict toward \(j\) later? This asks whether related targets receive conflict from similar sets of senders.

The model can be written as

\[ \eta_{ijt} = z_{ijt}^{\mathsf T}\theta + \left(A X_{t-1}B^{\mathsf T}\right)_{ij}, \qquad A=\sum_r\alpha_rW_r,\qquad B=\sum_s\beta_sW_s. \]

The term \(z_{ijt}^{\mathsf T}\theta\) contains the direct dyadic predictors. The matrix \(X_{t-1}\) records the previous month’s conflict. The matrix \(A\) describes which relationships connect senders, and \(B\) describes which relationships connect targets. Rather than estimating a separate influence weight for every possible country pair, SIR constructs these matrices from observed relationships such as alliances, trade, verbal cooperation, and geographic proximity.

11.2 What the Published Application Finds

The sender-side results show that states tend to direct more material-conflict events toward countries that their allies fought in the previous month. States also tend to align their targets with states with whom they exchange cooperative diplomatic statements. A country’s later conflict pattern therefore resembles the earlier conflict pattern of some of its allies and verbal-cooperation partners, even after the direct dyadic predictors enter.

On the receiver side, geographic proximity is especially important. Nearby states tend to receive conflict from similar sets of senders. This does not mean that proximity itself causes conflict. It means that the fitted model finds a recurring temporal pattern in which sources that target one country are more likely to target nearby countries as well.

These results show what SIR adds beyond a standard lagged dyadic regression. The direct terms describe persistence, retaliation, and measured features of the focal pair. The influence terms describe how activity involving other states helps predict the focal pair’s next outcome.

Important

In this model, influence means a conditional lagged association. The results do not establish that allies coordinated, that one country copied another, or that a relationship caused conflict to spread. Those stronger claims require a design or assumptions that address confounding and strategic anticipation.

11.3 Let Covariates Structure What Would Otherwise Be Hidden

SIR belongs to a broader family of strategies that do not leave every latent role, position, or influence weight unrestricted. They use observed variables to organize part of the hidden relational structure. The connection is useful, but the three examples below parameterize different objects.

NetMix and role membership. In the Day 9 specification, formula.monad = ~ 1 allowed each state-year’s mixed-membership weights to be learned from its partner list without nodal predictors. NetMix also permits a formula such as formula.monad = ~ democracy + log_gdp. Those coefficients shift the distribution over roles for a state-year. In applied language, they ask whether observed state traits help explain why a state-year’s alliance portfolio resembles one relational role more than another. The block matrix still describes how the roles connect, and residual uncertainty in role membership remains. Covariates do not hard-code the blocks.

Schematically, for role \(g\) relative to a reference role,

\[ \log\frac{\pi_{it,g}}{\pi_{it,G}} = x_{it}^{\mathsf T}\gamma_g, \]

where \(\boldsymbol\pi_{it}\) is the state-year’s vector of role weights and \(x_{it}\) contains its nodal covariates. The fitted network still helps infer \(\boldsymbol\pi_{it}\); the covariates make systematic differences in those weights interpretable.

Austin, Linkletter, and Wu’s covariate-defined latent space. Their model makes an actor’s continuous latent position partly predictable from nodal covariates:

\[ z_i = Bx_i+\varepsilon_i. \]

The term \(Bx_i\) is the position predicted by observed attributes, while \(\varepsilon_i\) is the residual position learned from the network. Ties then depend on distances among the \(z_i\) values, together with random effects. This construction helps explain which attributes organize the latent social space and permits a position, and therefore link probabilities, to be predicted for a new actor whose covariates are observed. It does not make the residual latent position disappear or turn the covariates into literal spatial axes.

SIR and covariate-defined influence channels. SIR does not regress an actor’s latent position on covariates. Instead, it writes the two directed influence matrices as

\[ A=\sum_r\alpha_r W_r, \qquad B=\sum_s\beta_s W_s. \]

Each observed matrix \(W_r\) encodes a proposed route through which prior activity could matter, such as alliance, geographic proximity, or shared institutional membership. The estimated \(\alpha_r\) and \(\beta_s\) weights say which source-side and target-side routes carry the lagged association in the fitted dynamic outcome model. The full influence contribution is \(A X_{t-1} B^\mathsf T\), so a path is strong when its source-side and target-side features jointly receive weight.

The common principle is structured relational heterogeneity. NetMix uses covariates to explain mixtures over discrete roles. Austin et al. use covariates to explain continuous latent positions. SIR uses observed relational matrices to construct directed influence channels. These are not three implementations of the same estimand: role membership describes partner-list types, latent distance describes tie propensity through proximity, and SIR describes how lagged network activity predicts later activity.

11.4 How SIR Is Estimated

For the Poisson fit below, \(\mu_{ijt}=\exp(\eta_{ijt})\), and the parameters maximize the following log-likelihood up to an additive term that depends only on the observed counts:

\[ \ell(\theta,\alpha,\beta) =\sum_{t,i\ne j} \left[ Y_{ijt}\eta_{ijt} -\exp(\eta_{ijt}) \right]. \]

The default fitting engine alternates between two familiar Poisson-regression jobs. First, it holds the receiver-side weights fixed and updates the direct and sender-side coefficients. Second, it holds the updated sender-side weights fixed and updates the direct and receiver-side coefficients. It repeats those steps until the likelihood and fitted coefficients stop changing meaningfully. Because the influence term multiplies the sender and receiver pieces, different starting values can reach different solutions, so compare their achieved likelihoods and fitted predictions.

The separate scales of \(\alpha\) and \(\beta\) are not identified: multiplying one side by \(c\) and dividing the other by \(c\) leaves \(A X_t B^{\mathsf T}\) unchanged. The software uses the reporting convention \(\alpha_1=1\). The channel matrix \(C=\alpha\beta^{\mathsf T}\) and the fitted means do not change under that rescaling, so those are the safer quantities to compare. Actor-clustered uncertainty assigns each dyad-period score contribution to both endpoint actors before constructing a sandwich covariance. Maximum likelihood still supplies the point estimates; the sandwich changes only their estimated covariance.

After fitting, check whether several starts reach the same likelihood and fitted values, inspect uncertainty for the channel coefficients, compare predictions on later held-out months, and examine whether important temporal or network patterns remain in the residuals.

Optional: Verify the Estimator on Simulated Data

The following example uses simulated data so that the fitted coefficients can be compared with known values. It checks the estimation code, not the substantive ICEWS findings described above.

library(sir)   # remotes::install_github("netify-dev/sir")

dat <- sim_sir(m = 20, T_len = 120, p = 2, q = 2,
               family = "poisson", gain = 0.9, seed = 42)
fit_sir <- cache_fit("fit_sir_frontier", sir(
  dat$Y, W = dat$W, X = dat$X, Z = dat$Z,
  family = "poisson", seed = 6886
))

ci <- confint(fit_sir)      # Cluster-robust by default
sir_tbl <- data.frame(
  term  = names(coef(fit_sir)),
  est   = round(unname(coef(fit_sir)), 3),
  lo    = round(ci[, 1], 3),
  hi    = round(ci[, 2], 3),
  truth = round(c(dat$theta, dat$alpha[-1], dat$beta), 3),
  row.names = NULL
)
sir_tbl

term	est	lo	hi	truth
(Z) Z1	0.237	0.231	0.243	0.237
(Z) Z2	-0.362	-0.372	-0.352	-0.365
(alphaW) W2	0.410	0.401	0.419	0.411
(betaW) W1	-0.654	-0.660	-0.648	-0.651
(betaW) W2	0.417	0.408	0.426	0.418

On simulated data where we know the truth, recovery is essentially exact, the largest gap between estimate and truth in that table is 0.003. Two reads matter. First, the fitted influence matrix produces actor-specific influence weights:

# The influence matrix's [i, k] entry records source k's fitted
# contribution to actor i's future outgoing ties.
fr_A <- fit_sir$A; diag(fr_A) <- NA
top <- base::order(abs(fr_A), decreasing = TRUE, na.last = NA)[1:5]
cells <- arrayInd(top, dim(fr_A))
data.frame(source = cells[, 2], influenced = cells[, 1],
           weight = round(fr_A[top], 2))

source	influenced	weight
14	5	-3.32
2	15	-2.90
1	15	-2.70
1	5	2.62
1	14	-2.57

Second, you can compare fitted predictions with the lagged channel present and with the lagged input set to zero. This is a fitted-prediction ablation:

mu_full <- predict(fit_sir)
mu_off  <- predict(fit_sir, newdata = list(
  W = dat$W, X = array(0, dim = dim(dat$X)), Z = dat$Z))

With the lagged input zeroed out, mean fitted activity falls from 1.55 to 1.1, a difference equal to about 29% of the full fitted mean here. This tells us how much the fitted predictions change when one input channel is removed and everything else is held at its estimated value. It is not a variance decomposition, proof that the lagged channel generated that share of the observed counts, or a causal intervention effect.

Where SIR Sits on the Map

SIR deliberately has no contemporaneous latent term. Influence enters through the lagged state. Actor-clustered standard errors allow score contributions sharing endpoint actors to be dependent, but they do not repair an omitted conditional-mean term for reciprocity, degree heterogeneity, or other network structure. Use this specification when the lagged channel is the target. Reach for a richer mean model when contemporaneous relational structure is part of the question.

The Day 10 Rhyme, in Miniature

The bilinear scale is not identified: \((c\alpha,\, \beta/c)\) gives the same likelihood for any \(c \neq 0\). The package therefore reports the normalized convention \(\alpha_1=1\). That normalization changes the separate coordinates, but it leaves \(C=\alpha\beta^{\mathsf T}\) and the fitted means unchanged. Compare \(C\) or fitted predictions across solutions rather than comparing unnormalized \(\alpha\) and \(\beta\) vectors.

The Full SIR Toolbox (Optional Self-Study)

None of this is part of the core path; all of it runs against the cached fit. sir also provides tidy summaries, alternative uncertainty checks, forecasts, and plots:

# The broom trio: tidy() rows are exactly the shape the comparison
# section's _compare_helpers.R pattern wants; augment() is the
# per-dyad frame (fitted values + Pearson residuals) behind any
# residual map you would draw
sir::tidy(fit_sir)

term	component	estimate	std.error	statistic
(Z) Z1	theta	0.2374254	0.0028342	83.77016
(Z) Z2	theta	-0.3617677	0.0047938	-75.46649
(alphaW) W2	alpha	0.4097317	0.0042275	96.91999
(betaW) W1	beta	-0.6540213	0.0029166	-224.24289
(betaW) W2	beta	0.4166543	0.0043674	95.40066

sir::glance(fit_sir)[, c("nobs", "n_periods", "AIC", "converged")]

nobs	n_periods	AIC	converged
45600	120	122295.1	TRUE

head(sir::augment(fit_sir), 3)

	sender	receiver	time	.observed	.fitted	.resid	.resid_pearson
2	2	1	1	2	0.6629741	1.3370259	1.6420694
3	3	1	1	2	2.1622159	-0.1622159	-0.1103173
4	4	1	1	3	0.6264633	2.3735367	2.9988043

# A naive-vs-cluster comparison, as one cbind: lower bounds
# side by side. On this clean simulation they barely differ; on the
# package's ICEWS application the cluster-robust intervals are much
# wider, so inspect the uncertainty choice rather than assuming it away
round(cbind(cluster = confint(fit_sir)[, 1],
            naive   = confint(fit_sir, se.type = "classical")[, 1]), 3)

#>             cluster  naive
#> (Z) Z1        0.231  0.230
#> (Z) Z2       -0.372 -0.369
#> (alphaW) W2   0.401  0.401
#> (betaW) W1   -0.660 -0.660
#> (betaW) W2    0.408  0.409

# A quick code-path check: resample whole periods and refit the
# pipeline. This treats sampled periods as exchangeable; it does not
# preserve serial dependence. R = 25 is too small for research inference.
bs <- cache_fit("sir_boot",
                boot_sir(fit_sir, R = 25, type = "block", seed = 6886))
round(cbind(analytic_se = sir::tidy(fit_sir)$std.error,
            period_boot_se = bs$se), 4)

#>             analytic_se period_boot_se
#> (Z) Z1           0.0028         0.0042
#> (Z) Z2           0.0048         0.0041
#> (alphaW) W2      0.0042         0.0045
#> (betaW) W1       0.0029         0.0030
#> (betaW) W2       0.0044         0.0039

# Out-of-sample check: Rolling-origin cross-validation against a
# carry-forward baseline, then a genuine one-step forecast
cv <- cv_sir(fit_sir, initial = 110, horizon = 1)
round(cv$aggregate[c("rmse", "naive_rmse")], 2)

#>       rmse naive_rmse 
#>       1.19       2.58

fc1 <- forecast(fit_sir, h = 1, Z_future = dat$Z[, , , 120, drop = FALSE])
round(mean(fc1, na.rm = TRUE), 2)   # Mean predicted activity, period 121

#> [1] 1.52

In this quick teaching run, the period-resampling and sandwich standard errors differ by at most 0.0014. This bootstrap draws whole periods but does not preserve serial dependence across them, so it is an additional uncertainty check rather than a time-series block bootstrap. Twenty-five resamples are enough to check that the code path works, not enough to establish bootstrap stability. A research analysis should use several hundred or more resamples, choose a resampling scheme that matches the dependence being claimed, and verify that the reported intervals stop moving. The rolling-origin comparison is more directly useful here: the fitted model’s RMSE is 1.19 against the carry-forward baseline’s 2.58. Unlike the switch-off ablation, that comparison evaluates predictions on later held-out periods.

set.seed(6886)
p_sir <- plot_sir_network(fit_sir, matrix = "A", threshold = 2)
# The package draws arrows, but its default arrowheads land underneath the
# destination nodes. An end cap stops each line before the node so the
# direction remains visible.
p_sir$layers[[1]]$geom_params$end_cap <- ggraph::circle(4, "mm")
p_sir$layers[[1]]$geom_params$arrow <- grid::arrow(
  length = grid::unit(3, "mm"), type = "closed"
)
p_sir +
  ggraph::scale_edge_colour_gradient2(
    low = "#535054", mid = "#D9D9D9", high = "#18453B", midpoint = 0,
    guide = "none"
  )

A directed network showing the strongest fitted sender-side influence weights in the simulated SIR example. Spartan-green arrows are positive and dark-gray arrows are negative. — The strongest fitted sender-side influence relationships. Arrows show direction, Spartan green is positive, and dark gray is negative. Only relationships with an absolute fitted weight above 2 are drawn; omitted lines are smaller, not necessarily zero. Setting matrix = ‘B’ draws the receiver side.

The default threshold of 0.1 would draw 345 of the 380 possible directed relationships, leaving an unreadable web. The higher display threshold keeps the 21 strongest relationships so we can see their direction and sign. This is a display filter, not a new model or a significance test. An isolated node in this picture can still have many smaller fitted influence relationships.

12 The Frontier: Multilayer Networks

Everything so far has usually treated one relation as one outcome matrix. A multilayer network keeps several relations among the same actors so we can study what is shared across them and what is specific to each kind of relationship.

12.1 When Theory Calls for Layers Rather Than Several Unrelated Outcomes

The motivation should come before the data structure. A multilayer representation is useful when different relations among the same actors are theoretically connected but not interchangeable. There are at least three distinct claims.

Shared causes. Cooperation and conflict may both increase when states interact frequently, when a crisis attracts attention, or when events are heavily reported. The layers co-vary because they respond to common conditions. Depending on its specification, a joint model can represent shared cross-layer mean structure, residual association, or both. None of these patterns by itself shows that one layer causes another.

Cross-layer dependence. A prior relation in one layer may help predict a later relation in another. In international relations, conflict may precede reduced trade. In American politics, cosponsorship may precede greater voting agreement. In education, friendship may precede a study-help tie. With repeated observations, a model can include both the layer’s own lag and lags from other layers:

\[ g\{\mathbb E(Y_{ij\ell t})\} = \alpha_\ell +\rho_{\ell\ell}Y_{ij\ell,t-1} +\sum_{\ell'\ne\ell}\rho_{\ell\ell'}Y_{ij\ell',t-1} +\text{covariates and network structure}. \]

The cross-layer coefficient \(\rho_{\ell\ell'}\) asks whether relation \(\ell'\) at the previous observation helps predict relation \(\ell\) now, conditional on the included history. It is a temporal conditional association, not automatically a causal effect.

Co-evolution and feedback. Two processes may change one another over time. Conflict can reduce trade, while trade dependence may discourage later conflict. Cosponsorship can build voting coalitions, while repeated voting agreement can make future cosponsorship easier. Friendship can create opportunities for study help, while repeated help can strengthen a friendship. A co-evolution claim therefore requires repeated measurements of every focal layer and a model that allows more than one direction of dependence. A single cross-section shows overlap but cannot identify which process moved first.

The practical question is not merely, “Do I have several network variables?” It is, “Would collapsing these relations erase a theoretical distinction, and does my theory predict shared causes, cross-layer dependence, or feedback?” If the answer is no, separate outcome models may be clearer. If the answer is yes, the layer index belongs in both the data and the model.

12.2 Keep the Extra Index Visible

For a multilayer network, \(Y_{ij\ell}\) records the relation from actor \(i\) to actor \(j\) in layer \(\ell\). The four ICEWS layers below distinguish verbal cooperation, material cooperation, verbal conflict, and material conflict among the same states in 2010. With repeated years, the object would be \(Y_{ij\ell t}\), adding time as a fourth mode.

data(icews, package = "netify")
ic10 <- icews[icews$year == 2010, ]
lay_names <- c("verbCoop", "matlCoop", "verbConf", "matlConf")
layer_labels <- c(
  verbCoop = "Verbal cooperation",
  matlCoop = "Material cooperation",
  verbConf = "Verbal conflict",
  matlConf = "Material conflict"
)
layers <- lapply(lay_names, function(v) {
  netify(ic10, actor1 = "i", actor2 = "j", symmetric = FALSE,
         weight = v)
})
ml <- layer_netify(layers, layer_labels = lay_names)
dim(ml)

#> [1] 152 152   4

summary(ml)[, c("layer", "num_actors", "density", "num_edges")]

layer	num_actors	density	num_edges
verbCoop	152	0.4346462	9976
matlCoop	152	0.0927588	2129
verbConf	152	0.1110579	2549
matlConf	152	0.1155891	2653

The object contains 4 relations on the same 152 states. Verbal cooperation has at least one event on 43.5% of eligible dyads, while material cooperation has at least one on 9.3%. Flattening the layers into “any interaction” would make a verbal statement and a material conflict count as the same outcome. Fitting each layer separately would preserve the outcomes but would not estimate what the layers share.

12.3 Use `netify` to Keep the Layers Aligned

The important software job is to preserve one actor ordering, one missingness rule, and explicit relation labels. layer_netify() combines already checked netify objects. The resulting object can be inspected, sliced, compared, plotted, or converted for a model without rebuilding the actor-by-actor matrices by hand.

is_multilayer(ml)

#> [1] TRUE

n_layers(ml)

#> [1] 4

summary(ml)[, c("layer", "num_actors", "density", "num_edges")]

layer	num_actors	density	num_edges
verbCoop	152	0.4346462	9976
matlCoop	152	0.0927588	2129
verbConf	152	0.1110579	2549
matlConf	152	0.1155891	2653

# Pull out one layer as a netify object or as an adjacency matrix.
verbal_conflict <- subset(ml, layers = "verbConf")
A_verbal_conflict <- as.matrix(ml, layer = "verbConf")
dim(verbal_conflict)

#> [1] 152 152

dim(A_verbal_conflict)

#> [1] 152 152

plot(ml) produces a faceted plot from the aligned multilayer object. For the teaching graphic below, the same basic plot data are restricted to 12 high-activity states and the 12 highest-volume dyads in each relation so the partner changes remain readable. compare_networks() can compare edge values or attributes across layers using correlations or binary overlap measures, and it can return the underlying comparison matrices and optional permutation tests.

Compare the four aligned relation layers

compare_networks(
  ml,
  method = "jaccard",
  edge_threshold = 0,
  return_details = TRUE
)$details$jaccard_matrix

#>           verbCoop  matlCoop  verbConf  matlConf
#> verbCoop 1.0000000 0.1976848 0.2391175 0.2376519
#> matlCoop 0.1976848 1.0000000 0.2969227 0.3257555
#> verbConf 0.2391175 0.2969227 1.0000000 0.3853529
#> matlConf 0.2376519 0.3257555 0.3853529 1.0000000

Construction and comparison do not choose the scientific model. A researcher can extract one layer for an ERGM, convert aligned matrices for a tensor model, or retain both layer and time indices for a longitudinal multilayer model. The theoretical claim determines which route is warranted.

12.4 Read the Four Layers Side by Side

Four network panels show the same twelve states in fixed positions, with different directed ties for verbal cooperation, material cooperation, verbal conflict, and material conflict. — Four ICEWS relation layers among the same twelve high-activity states in 2010. Node positions are held fixed across panels.

The display keeps the same twelve high-activity states in the same positions and shows the twelve most active directed state pairs within each layer. Thicker arrows indicate more coded events within that panel. Because ties and thickness are selected separately within each layer, use the figure to compare which partners appear and how the pattern changes, not to compare raw density across panels.

The United States is prominent in all four panels, but its partners and the surrounding structure change across cooperation and conflict. Some state pairs appear across several layers, while others are visible only for one kind of interaction. A single network labeled “international interaction” would hide those differences before any model was estimated.

This is a descriptive picture, not evidence that cooperation causes conflict or that one layer comes first. Its job is to show why the layer label belongs in the data and in the research question.

12.5 Benchmark Overlap Against the Layer Densities

Jaccard overlap alone can be misleading when one layer is much denser than another. We therefore report two quantities. Jaccard is the intersection divided by the union. Overlap lift is the observed joint-tie rate divided by the joint-tie rate expected if the two layer indicators were independent while keeping their observed densities. A lift of 1 is the density-only benchmark.

layer_overlap <- compare_networks(
  ml, method = "jaccard", edge_threshold = 0, return_details = TRUE
)
round(layer_overlap$details$jaccard_matrix, 2)

#>          verbCoop matlCoop verbConf matlConf
#> verbCoop     1.00     0.20     0.24     0.24
#> matlCoop     0.20     1.00     0.30     0.33
#> verbConf     0.24     0.30     1.00     0.39
#> matlConf     0.24     0.33     0.39     1.00

overlap_benchmark <- function(dat, vars, thresholds = c(1, 5)) {
  pairs <- combn(vars, 2, simplify = FALSE)
  out <- lapply(thresholds, function(cutoff) {
    do.call(rbind, lapply(pairs, function(pair) {
      left <- dat[[pair[1]]] >= cutoff
      right <- dat[[pair[2]]] >= cutoff
      p_left <- mean(left)
      p_right <- mean(right)
      p_both <- mean(left & right)
      data.frame(
        threshold = paste0("At least ", cutoff, " event",
                           ifelse(cutoff == 1, "", "s")),
        layer_pair = paste(layer_labels[pair], collapse = " with "),
        density_left = p_left,
        density_right = p_right,
        jaccard = p_both / mean(left | right),
        independence_jaccard =
          (p_left * p_right) /
          (p_left + p_right - p_left * p_right),
        overlap_lift = p_both / (p_left * p_right)
      )
    }))
  })
  do.call(rbind, out)
}

overlap_results <- overlap_benchmark(ic10, lay_names)
transform(
  overlap_results,
  density_left = round(density_left, 3),
  density_right = round(density_right, 3),
  jaccard = round(jaccard, 3),
  independence_jaccard = round(independence_jaccard, 3),
  overlap_lift = round(overlap_lift, 2)
)

threshold	layer_pair	density_left	density_right	jaccard	independence_jaccard	overlap_lift
At least 1 event	Verbal cooperation with Material cooperation	0.435	0.093	0.198	0.083	2.16
At least 1 event	Verbal cooperation with Verbal conflict	0.435	0.111	0.239	0.097	2.18
At least 1 event	Verbal cooperation with Material conflict	0.435	0.116	0.238	0.100	2.10
At least 1 event	Material cooperation with Verbal conflict	0.093	0.111	0.297	0.053	4.53
At least 1 event	Material cooperation with Material conflict	0.093	0.116	0.326	0.054	4.77
At least 1 event	Verbal conflict with Material conflict	0.111	0.116	0.385	0.060	4.91
At least 5 events	Verbal cooperation with Material cooperation	0.242	0.020	0.078	0.018	4.00
At least 5 events	Verbal cooperation with Verbal conflict	0.242	0.033	0.133	0.030	4.06
At least 5 events	Verbal cooperation with Material conflict	0.242	0.037	0.145	0.034	3.90
At least 5 events	Material cooperation with Verbal conflict	0.020	0.033	0.285	0.012	18.04
At least 5 events	Material cooperation with Material conflict	0.020	0.037	0.273	0.013	16.68
At least 5 events	Verbal conflict with Material conflict	0.033	0.037	0.402	0.018	16.39

The clearest pair is verbal and material conflict. At the one-event threshold their Jaccard overlap is about 0.39 and their overlap lift is about 4.9. At the five-event threshold the Jaccard remains about 0.40 while the lift rises above 16 because high-volume dyads are much rarer. The same state pairs therefore appear in both conflict layers far more often than their separate densities would predict. That can reflect common escalation, common reporting attention, or persistent high-activity dyads. It does not by itself show that verbal conflict causes material conflict. The threshold comparison matters because an overlap conclusion that disappears under a modest threshold change is fragile.

Translate the table into dyad percentages. At the one-event threshold, 38.5% of dyads that appear in either the verbal-conflict or material-conflict layer appear in both. The lift of 4.91 says that joint appearance is 4.91 times as common as it would be if the two layer indicators were independent while retaining their observed densities. At the five-event threshold, the overlapping share remains about 40%, but the lift becomes much larger because high-volume ties are rare in both layers.

The independence benchmark is deliberately simple. It does not preserve state activity, target exposure, geography, time-varying reporting attention, or network structure. After finding overlap, an applied analysis should take four steps: repeat the comparison under substantively meaningful coding thresholds; compare against a null that preserves obvious actor or temporal structure; use a joint factor model if the question concerns shared profiles or a lagged multilayer model if the question concerns temporal prediction; and check whether the fitted model reproduces residual cross-layer overlap as well as within-layer degree and triad patterns.

12.6 A Joint Multilinear Model Asks What the Layers Share

Start with one component \(r\) before reading the full equation. The number \(u_{ir}\) describes how strongly state \(i\) expresses the source side of that profile. The number \(v_{jr}\) describes how strongly state \(j\) expresses its target side. Their product \(u_{ir}v_{jr}\) is a score for the ordered pair \(i\rightarrow j\). The layer loading \(w_{\ell r}\) then changes that pair score for relation \(\ell\):

\[ \text{component contribution for }(i,j,\ell) =u_{ir}v_{jr}w_{\ell r}. \]

If \(u_{ir}v_{jr}\) is large and positive, then a positive \(w_{\ell r}\) raises the fitted outcome for that source-target pair in layer \(\ell\). A loading near zero means the shared pair pattern contributes little to that layer. A negative loading reverses the pattern: pairs pushed upward in a positive-loading layer are pushed downward in the negative-loading layer. The magnitude of \(w_{\ell r}\) is meaningful only after a normalization because scale can move among \(U\), \(V\), and \(W\) without changing their product.

Make that concrete with the fitted ICEWS example below. One estimated component concentrates on high-activity dyads involving the United States. Another concentrates on a Middle East pattern involving Iran, Israel, Lebanon, Egypt, and Iraq. The layer loadings ask whether those same sets of source-target pairs recur in verbal cooperation, material cooperation, verbal conflict, and material conflict. In this fit both components appear strongly across cooperation and conflict. The model therefore finds recurring dyadic profiles, not a clean axis that means cooperation at one end and conflict at the other.

The full model adds measured predictors, broad actor effects, and more than one component:

\[ \eta_{ij\ell} = x_{ij\ell}^{\mathsf T}\beta_{\ell} +a_{i\ell}+b_{j\ell} +\sum_{r=1}^{R}u_{ir}v_{jr}w_{\ell r}. \]

Read the equation from left to right. \(x_{ij\ell}^{\mathsf T}\beta_\ell\) allows measured covariates to have different associations in each relation. For example, geographic distance could be more strongly related to material conflict than to verbal cooperation. The terms \(a_{i\ell}\) and \(b_{j\ell}\) allow state \(i\)’s overall source activity and state \(j\)’s overall target exposure to differ by layer. Removing those broad source and target tendencies first means that any remaining United States-centered component refers to particular pairs that are unusually active, not simply to the United States having high activity overall.

The summation is the shared relational part. Each \(r\) indexes one recurring source-target profile. The \(U\) and \(V\) scores are shared across layers, which keeps “component 1” tied to the same dyadic pattern everywhere. The \(W\) loadings show where that component is expressed. With \(R=2\), every fitted source-target-layer value receives two contributions, \(u_{i1}v_{j1}w_{\ell1}+u_{i2}v_{j2}w_{\ell2}\), on top of the measured predictors and additive effects.

The linear predictor \(\eta_{ij\ell}\) still needs an outcome model. A Gaussian model can use it as the conditional mean of a continuous or transformed outcome. A binary model maps it through a logit or probit link to a tie probability. A count model maps it through a log link to a positive expected count. The shared-factor idea is the same, but the likelihood, estimation routine, and interpretation of a one-unit change differ.

Four separate factor fits preserve the layer labels but solve four unrelated coordinate problems. Their axes can rotate, change signs, or reorder independently, so visually similar “component 1” plots do not establish that the same component recurs. A joint tensor fit constrains the layers to share \(U\) and \(V\) and estimates \(W\) explicitly. That is the additional empirical claim purchased by the joint model.

The estimator chooses \(\beta_{\ell}\), the additive effects, and the three factor blocks so the joint fitted tensor is close to all observed layers at once. Under a Gaussian working model, that means minimizing the sum of squared residuals over eligible actor pairs and layers, which is the same as maximizing the Gaussian likelihood. Alternating least squares holds two factor blocks fixed, updates the third by least squares, updates the additive terms, and cycles until the objective stops improving. Binary and count versions use weighted, variational, or Bayesian updates for their likelihoods. Multiple starting values matter because the objective is not jointly convex.

This is the same broad computational idea used by the fast estimator in lame. Holding the other latent blocks fixed turns one difficult joint problem into an ordinary regression update. The algorithm then cycles through the blocks, using a small ridge term for numerical stability and short inner updates for the additive effects. These conditional updates are designed to reduce the joint squared-error objective, but the full objective can have several local optima. That is why we compare starts and evaluate rank on held-out data. It is different from ERGM maximum pseudolikelihood: pseudolikelihood replaces the graph likelihood with a product of conditional tie contributions, while ALS keeps one declared outcome loss and changes which parameter block is being optimized. It is also different from lame’s MCMC option, which explores a posterior rather than returning only a point solution.

Coordinates are not unique. Components can be permuted, signs can flip, and scale can move among \(u\), \(v\), and \(w\) without changing the product. Interpret the fitted surface \(\sum_r u_{ir}v_{jr}w_{\ell r}\), normalized within-component layer loadings, held-out predictions, and replicated network features. Do not attach meaning to an unaligned axis direction.

For applied work, the minimum audit is: compare several ranks out of sample; rerun several starts; add the covariates needed by the substantive question; inspect residual actor activity, triads, and cross-layer overlap; and use a binary or count likelihood when the raw outcome requires it. The fast Gaussian fit below is an interpretable teaching approximation. It is not a substitute for uncertainty intervals when the research claim depends on a particular coefficient, loading, or dyadic prediction.

12.7 A Reproducible Two-Component ICEWS Fit

netify preserves the common actor ordering, relation labels, and overlap diagnostics. Its bridge to lame currently fits extracted layers separately, so the joint fit below makes the alternating updates explicit rather than pretending that a layer-by-layer call estimates cross-layer dependence. The chunk fits the Gaussian least-squares version to \(\log(1+\text{event count})\) after standardizing each layer. It jointly updates layer-specific sender and receiver effects with a rank-2 multilinear surface and omits measured dyadic predictors so the new cross-layer term stays visible. Rank 2 is a compact teaching choice, not a claim that two components are optimal. A research analysis should add defensible covariates, compare ranks with held-out likelihood or prediction, inspect residual network features, and check more starting values.

Open the joint multilinear alternating-least-squares fit

Y_tensor <- simplify2array(lapply(layers, as.matrix))
dimnames(Y_tensor)[[3]] <- lay_names
Y_work <- log1p(Y_tensor)
for (ell in seq_along(lay_names)) {
  current <- Y_work[, , ell]
  Y_work[, , ell] <-
    (current - mean(current, na.rm = TRUE)) /
    sd(current, na.rm = TRUE)
}

fit_additive_layer <- function(M, maxit = 20) {
  n <- nrow(M)
  mu <- mean(M, na.rm = TRUE)
  a <- b <- numeric(n)
  for (iteration in seq_len(maxit)) {
    a <- rowMeans(sweep(M, 2, b, "-") - mu, na.rm = TRUE)
    shift <- mean(a)
    a <- a - shift
    mu <- mu + shift
    b <- colMeans(sweep(M, 1, a, "-") - mu, na.rm = TRUE)
    shift <- mean(b)
    b <- b - shift
    mu <- mu + shift
  }
  fitted <- mu + outer(a, rep(1, n)) + outer(rep(1, n), b)
  list(mu = mu, a = a, b = b, fitted = fitted)
}

fit_tensor_start <- function(Y, rank = 2, seed = 6886,
                             maxit = 120, tolerance = 1e-7) {
  n <- dim(Y)[1]
  n_layers <- dim(Y)[3]
  additive <- lapply(seq_len(n_layers), function(ell) {
    fit_additive_layer(Y[, , ell])
  })
  residual <- Y
  for (ell in seq_len(n_layers)) {
    residual[, , ell] <- Y[, , ell] - additive[[ell]]$fitted
  }
  additive_only_sse <- sum(residual^2, na.rm = TRUE)
  eligible <- !is.na(Y[, , 1])
  eligible_index <- which(eligible, arr.ind = TRUE)

  set.seed(seed)
  U <- matrix(rnorm(n * rank), n, rank)
  V <- matrix(rnorm(n * rank), n, rank)
  W <- matrix(rnorm(n_layers * rank), n_layers, rank)
  for (component in seq_len(rank)) {
    norm_u <- sqrt(sum(U[, component]^2))
    U[, component] <- U[, component] / norm_u
    W[, component] <- W[, component] * norm_u
    norm_v <- sqrt(sum(V[, component]^2))
    V[, component] <- V[, component] / norm_v
    W[, component] <- W[, component] * norm_v
  }

  previous_sse <- NA_real_
  for (iteration in seq_len(maxit)) {
    for (i in seq_len(n)) {
      targets <- seq_len(n)[-i]
      design <- do.call(rbind, lapply(seq_len(n_layers), function(ell) {
        sweep(V[targets, , drop = FALSE], 2, W[ell, ], "*")
      }))
      U[i, ] <- solve(crossprod(design) + diag(1e-8, rank),
                     crossprod(design, as.vector(residual[i, targets, ])))
    }
    for (j in seq_len(n)) {
      sources <- seq_len(n)[-j]
      design <- do.call(rbind, lapply(seq_len(n_layers), function(ell) {
        sweep(U[sources, , drop = FALSE], 2, W[ell, ], "*")
      }))
      V[j, ] <- solve(crossprod(design) + diag(1e-8, rank),
                     crossprod(design, as.vector(residual[sources, j, ])))
    }
    for (ell in seq_len(n_layers)) {
      design <-
        U[eligible_index[, 1], , drop = FALSE] *
        V[eligible_index[, 2], , drop = FALSE]
      W[ell, ] <- solve(crossprod(design) + diag(1e-8, rank),
                       crossprod(design, residual[, , ell][eligible]))
    }
    for (component in seq_len(rank)) {
      norm_u <- sqrt(sum(U[, component]^2))
      U[, component] <- U[, component] / norm_u
      W[, component] <- W[, component] * norm_u
      norm_v <- sqrt(sum(V[, component]^2))
      V[, component] <- V[, component] / norm_v
      W[, component] <- W[, component] * norm_v
    }

    multilinear <- array(NA_real_, dim(Y))
    for (ell in seq_len(n_layers)) {
      multilinear[, , ell] <-
        tcrossprod(sweep(U, 2, W[ell, ], "*"), V)
      additive[[ell]] <-
        fit_additive_layer(Y[, , ell] - multilinear[, , ell], maxit = 5)
      residual[, , ell] <- Y[, , ell] - additive[[ell]]$fitted
    }
    current_sse <- sum((residual - multilinear)^2, na.rm = TRUE)
    if (is.finite(previous_sse) &&
        abs(previous_sse - current_sse) / (1 + previous_sse) < tolerance) {
      break
    }
    previous_sse <- current_sse
  }

  component_order <- base::order(colSums(abs(W)), decreasing = TRUE)
  U <- U[, component_order, drop = FALSE]
  V <- V[, component_order, drop = FALSE]
  W <- W[, component_order, drop = FALSE]
  for (component in seq_len(rank)) {
    anchor <- which.max(abs(W[, component]))
    if (W[anchor, component] < 0) {
      U[, component] <- -U[, component]
      W[, component] <- -W[, component]
    }
  }
  rownames(U) <- rownames(Y)
  rownames(V) <- colnames(Y)
  rownames(W) <- dimnames(Y)[[3]]
  colnames(U) <- colnames(V) <- colnames(W) <-
    paste("Component", seq_len(rank))
  list(
    U = U, V = V, W = W, additive = additive,
    sse = current_sse,
    explained_after_additive = 1 - current_sse / additive_only_sse,
    iterations = iteration
  )
}

tensor_fit <- cache_fit("icews_2010_tensor_rank2_v2", {
  start_seeds <- c(6886, 6887, 6888, 6889, 6890, 6891, 6892, 6893)
  starts <- lapply(start_seeds, function(start_seed) {
    fit_tensor_start(Y_work, rank = 2, seed = start_seed)
  })
  start_sse <- vapply(starts, `[[`, numeric(1), "sse")
  best <- starts[[which.min(start_sse)]]
  best$start_sse <- setNames(start_sse, start_seeds)
  best
})

tensor_loadings <- sweep(
  tensor_fit$W, 2, apply(abs(tensor_fit$W), 2, max), "/"
)
rownames(tensor_loadings) <- layer_labels[rownames(tensor_loadings)]
round(tensor_loadings, 2)

#>                      Component 1 Component 2
#> Verbal cooperation          0.55        0.86
#> Material cooperation        0.97        0.86
#> Verbal conflict             1.00        0.98
#> Material conflict           0.94        1.00

round(tensor_fit$start_sse, 1)

#>    6886    6887    6888    6889    6890    6891    6892    6893 
#> 47853.6 47896.8 47853.7 47853.7 47853.4 47853.4 47857.0 47853.3

top_component_dyads <- do.call(rbind, lapply(seq_len(ncol(tensor_fit$U)),
  function(component) {
    surface <-
      tensor_fit$W["matlConf", component] *
      outer(tensor_fit$U[, component], tensor_fit$V[, component])
    diag(surface) <- NA
    top <- base::order(surface, decreasing = TRUE, na.last = NA)[1:8]
    cells <- arrayInd(top, dim(surface))
    data.frame(
      component = paste("Component", component),
      source = rownames(tensor_fit$U)[cells[, 1]],
      target = rownames(tensor_fit$V)[cells[, 2]],
      fitted_component = round(surface[top], 2)
    )
  }))
top_component_dyads

component	source	target	fitted_component
Component 1	United States	China	7.60
Component 1	China	United States	7.18
Component 1	United States	Pakistan	6.19
Component 1	United States	Iran, Islamic Republic Of	6.18
Component 1	Japan	United States	5.87
Component 1	United States	Afghanistan	5.77
Component 1	United States	Japan	5.39
Component 1	United States	India	5.23
Component 2	Iran, Islamic Republic Of	Israel	4.03
Component 2	Israel	Lebanon	3.26
Component 2	Israel	Iran, Islamic Republic Of	3.17
Component 2	Lebanon	Israel	3.06
Component 2	Iran, Islamic Republic Of	Lebanon	2.91
Component 2	Egypt	Israel	2.80
Component 2	Israel	Iraq	2.70
Component 2	Iraq	Israel	2.64

loading_plot_data <- data.frame(
  layer = rep(rownames(tensor_loadings), ncol(tensor_loadings)),
  component = rep(colnames(tensor_loadings),
                  each = nrow(tensor_loadings)),
  loading = as.vector(tensor_loadings)
)
ggplot(loading_plot_data,
       aes(layer, loading, fill = component)) +
  geom_hline(yintercept = 0, color = "grey70", linewidth = 0.4) +
  geom_col(position = "dodge", width = 0.72) +
  scale_fill_manual(values = c("#18453B", "#7BBD00")) +
  coord_cartesian(ylim = c(-1.05, 1.05)) +
  labs(x = NULL, y = "Normalized layer loading", fill = NULL) +
  theme(axis.text.x = element_text(angle = 20, hjust = 1),
        legend.position = "bottom")

A grouped bar chart of two normalized component loadings across verbal cooperation, material cooperation, verbal conflict, and material conflict. Both components have positive loadings in every layer. — Layer loadings from the rank-2 joint approximation, normalized so the largest absolute loading within each component equals one. Both shared relational profiles appear in all four layers rather than separating into a pure cooperation axis and a pure conflict axis.

The purpose of this fit is to answer a direct question: after removing which states are generally active as sources and targets within each relation, do the same state pairs remain unusually active across several relations? A component is one recurring source-target pattern. The layer loading says how strongly that same pair pattern appears in each relation. The code compares eight starting values and retains the fit with the smallest squared-error objective.

The answer is yes, although the fit is modest. The two-component surface accounts for about 23.3% of the transformed variation left after layer-specific sender and receiver effects. Component 1 concentrates on pairs involving the United States, including the United States with China, Pakistan, Iran, Afghanistan, and Japan. Component 2 concentrates on a Middle East set that includes Iran with Israel, Israel with Lebanon and Iran, and Egypt or Iraq with Israel.

Both pair patterns appear in verbal cooperation, material cooperation, verbal conflict, and material conflict. The model therefore does not find one cooperation component and one conflict component. It finds state pairs that remain unusually active across several forms of interaction after broad actor activity is removed. The loading table supports that conclusion because every relation expresses both patterns. The values are normalized within each component: 1.00 marks the strongest layer for that component, not a probability, and values should be compared down a column rather than across the two columns.

That result is in-sample and exploratory. Shared profiles can reflect common geopolitical engagement, common exposure to news coverage, or measurement carried across event categories. The next step would be a count model with covariates, rank sensitivity, held-out dyad-year prediction, and replicated layer-specific degree and overlap statistics. Flattening would hide the relation-specific loadings. Four separate fits would hide the shared profiles. The tensor formulation lets us ask about both.

12.8 Write the Model Choice as an Argument

Use the same six parts for any application. Question: name the quantity you want to learn. Data: state the actors, ties or outcomes, eligibility rule, direction, and time structure. Primary model: explain where it represents the dependence that matters. Runner-up: name the strongest alternative representation. Flip condition: identify a specific result or diagnostic that would make you change models. Claim: state what the fitted result would and would not establish.

12.9 Published Example: With Whom Do Syrian Rebels Cooperate?

Gade et al. (2019) ask why armed organizations chose some partners rather than others during the Syrian civil war. The paper is the published study behind today’s application. It uses claims of joint military operations from July 2012 through June 2015 and compares three explanations: organizations may cooperate with ideologically similar groups, groups of similar power, or groups backed by the same state. Shared location is important because joint operations require organizations to be able to operate in the same place.

The published outcome counts joint operations, while our classroom comparison simplifies the relationship to a binary cooperation network. In the paper’s main regression, AME represents broad organization-level activity and higher-order dependence while estimating associations with ideology, power, sponsorship, and location. The paper then uses activity-constrained simulations based on organizations’ observed joint-operation counts as a simpler check: do connected organizations look more ideologically similar than we would expect from their very different activity levels alone?

The published and classroom coefficient findings therefore do not have to match. The paper asks which characteristics are associated with how often pairs conducted joint operations. The capstone asks which characteristics are associated with whether any joint operation was recorded and uses a different specification. Changing the outcome and conditioning set changes the association being estimated.

Using the course framework, the model-choice argument can be reconstructed as follows:

Part of the argument	Application to Gade et al. (2019)
Question	Which differences between armed organizations are associated with more joint military operations?
Data	A valued, undirected cooperation network built from claims of joint operations, plus organization and dyad measures of ideology, power, sponsorship, and location.
Primary model	AME regression, which represents broad organization activity and higher-order dependence while estimating the measured associations.
Strongest model runner-up	An ERGM with the same dyadic covariates, organization heterogeneity, and any theoretically named configurations. The paper does not fit this model, but it would be the clearest alternative if the argument centered on an explicit graph configuration.
Published robustness check	Activity-constrained simulations based on observed joint-operation counts, which ask whether ideological similarity among partners exceeds what activity alone would generate.
Flip condition	Reconsider the ideology conclusion if it disappears after shared location enters or if the activity-constrained simulations reproduce the observed similarity without using ideology.
Claim	Ideologically closer organizations cooperated more often in these data. The fitted pattern does not by itself establish that ideology caused the partnerships.

12.9.1 What the Published Comparison Added

Ideological proximity had the clearest and most consistent association with cooperation across the paper’s AME specifications. Shared location also mattered. Evidence for power similarity was weaker and less consistent, while shared state sponsorship did not receive clear support. The activity-constrained simulations likewise found more ideological similarity among partners than differences in activity alone would generate.

This is a useful model-choice argument because the fitted AME and the robustness diagnostic ask connected but distinct questions. AME estimates the ideology association while representing several forms of network dependence. The simulations ask whether the same broad pattern survives a simpler benchmark tied to each organization’s observed activity. Agreement across them strengthens the descriptive conclusion, but the simulation is not a runner-up estimator and neither procedure identifies the behavioral mechanism that produced the partnerships.

The substantive conclusion is therefore straightforward: in this fragmented conflict, organizations with more similar ideologies cooperated more often. Power and sponsorship did less to distinguish partnerships in these analyses. The model and robustness check narrow the set of plausible descriptions, but they do not show that changing an organization’s ideology would change its partners.

13 The Frontier: Scale and the Adjacent Field

What breaks when \(n\) gets big. There is no universal cutoff. Runtime depends on density, the number of layers or periods, the model terms, the approximation, the implementation, and the hardware. The useful question is which object or operation grows fastest in the analysis you actually plan to run:

Family	Common bottleneck	A defensible scaling move	What that move gives up
Descriptives and community detection	Edge storage and repeated graph traversals	Sparse matrices and local algorithms	Some global summaries or exact optimization
Stochastic blockmodels	Searching over memberships and storing dense probability objects	Sparse implementations, variational approximations, or restricted block structure	Exact posterior exploration or some dependence
SRM, AME, and latent-space models	\(n\times n\) outcomes, latent coordinates, and repeated likelihood or posterior updates	Case-control likelihoods that keep ties and sample non-ties (Raftery et al. 2012)	An approximate likelihood whose quality depends on the sampling design
ERGMs	Simulation mixing, change-statistic cost, and degeneracy risk	Local-dependence or block-structured specifications	Dependence across the imposed blocks

A computational shortcut is also a statistical approximation. Before using one, state whether it changes the estimand, the dependence that can be represented, or only the numerical route to the same target.

Use the procedures from this course as anchors. lame’s ALS option and the tensor ALS fit use block-coordinate updates to find a fast point estimate of one nonconvex loss. They do not provide a posterior distribution, and multiple starts are needed. A variational blockmodel optimizes a lower bound rather than exploring the exact posterior. ERGM maximum pseudolikelihood changes the fitting target by multiplying conditional tie contributions; MCMC maximum likelihood instead simulates networks to approximate the graph-likelihood score. btergm combines MPLE with a temporal bootstrap, so its uncertainty depends on the observed network waves available for resampling. Node2vec uses sampled walks and negative examples to optimize a context-prediction objective rather than a tie likelihood.

The applied choice follows the deliverable. If the goal is fast exploratory structure, a point estimator plus held-out checks may be enough. If the claim concerns uncertainty in a coefficient or latent feature, use a procedure that propagates that uncertainty and verify its convergence. If the theory depends on global network configurations, a local or pseudolikelihood shortcut may omit exactly the dependence of interest. Always report what was optimized, what was approximated or sampled, and which diagnostic checks that approximation.

The adjacent field. Node embeddings and graph neural networks often prioritize link prediction or node classification and can exploit much larger sparse graphs than the fitted models used here. That does not mean they inherently lack uncertainty or interpretation, and it does not mean a classical model is automatically identified. It means their usual objective and reporting workflow are different. If prediction is the deliverable, evaluate held-out predictions against a serious baseline. If a parameter or causal contrast is the deliverable, state the assumptions that identify it and use an inferential procedure built for that target. In either case, treat coordinate axes and learned features as model-dependent objects rather than ready-made substantive dimensions.

14 The Frontier: The Econometrics of Networks

The econometrics of networks includes several distinct problems: estimating peer effects when outcomes respond through a network, modeling strategic network formation, and estimating dyadic outcome models with high-dimensional actor fixed effects. In the formation lane most relevant here, sender and receiver fixed effects can represent degree heterogeneity. Depending on the model, researchers condition them out, estimate and bias-correct them, or absorb them in two-way gravity-style specifications estimated by PPML. Those choices concern the conditional mean and identification. Dyadic or actor-clustered standard errors concern uncertainty around an estimator. They are complements, not substitutes for one another or for latent structure.

Where to start (all on the shelf below): de Paula (2017) is the doorway survey; Graham (2017, Econometrica) develops a fixed-effects approach for binary network formation; de Paula (2020, Annual Review of Economics) provides a shorter updated survey. The useful connection to this course is not a disciplinary label. It is the recurring need to decide whether actor heterogeneity belongs in fixed effects, random effects, latent structure, or the uncertainty calculation, and what assumptions each choice requires.

15 Close: The Map, Filled In

Day 8 began with a data object and a measurement decision. Each day since then added a tool for a different job. This is the complete toolbox in question-first form:

job	start with	day	first check
Define the network and risk set	netify and descriptives	8	direction, roster, eligibility, time, and zero versus missing
Represent broad actor heterogeneity	SRM	9	actor-level GOF, exogeneity, and enough actors
Learn discrete relational roles	SBM or NetMix	9	K, starts, role interpretation, and simulated GOF
Learn continuous residual profiles	latent distance, LFM, or AME	10	rank, starts or chains, held-out fit, and simulation GOF
Parameterize named graph configurations	ERGM	11	MCMC diagnostics and simulated GOF beyond fitted statistics
Model change across observed waves	TERGM, SAOM, or dynamic AME	10 through 12	risk set, number of transitions, and temporal GOF
Calibrate uncertainty for repeated actors	DCR	14	identical coefficients and term-specific interval changes
Define and identify spillover contrasts	exposure mappings or negative controls	13	support, exposure mapping, assignment probabilities, or proxy assumptions
Build scalable predictive node features	node2vec or related embeddings	14	held-out performance, simple baselines, and tuning stability
Represent several relation layers jointly	multilinear tensor or native multilayer model	14	layer definition, rank sensitivity, and held-out layers or dyads
Estimate directed temporal influence	SIR	14	lag definition, starts, uncertainty, and out-of-sample prediction

Three seams in this field remain open research questions, and you now know enough to read them as active arguments rather than settled doctrine.

ERGM versus latent variables. Day 11 versus Days 9–10. Both families can pass GOF on different networks, and there is no theorem that picks a family before the question and data are specified. The practical disagreement is whether named configurations should receive explicit conditional parameters or whether broader residual structure should be represented with latent effects. Cranmer and Desmarais (2011) and Minhas, Hoff, and Ward (2019) provide two entry points. I am an author in this debate, so treat my preferences as disclosed rather than decisive.
Modeling dependence versus adjusting the covariance. Day 10 versus today’s DCR section. AME changes the conditional mean by representing actor and pair-specific structure. DCR leaves that mean unchanged and adjusts uncertainty for dyads sharing actors. An open applied question is when those two choices lead to materially different conclusions under data-generating processes where the truth is known. A useful comparison would evaluate bias, interval coverage, and prediction rather than treating either strategy as a universal answer.
Node embeddings versus model-based latent factors. Day 10 versus today’s node2vec section. Both use low-dimensional dot products, but they reconstruct different objects and support different claims. The open practical question is how much a prediction-first representation and a dyadic outcome model agree when they are evaluated on the same held-out task, without pretending that predictive performance supplies coefficient interpretation or that a likelihood automatically supplies good prediction.

What to read next, per family: one entry point each, chosen so the next step is never “the literature”:

Family	Read next
SRM / AME	Hoff (2021), “Additive and multiplicative effects network models,” Statistical Science 36(1); the family’s own retrospective
Latent space (distance tradition)	Handcock, Raftery & Tantrum (2007), “Model-based clustering for social networks,” JRSS-A 170(2); the latentnet line Day 10 used as its foil; Krivitsky, Handcock, Raftery & Hoff (2009, Social Networks 31(3)) for the random-effects repair
Blockmodels	Peixoto (2019), “Bayesian stochastic blockmodeling,” in Advances in Network Clustering and Blockmodeling; the modern inference view
ERGM	Lusher, Koskinen & Robins (2013), Exponential Random Graph Models for Social Networks; the standard book-length treatment
SAOM	Snijders (2017), “Stochastic actor-oriented models for network dynamics,” Annual Review of Statistics and Its Application 4
TERGM	Taught on Day 11; the ancillary deep dive (`ancillary/tergm.qmd`) is the longer treatment, then Leifeld, Cranmer & Desmarais (2018, JSS 83(6))
Causal	Aronow & Samii (2017, AOAS) for what randomization buys; Ogburn et al. (2024, JASA) for the observational frontier
Multilayer networks	Hoff (2015, Annals of Applied Statistics) for multilinear tensor regression; Kivelä et al. (2014, J. Complex Networks) for the multilayer map
ML / embeddings	Perozzi, Al-Rfou & Skiena (2014) for DeepWalk; Grover & Leskovec (2016) for node2vec; Qiu et al. (2018) for the matrix-factorization connection; Hamilton, Ying & Leskovec (2017) for the broader map
Influence (SIR)	Minhas & Hoff (2026), “Decomposing network influence: social influence regression,” Political Analysis 34(2); the package’s `icews` help page documents its reduced illustrative ICEWS data
Econometrics of networks	de Paula (2017), “Econometrics of network models,” ch. 8 of Advances in Economics and Econometrics: Eleventh World Congress; the doorway; Graham (2017, Econometrica 85(4)) to see the FE approach executed end to end

What Today Deliberately Did Not Show, and Where It Lives

Each family’s native one-call GOF: lame’s gof_plot(fit, statistics = ...), ergm’s gof() + plot(), and RSiena’s sienaGOF(). Today ran everything through the shared calc_gof() yardstick so the lenses were comparable; the native versions are demonstrated on Days 9 through 12 and are what you would use inside a single family.
DCR’s coeftest() one-liner for a publication-style table from the robust variance; today’s tab_dcr build is the transparent version of the same arithmetic.
A fitted SAOM or TERGM on today’s data: this single snapshot cannot estimate change between waves. Your HW 4 sienaFit objects plug into row_saom() in _compare_helpers.R, with the estimand warning stated in that file.
The rest of the SIR surface: residual maps built from augment()’s per-dyad frame, alternative uncertainty estimators, cross-validation, and forecasts. The toolbox aside shows the code path, and the article’s replication archive contains the full ICEWS analysis.

15.1 After the Models Run: Write the Substantive Paragraph

The capstone is not finished when the comparison table is full. It is finished when a reader can understand why the estimates change the substantive argument. A strong final paragraph has five parts:

Question: State the quantity or mechanism the analysis is trying to learn about.
Model choice: Explain why the primary model represents the dependence, time, and actor roles that matter for that question.
Substantive result: Describe the meaningful pattern, actor comparison, predicted probability, or mechanism comparison. Do not march through every coefficient.
Diagnostic: Name the fit, sensitivity, convergence, or design evidence that makes the interpretation credible.
Limit and flip condition: State what the model does not establish and what empirical result would move you to the runner-up model or change the conclusion.

What the capstone comparison shows

“Among 30 armed organizations in Syria, ANF and ASIM cooperate with far more partners than the rest. A simple two-block model isolates those two as a small core but compresses the other 28 organizations into one heterogeneous periphery. In an initial ERGM, ties between organizations with shared partners receive a clear positive conditional association. After every organization receives its own baseline propensity to cooperate, the shared-partner interval includes zero and simulated networks cover the two selected scalar yardsticks. We should not compare 1.38 and 0.18 as though they were the same effect size under the same conditioning set. This saturated diagnostic suggests that the data do not clearly separate an additional shared-partner pattern from the hub structure. It does not show that closure never occurs, and it does not identify why ANF or ASIM became hubs.”

That paragraph is more useful than declaring one model the winner. It connects an estimand to an explanation, makes clear which patterns carry the conclusion, and gives the reader a concrete way to disagree. Advanced methods earn their place when they sharpen the substantive claim and its limits, not when they merely make the methods section longer.

16 Your Turn

Exercises

1. (Required: use the framework.) Choose a dataset or research question you care about. If you do not have one available, use today’s Syrian cooperation network. Walk it through Section 9: which branch matches the question? Write one paragraph justifying a primary model, naming the runner-up, and identifying the observation or diagnostic that would make you switch.

2. (Stretch: Two rivals, reconciled.) Actually fit your chosen model and the runner-up on your data (or on the rebel network: reproduce Section 7’s erg_soc from a cold cache and set it against the AME from Section 6). Extract one comparable claim from each: a coefficient verdict, a GOF band, a recovered structure, and write three sentences reconciling them in the comparison-table style: what question did each answer, and do the answers actually conflict once you say the questions out loud?

3. (Your data, forward-looking.) Which frontier method from today is closest to a dataset you expect to have within two years, and what is the first sentence of its entry-point reading? (Yes, this requires opening the reading. That is the point.)

Guidance: Open after you’ve tried them

## 1. ---------------------------------------------------------------
# There is no answer key for your data, but here is the paragraph
# pattern, worked for today's Syrian cooperation network:
#   "My question is whether shared location is associated with any
#    recorded cooperation after accounting for persistent differences
#    in organizations' overall activity. The network is binary,
#    symmetric, and one snapshot, so the SRM is my primary model.
#    Runner-up: an ERGM with sociality and shared-partner terms, because
#    an explicit closure parameter matters if the argument is about
#    cooperation spreading through common partners. The diagnostic
#    that would flip me is an SRM simulation that reproduces actor
#    heterogeneity but consistently underproduces the observed
#    shared-partner distribution."
# Grade yourself on: (a) Named branch, (b) Named runner-up,
# (c) A *specific, runnable* flip condition, not "if it fits badly."

## 2. ---------------------------------------------------------------
# On the Syrian network the final specification is already on screen:
# its GWESP estimate is about 0.18 with an interval that includes zero,
# and its simulations cover both yardstick statistics. Read that
# coefficient within the final model rather than comparing its size
# directly with the coefficient from another specification. Reproduce
# the fit from a cold cache and notice two things
# while you do:
#   (a) The spec drops the nodal covariates on purpose. sociality()
#       absorbs every nodal main effect, so keeping nodecov() terms
#       makes the model nonidentifiable; ERGM will warn about
#       singularities, and the warning is correct, not a bug.
#   (b) The three comparison-style sentences you are after: The AME
#       reproduced the observed transitivity without an explicit
#       closure parameter. In the ERGM with sociality, the residual
#       GWESP estimate is small and imprecise. The two families looked
#       more different before both represented actor heterogeneity,
#       but neither result proves the generating mechanism or the
#       absence of closure.

## 3. ---------------------------------------------------------------
# No code is required. Name the outcome, actors, data structure, and
# prediction or inference task the frontier method would need.

17 Reading

The day’s working citations, in order of appearance:

Gade, Gabbay, Hafez & Kelly (2019), “Networks of cooperation: rebel alliances in fragmented civil wars,” Journal of Conflict Resolution: the data every lens pointed at today.
Carlson, Incerti & Aronow (2024), “Dyadic clustering in international relations,” Political Analysis: the shared-actor uncertainty correction introduced in lens 0.
Leifeld & Cranmer (2019), Network Science; Block, Hollway, Stadtfeld, Koskinen & Snijders (2022) and Leifeld & Cranmer (2022), also Network Science: the TERGM–SAOM exchange that Day 12 introduced; the ancillary TERGM deep dive walks the papers piece by piece.
Minhas & Hoff (2026), “Decomposing network influence: social influence regression,” Political Analysis: the influence estimand and the covariate-parameterized bilinear structure that replaces an unrestricted influence matrix with a smaller set of estimable channel weights; sir is its software.
Kivelä, Arenas, Barthelemy, Gleeson, Moreno & Porter (2014), “Multilayer networks,” Journal of Complex Networks.
Hoff (2015), “Multilinear tensor regression for longitudinal relational data,” Annals of Applied Statistics: a joint low-rank model for multiway relational arrays.
Raftery, Niu, Hoff & Yeung (2012), “Fast inference for the latent space network model using a case-control approximate likelihood,” Journal of Computational and Graphical Statistics.
Perozzi, Al-Rfou & Skiena (2014), “DeepWalk: online learning of social representations,” KDD 2014: an early method for creating actor features from short trips through a graph.
Tang, Qu, Wang, Zhang, Yan & Mei (2015), “LINE: large-scale information network embedding,” WWW 2015: direct ties and similarity in neighbor patterns at very large scale.
Grover & Leskovec (2016), “Node2vec: scalable feature learning for networks,” KDD 2016: biased second-order walks and the node2vec training objective.
Qiu et al. (2018), “Network embedding as matrix factorization: unifying DeepWalk, LINE, PTE, and node2vec,” WSDM 2018: the careful matrix-factorization connection.
Hamilton, Ying & Leskovec (2017), “Inductive representation learning on large graphs,” NeurIPS 2017: GraphSAGE learns how to combine actor attributes with sampled neighbor information, including for actors not observed during training.
Hamilton, Ying & Leskovec (2017), “Representation learning on graphs: methods and applications,” IEEE Data Engineering Bulletin: the broader representation-learning map.
Brown et al. (2021), “Network embedding methods for large networks in political science”: an applied comparison using network links to predict ideological labels for Twitter users and internet domains.
de Paula (2017), “Econometrics of network models,” in Advances in Economics and Econometrics: Eleventh World Congress: the doorway into the econometrics-of-networks lane; de Paula (2020), “Econometric models of network formation,” Annual Review of Economics is the shorter, more current second pass.
Graham (2017), “An econometric model of network formation with degree heterogeneity,” Econometrica: fixed effects conditioned out of a binary network via tetrad logit; read it next to Day 10 and feel the two traditions pull.
Plus the per-family “read next” table in Section 15, which is the reading list that matters after this week.

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] sir_1.1.0         ggplot2_4.0.3     ergm_4.12.0       network_1.20.0   
#> [5] blockmodels_1.1.5 lame_1.3.5        netify_1.5.3     
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.1         viridisLite_0.4.3        dplyr_1.2.1             
#>  [4] farver_2.1.2             loo_2.9.0                viridis_0.6.5           
#>  [7] S7_0.2.2                 ggraph_2.2.2             fastmap_1.2.0           
#> [10] tensorA_0.36.2.1         tweenr_2.0.3             digest_0.6.39           
#> [13] lifecycle_1.0.5          rle_0.10.0               magrittr_2.0.5          
#> [16] posterior_1.6.1          compiler_4.3.3           rlang_1.2.0             
#> [19] tools_4.3.3              igraph_2.2.2             yaml_2.3.12             
#> [22] knitr_1.51               labeling_0.4.3           graphlayouts_1.2.2      
#> [25] htmlwidgets_1.6.4        multiwayvcov_1.2.3       RColorBrewer_1.1-3      
#> [28] abind_1.4-8              withr_3.0.2              purrr_1.2.2             
#> [31] grid_4.3.3               polyclip_1.10-7          scales_1.4.0            
#> [34] MASS_7.3-60.0.1          cli_3.6.6                rmarkdown_2.30          
#> [37] generics_0.1.4           otel_0.2.0               robustbase_0.99-7       
#> [40] biglm_0.9-3              DBI_1.3.0                cachem_1.1.0            
#> [43] ggforce_0.5.0            stringr_1.6.0            parallel_4.3.3          
#> [46] matrixStats_1.5.0        vctrs_0.7.3              boot_1.3-30             
#> [49] Matrix_1.6-5             sandwich_3.1-1           jsonlite_2.0.0          
#> [52] patchwork_1.3.2          ggrepel_0.9.6            trust_0.1-9             
#> [55] speedglm_0.3-5           ggnewscale_0.5.2         tidyr_1.3.2             
#> [58] glue_1.8.1               statnet.common_4.13.0    dcr_0.0.0.9000          
#> [61] DEoptimR_1.1-4           distributional_0.6.0     stringi_1.8.7           
#> [64] gtable_0.3.6             tibble_3.3.1             pillar_1.11.1           
#> [67] htmltools_0.5.9          R6_2.6.1                 Rdpack_2.6.4            
#> [70] tidygraph_1.3.1          lpSolveAPI_5.5.2.0-17.15 evaluate_1.0.5          
#> [73] lattice_0.22-5           rbibutils_2.4            backports_1.5.0         
#> [76] memoise_2.0.1            broom_1.0.11             Rcpp_1.1.1-1.1          
#> [79] coda_0.19-4.1            gridExtra_2.3            checkmate_2.3.4         
#> [82] xfun_0.55                zoo_1.8-15               pkgconfig_2.0.3