---
title: "Composition change in SAOMs"
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"))
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). We fitted `s50`, which is a tidy dataset: 50 girls, all present at all three waves. Your data will not be like that.

**The problem it solves:** actors who join or leave between waves. Students transfer schools, firms are founded and go bankrupt, states enter and exit the international system, patients are discharged. If you have panel network data of any real duration, you have this problem, and the default thing people do about it: drop everyone who was not there the whole time: throws away data and can bias what is left.

**What you need:** `RSiena`, plus the four `klas12p-*.dat` files shipped alongside this document in `data/`. Nothing to download.

:::

# 1. Why this is not just a missing-data problem

There is a distinction here that is easy to blur and expensive to get wrong.

**Missing data** means the tie exists as a question but you do not know the answer. Pupil 7 was off sick the day the survey ran, so you do not know whether 7 nominated 12 as a friend. The tie variable $y_{7,12}$ is well defined; you just did not observe it. RSiena handles this by imputing during simulation and excluding the cell from the target statistics.

**Composition change** means the tie does not exist as a question at all. Pupil 7 had not yet joined the class. There is no fact of the matter about whether 7 nominated 12, because 7 was not there to nominate anybody. $y_{7,12}$ is not unknown, it is **undefined**.

Treating the second as if it were the first inflates your sample with cells that were never real, and: more insidiously: it tells the model that a bunch of actors spent a wave making the deliberate choice to have zero friends. The rate function will then work hard to explain why those actors were so inactive. They were not inactive. They were absent.

::: {.callout-note collapse="true" title="Depth: why the actor-oriented formulation can do this at all"}
This is one of the genuinely nice consequences of writing the model at the level of actors rather than graphs.

A SAOM says: time runs continuously; every so often an actor gets an opportunity to change one tie; when they do, they pick the option that maximises their objective function plus noise. The model is a description of **what individual actors do**. If an actor is not present during some stretch of time, that is not a hole in the data: it simply means the actor had no opportunities during that stretch, which the model expresses directly by not giving them any.

