---
title: "A hard SAOM: diplomatic recognition (Duque 2018)"
engine: knitr
editor: source
subtitle: "Ancillary self-study · ICPSR Network Analysis: Advanced Topics"
---

```{r}
#| label: setup
#| include: false
source("_helpers.R")
require_pkgs(c("RSiena", "networkdata"))
set.seed(6886)
```

{{< include _ancillary_header.qmd >}}

::: {.callout-note title="What 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:

1. **The node set changes**: states enter and leave the international system across waves.
2. **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.

```{r}
#| label: prep
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)))
table(dipl, useNA = "ifany")
```

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 `NA`s (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:

```{r}
#| label: prep-check
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)
```

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.

::: {.callout-warning title="The 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:

1. **Positional indexing into a name-dimensioned array.** It looped `for (i in 1:nrow(d))` and wrote `dipl[i, j, t] <- val`, but `dipl` was dimensioned by the *pooled* 158-actor list while `d` had 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.

2. **Silent symmetrisation.** It wrote both `dipl[i, j, t] <- val` and `dipl[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.

3. **Coercion of the missing code without noticing.** `as.numeric(as.character(d))` turned the `"."` missing code into `NA` with a warning nobody read: 586 cells. We handle `"."` explicitly with `d[d == "."] <- NA` so 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.

```{r}
#| label: fit-setup
dv  <- sienaDependent(dipl, allowOnly = FALSE)
dat <- sienaDataCreate(dv)
eff <- getEffects(dat)
eff <- includeEffects(eff, transTies)
```

```{r}
#| label: fit-show
#| eval: false
# 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)
```

```{r}
#| label: fit-load
#| include: false
ans <- fit_or_load("saom_duque", {
  alg <- sienaAlgorithmCreate(projname = NULL, seed = 6886, n3 = 1000)
  siena07(alg, data = dat, effects = eff, batch = TRUE, verbose = FALSE, silent = TRUE)
})
```

The result below is precomputed from the code above.

```{r}
#| label: fit-result
data.frame(effect = ans$effects$effectName,
           est = round(ans$theta, 3),
           se  = round(sqrt(diag(ans$covtheta)), 3))
```

Substantively 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.**

```{r}
#| label: fit-conv
round(as.numeric(ans$tconv.max), 3)
```

::: {.callout-important title="This 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:

```{r}
#| label: repair-show
#| eval: false
# 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 = ans2
```

Three things are worth saying about this loop, because they are the practical knowledge the original materials hid:

1. **`prevAns` is the whole trick.** Each restart begins from the last run's estimates, so it starts near the solution and both converges faster and pushes `tconv.max` down. Re-running from scratch, or giving up, are the two things people do instead, and both are wrong.
2. **It can take several passes** on data this hard. That is normal. Budget for it.
3. **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.

::: {.callout-note collapse="true" title="Depth: 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](saom_composition.html) 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

::: {.panel-tabset}

## Core

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).

```{r}
#| eval: false
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?

## Solution

```{r}
#| label: ex-jaccard
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))
```

**The 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 `prevAns` procedure 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.

```{r}
#| label: versions
#| cache: false
#| code-fold: true
#| code-summary: "Versions: expand if your numbers differ from mine"
data.frame(package = c("RSiena", "networkdata"),
           version = c(as.character(packageVersion("RSiena")),
                       as.character(packageVersion("networkdata"))))
cat("R", as.character(getRversion()), "\n")
```
