---
title: "Day 12: Stochastic Actor-Oriented Models"
subtitle: "ICPSR Network Analysis: Advanced Topics · Day 12"
author: "Shahryar Minhas"
date: today
engine: knitr
editor: source
format:
  html:
    theme: [cosmo, _saom_session.scss]
    toc: true
    toc-depth: 3
    toc-location: left
    toc-title: "Run sheet"
    embed-resources: true
    code-tools: true
    code-link: false
    fig-width: 7
    fig-height: 4
    df-print: default
knitr:
  opts_chunk:
    message: false
    warning: false
    comment: "#>"
    collapse: true
    dpi: 140
execute:
  cache: false
  error: false
---

```{r}
#| label: setup
#| include: false
source("_saom_helpers.R")
# Isolate this session's fitted objects from any other document in this folder
formals(fit_or_load)$dir <- "cache_saom"
set.seed(6886)
options(width = 90, digits = 4)
```

## Before we start {#preflight}

This is a walkthrough, not a deck. The prose between the chunks is the lecture. If you are remote, or if you are reading this a week later on your own, you should be able to run every chunk in order and get the same numbers I got; every fit in this document is seeded.

[Open the Day 12 teaching deck](day12_saom_teaching_deck.html){target="_blank"}.

::: {.callout-note title="The three helper functions this document depends on"}
This document uses three small helpers defined in **`_saom_helpers.R`**: `require_pkgs()`, `fit_or_load()`, and `checkpoint()`. The helper file must sit in your working directory alongside the **`cache_saom/`** folder of pre-fitted objects. (The file also defines a fourth, `quiet_solver()`, which belongs to the ERGM day and is never called here.) **Source it first, before anything else runs:**

```r
source("_saom_helpers.R")
```

`require_pkgs()` fails early and loudly if a package is missing; `fit_or_load()` returns a shipped fit instantly (or refits from scratch when you set `options(session.refit = TRUE)`); `checkpoint()` prints the compact comparison lines you'll see throughout. The very first visible line below, `require_pkgs(c("RSiena", "netify"))`, is one of them, if you skipped the `source()` it will die with "could not find function."
:::

Today's path:

| Block | What |
|---|---|
| **A** | Is this even a SAOM problem? |
| **B** | Fit it, converge it, read it |
| **C** | Co-evolution, and what it does not buy you |
| Final | Wrap |