An ERGM or a TERGM has no comparable move available, because it is defined as a probability distribution **over graphs on a fixed node set**. A graph on 26 nodes and a graph on 22 nodes are elements of different sample spaces, and there is no natural way to write a single ERGM over both. This is why the TERGM companion document has to drop four pupils and take the loss, and it is a real, if narrow, advantage of the actor-oriented framing. See [Temporal ERGMs §6.1](tergm.html#a-composition-problem-and-why-it-is-not-a-nuisance).
:::

# 2. The data

Andrea Knecht's Dutch classroom study, class 12p: 33 pupils, four waves of friendship nominations, with pupils entering and leaving. The files ship with this document.

The coding scheme matters and it is exactly the distinction from §1:

| code | meaning | what to do with it |
|---|---|---|
| `0` | no friendship | leave as 0 |
| `1` | friendship | leave as 1 |
| `9` | **incidentally** missing nominator (absent that day) | → `NA`, genuine missing data |
| `10` | **structurally** missing (not in the class at that wave) | → composition change |

```{r}
#| label: read
suppressMessages(library(RSiena))

rd <- function(f) as.matrix(read.table(file.path("data", f)))
wave_files <- paste0("klas12p-friends-wave", c("A", "B", "C", "D"), ".dat")
raw <- lapply(wave_files, rd)

c(n_actors = nrow(raw[[1]]), n_waves = length(raw))
table(unlist(raw))
```

Nearly 400 cells coded `10`. Those are the ones this document is about.

# 3. Building the composition-change object

## 3.1 First, the two missing codes go to `NA`

Both `9` and `10` become `NA` in the network array. That may look like it defeats the purpose: we just said they are different things. The point is that the *distinction* is carried by a **separate object**, not by the tie values. RSiena needs to be told who was present when, and once it knows that, it can tell which `NA`s are structural and which are incidental.

```{r}
#| label: array
n <- nrow(raw[[1]])
arr <- array(NA_integer_, dim = c(n, n, 4))
for (i in 1:4) {
  m <- raw[[i]]
  m[m %in% c(9, 10)] <- NA
  diag(m) <- NA                      # self-ties are never defined
  arr[, , i] <- m
}
fr <- sienaDependent(arr)
fr
```

## 3.2 Then, who was present when

`klas12p-composition.dat` has one row per actor and two columns: the wave at which they **entered** and the wave at which they **left**, on a continuous scale where wave 1 is time 1 and wave 4 is time 4. A value of `2.5` means the actor arrived halfway between waves 2 and 3.

```{r}
#| label: comp
comp <- read.table(file.path("data", "klas12p-composition.dat"))
names(comp) <- c("enters", "leaves")
table(comp$enters, comp$leaves)
```

Thirty of the 33 pupils were present the whole time (`1` to `4`). Three were not: one **left** at 1.5 (partway through period 1), one **left** at 3.5 (partway through period 3), and one **arrived** at 2.5 (partway through period 2).

::: {.callout-tip collapse="true" title="Stuck? What the fractional times mean"}
The intervals are on the same continuous clock the SAOM uses internally. Because the model runs in continuous time between waves, it can accept "this actor joined partway through period 2" and simply give them opportunities to act from that moment onward.

If you do not know when someone joined, use the wave at which you first observed them. Being approximately right is much better than pretending they were there all along, and much better than deleting them.
:::

`sienaCompositionChange()` wants a **list**, one element per actor, each a numeric vector of the intervals during which that actor was present.

```{r}
#| label: compchange
cc_list <- Map(function(a, b) c(a, b), comp$enters, comp$leaves)

movers <- which(comp$enters != 1 | comp$leaves != 4)
str(cc_list[movers])                 # the three who moved
str(cc_list[1])                      # a typical pupil, present throughout
cc <- sienaCompositionChange(cc_list)
```

## 3.3 Covariates and the data object

```{r}
#| label: covars
dem <- read.table(file.path("data", "klas12p-demographics.dat"))
names(dem) <- c("sex", "age", "ethnicity", "religion")
dem[dem == 0] <- NA                  # 0 is the missing code for all four

sex <- coCovar(dem$sex)

dat <- sienaDataCreate(fr, sex, cc)
dat
```

Read the print-out. The line that matters is **`With composition change.`** That is your confirmation that RSiena took the object and will use it. If it is absent, `cc` did not make it into `sienaDataCreate()` and everything below is wrong.

```{r}
#| label: checkpoint-1
checkpoint(actors = n,
           waves = 4,
           movers = sum(comp$enters != 1 | comp$leaves != 4))
```

# 4. Fitting, and what changes

```{r}
#| label: fit
eff <- getEffects(dat)
eff <- includeEffects(eff, transTrip, cycle3)
eff <- includeEffects(eff, sameX, interaction1 = "sex")

alg <- sienaAlgorithmCreate(projname = NULL, seed = 6886, n3 = 1000)

ans <- fit_or_load("saom_comp",
  siena07(alg, data = dat, effects = eff, batch = TRUE, verbose = FALSE, silent = TRUE))

data.frame(effect = ans$effects$effectName,
           est = round(ans$theta, 3),
           se  = round(sqrt(diag(ans$covtheta)), 3))
```

Convergence first, as always:

```{r}
#| label: conv
round(as.numeric(ans$tconv.max), 4)
max(abs(ans$tstat))
```

`tconv.max` under 0.25 and all individual t-ratios under 0.1: this converged properly and you could report it.

Substantively: strong reciprocity, positive transitive triplets, negative 3-cycles (so closure is hierarchical rather than egalitarian: ties close into transitive triads but not into cycles), and same-sex friendship is clearly favoured. Standard adolescent-friendship findings, which is reassuring on a dataset this small.

::: {.callout-note collapse="true" title="Depth: what the rate parameters are telling you"}
First, a note on why they are visible at all. By default `siena07` uses **conditional estimation**, which conditions on the observed number of changes and does not report the rates as estimated parameters. Composition change makes conditional estimation impossible, so RSiena switches to unconditional estimation and the rate parameters appear in the output. You get them here as a side effect of the thing this document is about. (The Core exercise below turns this into the diagnostic it deserves to be.)

The three `constant fr rate` parameters are the estimated number of opportunities each actor gets to change a tie in each period. They are around 6, meaning the model thinks a typical pupil considered changing about six ties between consecutive waves.

Two things to notice.

First, these are **opportunities, not changes.** An actor who gets an opportunity is free to decide their current network is fine and change nothing: that is the "no change" option in the multinomial choice, and it is chosen often.

Second, this is precisely the parameter that composition change protects. Absent actors get **zero** opportunities in the periods when they are away, rather than being scored as actors who had six opportunities and used all of them to do nothing. Drop the composition-change object and the rate estimates absorb the absences, which distorts your reading of how fast the network is actually moving.
:::

# 5. Structural zeros, the alternative approach

You will see a different technique in the literature and in older teaching material: code the absent actors' ties as **structural zeros** (value `10` in the network array, with `allowOnly = FALSE`), which tells RSiena those cells are fixed and not to be modelled.

This is a supported approach for method-of-moments estimation, not a deprecated feature kept only for historical replication. The current RSiena manual presents structural-zero coding and composition-change directives as two ways to represent changing composition. It prefers a composition-change object when endowment or creation effects are needed, and it documents additional caveats for maximum-likelihood and Bayesian estimation. Those qualifications are more specific than a blanket rule against structural zeros.

Two reasons to prefer `sienaCompositionChange()`:

1. **A structural zero says "this tie is fixed at zero."** That is a statement about a tie that exists as a question. Composition change says "this actor was not here." Only the second one is true, and only the second one correctly removes the actor from the rate function.
2. **Composition directives can encode fractional entry and exit times directly.** Structural zero coding records the constraint at observation waves. The directive therefore uses more timing information when those entry and exit times are known.

Structural-zero coding persists because it remains useful and because a large body of published work uses it. If you are replicating a paper that used structural zeros, preserve the coding and say so. For a new method-of-moments analysis with known entry and exit times, a composition-change object is usually the clearer representation, especially when creation or endowment effects are part of the specification.

::: {.callout-warning title="What structural zeros are legitimately for"}
Do not conclude that structural zeros are obsolete generally. They are the right tool whenever a tie is genuinely **impossible** rather than merely absent: a directed network where some actors cannot by construction send ties, two departments that are forbidden to collaborate, a bipartite structure encoded in a one-mode array.

The objection here is narrow and specific: using them to represent *actors who were not there* is a workaround for a problem that now has a purpose-built solution.
:::

# 6. Exercises

::: {.panel-tabset}

## Core

Refit the §4 model **without** the composition-change object: that is, using only `sienaDataCreate(fr, sex)`, and compare the rate parameters with the ones you got in §4.

```{r}
#| eval: false
dat_nocc <- sienaDataCreate(fr, sex)          # note: no `cc`
eff2 <- getEffects(dat_nocc)
eff2 <- includeEffects(eff2, ________, ________)
eff2 <- includeEffects(eff2, sameX, interaction1 = "sex")
ans2 <- siena07(alg, data = dat_nocc, effects = eff2,
                batch = TRUE, verbose = FALSE, silent = TRUE)
```

In two sentences: which parameters move most, and does the direction of the change match what §4's Depth box predicted?

## Stretch

Only three of 33 pupils move in and out here, which is why the §4 vs Core comparison is undramatic. Construct a harsher test: build a second composition-change object that additionally treats the pupils with the **most incidentally missing data** (code `9`) as if they had been absent for the wave in which they are missing, refit, and see how much further the estimates shift.

Then argue: this is the real exercise: whether that second object is a *better* description of the data than the first, or a worse one. There is a defensible answer in each direction and what matters is that you can say which distinction from §1 you are relying on.

## Solution

```{r}
#| label: ex-solution
dat_nocc <- sienaDataCreate(fr, sex)
eff2 <- getEffects(dat_nocc)
eff2 <- includeEffects(eff2, transTrip, cycle3)
eff2 <- includeEffects(eff2, sameX, interaction1 = "sex")

ans2 <- fit_or_load("saom_nocomp",
  siena07(alg, data = dat_nocc, effects = eff2,
          batch = TRUE, verbose = FALSE, silent = TRUE))

c(with_comp = length(ans$theta), without_comp = length(ans2$theta))
```

The two fits do not have the same number of parameters, so you cannot line them up in a data frame. **That is the first finding, not an obstacle.** Here is why:

```{r}
#| label: ex-cond
c(with_comp = ans$cconditional, without_comp = ans2$cconditional)
```

Without composition change, RSiena used **conditional estimation**: it conditions on the observed number of tie changes in each period and does not estimate the rate parameters as free parameters at all: they are recovered afterwards and stored in `$rate`. With composition change, conditional estimation is **not available**, so RSiena silently switched to unconditional estimation and the rates became ordinary parameters in `$theta`.

Nothing warned you about this. It is a good habit to check `length(theta)` and `$cconditional` whenever two SAOM fits refuse to line up.

Now the actual comparison. Rates first, pulling them from wherever each fit put them:

```{r}
#| label: ex-rates
rate_with    <- ans$theta[grep("rate", ans$effects$effectName)]
rate_without <- ans2$rate

data.frame(period     = 1:3,
           with_comp  = round(rate_with, 3),
           without    = round(rate_without, 3),
           pct_diff   = round(100 * (rate_with - rate_without) / rate_without, 1))
```

And the structural effects, matched by name:

```{r}
#| label: ex-struct
nm <- c("outdegree (density)", "reciprocity", "transitive triplets", "3-cycles", "same sex")
data.frame(
  effect     = nm,
  with_comp  = round(ans$theta[match(nm, ans$effects$effectName)], 3),
  without    = round(ans2$theta[match(nm, ans2$effects$effectName)], 3),
  row.names  = NULL)
```

```{r}
#| label: ex-conv
c(with_comp = as.numeric(ans$tconv.max), without_comp = as.numeric(ans2$tconv.max))
```

**Core.** The **rate parameters move and the structural effects barely do**, which is what §4's Depth box predicted, and the pattern of movement is sharper than you might expect.

Look at *which* period moves. Periods 1 and 3 agree to within a fraction of a percent. **Period 2 is off by more than twenty percent**, and period 2 is the one containing the arrival, the pupil who joins at time 2.5.

The general mechanism is the one described earlier: without the composition-change object, an actor who was not in the class is scored as an actor who *was* there, was handed a full allocation of opportunities, and used every one of them to change nothing. The model has only one dial for that, the rate, so it concludes the network was moving more slowly than it really was. That is the direction we see.

I will not pretend to a complete account of why the *arrival* distorts period 2 so much more than the two departures distort periods 1 and 3. Plausibly an arriving actor is the harder case: their ties all become visible at once, and how the estimator reads that burst depends on whether it knows the actor is new. If you want to settle it, the experiment is right there: vary the composition list one mover at a time and refit. That is a genuinely good afternoon's work and I would rather flag the open question than paper over it.

The structural parameters: reciprocity, transitivity, homophily: hardly shift, because they are estimated from the actors who actually acted, and those are the same actors in both fits.

The honest summary of *this* comparison is that with three movers out of 33 the substantive conclusions do not change. **Do not generalise from that.** The distortion scales with the fraction of your panel that turns over, and it lands hardest on the periods where the turnover happens. With a third of the panel moving in and out: completely ordinary in organisational or clinical data: you get rate estimates that are badly wrong in specific periods, and if the timing of the turnover correlates with anything substantive, the structural estimates follow.

So the reason to use composition change here is not that it rescues this analysis. It is that it is the correct model of what happened, it costs three lines, and you will not remember to add it later when the stakes are higher.

**Stretch.** The argument against the harsher object: a pupil who was absent on survey day is a **genuine missing-data** problem, and RSiena's missing-data machinery is designed for it: it imputes during simulation and excludes those cells from the target statistics, which is the statistically correct treatment. Recoding them as composition change tells the model something false (they were not in the class) in order to avoid a problem the model already handles well.

The argument for it: if a pupil is missing at a wave, you have no information about their choices at that wave either way, and the rate function is arguably better off not attributing opportunities to them. This is more defensible when the missingness is heavy and long: a pupil missing three of four waves is functionally absent whatever the coding says.

The distinction from §1 that decides it: **was the tie well defined and unobserved, or undefined?** For a pupil who was enrolled but off sick, the tie was well defined. That points to missing data, not composition change. The line genuinely blurs when absences get long enough that "enrolled" stops describing anything real, and at that point you should say in your write-up which choice you made and why, because a reader cannot tell from the coefficients.
:::

# 7. Readings

- **Ripley, Snijders, Boda & Vörös**, *Manual for RSiena* (current edition: check [stats.ox.ac.uk/~snijders/siena](https://www.stats.ox.ac.uk/~snijders/siena/) for the latest). The composition-change section is short and worth reading in full. **Note the syllabus cites a 2012 edition with Preciado as an author**, which is roughly fourteen years stale, and the author line has changed. Use the current one.
- **Huisman & Snijders (2003)**, "Statistical analysis of longitudinal network data with changing composition," *Sociological Methods & Research* 32(2):253–287. The method underneath `sienaCompositionChange()`.
- **Huisman & Steglich (2008)**, "Treatment of non-response in longitudinal network studies," *Social Networks* 30(4):297–308. On the §1 distinction, done carefully, with simulations showing when it matters.

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