Day 10: Latent Geometry, From Distance to the Factor Model

Advanced Network Analysis · ICPSR

Author

Shahryar Minhas

Published

July 30, 2026

NoteHow to Use This Document

Open the Day 10 teaching deck.

Use the deck as the route through the class. Each link in the deck opens the exact walkthrough section we need, so we can run or inspect that part and return to the deck without searching through the file. After restarting R, run the setup chunk first. Every example is seeded at 6886, and the downloaded folder includes the data and saved model fits used to render this page.

The expensive model calls use cache_fit(), which stores their results in cache/. If the model code matches the saved version, the chunk loads the result quickly. If you change the call, the cache detects the change and refits. Delete an individual cache file only when you deliberately want to estimate that model again.

Most packages are on CRAN. Until lame reaches CRAN, install its current release from GitHub:

install.packages(c("latentnet", "sna", "netify", "ggplot2", "reshape2",
                   "patchwork", "ggrepel", "plyr", "gridExtra", "coda",
                   "digest", "posterior", "remotes"))
remotes::install_github(
  "netify-dev/lame",
  dependencies = TRUE,
  build_vignettes = FALSE,
  upgrade = "never"
)

Day 9 represented recurring actor patterns with the SRM and discrete relational roles with blockmodels. Day 10 keeps the same concern about connected observations but moves to continuous latent geometry. We begin with a distance map because it is intuitive, then use a factor model when similarity and complementary roles both matter.

lame is in the CRAN submission process and should be available there shortly. At the time this walkthrough was rendered, it was not yet listed on CRAN. Once it appears, the simplest installation will be install.packages("lame"). Until then, use one of the two routes below.

0.0.1 Option 1: Install From GitHub Source

This is the easiest route if your computer already has C++ compilation tools. Install the lighter remotes helper, then install lame from its current GitHub repository:

install.packages("remotes")
remotes::install_github(
  "netify-dev/lame",
  dependencies = TRUE,
  build_vignettes = FALSE,
  upgrade = "never"
)

The required compilation tools are Rtools on Windows, the Xcode command-line tools on macOS, and build-essential on Ubuntu or WSL. Use the version of Rtools that matches your version of R. On macOS, run xcode-select --install in Terminal. On Ubuntu or WSL, run sudo apt update followed by sudo apt install build-essential.

0.0.2 Option 2: Use a Compiled GitHub Release

If you do not have compilation tools, open the lame releases page, expand the assets for the latest release, and download the compiled file that matches your system. The current release provides Windows, Intel macOS, Apple silicon macOS, and x86_64 Linux builds.

Install the package dependencies from CRAN first:

install.packages(c(
  "Rcpp", "RcppArmadillo", "ggplot2", "ggrepel", "ggforce",
  "gridExtra", "coda", "patchwork", "cli", "abind", "netify"
))

For Windows, download the file ending in -windows.zip, then run:

install.packages(file.choose(), repos = NULL, type = "win.binary")

For an Apple silicon Mac, download the file ending in -macos-arm64.tgz. For an Intel Mac, download the file ending in -macos-intel.tgz. Then run:

install.packages(file.choose(), repos = NULL, type = "mac.binary")

For x86_64 Linux or WSL, choose the long filename containing _R_x86_64-pc-linux-gnu-linux.tar.gz, not the shorter source archive. Then run:

install.packages(file.choose(), repos = NULL, type = "source")

The Linux binary is the least portable option because it must match the system architecture and a compatible R setup. On Linux or WSL, installing build-essential and using the GitHub source route is usually more reliable.

After either route, verify the installation:

library(lame)
packageVersion("lame")
stopifnot(packageVersion("lame") >= "1.3.5")

1 Where Geometry Enters

Day 9 showed that some organizations work with many partners and that cooperation can be organized around a small tactical core. Day 10 asks what those summaries still miss. Can we represent which particular organizations work together, even when they have similar overall levels of activity? In the longitudinal application, can we distinguish a state becoming involved across many relationships from one particular directed state relationship changing?

The SRM asks whether some actors repeatedly appear in many relationships and whether the two directions of a pair remain connected. It does not directly say why one particular pairing is common while another pairing involving the same actor is rare. That distinction matters whenever broad involvement and bilateral relationships can change separately.

The blockmodel handles this by placing actors with similar partner lists into a few groups. Today replaces those boxes with continuous profiles, allowing two actors to resemble one another partly rather than forcing them into exactly the same role.

The day has three parts:

  1. The latent distance model gives every actor a position and makes ties more likely for nearby pairs.
  2. A simulated role network shows where a small Euclidean distance map struggles.
  3. The latent factor model, and its full AME form, represents both similarity and complementary roles through a continuous relational surface.

The main question is not whether a map looks appealing. It is whether the model helps us say something clearer about who relates to whom, what changed, and what the data still cannot explain.

1.1 What Might the Missing Pattern Represent?

Let’s begin with the relationships we will study, not the machinery. Among Syrian armed organizations, was tactical cooperation scattered across unrelated pairs or concentrated among organizations with overlapping sets of partners? In the ICEWS panel, did states simply become more involved in high-volume conflictual relationships overall, or did particular directed relationships change in distinctive ways? These are substantive questions about joint operations and coded state actions. The latent models provide different ways to summarize the partner patterns behind them.

A latent model asks whether a lower-dimensional relationship pattern helps account for what remains. Several possibilities should stay separate:

  • Proximity: armed organizations with similar unmeasured cooperation positions may be more likely to work together.
  • Complementary roles: two actors may recur together because their relationship profiles line up in a compatible way.
  • Adjustment: the latent term may help the model handle connected observations while we focus on a measured covariate.
  • Measurement: the latent relationship pattern itself may be the object we want to study.

Those uses lead to different interpretations. A useful predictive map is not automatically a causal adjustment, and a latent axis is not automatically ideology, power, or any other concept we recognize after seeing the labels.

1.2 The Applied Model Compass

Model Question We Want to Understand How the Model Helps What Makes It Distinct What We Would Say About the Result
Additive SRM Are some actors involved across many relationships, regardless of the particular partner? Represents recurring actor activity and, for directed networks, reciprocal dependence It cannot give one particular pairing a distinctive history “Syria can begin appearing in many more above-threshold relationships, but this alone does not tell us which directed relationships changed.”
Latent distance model Is tactical cooperation concentrated among organizations with overlapping sets of partners? Infers a map in which nearby actors have higher fitted tie probabilities It assumes one common notion of proximity “Organizations with similar overall cooperation patterns are placed near one another, but left and right on the map have no automatic meaning.”
Latent factor or eigenmodel Can actors play similar or complementary parts in the network without belonging to the same discrete group? Gives actors continuous relationship profiles and combines them pair by pair It can represent both similar partner lists and partnerships between unlike roles “Two states can repeatedly engage the same kinds of counterparts even when they do not engage one another.”
AME Is change spread across all of a state’s relationships, or concentrated in particular directed pairs? Combines broad state involvement with a separate adjustment for each recurring source-target pattern It keeps the broad and bilateral parts visible in one model “Iran toward Syria and Syria toward Iran can follow different paths even after accounting for how often each state appears with other partners.”

The model should follow the question. Use distance when proximity is a defensible picture of tie formation. Use factors when complementary roles matter. Use the full AME when we need observed covariates, broad actor tendencies, reciprocity, and pair-specific relationship profiles in the same model.

2 The Distance Idea: A Map You Can Predict From

Let’s start with the cooperation question. Was tactical cooperation scattered across unrelated pairs, or did organizations with overlapping sets of partners form a recognizable cluster? The distance model turns each organization’s complete partner list into a position on a map. Pairs placed closer together receive a higher fitted chance of a recorded joint operation.

The everyday analogy is a seating chart. If I tell you that nearby people are more likely to talk, you can predict many conversations from the arrangement. A latent distance model runs that reasoning backward: it uses the observed ties to estimate positions that make nearby pairs more likely to connect.

Let’s make the intuition concrete by running it forwards: put points in a space ourselves, and let the space generate a network. Keep this simulated network in mind; it comes back in Section 6 as one of the two test cases.

set.seed(6886)
np  <- 40
grp <- rep(c("group A", "group B"), each = np / 2)

# Every node gets a hidden 2-D position. Group A sits left, B sits right,
# but they overlap; this is not two clean clusters.
Z_true <- cbind(
  ifelse(grp == "group A", -1.1, 1.1) + rnorm(np, 0, 0.6),
  rnorm(np, 0, 0.8)
)

# Closer in the space -> higher tie probability. That is the whole model:
# a baseline minus the distance between the two points.
D  <- as.matrix(dist(Z_true))
P  <- plogis(2 - 1.6 * D)
# plogis() = inverse-logit: turns the log-odds into a tie probability in (0,1)
Y  <- (matrix(runif(np^2), np, np) < P) * 1  # Flip a weighted coin for each pair: tie with probability P
Y[lower.tri(Y, diag = TRUE)] <- 0          # Undirected: keep upper triangle
Y  <- Y + t(Y)

plot(Z_true, col = ifelse(grp == "group A", "#18453B", "#7BBD00"),
     pch = 19, xlab = "latent dim 1", ylab = "latent dim 2", asp = 1)
for (i in 1:np) for (j in 1:np) if (Y[i, j] == 1)
  segments(Z_true[i,1], Z_true[i,2], Z_true[j,1], Z_true[j,2],
           col = adjustcolor("grey40", 0.35))

A network generated from latent positions. Two loose groups on the left/right; points that sit close tie, points that sit far do not. The model we fit later runs this picture in reverse: ties in, positions out.

The picture is exactly the intuition: dense wiring inside each blob, thin wiring across the gap. We can put a number on it: tie density among pairs that are close in the latent space versus pairs that are far:

close <- D[upper.tri(D)] < median(D[upper.tri(D)])
ties  <- Y[upper.tri(Y)]
checkpoint(
  density_when_close = mean(ties[close]),
  density_when_far   = mean(ties[!close])
)
#> ------------------------------------------------------------------
#> CHECKPOINT: density_when_close = 0.523   |   density_when_far = 0.095
#> ------------------------------------------------------------------

Close pairs tie at 0.52; far pairs at 0.09. Distance in the latent space is doing all the work.

Now notice what this tends to buy us. Distances obey the triangle inequality:

‖ui−uk‖≤‖ui−uj‖+‖uj−uk‖. \lVert u_i-u_k\rVert \leq \lVert u_i-u_j\rVert + \lVert u_j-u_k\rVert.

That does not make closeness literally transitive. Two short legs can add up to a third leg that is longer than either one, and an observed tie is only a probabilistic clue about distance. What the inequality does is constrain the third distance: if ii is very near both jj and kk, then jj and kk cannot be arbitrarily far apart. Because tie probability falls with distance, latent homophily therefore raises the probability of clustered and transitive configurations. It makes triangles common in expectation; it does not force every open two-path to close. That probabilistic version is the structural reason this model exists, and Section 4 shows the trade-off that comes with it.

2.1 The Model, Written Down

The model just says what the picture said. For an undirected binary tie yijy_{ij}:

logitPr(yij=1)=α+β⊤xij⏟covariates−‖ui−uj‖⏟distance betweenlatent positions \text{logit}\, \Pr(y_{ij} = 1) \;=\; \alpha + \underbrace{\beta^\top x_{ij}}_{\text{covariates}} \;-\; \underbrace{\lVert u_i - u_j \rVert}_{\substack{\text{distance between}\\\text{latent positions}}}

Read it slowly. Each actor ii gets a position uiu_i in a dd-dimensional space (often d=2d = 2, so you can draw it). The further apart two actors sit, the more we subtract from the log-odds of a tie. The β⊤xij\beta^\top x_{ij} part is a regression on measured dyadic or nodal covariates; the geometry represents residual structure those covariates do not explain. This is Hoff, Raftery, and Handcock’s (2002) latent space model, and it is what latentnet::ergmm() fits.

2.2 What the Distance Likelihood Rewards

Conditional on the positions and covariates, the basic undirected binary LDM treats the unordered dyads as independent Bernoulli observations:

p(Y∣U,X,α,β)=∏i<jpijyij(1−pij)1−yij,logit(pij)=α+β⊤xij−‖ui−uj‖. p(Y\mid U,X,\alpha,\beta) = \prod_{i<j} p_{ij}^{y_{ij}}(1-p_{ij})^{1-y_{ij}}, \qquad \operatorname{logit}(p_{ij}) = \alpha+\beta^\top x_{ij}-\lVert u_i-u_j\rVert.

Each observed tie rewards arrangements that give the pair a higher fitted probability, while each observed non-tie rewards arrangements that give the pair a lower fitted probability. The algorithm must compromise across the entire network because moving one actor changes its distance from everyone else. One observed tie therefore does not tell us that two actors are close. Their complete partner lists determine where the model can place them.

Three details matter. First, the product is over i<ji<j, so an undirected dyad is counted once. Second, the prior prevents weakly informed positions from drifting without limit. Third, the basic distance term can create some degree variation indirectly because an actor near the middle of the cloud may be close to many others. It does not give every actor a free activity parameter. If residual degree heterogeneity remains, add a random sociality term such as rsociality() for an undirected ergmm() fit, or use AME’s additive ai+bja_i+b_j terms. Geometry and activity are different jobs, even when a fitted map partly performs both.

2.3 How the Distance Model Is Estimated

The model now has a clear substantive representation and a clear probability rule. Estimation runs that story backward: it searches for coefficients and actor positions that give high probabilities to observed ties and low probabilities to observed non-ties across the whole network.

Bayesian estimation combines the likelihood with priors on the coefficients, positions, and variance quantities:

p(α,β,Z,Σ∣Y,X)∝p(Y∣X,α,β,Z)p(α,β)p(Z∣Σ)p(Σ). p(\alpha,\beta,Z,\Sigma\mid Y,X) \propto p(Y\mid X,\alpha,\beta,Z)\, p(\alpha,\beta)\,p(Z\mid\Sigma)\,p(\Sigma).

The posterior is the statistical target. latentnet::ergmm() uses MCMC as the algorithm for exploring it. The sampler does not find one map and stop. It repeatedly updates the coefficients and positions and retains a collection of plausible explanations, so uncertainty in the map and coefficients is carried forward together.

2.4 How MCMC Represents Uncertainty in the Distance Model

Return to the seating-chart analogy. We observe the relationships but not the arrangement or coefficient values. Instead of pretending there is one certain answer, MCMC revisits many combinations that are plausible under the likelihood and priors.

The computer begins with a provisional explanation and then takes turns:

  1. Hold the current actor positions fixed and update the regression coefficients and variance quantities.
  2. Hold those updated pieces fixed and move one or more actors to positions supported by their complete relationship patterns.
  3. Calculate the likelihood and prior contribution for the revised state.
  4. Record the complete state and repeat the cycle.

One recorded state is a posterior draw. A sequence of draws is the chain. The chain is not collecting new actors or networks. It is exploring different explanations of the same observed network. Neighboring draws are correlated because each update begins from the preceding state.

We discard early burn-in draws that may still reflect the arbitrary starting configuration. We inspect trace plots to see whether retained draws move around a stable region. Effective sample size translates the number of correlated saved draws into the approximate amount of independent Monte Carlo information they contain. With multiple chains, R̂\widehat R asks whether different starting points reached the same posterior region.

Coordinate summaries require special care because the same fitted distances can be displayed after rotating or reflecting the map. The collapsed identification box later in the walkthrough explains this issue and the alignment procedure. For the main argument, interpret fitted distances and probabilities before trying to name axes.

Posterior means summarize the retained distribution, and credible intervals summarize uncertainty under the specified model. A longer chain can reduce Monte Carlo error. It cannot add actors, repair weak measurement, or turn a latent association into a causal effect.

NoteWhat the Priors Are Doing

The priors regularize weakly identified positions, their overall spread, and the regression coefficients. They matter most when the network gives limited information about a direction or scale. A prior is part of the model. It is not a convergence device or a substitute for evidence.

Compare this with the SRM from Day 9:

