---
title: "Day 9: The Social Relations Model and Blockmodels"
subtitle: "Advanced Network Analysis · ICPSR"
author: "Shahryar Minhas"
date: last-modified
engine: knitr
editor: source
format:
html:
theme: [cosmo, _theme.scss]
toc: true
toc-depth: 2
toc-location: left
toc-title: "Today"
number-sections: true
embed-resources: true
html-math-method: mathml
code-copy: true
code-tools: true
fig-width: 8
fig-height: 5
fig-align: center
df-print: kable
execute:
echo: true
warning: false
message: false
freeze: false
knitr:
opts_chunk:
comment: "#>"
---
```{r}
#| label: setup
#| include: false
source("_setup.R")
require_pkgs(c("lame", "netify", "blockmodels", "igraph", "NetMix", "ggplot2", "reshape2", "patchwork", "dplyr", "coda", "digest", "scales"))
if (utils::packageVersion("lame") < "1.3.5") {
stop("Day 9 requires lame 1.3.5 or later. Reinstall lame from GitHub, restart R, and render again.", call. = FALSE)
}
library(lame)
library(netify)
library(blockmodels)
library(igraph)
library(ggplot2)
library(reshape2)
library(patchwork)
theme_set(theme_bw(base_size = 13))
```
::: {.callout-note appearance="simple"}
## How to Use This Document
[Open the Day 9 teaching deck](day9_srm_blocks_teaching_deck.html){target="_blank"}.
Everything in this walkthrough runs. The expensive fits use `cache_fit()`, which saves the result in `cache/`. The first run estimates the model and later runs load the saved fit. Delete the corresponding cache file only when you deliberately want to refit a changed model.
```r
install.packages(c("netify", "blockmodels", "igraph", "NetMix", "ggplot2", "reshape2", "patchwork", "dplyr", "coda", "digest", "scales", "remotes"))
remotes::install_github("netify-dev/lame", dependencies = TRUE, build_vignettes = FALSE, upgrade = "never")
```
:::
::: {#lame-install .callout-note collapse="true" title="Detailed Install Instructions"}
`lame` is [in the CRAN submission process](https://github.com/netify-dev/lame/blob/main/cran-comments.md) 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.
### 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](https://github.com/netify-dev/lame):
```r
install.packages("remotes")
remotes::install_github(
"netify-dev/lame",
dependencies = TRUE,
build_vignettes = FALSE,
upgrade = "never"
)
```
The required compilation tools are [Rtools](https://cran.r-project.org/bin/windows/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`.
### Option 2: Use a Compiled GitHub Release
If you do not have compilation tools, open the [`lame` releases page](https://github.com/netify-dev/lame/releases), expand the assets for the latest release, and download the compiled file that matches your system. The release page currently provides Windows, Intel macOS, Apple silicon macOS, and x86_64 Linux builds.
Install the package dependencies from CRAN first:
```r
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:
```r
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:
```r
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:
```r
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:
```r
library(lame)
packageVersion("lame")
stopifnot(packageVersion("lame") >= "1.3.5")
```
:::
# What We Are Trying to Understand Today {#sec-today}
Day 8 was about constructing defensible network data. Day 9 asks what we can learn once the same actors appear in many relationships, which means those relationships are connected rather than independent rows. We will use the models to answer three concrete questions: which Syrian armed organizations worked with unusually many partners, whether cooperation was organized around a small tactical core, and whether states followed different defense-alliance paths after the Cold War.
1. **The social relations model (SRM)** asks whether some actors repeatedly appear in more relationships than their measured characteristics would lead us to expect.
2. **Blockmodels** ask whether actors follow a small number of recurring relationship patterns, first in one observed network and then across repeated networks.
The cross-sectional sections use claimed joint operations among armed organizations in the Syrian civil war collected by Gade et al. (2019). The same organizations and cooperation network carry us from the independent regression to the SRM and then to the blockmodel. The final sections use annual Integrated Crisis Early Warning System (ICEWS) material-conflict networks for the directed SRM and annual defense commitments for the longitudinal role model.
::: {.callout-important}
## The Thread That Connects the Models
We are not collecting methods for their own sake. We are asking what the measured variables still fail to explain. Do a few organizations take part in operations with many different partners? Is cooperation organized around a small tactical core? After the Cold War, did Poland, Romania, and Russia develop similar defense partners or follow different paths? The technical terms come after those questions are clear.
:::
## The Applied Model Compass {#sec-model-compass}
| Model | Research Question | What the Model Adds | What We Would Say About the Result |
|---|---|---|---|
| Independent dyadic regression | Which observed characteristics are associated with cooperation? | A familiar first comparison | "A useful first look, but it acts as if ANF-ASIM tells us nothing about ANF-ISIL." |
| SRM | Which organizations took part in more joint operations across many partners than the measured characteristics predict? | One adjusted cooperation score per organization | "ANF and ASIM took part in more joint operations across their partnerships than their measured characteristics would lead us to expect." |
| Cross-sectional SBM | Is cooperation spread broadly, divided into separate camps, or organized around a small set of widely connected organizations? | A role for each organization and a cooperation rate for each role pairing | "ANF and ASIM form a small tactical core with links across a much larger set of organizations that seldom worked with one another." |
| Longitudinal SRM | Which states appeared in many above-threshold ICEWS relationships, either as the coded source or the coded target, and when did that change? | Separate yearly scores for events coded from and toward each state, with bootstrap uncertainty | "Syria appeared in relatively few above-threshold relationships through 2010 and in many more after 2011, especially as the state toward which events were coded." |
| Repeated-network NetMix | Did Poland, Romania, and Russia develop similar sets of defense partners after the Cold War, or did their treaty networks take different paths? | A yearly summary of which recurring partner pattern each state's treaty list most resembles | "Poland's treaty partners begin to resemble the dense multilateral system in 1997, Romania remains mostly outside it in this panel, and Russia follows a different regional pattern." |
# Cooperation Among Armed Organizations in Syria {#sec-gade}
Gade et al. (2019) study armed opposition organizations in the Syrian civil war, often described as rebel groups or militias, and ask: **which organizations carried out joint tactical operations with one another?** They assembled claims of joint operations from July 2012 through June 2015 and measured ideological distance, power differences, shared operational location, and common state sponsorship. The prepared regression table contains 30 organizations with complete covariates; the binary network used for the block analysis contains the 31 organizations that recorded at least one cooperative relationship.
The dependent variable is the number of claimed joint operations for each pair. It is symmetric because a joint operation belongs to the pair rather than traveling from a sender to a receiver. We follow the published analysis and use the square root of the count. The transformation makes the very skewed counts more compatible with a Gaussian working model and keeps the classroom analysis close to the published application. Ideological distance runs from 0 to 4, power is recorded in thousands of estimated fighters, and shared location and sponsorship are binary. Current versions of `lame` also support an overdispersed Poisson family, which would be an appropriate sensitivity model for a new paper.
```{r}
#| label: load-gade
load("data/gade_app.rda")
gade <- data
rm(data)
gade$sqrt_coop <- sqrt(gade$coopActions)
checkpoint(actors = length(unique(gade$actor1)), directed_rows = nrow(gade), observed_pairs = nrow(gade) / 2, cooperative_pairs = sum(gade$coopActions > 0) / 2)
head(gade[, c("actor1", "actor2", "coopActions", "ideol_diff", "powerdiff", "loc", "spons")])
```
The duplicated directions are a storage convention. Because the relationship is symmetric, `A-B` and `B-A` carry the same outcome and covariates. The model counts each unordered pair once.
## Start With the Expectations
- If ideological compatibility helps organizations coordinate tactics and agree about long-term goals, greater `ideol_diff` should be associated with less cooperation.
- If organizations seek similarly powerful partners, greater `powerdiff` should be associated with less cooperation.
- Organizations operating in the same location have more opportunities to conduct joint operations, so `loc` should be positively associated with cooperation.
- If shared external patrons coordinate their clients, `spons` should be positively associated with cooperation. The published study does not find clear evidence for this last expectation.
These are associational expectations. The data are observational, and ideology, power, location, sponsorship, and cooperation may affect one another. A network model changes the dependence assumptions. It does not create random assignment.
# Put the Data Into the Form the Model Needs {#sec-netify}
```{r}
#| label: netify-gade
gade_net <- netify(
input = gade,
actor1 = "actor1", actor2 = "actor2",
symmetric = TRUE, weight = "sqrt_coop",
nodal_vars = c("averageId_actor1", "size_actor1", "spons_actor1"),
dyad_vars = c("ideol_diff", "powerdiff", "loc", "spons"),
dyad_vars_symmetric = rep(TRUE, 4),
missing_to_zero = TRUE
)
gade_lame <- to_lame(gade_net, family = "normal")
Yg <- gade_lame$Y
Xg <- gade_lame$Xdyad
Ng <- gade_lame$Xrow
dim(Yg)
dim(Xg)
colnames(Ng)
```
`Yg` is a symmetric 30 by 30 matrix. `Xg` is a 30 by 30 by 4 dyadic-covariate array. `Ng` contains one row per organization. The diagonal of `Yg` is missing because an organization cannot conduct a joint operation with itself.
::: {.callout-note}
## What Changes in a Directed Network?
This application is symmetric, so each organization receives one additive actor effect. In a directed network, such as one state directing a threat toward another, the SRM separates the source side from the target side:
$$
\begin{aligned}
y_{ij}&=\mathbf{x}_{ij}^{\mathsf T}\boldsymbol\beta+a_i+b_j+\epsilon_{ij},\\
\begin{pmatrix}a_i\\b_i\end{pmatrix}
&\sim\mathcal N\!\left[
\begin{pmatrix}0\\0\end{pmatrix},
\begin{pmatrix}\sigma_a^2&\sigma_{ab}\\\sigma_{ab}&\sigma_b^2\end{pmatrix}
\right],\qquad
\operatorname{Corr}(\epsilon_{ij},\epsilon_{ji})=\rho .
\end{aligned}
$$
Here $a_i$ records how often actor $i$ appears on the source side beyond the measured predictors, while $b_i$ records how often the same actor appears on the target side. The variances $\sigma_a^2$ and $\sigma_b^2$ describe how much those adjusted source-side and target-side scores differ across actors. The covariance $\sigma_{ab}$ asks whether actors that are often sources are also often targets. The residual correlation $\rho$ asks a different question: after the actor effects and predictors enter, is an unexpectedly high $i\!\rightarrow\!j$ outcome paired with an unexpectedly high $j\!\rightarrow\!i$ outcome? The label $\mathrm{iid}$ means independent and identically distributed across actors. The symmetric model is not a lesser model. It matches the meaning of a joint operation.
:::
# Why Independent Dyads Are the Wrong Benchmark {#sec-glm}
An independent dyadic regression says
$$
y_{ij} = \beta_0 + \mathbf{x}_{ij}^{\mathsf T}\boldsymbol\beta + \mathbf{w}_i^{\mathsf T}\boldsymbol\gamma + \mathbf{w}_j^{\mathsf T}\boldsymbol\gamma + \epsilon_{ij},
\qquad \epsilon_{ij} \stackrel{\text{iid}}{\sim} N(0,\sigma^2).
$$
The intercept $\beta_0$ is the baseline. The vector $\mathbf x_{ij}$ contains measured pair attributes and $\boldsymbol\beta$ their coefficients. The vectors $\mathbf w_i$ and $\mathbf w_j$ contain measured organization attributes and $\boldsymbol\gamma$ their coefficients. The error $\epsilon_{ij}$ is what remains unique to the pair. The label $\mathrm{iid}$ means independent and identically distributed: every pair error is assumed to be a separate draw from the same normal distribution.
The conditional mean is sensible, but the IID statement says that evidence from `ANF-ASIM` is unrelated to evidence from `ANF-ISIL`, even though both relationships involve ANF. If ANF is unusually cooperative for reasons we did not measure, all pairs involving ANF will lean in the same direction.
## Return the Residuals to the Network
```{r}
#| label: independent-lm
gade_pair_id <- apply(gade[, c("actor1", "actor2")], 1, function(x) paste(sort(x), collapse = "|"))
gade_unique <- gade[!duplicated(gade_pair_id), ]
gade_unique$ideology_sum <- gade_unique$averageId_actor1 + gade_unique$averageId_actor2
gade_unique$size_sum <- gade_unique$size_actor1 + gade_unique$size_actor2
gade_unique$sponsor_sum <- gade_unique$spons_actor1 + gade_unique$spons_actor2
glm_gade <- lm(sqrt_coop ~ ideology_sum + size_sum + sponsor_sum + ideol_diff + powerdiff + loc + spons, data = gade_unique)
summary(glm_gade)$coefficients[c("ideol_diff", "powerdiff", "loc", "spons"), ]
```
<div id="residual-actor"></div>
```{r}
#| label: residual-actor
Eg <- matrix(NA_real_, nrow(Yg), ncol(Yg), dimnames = dimnames(Yg))
for (ii in seq_len(nrow(gade_unique))) {
a <- gade_unique$actor1[ii]
b <- gade_unique$actor2[ii]
Eg[a, b] <- Eg[b, a] <- resid(glm_gade)[ii]
}
actor_resid <- sort(rowMeans(Eg, na.rm = TRUE), decreasing = TRUE)
head(actor_resid, 8)
sd(actor_resid)
```
If the independent model had absorbed the recurring-actor structure, these means would hover near zero. Instead, the same organizations repeatedly appear above or below the regression line.
# The Symmetric Social Relations Model {#sec-srm}
The SRM adds one latent effect for every organization:
$$
y_{ij} = \beta_0 + \mathbf{x}_{ij}^{\mathsf T}\boldsymbol\beta + \mathbf{w}_i^{\mathsf T}\boldsymbol\gamma + \mathbf{w}_j^{\mathsf T}\boldsymbol\gamma + a_i + a_j + \epsilon_{ij},
\qquad a_i \stackrel{\mathrm{iid}}{\sim} N(0,\sigma_a^2).
$$
Each organization gets one effect, $a_i$, and that same effect appears in every partnership involving the organization. A positive value means the organization took part in more joint operations across its partners than the model expected. The variance $\sigma_a^2$ summarizes how much organizations differ from one another in this overall adjusted activity. The distribution also partially pools the organization estimates toward zero, so they are not simply raw cooperation counts.
For a directed network, replace $a_i+a_j$ with $a_i+b_j$. The model then estimates sender variance, receiver variance, their covariance, and dyadic reciprocity. Back and Kenny (2010) provide the general SRM decomposition; Hoff (2005) develops the network regression version; Minhas et al. (2022) show why these components matter for network inference.
# The Inference Issue: Three Cases, Not One Slogan {#sec-inference}
A missing factor can move a coefficient up or down. In an SRM, this is not only a generic omitted-variable problem. The usual random-effects specification also assumes that the latent actor effects are unrelated to the included predictors after conditioning on the model. If organizations with large positive actor effects systematically have different observed characteristics, the coefficient can blend a predictor association with persistent unmeasured organization differences.
Let $x_{ij}$ be a measured predictor such as ideological distance, and let $z_{ij}$ be a missing process such as whether two organizations faced the same immediate battlefield threat. Suppose the true model is
$$
y_{ij} = \beta x_{ij} + \gamma z_{ij} + \epsilon_{ij},
$$
but we regress $y$ only on $x$. In the linear case, the omitted-variable calculation is
$$
\operatorname{plim}(\widehat\beta_{\text{naive}})
= \beta + \gamma\frac{\operatorname{Cov}(x,z)}{\operatorname{Var}(x)}.
$$
The probability limit, $\operatorname{plim}$, is the value an estimate approaches as the sample becomes very large. Here $\beta$ is the association we want for $x$, while $\gamma$ describes how the missing process $z$ relates to cooperation. The equation gives us three distinct cases.
| Relationship between measured $x$ and omitted $z$ | What happens to the point estimate? | What can still go wrong? |
|---|---|---|
| Positive covariance | If $\gamma>0$, the estimate for $x$ is pulled upward | Some of the missing process is mistakenly credited to $x$ |
| Negative covariance | If $\gamma>0$, the estimate for $x$ is pulled downward and can cross zero | A weak or negative estimate can hide a more positive direct association |
| Zero covariance | The estimate for $x$ can remain centered on $\beta$ | If $z$ contains recurring actor effects, the usual uncertainty calculation can still be wrong because pairs sharing actors remain connected |
The sign reverses if $\gamma<0$. With several measured predictors, the relevant overlap is what remains between $x$ and $z$ after accounting for the other included variables. "Networks make coefficients biased" is therefore too crude. A coefficient can be centered correctly even when the usual uncertainty calculation is wrong. Coefficient bias appears when the missing process both matters for the outcome and overlaps with the predictor.
## A Concrete Reading
Imagine that ideological distance is our measured $x$ and that facing the same immediate battlefield threat is the missing $z$. A common threat can make cooperation more likely. If ideologically similar organizations also tend to face the same opponents, the ideology coefficient can receive credit for some of that unmeasured threat. If common threats are what bring otherwise dissimilar organizations together, the bias can run in the opposite direction. The same logic applies to the SRM actor effects: if an observed organization characteristic is systematically related to the persistent activity stored in $a_i$, the random-effects independence assumption is not credible and the measured coefficient can absorb part of the actor pattern.
The algebra tells us the direction conditional on assumptions. The data do not reveal the sign of an unmeasured covariance by themselves.
::: {.callout-important}
## What an SRM Does Not Solve
The SRM accounts for patterns shared by pairs involving the same actor. It does not turn observational data into a causal design. If an unmeasured actor or pair process overlaps with a measured predictor, the coefficient can still blend the two explanations. In random-effects language, the key assumption is $E(a_i\mid X)=0$ in the symmetric model, or $E(a_i,b_i\mid X)=0$ in a directed model. The actor effects may be latent, but they are not automatically independent of the covariates.
:::
::: {.callout-note}
## How a Mundlak Specification Works
With repeated observations and a time-varying actor predictor $q_{it}$, calculate each actor's across-time mean $\bar q_i$ and its yearly deviation $q_{it}-\bar q_i$. Include both in the SRM. The deviation coefficient uses within-actor change, while the mean coefficient records the stable between-actor difference that may be related to the random effect. For a directed SRM, make this decomposition separately for source-side and target-side actor covariates and pass the resulting period-specific matrices through `Xrow` and `Xcol`. For a time-varying dyadic predictor, put its pair mean and within-pair deviation into `Xdyad`.
```{r}
#| eval: false
# q_it is one row per actor-year
q_it$q_mean <- ave(q_it$q, q_it$actor, FUN = mean)
q_it$q_within <- q_it$q - q_it$q_mean
# Build one actor-by-covariate matrix per year, then fit:
fit_mundlak <- lame(
Y = Y_list,
Xdyad = X_dyad_list,
Xrow = Xrow_list, # q_within and q_mean for source-side actors
Xcol = Xcol_list, # q_within and q_mean for target-side actors
family = "binary",
symmetric = FALSE,
R = 0
)
```
This is possible inside an SRM because it changes the observed part of the model while retaining the source and target random effects. It does not solve every kind of confounding, and one cross-section does not contain the within-actor time variation needed for this decomposition.
:::
Minhas et al. (2022) make the boundary especially clear. When an omitted network pattern is independent of $x$ and has a form the latent model can represent, the model can recover the measured coefficient well in their simulation. When the omitted quantity contains part of $x$, write $z=\alpha x+e$ and substitute it into the outcome equation:
$$
y=\beta x+\gamma z+\epsilon
=(\beta+\gamma\alpha)x+(\gamma e+\epsilon).
$$
The fitted data only reveal the combined association $\beta+\gamma\alpha$. In plain language, if $x$ and the missing network pattern move together, the model cannot know how much credit belongs to each one. If the shared part moves in the opposite direction, write $z=-\alpha x+e$ and the combined association becomes $\beta-\gamma\alpha$. The sign in front of the shared $x$ component controls the direction; changing the sign of the unrelated remainder $e$ does not. No SRM or other latent-variable model can separate $\beta$ without additional assumptions or information.
# How the SRM Is Estimated With MCMC {#sec-srm-mcmc}
At a high level, the model is trying to find combinations of coefficients and organization effects that make the observed pattern of joint operations plausible. A proposed positive effect for ANF, for example, is useful only if it helps account for ANF's outcomes across all of its partners without forcing the measured covariate coefficients or the remaining variance into implausible values.
Technically, `ame()` performs Bayesian estimation. For the observed network $Y$, coefficients $\boldsymbol\beta$, actor effects $\mathbf a$, actor-effect variance $\sigma_a^2$, and remaining variance $\sigma_\epsilon^2$, the target is the joint posterior
$$
p(\boldsymbol\beta,\mathbf a,\sigma_a^2,\sigma_\epsilon^2\mid Y,X)
\propto
p(Y\mid X,\boldsymbol\beta,\mathbf a,\sigma_\epsilon^2)\,
p(\boldsymbol\beta)\,p(\mathbf a\mid\sigma_a^2)\,
p(\sigma_a^2)\,p(\sigma_\epsilon^2).
$$
The first term is the likelihood: parameter settings receive more posterior weight when they make the observed cooperation counts more probable. The remaining terms are priors, which regularize coefficients and partially pool organization effects toward zero. The sampler does not maximize this expression and stop at one best value. It constructs draws whose long-run distribution is the posterior.
One MCMC sweep updates connected blocks of unknown quantities:
1. Start with provisional coefficients, organization effects, and variance parameters.
2. Draw the measured-variable coefficients conditional on the current actor effects and variances.
3. Draw every organization's actor effect conditional on its observed partnerships and the current coefficients.
4. Draw the actor-effect and remaining-error variances conditional on the current effects and residuals.
5. Repeat, discard the early burn-in sweeps, and retain later draws for posterior summaries.
A chain is a connected walk through plausible explanations of one observed network. It is not a collection of newly observed Syrian conflicts. The posterior mean is an average across retained explanations, a credible interval contains the central posterior mass under the model, and a trace plot checks whether the sampler moved through a stable region. These diagnostics assess computation. They do not establish that the model, priors, or exogeneity assumptions are substantively correct.
# Fit the Independence Model and SRM {#sec-srm-fit}
```{r}
#| label: fit-independence
gade_indep <- cache_fit("gade_indep_v2", ame(
Y = Yg, Xdyad = Xg, Xrow = Ng,
family = "normal", symmetric = TRUE,
nvar = FALSE, R = 0, seed = 6886,
nscan = 10000, burn = 10000, odens = 10,
verbose = FALSE
))
```
```{r}
#| label: fit-srm
gade_srm <- cache_fit("gade_srm_v2", ame(
Y = Yg, Xdyad = Xg, Xrow = Ng,
family = "normal", symmetric = TRUE,
nvar = TRUE, R = 0, seed = 6886,
nscan = 10000, burn = 10000, odens = 10,
verbose = FALSE
))
```
In these calls, `nscan` sets the number of MCMC updates, `burn` sets how many early updates are discarded while the chain settles, and `odens` saves every tenth later update. Saving every tenth draw reduces storage; it does not manufacture new information.
The only modeling switch is `nvar`: `FALSE` treats the pairs as independent conditional on the covariates, while `TRUE` estimates one persistent organization effect. `R = 0` says there are no multiplicative latent factors yet.
```{r}
#| label: srm-diagnostics
#| fig-height: 10
# Keep the classroom diagnostic focused enough that every trace and
# posterior density remains readable at the width of the walkthrough.
gade_srm_diagnostics <- gade_srm
gade_srm_diagnostics$BETA <- gade_srm$BETA[, c("ideol_diff_dyad", "powerdiff_dyad", "loc_dyad", "spons_dyad"), drop = FALSE]
colnames(gade_srm_diagnostics$BETA) <- c("Ideological distance", "Power difference", "Shared location", "Shared sponsorship")
gade_srm_diagnostics$VC <- gade_srm$VC[, c("va", "ve"), drop = FALSE]
colnames(gade_srm_diagnostics$VC) <- c("Actor-effect variance", "Remaining variance")
trace_plot(gade_srm_diagnostics, ncol = 2, title = "MCMC Diagnostics for the Main Results")
round(coda::effectiveSize(coda::mcmc(gade_srm$BETA)))
```
The diagnostic focuses on the four pair-level predictors we interpret below, the variance of the organization effects, and the remaining outcome variance. Each trace should move around a stable horizontal region without a sustained trend, while each density summarizes where that quantity spent most of its saved draws.
The `coda::effectiveSize()` line belongs with the settings in the model call. After discarding `burn = 10000` settling-in updates, `nscan = 10000` supplies the post-burn sampling updates and `odens = 10` saves every tenth one, leaving 1,000 stored draws. Consecutive MCMC draws can repeat much of the same information. ESS asks how many independent draws would carry about as much information as those 1,000 correlated saved draws, separately for each coefficient. An ESS near 1,000 means little information was lost to autocorrelation for that parameter; an ESS of 250 means the 1,000 saved draws carry roughly the information of 250 independent draws. Increasing `odens` only saves fewer draws. It does not improve the underlying chain.
## Read the Coefficients in Terms of Actions
<div id="srm-compare"></div>
```{r}
#| label: srm-compare
compare_fits(gade_indep, gade_srm, names = c("independent", "srm"))
```
```{r}
#| label: srm-numbers
b_ind <- colMeans(gade_indep$BETA)
b_srm <- colMeans(gade_srm$BETA)
ci_srm <- t(apply(gade_srm$BETA, 2, quantile, c(.025, .975)))
round(cbind(independent = b_ind, srm = b_srm, lo = ci_srm[, 1], hi = ci_srm[, 2]), 3)
```
The point estimates suggest less cooperation as ideological and power differences grow, more cooperation among organizations operating in the same location, and no clear shared-sponsorship pattern. The uncertainty is not equally decisive. In the SRM, the power-difference interval stays below zero, while the ideology, location, and sponsorship intervals cross zero. The estimates are on the square-root count scale, so they should not be narrated as percentage changes. The comparison shows which patterns survive persistent organizational activity rather than licensing a stronger causal story.
Do not reduce the comparison to "the standard errors got wider." Both point estimates and interval widths can move, and they need not move in the same direction. The new model changes the conditional mean by adding the actor effects and changes the covariance structure by recognizing repeated organizations.
## The Actor Effects Tell Us About Organizations {#sec-actor-effects}
```{r}
#| label: actor-effects
# APM stores the posterior mean additive actor effect for each organization.
gade_effects <- data.frame(actor = names(gade_srm$APM), effect = as.numeric(gade_srm$APM))
gade_effects <- gade_effects[order(gade_effects$effect), ]
tail(gade_effects, 8)
ggplot(gade_effects, aes(effect, reorder(actor, effect))) +
geom_vline(xintercept = 0, color = "grey70") +
geom_point(color = "#18453B", size = 2.5) +
labs(x = "Adjusted cooperative-activity effect", y = NULL)
```
Al-Nusrah Front (ANF) and Ahrar al-Sham Islamic Movement (ASIM) receive the two highest adjusted cooperation scores in this fit. ANF was an al-Qaeda-affiliated sectarian-jihadist organization during this period, while ASIM was a major Salafist-nationalist organization. Both recorded more joint operations across their partnerships than their measured ideology, size, sponsorship, and location would lead the model to expect. This does not mean either organization was inherently cooperative. The high values could reflect capacity, battlefield reach, brokerage, geography, reporting, or some combination of these. The plot shows posterior mean scores for this network and period, not uncertainty intervals for each organization.
Minhas et al. (2022) use the same logic in international conflict. Their reanalysis of Reiter and Stam shows that sender effects identify countries that initiate more or fewer disputes than the observed covariates predict. The latent effects become clues about missing explanations, not merely nuisance controls.
# Use the SRM as a Latent-Variable Tool {#sec-latent-tool}
The additive effect $a_i$ is already latent, meaning we infer it from the pattern of relationships rather than observe it directly. It can be used in two ways.
1. **A fuller model of recurring actors:** include $a_i$ so the measured variables do not have to explain every recurring difference among organizations. This helps with the dependence problem but does not guarantee that actor effects are unrelated to the predictors.
2. **A clue worth following:** inspect which organizations repeatedly sit above or below the model's prediction, compare that pattern with outside information, and ask what mechanism the measured variables missed.
That second use requires humility. The model shows us a repeated pattern that our measured variables missed. It cannot tell us by itself whether the explanation is capacity, brokerage, geography, reporting, or something else. Giving the effect a name requires evidence beyond the fitted network.
## What the SRM Still Misses
The SRM has now done its job: it accounts for organizations that work with many partners and shows us which measured associations remain after that adjustment. It still cannot tell us whether the same sets of organizations repeatedly choose the same partners. For example, one score cannot distinguish a collection of organizations that mostly cooperate with one another from a small set that connects otherwise separate organizations. We use blockmodels to ask that next question. Day 10 will develop the more flexible AME version carefully.
# Posterior Predictive Check: What the SRM Still Misses {#sec-gof}
```{r}
#| label: srm-gof
#| fig-height: 3.8
gof_plot(gade_srm, statistics = c("sd.rowmean", "trans.dep"))
```
The first check, `sd.rowmean`, is the standard deviation across organizations of their mean square-root cooperation count across partnerships. Put simply, it asks whether most organizations took part in similar numbers of joint operations or whether a few participated far more broadly than the rest. The second, `trans.dep`, is a standardized product of centered outcomes across connected triples. It asks whether high cooperation counts tend to appear among connected sets of three organizations, but it is not the percentage of triples that close.
The additive SRM reproduces the uneven actor activity but misses the network's triadic pattern. That is exactly what its formula implies: one actor effect changes every relationship involving that organization, but it cannot express a special affinity among a set of organizations.
::: {.callout-important}
## The SRM Inference Checklist
1. State whether the network is directed or symmetric.
2. Name the observational unit and the actors that repeat across dyads.
3. Separate point-estimate movement from interval-width movement.
4. Inspect the latent actor effects as quantities worth interpreting, but do not reify them.
5. Ask whether omitted structure might be correlated with the predictor of interest.
6. Check convergence and simulate from the fitted model before interpreting it.
:::
# From Continuous Latent Effects to Discrete Roles {#sec-blocks-bridge}
The SRM gives every organization one adjusted score for how broadly it cooperates. The fit check shows that broad participation is not the whole story. We next ask whether cooperation was spread evenly, divided into separate camps, or organized around a small set of organizations that worked across a much less connected field. A stochastic blockmodel summarizes that structure with a few roles, but those fitted roles are descriptions of partner lists rather than automatically meaningful real-world communities.
For the block analysis, we use the observed 31-organization cooperation network from the same study. One additional organization appears in the network matrix but not in the prepared regression table. The regression and block analyses therefore share a source and setting, but they are not nested fits to an identical analytic network. The regression uses square-root counts among 30 organizations with complete covariates. The blockmodel uses a binary indicator of any recorded cooperation among the 31 organizations with at least one observed cooperative relationship. Its periphery is the periphery among observed collaborators, not every armed organization active in Syria.
```{r}
#| label: load-block-data
load("data/gade_blocks.rda")
gade_block_net <- netify(
gadeData,
symmetric = TRUE,
diag_to_NA = TRUE,
missing_to_zero = FALSE
)
gade_binary_net <- binarize(gade_block_net, threshold = 0)
stopifnot(all(unlist(validate_netify(gade_binary_net))))
gade_A <- get_adjacency(gade_binary_net)
diag(gade_A) <- 0L
checkpoint(actors = nrow(gade_A), cooperative_pairs = sum(gade_A) / 2, density = round(mean(gade_A[upper.tri(gade_A)]), 3))
```
`netify()` preserves the organization labels and symmetric-network declaration, while `binarize()` makes the outcome transformation explicit. `blockmodels` requires a plain adjacency matrix, so `get_adjacency()` performs the final package-specific handoff. The binary outcome asks whether an unordered organization pair recorded at least one claimed joint operation at any point from July 2012 through June 2015. It discards how many operations the pair recorded. Every probability below refers to that whole-period relationship, not to a single operation, a future alliance, or the probability that two organizations share long-term goals.
# The Cross-Sectional Blockmodel {#sec-blockmodel}
The model is easiest to understand as a compressed description of partner lists. Suppose one role contains organizations that cooperate with many different organizations, while another contains organizations that mostly cooperate with the first role and rarely with one another. We would not need a separate probability for every organization pair. We would need each organization's role and three probabilities: core with core, core with periphery, and periphery with periphery.
That description is useful because it can represent more than community structure. Communities imply high within-role and low across-role cooperation. A core-periphery pattern instead implies high core-to-core and core-to-periphery cooperation but low periphery-to-periphery cooperation. The fitted role names come after we inspect those probabilities and the organizations associated with them.
Let $K$ be the number of roles, let $z_i\in\{1,\ldots,K\}$ be organization $i$'s role, and let $\theta_{gh}$ be the probability of a recorded tie between an organization in role $g$ and one in role $h$. Then
$$
Y_{ij} \mid z_i,z_j \sim \operatorname{Bernoulli}(\theta_{z_i z_j}).
$$
The central quantity is the matrix $\boldsymbol\theta$, which collects all role-to-role probabilities. Its diagonal describes cooperation within roles; its off-diagonal describes cooperation across roles. Conditional on the roles and this matrix, the basic SBM treats pair outcomes as independent Bernoulli observations.
### What Makes One Block Assignment Better Than Another
If the roles $\mathbf z$ were known, the likelihood would be
$$
L(\mathbf z,\boldsymbol\theta;Y)
=\prod_{i<j}
\theta_{z_i z_j}^{Y_{ij}}
(1-\theta_{z_i z_j})^{1-Y_{ij}}.
$$
Every recorded cooperation tie rewards assignments whose relevant $\theta_{z_i z_j}$ is large. Every recorded non-tie rewards assignments whose relevant $\theta_{z_i z_j}$ is small. Non-ties matter because a role is defined by the whole partner pattern, not only by an organization's observed partners. The likelihood therefore prefers a compression in which pairs assigned to the same role pairing behave consistently.
If we temporarily assigned ANF to a proposed core role, the model would evaluate all of ANF's ties and non-ties using the core-to-core and core-to-periphery rates. If ANF's partner list is much more probable under that assignment than under a peripheral assignment, the evidence moves ANF's role probability toward the core. The same calculation is made for every organization while the role-to-role rates are also unknown.
## How Variational Expectation-Maximization Learns the Roles {#sec-block-estimation}
The exact likelihood requires summing over $K^n$ possible role assignments, which is already impossible for modest networks. `blockmodels` therefore uses variational expectation-maximization, or VEM, to maximize an approximation to the observed-data likelihood.
1. Give each organization a tentative set of role probabilities. At first, ANF might be treated as 60% role 1 and 40% role 2.
2. In the variational E-step, update each organization's role probabilities using all of its ties and non-ties and the current role-to-role rates.
3. In the M-step, update each $\theta_{gh}$ as the expected number of ties between roles $g$ and $h$ divided by the expected number of possible pairs, where uncertain memberships provide the weights.
4. Recalculate memberships using the new rates, then rates using the new memberships.
5. Stop when another round barely improves the variational lower bound. Repeat from different starting guesses because the objective is not globally concave and the algorithm can settle on local solutions.
The optimized quantity is an evidence lower bound,
$$
\operatorname{ELBO}(q,\boldsymbol\theta)
=E_q[\log p(Y,\mathbf z\mid\boldsymbol\theta)]
-E_q[\log q(\mathbf z)],
$$
where $q(\mathbf z)$ is a tractable approximation to the unknown distribution over role assignments. The first term rewards assignments and tie probabilities that explain the observed network. The second term accounts for the uncertainty in the approximate memberships. Raising the ELBO improves the approximation, and at its optimum it supplies an approximate maximum-likelihood fit. It is not a posterior simulation, a goodness-of-fit test, or evidence that the selected roles are natural kinds.
For fixed soft memberships $\tau_{ig}=q(z_i=g)$, define the membership weight for an unordered pair as $w_{ij,gg}=\tau_{ig}\tau_{jg}$ for a within-role cell and $w_{ij,gh}=\tau_{ig}\tau_{jh}+\tau_{ih}\tau_{jg}$ for $g\neq h$. The Bernoulli rate update then has an especially useful interpretation:
$$
\widehat\theta_{gh}
=\frac{\sum_{i<j}w_{ij,gh}Y_{ij}}
{\sum_{i<j}w_{ij,gh}}.
$$
The numerator is the membership-weighted number of observed ties for that role pairing, and the denominator is the membership-weighted number of possible pairs. This is what the phrase “update the cooperation rates” means. Variational inference is fast, but it usually understates membership uncertainty, so we compare random starts, inspect nearby $K$, and simulate networks from the fitted roles.
<div id="fit-gade-blocks"></div>
```{r}
#| label: fit-gade-blocks
set.seed(6886)
gade_sbm <- BM_bernoulli("SBM_sym", gade_A, verbosity = 0, plotting = "")
gade_sbm$estimate()
plot(seq_along(gade_sbm$ICL), gade_sbm$ICL, type = "b", pch = 19, xaxt = "n", xlab = "Number of blocks (K)", ylab = "ICL", main = "Syrian Cooperation: Choosing K")
axis(1, at = seq_along(gade_sbm$ICL))
K_gade <- which.max(gade_sbm$ICL)
K_gade
round(gade_sbm$ICL, 1)
```
The integrated completed likelihood (ICL) rewards fit while penalizing unnecessary roles, and larger values are preferred here. ICL prefers two roles, but the three-role solution is close enough to inspect. The right conclusion is not that nature contains exactly two types. The two-role model gives the clearest balance between simplicity and fit, while the periphery can be divided more finely. As a stability check, ten additional random seeds all select two roles and recover the same ANF-ASIM core.
```{r}
#| label: block-seed-check
#| code-fold: true
#| code-summary: "Check ten additional random starts"
check_block_seed <- function(seed) {
set.seed(seed)
fit <- BM_bernoulli("SBM_sym", gade_A, verbosity = 0, plotting = "")
invisible(capture.output(fit$estimate()))
selected_k <- which.max(fit$ICL)
z <- apply(fit$memberships[[selected_k]]$Z, 1, which.max)
names(z) <- rownames(gade_A)
theta <- fit$model_parameters[[selected_k]]$pi
core <- which.max(rowMeans(theta))
data.frame(seed = seed, selected_k = selected_k, anf_asim_core = ifelse(identical(sort(names(z)[z == core]), c("ANF", "ASIM")), "yes", "no"))
}
block_seed_check <- do.call(rbind, lapply(1:10, check_block_seed))
block_seed_check
```
<div id="gade-block-membership"></div>
```{r}
#| label: gade-block-membership
z_gade <- apply(gade_sbm$memberships[[K_gade]]$Z, 1, which.max)
names(z_gade) <- rownames(gade_A)
split(names(z_gade), z_gade)
theta_gade <- gade_sbm$model_parameters[[K_gade]]$pi
core_block <- which.max(rowMeans(theta_gade))
periphery_block <- setdiff(seq_len(K_gade), core_block)
block_order <- c(core_block, periphery_block)
theta_ordered <- theta_gade[block_order, block_order, drop = FALSE]
dimnames(theta_ordered) <- list(c("core", "periphery"), c("core", "periphery"))
core_actors <- names(z_gade)[z_gade == core_block]
periphery_actors <- names(z_gade)[z_gade == periphery_block]
within_core <- gade_A[core_actors, core_actors, drop = FALSE]
within_core <- within_core[upper.tri(within_core)]
core_periphery <- gade_A[core_actors, periphery_actors, drop = FALSE]
within_periphery <- gade_A[periphery_actors, periphery_actors, drop = FALSE]
within_periphery <- within_periphery[upper.tri(within_periphery)]
cell_summary <- function(x, fitted) c(observed_ties = sum(x), possible_pairs = length(x), observed_rate = mean(x), fitted_probability = fitted)
block_cell_summary <- rbind(
"Core to core" = cell_summary(within_core, theta_ordered["core", "core"]),
"Core to periphery" = cell_summary(core_periphery, theta_ordered["core", "periphery"]),
"Periphery to periphery" = cell_summary(within_periphery, theta_ordered["periphery", "periphery"])
)
round(block_cell_summary, 3)
```
<div id="gade-block-plots"></div>
```{r}
#| label: gade-block-plots
#| fig-height: 5.2
#| fig-width: 10
#| fig-cap: "Left: the observed cooperation matrix reordered by the fitted roles. Right: the fitted probability of a tie for each role pairing. The small dense core connects widely into a much sparser periphery."
actor_order <- c(core_actors, periphery_actors)
adj_long <- melt(
gade_A[actor_order, actor_order],
varnames = c("organization_1", "organization_2"),
value.name = "tie"
)
adj_long$organization_1 <- factor(adj_long$organization_1, levels = rev(actor_order))
adj_long$organization_2 <- factor(adj_long$organization_2, levels = actor_order)
adj_plot <- ggplot(adj_long, aes(organization_2, organization_1, fill = factor(tie))) +
geom_tile(color = "white", linewidth = .15) +
geom_vline(xintercept = length(core_actors) + .5, color = "#7BBD00", linewidth = 1) +
geom_hline(yintercept = length(periphery_actors) + .5, color = "#7BBD00", linewidth = 1) +
scale_fill_manual(values = c("0" = "white", "1" = "#18453B"), guide = "none") +
coord_equal() +
labs(title = "Observed Ties, Reordered by Role", x = NULL, y = NULL) +
theme(axis.text = element_blank(), axis.ticks = element_blank())
theta_long <- melt(theta_ordered, varnames = c("source_role", "target_role"), value.name = "probability")
theta_plot <- ggplot(theta_long, aes(target_role, source_role, fill = probability)) +
geom_tile(color = "white", linewidth = 1) +
geom_text(aes(label = sprintf("%.2f", probability)), color = "white", fontface = "bold", size = 5) +
scale_fill_gradient(low = "#C7D4C8", high = "#18453B", limits = c(0, 1)) +
coord_equal() +
labs(title = "Fitted Role-to-Role Probabilities", x = "Partner Role", y = "Organization Role", fill = "Probability")
adj_plot + theta_plot
```
Al-Nusrah Front (ANF) and Ahrar al-Sham Islamic Movement (ASIM) form the two-organization core. The fitted within-core probability is about `r round(theta_ordered["core", "core"], 2)`, but that number is supported by only one possible within-core pair, and ANF-ASIM is tied. The stronger evidence is the much broader cross-cell pattern: the two core organizations connect to `r block_cell_summary["Core to periphery", "observed_ties"]` of `r block_cell_summary["Core to periphery", "possible_pairs"]` possible core-periphery pairs, while peripheral organizations connect to only `r block_cell_summary["Periphery to periphery", "observed_ties"]` of `r block_cell_summary["Periphery to periphery", "possible_pairs"]` possible pairs among themselves.
Al-Nusrah Front and Ahrar al-Sham Islamic Movement had different Islamist projects, but both worked across much of this observed cooperation network. They carried out an operation together and each recorded operations with many organizations that rarely cooperated among themselves. That does not make them a unified coalition or imply shared long-term goals. The model is describing a small tactical backbone, not two isolated ideological camps.
The standard SBM has no separate parameter for each organization's number of partners. It may therefore be grouping organizations mainly by how broadly they cooperate, turning the continuous activity difference from the SRM into two categories. ANF and ASIM appearing prominently in both models is not separate proof that they form a natural type. The useful cross-model finding is simpler: these two organizations have unusually broad tactical reach, while the exact description of that reach depends on the model.
## Check the Nearby Three-Block Solution {#sec-block-sensitivity}
```{r}
#| label: gade-k3
z_gade3 <- apply(gade_sbm$memberships[[3]]$Z, 1, which.max)
names(z_gade3) <- rownames(gade_A)
theta_gade3 <- gade_sbm$model_parameters[[3]]$pi
core_block3 <- which.max(rowMeans(theta_gade3))
list(core = names(z_gade3)[z_gade3 == core_block3], other_blocks = split(names(z_gade3)[z_gade3 != core_block3], z_gade3[z_gade3 != core_block3]))
round(theta_gade3, 2)
```
ANF and ASIM remain together in the core. The additional block divides the periphery into organizations with moderate access to the core and organizations with weaker access. That is the stable result: the two-organization core persists, while the exact number of peripheral types is less certain.
## Simulate Networks From the Fitted Roles {#sec-block-gof}
```{r}
#| label: block-gof
calc_block_gof <- function(M) c(degree_sd = sd(rowSums(M)), triangles = sum(diag(M %*% M %*% M)) / 6)
sim_gade_sbm <- function() {
P <- theta_gade[z_gade, z_gade]
M <- matrix(rbinom(length(P), 1, P), nrow(P), ncol(P))
M[lower.tri(M)] <- t(M)[lower.tri(M)]
diag(M) <- 0
M
}
set.seed(6886)
gade_obs_gof <- calc_block_gof(gade_A)
gade_sim_gof <- t(replicate(500, calc_block_gof(sim_gade_sbm())))
round(rbind(observed = gade_obs_gof, simulated_mean = colMeans(gade_sim_gof), standardized_gap = abs(gade_obs_gof - colMeans(gade_sim_gof)) / apply(gade_sim_gof, 2, sd)), 2)
```
We now ask a concrete question: if this two-role description were an adequate summary, what would networks generated from it look like? We keep each organization's fitted role and the three fitted core-to-core, core-to-periphery, and periphery-to-periphery probabilities. We then generate 500 new networks from those probabilities and compare their partner counts and triangles with the observed network.
The generated networks are a little too uniform in how many partners organizations have, and they contain far fewer triangles than the observed network. The model gets the large pattern right, a tiny well-connected core and a sparse periphery, but it treats all organizations inside a role as too interchangeable. There is still smaller-scale coalition structure that two boxes do not capture.
There is also a data-construction reason to be cautious about the triangle gap. If one recorded operation involved three organizations, the data create all three pairwise ties at once. That single event mechanically creates a triangle. The missing triangles can therefore reflect omitted coalition detail, differences among organizations in the same role, or the conversion of group operations into pairs. They are not clean evidence that one partnership caused another through friend-of-a-friend closure.
# One Binary Panel, Two Longitudinal Questions {#sec-longitudinal}
The cross-sectional Syrian analysis asks what organizes one observed network. We now use 13 annual Integrated Crisis Early Warning System (ICEWS) networks from 2002 through 2014 and fit two longitudinal models to the same binary outcome.
- The longitudinal SRM asks which states repeatedly appear on the source and target sides of high-volume coded material conflict after regime difference and year-to-year density changes enter.
- The repeated-network NetMix model asks how states' estimated relational profiles differ across the annual networks.
ICEWS records coded events with a source actor and a target actor (Boschee et al. 2015). Our directed outcome equals one when ICEWS records more than 20 material-conflict events with state $i$ as the source and state $j$ as the target in a year. This is a relatively high event-volume threshold for our classroom panel. It is not a standard definition of war intensity, and the source-target distinction does not always correspond to military initiation and victimization. Event volume can also reflect sustained activity, news visibility, and source coverage. We keep the threshold fixed across the two main models, then check nearby thresholds after the SRM fit.
```{r}
#| label: netmix-data
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)
states
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)
netmix_dyad <- icews_sub[, c("i", "j", "year", "matlConf", "high_event_volume", "polity_gap", "same_region")]
names(netmix_dyad)[1:4] <- c("source", "target", "year", "material_conflict_events")
netmix_dyad$year <- as.integer(netmix_dyad$year)
netmix_monad <- unique(icews_sub[, c("i", "year")])
names(netmix_monad) <- c("state", "year")
netmix_monad$year <- as.integer(netmix_monad$year)
checkpoint(dyad_years = nrow(netmix_dyad), state_years = nrow(netmix_monad), years = length(unique(netmix_dyad$year)), states = length(states), above_threshold_rate = round(mean(netmix_dyad$high_event_volume), 3))
```
The 18 states are the highest-volume ICEWS sources among states with complete Polity and GDP coverage for all 13 years. Polity runs from -10 for a strongly autocratic regime to 10 for a strongly democratic regime. We use the absolute difference between two states' scores and divide it by 20, so the model variable runs from 0 to 1. A 10-point difference on the original scale therefore enters the model as 0.5. We require complete Polity and GDP coverage to preserve one fixed case set for optional predictor extensions, although GDP does not enter the displayed SRM. GDP completeness is a classroom design choice, not an estimator requirement. This purposive, outcome-selected case set gives the dynamic SRM enough above-threshold relationships and changes to learn from. It is not a probability sample of the international system, and its event rate should not be generalized to all states. The NetMix application later uses a separate defense-alliance panel.
# A Fast Binary Longitudinal SRM {#sec-lame-long}
The directed panel lets the two actor scores change over time. For state $i$ acting as the source, state $j$ acting as the target, and year $t$, the binary probit SRM is
$$
Y^*_{ijt}=\alpha_t+\mathbf{x}_{ijt}^{\mathsf T}\boldsymbol\beta+a_{it}+b_{jt}+e_{ijt},
\qquad
Y_{ijt}=\mathbb{1}(Y^*_{ijt}>0).
$$
Read the pieces from left to right. $\alpha_t$ gives each year its own baseline event rate. $\mathbf{x}_{ijt}^{\mathsf T}\boldsymbol\beta$ contains the measured pair-year predictors and their coefficients. $a_{it}$ asks whether state $i$ appears as the source of high-volume event relationships more often than those measured pieces predict in year $t$. $b_{jt}$ asks the corresponding target-side question for state $j$. The indicator turns the unobserved probit score into the observed zero or one.
The yearly actor scores are connected rather than estimated as 13 unrelated snapshots:
$$
a_{it}=\rho_{ab}a_{i,t-1}+\nu_{it},
\qquad
b_{it}=\rho_{ab}b_{i,t-1}+\xi_{it}.
$$
The shared $\rho_{ab}$ term links each annual score to the previous year's score, while the innovations $\nu_{it}$ and $\xi_{it}$ allow the patterns of events coded from or toward a state to move. In the fast ALS fit below, $\rho_{ab}$ is a fixed smoothing control rather than a persistence coefficient estimated from these data. Put simply, the model lets a state's broad involvement change gradually across 2002 through 2014 instead of assigning one permanent score.
## Let `netify` Build the Longitudinal Matrices
`netify()` already knows how to turn a time-stamped directed edge table into aligned yearly matrices. We add the year indicators as ordinary dyadic predictors, build the object once, and let `to_lame()` produce the outcome and predictor lists.
```{r}
#| label: lame-long-data
years <- sort(unique(netmix_dyad$year))
year_predictors <- paste0("year_", years[-1])
for (yy in years[-1]) {
netmix_dyad[[paste0("year_", yy)]] <- as.integer(netmix_dyad$year == yy)
}
icews_predictors <- c("polity_gap", "same_region", year_predictors)
icews_net <- netify(
input = netmix_dyad,
actor1 = "source",
actor2 = "target",
time = "year",
symmetric = FALSE,
weight = "high_event_volume",
dyad_vars = icews_predictors,
dyad_vars_symmetric = rep(TRUE, length(icews_predictors)),
missing_to_zero = TRUE
)
icews_lame <- to_lame(
icews_net,
lame = TRUE,
family = "binary",
fit_method = "als",
bootstrap = 100
)
Y_icews <- icews_lame$Y
X_icews <- icews_lame$Xdyad
checkpoint(
yearly_outcome_matrices = length(Y_icews),
actors_per_matrix = nrow(Y_icews[[1]]),
predictors = dim(X_icews[[1]])[3]
)
```
## Use the Fast Point Estimator
<div id="lame-long-fit"></div>
```{r}
#| label: lame-long-fit
icews_srm_als <- cache_fit(
"icews_binary_dynamic_srm_boot100_v1",
{
fit <- lame(
Y = Y_icews,
Xdyad = X_icews,
family = "binary",
symmetric = FALSE,
R = 0,
dynamic_ab = TRUE,
method = "als",
als_max_iter = 500,
bootstrap = 100,
bootstrap_type = "parametric",
bootstrap_seed = 6886,
verbose = FALSE
)
fit$bootstrap_dynamic$a_paths <- simplify2array(lapply(
fit$bootstrap_dynamic$refits,
function(refit) refit$a_dynamic
))
fit$bootstrap_dynamic$b_paths <- simplify2array(lapply(
fit$bootstrap_dynamic$refits,
function(refit) refit$b_dynamic
))
fit$bootstrap_dynamic$refits <- NULL
fit
}
)
checkpoint(
converged = icews_srm_als$converged,
iterations = icews_srm_als$iterations,
successful_bootstrap_refits = icews_srm_als$bootstrap_dynamic$n_success
)
```
`method = "als"` is the fast point-estimation route. At a high level, the model first converts each binary outcome into a temporary continuous working score on the probit scale. It then finds coefficients and annual source and target paths that reconstruct those working scores while discouraging implausibly jagged year-to-year paths.
The simplified objective is
$$
\min_{\boldsymbol\beta,\mathbf a,\mathbf b}
\sum_{t}\sum_{i\neq j}w_{ijt}
\left[z_{ijt}-\left(\mathbf x_{ijt}^{\mathsf T}\boldsymbol\beta+a_{it}+b_{jt}\right)\right]^2
+\lambda_a\sum_{t>1}\lVert\mathbf a_t-\rho_{ab}\mathbf a_{t-1}\rVert^2
+\lambda_b\sum_{t>1}\lVert\mathbf b_t-\rho_{ab}\mathbf b_{t-1}\rVert^2.
$$
Here $z_{ijt}$ and $w_{ijt}$ are the current probit working response and weight. The first term rewards reconstruction of the observed zeroes and ones on that working scale. The last two terms are smoothing penalties that connect adjacent source and target scores. The fit alternates among updating the working response, solving for $\boldsymbol\beta$, updating the actor paths, and recomputing the weights until the penalized objective stops changing. With $R=0$, there is no multiplicative factor block. This is a penalized point estimate, not a posterior sampler and not the maximizer of the exact binary likelihood.
The 100 parametric bootstrap replicates add uncertainty intervals to the fast fit. Each replicate generates a new panel from the fitted SRM, refits the same model, and records the coefficients and actor paths. The intervals below show how much those estimates vary under repeated panels generated by this fitted model. They are not Bayesian credible intervals, they depend on the SRM being a reasonable data-generating approximation, and 100 replicates are a classroom compromise rather than a final publication run.
## Translate the Measured Associations Into Probabilities {#sec-lame-results}
```{r}
#| label: lame-long-coefs
coef_boot <- icews_srm_als$bootstrap_dynamic$coef_paths[, 1, , drop = FALSE]
coef_interval <- function(term) {
values <- coef_boot[term, 1, ]
c(
estimate = icews_srm_als$coefficients[term],
lower = quantile(values, .025),
upper = quantile(values, .975)
)
}
srm_coef_summary <- rbind(
"Polity gap, 0 to 20 points" = coef_interval("polity_gap_dyad"),
"Same region" = coef_interval("same_region_dyad")
)
round(srm_coef_summary, 3)
```
The intervals for both measured associations stay above zero in this selected panel. Larger regime differences and being in the same region are associated with a higher probability of crossing the event-volume threshold after yearly baselines and the changing source and target scores enter. They remain associations. The selected states, event coding, news coverage, and possible correlation between the actor effects and predictors all limit stronger claims.
```{r}
#| label: lame-long-probabilities
year_index <- match(netmix_dyad$year, years)
source_index <- match(netmix_dyad$source, rownames(icews_srm_als$a_dynamic))
target_index <- match(netmix_dyad$target, rownames(icews_srm_als$b_dynamic))
eta_actor_year <- icews_srm_als$coefficients["intercept"] +
icews_srm_als$a_dynamic[cbind(source_index, year_index)] +
icews_srm_als$b_dynamic[cbind(target_index, year_index)]
for (yy in years[-1]) {
eta_actor_year <- eta_actor_year +
icews_srm_als$coefficients[paste0("year_", yy, "_dyad")] *
(netmix_dyad$year == yy)
}
p_same0 <- pnorm(
eta_actor_year +
icews_srm_als$coefficients["polity_gap_dyad"] * netmix_dyad$polity_gap
)
p_same1 <- pnorm(
eta_actor_year +
icews_srm_als$coefficients["polity_gap_dyad"] * netmix_dyad$polity_gap +
icews_srm_als$coefficients["same_region_dyad"]
)
p_gap0 <- pnorm(
eta_actor_year +
icews_srm_als$coefficients["same_region_dyad"] * netmix_dyad$same_region
)
p_gap10 <- pnorm(
eta_actor_year +
.5 * icews_srm_als$coefficients["polity_gap_dyad"] +
icews_srm_als$coefficients["same_region_dyad"] * netmix_dyad$same_region
)
srm_probability_contrasts <- data.frame(
comparison = c(
"Different region to same region",
"Polity gap of 0 to 10 points"
),
probability_before = c(mean(p_same0), mean(p_gap0)),
probability_after = c(mean(p_same1), mean(p_gap10))
)
srm_probability_contrasts$difference <-
srm_probability_contrasts$probability_after -
srm_probability_contrasts$probability_before
srm_probability_contrasts[-1] <- round(100 * srm_probability_contrasts[-1], 1)
srm_probability_contrasts
```
These comparisons score the same state-pair-years twice while keeping their fitted year, source, and target pieces fixed. They translate a probit coefficient into an average fitted probability difference. They are not observed before-and-after changes or causal effects.
## Show the Changing Source and Target Scores {#sec-lame-trajectories}
```{r}
#| label: lame-long-trajectories
#| fig-height: 5.8
#| fig-width: 12
#| fig-cap: "Yearly source scores and target scores for four states. The lines are the dynamic ALS estimates drawn with `lame::ab_plot()`; the shaded areas are 95% parametric-bootstrap intervals from 100 model refits."
#| code-fold: true
a_boot <- icews_srm_als$bootstrap_dynamic$a_paths
b_boot <- icews_srm_als$bootstrap_dynamic$b_paths
focus_states <- c(
"Syrian Arab Republic",
"Russian Federation",
"Iran, Islamic Republic Of",
"United States"
)
focus_labels <- c(
"Syrian Arab Republic" = "Syria",
"Russian Federation" = "Russia",
"Iran, Islamic Republic Of" = "Iran",
"United States" = "United States"
)
actor_colors <- c(
"Syria" = "#18453B",
"Russia" = "#008208",
"Iran" = "#7BBD00",
"United States" = "#535054"
)
actor_linetypes <- c(
"Syria" = "solid",
"Russia" = "dashed",
"Iran" = "dotdash",
"United States" = "longdash"
)
dynamic_ab_plot <- function(fit, effect, boot, title, y_label) {
plot_fit <- fit
class(plot_fit) <- unique(c(
"lame",
setdiff(class(plot_fit), c("ame_als", "lame_als"))
))
p <- lame::ab_plot(
plot_fit,
effect = effect,
plot_type = "trajectory",
show_actors = focus_states,
title = title
)
p$data$actor <- factor(
unname(focus_labels[as.character(p$data$actor)]),
levels = unname(focus_labels[focus_states])
)
interval_data <- expand.grid(
actor = rownames(if (effect == "sender") fit$a_dynamic else fit$b_dynamic),
time = seq_along(years),
KEEP.OUT.ATTRS = FALSE
)
interval_data$lower <- c(apply(boot, c(1, 2), quantile, .025))
interval_data$upper <- c(apply(boot, c(1, 2), quantile, .975))
interval_data <- interval_data[interval_data$actor %in% focus_states, ]
interval_data$actor <- factor(
unname(focus_labels[as.character(interval_data$actor)]),
levels = unname(focus_labels[focus_states])
)
interval_data <- interval_data[order(interval_data$actor, interval_data$time), ]
ribbon_layer <- geom_ribbon(
data = interval_data,
aes(
x = time,
ymin = lower,
ymax = upper,
group = actor,
fill = actor
),
inherit.aes = FALSE,
stat = "identity",
alpha = .12,
color = NA,
show.legend = FALSE
)
p$layers <- append(p$layers, list(ribbon_layer), after = 0L)
p +
aes(linetype = actor) +
scale_x_continuous(
breaks = seq(1, length(years), by = 3),
labels = years[seq(1, length(years), by = 3)]
) +
scale_color_manual(values = actor_colors, name = NULL) +
scale_fill_manual(values = actor_colors, name = NULL) +
scale_linetype_manual(values = actor_linetypes, name = NULL) +
labs(x = NULL, y = y_label) +
theme(
legend.position = "bottom",
panel.grid.minor.x = element_blank(),
plot.title = element_text(size = 13, face = "bold")
)
}
from_state_plot <- dynamic_ab_plot(
icews_srm_als,
effect = "sender",
boot = a_boot,
title = "Events coded from each state",
y_label = "Source score (a)"
)
toward_state_plot <- dynamic_ab_plot(
icews_srm_als,
effect = "receiver",
boot = b_boot,
title = "Events coded toward each state",
y_label = "Target score (b)"
)
from_state_plot + toward_state_plot +
plot_layout(guides = "collect") &
theme(legend.position = "bottom")
```
Read the left panel as follows: a higher source score means that events above the threshold were coded from that state toward more other states than the measured predictors, the year, and the other states' target scores would lead us to expect. Read the right panel the same way for events coded toward the state. Zero is the model's average after those other pieces are taken into account. The shaded interval shows how much the estimated path changes across 100 new networks generated from the fitted SRM and then refitted. These paths describe where coded events were directed. They do not by themselves tell us who initiated a conflict, who was harmed, or how severe the events were.
::: {.callout-warning}
## Scope of the Dynamic SRM
The source and target scores follow smoothed yearly paths, but the model still has no dyad-specific long-run baseline, lagged outcome, or leftover serial-dependence process for the same pair. The bootstrap reflects uncertainty under the fitted SRM and does not repair case selection, measurement error, or correlation between the random effects and observed predictors. A Mundlak specification such as the one in @sec-inference addresses one particular random-effects correlation problem when suitable time-varying covariates are available.
:::
## Does the Event-Volume Threshold Drive the Result? {#sec-threshold-sensitivity}
```{r}
#| label: threshold-sensitivity
fit_threshold_srm <- function(cutoff) {
threshold_dyad <- netmix_dyad
threshold_dyad$threshold_outcome <- as.integer(
threshold_dyad$material_conflict_events > cutoff
)
threshold_net <- netify(
input = threshold_dyad,
actor1 = "source",
actor2 = "target",
time = "year",
symmetric = FALSE,
weight = "threshold_outcome",
dyad_vars = icews_predictors,
dyad_vars_symmetric = rep(TRUE, length(icews_predictors)),
missing_to_zero = TRUE
)
stopifnot(all(unlist(validate_netify(threshold_net))))
threshold_lame <- to_lame(
threshold_net,
lame = TRUE,
family = "binary",
fit_method = "als"
)
cache_fit(paste0("icews_binary_dynamic_srm_volume", cutoff, "_v1"), lame(
Y = threshold_lame$Y,
Xdyad = threshold_lame$Xdyad,
family = "binary", symmetric = FALSE, R = 0,
dynamic_ab = TRUE,
method = "als", als_max_iter = 500,
verbose = FALSE
))
}
threshold_fits <- list("10" = fit_threshold_srm(10), "20" = icews_srm_als, "30" = fit_threshold_srm(30))
threshold_summary <- do.call(rbind, lapply(names(threshold_fits), function(cutoff) {
fit <- threshold_fits[[cutoff]]
data.frame(
threshold = as.integer(cutoff),
outcome_rate = mean(netmix_dyad$material_conflict_events > as.integer(cutoff)),
polity_gap = unname(fit$coefficients["polity_gap_dyad"]),
same_region = unname(fit$coefficients["same_region_dyad"]),
actor_correlation = unname(fit$VC["cab"] / sqrt(fit$VC["va"] * fit$VC["vb"])),
row.names = NULL
)
}))
threshold_summary[-1] <- round(threshold_summary[-1], 3)
threshold_summary
```
Moving the threshold from 10 to 20 to 30 coded events changes the number of positive dyad-years and the magnitude of the point estimates, as it should. The positive Polity-gap and same-region patterns, along with a strong source-target actor correlation, survive all three choices. This does not make the threshold correct or the associations causal. It shows that the main classroom interpretation is not an artifact of exactly 20 events.
# Did Post-Cold War Defense Partners Follow the Same Path? {#sec-netmix-long}
The longitudinal question now changes. The SRM asked which states appeared in many high-volume ICEWS relationships. A role model should tell us more than who has many ties, so we switch to defense alliances and ask: did Poland, Romania, and Russia develop similar sets of defense partners after the Cold War, or did their treaty networks take different paths? We also ask which states linked partner systems that would otherwise remain separate.
## The Defense-Alliance Panel
The panel contains 50 states observed yearly from 1991 through 2000. The outcome is one when a pair had a defense commitment in that year and zero otherwise. The data come from the international alliance data assembled for Cranmer, Desmarais, and Menninga (2012) and Cranmer, Desmarais, and Kirkland (2012). This is an undirected network because a recorded defense commitment connects the pair rather than running from one state toward the other.
```{r}
#| label: netmix-alliance-data
alliance_panel <- readRDS("data/alliance_panel_1991_2000.rds")
alliance_dyad <- alliance_panel$dyads[, c(
"state1", "state2", "year", "defense_alliance"
)]
names(alliance_dyad)[1:2] <- c("state_a", "state_b")
alliance_monad <- expand.grid(
state = alliance_panel$state_order,
year = alliance_panel$years,
stringsAsFactors = FALSE
)
alliance_net <- netify(
input = alliance_dyad,
actor1 = "state_a",
actor2 = "state_b",
time = "year",
symmetric = TRUE,
weight = "defense_alliance",
missing_to_zero = FALSE
)
stopifnot(all(unlist(validate_netify(alliance_net))))
checkpoint(
pair_years = nrow(alliance_dyad),
years = n_periods(alliance_net),
states = n_actors(alliance_net),
alliance_rate = round(mean(alliance_dyad$defense_alliance), 3)
)
```
One row is an unordered state-pair-year. `netify()` verifies the actor roster, annual slices, symmetry, and missingness before estimation. NetMix requires its own long dyad and monad tables, so we pass the same validated source tables directly rather than rebuilding matrices by hand. A tie records a formal defense commitment, not how credible it was, how closely the states coordinated, or whether the treaty would have been honored in a crisis.
## What NetMix Is Trying to Learn {#sec-netmix-fit}
Start with three blank partner-list patterns. NetMix learns which patterns tend to contain alliances within or across them and how much each state-year's observed set of defense partners resembles each pattern. The result is not a permanent assignment. Every state-year receives three weights, $\boldsymbol\pi_{it}$, that add to one.
For a pair of states $i$ and $j$ in year $t$, imagine that the model temporarily assigns each endpoint a role using those weights:
$$
g_{ijt}\sim\operatorname{Categorical}(\boldsymbol\pi_{it}),
\qquad
h_{ijt}\sim\operatorname{Categorical}(\boldsymbol\pi_{jt}),
\qquad
Y_{ijt}\mid g_{ijt},h_{ijt}\sim\operatorname{Bernoulli}(\theta_{g_{ijt},h_{ijt}}).
$$
`Categorical` simply means “choose one of the three roles using the state's three weights.” The number $\theta_{gh}$ is the fitted chance of an alliance for a pair whose endpoints fully match roles $g$ and $h$. Because observed states can mix roles, their fitted alliance probability averages all nine role pairings:
$$
\Pr(Y_{ijt}=1\mid\boldsymbol\pi_{it},\boldsymbol\pi_{jt},\boldsymbol\Theta)
=\sum_g\sum_h\pi_{it,g}\pi_{jt,h}\theta_{gh}
=\boldsymbol\pi_{it}^{\mathsf T}\boldsymbol\Theta\boldsymbol\pi_{jt}.
$$
In plain language, the model asks how much each state resembles each alliance pattern, looks up how strongly those patterns connect, and combines the answers.
```{r}
#| label: netmix-alliance-load
# The downloaded folder contains the fitted object and start-value audit.
load("cache/netmix_alliance_k3_v1.rda")
netmix_long <- netmix_alliance
c(
converged = netmix_long$converged,
iterations = netmix_long$niter,
final_elbo = netmix_long$LowerBound
)
summary(netmix_alliance_starts$elbo)
```
The code below shows the start-value audit that produced the saved fit. It is not evaluated during rendering because the fitted objects are included with the handout.
```{r}
#| label: netmix-alliance-call
#| eval: false
set.seed(6886)
netmix_candidates <- lapply(1:12, function(start) {
mm_start <- netmix_alliance_init_mm +
matrix(
runif(length(netmix_alliance_init_mm), 0, .03),
nrow = nrow(netmix_alliance_init_mm)
)
mm_start <- sweep(mm_start, 2, colSums(mm_start), "/")
block_start <- netmix_alliance_init_block +
matrix(rnorm(9, 0, .08), 3, 3)
block_start <- (block_start + t(block_start)) / 2
NetMix::mmsbm(
defense_alliance ~ 1,
~ 1,
senderID = "state_a",
receiverID = "state_b",
nodeID = "state",
timeID = "year",
data.dyad = alliance_dyad,
data.monad = alliance_monad,
n.blocks = 3,
n.hmmstates = 1,
directed = FALSE,
mmsbm.control = list(
seed = 6800 + start,
nstart = 1,
svi = FALSE,
vi_iter = 15000,
conv_tol = 1e-4,
hessian = FALSE,
assortative = FALSE,
mu_block = c(-2.5, -1.5),
mm_init_t = mm_start,
b_init_t = block_start,
verbose = FALSE
)
)
})
netmix_long <- netmix_candidates[[
which.max(vapply(
netmix_candidates,
function(fit) fit$LowerBound,
numeric(1)
))
]]
```
We set the number of roles to three, fit without predictors, and use full-batch variational updates. Twelve nearby starting values recover the same three-role pattern, although broader starts can reach less useful local solutions. That is a reason to inspect the role matrix and rerun the optimizer, not a reason to report whichever run gives the easiest story. `hessian = FALSE` makes this a descriptive point fit without standard errors.
The argument `n.hmmstates = 1` means this classroom fit does not estimate a hidden transition process or smooth memberships across years. It estimates a separate role mixture for every state-year and one role matrix shared across the decade. We can compare annual mixtures, but we should not call them estimated transitions.
## How the Fast Estimator Works
At a high level, the estimator asks whether each state-year's entire treaty list looks more like the regional, dense-multilateral, or sparse pattern, then revises the alliance rates among those patterns. Poland's 1997 role mixture is supported by every alliance and non-alliance in Poland's 1997 row, not by one treaty or its degree alone.
Technically, NetMix uses variational inference for the mixed-membership blockmodel. The latent quantities include each dyad endpoint's temporary role assignment and the state-year membership vectors that generate those assignments. The algorithm chooses a tractable distribution $q$ and maximizes
$$
\operatorname{ELBO}(q,\boldsymbol\Theta)
=E_q[\log p(Y,\text{roles},\boldsymbol\pi\mid\boldsymbol\Theta)]
-E_q[\log q(\text{roles},\boldsymbol\pi)].
$$
Local updates revise the approximate role assignment for each dyad endpoint and the membership vector for each state-year. Global updates revise the role-pair log-odds in $\boldsymbol\Theta$ and any regression coefficients. This fit uses full-batch updates, so every iteration uses the whole panel rather than a stochastic subset. It keeps alternating until the ELBO changes by less than the declared tolerance.
The evidence lower bound, or ELBO, is the score used to track that process. A larger ELBO is better when comparing runs of the same specification. It is not a probability, a goodness-of-fit test, or proof that three roles are the truth. We will discuss the related back-and-forth logic of ALS in more detail on Day 10.
## Read the Relationship Pattern Before Naming the Roles {#sec-netmix-results}
```{r}
#| label: netmix-alliance-theta
#| fig-height: 5
#| fig-width: 7
#| fig-cap: "Estimated alliance probability for each pair of roles. The values, not color alone, show one densely connected pattern, one moderately connected regional pattern, and a sparse pattern."
#| fig-alt: "A three-by-three heat map. The regional-pattern diagonal cell is about 0.48, the dense-multilateral diagonal cell is about 0.94, and all cells involving the sparse-or-outside role are near zero."
netmix_role_names <- c(
"Regional alliance pattern",
"Dense multilateral pattern",
"Sparse or outside pattern"
)
netmix_theta <- plogis(netmix_long$BlockModel)
dimnames(netmix_theta) <- list(
role_a = netmix_role_names,
role_b = netmix_role_names
)
theta_plot <- as.data.frame(as.table(netmix_theta))
names(theta_plot) <- c("role_a", "role_b", "probability")
ggplot(theta_plot, aes(role_b, role_a, fill = probability)) +
geom_tile(color = "white", linewidth = 1) +
geom_text(
aes(label = scales::percent(probability, accuracy = 1)),
color = ifelse(theta_plot$probability > .35, "white", "black"),
fontface = "bold"
) +
scale_fill_gradient(
low = "#F2F2F2",
high = "#18453B",
labels = scales::label_percent(),
limits = c(0, 1)
) +
labs(
x = NULL,
y = NULL,
fill = "Alliance chance"
) +
theme(
axis.text.x = element_text(angle = 25, hjust = 1),
panel.grid = element_blank()
)
```
The first thing to learn is where the recorded treaties go. State-years whose partner lists match the dense multilateral pattern are almost always allied with one another. State-years matching the regional pattern are also often allied with one another, but less consistently. Alliances across these patterns are rare in this panel. The sparse pattern contains state-years with few recorded defense partners among the 50 states included here.
We gave the roles readable names only after inspecting the fitted matrix and the states with the largest weights. The labels describe alliance portfolios in this panel. They are not claims about common ideology, shared interests, or a state's entire foreign policy.
## See Mixtures and Change Directly {#sec-netmix-profiles}
```{r}
#| label: netmix-alliance-memberships
#| fig-height: 7
#| fig-width: 10
#| fig-cap: "Yearly NetMix role weights for four states. Color, line type, and point shape all distinguish the three fitted alliance patterns."
#| fig-alt: "Four small line charts for the United States, Poland, Romania, and Russia from 1991 to 2000. The United States mixes regional and dense-multilateral roles early and becomes mostly dense-multilateral. Poland moves sharply from the sparse role to the dense-multilateral role in 1997. Romania remains mostly sparse. Russia moves from sparse to a mix dominated by the regional role."
membership_long <- data.frame(
key = colnames(netmix_long$MixedMembership),
t(netmix_long$MixedMembership),
check.names = FALSE
)
names(membership_long)[2:4] <- netmix_role_names
membership_long$state <- sub("@[0-9]+$", "", membership_long$key)
membership_long$year <- as.integer(sub("^.*@", "", membership_long$key))
membership_long$state_name <- unname(
alliance_panel$state_labels[membership_long$state]
)
alliance_degree <- do.call(
rbind,
lapply(alliance_panel$years, function(yy) {
tied <- alliance_dyad[
alliance_dyad$year == yy &
alliance_dyad$defense_alliance == 1,
]
degree <- table(factor(
c(tied$state_a, tied$state_b),
levels = alliance_panel$state_order
))
data.frame(
state = alliance_panel$state_order,
year = yy,
alliance_count = as.numeric(degree)
)
})
)
membership_long <- merge(
membership_long,
alliance_degree,
by = c("state", "year"),
all.x = TRUE
)
dense_degree_correlation <- cor(
membership_long[["Dense multilateral pattern"]],
membership_long$alliance_count
)
focus_states <- c(
"United States", "Poland", "Romania", "Russia"
)
profile_long <- melt(
membership_long[
membership_long$state_name %in% focus_states,
],
id.vars = c("state_name", "year"),
measure.vars = netmix_role_names,
variable.name = "role",
value.name = "weight"
)
profile_long$state_name <- factor(
profile_long$state_name,
levels = focus_states
)
role_colors <- c(
"Regional alliance pattern" = "#18453B",
"Dense multilateral pattern" = "#7BBD00",
"Sparse or outside pattern" = "#535054"
)
role_linetypes <- c(
"Regional alliance pattern" = "solid",
"Dense multilateral pattern" = "longdash",
"Sparse or outside pattern" = "dotted"
)
role_shapes <- c(
"Regional alliance pattern" = 16,
"Dense multilateral pattern" = 17,
"Sparse or outside pattern" = 15
)
ggplot(
profile_long,
aes(
year,
weight,
color = role,
linetype = role,
shape = role
)
) +
geom_line(linewidth = .9) +
geom_point(size = 1.8) +
facet_wrap(~ state_name, ncol = 2) +
scale_color_manual(values = role_colors) +
scale_linetype_manual(values = role_linetypes) +
scale_shape_manual(values = role_shapes) +
scale_x_continuous(
breaks = c(1991, 1994, 1997, 2000)
) +
scale_y_continuous(
limits = c(0, 1),
labels = scales::label_percent()
) +
labs(
x = NULL,
y = "Share of fitted role profile",
color = NULL,
linetype = NULL,
shape = NULL
) +
theme(legend.position = "bottom")
```
The United States begins with a partner list that partly resembles both connected patterns. From 1992 onward, its recorded defense partners overwhelmingly match the dense multilateral system. Poland provides the clearest later change: through 1996 it has few ties into either connected pattern, while from 1997 onward its partners are mainly states already tied into the dense multilateral system. Romania does not make the same move during this window. Russia instead develops a partner list concentrated in the more moderately connected regional system.
The result is partly, but not entirely, a summary of alliance counts. Dense-pattern membership correlates about `r round(dense_degree_correlation, 2)` with the number of alliances a state-year has. The partner identities still matter: Poland and Romania can have similarly short early partner lists but move differently when the identities of their partners diverge. The finding supports a comparison of those paths, not a claim that the fitted roles are permanent types of states.
## Does the Role Model Reconstruct the Alliances? {#sec-netmix-calibration}
```{r}
#| label: netmix-alliance-calibration
netmix_mm <- t(netmix_long$MixedMembership)
netmix_dd <- netmix_long$dyadic.data
state_a_key <- paste(
netmix_dd[["(sid)"]],
netmix_dd[["(tid)"]],
sep = "@"
)
state_b_key <- paste(
netmix_dd[["(rid)"]],
netmix_dd[["(tid)"]],
sep = "@"
)
state_a_mm <- netmix_mm[
match(state_a_key, rownames(netmix_mm)),
,
drop = FALSE
]
state_b_mm <- netmix_mm[
match(state_b_key, rownames(netmix_mm)),
,
drop = FALSE
]
netmix_pred <- vapply(
seq_len(nrow(netmix_dd)),
function(ii) {
sum(
outer(state_a_mm[ii, ], state_b_mm[ii, ]) *
plogis(netmix_long$BlockModel)
)
},
numeric(1)
)
brier_model <- mean((netmix_long$Y - netmix_pred)^2)
brier_null <- mean(
(netmix_long$Y - mean(netmix_long$Y))^2
)
round(
c(
role_model = brier_model,
constant_probability = brier_null
),
3
)
```
The Brier score is the average squared gap between the fitted probability and the observed zero or one, so lower is better. The three-role fit scores about `r round(brier_model, 3)`, compared with `r round(brier_null, 3)` if every pair-year receives the same alliance probability. This is a large in-sample reconstruction improvement, but it is not a forecast. The same observations taught us the roles and were used to check them.
::: {.callout-important}
## What the Alliance Records Show
The three countries do not follow one common post-Cold War path in these records. Poland's defense partners shift sharply in 1997 toward states already tied into a dense multilateral alliance system. Romania remains connected to relatively few states in this 50-state panel, while Russia's partners concentrate in a different, more regional system. NetMix makes those differences in partner lists easy to compare from year to year. Much of the dense-pattern score still reflects how many alliances a state has, and the model does not tell us why the treaty lists changed or whether every commitment was equally credible.
:::
::: {.callout-warning}
## What We Can and Cannot Claim
We can describe how each state-year's observed alliance portfolio resembles the fitted roles. We cannot treat the roles as natural kinds, read a weight as the probability that a state “belongs to the West,” claim that this one-hidden-state fit estimated transitions, or infer that a treaty caused a broader realignment. A research analysis should compare role counts and starting values, examine held-out years or pairs, add relevant predictors, and investigate why the portfolio changed.
:::
# How the Three Representations Handle the Missing Pattern {#sec-spine}
| Representation | Latent quantity | Best at | Main limitation |
|---|---|---|---|
| Social relations model (SRM) | One activity score per actor, or separate source-side and target-side scores | Recurring broad actor involvement and reciprocity | Cannot represent special affinity among particular groups of actors |
| Stochastic blockmodel (SBM) | One discrete relational role per actor | A readable role-to-role cooperation pattern | Actors in one role are treated as interchangeable unless the model is extended |
| Repeated-network mixed-membership stochastic blockmodel, fitted with NetMix | A separate blend of roles for every actor-year, with one shared role matrix | Comparing how alliance portfolios line up with recurring relational patterns | With one hidden state, it has no transition process or smoothing across years |
# What Did the Models Show? {#sec-after-fit}
## Cooperation Among Armed Organizations in Syria
> These data record claimed joint operations among armed opposition organizations in Syria from July 2012 through June 2015. Two findings stand out. First, organizations closer to one another in estimated size tended to carry out more joint operations together. Second, some organizations took part in joint operations with far more partners than their measured ideology, size, location, and sponsorship would lead us to expect.
>
> Al-Nusrah Front (ANF) and Ahrar al-Sham Islamic Movement (ASIM) stand out most clearly. They carried out an operation together, and each also carried out operations with many organizations that seldom worked with one another. The blockmodel summarizes this as a two-organization core connected to a much larger set of organizations with few ties among themselves. That is a description of the cooperation network. It does not mean ANF and ASIM had the same goals or formed one lasting coalition.
## Conflictual Events Among States
> Here one case is an ordered pair of states in one year. The outcome equals one when ICEWS records more than 20 materially conflictual events directed from the first state toward the second. It does not indicate whether a war occurred or how severe the conflict was.
>
> After accounting for the year and for each state's recurring place on both sides of these relationships, pairs in the same region crossed the threshold more often than pairs in different regions. Pairs with larger Polity-score differences also crossed it more often. These are adjusted comparisons within this dataset, not evidence that moving a state to another region or changing its regime score would cause more events.
>
> The United States appeared in many above-threshold relationships in both directions throughout 2002 to 2014. Syria appeared in relatively few through 2010. After 2011, Syria appeared in many more, especially as the state toward which events were directed. The bootstrap intervals show that this change is large relative to the uncertainty in the fitted model.
>
> ICEWS is built from events coded from news reports. The patterns can reflect state behavior, what news sources covered, and how the events were coded. They do not measure battlefield severity and they do not establish that region or regime difference caused the events.
## Defense-Alliance Portfolios
> The NetMix application uses a different panel because the ICEWS role fit mostly renamed how many above-threshold partners a state had. The defense-alliance panel gives the role model a clearer job: compare which states are treaty partners, not just how many treaties each state has.
>
> The three countries do not follow one common post-Cold War path in these records. Poland's defense partners shift sharply in 1997 toward states already tied into a dense multilateral alliance system. Romania remains connected to relatively few states in this panel, while Russia's partners concentrate in a different regional system.
>
> These are summaries of observed treaty portfolios. The model does not explain the realignment, judge treaty credibility, or estimate a transition process in this one-hidden-state specification.
# Exercises {#sec-exercises}
1. **Inference:** Let $\beta=0.5$, $\gamma=1$, and $z=\alpha x+e$, where $e$ contains recurring actor effects. Compute the naive target when $\alpha=0.4$, $-0.4$, and $0$. Explain why the zero-correlation case can still have incorrect independent-row uncertainty.
2. **SRM:** Inspect the eight largest and smallest `gade_srm$APM` values. Choose two organizations and write an interpretation that does not turn the latent effect into a fixed trait.
3. **SRM Fit:** Explain in one sentence what the posterior predictive check reproduces and what third-order pattern it misses.
4. **Blocks:** Use the observed numerators and denominators to explain why the cross-cell evidence is stronger than the fitted core-core probability. Then compare the two- and three-block solutions.
5. **Longitudinal SRM:** Choose one state and compare its source-side and target-side paths. Identify a change that is clear relative to its bootstrap interval and one that is not.
6. **NetMix:** Compare the United States, Poland, Romania, and Russia in the membership plot. Explain what the role weights add beyond counting alliances, then identify one claim the one-hidden-state fit cannot support.
# Reading {#sec-reading}
- Back and Kenny (2010), "The Social Relations Model: How to Understand Dyadic Processes," for the core actor, partner, and relationship decomposition.
- Hoff (2005), "Bilinear Mixed-Effects Models for Dyadic Data," for the additive social-relations regression framework introduced here. Day 10 adds the multiplicative term.
- Gade, Gabbay, Hafez, and Kelly (2019), "Networks of Cooperation: Rebel Alliances in Fragmented Civil Wars," for the Syrian cooperation application used here.
- Minhas et al. (2022), "Taking Dyads Seriously," for the simulation evidence, applied reanalyses, latent-effect interpretation, and the correlated-omitted-variable limitation.
- Holland, Laskey, and Leinhardt (1983), "Stochastic Blockmodels: First Steps," for the generative blockmodel.
- Karrer and Newman (2011), "Stochastic Blockmodels and Community Structure in Networks," for degree-corrected blocks.
- Boschee et al. (2015), "ICEWS Coded Event Data," for the source-target event records used in the longitudinal SRM panel.
- Cranmer, Desmarais, and Menninga (2012) and Cranmer, Desmarais, and Kirkland (2012), for the international alliance data used in the NetMix application.
- Olivella, Pratt, and Imai (2022), "Dynamic Stochastic Blockmodel Regression for Network Data," for the broader dynamic mixed-membership framework behind `NetMix`. Our one-hidden-state classroom fit deliberately does not estimate its HMM transition layer.
```{r}
#| label: session
#| code-fold: true
#| code-summary: "Session info"
sessionInfo()
```