::: {.callout-note collapse="true" title="Already comfortable with parts of this? Where to jump"}
* You've fit SAOMs before and want the diagnostics and limitations → [Block C](#block-c)
* You want the convergence machinery specifically → [B2](#b2)
* You want to know how to *referee* a SAOM paper without fitting one → [the referee checklist](#referee)
* You want the mathematics the core path skips → [Appendix A](#appendix-a)
:::

### Setup check: Run this the night before {#setup-check}

If this chunk does not print `SETUP OK`, post in the shared doc *before* the session rather than during the session.

```{r}
#| label: setup-check
require_pkgs(c("RSiena", "netify"))
suppressMessages({library(RSiena)})
stopifnot(exists("s501"), exists("s50a"), exists("s50s"))
cat("RSiena  ", as.character(packageVersion("RSiena")), "\n")
cat("R       ", as.character(getRversion()), "\n")
cat("SETUP OK\n")
```

::: {.callout-important title="A version note you need before you read anything else"}
This document is written against **RSiena 1.5.0**, which is what is installed here. CRAN is further along, and **RSiena ≥ 1.6.3 requires R ≥ 4.5.0**, so on an older R, `install.packages("RSiena")` will silently offer you nothing rather than telling you why. That failure mode is easy to encounter and poorly explained, so it is worth knowing about before you hit it.

More importantly: **RSiena 1.6 renamed essentially the whole user-facing API.** The old names still work, nothing you learn today stops working, but the current manual uses the new ones exclusively. The full translation table is in a collapsed box in [A1](#a1), reference on demand, not something to memorise.
:::

---

# Block A: Is this even a SAOM problem? {#block-a}

## A1 · The motivating question {#a1}

Start with a question the earlier material left open: *suppose I want to study how smoking and friendships develop together. Can I answer that by comparing the network at time 1 with the network at time 2?*

No, not with a simple before-and-after comparison. Yesterday, with cross-sectional ERGMs, we could write down a transitivity parameter and test it, but from a single cross-section we could not distinguish *triadic closure* (I befriend my friends' friends) from *homophily on something we did not measure*. That is the problem raised by Shalizi and Thomas, which you met on Day 10. One snapshot is not enough information.

Today's answer is: **use time**. Yesterday's C4 already made that move once, tie-oriented, the TERGM conditioned each wave on the last, with no actor anywhere in the machinery. Today the actor takes the wheel, and whether that is an upgrade or just a different bet on the same waves is [C3](#c3)'s fight. If you observe the same people at several waves, you can ask whether ties changed to match behavior (selection) or behavior changed to match ties (influence). The temporal ordering helps, but the separation still depends on the model's continuous-time micro-step assumptions.

That is the central promise of the stochastic actor-oriented model, and by the end of Block C you should also understand its limits. It provides a useful separation of mechanisms under a declared process model, but it does not deliver as much as some abstracts claim.

### Did Similarity Produce the Relationship, or Did the Relationship Produce Similarity? {#sec-political-question}

The same question appears in many settings. Do states form alliances because their security policies already align, or do alliances bring those policies closer together? Do legislators collaborate because they already hold similar positions, or does collaboration bring their positions closer? Do armed groups cooperate because their strategies already align, or does cooperation lead their strategies to converge? Repeated observations help us see what changed first, but they do not record every decision or change between waves. A SAOM fills in that missing sequence using an assumed actor-oriented process, so its results show patterns consistent with selection or influence rather than proving either process.

A SAOM is plausible only when the actor-oriented story makes substantive sense. Actors must have meaningful opportunities to change an outgoing tie or behavior, and the observed waves must be snapshots of a process that plausibly continued between observations. The model does not require actors to solve a literal optimization problem consciously, but its stochastic choice process still needs to be a defensible approximation to how change occurred.

Before fitting, ask:

1. Who has agency over the outgoing relation?
2. What changes were possible between waves?
3. Were actors present and eligible throughout each period?
4. Which network mechanisms could produce change, such as reciprocity, closure, or homophily?
5. Which behavior process could produce change, and how could the network enter it?
6. What unmeasured process could still generate both selection and influence patterns?

### The Applied Model Compass {#sec-model-compass}

| Model or specification | Question we want to understand | How it helps | What makes it distinct | What we should be able to say after fitting |
|---|---|---|---|---|
| Network-evolution SAOM | Which local network changes are more attractive when actors receive opportunities to alter outgoing ties? | Represents the unobserved sequence between panel waves as stochastic actor-level micro-steps | Places agency and the opportunity to change at the actor level | "When an adolescent has an opportunity to change a friendship, reciprocating an existing nomination makes the resulting network state more attractive under the fitted objective function." |
| Network-behavior co-evolution SAOM | Are similarity patterns more consistent with selecting similar peers, becoming similar to peers, or both? | Estimates a network equation and a behavior equation within one simulated change process | Separately parameterizes selection and influence while allowing them to evolve together | "Students are more likely to choose peers with similar drinking behavior, while the estimated peer-influence component is weaker and more uncertain." |
| TERGM comparator | Is the theory better expressed as a distribution over network transitions rather than actor micro-steps? | Conditions each network on its observed past without requiring actor-oriented opportunities | Makes the transition, rather than the actor's choice occasion, the modeling unit | "The lagged network predicts the next friendship configuration, but the model does not claim that one actor made each intermediate change." |

The interpretation should always name the candidate change being compared. A positive SAOM evaluation parameter is not an unconditional probability that a tie forms, and a fitted selection or influence term is not automatically a causal effect.

### One note on names, then we move on {#names}

We use the installed **1.5** function names throughout today, they're what's installed, they still work in 1.6, and every published replication script and every tutorial online uses them. The current manual uses the renamed 1.6 API, so the full 1.5→1.6 map lives in the collapsed box below (and in [Appendix A](#appendix-a)); you never have to memorise it. With that parked, the motivating question runs straight into the pre-flight checks and the micro-step.

::: {.callout-note collapse="true" title="Reference on demand: The full 1.5 → 1.6 name-translation table"}
In version 1.6, RSiena renamed most of its functions. The old names still work, nothing you learn today stops working, but the manual you'll download uses the new ones:

| What we type today (1.5) | What the current manual calls it (1.6+) |
|---|---|
| `siena07()` | `siena()` |
| `sienaDependent()` | `as_dependent_rsiena()` |
| `coCovar()`, `varCovar()` | `as_covariate_rsiena()` |
| `sienaDataCreate()` | `make_data_rsiena()` |
| `getEffects()` | `make_specification()` |
| `includeEffects()`, `setEffect()` | `set_effect()` |
| `sienaGOF()` | `test_gof()` |
| `sienaTimeTest()` | `test_time()` |
| `selectionTable()` / `influenceTable()` | `interpret_selection()` / `interpret_influence()` |
| `print01Report()` | `write_report()` |
| Effect arguments `interaction1=` / `interaction2=` | `covar1=` / `covar2=` |

There is also `transformScript()`, which converts a whole file for you.

**Why we're not just teaching the new names:** switching would require an R upgrade in the middle of the workshop, even though the installed names remain supported aliases. That interruption would add setup work without changing the model or interpretation, so we use the installed API and provide the translation table for current documentation.
:::

::: {.callout-note collapse="true" title="Depth: A correction to the syllabus citation"}
The syllabus lists *Ripley, Snijders & Preciado, RSiena manual (2012)*. The current manual is **Ripley, Snijders, Boda & Vörös, "Manual for RSiena"**, Oxford/Groningen, revised 2026, with contributions from Huisman, Schweinberger, Steglich and Amati. Preciado is no longer on the author line and the function names have all changed. Use the current one.

Worth adding to your reading, and not on the syllabus: **Snijders (2017), "Stochastic Actor-Oriented Models for Network Dynamics," *Annual Review of Statistics and Its Application* 4:343–363**, the canonical modern overview, and much shorter than the manual.
:::

## A2 · Pre-flight: Do not fit before you check this {#a2}

Keep the same rule we have used throughout the course: *diagnose before you fit*. Applied to SAOMs, that means checking whether the observed network changes could plausibly arise from the process the model assumes. This is one of the most important checks in the entire workflow.

The data: **s50**, the "teenage friends and lifestyle study", 50 Scottish girls, three waves, with friendship, alcohol and smoking measured at each wave. It ships with RSiena, so there are no files to download and nothing to install. I use it because friendship nominations and behavior changes require little background knowledge, which lets us concentrate on the model's assumptions and interpretation.

```{r}
#| label: a2-data
# These objects are lazy-loaded by RSiena; no data() call is needed, and
# data(s50) will actually warn at you, because the objects are named individually
dim(s501); dim(s502); dim(s503)
table(s50a, useNA = "ifany")   # Alcohol, 1-5
table(s50s, useNA = "ifany")   # Smoking, 1-3
s50_net <- netify::netify(
  list(s501, s502, s503),
  symmetric = FALSE,
  diag_to_NA = TRUE,
  missing_to_zero = FALSE
)
stopifnot(all(unlist(netify::validate_netify(s50_net))))
s50_friend_array <- simplify2array(netify::get_raw(s50_net))
```

Three 50×50 adjacency matrices, and two 50×3 behavior matrices. `netify()` verifies the common actor order, direction, and dimensions before `simplify2array()` creates the 50×50×3 array that RSiena requires. The behavior matrices remain separate because they are actor outcomes rather than networks.

Now the check. SAOMs model a network *evolving*, the model assumes wave 2 is reachable from wave 1 by a sequence of small steps. If the network turns over almost completely between waves, that assumption is false and the model is not applicable, no matter how nicely it converges.

The diagnostic is the **Jaccard index** of successive waves: of the ties present in either wave, what fraction is present in both.

```{r}
#| label: a2-jaccard
jaccard <- function(a, b) {
  keep <- !is.na(a) & !is.na(b)
  n11 <- sum(a[keep] == 1 & b[keep] == 1)
  n10 <- sum(a[keep] == 1 & b[keep] == 0)
  n01 <- sum(a[keep] == 0 & b[keep] == 1)
  n11 / (n11 + n10 + n01)
}
j12 <- jaccard(s501, s502)
j23 <- jaccard(s502, s503)
deg <- c(mean(rowSums(s501)), mean(rowSums(s502)), mean(rowSums(s503)))

checkpoint(`jaccard 1->2` = j12, `jaccard 2->3` = j23, `avg degree` = deg)
```

The manual's thresholds, quoted:

> *"Jaccard values of .3 and higher are good; values lower than .2 indicate that there might be difficulties in estimation; values lower than .1 are quite low indeed."*

We're at `r round(j12,3)` and `r round(j23,3)`. Comfortably fine.

The average degree matters too, because low Jaccard is *acceptable* if the network is growing or shrinking monotonically, since we would expect low overlap then. Ours is roughly flat (`r paste(round(deg,2), collapse=" → ")`), so the Jaccard provides useful information. The manual's sharper version:

> *"If Jaccard indices are very low while the average degree is not strongly increasing, this indicates that the turnover in the network may be too high to consider the data as an evolving network, and perhaps the SIENA method is not suitable for the data set."*

That is four lines of code and it is the first question you should ask of any SAOM paper you referee.

::: {.callout-note collapse="true" title="Depth: What else to check before fitting"}
Beyond Jaccard: **number of waves** (two waves gets you tie dynamics but influence is essentially unestimable, you need three as a working minimum, four before it behaves); **actor set stability** (SAOM assumes a roughly fixed actor set; joiners and leavers need `sienaCompositionChange()`); and **missingness pattern**, RSiena handles missing ties by treating them as non-informative for the simulation, which is fine when missingness is sparse and structural, and quietly harmful when it isn't.

`s50` has no missing data at all, which is why it's the teaching dataset rather than the realistic one.

::: {.callout-note collapse="true" title="More Detail: Structural Zeros Versus Actors Entering or Leaving"}
Two different data situations are easy to confuse.

**A structural zero concerns a particular tie.** In an RSiena network array, `0` means a tie is absent but possible, while `10` means the tie is structurally impossible. For a directed network, `friend_array[1, 2, ] <- 10` prohibits only $1\rightarrow2$. Prohibiting both directions requires setting both `[1, 2, ]` and `[2, 1, ]` to 10. RSiena removes a structurally fixed tie from the meaningful choices rather than treating it as an ordinary tie that the actor declined to make. The value `11` represents a structurally required tie.

Use this coding only when the dyad itself cannot take another value, such as a tie prohibited by the study design or institution. The distinction is:

```r
friend_array[1, 2, 1] <- 0   # absent at wave 1, but possible
friend_array[1, 2, 1] <- 10  # impossible at wave 1
```

**Composition change concerns whether an actor exists in the modeled population.** If a student joins or leaves, coding every associated dyad as a structural zero is a crude workaround. It can wrongly imply that the student was unavailable for an entire interval even when entry or exit occurred between waves. The preferred approach is to keep a common superset containing everyone observed at any wave, preserve the same actor order in every adjacency matrix, and tell RSiena when each actor was present.

```r
n_actors <- dim(friend_array)[1]
n_waves  <- dim(friend_array)[3]
composition_list <- rep(list(c(1, n_waves)), n_actors)

# Actor 12 joins at wave 2.
composition_list[[12]] <- c(2, n_waves)

# Actor 27 leaves after wave 2.
composition_list[[27]] <- c(1, 2)

# Actor 34 enters halfway between waves 1 and 2.
composition_list[[34]] <- c(1.5, n_waves)

composition <- sienaCompositionChange(
  composition_list,
  nodeSet = "Actors",
  option = 3
)
```

Each pair gives the beginning and end of an actor's presence. The times can fall between integer-valued observation waves. Someone who leaves and returns can have more than one pair, such as `c(1, 1.8, 2.4, 3)`. With `option = 3`, values are treated as missing whenever the actor is absent, with previous values used internally when available but still treated as missing. Options 1 and 2 process pre-entry and post-exit values differently, so choose among them based on what is observed and what absence means in the application.

While an actor is absent, a clear data setup is to mark their incoming and outgoing network observations as `NA` and let the composition object identify those cells as absence from the population:

```r
friend_array[12, , 1] <- NA
friend_array[, 12, 1] <- NA

actors <- sienaNodeSet(n = n_actors, nodeSetName = "Actors")
friend <- sienaDependent(friend_array, nodeSet = "Actors")
dat <- sienaDataCreate(friend, composition, nodeSets = list(actors))
```

The node-set name on the composition object must match the dependent network. Explicit `nodeSets` can be unnecessary for a simple one-mode network, but when supplied it is a list, as shown here, and makes that match visible.

This matters even more in a co-evolution model. Network structural zeros stop the corresponding ties from changing, but they do not automatically stop an absent actor's behavior from changing. A composition-change object makes the actor active only during the stated intervals and handles participation more coherently across the network and behavior processes.

Statistically, adjusting for changing composition does not mean adding a joiner or leaver control variable. It means changing the process's risk set. Before entry, the actor cannot receive opportunities or participate in feasible tie changes. While present, the actor participates normally. After exit, the actor stops participating. This prevents absence from being interpreted as mass tie dissolution, social isolation, avoidance, or ordinary missing data.

The practical rule is simple: use `10` when a **specific dyad is structurally impossible**; use `sienaCompositionChange()` when an **actor enters or leaves the population**.
:::
:::

::: {.callout-note collapse="true" title="The full toolbox: `print01Report()`, the canonical pre-fit artifact (optional self-study)"}
We hand-computed the Jaccard above because writing the four lines yourself is how you learn what the number is. RSiena will also compute it *for* you, along with everything else you should look at before fitting, via `print01Report()`, the same "initial report" that `siena07` writes when you give it a `projname`. It needs a siena data object, which we don't build until [B1](#b1), so this aside builds a minimal one inline:

```{r}
#| label: tb-a-report
frnd0 <- sienaDependent(s50_friend_array)
dat0  <- sienaDataCreate(frnd0)
print01Report(dat0, modelname = "s50_prefit_report")
rl <- readLines("s50_prefit_report.txt")
writeLines(rl[grep("Tie changes", rl) + (0:3)])
```

That writes **`s50_prefit_report.txt`** into your working directory, `r length(rl)` lines of densities, degrees, missingness, and the change table shown above. Its `Jaccard` column is exactly the `r round(j12, 3)` and `r round(j23, 3)` we computed by hand, which is the point: the hand version teaches you the diagnostic, the report is what you run on your own data, and when the two agree you know you understood it. Read the whole file once in your life; skim it for every dataset you ever fit.
:::

## A3 · The micro-step {#a3}

This is the main new estimation idea, and it is worth slowing down for.

Everything you've fit so far models a network **state**. An ERGM puts a probability on the whole graph. SAOM does something different: it models the **transition** from one observed wave to the next, as the accumulation of many tiny unobserved changes.

The story runs like this. Between wave 1 and wave 2 time runs continuously. At random moments, one randomly chosen actor gets an opportunity to change **one** thing, a single outgoing tie, or their own behavior score, never both, never two at once. They choose among the options, favouring the ones they prefer (probabilistically, piece 3 makes this precise). Repeat until you arrive at wave 2.

Three pieces:

**1. The rate function**, how often does actor $i$ get an opportunity? Waiting times are exponential; the rate parameter $\lambda$ is estimated from how much change there is between waves. More change between waves ⇒ higher estimated rate. This is why the rate parameters are usually not substantively interesting: they're absorbing the wave spacing.

**2. The objective function**, when actor $i$ *does* get an opportunity, what do they want? This is a linear combination of network statistics evaluated from $i$'s point of view:

$$f_i(x) \;=\; \sum_k \beta_k \, s_{ik}(x)$$

Where $s_{ik}(x)$ might be $i$'s outdegree, or the number of reciprocated ties $i$ has, or the number of transitive triplets $i$ is embedded in. The $\beta_k$ are what we estimate.

**3. The choice**, among all the networks reachable in one step (change one tie, or change nothing), actor $i$ picks probabilistically, with

$$\Pr(i \text{ moves to } x') \;=\; \frac{\exp\{f_i(x')\}}{\sum_{x''} \exp\{f_i(x'')\}}$$

**That is a multinomial logit.** The choice set is "which one tie do I change, or do I make no change?" If you understand a discrete-choice model, you already understand the core of the SAOM. The remaining work is defining which network and actor statistics enter $s_{ik}$ and simulating the unobserved sequence of choices between waves.

::: {.callout-important title="The one interpretive rule that matters"}
A SAOM coefficient is a contribution to the **score assigned to a candidate state at a modeled micro-step**. It is *not* a probability, *not* a marginal effect on the observed network, and *not* comparable to an ERGM coefficient even when the term has the same name. The micro-step is an unobserved part of the model, so phrases such as "students chose reciprocity" overstate what the panel records. We return to the comparison in [B3](#b3) with an actual number.
:::

### How SAOM Estimation Works {#sec-estimation-roadmap}

Default SAOM estimation uses **simulation-based method of moments**. In one sentence, RSiena chooses parameter values that make important summaries from simulated network evolution match the corresponding summaries in the observed panel. Open the box below for the complete logic, including a numerical example, Robbins-Monro updating, convergence, and the difference from maximum likelihood.

::: {.callout-note collapse="true" title="More Detail: What Method of Moments Is Doing"}
#### Start with the basic idea

**Method of moments estimation** means choosing parameter values so that important features produced by the model match the corresponding features observed in the data.

A **moment** here does not mean a moment in time. It means a numerical summary, such as:

- The number of ties
- The amount of reciprocity
- The number of transitive triplets
- The tendency for similar actors to form ties
- The amount of network or behavior change between waves

#### Begin with a non-network example

Suppose observations come from a distribution with an unknown mean $\mu$, and the sample mean is 10. The model implies

$$
E_{\mu}[Y]=\mu.
$$

Method of moments sets the model's expected mean equal to the observed mean:

$$
E_{\mu}[Y]-\bar Y=0.
$$

Because $\bar Y=10$, the solution is $\hat\mu=10$. We choose the parameter value that makes the model reproduce the observed summary.

#### Apply the same idea to a SAOM

For a panel network, we observe the network at two waves,

$$
X(t_1) \qquad \text{and} \qquad X(t_2),
$$

but we do not observe every individual tie decision between them. RSiena starts with a parameter vector such as

$$
\theta=
\begin{bmatrix}
\beta_{\text{outdegree}}\\
\beta_{\text{reciprocity}}\\
\beta_{\text{transitivity}}\\
\beta_{\text{similarity}}
\end{bmatrix},
$$

then simulates many possible sequences of actor-level micro-steps between the observed waves. From each simulated path it calculates the same target statistics constructed from the observed panel, including the amount of change, the number of ties, reciprocity, transitive triplets, and similarity-related configurations.

Let $S_{\mathrm{obs}}$ collect the observed target statistics and let $S_{\mathrm{sim}}(\theta)$ collect their simulated counterparts. The estimator seeks parameter values satisfying

$$
E_{\theta}\!\left[S_{\mathrm{sim}}(\theta)\right]-S_{\mathrm{obs}}=0.
$$

This equation says: find parameter values under which the average simulated statistics equal the observed statistics.

#### Follow one illustrative reciprocity update

Suppose, only for illustration, that the observed reciprocity target is 40 but paths simulated under the current parameters average 25:

$$
S_{\mathrm{sim,reciprocity}}-S_{\mathrm{obs,reciprocity}}=25-40=-15.
$$

The current simulations are producing too little reciprocity. Using the joint relationship between all parameters and all target statistics, the algorithm adjusts the parameters in a direction that should reduce that discrepancy. If the next simulations average 35 reciprocal configurations, the model is closer but not finished. It continues until the simulated reciprocity target is approximately 40 while also bringing the other simulated targets close to their observed values.

The parameters are estimated jointly. Changing the reciprocity parameter can also change tie count, transitivity, and other statistics, so RSiena cannot repair one discrepancy at a time while ignoring the rest.

#### See why changing a coefficient changes the simulations

The evaluation function assigns each candidate next state a score:

$$
f_i(x)=\sum_k\beta_k s_{ik}(x).
$$

Suppose a proposed friendship increases outdegree by 1, creates reciprocity, and completes two transitive triplets. With

$$
\beta_{\text{outdegree}}=-2,\qquad
\beta_{\text{reciprocity}}=1.5,\qquad
\beta_{\text{transitivity}}=0.4,
$$

the proposed friendship changes the evaluation score by

$$
\Delta f_i=-2(1)+1.5(1)+0.4(2)=0.3.
$$

That score helps determine the action's probability relative to the actor's other available micro-steps. Changing the coefficients changes those relative choice probabilities. Different probabilities produce different simulated paths, which produce different target statistics.

#### Understand what Robbins-Monro contributes

The expected simulated statistics cannot usually be calculated exactly, so RSiena estimates them with simulation. This makes the discrepancy $S_{\mathrm{sim}}(\theta)-S_{\mathrm{obs}}$ noisy. One simulated set of paths can produce too much reciprocity even when the parameter is close to its final value.

Robbins-Monro stochastic approximation handles that noise by updating the parameters gradually:

$$
\theta_{r+1}
=
\theta_r
-
a_rD_r^{-1}
\left[
S_{\mathrm{sim}}(\theta_r)-S_{\mathrm{obs}}
\right].
$$

Here $D_r$ describes how the expected statistics change when the parameters change, and $a_r$ controls the update size. At a high level, RSiena:

1. Starts with parameter guesses
2. Simulates network evolution
3. Calculates the simulated target statistics
4. Compares them with the observed targets
5. Adjusts the parameters in a direction that reduces the joint mismatch
6. Repeats with progressively smaller adjustments

The shrinking adjustment size helps keep random simulation variation from making the estimates bounce indefinitely.

RSiena organizes this work into three phases. Phase 1 learns the local relationship between parameters and target statistics. Phase 2 carries out the repeated parameter updates with decreasing step sizes. Phase 3 holds the estimates fixed and uses additional simulations to calculate standard errors and convergence diagnostics.

#### Read convergence correctly

Good convergence means that, at the estimated parameters, the average simulated target statistics are sufficiently close to the observed target statistics. An effect-specific convergence t-ratio is roughly

$$
\frac{
E_{\hat\theta}[S_{\mathrm{sim}}]-S_{\mathrm{obs}}
}{
\text{Monte Carlo standard deviation}
}.
$$

A value near zero means the estimator is matching that targeted statistic relative to the remaining simulation noise.

Convergence does **not** mean:

- The model is substantively true
- Actors literally follow the proposed decision process
- The model reproduces every feature of the network
- A selection or influence estimate is causal
- No important variables or network mechanisms were omitted

It means only that the estimation algorithm has approximately solved the chosen moment equations. This is why convergence must be checked before interpreting coefficients, but convergence alone is not evidence of good overall fit.

#### Distinguish method of moments from maximum likelihood

Maximum likelihood asks which parameter values make the exact observed data most probable:

$$
\hat\theta_{\mathrm{ML}}
=
\arg\max_{\theta}L(\theta\mid X_{\mathrm{obs}}).
$$

Method of moments instead asks which parameter values make selected model-generated summaries match the corresponding observed summaries:

$$
E_{\theta}[S(X)]=S(X_{\mathrm{obs}}).
$$

Default RSiena estimation does not maximize the full likelihood of the observed panel. It simulates the unobserved paths and searches for parameters that match the selected target statistics. RSiena also offers a separate maximum-likelihood estimator, but that is not what the standard `siena07()` workflow is doing.

#### See why the choice of moments matters

In the standard setup, each parameter is paired with a target statistic:

| Parameter | Main observed target |
|---|---|
| Rate | Amount of network or behavior change between waves |
| Outdegree | Overall number of ties |
| Reciprocity | Reciprocal tie configurations |
| Transitive triplets | Transitive configurations |
| Similarity effect | Ties between actors with similar covariate or behavior values |

These targets are statistically connected, so RSiena solves their moment equations jointly. Estimation forces the model to match the targeted statistics, within simulation error. It does not directly force the model to reproduce the full degree distribution, every triad-census category, geodesic distances, subgroup-specific clustering, or other omitted features.

That is why goodness-of-fit should examine meaningful statistics that were **not directly targeted** during estimation. Otherwise, we would mostly be checking whether the estimator solved the equations it was designed to solve.

The most intuitive summary is:

> Suppose actors behaved according to these parameter values. If we replayed the period between the observed waves many times, would the resulting networks resemble the changes that were actually observed?

If the answer is no, RSiena revises the parameters and simulates again.
:::

::: {.callout-note collapse="true" title="Depth: The manual's own candour about identification"}
The manual states the limitation directly, and this is worth reading carefully:

> *"The Stochastic Actor-Oriented Model is so complicated mathematically that it is hard to give proofs of any properties... Mathematical proofs of desirable properties are so hard that they should not be prerequisites for the use of the model."*

And on identifiability specifically:

> *"The fact that there is convergence strongly suggests that there is identifiability."*

That is an argument from behavior, not from theory. I point it out not to undermine the model but because it changes what the convergence check in [B2](#b2) *means*: it isn't a box you tick before reporting, it's a substantial part of the evidence that the thing you fitted is estimable at all.
:::

```{r}
#| label: a3-checkpoint
checkpoint(`waves` = 3, `actors` = nrow(s501),
           `ties w1/w2/w3` = c(sum(s501), sum(s502), sum(s503)),
           label = "END BLOCK A")
```

---

# Block B: Fit it, converge it, read it {#block-b}

## B1 · The first fitted model {#b1}

We start with network dynamics only, no behavior yet. The pipeline is four objects and then a fit.

```{r}
#| label: b1-build
friend  <- sienaDependent(s50_friend_array)
smoke1  <- coCovar(s50s[, 1])       # Smoking at wave 1, treated as fixed
dat1    <- sienaDataCreate(friend, smoke1)
dat1
```

Read that summary, it tells you the number of observations, the actor set, and the density per wave. If `sienaDataCreate` reports a different number of waves than you expect, stop there; something is wrong with your array. Drinking is intentionally absent from this first model. Friendship is the only changing dependent variable, while smoking is a fixed wave-1 covariate. Drinking enters later, when we build the joint co-evolution model.

Now the effects. `getEffects()` gives you a default specification, outdegree and reciprocity, and you add to it.

```{r}
#| label: b1-effects
eff1 <- getEffects(dat1)
eff1 <- includeEffects(eff1, recip, transTrip, cycle3)
eff1 <- includeEffects(eff1, sameX, interaction1 = "smoke1")
```

`transTrip` is transitive triplets, closure. `cycle3` is three-cycles, which is generalized exchange rather than hierarchy; including both lets closure and non-hierarchical reciprocation compete rather than loading onto one term. `sameX` with `interaction1="smoke1"` is an exact match on the fixed wave-1 smoking category.

### The Complete First Specification {#b1-specification}

This is a network-evolution SAOM, not a co-evolution model. When student \(i\) receives a friendship-change opportunity, the model scores every network that could result from adding one outgoing nomination, deleting one outgoing nomination, or making no change. Schematically, the evaluation function is

$$
f_i^{\text{friend}}(x)
=
\beta_{\text{out}}s_{\text{out},i}(x)
+\beta_{\text{rec}}s_{\text{rec},i}(x)
+\beta_{\text{trans}}s_{\text{trans},i}(x)
+\beta_{\text{cycle}}s_{\text{cycle},i}(x)
+\beta_{\text{smoke}}s_{\text{same smoke},i}(x).
$$

Here \(x\) is a candidate friendship network and each \(s\) records what that candidate network would contain for student \(i\). The same-smoking statistic counts outgoing nominations to students in the same wave-1 smoking category. It does not allow smoking to change, and it does not include a drinking effect.

| Part of the model | What is included |
|---|---|
| Changing dependent outcome | Directed friendship nominations at three observed waves |
| Opportunity process | Separate change rates for waves 1 to 2 and waves 2 to 3 |
| Structural evaluation terms | Outdegree, reciprocity, transitive triplets, and three-cycles |
| Actor-covariate evaluation term | Exact same-smoking status, measured at wave 1 and held fixed |

The first four structural coefficients ask how a proposed friendship change alters network structure. The same-smoking coefficient asks whether a candidate state receives more weight when it contains more nominations between students in the same baseline smoking category, conditional on those structural terms.

::: {.callout-tip collapse="true" title="Stuck? `interaction1` is not an interaction"}
This argument name is misleading and it trips people up. `interaction1=` does not specify a statistical interaction, it names **which covariate the effect applies to**. `sameX` needs to know *same on what*. In RSiena 1.6 it was renamed `covar1=`, which is a clearer name.
:::

One more thing before we fit, and it is the single most useful line in this section for your own work. Today I dictate every shortName, `recip`, `transTrip`, `cycle3`, `sameX` - but on your data nobody dictates them, and guessing shortNames from memory is how people end up asking "what's the name for indegree popularity again?" At midnight. The answer is that the catalogue is queryable:

```{r}
#| label: b1-effectsdoc
#| eval: false
effectsDocumentation(eff1, type = "html", display = FALSE, filename = "eff1_doc")
```

That writes **`eff1_doc.html`**, a browsable, searchable table of every effect *this particular data object* supports (`r nrow(eff1)` rows for `dat1`), with the `shortName` column being exactly what you feed to `includeEffects()`, plus which `interaction1` covariate each one needs. Called with no effects object at all, `effectsDocumentation()` documents RSiena's entire built-in catalogue. We don't execute it in the render (it writes a file, not console output), run it once tonight and keep the tab open all week. It is the manual's effect chapter, restricted to your data.

Now fit. The fitted object ships with this document, so on your machine it returns instantly; when you refit it from scratch (`options(session.refit = TRUE)`) `siena07` runs its three phases and takes a few seconds. Either way the estimates are identical, the fit is seeded.

```{r}
#| label: b1-fit
alg <- sienaAlgorithmCreate(projname = NULL, seed = 6886, n3 = 1000)
m1 <- fit_or_load("m1_network_only",
  siena07(alg, data = dat1, effects = eff1,
          batch = TRUE, verbose = FALSE, returnDeps = TRUE))
```



```{r}
#| label: b1-print
m1
```

::: {.callout-tip title="Reading that output: Which number is which"}
This is the first `siena07` table, so here is a key. The columns are, left to right, the **estimate**, its **standard error**, and a **convergence** t-ratio (that third column is *not* a significance test, it is the [B2](#b2) diagnostic; you want it small). To assess conventional statistical significance yourself, take estimate ÷ SE: larger than about 2 in absolute value is conventionally "significant." The row to find first is **`reciprocity`**. The later co-evolution model is shown through a tidy table instead ([C1](#c1)); this one is the raw console form on purpose, so you have seen it once.
:::

::: {.callout-important title="Set a seed on every fit. Every one."}
`sienaAlgorithmCreate(seed = 6886)`. This is a method-of-moments estimator with simulation inside. Without a seed, you get different numbers every run and cannot tell "my model is unstable" apart from "Monte Carlo noise." Two seeded runs of this model differ by exactly zero. If you are rerunning this alone next week, that is the difference between knowing you got it right and guessing.

`returnDeps = TRUE` is also not optional, it stores the simulated networks, and [sienaGOF in C2](#c2) needs them. Asking for it after the fact means refitting.
:::

Reciprocity is large and transitivity is clearly positive, friendship is reciprocated and closes triangles, which is what friendship does. Three-cycles is indistinguishable from zero. And smoking homophily, as a fixed wave-1 covariate, is **not** significant.

Hold onto that last one. It is not evidence that smoking doesn't matter for friendship. It's evidence that treating smoking as a fixed characteristic, measured once, never changing, is the wrong model of it. That's exactly what Block C fixes.

::: {.callout-note collapse="true" title="The full toolbox: Covariates your own data will have, interactions, and parallel fits (optional self-study)" #toolbox-b}
The core path used exactly one covariate type, `coCovar`, a fixed actor attribute - because s50 lets us get away with it. Your data won't. The other two covariate constructors, plus the interaction idiom and the parallel-fit flags, in one place. All of this is runnable; none of it is part of the core path.

**A changing actor covariate: `varCovar()`.** An attribute measured at every wave but treated as *exogenous*, it moves, but you're not modeling why. The shape rule is the thing to memorise: **n × (M−1)**, one column per *period*, not per wave, because the covariate's wave-3 value has no period left to act on.

```{r}
#| label: tb-b-varcovar
drink_vc <- varCovar(s50a[, 1:2])   # Waves 1 and 2 only: 2 periods, 2 columns
dim(drink_vc)
```

(This is the halfway house between Block B, drinking frozen at wave 1, and Block C - drinking fully endogenous. Use it when a covariate changes but you have no theory of its change. RSiena centers it for you, same as `coCovar`.)

**A dyadic covariate: `coDyadCovar()`.** A fixed n × n matrix scoring each *pair* - distance, shared classroom, prior co-membership. Here, dissimilarity in wave-1 drinking:

```{r}
#| label: tb-b-dyadcovar
drinkdist <- coDyadCovar(abs(outer(s50a[, 1], s50a[, 1], "-")))
datB <- sienaDataCreate(friend, drink_vc, drinkdist)
effB <- getEffects(datB)
effB <- includeEffects(effB, recip, transTrip)
effB <- includeEffects(effB, egoX, altX, simX, interaction1 = "drink_vc")
effB <- includeEffects(effB, X, interaction1 = "drinkdist")
```

Note the naming asymmetry: actor covariates enter through `egoX`/`altX`/`simX`/`sameX` (which side of the dyad, and how), but a dyadic covariate is already a statement about the pair, so it enters through the single effect **`X`**. (Its time-varying sibling `varDyadCovar()` exists too, pointer in [the closing list](#not-shown).)

**A genuine interaction: `includeInteraction()`.** [B1's stuck-box](#b1) warned you that `interaction1=` is *not* an interaction. This is the function that actually is one - here, "does the reciprocity preference depend on ego's drinking?":

```{r}
#| label: tb-b-interaction
effB <- includeInteraction(effB, egoX, recip, interaction1 = c("drink_vc", ""))
```

The quirk to remember is that `interaction1=` is now a **two-slot vector**, one slot per parent effect in order: `egoX` needs to know its covariate (`"drink_vc"`), `recip` needs nothing (`""`). Get the slots misaligned and you'll create a different interaction than you meant, silently. The fit call itself is unchanged, `siena07(alg, data = datB, effects = effB, ...)` Exactly as on the core path, so we don't re-run it here.

**Parallel estimation.** S50 fits in seconds; your 300-actor, 5-wave dataset will not. `siena07` parallelizes its simulations across cores:

```{r}
#| label: tb-b-parallel
#| eval: false
siena07(alg, data = datB, effects = effB, batch = TRUE, verbose = FALSE,
        useCluster = TRUE, nbrNodes = parallel::detectCores() - 1)
```

Keep the seed in `sienaAlgorithmCreate`, but treat `nbrNodes` as part of the reproducibility recipe and report it, because the simulations are distributed across nodes and the run is only guaranteed to replicate with the same node count. (Expect one byproduct: the worker nodes write a small `cluster.out` log file into your working directory, harmless, delete at will.)
:::

## B2 · Convergence, and how to repair it {#b2}

This section is the one that most changes how you read other people's papers.

There are **two** convergence diagnostics in that output and students routinely conflate them.

```{r}
#| label: b2-diagnostics
cat("overall maximum convergence ratio (tconv.max):", round(m1$tconv.max, 4), "\n")
cat("largest individual t-ratio            :", round(max(abs(m1$tconv)), 4), "\n")
```

**Individual t-ratios** (`m1$tconv`), one per parameter: how far is the simulated mean of this statistic from its observed target, in standard errors. You want every one below **0.1** in absolute value.

**The overall maximum convergence ratio** (`m1$tconv.max`): the maximum over *any linear combination* of the targets. It is strictly harder to satisfy than the individual ones, which is exactly why it exists, parameters can each look fine while a combination of them is badly off.

The manual's thresholds:

> *"Convergence is excellent when the overall maximum convergence ratio is less than 0.2, and for all the individual parameters the t-ratios for convergence all are less than 0.1 in absolute value; convergence is reasonable when the former is less than 0.30. **For published results, it is suggested that estimates presented come from runs in which the overall maximum convergence ratio is less than 0.25.**"*

Our model is at `r round(m1$tconv.max,3)`. That's excellent, and it's easy because the model is small.

### Now let's break it on purpose

A converged model teaches you little about diagnosing convergence. The next example is a deliberately starved run with too few iterations and too few subphases, so you can see what a model that has not converged looks like.

```{r}
#| label: b2-bad
starved <- sienaAlgorithmCreate(projname = NULL, seed = 42, n3 = 200, nsub = 2)
bad <- fit_or_load("m1_starved",
  siena07(starved, data = dat1, effects = eff1, batch = TRUE, verbose = FALSE))
checkpoint(`tconv.max` = bad$tconv.max,
           `max |t|`   = max(abs(bad$tconv)), label = "STARVED RUN")
```



Both diagnostics fail: `tconv.max` is `r round(bad$tconv.max,3)`, above the 0.25 publication threshold, and the worst individual t-ratio is `r round(max(abs(bad$tconv)),3)`, above 0.1. This is a model you must not report.

**The repair is `prevAns`.** You hand the failed fit back to `siena07` as a starting point, it carries forward both the estimates and the derivative matrix, so the new run starts near the answer instead of at the defaults.

```{r}
#| label: b2-repair1
full <- sienaAlgorithmCreate(projname = NULL, seed = 42, n3 = 1000)
rep1 <- fit_or_load("m1_repair1",
  siena07(full, data = dat1, effects = eff1,
          batch = TRUE, verbose = FALSE, prevAns = bad))
checkpoint(`tconv.max` = rep1$tconv.max,
           `max |t|`   = max(abs(rep1$tconv)), label = "REPAIR ROUND 1")
```



One round fixes it: `tconv.max` `r round(rep1$tconv.max,3)`, worst individual t-ratio `r round(max(abs(rep1$tconv)),3)`. Both now pass, though `tconv.max` is only just under 0.25, close enough to the line that I'd run it once more before reporting.

```{r}
#| label: b2-repair2
rep2 <- fit_or_load("m1_repair2",
  siena07(full, data = dat1, effects = eff1,
          batch = TRUE, verbose = FALSE, prevAns = rep1))
checkpoint(`tconv.max` = rep2$tconv.max,
           `max |t|`   = max(abs(rep2$tconv)), label = "REPAIR ROUND 2")
```



And now it's comfortable: `r round(rep2$tconv.max,3)`. The manual's instruction is exactly this loop:

> *"The practical procedure therefore, in usual cases, is that one runs siena, if necessary repeatedly with subsequent runs using the prevAns argument as above, until the Overall maximum convergence ratio is less than 0.25 and all convergence t ratios are less than 0.1."*

::: {.callout-important title="Repair rounds do not monotonically improve, and this surprises people"}
I ran that exact three-step sequence under five different seeds before writing this section. This is an **off-document experiment**: the render only ships the seed-42 fits, so only the seed-42 row below is reproducible from `cache_saom/`. Those are the `m1_starved → m1_repair1 → m1_repair2` numbers displayed in the cached output. To regenerate the other four seeds yourself later, run this exact loop; it refits from scratch and will take some time:

```r
for (s in c(42, 99, 111, 2024, 7)) {
  starved_s <- sienaAlgorithmCreate(projname = NULL, seed = s, n3 = 200, nsub = 2)
  full_s    <- sienaAlgorithmCreate(projname = NULL, seed = s, n3 = 1000)
  b  <- siena07(starved_s, data = dat1, effects = eff1, batch = TRUE, verbose = FALSE)
  r1 <- siena07(full_s,   data = dat1, effects = eff1, batch = TRUE, verbose = FALSE, prevAns = b)
  r2 <- siena07(full_s,   data = dat1, effects = eff1, batch = TRUE, verbose = FALSE, prevAns = r1)
  cat(sprintf("seed %-4d  starved %.3f  repair1 %.3f  repair2 %.3f\n",
              s, b$tconv.max, r1$tconv.max, r2$tconv.max))
}
```

The results I got:

| Seed | Starved | After repair 1 | After repair 2 |
|---|---|---|---|
| 42 (Used above) | 0.343 | 0.236 | **0.128** |
| 99 | 0.328 | 0.207 | **0.253** ↑ |
| 111 | 0.374 | 0.173 | **0.269** ↑ |
| 2024 | 0.149 | 0.096 | **0.355** ↑ |
| 7 | 0.130 | 0.294 ↑ | 0.213 |

(All `tconv.max`.) In three of five, the *second* repair round made things worse, and in one the first round did.

This is not a bug and it is not your fault. `tconv.max` is itself estimated from a finite simulation, so it wanders. The rule is **"run until you are under 0.25,"** not "run until it stops improving", and if a rerun bounces you back over the line, rerun again rather than concluding the model is broken.

Note also seed 7, where round 1 landed at `tconv.max` 0.294 while its worst individual t-ratio was 0.045. That is the two diagnostics disagreeing exactly as described above: every parameter looks converged, and a linear combination of them is not. If you only read the individual t-ratios you would have reported that model.
:::

::: {.callout-warning title="What `prevAns` is and isn't"}
`prevAns` is a **convergence repair**, not a speed trick. Re-running an already-converged model with `prevAns` is not meaningfully faster, I measured 25.9s vs 22.4s cold on the Block C model, which is noise. What it buys you is a better starting point, which is what gets a *failing* model to converge.

If you've seen it described as a big speedup, that's a claim about a specific case, not a general property.
:::

::: {.callout-note collapse="true" title="Depth: What to do when repeated `prevAns` doesn't fix it"}
If you loop three or four times and `tconv.max` won't come down, the problem is usually not the algorithm, it's the specification. The usual culprits, in order:

- **Too many parameters for the data.** Especially on the behavior side; see [C1](#c1), where this is the deck's own historical bug.
- **Near-collinear effects.** E.G. Including `transTrip`, `transTies` and `gwespFF` together; they're competing for the same variance.
- **An effect the data can't inform.** A covariate with almost no variance, or a homophily term on a nearly-constant attribute.
- **The model does not apply**, go back to [A2](#a2) and re-check Jaccard.

Raising `n3` gives phase 3 more simulations and can stabilize the Monte Carlo estimates of standard errors and convergence diagnostics. It does not make the data more informative, does not guarantee smaller standard errors, and does not fix a model that's wrong.
:::

## B3 · Reading a coefficient, on the model we actually fit {#b3}

Most treatments of SAOM interpretation work through hypothetical coefficients. That skips the step people actually get stuck on, so let's do it on `m1`.

Recall the objective function is $f_i(x) = \sum_k \beta_k s_{ik}(x)$. A coefficient tells us how much a network feature raises or lowers the evaluation score of one candidate micro-change, holding the other modeled features constant. A positive coefficient makes changes that increase the statistic more attractive. A negative coefficient makes them less attractive.

Let's take reciprocity. Suppose actor $i$ is deciding whether to add a tie to $j$, and compare two situations: $j$ already nominates $i$, versus $j$ doesn't.

```{r}
#| label: b3-hand
th <- m1$theta
names(th) <- m1$effects$effectName
b_dens  <- th[["outdegree (density)"]]
b_recip <- th[["reciprocity"]]

# Adding a tie always changes outdegree by 1.
# If j already nominates i, it ALSO creates one reciprocated pair.
f_unreciprocated <- b_dens
f_reciprocated   <- b_dens + b_recip

cat("objective change, tie to a non-nominator:", round(f_unreciprocated, 3), "\n")
cat("objective change, tie to a nominator    :", round(f_reciprocated, 3), "\n")
cat("relative choice weight                  :", round(exp(b_recip), 2), "\n")
```

Read this in three steps. First, reciprocity is positive, so a proposed change is more attractive when it creates a returned nomination. Second, name the change: if $j$ already nominates $i$, adding $i\rightarrow j$ raises $i$'s evaluation score by `r round(b_recip,2)`. Third, exponentiate only if the relative weight helps: holding the other modeled features constant, that reciprocating move receives about `r round(exp(b_recip),1)` times the choice weight of an otherwise comparable move that does not create reciprocity.

Now the same logic for closure:

```{r}
#| label: b3-trans
b_tt <- th[["transitive triplets"]]
cat("one new transitive triplet :", round(b_tt, 3), "\n")
cat("closing 3 triplets at once :", round(3 * b_tt, 3),
    "  (relative weight", round(exp(3 * b_tt), 2), ")\n")
```

Again, start with the sign. The positive coefficient says proposed nominations are more attractive when they close friend-of-friend paths. If adding $i\rightarrow j$ closes one path, the score rises by `r round(b_tt,2)` and the move receives about `r round(exp(b_tt),2)` times the relative weight of an otherwise comparable move that closes none. The contribution accumulates: a tie that closes three paths adds $3(`r round(b_tt,2)`)=`r round(3*b_tt,2)`$ to the score. This is the same "count the configurations you create" logic as an ERGM change statistic, with the statistic evaluated from the selected actor's perspective.

### One proposed nomination can change several statistics at once

The previous illustration contained one path from Alice through Ben to Cara. Now expand the network so that it contains two separate paths from Alice to Cara:

$$
\text{Alice}\rightarrow\text{Ben}\rightarrow\text{Cara}
\qquad\text{and}\qquad
\text{Alice}\rightarrow\text{Elena}\rightarrow\text{Cara}.
$$

Also suppose Cara already nominates Alice:

$$
\text{Cara}\rightarrow\text{Alice}.
$$

Alice now considers adding $\text{Alice}\rightarrow\text{Cara}$. This one proposed move adds one outgoing nomination, returns Cara's nomination, and closes both friend-of-friend paths, one through Ben and one through Elena. It therefore creates two transitive triplets. Using the fitted estimates, its combined change in Alice's evaluation score is

```{r}
#| label: b3-combined
combined_move <- c(
  "one more outgoing nomination" = b_dens,
  "one returned nomination" = b_recip,
  "two friend-of-friend paths closed" = 2 * b_tt
)
round(combined_move, 3)
combined_score <- sum(combined_move)
cat("combined change in evaluation score:", round(combined_score, 3), "\n")
cat("relative choice weight versus no change:", round(exp(combined_score), 2), "\n")
```

The negative outdegree contribution says that adding another nomination is unattractive in the fitted process. Returning Cara's nomination and closing the two paths more than offset that cost. The proposed move gains `r round(combined_score,2)` points in the evaluation score and receives about `r round(exp(combined_score),2)` times the choice weight of making no change, assuming the cycle and smoking statistics do not change.

The example also shows why we should not ask whether reciprocity or transitivity "caused" a particular observed friendship. One move can activate both statistics. The reciprocity coefficient compares otherwise similar options that differ in whether the other student already nominated Alice. The transitivity coefficient compares otherwise similar options that differ in how many friend-of-friend paths would close. The model separates their contributions to its choice score, not the private historical reason for one student's decision. If reciprocating moves almost always close paths in these data, the two comparisons contain little independent information and their separate estimates become less precise.

The important warning is that these exponentiated values are not ordinary odds ratios for whether a tie exists. At an opportunity, the actor can consider all eligible additions, all eligible deletions, and no change. The probability assigned to one action depends on that whole choice set. Interpret the exponentiated coefficient as a change in the **relative attractiveness of a candidate micro-step, not the unconditional probability that a tie exists.**

The reusable sentence is: **Holding the other modeled features constant, a one-unit increase in [statistic] changes the evaluation score of this proposed action by [coefficient], making the action more or less attractive.** Add the exponentiated relative weight only when it helps the reader.

### Interpret the same-smoking coefficient as a 0-to-1 change {#smoke-coefficient}

Start with the variable, not the coefficient. The `sameX` statistic is an indicator for the proposed nomination:

| Same-smoking indicator | Meaning |
|---|---|
| \(0\) | The nominating student and possible friend report different smoking categories |
| \(1\) | They report the same smoking category |

Moving from different smoking categories to the same category is therefore a one-unit increase in the predictor. The fitted coefficient is 0.14. Exponentiating it converts the change in the evaluation function into the conditional probability ratio described in the RSiena manual (Ripley et al. 2026):

```{r}
#| label: b3-smoke-coefficient
beta_same_smoking <- unname(
  m1$theta[m1$effects$effectName == "same smoke1"]
)
same_smoking_probability_ratio <- exp(beta_same_smoking)

cat("Same-smoking coefficient:", round(beta_same_smoking, 2), "\n")
cat(
  "Conditional probability ratio:",
  round(same_smoking_probability_ratio, 2), "\n"
)
cat(
  "Percent difference:",
  round(100 * (same_smoking_probability_ratio - 1)), "%\n"
)
```

> At a modeled friendship-change opportunity, a nomination to a student with the same smoking status is about `r round(100 * (same_smoking_probability_ratio - 1))`% more likely to be chosen than an otherwise identical nomination to a student with a different smoking status.

This is the probability statement attached directly to the 0-to-1 predictor change. It is a **conditional probability ratio** between two possible actions in the same choice situation. It is not a 15 percentage-point increase. An absolute percentage-point change would depend on the other nominations, deletions, and the no-change option available to the student at that opportunity.

The estimate is 0.14 with a standard error of 0.15. The point estimate favors a same-smoking nomination, but the uncertainty is large enough that these data do not show a clear difference.

::: {.callout-note collapse="true" title="Optional Detail: Removing the Entire Same-Smoking Coefficient"}

### A different question: remove the entire coefficient {#smoke-first-difference}

This summary compares two parameter settings within the complete network-only specification above. It does not fit a bivariate smoking model and does not change any student's smoking category.

| Model component | Same in both scenarios | What differs |
|---|---|---|
| Starting state | The observed friendship network at the beginning of each interval and every student's wave-1 smoking category | Nothing |
| Opportunity process | The estimated interval-specific friendship-change rates | Nothing |
| Structural evaluation | Outdegree, reciprocity, transitive triplets, and three-cycles at their fitted values | Nothing |
| Smoking evaluation | The same exact-match statistic | \(\beta_{\text{smoke}}=0.14\) in the fitted scenario and \(0\) in the term-removed scenario |

For each scenario, simulate the full process and ask one outcome question: among directed nominations that were absent at the beginning of an interval and joined students in the same baseline smoking category, what proportion are present in the simulated next wave? The difference is a model-based term-ablation summary. It describes what the fitted model produces with and without the same-smoking term. It is not the effect of changing smoking and is not a causal effect.

The simulation follows the fixed-parameter procedure in the [current RSiena manual](https://www.stats.ox.ac.uk/~snijders/siena/RSiena_Manual.pdf) (Ripley et al. 2026). Read the code through five verbs: **configure, copy, simulate, define, average**. We will keep those jobs in separate chunks so that the bookkeeping does not hide the idea.

#### Step 1: Configure a Simulation, Not a New Estimation

`simOnly = TRUE` and `nsub = 0` tell RSiena to generate paths without estimating parameters. `cond = FALSE` lets the interval rates generate the amount of change rather than forcing every simulation to reproduce the observed number of changes. `n3 = 1000` requests 1,000 simulated paths, and the seed makes the two scenarios use the same random-number stream as closely as RSiena permits.

```{r}
#| label: b3-smoke-algorithm
smoke_sim_algorithm <- sienaAlgorithmCreate(
  projname = NULL,
  cond = FALSE,
  useStdInits = FALSE,
  nsub = 0,
  n3 = 1000,
  simOnly = TRUE,
  seed = 9321,
  silent = TRUE
)
```

This object contains instructions for the simulation. It does not contain a fitted model or simulated networks yet.

#### Step 2: Copy the Fitted Coefficients and Change One Number

`updateTheta()` places the fitted evaluation coefficients and interval rates into a simulation-ready effects object. The first object keeps the fitted same-smoking coefficient. The second starts as an exact copy and then replaces only that coefficient with zero.

```{r}
#| label: b3-smoke-scenarios
eff1_smoke_on <- updateTheta(eff1, m1)
eff1_smoke_zero <- setEffect(
  eff1_smoke_on,
  sameX,
  interaction1 = "smoke1",
  initialValue = 0
)
```

The two objects differ in one place: `same smoke1` is 0.14 in `eff1_smoke_on` and 0 in `eff1_smoke_zero`. The starting networks, rates, and four structural coefficients are identical.

#### Step 3: Simulate the Two Versions of the Same Model

Each call begins from the observed start of each interval and stores the 1,000 simulated end-of-interval networks because `returnDeps = TRUE`. The fitted objects are shipped in `cache_saom`, so students reading the rendered walkthrough do not need to rerun the simulations.

```{r}
#| label: b3-smoke-simulations
m1_smoke_on <- fit_or_load(
  "m1_smoke_on_sim",
  siena07(
    smoke_sim_algorithm,
    data = dat1,
    effects = eff1_smoke_on,
    returnDeps = TRUE,
    batch = TRUE,
    verbose = FALSE,
    silent = TRUE
  )
)

m1_smoke_zero <- fit_or_load(
  "m1_smoke_zero_sim",
  siena07(
    smoke_sim_algorithm,
    data = dat1,
    effects = eff1_smoke_zero,
    returnDeps = TRUE,
    batch = TRUE,
    verbose = FALSE,
    silent = TRUE
  )
)
```

#### Step 4: Turn Simulated Networks Into Dyad-Level Probabilities

RSiena stores each simulated network as an edge list. `edge_list_to_matrix()` turns one edge list into a familiar adjacency matrix. `mean_simulated_network()` repeats that conversion for all 1,000 paths and averages the matrices. The resulting cell \([i,j]\) is the proportion of simulated next-wave networks in which \(i\rightarrow j\) exists.

```{r}
#| label: b3-simulation-helpers
edge_list_to_matrix <- function(edges, n) {
  out <- matrix(0L, n, n)
  if (!is.null(edges) && length(edges) > 0L && nrow(edges) > 0L) {
    out[cbind(edges[, 1], edges[, 2])] <- edges[, 3]
  }
  out
}

mean_simulated_network <- function(fit, period, n, variable = "friend") {
  simulated_networks <- lapply(
    fit$sims,
    function(run) edge_list_to_matrix(run[[1]][[variable]][[period]], n)
  )
  Reduce(`+`, simulated_networks) / length(simulated_networks)
}
```

For example, `mean_simulated_network(m1_smoke_on, 1, 50)` returns a \(50\times50\) matrix for waves 1 to 2. It is not one simulated network. It is a matrix of simulation frequencies, which we use as model-implied next-wave probabilities.

#### Step 5: Define the Risk Set, Subtract, and Average

The risk set contains only ordered pairs that satisfy three conditions: the pair is not a self-tie, the nomination was absent at the start of the interval, and the two students shared the same wave-1 smoking category. Existing nominations are excluded because they answer a persistence question rather than a formation question.

The loop below handles one interval at a time. For each interval it creates the two probability matrices, selects the eligible cells, and averages the cell-by-cell difference.

```{r}
#| label: b3-smoke-first-difference

n_students <- dim(s50_friend_array)[1]
off_diagonal <- row(s50_friend_array[, , 1]) !=
  col(s50_friend_array[, , 1])
same_smoking <- outer(s50s[, 1], s50s[, 1], `==`)

smoke_rows <- vector("list", 2)
for (period in 1:2) {
  p_on <- mean_simulated_network(m1_smoke_on, period, n_students)
  p_zero <- mean_simulated_network(m1_smoke_zero, period, n_students)

  start_network <- s50_friend_array[, , period]
  risk_set <- off_diagonal & start_network == 0 & same_smoking

  smoke_rows[[period]] <- data.frame(
    interval = paste0("Waves ", period, " to ", period + 1L),
    eligible_dyads = sum(risk_set),
    smoke_term_zero = mean(p_zero[risk_set]),
    smoke_term_fitted = mean(p_on[risk_set]),
    first_difference = mean(p_on[risk_set] - p_zero[risk_set])
  )
}
smoke_first_difference <- do.call(rbind, smoke_rows)
```

The last line stacks the two interval summaries. We then pool them using the number of eligible nominations as weights, so an interval with more eligible same-smoking nominations contributes proportionally more to the overall result.

```{r}
#| label: b3-smoke-pool

smoke_pooled <- data.frame(
  interval = "Pooled across intervals",
  eligible_dyads = sum(smoke_first_difference$eligible_dyads),
  smoke_term_zero = weighted.mean(
    smoke_first_difference$smoke_term_zero,
    smoke_first_difference$eligible_dyads
  ),
  smoke_term_fitted = weighted.mean(
    smoke_first_difference$smoke_term_fitted,
    smoke_first_difference$eligible_dyads
  ),
  first_difference = weighted.mean(
    smoke_first_difference$first_difference,
    smoke_first_difference$eligible_dyads
  )
)

smoke_summary <- rbind(smoke_first_difference, smoke_pooled)
knitr::kable(
  transform(
    smoke_summary,
    smoke_term_zero = scales::percent(smoke_term_zero, accuracy = 0.01),
    smoke_term_fitted = scales::percent(smoke_term_fitted, accuracy = 0.01),
    first_difference = paste0(round(100 * first_difference, 2), " pp")
  ),
  row.names = FALSE,
  caption = "Model-based first difference for the same-smoking term"
)
```

The reusable interpretation is:

> At the fitted point estimates, including the same-smoking term raises the model-implied probability that an initially absent same-smoking nomination exists at the next wave from `r scales::percent(smoke_pooled$smoke_term_zero, accuracy=.01)` to `r scales::percent(smoke_pooled$smoke_term_fitted, accuracy=.01)`, an increase of `r round(100 * smoke_pooled$first_difference, 2)` percentage points, with the starting networks and all other coefficients held fixed.

That is much easier to read than the coefficient alone. It is also small, which matches the statistical result: the same-smoking estimate is 0.14 with a standard error of 0.15, so the data do not show a clear smoking-similarity pattern. This table uses the point estimates; a research paper should repeat the simulations over parameter draws to attach an uncertainty interval to the first difference.

:::

::: {.callout-note collapse="true" title="More Detail: Why This Is a Model-Based First Difference, Not an Ordinary AME" #saom-first-differences}
A coefficient interpretation and a simulation-based first difference answer different questions. The coefficient tells us how a feature changes the score of one proposed micro-step. A simulated first difference asks how a clearly defined change alters a recognizable outcome after the modeled network process unfolds. For substantive magnitude, the second quantity is often more useful.

There is no single generic "probability of a tie" in a SAOM. Before simulating anything, define the estimand:

1. **Local micro-step choice:** when actor $i$ receives an opportunity, how does the focal contrast change the relative attractiveness or choice probability of one candidate action?
2. **Tie formation over an interval:** among ties absent at wave $t$, how does the contrast change the simulated chance that $i\rightarrow j$ exists at wave $t+1$?
3. **Expected network or behavior outcome:** how does the contrast change expected outdegree, reciprocated ties, same-attribute ties, or an actor's behavior by the next wave?

For most applied questions, the second or third estimand is easier to communicate because it incorporates the full sequence of modeled opportunities and choices between waves. It also avoids pretending that one coefficient is itself a probability.

The workflow is:

1. Begin from the same observed network and behavior state at wave $t$.
2. Simulate many paths to wave $t+1$ under the baseline values.
3. Make one substantively defensible change while holding the other starting information fixed.
4. Simulate the paths again, ideally using common random-number streams so the paired difference contains less Monte Carlo noise.
5. Calculate the outcome in both scenarios, subtract baseline from the counterfactual, and average over the actors, dyads, and intervals that answer the research question.

For friendship formation, average only over ties that are absent and eligible at the start of the interval. Already-existing ties face persistence or dissolution, which is a different outcome and should be summarized separately. For example, an interval-level formation estimand is

$$
\frac{1}{|\mathcal D_0|}
\sum_{(i,j)\in\mathcal D_0}
\left[
\widehat{\Pr}\!\left(Y_{ij,t+1}=1\mid\text{counterfactual}\right)
-
\widehat{\Pr}\!\left(Y_{ij,t+1}=1\mid\text{baseline}\right)
\right],
$$

where $\mathcal D_0$ contains eligible dyads with $Y_{ij,t}=0$. A useful interpretation would be: "Among friendships absent at the start of the interval, the specified contrast changed the simulated chance of a friendship at the next wave by an average of X percentage points."

The counterfactual must make sense for the variable. If similarity is constructed from actor attributes, changing one actor's value can affect that actor's similarity with several possible partners. We cannot generally flip one dyad's similarity independently while pretending the actor attributes stayed unchanged. Setting a coefficient to zero is another possible simulation contrast, but it answers a different question: what changes when the fitted mechanism is removed, not what happens when an observed attribute changes.

For an actor covariate, a more natural outcome may be expected outdegree or the expected number of same-attribute friends. For a behavior influence effect, one could compare the simulated probability that an actor moves up a behavior category, the expected behavior level, or the probability of remaining unchanged. The outcome should follow the research question rather than defaulting to a probability for one arbitrary dyad.

Uncertainty requires more than repeated paths at $\hat\theta$. Repeat the paired baseline and counterfactual simulations over parameter draws, using the fitted covariance matrix as an approximation when appropriate, and summarize the resulting first differences with an interval. This incorporates parameter uncertainty as well as dynamic simulation uncertainty.

The reporting pair is:

- **Local mechanism:** "Holding the other modeled changes constant, [feature] raises or lowers the evaluation score of [specific candidate action] by [coefficient]."
- **Simulation-based magnitude:** "Across [defined risk set], changing [baseline] to [counterfactual] changed the simulated [outcome] by an average of [difference], with [uncertainty interval]."

The second quantity is model-based, not automatically causal. A causal interpretation additionally requires a coherent intervention, defensible adjustment for confounding, correct temporal and network-process assumptions, and a population over which the average is meaningful.
:::

::: {.callout-note collapse="true" title="Depth: Why these are not ERGM coefficients"}
Yesterday's ERGM reciprocity coefficient and today's SAOM reciprocity coefficient have the same name and are not the same quantity. The ERGM one is a contribution to the conditional log-odds of a tie **given the entire rest of the graph** (a static distribution over graphs, no dynamics, no equilibrium enter into it). The SAOM one is a contribution to **one actor's objective function at one micro-step**, where the actor is choosing among their own outgoing ties only.

Block, Stadtfeld & Snijders (2019, *SMR* 48(1):202–239) make the general version of the point: ERGM is specified at the **graph** level, SAOM at the **transition** level, and SAOM's statistics nest within actors, which means its dependence assumptions are *stronger*, not weaker. Do not put the two coefficients in the same table and compare them.
:::

### Quick exercise {#exercise-b}

::: {.panel-tabset}

## Core

The model has a `same smoke1` effect, homophily on smoking. Convert it to a relative choice weight the way we just did for reciprocity, and say in one sentence what it would mean *if* it were distinguishable from zero.

```{r}
#| label: ex-b-core
#| eval: false
b_smoke <- th[["same smoke1"]]
exp(______)
```

## Stretch

`transTrip` counts transitive triplets, so its contribution scales with how many triplets a tie closes. Compute the objective-function contribution of a tie that closes 0, 1, 2 and 5 triplets, and plot it. Then explain why a modeller might prefer `gwespFF`, the geometrically weighted version, and what specifically it changes about the shape you just plotted.

## Solution

```{r}
#| label: ex-b-solution
#| fig-height: 3.4
b_smoke <- th[["same smoke1"]]
cat("same-smoking relative choice weight:", round(exp(b_smoke), 2),
    " (t =", round(b_smoke / sqrt(diag(m1$covtheta))[m1$effects$effectName == "same smoke1"], 2), ")\n")

k <- 0:5
contrib <- b_tt * k
plot(k, contrib, type = "b", pch = 19,
     xlab = "transitive triplets closed by the tie",
     ylab = "contribution to objective function",
     main = "transTrip is linear in the number of triplets closed")
```

The relative weight is a little above 1, so the point estimate leans toward nominations of students with the same baseline smoking status receiving more weight. With a t-ratio well under 2, however, we cannot distinguish it from no effect. The right sentence is "these data do not show a clear baseline smoking-similarity pattern," not "smoking does not matter." As [C1](#c1) shows, treating smoking as fixed also cannot answer an influence question.

On the plot: `transTrip` is **linear** in triplets closed, so the 5th triplet is worth as much as the 1st. `gwespFF` bends that curve so additional closure has diminishing returns, which is usually more plausible and, exactly as with `gwesp` in yesterday's ERGM, much better behaved, because a linear term rewards runaway clustering.
:::

```{r}
#| label: b3-checkpoint
checkpoint(`recip` = b_recip, `transTrip` = b_tt,
           `tconv.max` = m1$tconv.max, label = "END BLOCK B")
```

```{r}
#| label: offroad-1
# Optional extension space for a specification question.
# Any exploratory fit belongs here so the main fitted objects remain unchanged.

```

---

# Block C: Co-evolution, and what it does not buy you {#block-c}

## C1 · Selection and influence, jointly {#c1}

That distinction is central to Steglich, Snijders, and Pearson (2010).

The problem: we observe that friends resemble each other in drinking. Two mechanisms produce that, and they are observationally identical in a cross-section.

- **Selection**, people choose friends who already drink like them.
- **Influence**, people come to drink like the friends they already have.

The distinction is easiest to see as two possible histories for the same endpoint. Under selection, Alice and Cara report similar drinking at an earlier interview and become friends afterward. Under influence, they are already friends and Alice's later drinking moves toward Cara's. The dependent outcome is therefore different in the two equations:

| Process | Earlier information | Later change being modeled |
|---|---|---|
| Selection | Alice's and Cara's drinking levels | Whether Alice adds or keeps a friendship nomination to Cara |
| Influence | Alice's existing friends and their drinking levels | Whether Alice's own reported drinking moves up, down, or stays unchanged |

The waves help establish which variables were observed earlier, but they do not reveal the exact order of every event between interviews. The continuous-time micro-step process is the model's account of that missing path.

If you only have one wave you cannot tell these apart. With several waves, a network-only model can estimate whether similarity predicts tie change but cannot represent subsequent behavior convergence in the same process. A behavior-only model can estimate whether friends' prior behavior predicts change but cannot represent the selection of similar friends. Steglich et al.'s argument is that the two equations should be estimated **jointly** because each process can generate the pattern the other equation is trying to explain.

So now drinking becomes a *second dependent variable*, evolving by the same micro-step logic as the network.

```{r}
#| label: c1-build
friend <- sienaDependent(s50_friend_array)
drink  <- sienaDependent(array(s50a, dim = c(50, 1, 3)), type = "behavior")
smoke1 <- coCovar(s50s[, 1])
dat2   <- sienaDataCreate(friend, drink, smoke1)
dat2
```

Note `type = "behavior"` and the `c(50, 1, 3)` shape, 50 actors, 1 variable, 3 waves. Getting that array shape wrong is the most common error in setting up a co-evolution model.

The effects come in two groups, and the grouping *is* the identification story:

```{r}
#| label: c1-effects
eff2 <- getEffects(dat2)

#: Network equation: structure, then SELECTION on drinking ---
eff2 <- includeEffects(eff2, recip, transTrip, cycle3)
eff2 <- includeEffects(eff2, egoX, altX, simX, interaction1 = "drink")

#: Behavior equation: shape, then INFLUENCE from friends ---
eff2 <- includeEffects(eff2, name = "drink", avAlt, interaction1 = "friend")
```

### The Complete Co-Evolution Specification {#c1-specification}

This is a different model from the fixed-smoking specification in Block B. Both friendship \(x\) and drinking \(z\) are now changing dependent variables. Each process has separate period-specific rate parameters and its own evaluation function. Schematically,

$$
\begin{aligned}
f_i^{\text{friend}}(x,z)
&=
\beta_{\text{out}}s_{\text{out},i}
+\beta_{\text{rec}}s_{\text{rec},i}
+\beta_{\text{trans}}s_{\text{trans},i}
+\beta_{\text{cycle}}s_{\text{cycle},i} \\
&\quad
+\beta_{\text{ego}}s_{\text{drink ego},i}
+\beta_{\text{alter}}s_{\text{drink alter},i}
+\beta_{\text{sim}}s_{\text{drink similarity},i},
\\[4pt]
f_i^{\text{drink}}(z,x)
&=
\gamma_{\text{linear}}s_{\text{linear},i}
+\gamma_{\text{quadratic}}s_{\text{quadratic},i}
+\gamma_{\text{avAlt}}s_{\text{average friend drinking},i}.
\end{aligned}
$$

The notation is schematic: every statistic is evaluated for a candidate next state. At a friendship opportunity, student \(i\) can add one outgoing nomination, delete one, or make no change. At a drinking opportunity, student \(i\) can move one response category up, move one category down, or make no change, subject to the scale boundaries.

| Equation | Terms included before the focal interpretation |
|---|---|
| Friendship change | Outdegree, reciprocity, transitive triplets, three-cycles, drinking ego, drinking alter, and drinking similarity |
| Drinking change | Linear shape, quadratic shape, and average friends' drinking |

On the network side, `simX` on drinking is **selection**: do actors prefer friends whose drinking resembles their own? In Block B we used `sameX` for smoking, not `simX`, because `sameX` is an exact match for a categorical attribute while `simX` is graded similarity for a numeric variable such as drinking, which runs from 1 to 5. `egoX` and `altX` go alongside it because otherwise a similarity term can pick up the fact that students at some drinking levels simply nominate more peers or receive more nominations.

On the behavior side, `avAlt` is **influence**: does an actor move their drinking toward the average of their friends'? The two shape terms provide the baseline for drinking change before we ask about friends. The **linear shape** term captures an overall tendency to move toward higher rather than lower drinking categories. The **quadratic shape** term captures whether movement tends toward the middle of the scale or toward its extremes. These are not linear and quadratic trends over time. (`avAlt` has siblings: `totAlt` uses the *total* rather than the average of friends' values, so more friends means more exposure, and `avSim`/`totSim` reward *resembling* friends rather than moving toward their level. The av-vs-tot choice is about whether influence scales with degree; the alt-vs-sim choice is the [stretch exercise](#exercise) below.)

```{r}
#| label: c1-fit
m2 <- fit_or_load("m2_coevolution",
  siena07(alg, data = dat2, effects = eff2,
          batch = TRUE, verbose = FALSE, returnDeps = TRUE))
```



```{r}
#| label: c1-print
m2
```

```{r}
#| label: c1-table
th2 <- m2$theta; se2 <- sqrt(diag(m2$covtheta))
res <- data.frame(effect = m2$effects$effectName,
                  est = round(th2, 4), se = round(se2, 4),
                  t = round(th2 / se2, 2))
key <- res[grepl("similarity|average alter|reciprocity|transitive", res$effect), ]
knitr::kable(key, row.names = FALSE,
             caption = "the two rows this whole session is about, plus structure for context")
```

(We built that table by hand from `$theta` and `$covtheta` because you should know where the numbers live. RSiena will also format a publication table for you, `siena.table()`, demonstrated in [the C3 toolbox](#toolbox-c).)

The fitted results are clearer if we keep the two dependent outcomes visible:

| Equation | Question | Estimate | SE | What these data show |
|---|---|---:|---:|---|
| Friendship | Are students with more similar drinking more attractive friendship choices? | `r round(th2[m2$effects$effectName=="drink similarity"],2)` | `r round(se2[m2$effects$effectName=="drink similarity"],2)` | A positive selection-consistent pattern |
| Drinking | Does a student's drinking move toward the average among current friends? | `r round(th2[m2$effects$effectName=="drink average alter"],2)` | `r round(se2[m2$effects$effectName=="drink average alter"],2)` | A positive point estimate with enough uncertainty to include no influence |

The first row models friendship change. Candidate friendship states receive greater weight when two students report more similar drinking. The second row models behavior change. Its estimate leans toward students moving in the direction of their friends' average drinking, but the uncertainty is too wide to distinguish that pattern from no influence in this sample.

The tempting sentence is "so it's selection, not influence." **Do not write that sentence.** Here is why not.

### What the Result Does and Does Not Establish {#selection-influence-limits}

**1. Joint estimation is an identifying assumption, not a proof.** Estimating both equations at once means selection and influence are separately *parameterized*. It does not by itself mean the data cleanly separate them. The distinction relies on the continuous-time micro-step structure: the model represents one actor changing one tie or one behavior step at a time. If Alice's friendship and drinking both change between two interviews, the panel does not reveal which happened first, so the model must allocate that pattern through its assumed process.

**2. Both mechanisms are badly underpowered at realistic sample sizes.** We have 50 actors and 3 waves, which is a *normal* size for this literature. Look at the influence standard error relative to the estimate. A null result here is weak evidence of absence.

::: {.callout-note collapse="true" title="Depth: Why the selection/influence balance is fragile at this N"}
This is worth knowing before you over-read our result. We got selection significant and influence not, but at n = 50 that balance is not stable. Flip the datasets, the waves, or the behavior coding and the two can swap which one clears significance, precisely because both equations are estimated from very little information. The other canonical teaching dataset here, Knecht's classroom data (26 pupils, 4 waves), is smaller still and is exactly where you would expect the estimates to be even noisier.

The defensible reading is the identification point, not either specific verdict: "SAOM separates selection from influence" means the model separately *parameterizes* them, not that either is reliably *estimable* at these sample sizes. If you want to see this for yourself, the ancillary TERGM deep dive, the long-form companion to yesterday's C4, fits the Knecht data (it ships with `btergm` as `data(knecht)`, not with RSiena) and you can add the behavior equation there and watch the standard errors.
:::

**3. The behavior equation is where people over-fit, and it matters here.** Our behavior equation has five parameters: two rate parameters plus linear shape, quadratic shape, and `avAlt`. That is close to the maximum a 1–5 ordinal outcome on 50 actors will support. A common failure, and it appears in widely-circulated teaching material for this exact dataset - is piling on more behavior effects and getting an influence estimate with a **standard error larger than the estimate itself**, then reporting it without comment. Keeping the behavior side minimal is not laziness; it's what makes the influence coefficient readable at all.

**4. Unobserved time-varying confounders are untouched.** If joining a new peer group or experiencing a family change affects both friendship formation and drinking, that omitted common cause can appear as selection or influence. This is Shalizi & Thomas again, and having panel data narrows the problem but does not close it.

**5. Wave spacing is a modeling assumption.** The rate function absorbs it. If the meaningful friendship and drinking changes happen much faster than the interval between interviews, we do not observe their order. If they happen much more slowly, the waves may contain little useful change.

**6. Measurement still matters.** A student can omit a friendship nomination or report drinking imprecisely. Joint estimation does not correct those errors automatically. Separate equations make the two mechanisms explicit, but they do not make these alternative explanations disappear.

### Interpret the drinking-similarity coefficient step by step {#selection-interpretation}

Start by naming the coefficient and the outcome. The estimate of `r round(th2[m2$effects$effectName=="drink similarity"],2)` belongs to the **friendship equation**. It describes how drinking similarity changes the relative attractiveness of one possible friendship nomination. It is not a coefficient in the drinking equation.

**Step 1: Put the predictor on a scale we can see.** Reported drinking is coded 1 for none, 2 for once or twice a year, 3 for once a month, 4 for once a week, and 5 for more than once a week. RSiena turns the absolute gap between two students into a similarity score running from 0 to 1. Students in the same category have similarity 1. Students at opposite ends of the scale have similarity 0.

```{r}
#| label: c1-similarity-scale
similarity_scale <- data.frame(
  drinking_gap = 0:4,
  similarity = 1 - (0:4 / 4)
)
knitr::kable(
  similarity_scale,
  row.names = FALSE,
  col.names = c(
    "Difference in reported drinking categories",
    "Drinking-similarity score"
  )
)
```

**Step 2: Define a one-category change in the predictor.** Moving a possible friend one drinking category closer reduces the gap by one and raises the similarity score by \(1/4=0.25\). This is the predictor change we want to interpret.

**Step 3: Multiply that predictor change by the coefficient.** The drinking-similarity coefficient is `r round(th2[m2$effects$effectName=="drink similarity"],2)`, so moving one category closer changes the evaluation score of that possible nomination by

```{r}
#| label: c1-selection-one-category
beta_similarity <- unname(
  th2[m2$effects$effectName == "drink similarity"]
)
one_category_similarity_change <- 1 / 4
score_change_one_category <-
  beta_similarity * one_category_similarity_change
relative_choice_weight_one_category <-
  exp(score_change_one_category)

cat(
  "Change in similarity score:",
  round(one_category_similarity_change, 2), "\n"
)
cat(
  "Change in nomination's evaluation score:",
  round(score_change_one_category, 2), "\n"
)
cat(
  "Relative choice weight:",
  round(relative_choice_weight_one_category, 2), "\n"
)
```

**Step 4: State what the model says about the application.**

> At a modeled friendship-change opportunity, a student who is one drinking category closer is about `r round(100 * (relative_choice_weight_one_category - 1))`% more likely to be chosen than an otherwise comparable student.

This is the closest SAOM analogue to the familiar regression sentence, “a one-unit increase in X changes Y by the coefficient.” Here the one-unit substantive change is **one drinking category closer**, and “43% more likely” is a conditional probability ratio of 1.43 between otherwise comparable possible nominations. It is not a 43 percentage-point increase. An absolute percentage-point change would depend on every other available nomination, deletion, and the option of making no change.

**Step 5: Check the pattern across the whole drinking scale.** Raw `simX` contributions are centered inside RSiena. Centering changes the common reference point but not the one-category comparison above. `selectionTable()` combines the relevant drinking terms and shows the fitted pattern for every ego-alter pairing:

```{r}
#| label: c1-seltable
st <- selectionTable(m2, xd = dat2, name = "friend", vname = "drink",
                     levls = 1:5, as.matrix = TRUE, silent = TRUE)
round(st, 3)
```

Rows are the nominating student's drinking level and columns are the possible friend's drinking level. Read across one row while holding the nominating student fixed. Higher cells mark candidates who receive more weight. The diagonal cells describe candidates who report the same drinking level.

```{r}
#| label: c1-inftable
it <- influenceTable(m2, xd = dat2, netname = "friend", behname = "drink",
                     levls = 1:5, as.matrix = TRUE, silent = TRUE)
round(it, 3)
```

Here rows are the average drinking of ego's friends and columns are ego's own possible next value. Read across the top row, friends who barely drink, and ego's objective is maximized at a low value. Read across the bottom row and it's maximized at a high value. That is the influence effect made visible; remember it's the one that *wasn't* significant, so this is the shape of an effect we cannot confidently distinguish from zero.

These interpretation tables display contributions to the evaluation function. Read higher cells as candidate states the fitted model scores as more attractive and lower cells as states it scores as less attractive. They do not describe how often those states occur in the observed panel.

### Optional model-level summaries {#average-first-differences}

The one-category comparison above is the direct interpretation of the drinking-similarity coefficient. The simulations below answer a different question: how do entire next-wave outcomes change if we remove one fitted mechanism from the model? These are useful summaries of the fitted process, but they are not one-unit changes in a predictor and should not be described as marginal effects.

| Scenario | Friendship equation | Drinking equation |
|---|---|---|
| Complete fitted model | All structural, drinking ego, drinking alter, and drinking-similarity terms at their fitted values | Shape and average-friend terms at their fitted values |
| Selection term removed | Same as the complete model except \(\beta_{\text{sim}}=0\) | Unchanged |
| Influence term removed | Unchanged | Same as the complete model except \(\gamma_{\text{avAlt}}=0\) |

Every scenario begins from the same observed friendship and drinking values and retains the same period-specific rate parameters. A **model-based average first difference** averages the simulated outcome difference between the complete model and one term-removed version over a clearly defined set of students or nominations. This is a term-ablation summary. It asks what the fitted process produces with and without one coefficient. It does not change an observed characteristic, so it is not an ordinary average marginal effect, an inverse-logit transformation, or a causal effect.

::: {.callout-note collapse="true" title="More Detail: Removing the Entire Selection Coefficient"}

The code follows the same five-step recipe as the smoking example. The only difference is that we now need three scenarios because the joint model contains two focal mechanisms.

#### Step 1: Create the Complete and Term-Removed Specifications

`updateTheta()` copies every fitted coefficient and rate into `eff2_fitted`. We then make two copies. In one, only drinking-based friendship selection is set to zero. In the other, only movement toward friends' average drinking is set to zero.

```{r}
#| label: c1-first-difference-scenarios
coevolution_sim_algorithm <- sienaAlgorithmCreate(
  projname = NULL,
  cond = FALSE,
  useStdInits = FALSE,
  nsub = 0,
  n3 = 1000,
  simOnly = TRUE,
  seed = 7312,
  silent = TRUE
)

eff2_fitted <- updateTheta(eff2, m2)
eff2_no_selection <- setEffect(
  eff2_fitted,
  simX,
  interaction1 = "drink",
  initialValue = 0
)
eff2_no_influence <- setEffect(
  eff2_fitted,
  name = "drink",
  avAlt,
  interaction1 = "friend",
  initialValue = 0
)
```

Read those objects as a controlled comparison:

- `eff2_fitted`: complete fitted friendship and drinking equations.
- `eff2_no_selection`: same model, but the **coefficient on drinking similarity** is set from its fitted value to 0.
- `eff2_no_influence`: same model, but `drink average alter = 0`.

No starting observation is edited. No other coefficient is re-estimated.

#### Step 2: Simulate All Three Specifications

Each call stores 1,000 simulated paths from the observed start of each interval. The same algorithm object and seed keep the simulation settings aligned across the three scenarios.

```{r}
#| label: c1-first-difference-simulations
m2_full_sim <- fit_or_load(
  "m2_full_sim",
  siena07(
    coevolution_sim_algorithm,
    data = dat2,
    effects = eff2_fitted,
    returnDeps = TRUE,
    batch = TRUE,
    verbose = FALSE,
    silent = TRUE
  )
)
m2_no_selection_sim <- fit_or_load(
  "m2_no_selection_sim",
  siena07(
    coevolution_sim_algorithm,
    data = dat2,
    effects = eff2_no_selection,
    returnDeps = TRUE,
    batch = TRUE,
    verbose = FALSE,
    silent = TRUE
  )
)
m2_no_influence_sim <- fit_or_load(
  "m2_no_influence_sim",
  siena07(
    coevolution_sim_algorithm,
    data = dat2,
    effects = eff2_no_influence,
    returnDeps = TRUE,
    batch = TRUE,
    verbose = FALSE,
    silent = TRUE
  )
)
```

#### Step 3: Define the Selection Outcome and Its Risk Sets

For selection, the outcome is whether a nomination absent at the start exists in the simulated next wave. We sort these initially absent nominations into three groups using the students' drinking values at the beginning of that interval: same level, one category apart, or at least two categories apart.

We keep the groups separate because the centered similarity statistic raises the score for similar candidates and lowers it for dissimilar candidates. A single average across every absent nomination could hide both changes through cancellation.

```{r}
#| label: c1-selection-first-differences
gap_labels <- c(
  "Same reported drinking",
  "One category apart",
  "Two or more categories apart"
)
selection_rows <- list()
row_number <- 1

for (period in 1:2) {
  p_fitted <- mean_simulated_network(m2_full_sim, period, n_students)
  p_no_selection <- mean_simulated_network(
    m2_no_selection_sim, period, n_students
  )

  start_network <- s50_friend_array[, , period]
  drinking_gap <- abs(outer(s50a[, period], s50a[, period], "-"))
  initially_absent <- off_diagonal & !is.na(start_network) &
    start_network == 0

  gap_groups <- list(
    drinking_gap == 0,
    drinking_gap == 1,
    drinking_gap >= 2
  )

  for (group_number in seq_along(gap_groups)) {
    risk_set <- initially_absent & gap_groups[[group_number]]

    selection_rows[[row_number]] <- data.frame(
      interval = period,
      starting_gap = gap_labels[group_number],
      nominations = sum(risk_set),
      selection_term_zero = mean(p_no_selection[risk_set]),
      fitted_model = mean(p_fitted[risk_set]),
      first_difference = mean(
        p_fitted[risk_set] - p_no_selection[risk_set]
      )
    )
    row_number <- row_number + 1
  }
}
selection_first_differences <- do.call(rbind, selection_rows)
```

The outer loop moves through the two intervals. The inner loop moves through the three starting drinking gaps. Inside the loop, `risk_set` is a logical matrix marking exactly which ordered pairs belong in that comparison.

#### Step 4: Pool Like With Like Across Intervals

We combine the same drinking-gap group across the two intervals. The weights are the numbers of initially absent nominations in each interval, not the numbers of simulations.

```{r}
#| label: c1-selection-pool
selection_pooled_rows <- vector("list", length(gap_labels))
for (group_number in seq_along(gap_labels)) {
  interval_rows <- selection_first_differences[
    selection_first_differences$starting_gap == gap_labels[group_number],
  ]

  selection_pooled_rows[[group_number]] <- data.frame(
    starting_gap = gap_labels[group_number],
    nominations = sum(interval_rows$nominations),
    selection_term_zero = weighted.mean(
      interval_rows$selection_term_zero,
      interval_rows$nominations
    ),
    fitted_model = weighted.mean(
      interval_rows$fitted_model,
      interval_rows$nominations
    ),
    first_difference = weighted.mean(
      interval_rows$first_difference,
      interval_rows$nominations
    )
  )
}
selection_pooled <- do.call(rbind, selection_pooled_rows)

knitr::kable(
  transform(
    selection_pooled,
    selection_term_zero = scales::percent(
      selection_term_zero, accuracy = 0.01
    ),
    fitted_model = scales::percent(fitted_model, accuracy = 0.01),
    first_difference = scales::percent(
      first_difference, accuracy = 0.01
    )
  ),
  row.names = FALSE,
  col.names = c(
    "Starting drinking gap",
    "Initially absent nominations",
    "Drinking-similarity coefficient set to zero",
    "Fitted model",
    "Average first difference"
  ),
  caption = "What the drinking-similarity selection term changes"
)
```

This is a whole-coefficient comparison, not the interpretation of a one-category change in drinking similarity. Removing the 1.44 coefficient lowers the model-implied next-wave probability for initially absent nominations between same-drinking students from `r scales::percent(selection_pooled$fitted_model[selection_pooled$starting_gap=="Same reported drinking"], accuracy=.01)` to `r scales::percent(selection_pooled$selection_term_zero[selection_pooled$starting_gap=="Same reported drinking"], accuracy=.01)`. Use the one-category comparison in the main section when interpreting the coefficient itself.

:::

#### Step 5: Define and Average the Influence Outcome

For influence, “reports higher drinking” is the wrong outcome because movement toward friends can mean moving either up or down. Instead, take students who begin an interval with at least one outgoing friend and whose drinking differs from those friends' average. For every simulated next wave, calculate how far the student's reported drinking ends from that fixed starting friend average. A smaller distance means the simulated behavior ends closer to where the student's friends began.

The first helper calculates each student's average friend drinking at the observed start of an interval. The second extracts the simulated next-wave drinking values and averages their absolute distance from that starting friend average.

```{r}
#| label: c1-influence-first-difference
mean_starting_friend_drinking <- function(network, behavior) {
  vapply(
    seq_len(nrow(network)),
    function(i) {
      friends <- !is.na(network[i, ]) & network[i, ] == 1
      if (any(friends)) mean(behavior[friends]) else NA_real_
    },
    numeric(1)
  )
}

mean_distance_from_starting_friends <- function(
  fit, period, friend_average, risk_set
) {
  simulated_end <- sapply(
    fit$sims,
    function(run) run[[1]][["drink"]][[period]]
  )
  mean(abs(
    simulated_end[risk_set, , drop = FALSE] - friend_average[risk_set]
  ))
}

influence_rows <- vector("list", 2)
for (period in 1:2) {
  friend_average <- mean_starting_friend_drinking(
    s50_friend_array[, , period],
    s50a[, period]
  )

  risk_set <- is.finite(friend_average) &
    abs(s50a[, period] - friend_average) > sqrt(.Machine$double.eps)

  no_influence_distance <- mean_distance_from_starting_friends(
    m2_no_influence_sim, period, friend_average, risk_set
  )
  fitted_distance <- mean_distance_from_starting_friends(
    m2_full_sim, period, friend_average, risk_set
  )

  influence_rows[[period]] <- data.frame(
    interval = paste0("Waves ", period, " to ", period + 1L),
    student_intervals = sum(risk_set),
    influence_term_zero = no_influence_distance,
    fitted_model = fitted_distance,
    first_difference = fitted_distance - no_influence_distance
  )
}
influence_first_differences <- do.call(rbind, influence_rows)
```

The risk set now contains students rather than dyads. It keeps students who start with at least one outgoing friend and are not already exactly at their friends' average. The interval difference is “distance with influence included” minus “distance with influence removed,” so a negative number means the fitted influence term moves simulated drinking closer to the friends' starting average.

Pool the two interval summaries using the number of eligible student-intervals as weights:

```{r}
#| label: c1-influence-pool
influence_pooled <- data.frame(
  interval = "Pooled across intervals",
  student_intervals = sum(influence_first_differences$student_intervals),
  influence_term_zero = weighted.mean(
    influence_first_differences$influence_term_zero,
    influence_first_differences$student_intervals
  ),
  fitted_model = weighted.mean(
    influence_first_differences$fitted_model,
    influence_first_differences$student_intervals
  ),
  first_difference = weighted.mean(
    influence_first_differences$first_difference,
    influence_first_differences$student_intervals
  )
)

knitr::kable(
  transform(
    influence_pooled,
    influence_term_zero = round(influence_term_zero, 2),
    fitted_model = round(fitted_model, 2),
    first_difference = round(first_difference, 2)
  ),
  row.names = FALSE,
  col.names = c(
    "Summary",
    "Student-intervals",
    "Influence term set to zero",
    "Fitted model",
    "Average first difference"
  ),
  caption = "Average distance from the starting drinking level of current friends"
)
```

With the influence term set to zero, students finish about `r sprintf("%.2f", influence_pooled$influence_term_zero)` drinking categories from their starting friends' average. With the fitted influence term, that distance is about `r sprintf("%.2f", influence_pooled$fitted_model)` categories, a reduction of `r sprintf("%.2f", abs(influence_pooled$first_difference))` categories. That is the model's plain-language influence summary: including the fitted term reduces the simulated distance to where current friends began. The estimated influence coefficient is `r round(th2[m2$effects$effectName=="drink average alter"],2)` with SE `r round(se2[m2$effects$effectName=="drink average alter"],2)`, so this small point-estimate contrast should not be presented as clear evidence of influence.

These first differences average over the observed starting configurations and over the model's simulated paths. They hold the other fitted coefficients fixed, but they do not include parameter uncertainty and they are not holdout forecasts. For a research analysis, repeat the calculation across parameter draws or bootstrap refits and report an interval. Use the goodness-of-fit checks in [C2](#c2) to ask whether the fitted process reproduces network and behavior features beyond the statistics targeted in estimation.

::: {.callout-tip collapse="true" title="Stuck? Two quirks in these functions"}
**First**, by default `selectionTable()` returns a long-format *data frame* (`ego`, `vego`, `valter`, `select`), not the grid above, so `round(st, 3)` errors with `non-numeric-alike variable(s) in data frame: ego`. Pass **`as.matrix = TRUE`** to get the grid. `silent = TRUE` suppresses the parameter-echo it otherwise prints.

**Second**, never `head()` the default (data-frame) form and print it. `print(head(st, 3))` errors with a `dimnames` complaint, because the print method assigns row names from a fixed-length attribute regardless of how you subset. `print(st)` on the whole object is fine. That quirk is present in 1.5 and 1.6 alike, you did nothing wrong.
:::

### Exercise {#exercise}

::: {.panel-tabset}

## Core

The motivating question was about **smoking**, not drinking. Swap it in: make `s50s` the behavior dependent variable instead of `s50a`, refit, and answer the question in three sentences.

```{r}
#| label: ex-core
#| eval: false
smokeDep <- sienaDependent(array(______, dim = c(50, 1, 3)), type = "behavior")
dat3 <- sienaDataCreate(friend, smokeDep)
eff3 <- getEffects(dat3)
eff3 <- includeEffects(eff3, recip, transTrip)
eff3 <- includeEffects(eff3, egoX, altX, simX, interaction1 = "smokeDep")
eff3 <- includeEffects(eff3, name = "smokeDep", avAlt, interaction1 = "friend")
m3 <- siena07(alg, data = dat3, effects = eff3, batch = TRUE, verbose = FALSE)
```

## Stretch

Replace `avAlt` with `avSim` (average similarity) in the behavior equation. These encode different mechanisms: `avAlt` says ego moves toward the *mean level* of their friends; `avSim` says ego is rewarded for *resembling* friends, which is symmetric around each friend's value.

Fit both, and decide which the data prefer. State how you decided, and note that you cannot simply compare t-ratios on differently-parameterized terms, so say what you actually used.

## Solution

```{r}
#| label: ex-solution
smokeDep <- sienaDependent(array(s50s, dim = c(50, 1, 3)), type = "behavior")
dat3 <- sienaDataCreate(friend, smokeDep)
eff3 <- getEffects(dat3)
eff3 <- includeEffects(eff3, recip, transTrip)
eff3 <- includeEffects(eff3, egoX, altX, simX, interaction1 = "smokeDep")
eff3 <- includeEffects(eff3, name = "smokeDep", avAlt, interaction1 = "friend")

m3 <- fit_or_load("m3_smoking",
  siena07(alg, data = dat3, effects = eff3, batch = TRUE, verbose = FALSE))

th3 <- m3$theta; se3 <- sqrt(diag(m3$covtheta))
out <- data.frame(effect = m3$effects$effectName,
                  est = round(th3, 3), se = round(se3, 3), t = round(th3/se3, 2))
knitr::kable(out[grepl("similarity|average alter", out$effect), ], row.names = FALSE)
cat("tconv.max:", round(m3$tconv.max, 3), "\n")
```

```{r}
#| label: ex-solution-t
tt3 <- setNames(th3 / se3, m3$effects$effectName)
round(tt3[c("smokeDep similarity", "smokeDep average alter")], 2)
```

**Answering the question.** No, a simple comparison of wave 1 with wave 2 cannot separate the two mechanisms. Two friends who smoke alike at wave 2 could have chosen each other partly because they were already similar, or their smoking could have become more similar after the friendship formed. A joint network-behavior model gives those processes separate parameters, but the separation still depends on the model's assumptions about change between waves.

On these data neither reaches conventional significance: smoking selection has a t-ratio of about `r round(tt3[["smokeDep similarity"]],2)` and smoking influence about `r round(tt3[["smokeDep average alter"]],2)`. The three-sentence summary is that her comparison cannot work even in principle, that the model which *can* work gives suggestive but inconclusive evidence for both mechanisms at n = 50, and that the right response is more waves or more actors rather than a stronger claim from these data.

(Check `tconv.max` before you read anything above. If it's over 0.25, refit with `prevAns = m3` before interpreting.)

**Stretch solution, `avSim` instead of `avAlt`.** Same data, same network side, one line changed in the behavior equation:

```{r}
#| label: ex-solution-stretch
eff3s <- getEffects(dat3)
eff3s <- includeEffects(eff3s, recip, transTrip)
eff3s <- includeEffects(eff3s, egoX, altX, simX, interaction1 = "smokeDep")
eff3s <- includeEffects(eff3s, name = "smokeDep", avSim, interaction1 = "friend")

m3s <- fit_or_load("m3s_smoking_avsim",
  siena07(alg, data = dat3, effects = eff3s, batch = TRUE, verbose = FALSE))

th3s <- m3s$theta; se3s <- sqrt(diag(m3s$covtheta))
tt3s <- setNames(th3s / se3s, m3s$effects$effectName)
c(avAlt_influence_t = round(tt3[["smokeDep average alter"]], 2),
  avSim_influence_t = round(tt3s[["smokeDep average similarity"]], 2),
  avSim_tconv       = round(m3s$tconv.max, 3))
```

`avAlt` and `avSim` are **not nested**, so you cannot ask "which term is significant" and read off a winner, they are two different single-parameter influence specifications. The defensible ways to choose are (a) theory, do you believe actors move toward the *mean level* of their friends (`avAlt`) or toward *resembling each friend* (`avSim`, symmetric around each friend's value)?, and (b) fit, via `sienaGOF` on the `BehaviorDistribution` auxiliary for each. What you must **not** do is compare the two influence t-ratios as if the larger one wins. Here influence is `r round(tt3s[["smokeDep average similarity"]],2)` under `avSim` versus `r round(tt3[["smokeDep average alter"]],2)` under `avAlt`, both short of significance, so at n = 50 these data do not distinguish between the two mechanisms.
:::

## C2 · Goodness of fit {#c2}

Converging is not fitting. `tconv.max` says the estimator found the parameters that match the statistics you asked it to match. It says nothing about whether the model reproduces features of the data you *didn't* target.

This is the same epistemic move as `ergm::gof()` yesterday, and as the posterior predictive checks on the SRM day: **simulate from the fitted model, compare to what you observed on statistics you did not fit.** Third time this course has done it, different implementation each time.

`sienaGOF` needs those simulated networks, which is why we set `returnDeps = TRUE`.

```{r}
#| label: c2-gof
g_indeg <- sienaGOF(m2, IndegreeDistribution, varName = "friend",
                    verbose = FALSE, join = TRUE)
g_behav <- sienaGOF(m2, BehaviorDistribution, varName = "drink",
                    verbose = FALSE, join = TRUE)
g_triad <- sienaGOF(m2, TriadCensus, varName = "friend",
                    verbose = FALSE, join = TRUE)

checkpoint(`indegree p`   = g_indeg[[1]]$p,
           `behavior p`   = g_behav[[1]]$p,
           `triadcensus p`= g_triad[[1]]$p, label = "GOF")
```



Two checks are compatible with the observed data and one is not. Read each result as a diagnostic rather than treating the p-value as a grade:

| Check | Result | What a miss would lead us to examine |
|---|---:|---|
| In-degree distribution | $p=$ `r round(g_indeg[[1]]$p,3)` | Unmodeled differences in how popular students are |
| Behavior distribution | $p=$ `r round(g_behav[[1]]$p,3)` | A poor description of how drinking changes |
| Triad census | **$p=$ `r round(g_triad[[1]]$p,3)`** | Missing closure, cycles, or local subgroup structure |

One isolated miss points toward a targeted specification question. Broad misses across several different checks would raise a more basic concern about whether the micro-step account is a useful description of the process. Here the degree and behavior checks are compatible with the observed data, while the triad census is not, so the immediate problem is local three-person structure.

A low p-value here means the observed statistic is extreme relative to what the model simulates. So our model gets degree and behavior right and local triadic structure wrong, despite having `transTrip` and `cycle3` in it. (These three are not the whole GOF menu - `OutdegreeDistribution` and the numbers-behind-the-plot function `descriptives.sienaGOF()` are in [the C3 toolbox](#toolbox-c), along with the pointer for writing your own auxiliary.)

```{r}
#| label: c2-plot
#| fig-height: 4.2
plot(g_indeg)
```

The plot shows the simulated distribution as a band with the observed values overlaid. When the observed line sits inside the band, that statistic fits.

**Now act on it.** A GOF failure you don't respond to is decoration. The triad census misfit says we're missing something about local closure, so let's add degree-based popularity, a common source of triadic structure that pure closure terms don't capture.

```{r}
#| label: c2-fix
eff2b <- includeEffects(eff2, inPopSqrt)
m2b <- fit_or_load("m2b_inpop",
  siena07(alg, data = dat2, effects = eff2b,
          batch = TRUE, verbose = FALSE, returnDeps = TRUE, prevAns = m2))
g_triad_b <- sienaGOF(m2b, TriadCensus, varName = "friend",
                      verbose = FALSE, join = TRUE)
checkpoint(`tconv.max` = m2b$tconv.max,
           `triad p, before` = g_triad[[1]]$p,
           `triad p, after`  = g_triad_b[[1]]$p, label = "AFTER ADDING inPopSqrt")
```



**It didn't work, and that is the lesson.** The triad census p-value went from `r round(g_triad[[1]]$p,3)` to `r round(g_triad_b[[1]]$p,3)`, no better, arguably worse - *and* `tconv.max` came out at `r round(m2b$tconv.max,3)`, which is over the 0.25 threshold. So the model with the extra effect both fails GOF on the same statistic and is less well converged than the one we started with.

Two things follow, and they're the reason I didn't script this to succeed.

First, **check convergence before you read any GOF number.** A GOF p-value from a model that hasn't converged is not telling you about the model; it's telling you about wherever the estimator happened to stop. The correct next move here is `prevAns` again until `tconv.max` is under 0.25, *then* re-read the GOF.

Second, this unconverged run cannot tell us whether the proposed mechanism helped. Degree-based popularity was a reasonable hypothesis about what generates the extra triadic structure, but we must first refit to convergence and then rerun GOF. Only if the converged model still misses the triad census would this specification count against that hypothesis. Candidates we haven't tried include `transTies` (closure counted per alter rather than per triplet), `gwespFF` (geometrically weighted closure, which downweights repeated closure the way `gwesp` did yesterday), or the possibility that the misfit is driven by a handful of dense subgroups no dyadic effect will capture.

What you should *not* do is add effects until the p-value crosses 0.05 and then report that specification. That's the same pathology as stepwise regression, and with a triad census of 16 correlated statistics it will find something eventually.

::: {.callout-note collapse="true" title="Depth: What `sienaGOF()` is actually doing"}
It is a Monte Carlo test using the Mahalanobis distance between the observed auxiliary-statistic vector and the distribution of simulated ones. That calculation accounts for the fact that components such as counts of different triad types are heavily correlated. `join = TRUE` pools across waves rather than testing each period separately.

The reference is **Lospinoso & Snijders (2019), "Goodness of fit for stochastic actor-oriented models," *Methodological Innovations* 12(3)**. Worth reading if you plan to report GOF, because it's specific about which auxiliary statistics are informative, and about the fact that with enough simulations you can reject any model, so the p-value is a guide rather than a verdict.
:::

## C3 · Testing without fitting, and the SAOM–TERGM dispute {#c3}

### The score test

A score test asks a narrow question: **if we keep a proposed new effect fixed at zero, does the current model show evidence that its coefficient wants to move away from zero?** This lets us screen an omitted effect without fully estimating the larger model.

Here the proposed addition is `inPopSqrt`, an indegree-popularity effect. It asks whether students who already receive many friendship nominations are especially likely to receive additional nominations. We add the effect to the calculations but do not let RSiena estimate its coefficient:

```{r}
#| label: c3-score
eff_sc <- setEffect(eff2, inPopSqrt, fix = TRUE, test = TRUE, initialValue = 0)
m_sc <- fit_or_load("m2_scoretest",
  siena07(alg, data = dat2, effects = eff_sc,
          batch = TRUE, verbose = FALSE, prevAns = m2))
m_sc
```



The three arguments do different jobs:

| Argument | What it tells RSiena to do |
|---|---|
| `fix = TRUE` | Keep the popularity coefficient fixed rather than estimate it |
| `initialValue = 0` | Fix it at zero |
| `test = TRUE` | Test whether the restriction at zero is consistent with the observed moment |

Look at the `inPopSqrt` row in the output. Its estimate is exactly 0.0000 and its standard error is `NA` because we forced the coefficient to stay at zero. Those values are not a fitted null effect and a failed uncertainty calculation. The coefficient was never freely estimated.

### What the score measures

At the restricted value of zero, RSiena compares the popularity-related change observed in the data with the corresponding change in networks simulated from the current model. If the simulations systematically contain too little or too much popularity concentration, there is pressure for the coefficient to move away from zero. The score summarizes that local pressure.

A regression analogy helps. Suppose we have fitted $Y=\beta_0+\beta_1X_1$ and wonder whether $X_2$ belongs in the model. A score test does not travel all the way to the best-fitting value of $\beta_2$. It stands at $\beta_2=0$ and checks the local slope: flat ground provides little evidence for moving, while a steep slope suggests that zero is a poor restriction.

This is also why the score test is usually cheaper than a full refit. The larger model would repeatedly update the new coefficient jointly with every existing coefficient. The score test leaves the proposed coefficient at zero and evaluates the evidence at that one location.

That number is buried in console output, and papers report it. `score.Test()` extracts it as an object:

```{r}
#| label: c3-scoretest-fn
st_inpop <- score.Test(m_sc)
st_inpop
```

Same chi-square, same p, but now `st_inpop$pvalue` is `r round(st_inpop$pvalue, 3)`, a value you can put in a table or an inline sentence instead of retyping it from the console. The result is $\chi^2=`r round(st_inpop$chisquare,2)`$ with $p=`r round(st_inpop$pvalue,3)`. This is evidence against the restriction $\beta_{\text{popularity}}=0$. In practical terms, the observed friendship changes contain more popularity-related structure than the restricted model reproduces, so a freely estimated popularity effect is worth investigating.

The one-sided Z is `r round(st_inpop$onesided,2)`, but its sign follows RSiena's score and parameterization conventions. Do not use that sign to describe the direction of the substantive effect. To learn direction and magnitude, estimate the expanded model.

### What the test does not tell us

The score test does not estimate the popularity coefficient, its standard error, or its effect on predictions. It also cannot tell us how outdegree, reciprocity, transitivity, or drinking selection will change when popularity is freely estimated. Some apparent closure, for example, may be reallocated to popular students attracting many nominations once both processes compete in the fitted model.

It also does not show that the expanded model will converge or improve goodness-of-fit. We already saw that distinction in [C2](#c2): the freely estimated popularity model did not converge adequately and did not improve the triad-census check. The score test correctly detected local evidence against zero, but that did not guarantee a successful larger model.

The proper sequence is therefore:

1. Use the score test to screen a theoretically motivated omitted effect.
2. If zero looks inadequate, add the effect as a freely estimated parameter.
3. Refit the expanded model.
4. Check convergence before interpreting it.
5. Examine the new estimate, changes in the other coefficients, and goodness-of-fit.

::: {.callout-tip collapse="true" title="Stuck? `includeEffects()` will not fix a parameter"}
If you try `includeEffects(eff, inPopSqrt, initialValue = 0)` RSiena warns `argument 'initialValue' has no effect in includeEffects; use setEffect`. The warning is easy to miss and you'll end up estimating the effect rather than testing it. The correct idiom is the `setEffect(..., fix = TRUE, test = TRUE, initialValue = 0)` above.
:::

Why bother, when we just fit the bigger model in C2? With several theoretically plausible candidate effects, score tests can help prioritize which expanded models deserve the cost of a full fit. They are diagnostic screens, not substitutes for those fits. A nonsignificant score test is also not proof that an effect is absent, especially when the data contain little information about it.

### Time heterogeneity

We've been assuming one set of parameters governs both periods. That's testable.

```{r}
#| label: c3-timetest
tt <- sienaTimeTest(m2)
summary(tt)
```



The joint test at the top is the one to read: it asks whether adding period dummies to every effect improves fit. If it's non-significant, a single parameter set across periods is defensible and you can report the model you have. If it's significant, your parameters are changing between waves and pooling them is hiding that, you'd add period interactions for the offending effects with `includeTimeDummy()`, which [the toolbox below](#toolbox-c) demonstrates so you have the idiom when your own data needs it.

::: {.callout-note collapse="true" title="The full toolbox: Joint tests, the time-dummy remedy, more GOF, and the table your paper needs (optional self-study)" #toolbox-c}
Four things a finished SAOM paper needs that the core path only pointed at. All runnable; the one refit is cached, so it loads instantly.

**Joint tests.** The score test above tests one effect. Referees routinely ask for joint tests, "are the structural terms needed *as a block*?", and `Multipar.RSiena()` is a joint Wald test on any set of parameters, addressed by position (read the positions off `m2$effects$effectName`):

```{r}
#| label: tb-c-joint
Multipar.RSiena(m2, 4, 5, 6)   # H0: recip = transTrip = cycle3 = 0
testSame.RSiena(m2, 7, 8)      # H0: drink alter effect == drink ego effect
```

`testSame.RSiena()` asks a different question, not "are these zero" but "are these *equal*", and here it cannot distinguish the alter effect of drinking from the ego effect, which is worth a sentence in a paper all by itself. (For arbitrary linear combinations there's `Wald.RSiena(A, m2)` with your own contrast matrix `A`; `Multipar` is the convenience wrapper you'll actually use.)

**The time-dummy remedy.** Our joint `sienaTimeTest` came back non-significant (p = `r round(tt$JointTest[1, 1], 3)`), so *this* model needs no repair, but when yours comes back significant, the remedy is `includeTimeDummy()`, and you should have run it once before you need it. Here it is on reciprocity, period 2:

```{r}
#| label: tb-c-timedummy
# capture.output() only mutes includeTimeDummy's verbose row-echo
invisible(capture.output(
  eff2t <- includeTimeDummy(eff2, recip, timeDummy = "2")))
m2t <- fit_or_load("m2_timedummy",
  siena07(alg, data = dat2, effects = eff2t,
          batch = TRUE, verbose = FALSE, prevAns = m2))
i_d <- grep("Dummy2.*reciprocity", m2t$effects$effectName)
checkpoint(`dummy est` = m2t$theta[i_d],
           `dummy t`   = m2t$theta[i_d] / sqrt(diag(m2t$covtheta))[i_d],
           `tconv.max` = m2t$tconv.max, label = "PERIOD-2 RECIPROCITY DUMMY")
```

The dummy is the *period-2 shift* in reciprocity, and its t-ratio (`r round(m2t$theta[i_d] / sqrt(diag(m2t$covtheta))[i_d], 2)`) agrees with the joint test's verdict: no time heterogeneity worth modeling here. On your own data the loop is: `sienaTimeTest` → dummy the flagged effects → refit → re-check. One 1.5.0 quirk found while writing this: calling `sienaTimeTest()` *again* on this dummy-containing fit dies with a `dims` error on this specification, so judge the remedy by the dummy's own t-ratio, which is what it's for.

**More GOF, and the numbers behind the band.** [C2](#c2) used three auxiliaries; the fourth stock one is `OutdegreeDistribution`, and `descriptives.sienaGOF()` prints the simulated envelope the GOF plot draws:

```{r}
#| label: tb-c-gofplus
g_out <- sienaGOF(m2, OutdegreeDistribution, varName = "friend",
                  verbose = FALSE, join = TRUE)
round(descriptives.sienaGOF(g_out), 2)
```

Outdegree fits too (p = `r round(g_out[[1]]$p, 3)`), unsurprising, since outdegree is close to what the rate and density parameters target. The `obs` row against the `perc.lower`/`perc.upper` rows is the plot's band in numeric form, which is what you report when a referee asks for more than a picture. For auxiliaries beyond the stock four, geodesic distributions, mixed triads, see `?"sienaGOF-auxiliary"`, which documents how to write your own.

**The table your paper needs.** [C1](#c1) built its table by hand so you'd know where `$theta` and `$covtheta` live. For the manuscript, one line:

```{r}
#| label: tb-c-table
siena.table(m2, type = "html", sig = TRUE, d = 3)
```

That writes **`m2.html`**, estimates, standard errors, significance daggers, network and behavior equations grouped, next to this document; `type = "tex"` gives you the LaTeX version. It is deliberately plain, which is correct for a first draft you'll edit anyway.
:::

### The dispute you should know about

You are in an unusually good seat for this one, because as of yesterday you have fit both contenders on the very same data. Day 11's C4 put a TERGM on this s50 panel: each wave conditioned on the last, ties updating with no actor anywhere in the machinery, estimated by pooled MPLE with bootstrapped intervals. Today you put a SAOM on it: one actor, one micro-step at a time, in the unobserved interval between photographs. The model you fit yesterday assumes the waves *are* the process; the model you fit today assumes the waves are censored glimpses of a continuous one. Same panel, two stories about what happened between snapshots, so it should not surprise you that people fought about it in print.

Here is the single best story in this literature, and you'll meet it if you ever submit a paper using either model. In 2019 a paper claimed TERGMs beat SAOMs at out-of-sample prediction. In 2022 someone re-read the code and found that **3 of the model's 13** predictors had been computed from the very networks it was trying to predict, it was forecasting the future using the future. The original authors published a formal correction. The methodological question underneath it all, though, is still open. Now the same three movements with the citations attached:

**Leifeld & Cranmer (2019, *Network Science* 7(1):20–51)** compared TERGM and SAOM theoretically and by out-of-sample prediction. Their own conclusion was not triumphal:

> *"Though we do not aim to make a general claim about either being superior to the other across all specifications... Each model out-predicts the other one the more the specific assumptions of the respective model are met."*

**Block, Hollway, Stadtfeld, Koskinen & Snijders (2022, *Network Science* 10(1))** are the ones who re-read the code: the TERGM specification used nodal attributes computed from the *outcome* networks' observed degrees, so 3 of its 13 terms carried information from the held-out wave. The out-of-sample prediction was partly predicting the future using the future. They also argued tie-level predictive accuracy is the wrong criterion, since these models are about endogenous processes rather than tie location.

**Leifeld & Cranmer (2022, *Network Science* 10(1):15–19)** conceded the coding error, the formal acknowledgement is their corrigendum, *Network Science* 10(1):111, and pressed a deeper point that does not go away: **the SAOM is itself a theory**, not a neutral estimator. It hard-codes sequential one-tie-at-a-time micro-steps, actors with full network knowledge, and one-period look-ahead. Embedded theory ought to be falsifiable, and they argue held-out prediction is the only feasible way to test it.

**Where that leaves you: the factual dispute is settled, the methodological one is not.** The coding error occurred and was acknowledged. Whether predictive fit is a legitimate criterion for comparing generative network models remains an open methodological debate.

The useful framing is Block, Stadtfeld & Snijders (2019) again: ERGM is defined at the **graph** level, SAOM at the **transition** level. SAOM's statistics nest within actors, so its dependence assumptions are *stronger*. That structure provides a route to identification, but only if the micro-step and specification assumptions are defensible. You should know what you paid.

```{r}
#| label: offroad-2
# Optional extension space for a co-evolution question.
# Keep exploratory terms here rather than changing the core specification.

```

```{r}
#| label: c3-checkpoint
checkpoint(`selection t` = (th2/se2)[m2$effects$effectName == "drink similarity"],
           `influence t` = (th2/se2)[m2$effects$effectName == "drink average alter"],
           `triad GOF p` = g_triad[[1]]$p, label = "END BLOCK C")
```

---

# Wrap {#wrap}

Three days, one sentence: **absorb it (Day 10) → test it (Day 11) → parameterize temporal mechanisms with time (Day 12, today).**

Day 10's latent factor models represent residual dependence so we can see how covariate estimates and predictions change after adding latent structure. They do not automatically make the coefficients trustworthy: if an omitted factor is correlated with an observed predictor, latent adjustment cannot recover a separate direct effect without additional assumptions. Day 11's ERGMs treat dependence as a hypothesis to parameterize and test. Today, SAOMs use repeated observation to go after something neither can get from one cross-section, the relative contribution of two mechanisms that are observationally identical in a single wave, and pay for it with a strong, untestable assumption about how change happens (one tie or one behaviour step at a time).

## After the Model Runs: Describe the Change Process {#sec-after-fit}

A SAOM result should sound like an account of the fitted change process between observed waves, not like a cross-sectional regression with unfamiliar coefficients and not like a record of choices we directly observed. Work through the fitted model in this order:

1. **Separate opportunities from evaluations.** Rate parameters describe how many opportunities for change the model assigns to actors between waves. Evaluation parameters describe which candidate microsteps are relatively attractive when an opportunity occurs.
2. **Make the comparison explicit.** An evaluation coefficient compares two possible next network or behavior states that differ in a statistic, holding the other evaluation terms fixed. Exponentiating a one-unit difference gives a relative weight for those candidate microsteps, not an unconditional tie probability.
3. **Keep selection and influence in their own equations.** A network similarity effect is evidence consistent with selection into ties with similar actors. A behavior average-alter effect is evidence consistent with influence toward alters. Neither label alone establishes a causal mechanism.
4. **Earn the story with diagnostics.** Report convergence for all effects, the overall maximum convergence ratio, repeated-run stability, and goodness-of-fit for the network and behavior features relevant to the claim.

::: {.callout-tip title="What the fitted s50 result means"}
"The panel follows 50 Scottish adolescents across three friendship and drinking observations. In the fitted change process, candidate friendship states with a returned nomination or a closed friendship path receive substantially greater weight than otherwise similar states. The drinking-similarity estimate in the friendship equation is 1.44 with a standard error of 0.61, which is consistent with students sorting into friendships with peers who report similar drinking. The peer-drinking influence estimate is 1.36 with a standard error of 0.95, so these data do not distinguish that influence pattern from zero. This is evidence about the process represented by the model, not proof that students consciously optimized these features or that friends caused one another's drinking."
:::

In another application, replace students, friendship nominations, and drinking with the actual actors, tie, behavior, and observation intervals. State which process is consistent with the estimates, which result is uncertain, whether simulations reproduce the observed changes, and which unobserved traits could still generate the same pattern.

### The referee checklist {#referee}

Realistically, most of you will read more SAOMs than you fit. Every diagnostic today doubles as a question to ask of a paper:

1. Did they report the **Jaccard index**, or otherwise establish the model applies? ([A2](#a2))
2. Did they report the **overall maximum convergence ratio**, not just individual t-ratios, and is it under 0.25? ([B2](#b2))
3. Is the **influence standard error** larger than the estimate? ([C1](#c1))
4. Did they run **GOF**, and on which auxiliary statistics? Did any fail, and did they act on it? ([C2](#c2))
5. Do they interpret coefficients as probabilities or marginal effects? (They shouldn't - [B3](#b3))
6. Do they claim selection *versus* influence from a null result on one of them? ([C1](#c1))

I have seen every one of those failures in published work.

### What we did not show {#not-shown}

RSiena is a large package. The main walkthrough and toolboxes cover the working core; the items below are useful next steps when a research question requires them:

- **`varDyadCovar()`**, dyadic covariates that change between waves (trade flows per period, shifting alliances). Same idea as `coDyadCovar()` in [the B1 toolbox](#toolbox-b), one array dimension bigger. `?varDyadCovar`.
- **`sienaDataConstraint()`**, structural zeros and ones: forbidding or freezing specific ties (siblings can't be "not friends", the teacher can't be nominated). `?sienaDataConstraint`.
- **Multi-group designs**, `sienaGroupCreate()` stacks many small networks (classrooms, villages) into one estimation; `siena08()` meta-analyzes separately fitted groups and `funnelPlot()` displays them. The standard design when no single network is big enough. `?sienaGroupCreate`, `?siena08`.
- **Generalized method of moments**, `includeGMoMStatistics()` adds extra target statistics beyond one-per-parameter, for efficiency. Current research frontier, not standard practice yet. `?includeGMoMStatistics`.
- **Two-mode (bipartite) networks**, SAOMs for actors choosing among clubs, boards, or events, via `sienaNodeSet()` and a `bipartite` dependent variable. `?sienaNodeSet`.
- **`updateTheta()` / `updateSpecification()`**, carrying estimates from one fit into a *differently specified* effects object, for when [B2](#b2)'s `prevAns` doesn't apply because the specification changed too much. `?updateTheta`.
- **Maximum-likelihood estimation**, beyond the comparison in the [method-of-moments detail box](#sec-estimation-roadmap): `sienaAlgorithmCreate(maxlike = TRUE)` is slower but more efficient with substantial missingness or very small networks. The manual's estimation chapter covers when it earns its runtime.

### Homework 4

Homework 4 offers two choices. You can interpret a supplied longitudinal-network result, focusing on what the estimates say about friendship change, selection, influence, uncertainty, and the limits of the evidence. Alternatively, you can share a piece of research connected to the course and meet with me to work through its research question, network measurement, model choice, or interpretation. Neither option requires fitting a new model.

### Where we go next: Day 13

Today closes the modeling arc of the course: every day so far has answered dependence by changing the *model*. Tomorrow (Day 13) leaves the model alone and interrogates the *design*: what a causal claim on a network even means, why interference breaks the usual potential-outcomes machinery, what randomization buys you, and which claims an observational design can support. Then Day 14 puts every model from this week head-to-head on one dataset and asks the question the syllabus opened with: which model, when?

---

# Appendix A: Depth track {#appendix-a}

::: {.callout-note collapse="true" title="A.1: Why the behavior equation needs shape terms"}
The linear and quadratic shape terms are included by default because a bounded ordinal behavior needs its own baseline change pattern. With only a linear term, the baseline evaluation can tilt toward higher or lower values but cannot represent a preference for the middle rather than the extremes, or the reverse. Omitting that curvature can force the friends' drinking term to absorb a baseline pattern that has nothing to do with friends.

The quadratic shape also has a substantive reading: negative means the behavior is self-correcting toward the middle, positive means it's self-reinforcing toward the extremes (polarizing). That's often more interesting than the influence term and it usually goes unreported.
:::

::: {.callout-note collapse="true" title="A.2: What to read next"}
- **Snijders (2017)**, *Annu Rev Stat Appl* 4:343–363, the overview. Start here.
- **Steglich, Snijders & Pearson (2010)**, *Sociological Methodology* 40(1):329–393 - the assigned reading; the selection/influence argument in full.
- **Lospinoso & Snijders (2019)**, *Methodological Innovations* 12(3), GOF.
- **Block, Stadtfeld & Snijders (2019)**, *SMR* 48(1):202–239, the graph-level vs transition-level framing.
- **Ripley, Snijders, Boda & Vörös**, *Manual for RSiena*, current edition, the reference manual; note the 1.6 function names.
- The **ancillary TERGM deep dive** (`ancillary/tergm.qmd`), the long-form companion to the TERGM you fit yesterday in Day 11's C4, fits the Knecht data both ways, so you can see the comparison in this section made empirical.
:::

---

# Appendix B: Session info {#appendix-b}

```{r}
#| label: sessioninfo
sessionInfo()
```