yij≈β⊤xij+ai+bj+γ(ui,uj)⏟latent pattern to estimate+εij y_{ij} \approx \beta^\top x_{ij} + a_i + b_j + \underbrace{\gamma(u_i, u_j)}_{\text{latent pattern to estimate}} + \varepsilon_{ij}

The open term is γ(ui,uj)\gamma(u_i,u_j), the relationship pattern left after measured predictors and broad actor effects. The distance model uses γ(ui,uj)=−‖ui−uj‖\gamma(u_i,u_j) = -\lVert u_i - u_j\rVert. That choice is intuitive and useful, but it also imposes a specific idea: similarity raises tie probability.

ImportantThe One Sentence for This Section

The model uses an actor’s complete partner list to decide where to place it. Actors with similar relationship patterns tend to appear closer together, and nearby pairs receive a higher fitted chance of a tie. This makes clustered relationships more common without claiming that one observed tie caused another.

3 A Latent-Distance Example: Syrian Armed Organizations

Latent-distance maps appear in studies of alliances, legislative collaboration, friendship, trade, and armed-group relationships. We will fit one with latentnet::ergmm() and practice reading what the map does and does not say.

We return to the cooperation data from Gade et al. (2019), covering July 2012 through June 2015. The observed network contains 31 Syrian armed organizations. An undirected tie means the pair took part in at least one recorded joint operation during the study period. When an operation involved more than two organizations, the data construction records each pair, so one operation can create several ties and a triangle. This whole-period binary projection is a pedagogical reanalysis, not a replication of the article’s square-root count model. It describes recorded tactical cooperation, not ideological agreement or a stable coalition.

load("data/gade_blocks.rda")
gade_weighted_net <- netify(
  gadeData,
  symmetric = TRUE,
  diag_to_NA = TRUE,
  missing_to_zero = FALSE
)
gade_binary_net <- binarize(gade_weighted_net, threshold = 0)
stopifnot(all(unlist(validate_netify(gade_binary_net))))
gade_binary <- get_adjacency(gade_binary_net)
diag(gade_binary) <- 0L
CooperationNet <- to_statnet(gade_binary_net)

checkpoint(
  organizations = network.size(CooperationNet),
  cooperative_pairs = network.edgecount(CooperationNet),
  density = network.density(CooperationNet),
  weak_transitivity = gtrans(gade_binary, mode = "graph")
)
#> ------------------------------------------------------------------
#> CHECKPOINT: organizations = 31   |   cooperative_pairs = 86   |   density = 0.185   |   weak_transitivity = 0.357
#> ------------------------------------------------------------------

netify() now performs the same label, symmetry, diagonal, and missingness checks used elsewhere in the course. binarize() records the change from operation counts to any recorded operation, and to_statnet() creates the network object required by latentnet. The weak transitivity score is 0.357. It summarizes how often the observed binary ties satisfy a weak closure condition among triples that could violate it. It is a descriptive network statistic, not an estimate that one partnership caused another.

The substantive question is narrow: can the observed cooperation network be summarized as clusters of organizations with overlapping partners? We deliberately fit a simple two-dimensional map with no measured covariates or separate adjustment for organizations that cooperate broadly. Its job is to make the distance idea visible, not to provide a new explanation of the Syrian conflict.

ls2d <- cache_fit("gade_ldm_2d_v1", ergmm(
  CooperationNet ~ euclidean(d = 2),
  control = ergmm.control(sample.size = 4000, interval = 10,
                          burnin = 10000),
  seed = 6886, verbose = FALSE
))

This is MCMC, so we inspect the chain before reading the map. The package’s full diagnostic display is useful on a large screen but becomes cramped in a handout. The focused display below shows the trace and density for the intercept and the latent-position variance. We do not diagnose a coordinate trace because an equivalent map can rotate or reflect while the fitted distances stay the same.

ldm_diag <- data.frame(
  iteration = seq_len(nrow(ls2d$sample$beta)),
  Intercept = ls2d$sample$beta[, 1],
  `Latent-position variance` = ls2d$sample$Z.var[, 1],
  check.names = FALSE
)
ldm_diag_long <- melt(
  ldm_diag,
  id.vars = "iteration",
  variable.name = "quantity",
  value.name = "value"
)
ldm_trace <- ggplot(ldm_diag_long, aes(iteration, value)) +
  geom_line(color = "#18453B", linewidth = .35) +
  facet_wrap(~ quantity, scales = "free_y", ncol = 1) +
  labs(x = "Retained draw", y = NULL)
ldm_density <- ggplot(ldm_diag_long, aes(value)) +
  geom_density(fill = "#7BBD00", color = "#18453B", alpha = .55, linewidth = .7) +
  facet_wrap(~ quantity, scales = "free", ncol = 1) +
  labs(x = "Posterior draw", y = "Density")
ldm_trace + ldm_density

Focused MCMC diagnostics for the latent-distance fit. Each trace should move around a stable level, and the density shows where the retained draws spent their time.
round(coda::effectiveSize(coda::mcmc(ldm_diag[, -1])))
#>                Intercept Latent-position variance 
#>                     1512                     2130
gade_positions <- data.frame(
  organization = rownames(gade_binary),
  x = ls2d$mkl$Z[, 1],
  y = ls2d$mkl$Z[, 2]
)
gade_edge_index <- which(gade_binary == 1 & upper.tri(gade_binary), arr.ind = TRUE)
gade_edges <- data.frame(
  x = gade_positions$x[gade_edge_index[, 1]],
  y = gade_positions$y[gade_edge_index[, 1]],
  xend = gade_positions$x[gade_edge_index[, 2]],
  yend = gade_positions$y[gade_edge_index[, 2]]
)
ggplot(gade_positions, aes(x, y, label = organization)) +
  geom_segment(
    data = gade_edges,
    aes(x = x, y = y, xend = xend, yend = yend),
    inherit.aes = FALSE,
    color = "grey72",
    linewidth = .45
  ) +
  geom_point(color = "#18453B", fill = "#7BBD00", shape = 21, size = 3.1, stroke = .8) +
  ggrepel::geom_text_repel(
    color = "#18453B",
    size = 3.4,
    box.padding = .35,
    point.padding = .2,
    max.overlaps = Inf,
    seed = 6886
  ) +
  coord_equal() +
  labs(x = "Latent coordinate 1", y = "Latent coordinate 2")

The fitted two-dimensional latent-distance map for recorded tactical cooperation. Organizations with similar whole-network cooperation patterns tend to be placed near one another. The axis directions have no substantive meaning.

Read the map through partner lists, not as a recovered historical geography. Al-Nusrah Front and Ahrar al-Sham Islamic Movement sit in the same tightly connected part of the map. They recorded an operation with one another, and each also recorded operations with 18 of the same other organizations. That reinforces the broad tactical-core result from Day 9 rather than creating a new historical explanation. A single operation does not determine a position, and the horizontal and vertical directions are arbitrary. The intercept is 1, the fitted log-odds when latent distance is zero. The distance term is then subtracted as pairs move apart.

summary(ls2d)
#> 
#> ==========================
#> Summary of model fit
#> ==========================
#> 
#> Formula:   CooperationNet ~ euclidean(d = 2)
#> Attribute: edges
#> Model:     Bernoulli 
#> MCMC sample of size 4000, draws are 10 iterations apart, after burnin of 10000 iterations.
#> Covariate coefficients posterior means:
#>             Estimate    2.5%  97.5% 2*min(Pr(>0),Pr(<0))   
#> (Intercept)  0.99617 0.35122 1.6923               0.0015 **
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Overall BIC:        507.9192 
#> Likelihood BIC:     354.464 
#> Latent space/clustering BIC:     153.4552 
#> 
#> Covariate coefficients MKL:
#>               Estimate
#> (Intercept) -0.4155834

The useful result is that the model can summarize much of the clustering in the recorded cooperation network. It does not give us a substantively new account of why these organizations worked together, so we use it as a bridge and move on. The optional identification box explains why the picture’s orientation itself should not be interpreted.

The distance geometry should make clustered triples more common. We can check that claim by simulating one network from each of 100 retained posterior draws and recalculating the same weak transitivity statistic. We build the simulation directly from plogis(intercept - distance) so the relationship between a posterior draw and a replicated network stays visible.

set.seed(6886)
draw_ids <- round(seq(1, dim(ls2d$sample$Z)[1], length.out = 100))
tr_sim <- sapply(draw_ids, function(s) {
  Pm <- plogis(ls2d$sample$beta[s] - as.matrix(dist(ls2d$sample$Z[s, , ])))
  n  <- nrow(Pm)
  Ys <- (matrix(runif(n^2), n, n) < Pm) * 1
  Ys[lower.tri(Ys, diag = TRUE)] <- 0
  Ys <- Ys + t(Ys)
  gtrans(Ys, mode = "graph")
})
checkpoint(
  observed_transitivity = gtrans(gade_binary, mode = "graph"),
  ldm_simulated_mean = mean(tr_sim),
  independent_tie_baseline = mean(gade_binary[upper.tri(gade_binary)])
)
#> ------------------------------------------------------------------
#> CHECKPOINT: observed_transitivity = 0.357   |   ldm_simulated_mean = 0.287   |   independent_tie_baseline = 0.185
#> ------------------------------------------------------------------

An independent-tie model at the observed density would have weak transitivity near 0.18. The latent-distance simulations average 0.29, compared with 0.36 in the observed network. The geometry moves the replicated networks toward the observed clustering, although the posterior predictive interval from 0.17 to 0.39 shows substantial variation. This check says the model can reproduce much of the triadic pattern. It does not show that social closure produced the observed operations.

The core path stops here. If the distance model is useful for your project, these commands expose the dimension-selection summary and fitted tie probabilities:

bic.ergmm(ls2d)$overall
#> [1] 507.9192
round(predict(ls2d)[1:4, 1:4], 2)
#>      [,1] [,2] [,3] [,4]
#> [1,] 0.73 0.00 0.00 0.00
#> [2,] 0.42 0.73 0.00 0.00
#> [3,] 0.09 0.09 0.73 0.00
#> [4,] 0.10 0.10 0.43 0.73

Dyadic and nodal covariates enter through edgecov() and nodecov(). Model-based clustering on top of the space uses euclidean(d = 2, G = k) (Handcock, Raftery, and Tantrum 2007). latentnet also provides a bilinear() term, which connects directly to the factor model later in the walkthrough. Krivitsky and Handcock (2008) provide the detailed software reference.

latentnet also ships the Sampson monastery data. Running data(sampson) creates an object called samplike, a directed network of 18 monks and their positive nominations. The corresponding call is ergmm(samplike ~ euclidean(d = 2)). The same estimation and identification lessons apply, but sender-receiver direction now matters.

4 Where the Map Struggles

Now the limitation. The same geometry that favors clustered ties makes a low-dimensional distance model an awkward language for sharp disassortative roles. If two actors have nearly identical tie profiles but do not tie to each other, a one-map distance model faces a tension: similar profiles encourage similar locations, while similar locations also raise their direct tie probability. A factor or block model can separate those two ideas cleanly.

This is a representational and parsimony claim, not an absolute impossibility claim. With enough dimensions, a finite Euclidean configuration can mimic more structure than a two-dimensional picture suggests. Hoff (2007) shows that doing so for role structure can require dimension that grows with the pattern being represented. The practical question is whether the small dd we can estimate and interpret represents the pattern compactly. Rather than answer in the abstract, build a network that puts that small model under stress: two groups, ties almost always between, almost never within, as in a buyer-seller or patron-client system.

set.seed(6886)
same  <- outer(grp, grp, "==")
Pdis  <- ifelse(same, 0.10, 0.85)
Ydis  <- (matrix(runif(np^2), np, np) < Pdis) * 1
Ydis[lower.tri(Ydis, diag = TRUE)] <- 0
Ydis  <- Ydis + t(Ydis)
diag(same) <- NA
image(1:np, 1:np, Ydis,
      breaks = c(-0.5, 0.5, 1.5),
      col = c("white", "#18453B"),
      xlab = "actor", ylab = "actor",
      main = "ties almost only BETWEEN groups (disassortative)")

The disassortative test network. Ties (dark cells) run almost entirely BETWEEN the two groups. To make between-group pairs likely, a distance model must put them close, but then within-group pairs get dragged close too, and it predicts within-group ties that do not exist.
obs_within  <- mean(Ydis[same],  na.rm = TRUE)
obs_between <- mean(Ydis[!same], na.rm = TRUE)
checkpoint(
  observed_within_group_density  = obs_within,
  observed_between_group_density = obs_between
)
#> ------------------------------------------------------------------
#> CHECKPOINT: observed_within_group_density = 0.097   |   observed_between_group_density = 0.815
#> ------------------------------------------------------------------

The observed contrast is sharp: within-group density is 0.1, while between-group density is 0.81. We now fit the same two-dimensional distance model used for the cooperation network and calculate the fitted within-group and between-group probabilities from each posterior draw.

ldm_dis <- cache_fit("ldm_dis", ergmm(
  network(Ydis, directed = FALSE) ~ euclidean(d = 2),
  control = ergmm.control(sample.size = 2000, interval = 10,
                          burnin = 5000),
  seed = 6886, verbose = FALSE
))
pp <- sapply(seq_len(dim(ldm_dis$sample$Z)[1]), function(s) {
  Pm <- plogis(ldm_dis$sample$beta[s] -
               as.matrix(dist(ldm_dis$sample$Z[s, , ])))
  c(within  = mean(Pm[same],  na.rm = TRUE),
    between = mean(Pm[!same], na.rm = TRUE),
    overall = mean(Pm[row(Pm) != col(Pm)]))
})
checkpoint(
  ldm_predicted_within  = mean(pp["within", ]),
  ldm_predicted_between = mean(pp["between", ]),
  observed_within       = obs_within,
  observed_between      = obs_between
)
#> ------------------------------------------------------------------
#> CHECKPOINT: ldm_predicted_within = 0.461   |   ldm_predicted_between = 0.47   |   observed_within = 0.097   |   observed_between = 0.815
#> ------------------------------------------------------------------

The fitted distance model predicts 0.46 within-group and 0.47 between-group tie probabilities. It nearly removes the contrast that generated the data. To make cross-group ties likely, the model needs actors from the two groups near one another. Once they are close, the same metric geometry also raises probabilities within each group.

The average fitted probability is 0.47, close to the observed overall density of 0.47. The model matches the average while missing the defining within-between pattern. This is evidence that this two-dimensional Euclidean distance model does not represent this role structure well. It is not a claim that every Euclidean model at every dimension must fail.

A poor run could also flatten a meaningful contrast. The following check repeats the specification with a longer warm-up and recomputes the posterior probabilities.

ldm_dis_long <- cache_fit("ldm_dis_long", ergmm(
  network(Ydis, directed = FALSE) ~ euclidean(d = 2),
  control = ergmm.control(sample.size = 2000, interval = 10,
                          burnin = 30000),
  seed = 6886, verbose = FALSE
))
pp_long <- sapply(seq_len(dim(ldm_dis_long$sample$Z)[1]), function(s) {
  Pm <- plogis(ldm_dis_long$sample$beta[s] -
               as.matrix(dist(ldm_dis_long$sample$Z[s, , ])))
  c(within  = mean(Pm[same],  na.rm = TRUE),
    between = mean(Pm[!same], na.rm = TRUE))
})
checkpoint(
  long_burnin_within  = mean(pp_long["within", ]),
  long_burnin_between = mean(pp_long["between", ])
)
#> ------------------------------------------------------------------
#> CHECKPOINT: long_burnin_within = 0.46   |   long_burnin_between = 0.469
#> ------------------------------------------------------------------

The longer run gives 0.46 within and 0.47 between. The result is unchanged to two decimals. This check does not prove global convergence, but it does show that simply extending this run does not recover the missing contrast. Changing the dimension or model family is a separate decision.

If these point estimates recovered the group structure, a simple two-cluster summary of the positions should separate the known groups. It does not:

