suppressMessages({
library(RSiena)
library(networkdata)
})
data(duqueData)
W <- 3L
actors <- sort(unique(unlist(lapply(dipl_ties[1:W], names))))
n <- length(actors)
dipl <- array(10L, dim = c(n, n, W), dimnames = list(actors, actors, NULL))
for (t in 1:W) {
d <- as.matrix(dipl_ties[[t]])
d[d == "."] <- NA # "." is this dataset's missing code
storage.mode(d) <- "integer"
ids <- names(dipl_ties[[t]])
dimnames(d) <- list(ids, ids)
dipl[ids, ids, t] <- d # NAME-based assignment into the pooled array
}
for (t in 1:W) diag(dipl[, , t]) <- 10L # structural zeros on the diagonal
c(n_actors = n, na_cells = sum(is.na(dipl)))
#> n_actors na_cells
#> 158 586
table(dipl, useNA = "ifany")
#> dipl
#> 0 1 10 <NA>
#> 49948 13536 10822 586A hard SAOM: diplomatic recognition (Duque 2018)
Ancillary self-study · ICPSR Network Analysis: Advanced Topics
ImportantOptional: not covered in the 2026 course
This document is not part of any taught session. Nobody is going to lecture it at you, nobody is going to grade you on it, and you can ignore it completely without falling behind.
It exists because it is a thing you will plausibly need next, and because it is easier to hand you something that works than to hand you a citation and wish you luck. Everything here runs end to end on the same machine setup you used in class. Every number and every figure on this page was produced by a chunk in this file: there are no screenshots of R output anywhere in the ancillary package.
How to work through it alone. Get the source (click the </> Code button at the top right, then View Source, and copy it into a .qmd file: or open the .qmd you were given directly), run the chunks in order, and read the prose between them: the prose is where the argument lives. Slow fits are identified in the prose, and precomputed fits load a shipped result while still showing you the code that produced it. Folded Depth boxes are optional extras; folded Stuck? boxes are hints. Open whichever you need and ignore the other. Exercises have a Solution tab: you are working alone, so use it.
NoteWhat this extends, and why it is here
Extends: the SAOM session (day 12, W3 Tue). The taught session runs on s50, which is a well-behaved dataset: it fits quickly and converges cleanly. This document is the opposite experience on purpose: a larger, directed, IR network that is genuinely hard to fit, and does not fully converge in one pass.
Why that is worth your time. Most of what you will meet in the wild looks like this, not like s50. A document that only ever shows you a SAOM that works teaches you nothing about the far more common situation where it does not, and where the skill is reading the diagnostics honestly rather than reporting a number you should not trust. This is the “what a hard SAOM actually looks like” document.
What you need: RSiena, networkdata (carries the Duque data; no install_github runs here). The fit is slow, so it is precomputed and shipped. The code is shown, and you can run it yourself with options(session.refit = TRUE).
1. The setting
Duque (2018) studies status recognition between states. A state “recognises” another by establishing a diplomatic mission: an embassy. The claim is that recognition is social: you get recognised because of your position in the network of recognition, not only because of your material attributes. That is a dynamic, endogenous story about who recognises whom over time, which is SAOM territory.
The data: dipl_ties, eight waves (1970–2005, every five years) of directed diplomatic ties among states, plus contiguity, alliances, and regime type. We will use the first three waves to keep the fit manageable.
Two features make this hard in ways s50 was not:
- The node set changes: states enter and leave the international system across waves.
- It is directed and dense: 158 states, thousands of ties, and the density plus the directedness make the estimator work much harder than 50 sparse friendships did.
2. Data preparation: where most of the danger lives
This is the part that has to be right, and it is the part that is easy to get silently wrong. The original version of this construction was wrong in three independent ways, none of which threw an error. We will do it correctly and then, in the warning box, name exactly what the three errors were, because they are errors you are likely to make yourself.
The 10 code marks structural zeros: a state that is not in the system at a given wave has all its ties coded 10, which tells RSiena “this actor is not available here” rather than “this actor chose to have no ties.” The scattered NAs (586 of them) are genuine missing data from the "." code.
Confirm the network is still directed after all that reshaping: the single most important check, because the substantive question is entirely about asymmetry:
asym <- sapply(1:W, function(t) {
m <- dipl[, , t]; m[m == 10 | is.na(m)] <- NA
sum(m != t(m), na.rm = TRUE)
})
checkpoint(actors = n, waves = W, asymmetric_cells_per_wave = asym)
#> ------------------------------------------------------------------
#> CHECKPOINT: actors = 158 | waves = 3 | asymmetric_cells_per_wave = 2070, 2050, 2264
#> ------------------------------------------------------------------Thousands of asymmetric cells per wave. If that had come back near zero you would know you had accidentally symmetrised the network, which would silently destroy the “who recognises whom” question the whole study is about.
WarningThe three bugs this construction replaces
The version that circulated in the original teaching materials was eval = FALSE, so it never ran and never errored, and was handed to students as a template. It was wrong three ways:
Positional indexing into a name-dimensioned array. It looped
for (i in 1:nrow(d))and wrotedipl[i, j, t] <- val, butdiplwas dimensioned by the pooled 158-actor list whiledhad only that wave’s actors (134 in wave 1). So from about the third actor onward, every value was written to the wrong cell. The fix is name-based assignment:dipl[ids, ids, t] <- d, which cannot misalign because R matches the names.Silent symmetrisation. It wrote both
dipl[i, j, t] <- valanddipl[j, i, t] <- val, forcing symmetry on a directed network. That corrupts roughly a thousand cells per wave and erases the asymmetry the study depends on. We never write the transpose.Coercion of the missing code without noticing.
as.numeric(as.character(d))turned the"."missing code intoNAwith a warning nobody read: 586 cells. We handle"."explicitly withd[d == "."] <- NAso the missingness is deliberate and countable.
The answer key that shipped alongside the slides actually had the correct name-based form. The slide was a degraded copy of working code. The lesson is the one from the ERGM applications document, restated: index relational data by name, and check dim(), table(), and symmetry after every reshape.
3. The fit, and reading a convergence failure honestly
Build the RSiena data object and a minimal three-effect specification: density, reciprocity, and transitive ties.
dv <- sienaDependent(dipl, allowOnly = FALSE)
dat <- sienaDataCreate(dv)
eff <- getEffects(dat)
eff <- includeEffects(eff, transTies)
#> effectNumber effectName shortName include fix test initialValue parm
#> 1 44 transitive ties transTies TRUE FALSE FALSE 0 0# This is a slow fit. A precomputed result appears below; run it yourself with
# options(session.refit = TRUE).
alg <- sienaAlgorithmCreate(projname = NULL, seed = 6886, n3 = 1000)
ans <- siena07(alg, data = dat, effects = eff,
batch = TRUE, verbose = FALSE, silent = TRUE)The result below is precomputed from the code above.
data.frame(effect = ans$effects$effectName,
est = round(ans$theta, 3),
se = round(sqrt(diag(ans$covtheta)), 3))
#> effect est se
#> 1 outdegree (density) -1.989 0.202
#> 2 reciprocity 2.179 0.035
#> 3 transitive ties 0.891 0.203Substantively this is what the theory predicts: strong reciprocity (recognition is returned) and positive transitivity (states recognise the states their partners recognise). If you stopped here you would report a clean-looking result. Do not stop here.
round(as.numeric(ans$tconv.max), 3)
#> [1] 0.653
ImportantThis did not converge, and that is the point of the document
The overall maximum convergence ratio is well above the 0.25 threshold the RSiena manual sets for publishable results, and above even the loose 0.35 “nearly converged” bar. By the manual’s own standard, you may not report these estimates.
The honest description of what happened: a clean three-effect run on this data does not converge in a single pass. The stored result that circulated with the original materials showed a tconv.max around 0.10: a converged fit, but reaching that took undocumented restarts with prevAns, chaining one run’s ending values into the next run’s starting values, repeated until the ratio came down. The single-pass number you see here is the truthful one.
This is not a defect in the data or the method. It is what a hard SAOM looks like, and the skill this document is really teaching is the discipline to look at tconv.max before you believe the coefficients, and to keep restarting until it passes, rather than reporting the first fit that ran without erroring.
4. What you would actually do next
Not report the table above. Instead, the standard repair loop from the taught session, made concrete:
# Feed the previous run's estimates in as starting values, and refit.
# Repeat until tconv.max < 0.25. Each restart is faster than the first
# because you begin much closer to the solution.
alg <- sienaAlgorithmCreate(projname = NULL, seed = 6886, n3 = 1000)
ans2 <- siena07(alg, data = dat, effects = eff, prevAns = ans,
batch = TRUE, verbose = FALSE, silent = TRUE)
ans2$tconv.max # check; if still > 0.25, do it again with prevAns = ans2Three things are worth saying about this loop, because they are the practical knowledge the original materials hid:
prevAnsis the whole trick. Each restart begins from the last run’s estimates, so it starts near the solution and both converges faster and pushestconv.maxdown. Re-running from scratch, or giving up, are the two things people do instead, and both are wrong.- It can take several passes on data this hard. That is normal. Budget for it.
- If it will not come down after several restarts, the message is about your model, not your patience. A specification the data cannot support will not converge no matter how many times you restart it: at which point you simplify the model, not the threshold.
NoteDepth: why this network is so much harder than s50
Two reasons, both structural.
Density and direction. s50 is 50 nodes at about 5% density, undirected in effect for the structural terms that matter. Duque is 158 nodes at far higher density, fully directed. The simulation the method of moments relies on has to explore a vastly larger and more constrained space, and the targets it is matching are correspondingly noisier.
Composition change via structural zeros. The 10 coding used for states that enter and leave is one supported method-of-moments representation, not an obsolete error. A sienaCompositionChange() object can use known entry and exit times more directly and is needed when creation or endowment effects are part of the model. See the composition-change document for the tradeoffs. This document preserves the structural-zero approach because it is what the published replication files use, and recognizing that coding is part of reading the literature.
5. Exercise
Before touching prevAns, ask whether this data should be modelled as an evolving network at all. Compute the Jaccard index between consecutive waves: the pre-flight check from the taught session, and compare against the manual’s thresholds (≥ 0.3 good, < 0.2 trouble, < 0.1 quite low).
jac <- function(a, b) {
a <- a == 1; b <- b == 1 # ties only; ignore NA / structural zeros
sum(a & b, na.rm = TRUE) / sum(______, na.rm = TRUE)
}Then answer: if the Jaccard indices are healthy but the model still will not converge, what does that tell you about where the problem is: the data, or the specification?
jac <- function(a, b) {
a <- (a == 1); b <- (b == 1) # TRUE where a tie exists
sum(a & b, na.rm = TRUE) / sum(a | b, na.rm = TRUE)
}
data.frame(
transition = c("wave 1 -> 2", "wave 2 -> 3"),
jaccard = round(c(jac(dipl[, , 1], dipl[, , 2]),
jac(dipl[, , 2], dipl[, , 3])), 3))
#> transition jaccard
#> 1 wave 1 -> 2 0.615
#> 2 wave 2 -> 3 0.679The blank is a | b: the Jaccard index is the size of the intersection (ties present at both waves) over the size of the union (ties present at either). Using the union in the denominator is what makes it a measure of stability rather than of density.
The interpretation. These Jaccard values are high: diplomatic ties are extremely stable, far more so than adolescent friendships. So the applicability check passes: this genuinely is an evolving network with enough overlap between waves for SAOM to have something to model. The turnover is not too fast.
That is the useful diagnostic split. Healthy Jaccard plus a convergence failure points at the specification, not the data. The data are fine to model; the three-effect model just is not converging on them in one pass, and the fix is prevAns restarts (and, if those fail, a simpler or better-specified model), not a different dataset and not a relaxed threshold. If the Jaccard had come back below 0.1 the diagnosis would be the opposite: the network turns over too fast to be treated as evolving, and no amount of restarting would save it.
6. Readings
- Duque (2018), “Recognizing international status: A relational approach,” International Studies Quarterly 62(3):577–592. The study. Assigned for the ERGM session as a referee exercise: read it and ask whether it reports the diagnostics this document is about.
- Ripley, Snijders, Boda & Vörös, Manual for RSiena (current edition). The convergence thresholds and the
prevAnsprocedure are in the estimation section. The syllabus’s 2012 citation is badly stale: use the current one. - Snijders (2017), “Stochastic actor-oriented models for network dynamics,” Annual Review of Statistics and Its Application 4:343–363. The canonical modern overview, and the reading the syllabus is missing.
Versions: expand if your numbers differ from mine
data.frame(package = c("RSiena", "networkdata"),
version = c(as.character(packageVersion("RSiena")),
as.character(packageVersion("networkdata"))))
#> package version
#> 1 RSiena 1.5.0
#> 2 networkdata 0.1
cat("R", as.character(getRversion()), "\n")
#> R 4.3.3