truth   <- as.integer(factor(grp))
km_ldm  <- kmeans(ldm_dis$mkl$Z, 2, nstart = 50)
acc_ldm <- max(mean(km_ldm$cluster == truth), mean(km_ldm$cluster == 3 - truth))
checkpoint(kmeans_on_LDM_positions_accuracy = acc_ldm)
#> ------------------------------------------------------------------
#> CHECKPOINT: kmeans_on_LDM_positions_accuracy = 0.55
#> ------------------------------------------------------------------

The classification accuracy is 55%, approximately chance for this balanced example. This is not a general test of latent-distance models. It confirms that these two-dimensional point estimates did not retain the simulated role labels, while Section 6 shows that a factor surface separates the same roles much more clearly.

The k-means check is deliberately simple. latentnet also provides model-based clustering that fits a mixture over the latent positions (Handcock, Raftery, and Tantrum 2007). The mixture is valuable for assortative clusters, but it still operates inside the same Euclidean distance geometry. We can check it on the simulated role network:

hrt_dis <- cache_fit("hrt_dis", ergmm(
  network(Ydis, directed = FALSE) ~ euclidean(d = 2, G = 2),
  seed = 6886, verbose = FALSE
))
cl_hrt  <- hrt_dis$mkl$Z.K            # The model-based cluster assignments
acc_hrt <- max(mean(cl_hrt == truth), mean(cl_hrt == 3 - truth))
checkpoint(clusters_found = length(unique(cl_hrt)),
           hrt_cluster_accuracy = acc_hrt)
#> ------------------------------------------------------------------
#> CHECKPOINT: clusters_found = 1   |   hrt_cluster_accuracy = 0.5
#> ------------------------------------------------------------------

The fitted mixture’s classification accuracy is 50%. In this run it assigns every actor to one cluster. The result is specific to this data-generating pattern and specification, but it reinforces the representational point: adding a mixture does not make a low-dimensional distance geometry equivalent to a model built for across-role ties.

ImportantThe Midpoint Summary

The distance model is built for latent homophily and tends to generate clustered ties. A small Euclidean space can represent sharp stochastic equivalence or disassortative mixing inefficiently, and this two-dimensional fit missed the simulated contrast. Day 9’s blockmodel handled roles directly. The eigenmodel in Section 5 gives us a continuous low-rank representation that can handle both within-role and across-role patterns.

Open this section when you want the formal explanation of why latent maps can face different directions and how we compare equivalent configurations. The exact invariances here belong to a distance map. The factor model carries the same broad warning, “the coordinates are not the estimand,” but its directed version has a larger invariance than an ordinary rotation. We handle that separately in Section 7 and Section 7.7.

Here is the problem, and it is not a technicality. The model is only identified up to rotation, reflection, and translation of the latent positions. The likelihood depends on the positions only through the pairwise distances ‖ui−uj‖\lVert u_i - u_j \rVert. And distances do not change if you spin the whole configuration, flip it in a mirror, or slide it across the page.

So there is not one set of positions that fits your network. There are infinitely many: every rotated, reflected, shifted copy fits exactly as well. The plot latentnet drew you in Section 3 is one arbitrary member of that infinite family.

We can see the problem directly. Take the cooperation positions, rotate them by 60 degrees, reflect them, and move them across the page. Every pairwise distance, fitted probability, and likelihood contribution stays unchanged.

# ls2d$mkl$Z is latentnet's single best point-estimate of the positions; the
# MKL (minimum Kullback-Leibler) configuration, the one map it picks to
# summarize the whole posterior.
Z <- ls2d$mkl$Z

theta  <- pi / 3                          # 60 degrees
Rot    <- matrix(c(cos(theta), sin(theta), -sin(theta), cos(theta)), 2)
Reflect<- diag(c(1, -1))                  # Mirror the 2nd axis
Z_alt  <- Z %*% Rot %*% Reflect
Z_alt  <- sweep(Z_alt, 2, c(4, -3), "+")  # Translate

par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(Z, pch = 19, col = "#18453B", asp = 1, main = "reported",
     xlab = "dim 1", ylab = "dim 2")
plot(Z_alt, pch = 19, col = "#7BBD00", asp = 1, main = "rotated + reflected + shifted",
     xlab = "dim 1", ylab = "dim 2")

Left: the positions latentnet reported. Right: the SAME positions rotated 60 degrees, reflected, and translated. Different-looking map, identical model. Anything you ‘read off’ the axes on the left is an artifact.

The two point clouds look nothing alike. But the pairwise distance matrices are, to numerical precision, the same object:

max_dist_diff <- max(abs(as.matrix(dist(Z)) - as.matrix(dist(Z_alt))))
checkpoint(
  max_abs_difference_in_all_pairwise_distances = max_dist_diff
)
#> ------------------------------------------------------------------
#> CHECKPOINT: max_abs_difference_in_all_pairwise_distances = 0
#> ------------------------------------------------------------------

The maximum difference is 8.9e-16, which is zero to numerical precision. The model cannot distinguish these two maps because they encode the same fitted distances.

ImportantWhat You May and May Not Read From a Latent-Space Plot

Defensible: which actors are close, which sets cluster, and how distances compare. These quantities are unchanged by rotation, reflection, and translation.

Not identified without outside information: “the x-axis is ideology,” “this actor is on the left,” or “there is a north-south split.” Absolute direction and orientation come from the displayed solution. If an interpretation changes after a mirror reflection, it is not an invariant model result.

This box applies exactly to distance models. The reliable practice travels to every latent model: identify which functions of the latent variables leave the likelihood unchanged, and interpret only quantities that survive those transformations.

4.1 What Is Procrustes Alignment?

If two maps that differ only by rotation/reflection/translation are the same model, then to compare two fitted maps across seeds, across models, or across your fit and a coauthor’s, you first have to strip away that arbitrary difference. That alignment step is called a Procrustes rotation: rotate, reflect, and shift one configuration to line up as closely as possible with a target, then look at whatever is left.

(The name comes from a Greek myth: Procrustes stretched or lopped his guests to fit an iron bed. Here we are allowed to rotate and flip a configuration but not stretch it; the shape is sacred, only the pose is negotiable.)

And here is the useful connection to the second half of the day: the optimal Procrustes rotation is computed with a singular value decomposition. This is the same matrix factorization we use in the movie example in Section 5.2. To align XX to a target YY after centering both matrices, take the SVD of X⊤Y=UΣV⊤X^\top Y = U \Sigma V^\top; the best rotation is R=UV⊤R = U V^\top. That is the whole algorithm.

# Align X onto target, allowing rotation + reflection + translation
# (but NOT rescaling: distances are the model, so we may not stretch them).
align_map_to_target <- function(X, target) {
  Xc <- scale(X,      center = TRUE, scale = FALSE)
  Tc <- scale(target, center = TRUE, scale = FALSE)
  s  <- svd(t(Xc) %*% Tc)          # <- The SVD doing the work
  R  <- s$u %*% t(s$v)             # Optimal orthogonal (rotation/reflection)
  aligned <- Xc %*% R
  list(aligned = aligned,
       mse = mean((Tc - aligned)^2))
}

Run it on our deliberately-mangled map. It should recover the original exactly, because the mangling was only rotation/reflection/translation; nothing Procrustes cannot undo.

rec <- align_map_to_target(Z_alt, Z)
mse_before <- mean((scale(Z, scale = FALSE) -
                    scale(Z_alt, scale = FALSE))^2)
checkpoint(
  mse_before_alignment = mse_before,
  mse_after_alignment  = rec$mse
)
#> ------------------------------------------------------------------
#> CHECKPOINT: mse_before_alignment = 1.015   |   mse_after_alignment = 0
#> ------------------------------------------------------------------

Before alignment the two maps disagree by 1.01; after Procrustes they agree to 5.0e-32, effectively machine zero. The difference between the two plots was entirely the arbitrary pose, and Procrustes removed all of it.

The demo above rotated the same positions, so Procrustes recovers them perfectly. Independent runs are messier: start the sampler twice from different seeds and you get two potentially different point estimates, with contributions from posterior uncertainty, Monte Carlo error, weak identification, and possibly different local modes. Procrustes removes only the orthogonal pose. A residual after alignment tells you that the two summaries differ beyond that pose; it does not diagnose which of those sources caused the difference. This is exactly why latentnet aligns internally: plot() and summary() align every posterior draw to a common reference (the MKL configuration) before averaging, which is the only reason the plotted positions are stable at all. Run the two-chain version yourself as exercise 3.

Rotation and reflection are the easy part of non-identifiability; they are inside the model and Procrustes fixes them. There is a harder layer underneath.

We chose d=2d = 2. We chose Euclidean distance: flat space. We chose a space with no curvature. None of these came from the data; we imposed them before the model saw a single tie. Lubold, Chandrasekhar, and McCormick (2023) make the uncomfortable case precise: the dimension, the geometry, and the curvature of the latent space are a priori modeling assumptions, and different manifolds can fit the same network comparably well while telling different substantive stories. A network that looks like tight clustering in flat 2-D space might be perfectly ordinary points on a sphere.

The practical discipline that follows:

  • Report the dimension you chose and why (fit? Interpretability? A theory?).
  • Do not treat a 2-D picture as “the” geometry of your network; it is the geometry conditional on your assuming 2-D flat Euclidean space.
  • When curvature could plausibly matter (hierarchies, cores-and-peripheries), say so, and know that the tools to test it are recent and still developing.

This does not make these models useless; it makes their assumptions explicit. You are fitting a map under assumptions, and naming the assumptions is the whole job. The factor model we are about to build also carries a chosen rank RR and an identification convention. Its directed invariance is broader than rotation, so nothing this afternoon is exempt from the habit even though the exact algebra changes.

5 Two Kinds of Similarity and the Factor Model

The fix starts by naming precisely what the distance model could not say. There are two classic ways third-order structure shows up, and they imply different theories of what makes two actors “alike.”

hom_nodes <- data.frame(
  id = 1:8,
  x = c(-1.4, -1, -.6, -1, .6, 1, 1.4, 1),
  y = c(0, .6, 0, -.6, 0, .6, 0, -.6),
  group = rep(c("Group 1", "Group 2"), each = 4)
)
hom_pairs <- rbind(
  t(combn(1:4, 2)),
  t(combn(5:8, 2))
)
hom_edges <- data.frame(
  x = hom_nodes$x[hom_pairs[, 1]],
  y = hom_nodes$y[hom_pairs[, 1]],
  xend = hom_nodes$x[hom_pairs[, 2]],
  yend = hom_nodes$y[hom_pairs[, 2]]
)
hom_plot <- ggplot(hom_nodes, aes(x, y)) +
  geom_segment(
    data = hom_edges,
    aes(x = x, y = y, xend = xend, yend = yend),
    inherit.aes = FALSE,
    color = "#535054",
    linewidth = .8
  ) +
  geom_point(aes(fill = group), shape = 21, color = "#18453B", size = 5, stroke = 1) +
  scale_fill_manual(values = c("Group 1" = "#18453B", "Group 2" = "#7BBD00")) +
  coord_equal(xlim = c(-1.8, 1.8), ylim = c(-1, 1)) +
  labs(title = "Homophily: Similar Actors Connect", x = NULL, y = NULL, fill = NULL) +
  theme_void(base_size = 13) +
  theme(legend.position = "bottom", plot.title = element_text(color = "#18453B", face = "bold", hjust = .5))

role_nodes <- data.frame(
  id = c("A", "B", "1", "2", "3"),
  x = c(-1, -1, 1, 1, 1),
  y = c(.65, -.65, .9, 0, -.9),
  type = c("Same role", "Same role", "Shared partners", "Shared partners", "Shared partners")
)
role_pairs <- expand.grid(from = c("A", "B"), to = c("1", "2", "3"))
role_edges <- merge(role_pairs, role_nodes[, c("id", "x", "y")], by.x = "from", by.y = "id")
names(role_edges)[3:4] <- c("x", "y")
role_edges <- merge(role_edges, role_nodes[, c("id", "x", "y")], by.x = "to", by.y = "id")
names(role_edges)[names(role_edges) == "x.x"] <- "x"
names(role_edges)[names(role_edges) == "y.x"] <- "y"
names(role_edges)[names(role_edges) == "x.y"] <- "xend"
names(role_edges)[names(role_edges) == "y.y"] <- "yend"
role_plot <- ggplot(role_nodes, aes(x, y)) +
  geom_segment(
    data = role_edges,
    aes(x = x, y = y, xend = xend, yend = yend),
    inherit.aes = FALSE,
    color = "#535054",
    linewidth = .8
  ) +
  geom_point(aes(fill = type), shape = 21, color = "#18453B", size = 5, stroke = 1) +
  geom_text(aes(label = id, color = type), fontface = "bold", size = 4) +
  scale_fill_manual(values = c("Same role" = "#18453B", "Shared partners" = "#7BBD00")) +
  scale_color_manual(values = c("Same role" = "white", "Shared partners" = "#18453B"), guide = "none") +
  coord_equal(xlim = c(-1.5, 1.5), ylim = c(-1.2, 1.2)) +
  labs(title = "Stochastic Equivalence: Same Partners", x = NULL, y = NULL, fill = NULL) +
  theme_void(base_size = 13) +
  theme(legend.position = "bottom", plot.title = element_text(color = "#18453B", face = "bold", hjust = .5))

hom_plot + role_plot

Two meanings of similarity. In the homophily panel, actors connect within the same group. In the stochastic-equivalence panel, A and B occupy the same role because they connect to the same partners, even though A and B do not connect to one another.

Homophily: birds of a feather. Actors are more likely to tie to actors who are like them. This is the distance story. If ii is very near both jj and kk, the triangle inequality constrains how far apart jj and kk can be, so a decreasing-distance tie probability tends to produce clustered triads. It raises a probability; it does not turn two observed ties into a logical proof of the third.

Stochastic equivalence: actors occupy the same role: they tie to the same others without necessarily tying to each other. Two junior faculty members may both work with senior colleagues and graduate students without collaborating with one another. In an alliance network, two states can have similar sets of allies without being allied to each other. Day 9’s blockmodel represented this idea with discrete groups; Section 4’s disassortative network is its extreme case.

Important

Homophily says similar actors connect to each other. Stochastic equivalence says similar actors connect to the same others. The distance model is a compact language for the first; a small distance space can need many more dimensions to approximate the second. Section 4 was a concrete failure at d=2d=2, not a universal impossibility proof. The multiplicative term represents both patterns compactly, and that is the whole reason it exists.

5.1 Build the Factor Model From a Regression

The latent factor model is easiest to understand as a regression with an unusual interaction term. The outcome and link function do not disappear. We keep the parts of regression that are already familiar and add one new way for a source and target to fit together.

Let μij=E(Yij∣⋅)\mu_{ij}=E(Y_{ij}\mid\cdot) and let g(⋅)g(\cdot) be the link function. For a continuous outcome, gg can be the identity link. For a binary outcome, it can be a logit or probit link. An ordinary dyadic regression starts with

g(μij)=ηij=β0+𝐱ij𝖳𝛃. g(\mu_{ij})=\eta_{ij}=\beta_0+\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}.

The linear predictor ηij\eta_{ij} is the model’s score for the ordered relationship from ii to jj. The intercept supplies a baseline, and the measured covariates move that score up or down.

Day 9 added broad source and target differences:

ηij=β0+𝐱ij𝖳𝛃+ai+bj. \eta_{ij}=\beta_0+\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+a_i+b_j.

The source effect aia_i shifts every relationship that begins with actor ii. The target effect bjb_j shifts every relationship aimed at actor jj. These terms explain why some rows or columns are broadly high, but they cannot explain why one particular source-target combination is unexpectedly high while another combination involving the same source is unexpectedly low.

5.1.1 Start With an Ordinary Measured Interaction

Suppose we observe a source characteristic cic_i, such as military capability, and a target characteristic rjr_j, such as vulnerability. Assume their separate source and target effects are already included in 𝐱ij\mathbf{x}_{ij}. A familiar regression can then include their product:

ηij=β0+𝐱ij𝖳𝛃+ai+bj+δcirj. \eta_{ij}=\beta_0+\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+a_i+b_j+\delta\,c_i r_j.

Nothing mysterious has happened. We create the dyadic predictor cirjc_i r_j, estimate its coefficient δ\delta, and ask whether the combination of source capability and target vulnerability raises or lowers the relationship score. The same value of cic_i can matter differently across targets because it is multiplied by rjr_j.

5.1.2 Now Let the Interacting Characteristics Be Unobserved

The factor model keeps that multiplication but admits that the relevant source and target characteristics may not be measured:

ηij=β0+𝐱ij𝖳𝛃+ai+bj+δuivj. \eta_{ij}=\beta_0+\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+a_i+b_j+\delta\,u_i v_j.

This is a rank-one latent factor model. The source score uiu_i and target score vjv_j play the same algebraic roles as cic_i and rjr_j, but they are not columns in a dataset. The model estimates them from the repeated pattern across the whole relationship matrix. It chooses source scores and target scores whose products reconstruct systematic structure left after the measured predictors and additive effects enter.

The distinction from an ordinary interaction is therefore simple but important:

Ordinary interaction Latent interaction
cic_i and rjr_j are observed uiu_i and vjv_j must be estimated
cirjc_i r_j is a known predictor uivju_i v_j is learned from the relationship matrix
The coefficient is the main unknown The regression coefficients and both sets of profiles are unknown
The variables arrive with substantive names A latent direction needs outside evidence before it receives a substantive name

5.1.3 Rank Means More Than One Hidden Interaction

One source-target product can represent only one repeated pattern. With RR dimensions, the model adds RR products:

ηij=β0+𝐱ij𝖳𝛃+ai+bj+∑r=1Rδruirvjr=β0+𝐱ij𝖳𝛃+ai+bj+𝐮i𝖳D𝐯j, \eta_{ij} =\beta_0+\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+a_i+b_j +\sum_{r=1}^{R}\delta_r u_{ir}v_{jr} =\beta_0+\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+a_i+b_j +\mathbf{u}_i^{\mathsf T}D\mathbf{v}_j,

where D=diag(δ1,…,δR)D=\operatorname{diag}(\delta_1,\ldots,\delta_R). Each dimension is another latent interaction. For directed models, software often absorbs the weights in DD into the scales of UU or VV and reports the equivalent product UV𝖳UV^{\mathsf T}. For symmetric eigenmodels, keeping the diagonal weights visible is especially useful because positive and negative eigenvalues distinguish assortative and disassortative components.

Here is a rank-one numerical example. Suppose actor A has source score uA=2u_A=2, actor B has uB=−1u_B=-1, target C has vC=1.5v_C=1.5, and target D has vD=−2v_D=-2. With δ=1\delta=1, the latent contributions are:

Ordered relationship Calculation Latent contribution
A →\rightarrow C 2(1.5)2(1.5) 33
A →\rightarrow D 2(−2)2(-2) −4-4
B →\rightarrow C −1(1.5)-1(1.5) −1.5-1.5
B →\rightarrow D −1(−2)-1(-2) 22

Actor A does not receive one universal positive factor effect. Its factor contribution is positive with C and negative with D. That is why this term represents pair-specific compatibility rather than broad source activity.

ImportantThe Regression Connection

AME is not a replacement for regression. It extends the linear predictor with source effects, target effects, and a small number of latent interactions, while the residual structure can also represent reciprocity within directed pairs. If UU and VV were known, their products could enter a regression like any other interaction predictors. The hard part is that they are unknown and must be estimated jointly with the rest of the model.

This regression view also makes the estimation logic less intimidating. Holding the target profiles VV fixed makes the model regression-like in the source profiles UU; holding UU fixed makes it regression-like in VV. ALS alternates between those jobs while updating the other parameters. MCMC instead samples each block conditional on the current values of the others. The movie example now gives us a concrete way to see how repeated row-column patterns reveal those profiles.

5.2 Start With a Movie Preference Matrix

The factor idea is easier to see in a rectangular matrix before we return to a square network. Rows are viewers, columns are movies, and each cell is a recorded preference score from 0 to 5. A zero here means a recorded score of zero. It does not mean that the movie was unseen or that the value is missing. That distinction matters because an ordinary SVD treats every number in the matrix as observed.

movie_scores <- matrix(
  c(
    1, 1, 1, 0, 0,
    3, 3, 3, 0, 0,
    4, 4, 4, 0, 0,
    5, 5, 5, 0, 0,
    0, 2, 0, 4, 4,
    0, 0, 0, 5, 5,
    0, 1, 0, 2, 2
  ),
  ncol = 5,
  byrow = TRUE,
  dimnames = list(
    c("Mike", "Cindy", "Hyerin", "Emily", "Cassy", "Juan", "Max"),
    c("Star Wars", "Alien", "Blade Runner", "Casablanca", "Pretty Woman")
  )
)
movie_scores
#>        Star Wars Alien Blade Runner Casablanca Pretty Woman
#> Mike           1     1            1          0            0
#> Cindy          3     3            3          0            0
#> Hyerin         4     4            4          0            0
#> Emily          5     5            5          0            0
#> Cassy          0     2            0          4            4
#> Juan           0     0            0          5            5
#> Max            0     1            0          2            2

We could describe 35 cells one at a time, but the repetitions are obvious. Mike, Cindy, Hyerin, and Emily differ mainly in how strongly they prefer the same three science-fiction movies. Juan favors the two romantic films. Cassy and Max mix the two patterns. The question behind SVD is: can a small number of repeated row and column patterns reconstruct almost all of this matrix?

5.3 What the SVD Separates

For any complete numeric matrix MM with nn rows and mm columns, the singular value decomposition writes

M=UDV𝖳. M=U D V^{\mathsf T}.

Read the three pieces by the job they do:

  • The rows of UU describe how each row unit loads on the shared directions. Here the row units are viewers.
  • The rows of VV describe how each column unit loads on those same directions. Here the column units are movies.
  • The diagonal entries of DD are nonnegative singular values. They order the directions by how much squared matrix magnitude they reconstruct.

Keeping only the first RR directions gives

MR=URDRVR𝖳. M_R=U_R D_R V_R^{\mathsf T}.

It is often clearer to split each singular value evenly across the two sides. Define P=URDR1/2P=U_R D_R^{1/2} and Q=VRDR1/2Q=V_R D_R^{1/2}. Then

MR=PQ𝖳,m̂ij=𝐩i𝖳𝐪j=∑r=1Rpirqjr. M_R=P Q^{\mathsf T}, \qquad \widehat m_{ij}=\mathbf p_i^{\mathsf T}\mathbf q_j =\sum_{r=1}^{R}p_{ir}q_{jr}.

That final expression is the factor-model intuition. Viewer ii does not carry one universal “likes movies” score, and movie jj does not carry one universal “good movie” score. Their fitted cell depends on how the viewer profile and movie profile line up, dimension by dimension.

movie_svd <- svd(movie_scores)
movie_rank <- 2L
movie_P <- sweep(
  movie_svd$u[, seq_len(movie_rank), drop = FALSE],
  2,
  sqrt(movie_svd$d[seq_len(movie_rank)]),
  `*`
)
movie_Q <- sweep(
  movie_svd$v[, seq_len(movie_rank), drop = FALSE],
  2,
  sqrt(movie_svd$d[seq_len(movie_rank)]),
  `*`
)
rownames(movie_P) <- rownames(movie_scores)
rownames(movie_Q) <- colnames(movie_scores)
colnames(movie_P) <- colnames(movie_Q) <- paste0("Direction ", seq_len(movie_rank))
movie_rank2 <- movie_P %*% t(movie_Q)
movie_svd_summary <- data.frame(
  direction = seq_along(movie_svd$d),
  singular_value = movie_svd$d,
  cumulative_squared_magnitude = cumsum(movie_svd$d^2) / sum(movie_svd$d^2)
)
round(movie_svd_summary, 3)
direction singular_value cumulative_squared_magnitude
1 12.481 0.628
2 9.509 0.993
3 1.346 1.000
4 0.000 1.000
5 0.000 1.000 
checkpoint(
  rank_two_share_of_squared_matrix_magnitude =
    round(sum(movie_svd$d[1:2]^2) / sum(movie_svd$d^2), 3),
  rank_two_rmse =
    round(sqrt(mean((movie_scores - movie_rank2)^2)), 3)
)
#> ------------------------------------------------------------------
#> CHECKPOINT: rank_two_share_of_squared_matrix_magnitude = 0.993   |   rank_two_rmse = 0.227
#> ------------------------------------------------------------------

The first two directions retain about 99.3% of the squared magnitude of this uncentered toy matrix. That phrase is more accurate than “variance explained” because we did not center the matrix first. It does not mean that two true psychological traits generated the ratings. It means that a rank-two matrix nearly reproduces these 35 recorded scores under squared-error loss.

movie_matrix_long <- function(M, label) {
  out <- melt(M, varnames = c("viewer", "movie"), value.name = "score")
  out$matrix <- label
  out
}
movie_heat <- rbind(
  movie_matrix_long(movie_scores, "Observed Scores"),
  movie_matrix_long(movie_rank2, "Rank-Two Reconstruction")
)
movie_heat$score_label <- ifelse(
  abs(movie_heat$score) < .05,
  "0.0",
  sprintf("%.1f", movie_heat$score)
)
movie_heat$viewer <- factor(movie_heat$viewer, levels = rev(rownames(movie_scores)))
movie_heat$movie <- factor(movie_heat$movie, levels = colnames(movie_scores))
ggplot(movie_heat, aes(movie, viewer, fill = score)) +
  geom_tile(color = "white", linewidth = .5) +
  geom_text(aes(label = score_label), size = 3.4) +
  facet_wrap(~ matrix, nrow = 1) +
  scale_fill_gradient2(
    low = "#C7D4C8",
    mid = "white",
    high = "#18453B",
    midpoint = 2.5,
    limits = c(-.5, 5)
  ) +
  labs(x = NULL, y = NULL, fill = "Score") +
  theme(
    axis.text.x = element_text(angle = 35, hjust = 1),
    panel.grid = element_blank()
  )

The observed preference matrix and its rank-two SVD reconstruction. Two repeated row-column patterns reproduce almost all of the toy matrix, including the broad science-fiction and romantic-film split.

5.4 Read a Fitted Cell as a Product

The profile plot below uses the balanced coordinates PP and QQ. The orientation is arbitrary, so do not call the horizontal direction “science fiction” or the vertical direction “romance.” What matters is how viewer and movie vectors line up with one another and the fitted cells their inner products produce.

movie_profiles <- rbind(
  data.frame(name = rownames(movie_P), movie_P, side = "Viewer Profiles", check.names = FALSE),
  data.frame(name = rownames(movie_Q), movie_Q, side = "Movie Profiles", check.names = FALSE)
)
names(movie_profiles)[2:3] <- c("direction_1", "direction_2")
ggplot(movie_profiles, aes(direction_1, direction_2, color = side)) +
  geom_hline(yintercept = 0, color = "grey80") +
  geom_vline(xintercept = 0, color = "grey80") +
  geom_segment(
    aes(x = 0, y = 0, xend = direction_1, yend = direction_2),
    arrow = arrow(length = unit(.12, "inches")),
    linewidth = .7
  ) +
  geom_point(size = 2.5) +
  ggrepel::geom_text_repel(aes(label = name), size = 3.4, show.legend = FALSE) +
  facet_wrap(~ side, nrow = 1) +
  scale_color_manual(values = c("Viewer Profiles" = "#18453B", "Movie Profiles" = "#7BBD00")) +
  coord_equal() +
  labs(x = "Direction 1", y = "Direction 2", color = NULL) +
  theme(legend.position = "bottom")

Balanced rank-two SVD profiles. Viewers and movies that point in compatible directions receive larger reconstructed scores. The axes may flip or rotate without changing the reconstructed matrix.

The multiplication is literal. For a rank-two fit, the Emily and Star Wars reconstruction is

m̂Emily, Star Wars=pEmily,1qStar Wars,1+pEmily,2qStar Wars,2. \widehat m_{\text{Emily, Star Wars}} =p_{\text{Emily},1}q_{\text{Star Wars},1} +p_{\text{Emily},2}q_{\text{Star Wars},2}.

The same viewer coordinate can raise the fitted score for one movie and do little for another because it is always multiplied by the movie’s coordinate. That pair-specific interaction is what an additive row effect or additive column effect cannot provide.

movie_cell_example <- function(viewer, movie) {
  pieces <- movie_P[viewer, ] * movie_Q[movie, ]
  data.frame(
    viewer = viewer,
    movie = movie,
    direction_1 = pieces[1],
    direction_2 = pieces[2],
    fitted_score = sum(pieces),
    observed_score = movie_scores[viewer, movie],
    row.names = NULL
  )
}
movie_cell_examples <- rbind(
  movie_cell_example("Emily", "Star Wars"),
  movie_cell_example("Juan", "Casablanca"),
  movie_cell_example("Cassy", "Alien")
)
movie_cell_examples[-c(1, 2)] <- round(
  movie_cell_examples[-c(1, 2)],
  2
)
movie_cell_examples
viewer movie direction_1 direction_2 fitted_score observed_score
Emily Star Wars 4.83 0.14 4.97 5
Juan Casablanca 0.08 4.84 4.92 5
Cassy Alien 1.13 0.16 1.29 2

5.5 What SVD Is Actually Doing

The truncated SVD is not searching for a label such as genre. It is solving a reconstruction problem. Among all matrices with rank at most RR, MRM_R is the one that minimizes the sum of squared cell-by-cell errors. The first direction captures the largest remaining squared pattern, the second captures the largest pattern orthogonal to the first, and so on.

That result gives us a disciplined interpretation:

  1. Rows with similar profiles receive similar row coordinates. They tend to have similar fitted values across columns.
  2. Columns with similar profiles receive similar column coordinates. They tend to have similar fitted values across rows.
  3. A cell comes from compatibility. The row coordinates and column coordinates are multiplied and added.
  4. Rank controls compression. A small RR forces many cells to share a small number of repeated patterns.
  5. The axes are not named concepts by themselves. Genres are visible here because we already know the movie labels. In an application, a substantive label requires outside evidence.
WarningDo Not Confuse Zero With Missing

The movie scores are complete by construction. If zero meant “not observed,” an ordinary svd(movie_scores) call would incorrectly treat an unobserved score as dislike. Practical recommendation systems use methods designed for incomplete matrices or fit only the observed entries. The same discipline carries into networks: the risk set, structural zeros, and unobserved dyads must be defined before fitting the factors.

5.6 Carry the Same Logic Into Directed Relations

A movie matrix is rectangular because viewers and movies are different kinds of units. A directed network is square because the same actors appear in the rows and columns, but the two sides still have different jobs. In a state-event network, row ii describes state ii as a source and column jj describes state jj as a target. We therefore keep separate source profiles uiu_i and target profiles vjv_j:

zij≈ui𝖳vj=∑r=1Ruirvjr. z_{ij}\approx u_i^{\mathsf T}v_j =\sum_{r=1}^{R}u_{ir}v_{jr}.

The product does not say that ii is generally active or that jj is generally exposed. The additive effects from the SRM handle those broad differences. The product says that this source profile lines up especially well or badly with this target profile.

Data object Row profile describes Column profile describes What a large product means
Viewer-by-movie preference matrix What kinds of movies a viewer tends to prefer Which viewers tend to prefer a movie The viewer and movie patterns are compatible
Country-by-resolution voting matrix A country’s pattern across roll calls The coalition or issue pattern of a resolution The country profile aligns with that resolution
Donor-by-recipient aid matrix Which recipients a donor tends to fund Which donors tend to fund a recipient The donor and recipient patterns fit beyond their broad activity
Source-by-target conflict-event matrix Which targets appear in a state’s outgoing event pattern Which sources appear in a state’s incoming event pattern That source-target combination recurs beyond broad source and target involvement

These examples are ways to read the algebra, not labels produced automatically by a factor column. A fitted direction might mix ideology, region, institutional role, measurement practices, and other processes. The identified result is the pair surface and the probabilities built from it. Naming a direction requires validation outside the same matrix.

5.7 SVD Is the Scaffold, Not the Network Estimator

The phrase “similar to an SVD” is useful, but only if we keep the statistical differences visible.

Truncated SVD of a Fixed Matrix Latent Factor or AME Model
Minimizes squared reconstruction error Fits a likelihood or posterior appropriate to the outcome
Treats every supplied cell as an observed number Can define missing dyads, structural zeros, and a risk set
Has no measured predictors or additive actor effects Can include predictors, source effects, target effects, and reciprocity
Produces a point decomposition Can quantify uncertainty with MCMC or bootstrap refits
Decomposes the matrix you hand it Estimates a low-rank surface jointly with the other model pieces

For a binary probit AME model, the multiplicative surface helps explain a latent continuous relationship score whose threshold generates the observed zero or one. We do not run svd() once on the raw adjacency matrix and call the result an AME fit. ALS resembles repeated low-rank least-squares updates inside a larger estimation loop, while MCMC repeatedly samples coefficients, additive effects, factors, and variance quantities from their conditional distributions.

The connection is still fundamental. SVD shows why a low-rank product compresses repeated row-column patterns. The statistical model puts that product on the correct outcome scale, estimates it together with the other terms, and carries uncertainty. In both cases, the reconstructed surface UV𝖳UV^{\mathsf T} is more directly identified than one particular orientation or scaling of the raw factors.

5.8 The Multiplicative Term

The movie example gave us the operation. Replace “distance between positions” with an inner product of source and target profiles. Multiply matching coordinates and add the results. The sum is large and positive when the profiles line up in a relationship-enhancing way, near zero when they contribute little, and negative when their alignment lowers the fitted relationship score. A high source coordinate has no universal meaning because its contribution always depends on the corresponding target coordinate.

Start with the undirected case, since both of our simulated test networks are undirected. Each actor gets one vector uiu_i of length RR (the number of dimensions you choose), and the geometry enters as the symmetric eigenmodel (Hoff 2007):

γ(ui,uj)=ui⊤Λuj=∑r=1Rλrui,ruj,r,Λ=diag(λ1,…,λR), \gamma(u_i, u_j) = u_i^\top \Lambda\, u_j = \sum_{r=1}^{R} \lambda_r\, u_{i,r}\, u_{j,r}, \qquad \Lambda = \text{diag}(\lambda_1, \dots, \lambda_R),

Where the λr\lambda_r are (estimated) eigenvalue-like weights that can be positive or negative. Read the two moving parts separately:

  • The vectors carry stochastic equivalence. If ui≈uku_i \approx u_k then ii and kk relate to everyone else the same way; for any third actor jj, ui⊤Λuj≈uk⊤Λuju_i^\top \Lambda u_j \approx u_k^\top \Lambda u_j. Same role, similar rows, whether or not ii and kk tie to each other.
  • The sign of λr\lambda_r carries homophily vs heterophily. When λr>0\lambda_r > 0, actors that match on dimension rr (both high, or both low) get a positive boost: homophily, the distance-like behavior. When λr<0\lambda_r < 0, opposites attract on that dimension, which is exactly the disassortative structure the distance model choked on.

That sign is the symmetric model’s argument compressed into one symbol:

  • λr>0\lambda_r > 0: dimension rr contributes an assortative, positive semidefinite component. This is distance-like, but not literally the same likelihood as subtracting Euclidean distance.
  • λr<0\lambda_r < 0: dimension rr is disassortative: similar scores repel, opposite scores attract. A low-dimensional distance model cannot represent that component compactly, which is precisely what Section 4’s network exposed.

Hoff (2007) proves a careful version of containment. For a finite symmetric network, the eigenmodel can weakly represent the ordering of tie propensities generated by latent distance and latent class models, sometimes using a larger rank. Weak representation is not equality of likelihoods at the same dimension. One useful identity shows the connection for squared distance:

−‖ui−uj‖2=2ui⊤uj−‖ui‖2−‖uj‖2. -\lVert u_i-u_j\rVert^2 =2u_i^\top u_j-\lVert u_i\rVert^2-\lVert u_j\rVert^2.

The inner product carries the positive semidefinite geometry, while the two norm terms act like actor-specific sociality terms. The unrestricted eigenmodel also permits negative eigenvalues, which is why it can represent role patterns much more economically. Section 6 illustrates that distinction on two simulations; it does not prove the theorem from two k-means scores.

The movie example in Section 5.2 used an SVD to reconstruct a fixed matrix. The factor model estimates the same kind of low-rank product inside a probability model, together with measured predictors and actor effects. The application in Section 7 returns to the symmetric form first; the directed source-profile and target-profile version is shown alongside it so the distinction stays clear.

5.9 What Factor-Model Estimation Is Trying to Learn

The model must learn several things at once: the measured associations in β\beta, broad source and target differences in aa and bb, the residual relationship surface UV𝖳UV^{\mathsf T}, and any variance or reciprocity parameters. A good fit gives high probabilities to observed ties and low probabilities to observed non-ties without asking the low-rank surface to perform jobs already handled by the measured or additive terms.

For a probit AME fit, the sampler introduces a latent continuous score Yij*Y_{ij}^{*}. Scores for observed ties are sampled above zero, and scores for observed non-ties are sampled below zero. Conditional on those scores, the remaining updates become familiar regression and matrix problems:

  1. Update β\beta, aa, and bb given the current factor surface.
  2. Update the source profiles UU given the current target profiles VV.
  3. Update VV given the revised UU.
  4. Update reciprocity and variance parameters.
  5. Save the complete state and repeat.

The statistical target is

p(β,a,b,U,V,Σ∣Y,X)∝p(Y∣X,β,a,b,U,V,Σ)p(β,a,b,U,V,Σ). p(\beta,a,b,U,V,\Sigma\mid Y,X) \propto p(Y\mid X,\beta,a,b,U,V,\Sigma)\, p(\beta,a,b,U,V,\Sigma).

This is the joint posterior, not the squared error from one SVD of the adjacency matrix. The SVD explains the low-rank geometry and can help initialize or solve a conditional update. The probability model decides how observed zeros and ones contribute, includes measured predictors and actor effects, and propagates uncertainty through the other unknowns.

6 The Eigenvalue-Sign Demo: Containment Illustrated

The theorem is broader than any classroom simulation, so this section has a narrower job: show what positive and negative eigenvalues look like when we know the generating pattern. We have the distance-generated network from Section 2 and the disassortative-block network from Section 4. We fit the same symmetric AME specification (R=2R = 2, no measured covariates) to both. The fit still includes an additive actor-activity effect, so the eigen component represents the lower-rank pattern that remains after broad activity enters. We then read the signs of the posterior-mean Λ\Lambda. Each cold fit takes about 7 seconds, and both ship cached.

Ya <- Y;    diag(Ya) <- NA    # The distance-generated network (sec 2)
Yb <- Ydis; diag(Yb) <- NA    # The disassortative network (sec 4)

eig_assort <- cache_fit("eig_assort", ame(
  Y = Ya, family = "binary", symmetric = TRUE, R = 2, seed = 6886,
  nscan = 4000, burn = 2000, odens = 10,
  plot = FALSE, verbose = FALSE, print = FALSE
))
eig_dis <- cache_fit("eig_dis", ame(
  Y = Yb, family = "binary", symmetric = TRUE, R = 2, seed = 6886,
  nscan = 4000, burn = 2000, odens = 10,
  plot = FALSE, verbose = FALSE, print = FALSE
))

# Display dominant weight first: sort by decreasing |lambda|
lam_assort <- diag(eig_assort$L)[order(-abs(diag(eig_assort$L)))]
lam_dis    <- diag(eig_dis$L)[order(-abs(diag(eig_dis$L)))]
checkpoint(
  lambda_distance_like_network = lam_assort,
  lambda_disassortative_network = lam_dis
)
#> ------------------------------------------------------------------
#> CHECKPOINT: lambda_distance_like_network = 33.682,  3.737   |   lambda_disassortative_network = -51.048,  -1.986
#> ------------------------------------------------------------------

Read the dominant signs as a demonstration, not as draw-level inference. On the network that a distance actually generated, the posterior-mean multiplicative surface has Λ\Lambda = (+33.7, +3.7), with a dominant positive weight consistent with the assortative generating pattern. On the disassortative network, it has Λ\Lambda = (-51.0, -2.0), with a dominant negative weight: opposite-signed scores receive a positive contribution on that dimension, which compactly represents the structure the d=2d=2 distance fit missed. The smaller weight in either fit is not evidence for a second stable dimension without draw-level diagnostics. On observational data, even a stable sign describes the fitted surface rather than identifying the process that generated it.

And the positions are not just diagnostics; they carry the structure the LDM lost. Cluster the factor model’s UU on the disassortative network, where K-means on the LDM’s positions was a coin flip (Section 4):

km_fac  <- kmeans(eig_dis$U, 2, nstart = 50)
acc_fac <- max(mean(km_fac$cluster == truth), mean(km_fac$cluster == 3 - truth))
km_facA  <- kmeans(eig_assort$U, 2, nstart = 50)
acc_facA <- max(mean(km_facA$cluster == truth), mean(km_facA$cluster == 3 - truth))
checkpoint(
  kmeans_on_factor_U_disassortative = acc_fac,
  kmeans_on_factor_U_distance_like  = acc_facA,
  kmeans_on_LDM_positions_was       = acc_ldm
)
#> ------------------------------------------------------------------
#> CHECKPOINT: kmeans_on_factor_U_disassortative = 1   |   kmeans_on_factor_U_distance_like = 0.95   |   kmeans_on_LDM_positions_was = 0.55
#> ------------------------------------------------------------------

And to grade both models on both pitches, not just the factor model on the LDM’s failure case: fit the distance model on its own home turf, the network a distance actually generated. The cold fit takes about 5 seconds and ships cached:

ldm_home <- cache_fit("ldm_home", ergmm(
  network(Y, directed = FALSE) ~ euclidean(d = 2),
  control = ergmm.control(sample.size = 2000, interval = 10,
                          burnin = 5000),
  seed = 6886, verbose = FALSE
))
km_home  <- kmeans(ldm_home$mkl$Z, 2, nstart = 50)
acc_home <- max(mean(km_home$cluster == truth), mean(km_home$cluster == 3 - truth))
checkpoint(kmeans_on_LDM_positions_home_turf = acc_home)
#> ------------------------------------------------------------------
#> CHECKPOINT: kmeans_on_LDM_positions_home_turf = 0.95
#> ------------------------------------------------------------------
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(eig_assort$U, col = ifelse(truth == 1, "#18453B", "#7BBD00"), pch = 19,
     xlab = "u1", ylab = "u2", main = "distance-like net", asp = 1)
plot(eig_dis$U, col = ifelse(truth == 1, "#18453B", "#7BBD00"), pch = 19,
     xlab = "u1", ylab = "u2", main = "disassortative net", asp = 1)

The factor model’s estimated positions (U) for the two test networks, colored by TRUE group. Left: the distance-generated network, with positive weights. Right: the disassortative network, with separation on the dimension carrying the large negative weight. The chosen two-dimensional distance fit missed the latter pattern. Orientation is arbitrary; fitted pairwise structure is not.

Now the full two-by-two. On the distance model’s home turf, the two models tie: 95% for the LDM, 95% for the factor model. Off it, the factor model recovers the true blocks at 100% where the distance model’s positions scored 55%. This is an illustration consistent with Hoff’s containment result: the factor fit retained the coarse group signal on the distance-generated example and represented the disassortative example much better at the same rank. K-means accuracy on two draws is not a proof, and it does not establish that one family will dominate on every dataset.

ImportantThe One Sentence for This Section

In a symmetric eigenmodel, the sign of λr\lambda_r describes one component of the fitted low-rank surface: positive is assortative on that dimension; negative is disassortative. That flexibility is the eigenmodel’s advantage. When a network is assortative, the LDM’s restriction can be useful parsimony and its distance map is often easier to explain. When role mixing is plausible, compare it with an eigenmodel rather than assuming the geometry in advance. Prediction, posterior predictive checks, stability, and the research question should decide between them, not this one simulation.

7 From Broad State Involvement to Particular Directed Relationships

Day 9 ended with a longitudinal SRM that tracked whether ICEWS coded events from or toward each state across many partners. That model told us which states appeared in many above-threshold relationships, but not which bilateral relationships were driving the change. We now ask: when high-volume conflictual relationships change, does a state become more involved across many relationships, or do particular directed relationships change in their own ways?

For source state ii, target state jj, and year tt, the dynamic binary AME model is

Yijt*=αt+𝐱ijt𝖳𝛃+ait+bjt+𝐮it𝖳𝐯jt+eijt,Yijt=𝟙(Yijt*>0). Y^*_{ijt} =\alpha_t+\mathbf{x}_{ijt}^{\mathsf T}\boldsymbol\beta +a_{it}+b_{jt} +\mathbf{u}_{it}^{\mathsf T}\mathbf{v}_{jt} +e_{ijt}, \qquad Y_{ijt}=\mathbb{1}(Y^*_{ijt}>0).

Read the equation as a sequence of adjustments. The year intercept αt\alpha_t lets the whole network become busier or quieter. The measured part 𝐱ijt𝖳𝛃\mathbf{x}_{ijt}^{\mathsf T}\boldsymbol\beta contains the Polity-score gap and same-region indicator. The additive effects aita_{it} and bjtb_{jt} track whether events are coded from state ii or toward state jj across many partners. The inner product 𝐮it𝖳𝐯jt\mathbf{u}_{it}^{\mathsf T}\mathbf{v}_{jt} is the new part. It asks whether this exact directed pair is more or less likely to cross the cutoff than those broad state patterns would suggest.

ImportantThe Main Distinction

The additive terms answer, “Are high-volume events coded from or toward this state across many partners?” The multiplicative term answers, “Does this exact directed state pair have a different history than those broad patterns would suggest?” A dynamic multiplicative term lets that bilateral answer change across years.

7.1 Return to the ICEWS Panel

We use the same Integrated Crisis Early Warning System panel as Day 9 so the change in the model is easy to see. ICEWS records coded events with a source actor and a target actor (Boschee et al. 2015). Our directed outcome equals one when the data contain more than 20 material-conflict events from state ii toward state jj in a year. Crossing that threshold means the relationship generated a high volume of coded events. It does not by itself identify military initiation, victimization, conflict severity, or a causal effect. Event volume can also reflect sustained activity, media visibility, and source coverage.

The classroom panel contains 18 states from 2002 through 2014. They are the highest-volume ICEWS source states among those with complete Polity and GDP coverage in all 13 years. This purposive case set gives the model enough above-threshold relationships and enough changes to learn from. It is not a probability sample of the international system.

icews <- readRDS("data/icews.rds")
monadic_coverage <- unique(
  icews[, c("i", "year", "i_polity2", "i_log_gdp")]
)
coverage_counts <- aggregate(
  cbind(
    polity_years = as.integer(!is.na(i_polity2)),
    gdp_years = as.integer(!is.na(i_log_gdp))
  ) ~ i,
  monadic_coverage,
  sum
)
eligible_states <- coverage_counts$i[
  coverage_counts$polity_years == 13 &
    coverage_counts$gdp_years == 13
]
conflict_totals <- aggregate(matlConf ~ i, icews, sum)
conflict_totals <- conflict_totals[
  conflict_totals$i %in% eligible_states,
]
states <- head(
  conflict_totals$i[
    order(conflict_totals$matlConf, decreasing = TRUE)
  ],
  18
)

icews_sub <- icews[
  icews$i %in% states & icews$j %in% states,
]
icews_sub$high_event_volume <-
  as.integer(icews_sub$matlConf > 20)
icews_sub$polity_gap <-
  abs(icews_sub$i_polity2 - icews_sub$j_polity2) / 20
icews_sub$same_region <-
  as.integer(icews_sub$i_region == icews_sub$j_region)
icews_dyad <- icews_sub[, c(
  "i", "j", "year", "matlConf", "high_event_volume",
  "polity_gap", "same_region"
)]
names(icews_dyad)[1:4] <- c(
  "source", "target", "year", "material_conflict_events"
)
icews_dyad$year <- as.integer(icews_dyad$year)
years <- sort(unique(icews_dyad$year))

checkpoint(
  state_pair_years = nrow(icews_dyad),
  years = length(years),
  states = length(states),
  high_volume_rate =
    round(mean(icews_dyad$high_event_volume), 3)
)
#> ------------------------------------------------------------------
#> CHECKPOINT: state_pair_years = 3978   |   years = 13   |   states = 18   |   high_volume_rate = 0.218
#> ------------------------------------------------------------------
states
#>  [1] "United States"             "Israel"                   
#>  [3] "Russian Federation"        "Pakistan"                 
#>  [5] "India"                     "United Kingdom"           
#>  [7] "Iran, Islamic Republic Of" "China"                    
#>  [9] "France"                    "Japan"                    
#> [11] "Syrian Arab Republic"      "Thailand"                 
#> [13] "Australia"                 "Egypt"                    
#> [15] "Korea, Republic Of"        "Spain"                    
#> [17] "Indonesia"                 "Germany"

The unit is an ordered state-pair-year. Reversing a pair produces a different observation because a high event volume from Iran toward Syria is not the same record as a high event volume from Syria toward Iran. The two measured predictors have simple units. polity_gap is the absolute difference between the two states’ Polity scores divided by 20, so it runs from 0 to 1. same_region equals one when the states share a World Bank region.

The outcome changes 462 times across adjacent years: 244 relationships cross above the threshold and 218 fall below it. That movement is the practical reason this example can support dynamic effects. A model cannot learn meaningful movement merely because we give it time-varying parameters. The observed network has to move too.

icews_change <- do.call(
  rbind,
  lapply(years[-1], function(yy) {
    now <- icews_dyad[
      icews_dyad$year == yy,
      c("source", "target", "high_event_volume")
    ]
    before <- icews_dyad[
      icews_dyad$year == yy - 1,
      c("source", "target", "high_event_volume")
    ]
    names(now)[3] <- "now"
    names(before)[3] <- "before"
    both <- merge(
      before,
      now,
      by = c("source", "target"),
      sort = FALSE
    )
    data.frame(
      year = yy,
      changes = sum(both$now != both$before),
      onsets = sum(both$now == 1 & both$before == 0),
      dissolutions = sum(both$now == 0 & both$before == 1)
    )
  })
)
icews_change
year changes onsets dissolutions
2003 32 25 7
2004 29 11 18
2005 23 15 8
2006 36 22 14
2007 43 20 23
2008 41 18 23
2009 37 18 19
2010 47 32 15
2011 57 22 35
2012 44 28 16
2013 33 12 21
2014 40 21 19

7.2 Let netify Build the Matrices

The model needs one outcome matrix and one predictor array for each year, with exactly the same state order in every object. netify() handles that bookkeeping and to_lame() produces the lists that lame() expects.

icews_net <- netify(
  input = icews_dyad,
  actor1 = "source",
  actor2 = "target",
  time = "year",
  symmetric = FALSE,
  weight = "high_event_volume",
  dyad_vars = c("polity_gap", "same_region"),
  dyad_vars_symmetric = c(TRUE, TRUE),
  missing_to_zero = FALSE
)
stopifnot(all(unlist(validate_netify(icews_net))))

icews_lame <- to_lame(
  icews_net,
  lame = TRUE,
  family = "binary",
  fit_method = "als"
)
Y_icews <- icews_lame$Y
X_icews <- icews_lame$Xdyad

checkpoint(
  outcome_matrices = length(Y_icews),
  actors_per_matrix = nrow(Y_icews[[1]]),
  predictors = dim(X_icews[[1]])[3]
)
#> ------------------------------------------------------------------
#> CHECKPOINT: outcome_matrices = 13   |   actors_per_matrix = 18   |   predictors = 2
#> ------------------------------------------------------------------

validate_netify() checks the actor labels, time slices, direction, missingness, and agreement among the stored representations before estimation begins. The diagonal is NA because a state cannot have a relationship with itself. Every off-diagonal zero is observed because the source table contains every ordered pair in every year. We therefore use missing_to_zero = FALSE: no unrecorded pair is silently converted into data.

7.3 What Is Dynamic in This Model?

Three parts move over time, and each has a different job.

Dynamic Part What Changes Question It Answers
αt\alpha_t The baseline intercept Did above-threshold relationships become more or less common across the whole 18-state panel this year?
aita_{it} and bjtb_{jt} Broad source-side and target-side effects Were events coded from or toward this state across many different partners?
𝐮it\mathbf{u}_{it} and 𝐯jt\mathbf{v}_{jt} Directed relationship profiles After accounting for those broad patterns, did particular directed state relationships change?

The paths are smoothed across adjacent years. In the ALS implementation, the model does not estimate 13 unrelated versions of every effect. It penalizes large jumps from one year to the next, which lets neighboring years share information while still allowing sustained change.

NoteWhat dynamic_beta = "intercept" Does

dynamic_beta = "intercept" lets the baseline intercept change by year. dynamic_beta_kind = "rw1" uses a first-order random-walk penalty, which treats a smooth drift as more plausible than a path that repeatedly jumps back and forth. The Polity-gap and same-region coefficients remain fixed across years. That choice keeps the application focused: the volume of the system, broad state involvement, and partner profiles can change, while the two measured associations are summarized over the full panel.

7.4 How Dynamic ALS Estimates the Model

The dynamic model is trying to reconstruct every ordered state-pair-year while keeping neighboring years connected. It rewards a close match between the model’s working relationship scores and the data, but it also penalizes paths that jump sharply from one year to the next. This is how it separates sustained movement from isolated noise without forcing every year to look identical.

For the binary model, an outer iteratively reweighted least-squares step converts the current fitted probabilities into a working response zijtz_{ijt} and weights wijtw_{ijt}. Conditional on that working response, ALS approximately minimizes an objective of the form

Q=∑t,i≠jwijt[zijt−αt−xijt𝖳β−ait−bjt−uit𝖳vjt]2+λab∑t(‖at−at−1‖2+‖bt−bt−1‖2)+λuv∑t(‖Ut−Ut−1‖F2+‖Vt−Vt−1‖F2)+λα∑t(αt−αt−1)2. \begin{aligned} Q={}&\sum_{t,i\ne j}w_{ijt} \left[ z_{ijt} -\alpha_t-x_{ijt}^{\mathsf T}\beta-a_{it}-b_{jt} -u_{it}^{\mathsf T}v_{jt} \right]^2\\ &+\lambda_{ab}\sum_t\left(\lVert a_t-a_{t-1}\rVert^2+\lVert b_t-b_{t-1}\rVert^2\right) +\lambda_{uv}\sum_t\left(\lVert U_t-U_{t-1}\rVert_F^2+\lVert V_t-V_{t-1}\rVert_F^2\right) +\lambda_{\alpha}\sum_t(\alpha_t-\alpha_{t-1})^2. \end{aligned}

The first line is the reconstruction job. It asks the baseline, measured predictors, broad actor paths, and factor surface to reproduce the current working scores. The remaining lines are the smoothing job. They charge the model for large year-to-year changes in the actor effects, factor profiles, and intercept. The λ\lambda values control how expensive those changes are.

  1. Start with provisional fitted probabilities for every ordered pair in every year.
  2. Construct the working responses and weights from those probabilities.
  3. Update the yearly baseline and the two measured coefficients while temporarily holding the actor and factor paths fixed.
  4. Update each state’s broad source and target paths while temporarily holding the other pieces fixed.
  5. Update the changing source and target factor profiles, using low-rank matrix updates to revise the pair-specific surface.
  6. Recalculate the fitted probabilities and repeat until the objective and fitted surface barely change.

This is fast penalized point estimation. For a binary model, it solves a sequence of weighted least-squares approximations rather than directly maximizing the exact Bernoulli likelihood in one step. With rank greater than zero, the problem is nonconvex, so different starts can reach different local solutions. The rho_ab, rho_uv, and rho_beta settings determine the smoothing penalties. They are fixed tuning values, not persistence parameters estimated from these data. Bootstrap refits describe how the selected fitting procedure varies across replicated data, while MCMC would instead represent a posterior distribution under a fully specified Bayesian model.

7.5 Fit the Dynamic SRM and Dynamic AME

The first model is the Day 9 structure with a changing intercept and additive effects. The second adds two changing multiplicative dimensions.

# The downloaded folder contains these fitted objects.
# Run _build_icews_dynamic_cache.R only when you intend to refit them.
load("cache/icews_dynamic_srm_r0_v1.rda")
load("cache/icews_dynamic_ame_r2_boot100_v1.rda")

dynamic_fit_summary <- data.frame(
  model = c("Dynamic SRM, rank 0", "Dynamic AME, rank 2"),
  converged = c(fit_r0$converged, fit_r2$converged),
  iterations = c(fit_r0$iterations, fit_r2$iterations),
  deviance = round(c(fit_r0$deviance, fit_r2$deviance), 1)
)
dynamic_fit_summary
model converged iterations deviance
Dynamic SRM, rank 0 TRUE 28 3127.9
Dynamic AME, rank 2 TRUE 48 2085.1

The call that produced the rank-two fit is shown below. It is not evaluated during rendering because the saved fit and its 100 bootstrap refits are included with the handout.

fit_r2 <- lame(
  Y = Y_icews,
  Xdyad = X_icews,
  family = "binary",
  symmetric = FALSE,
  R = 2,
  dynamic_ab = TRUE,
  dynamic_uv = TRUE,
  dynamic_beta = "intercept",
  dynamic_beta_kind = "rw1",
  method = "als",
  als_max_iter = 1000,
  als_tol = 1e-5,
  als_stability = "validation",
  bootstrap = 100,
  bootstrap_type = "parametric",
  bootstrap_seed = 6886,
  seed = 6886,
  verbose = FALSE
)

dynamic_ab = TRUE gives each state a changing broad source and target path. dynamic_uv = TRUE gives each state changing source and target factor profiles. als_stability = "validation" repeats the optimization from four jittered starts. All four starts converged, their fitted probability surfaces correlated at essentially 1.00, and the largest root-mean-squared difference from the main surface was below 0.0001. That is strong evidence that this point solution is not an accident of one starting value.

7.6 Read Broad State Movement With ab_plot()

The additive paths describe broad involvement across many partners. They do not tell us which partner accounts for the movement.

fit_r2_display <- fit_r2
display_names <- fit_r2$row_names
display_names[display_names == "Iran, Islamic Republic Of"] <- "Iran"
display_names[display_names == "Syrian Arab Republic"] <- "Syria"
display_names[display_names == "Korea, Republic Of"] <- "South Korea"
rownames(fit_r2_display$a_dynamic) <- display_names
rownames(fit_r2_display$b_dynamic) <- display_names
dimnames(fit_r2_display$U)[[1]] <- display_names
dimnames(fit_r2_display$V)[[1]] <- display_names
fit_r2_display$row_names <- display_names
fit_r2_display$col_names <- display_names
focus_states <- c(
  "United States", "Syria",
  "Iran", "Russian Federation"
)
state_colors <- c(
  "Syria" = "#18453B",
  "Russian Federation" = "#008208",
  "Iran" = "#7BBD00",
  "United States" = "#535054"
)
p_source <- lame::ab_plot(
  fit_r2_display,
  effect = "sender",
  plot_type = "trajectory",
  show_actors = focus_states,
  title = "Broad Source-Side Activity"
) +
  aes(linetype = actor, shape = actor) +
  scale_color_manual(values = state_colors, name = NULL) +
  scale_linetype_discrete(name = NULL) +
  scale_shape_discrete(name = NULL) +
  labs(x = "Year", y = "Source-side adjustment")
p_target <- lame::ab_plot(
  fit_r2_display,
  effect = "receiver",
  plot_type = "trajectory",
  show_actors = focus_states,
  title = "Broad Target-Side Exposure"
) +
  aes(linetype = actor, shape = actor) +
  scale_color_manual(values = state_colors, name = NULL) +
  scale_linetype_discrete(name = NULL) +
  scale_shape_discrete(name = NULL) +
  labs(x = "Year", y = "Target-side adjustment")
p_source / p_target

Two stacked line charts show broad source-side and target-side adjustments for Iran, Russia, Syria, and the United States from 2002 to 2014. The states are distinguished by color, line type, and point shape.

Changing broad source-side activity and target-side exposure in the dynamic rank-two fit. These are additive effects, so each path summarizes involvement across all possible partners.

Read one line at a time. A rising source-side path means that the state appears as the coded source with more partners than the measured variables and pair-specific pattern would otherwise predict. A rising target-side path means that the state appears as the coded target with more partners. Neither line says that a state started more conflicts or suffered greater harm. The outcome is high coded event volume, and each additive path summarizes all possible partners.

The package’s ab_plot() function is doing the right job here because it reads the dynamic additive-effect arrays directly. The bootstrap object also contains refitted aita_{it} and bjtb_{jt} paths. Those refits describe uncertainty under the fitted model and the chosen smoothing specification. They do not add uncertainty about whether the ICEWS threshold or case selection was the right measurement choice.

7.7 A Bird’s-Eye View Before We Follow Named Pairs

The multiplicative paths are the relationship-specific part. We align the yearly source and target profiles jointly, verify that their fitted inner products do not change, and then use the package’s uv_plot() function to inspect selected years. The optional identification box explains why this display step is necessary.

fit_r2_aligned <- align_directed_uv(fit_r2_display)
uv_snapshot <- function(tt) {
  p <- lame::uv_plot(
    fit_r2_aligned,
    layout = "biplot",
    plot_type = "snapshot",
    time_point = tt,
    label.nodes = FALSE,
    show_arrows = FALSE,
    title = years[tt]
  )
  focus_data <- p$data[p$data$name %in% c("Iran", "Syria"), ]
  focus_data$label <- paste0(
    focus_data$name,
    ifelse(focus_data$type == "Sender", " (source)", " (target)")
  )
  p +
    ggrepel::geom_text_repel(
      data = focus_data,
      aes(label = label),
      size = 3,
      max.overlaps = Inf,
      show.legend = FALSE
    ) +
    scale_color_manual(
      values = c(Sender = "#18453B", Receiver = "#7BBD00"),
      labels = c(
        Sender = "Source profile",
        Receiver = "Target profile"
      ),
      name = NULL
    ) +
    scale_shape_manual(
      values = c(Sender = 16, Receiver = 1),
      labels = c(
        Sender = "Source profile",
        Receiver = "Target profile"
      ),
      name = NULL
    )
}
wrap_plots(
  lapply(c(1, 13), uv_snapshot),
  ncol = 2,
  guides = "collect"
) & theme(legend.position = "bottom")

Two biplots compare source and target profiles in 2002 and 2014. Circles and colors distinguish source from target profiles, while labels identify Iran and Syria. The axes are unlabeled because their orientation has no unique substantive meaning.

Aligned source and target profiles at the beginning and end of the panel. Labels identify the Iran and Syria profiles. The axes themselves do not have substantive names.

The map is a compact picture of many fitted pair adjustments, not the final interpretation. When a source profile and a target profile point in similar directions, their multiplicative term raises that ordered pair’s fitted score. When they point in opposite directions, it lowers the score. The horizontal and vertical axes can rotate or reflect without changing any fitted pair adjustment, so we do not name them “democracy,” “region,” or “conflict.” For a substantive claim, return to a named pair and its fitted probability.

WarningDo Not Read This as Literal Geopolitical Movement

If Syria’s fitted profile moves toward Iran’s target profile, the model is not saying that the countries literally became closer in a geographic or diplomatic space. It is saying that the pattern of above-threshold source-target relationships involving those states changed in a way the year baseline and broad additive paths could not capture.

7.8 Choose Rank With More Than In-Sample Fit

Adding dimensions will almost always improve reconstruction of the observations used for estimation. We therefore hide complete unordered state-pair histories: both directions and all 13 years for a pair are withheld together. The model can learn each state’s profile from its relationships with other states, but it cannot see the held pair in another year or in reverse.

load("cache/icews_dynamic_rank_cv_v1.rda")
rank_cv_summary <- rank_cv$summary
rank_cv_summary$auroc <- round(rank_cv_summary$auroc, 3)
rank_cv_summary$auprc <- round(rank_cv_summary$auprc, 3)
rank_cv_summary$brier <- round(rank_cv_summary$brier, 3)
rank_cv_summary$logloss <- round(rank_cv_summary$logloss, 3)
names(rank_cv_summary)[
  names(rank_cv_summary) == "auprc"
] <- "precision_recall_auc"
rank_cv_summary
rank n_eval auroc precision_recall_auc brier logloss converged
0 795.6 0.825 0.676 0.118 0.414 5
1 795.6 0.837 0.694 0.113 0.400 5
2 795.6 0.837 0.695 0.113 0.404 5
3 795.6 0.833 0.686 0.116 0.423 5
4 795.6 0.839 0.687 0.115 0.410 5

Ranks one and two are essentially tied on held-pair prediction. Rank two has a precision-recall AUC of about 0.695 and a Brier score of about 0.113. Ranks three and four do not produce a consistent predictive gain. We keep rank two because it also improves the temporal reproduction check and gives a two-dimensional profile display that we can inspect. This is a judgment across prediction, temporal fit, and interpretability, not a declaration that rank two is the uniquely true dimension.

7.9 Put the Result Back Into State-Pair-Years

The easiest way to understand the model is to follow specific ordered pairs. For each year, we show the observed threshold outcome, the fitted probability, and the multiplicative contribution.

p_r2 <- fit_r2$fitted
uv_r2 <- lapply(
  seq_along(years),
  function(tt) fit_r2$Oarr[, , tt]
)

pair_path <- function(source, target, source_label, target_label) {
  ii <- match(source, rownames(Y_icews[[1]]))
  jj <- match(target, colnames(Y_icews[[1]]))
  data.frame(
    source = source_label,
    target = target_label,
    year = years,
    observed = vapply(
      Y_icews,
      function(M) M[ii, jj],
      numeric(1)
    ),
    fitted_probability = vapply(
      p_r2,
      function(M) M[ii, jj],
      numeric(1)
    ),
    lower = vapply(
      seq_along(years),
      function(tt) {
        quantile(
          vapply(
            fit_r2$bootstrap_dynamic$refits,
            function(refit) {
              refit$fitted[[tt]][ii, jj]
            },
            numeric(1)
          ),
          .025
        )
      },
      numeric(1)
    ),
    upper = vapply(
      seq_along(years),
      function(tt) {
        quantile(
          vapply(
            fit_r2$bootstrap_dynamic$refits,
            function(refit) {
              refit$fitted[[tt]][ii, jj]
            },
            numeric(1)
          ),
          .975
        )
      },
      numeric(1)
    ),
    multiplicative_contribution = vapply(
      uv_r2,
      function(M) M[ii, jj],
      numeric(1)
    )
  )
}

pair_paths <- rbind(
  pair_path(
    "Iran, Islamic Republic Of", "Syrian Arab Republic",
    "Iran", "Syria"
  ),
  pair_path(
    "Syrian Arab Republic", "Iran, Islamic Republic Of",
    "Syria", "Iran"
  ),
  pair_path(
    "United States", "Syrian Arab Republic",
    "United States", "Syria"
  ),
  pair_path(
    "Syrian Arab Republic", "United States",
    "Syria", "United States"
  )
)
pair_paths$pair <- paste(
  pair_paths$source,
  "toward",
  pair_paths$target
)

ggplot(pair_paths, aes(year, fitted_probability)) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    fill = "#7BBD00",
    alpha = .18
  ) +
  geom_line(color = "#18453B", linewidth = 1) +
  geom_point(
    aes(shape = factor(observed)),
    size = 2.3
  ) +
  facet_wrap(~ pair, ncol = 2) +
  scale_shape_manual(
    values = c("0" = 1, "1" = 19),
    name = "Observed threshold"
  ) +
  scale_y_continuous(
    labels = scales::label_percent(),
    limits = c(0, 1)
  ) +
  labs(
    x = NULL,
    y = "Estimated chance of more than 20 coded events",
    shape = "Observed threshold"
  ) +
  theme(
    legend.position = "bottom"
  )

Four small line charts show Iran toward Syria, Syria toward Iran, the United States toward Syria, and Syria toward the United States from 2002 to 2014. Iran-Syria probabilities rise sharply around 2011 and remain high. United States-Syria probabilities are high throughout. Ribbons show bootstrap uncertainty, and point shape shows the observed binary outcome.

Estimated probability paths for four ordered relationships. Ribbons are 95% intervals from 100 parametric-bootstrap refits, and filled points mark years that actually crossed the 20-event threshold.

Iran toward Syria stays below the threshold through 2011, then crosses it in 2012, 2013, and 2014. Its fitted probability rises from about 2% in 2002 to 92% in 2012, with a 95% bootstrap interval of roughly 57% to 98% in 2012. Syria toward Iran crosses one year earlier. Its fitted probability rises from about 1% in 2002 to 69% in 2011 and remains elevated through 2014. The two directions are connected, but the model does not force them to have the same timing or probability.

The United States toward Syria is above the threshold in nearly every year after 2002, and the model assigns it a high probability throughout. Syria toward the United States is less regular early in the panel and becomes more consistently above the threshold later. This is what the multiplicative term adds: states can both appear with many partners, yet a particular ordered relationship can have its own timing.

These paths show that the rise involving Iran and Syria is not simply a uniform increase across every relationship involving either state. The timing coincides with the escalation of the Syrian conflict, but the coded events and fitted model do not establish why the records changed or measure battlefield severity. A useful next step would connect the probability paths to dated events and alternative measurements, then ask whether the interpretation survives those checks.

7.10 Check the Annual Networks and Their Changes

We now simulate 500 complete panels while holding the fitted ALS solution fixed. For each year, we compare network density, variation in source and target rates, reciprocity, and two normalized triadic-dependence statistics. We also count how many ordered pairs change status from one simulated year to the next.

load("cache/icews_dynamic_gof_v1.rda")
gof_coverage <- data.frame(
  statistic = c(
    "Density",
    "Source-rate variation",
    "Target-rate variation",
    "Reciprocity",
    "Cyclic triadic dependence",
    "Transitive triadic dependence",
    "Adjacent-year changes"
  ),
  dynamic_srm = c(gof$r0$coverage, gof$r0$change_cover),
  dynamic_ame_rank2 = c(gof$r2$coverage, gof$r2$change_cover),
  periods = c(rep(13, 6), 12)
)
gof_coverage
statistic dynamic_srm dynamic_ame_rank2 periods
density Density 13 13 13
sd.rowmean Source-rate variation 13 13 13
sd.colmean Target-rate variation 13 13 13
dyad.dep Reciprocity 5 6 13
cycle.dep Cyclic triadic dependence 11 12 13
trans.dep Transitive triadic dependence 10 12 13
Adjacent-year changes 2 8 12

The dynamic rank-two model covers the observed density and source- and target-rate variation in all 13 years. It covers the cyclic and transitive summaries in 12 of 13 years. Reciprocity is the main annual weakness, with the observed value inside the simulation interval in 6 of 13 years. The multiplicative surface captures much of the recurring relational structure, but it does not perfectly reproduce the connection between opposite directions of a pair.

The transition check is especially important because this is a dynamic application. The additive model covers the observed number of changes in only 2 of 12 adjacent-year transitions. Rank two covers 8 of 12. The factor paths therefore make year-to-year redraws much more realistic, but four transitions remain outside the simulation interval.

transition_plot_data <- do.call(
  rbind,
  lapply(c("r0", "r2"), function(model_name) {
    change_matrix <- gof[[model_name]]$changes
    data.frame(
      year = years[-1],
      model = ifelse(
        model_name == "r0",
        "Dynamic SRM, rank 0",
        "Dynamic AME, rank 2"
      ),
      observed = icews_change$changes,
      median = apply(change_matrix, 1, median),
      lower = apply(change_matrix, 1, quantile, 0.025),
      upper = apply(change_matrix, 1, quantile, 0.975)
    )
  })
)
transition_plot_data$model <- factor(
  transition_plot_data$model,
  levels = c("Dynamic SRM, rank 0", "Dynamic AME, rank 2")
)

ggplot(
  transition_plot_data,
  aes(year, median, color = model, fill = model)
) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    alpha = 0.14,
    color = NA
  ) +
  geom_line(linewidth = 0.9) +
  geom_point(
    aes(y = observed),
    color = "#111111",
    size = 2
  ) +
  facet_wrap(~ model, ncol = 1) +
  scale_color_manual(values = c(
    "Dynamic SRM, rank 0" = "#535054",
    "Dynamic AME, rank 2" = "#18453B"
  )) +
  scale_fill_manual(values = c(
    "Dynamic SRM, rank 0" = "#535054",
    "Dynamic AME, rank 2" = "#18453B"
  )) +
  labs(
    x = NULL,
    y = "Ordered pairs that changed threshold status",
    color = NULL,
    fill = NULL
  ) +
  guides(color = "none", fill = "none")

Two stacked panels compare observed year-to-year tie changes with simulated intervals for the dynamic SRM and rank-two AME. The rank-two intervals contain eight of twelve observed changes, compared with two of twelve for the SRM.

Observed year-to-year changes and point-fit simulation intervals. Rank two closes much of the additive model’s temporal gap, although four observed transitions remain outside its 95% interval.

The black points are the observed counts. The lines are the typical counts from simulated panels, and the shaded areas contain the middle 95% of simulated counts. If a black point falls outside the shaded area, the model is generating too much or too little change for that transition. The rank-two panel is visibly closer to the observations than rank zero, but it is not perfect.

ImportantWhat the GOF Checks Show

The dynamic rank-two model gets the annual network shape largely right and makes year-to-year redraws much more realistic than the dynamic SRM. It still misses reciprocity in several years and does not reproduce every transition. Those failures limit the interpretation. They do not erase the model’s useful improvement.

7.11 What the Bootstrap Adds

The 100 parametric bootstrap refits simulate a panel from the fitted model and estimate the same dynamic AME again. The spread across successful refits tells us how much the point estimates would vary if the fitted data-generating process were repeated. This is useful for intervals around coefficient paths, additive paths, and derived probabilities.

It is still model-based uncertainty. The bootstrap assumes the fitted rank, smoothing setup, binary threshold, case selection, and measurement process are correct. It does not turn the Polity-gap or same-region coefficient into a causal effect, and it does not account for alternative ways of coding the outcome.

7.12 What to Do After a Dynamic AME Runs

  1. Check the optimization. Confirm convergence, inspect the objective path, and compare multiple starts.
  2. Check the prediction target. Hold complete pair histories together when the research goal involves relationships that were absent from estimation.
  3. Check annual structure and transitions. A model can reconstruct each year while still producing unrealistic movement between years.
  4. Interpret broad paths and partner paths separately. Use aita_{it} and bjtb_{jt} for broad involvement, then use 𝐮it𝖳𝐯jt\mathbf{u}_{it}^{\mathsf T}\mathbf{v}_{jt} and fitted probabilities for particular ordered relationships.
  5. Use the package plots as entry points. ab_plot() displays broad actor paths, while uv_plot() and latent_positions() help inspect multiplicative profiles. The fitted inner products and probabilities carry the relationship interpretation.
  6. Return to the historical record and measurement process. Use external evidence to explain a path, test nearby outcome definitions, and say clearly what the data cannot establish.
TipWhat the Coded Event Records Show

Across these 18 states, the model distinguishes a general rise in how often a state appears in above-threshold relationships from a rise concentrated in one directed pair. The estimated chance for Iran toward Syria rises from about 2% in 2002 to 92% in 2012. Syria toward Iran rises on a different timetable, reaching about 69% in 2011. These changes occur around the escalation of the Syrian conflict, but the event records and model do not establish why they occurred or measure battlefield severity. The richer model reproduces much more of the observed year-to-year change than the broad state-level model, although it still misses reciprocity in several years and four of the twelve annual transitions.

8 What Latent Models Did Not Buy Us

Latent models can give us a better representation of connected observations, better prediction, and a useful summary of recurring relationships. They do not, by themselves, tell us what would happen under an intervention. That boundary matters because the model can represent an omitted pattern without separating the part of that pattern that is also associated with a measured predictor.

Return to the ICEWS example. Suppose states with a particular regime profile also differ in an unmeasured, time-varying feature such as international news visibility. That feature can raise the number of coded events involving the state and can be correlated with the Polity-gap predictor. A random-effects model that omits it can give the measured predictor credit for both patterns. Additive and multiplicative effects can represent the leftover relational pattern, but the same event records may still be unable to tell us how much credit belongs to regime difference and how much belongs to the omitted feature.

Minhas et al. (2022) make this boundary especially clear in simulation. When an omitted low-rank term WW is independent of the observed predictor XX, AME can absorb the patterned residual dependence and recover the measured association much better than an independence model. When WW and XX are correlated, the network alone does not tell the model how much of their shared pattern belongs to XX and how much belongs to WW. Even a well-fitting AME coefficient can then remain biased.

Here is the basic version in notation. Suppose the latent outcome is

Y*=βX+γW+e, Y^*=\beta X+\gamma W+e,

where WW is an omitted relational pattern. If that pattern is related to the measured predictor so that

W=αX+z, W=\alpha X+z,

then substituting the second equation into the first gives

Y*=(β+γα)X+γz+e. Y^*=(\beta+\gamma\alpha)X+\gamma z+e.

An additive model that leaves out WW can therefore give XX credit for both its own association β\beta and the portion of the relational pattern that moves with it, γα\gamma\alpha. Adding an AME surface can represent WW, but the same network still may not contain enough information to separate β\beta from the part of WW that overlaps XX. Put simply, the factor model can notice that some relationships are systematically unusual. It cannot automatically decide whether a measured predictor caused that pattern or merely travels with an unmeasured cause.

This is the same random-effects issue we met in the SRM. The usual random-effects SRM assumes that the additive actor effects are unrelated to the observed predictors after conditioning on the model. If states with large unobserved source-side or target-side tendencies also systematically have particular observed characteristics, a coefficient can mix the predictor association with those state tendencies. A correlated-random-effects or Mundlak specification addresses one part of that problem. For a time-varying actor predictor xitx_{it}, include both the actor’s time average x‾i\bar{x}_i and the deviation xit−x‾ix_{it}-\bar{x}_i in the SRM. The deviation coefficient then compares a state with itself across times when its predictor differs from its usual level, while the average captures stable differences between states. This can be done in lame() by constructing those two variables before netify() and passing them as actor or dyadic predictors. It does not automatically solve dyad-level confounding, reverse causation, or causal identification.

ImportantDependence Adjustment Is Not Confounding Control

Adding rank does not mechanically remove confounding. It represents residual relational structure under a specified model. Causal identification still requires a defensible design or assumptions about assignment, timing, overlap, interference, and the relationship between measured and unmeasured causes.

The usual omitted-variable formula is a useful directional heuristic. Suppose, only for this table, that a larger omitted quantity WW raises the outcome and that we are thinking about a linear model after the other measured predictors have been removed.

Relationship Between the Residualized XX and Omitted WW Likely Direction of the Simple Linear Bias Plain-Language Reading
Positive covariance Upward High values of XX tend to occur where the omitted outcome-raising quantity is also high, so the fitted slope can give XX credit for both.
Negative covariance Downward High values of XX tend to occur where the omitted outcome-raising quantity is low, so the fitted slope can hide part of the association with XX.
Zero covariance No simple slope bias from this omitted term under the linear assumptions The omitted pattern may still create dependent errors and bad uncertainty estimates, but this particular covariance channel does not shift the slope.

If WW lowers the outcome, reverse the upward and downward labels. In a probit AME, an LDM with jointly estimated positions, or a blockmodel with estimated memberships, this table is intuition rather than an exact bias formula. It helps us ask the right question: could the measured predictor and the unmeasured relational pattern be carrying the same information?

8.1 What β\beta Means in Each Model

Across these models, β\beta lives in a conditional linear predictor. The thing held fixed changes, so the substantive interpretation changes too.

Model What the Coefficient Conditions On A Defensible Reading
SRM Source effects, target effects, and reciprocal dependence in the observation model Association on the link scale among observations with the same modeled source-side and target-side tendencies
SBM The relevant block-pair baseline, or the mixed-membership average over block pairs Block-adjusted association, provided the predictor still varies within the comparisons doing the estimation
LDM Fitted latent distance and any activity terms in the specification Association at a fixed fitted distance, not a treatment effect
AME Additive actor effects and the fitted multiplicative surface Association at a fixed modeled source, target, and recurring relationship profile

A probit or logit coefficient is not a percentage-point change. For a binary predictor, the conditional probability contrast for one pair is

g−1(ηij+β)−g−1(ηij), g^{-1}(\eta_{ij}+\beta)-g^{-1}(\eta_{ij}),

which depends on the pair’s baseline score ηij\eta_{ij}. That is why Section 7.9 translates the dynamic fit into state-pair-year probabilities while holding the other fitted pieces fixed. An average fitted contrast is still a model-based comparison, not an observed before-and-after change or a causal effect.

Every conditional comparison also needs support. If a predictor is almost fixed within the source-target profiles doing the estimation, the model has little information for comparing otherwise similar observations at different predictor values. Adding a latent surface cannot create a comparison the data do not contain. If changing XX would itself change the latent profile, “hold the profile fixed” describes the fitted conditional association rather than the total effect of changing XX.

Three ideas are easy to blur:

  1. A random-effects assumption. A standard latent model often assumes that latent effects have zero conditional mean given XX, or are independent of XX. A different relationship can be modeled, but it must come from an explicit specification or additional information.
  2. An omitted-variable calculation. The covariance table above provides linear intuition. It is not an exact formula for nonlinear latent models.
  3. A causal identification claim. Causal language additionally needs a defensible assignment and timing story, consistency, overlap, and an interference assumption that matches the network.

A path of β̂\hat\beta across ranks or numbers of blocks is a sensitivity analysis, not an exogeneity test. A stable path can be stably biased, and a moving path does not separate selection from treatment. Posterior predictive checks assess what the fitted model reproduces. Cross-validation assesses a stated prediction target. Neither establishes causal exchangeability. Regressing estimated factors on XX can describe overlap in the fitted representation, but it cannot test whether the unobserved true factor was independent of XX.

Panel data raise an additional distinction. Strict exogeneity conditions the current disturbance on past, present, and future covariates. Sequential exogeneity conditions it on information available through the current period. Lagging a predictor does not establish either condition when ties and attributes can affect one another over time.

The practical response is to measure important pre-treatment confounders, use actor fixed effects or correlated random-effects strategies when the needed within-actor variation exists, defend a design when one is available, and report sensitivity across plausible relational specifications. Actor fixed effects cannot separately estimate a time-invariant actor predictor. Dyadic cluster-robust standard errors can repair an uncertainty calculation under their assumptions, but they do not remove omitted-variable bias. Day 13 develops those distinctions.

In this course, read every AME, SRM, LDM, or blockmodel coefficient as a dependence-adjusted association unless a separate design supports stronger language. Its interval and structural summaries are model-based and inherit the model’s assumptions.

8.2 After the Model Runs: Interpret the Relational Surface, Not the Axes

A latent model gives us at least two kinds of results. It estimates measured associations conditional on the fitted latent structure, and it estimates a relational surface that summarizes patterned residual affinity. We should report both without turning an arbitrary axis into a named social dimension.

  1. Translate the measured association. Give its sign, uncertainty, and a probability or outcome contrast in meaningful units. Explain what was held fixed and call it dependence-adjusted unless a separate design supports causal language.
  2. Interpret invariant relationship quantities. Discuss fitted probabilities, pairwise distances, inner products, or entries of UV𝖳UV^{\mathsf T}. Do not name raw coordinate axes or separate AME factor columns as if their orientation and scaling were uniquely identified.
  3. Say what the latent term could contain. Strategy, institutions, geography, measurement, media coverage, and shared exposure may all contribute. The model does not label that mixture.
  4. Show whether the model earned the interpretation. Report sampling diagnostics, posterior predictive checks, the exact held-out prediction target, sensitivity to dimension or rank, and any important failures.
TipWhat We Can Say After the Applications

Syrian armed-organization cooperation: The fitted map places Al-Nusrah Front and Ahrar al-Sham Islamic Movement in the same tightly connected part of the cooperation network. They recorded an operation with one another, and each also recorded operations with 18 of the same other organizations. Simulated networks from the model reproduce much of the observed clustering. This largely reinforces the tactical-core finding from Day 9. It does not tell us why the organizations cooperated, and some clustering is created mechanically when one multi-organization operation is recorded as several pairwise ties.

Changing ICEWS state-pair relationships: The estimated chance for Iran toward Syria rises from about 2% in 2002 to 92% in 2012, while Syria toward Iran rises on a different timetable and reaches about 69% in 2011. The model shows that these changes are not just a by-product of Iran or Syria appearing in more relationships with every partner. It predicts unseen state-pair histories slightly better and reproduces much more of the observed year-to-year movement than the broad state-level model. It still misses reciprocity in several years and four of twelve transition intervals. These are descriptions of coded event volume, not estimates of conflict initiation, severity, or causal effects.

Nigerian armed conflict: Dorff, Gallop, and Minhas (2020) show how directed additive and multiplicative terms can identify broad initiators, broad targets, retaliation, and recurring opponent profiles in one model. Their predictive comparison is strong, but the post-2009 and civilian-targeting results remain conditional associations.

If a close or highly compatible pair makes sense, explain why using information outside the fitted coordinates. If it does not, treat that as a prompt to inspect measurement and fit rather than forcing a story onto the map.

9 Your Turn

TipExercises

Required.

1. Use ls2d$mkl$Z to find the two Syrian organizations that sit closest together and the two that sit farthest apart. Check gade_binary to see whether the close pair has an observed cooperation tie. Write one sentence that treats distance as a summary of the full tie pattern rather than proof of a shared ideology.

2. Compare Iran toward Syria with Syria toward Iran in pair_paths. Report the first year each direction crosses the outcome threshold, describe the fitted probability paths, and explain why the two directions should not be collapsed.

3. Use rank_cv_summary to compare ranks zero through four. State exactly what was held out, choose a rank, and explain why convergence alone does not make a higher rank useful.

Stretch.

4. Use gof_coverage and transition_plot_data to compare annual network reproduction with adjacent-year changes. Identify the largest remaining weakness and explain why a good annual fit does not guarantee a good transition fit.

5. Fit the Gade distance model again with a different seed. Plot the two maps before and after align_map_to_target(). Explain what translation, rotation, and reflection Procrustes removes and what a remaining difference could represent.

6. Fit a rank-two AME to a network of your own. Examine the invariant multiplicative surface, compare at least two starts, and state a held-out or posterior-predictive target before interpreting the factors.

Solutions: open after you have tried them 
## 1. Closest and farthest organizations
Zsol <- ls2d$mkl$Z
rownames(Zsol) <- rownames(gade_binary)
Dsol <- as.matrix(dist(Zsol))
Dsol[lower.tri(Dsol, diag = TRUE)] <- NA
close_pair <- which(Dsol == min(Dsol, na.rm = TRUE), arr.ind = TRUE)[1, ]
far_pair <- which(Dsol == max(Dsol, na.rm = TRUE), arr.ind = TRUE)[1, ]
close_names <- rownames(Zsol)[close_pair]
far_names <- rownames(Zsol)[far_pair]
close_names
far_names
gade_binary[close_names[1], close_names[2]]

## 2. What the changing partner surface adds
pair_paths[
  pair_paths$pair %in% c(
    "Iran toward Syria",
    "Syria toward Iran"
  ),
]

## 3. Rank and held-pair validation
rank_cv_summary

## 4. Fixed-fit simulation coverage by year
gof_coverage
transition_plot_data

## 5. Two seeds and Procrustes alignment
ls2d_b <- cache_fit("gade_ldm_2d_seed2", ergmm(
  CooperationNet ~ euclidean(d = 2),
  control = ergmm.control(
    sample.size = 4000,
    interval = 10,
    burnin = 10000
  ),
  seed = 12345,
  verbose = FALSE
))
Za <- ls2d$mkl$Z
Zb <- ls2d_b$mkl$Z
aligned_b <- align_map_to_target(Zb, Za)$aligned
par(mfrow = c(1, 2))
plot(Za, pch = 19, col = "#18453B", main = "Seed 6886")
plot(Zb, pch = 19, col = "#18453B", main = "Seed 12345")
plot(scale(Za, scale = FALSE), pch = 19, col = "#18453B", main = "After Alignment")
points(aligned_b, pch = 1, col = "#7BBD00")

## 6. Your own network
# fit_a <- lame(
#   Y = <your matrix, with diagonal set to NA>,
#   family = "binary",
#   symmetric = FALSE,
#   R = 2,
#   seed = 6886,
#   nscan = 6000,
#   burn = 6000,
#   odens = 10,
#   plot = FALSE,
#   verbose = FALSE,
#   print = FALSE
# )
# fit_b <- update(fit_a, seed = 12345)
# cor(c(fit_a$UVPM), c(fit_b$UVPM))

10 What Comes Next?

AME represents higher-order patterns through shared multiplicative profiles. It does not estimate a closure coefficient, and reproducing a transitive statistic does not identify closure as the process that generated the network. In the ICEWS application, rank two is useful because it improves a clearly defined held-pair target, reproduces more of the observed movement, and gives us a readable changing partner surface. Those successes still do not identify the process that created the coded event relationships.

The ERGM on Day 11 asks a different question by defining a joint probability model for the whole graph,

p(Y∣θ)∝exp{θ⊤g(Y)}. p(Y \mid \theta) \propto \exp\{\theta^\top g(Y)\}.

The chosen statistics g(Y)g(Y) can include edges, mutual ties, geometrically weighted shared partners, and other network features. An ERGM coefficient is interpreted through the change in those statistics when a tie is toggled, conditional on the rest of the graph. It is not automatically a causal effect of a triangle. AME and ERGM therefore represent dependence differently: AME uses latent actor and dyadic terms in the link-scale mean, while an ERGM uses chosen graph statistics in a joint graph distribution.

The course map now looks like this:

  • Describe the network: centrality and community.
  • Represent recurring actor and relationship patterns: SRM, blockmodels, latent distance, and AME.
  • Model graph structure directly: ERGM on Day 11.
  • Model network change: SAOM on Day 12.
  • Separate dependence, uncertainty, and causal identification: Day 13.
  • Compare the approaches on shared problems: Day 14.

11 Reading

  • Hoff, Raftery, and Handcock (2002), “Latent Space Approaches to Social Network Analysis,” JASA. The foundation for the distance model in Section 2.
  • Hoff (2007), “Modeling Homophily and Stochastic Equivalence in Symmetric Relational Data,” NeurIPS. The formal argument that an eigenmodel can weakly represent tie-propensity orderings from distance and class models, sometimes at a larger rank. This is not a claim of equal likelihoods at the same dimension.
  • Minhas, Hoff, and Ward (2019), “Inferential Approaches for Network Analysis: AMEN for Latent Factor Models,” Political Analysis. The AME framework developed in Section 7 through Section 7.7.
  • Gade, Gabbay, Hafez, and Kelly (2019), “Networks of Cooperation: Rebel Alliances in Fragmented Civil Wars,” Journal of Conflict Resolution. The source for the Syrian cooperation data. Their published model uses a square-root-transformed count outcome; our binary whole-period distance model is a teaching reanalysis.
  • Boschee, Lautenschlager, O’Brien, Shellman, Starz, and Ward (2015), “ICEWS Coded Event Data.” The event-data source used in the dynamic state-pair application.
  • Dorff, Gallop, and Minhas (2020), “Networks of Violence: Predicting Conflict in Nigeria,” The Journal of Politics. The published directed AME illustration in Section 7.
  • Raleigh, Linke, Hegre, and Karlsen (2010), “Introducing ACLED,” Journal of Peace Research. The event-data source used in the Nigeria application.
  • Lubold, Chandrasekhar, and McCormick (2023), “Identifying the Latent Space Geometry of Network Models Through Analysis of Curvature,” JRSS-B. The deeper warning that dimension and geometry are modeling choices rather than labels recovered automatically from a two-dimensional picture.
  • Krivitsky and Handcock (2008), “Fitting Position Latent Cluster Models for Social Networks With latentnet,” Journal of Statistical Software. The software reference for the Gade distance model and model-based clustering extension.
  • Weschle (2018), “Quantifying Political Relationships,” American Political Science Review, and Cheng and Minhas (2020), “Keeping Friends Close, but Enemies Closer: Foreign Aid Responses to Natural Disasters,” British Journal of Political Science. Two applied examples of using latent relational structure to study international relationships while connecting the model back to a substantive question.
  • Minhas, Dorff, Gallop, Foster, Liu, Tellez, and Ward (2022), “Taking Dyads Seriously,” Political Science Research and Methods. The simulation evidence for dependence adjustment and the correlated-omitted-variable boundary in Section 8.
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] patchwork_1.3.2       reshape2_1.4.5        ggplot2_4.0.3        
#>  [4] netify_1.5.3          lame_1.3.5            sna_2.8              
#>  [7] statnet.common_4.13.0 latentnet_2.12.0      ergm_4.12.0          
#> [10] network_1.20.0       
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6             tensorA_0.36.2.1         xfun_0.55               
#>  [4] htmlwidgets_1.6.4        ggrepel_0.9.6            rle_0.10.0              
#>  [7] lattice_0.22-5           vctrs_0.7.3              tools_4.3.3             
#> [10] Rdpack_2.6.4             generics_0.1.4           parallel_4.3.3          
#> [13] tibble_3.3.1             DEoptimR_1.1-4           pkgconfig_2.0.3         
#> [16] Matrix_1.6-5             data.table_1.18.0        checkmate_2.3.4         
#> [19] ggnewscale_0.5.2         RColorBrewer_1.1-3       S7_0.2.2                
#> [22] distributional_0.6.0     lifecycle_1.0.5          compiler_4.3.3          
#> [25] farver_2.1.2             stringr_1.6.0            precrec_0.14.5          
#> [28] htmltools_0.5.9          yaml_2.3.12              pillar_1.11.1           
#> [31] tidyr_1.3.2              MASS_7.3-60.0.1          cachem_1.1.0            
#> [34] trust_0.1-9              abind_1.4-8              robustbase_0.99-7       
#> [37] posterior_1.6.1          tidyselect_1.2.1         digest_0.6.39           
#> [40] mvtnorm_1.3-3            stringi_1.8.7            dplyr_1.2.1             
#> [43] purrr_1.2.2              labeling_0.4.3           fastmap_1.2.0           
#> [46] grid_4.3.3               cli_3.6.6                magrittr_2.0.5          
#> [49] loo_2.9.0                broom_1.0.11             withr_3.0.2             
#> [52] scales_1.4.0             backports_1.5.0          rmarkdown_2.30          
#> [55] matrixStats_1.5.0        igraph_2.2.2             otel_0.2.0              
#> [58] gridExtra_2.3            coda_0.19-4.1            memoise_2.0.1           
#> [61] evaluate_1.0.5           lpSolveAPI_5.5.2.0-17.15 knitr_1.51              
#> [64] rbibutils_2.4            rlang_1.2.0              Rcpp_1.1.1-1.1          
#> [67] glue_1.8.1               jsonlite_2.0.0           plyr_1.8.9              
#> [70] R6_2.6.1