Day 12: Stochastic Actor-Oriented Models

ICPSR Network Analysis: Advanced Topics · Day 12

Author

Shahryar Minhas

Published

July 30, 2026

Before we start

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.

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:

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

Already comfortable with parts of this? Where to jump

You’ve fit SAOMs before and want the diagnostics and limitations → Block C
You want the convergence machinery specifically → B2
You want to know how to referee a SAOM paper without fitting one → the referee checklist
You want the mathematics the core path skips → Appendix A

Setup check: Run this the night before

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

require_pkgs(c("RSiena", "netify"))
suppressMessages({library(RSiena)})
stopifnot(exists("s501"), exists("s50a"), exists("s50s"))
cat("RSiena  ", as.character(packageVersion("RSiena")), "\n")
#> RSiena   1.5.0
cat("R       ", as.character(getRversion()), "\n")
#> R        4.3.3
cat("SETUP OK\n")
#> SETUP OK

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, reference on demand, not something to memorise.

Block A: Is this even a SAOM problem?

A1 · The motivating question

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’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?

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:

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

The Applied 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

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); you never have to memorise it. With that parked, the motivating question runs straight into the pre-flight checks and the micro-step.

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.

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

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.

# 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)
#> [1] 50 50
#> [1] 50 50
#> [1] 50 50
table(s50a, useNA = "ifany")   # Alcohol, 1-5
#> s50a
#>  1  2  3  4  5 
#>  9 41 41 42 17
table(s50s, useNA = "ifany")   # Smoking, 1-3
#> s50s
#>   1   2   3 
#> 102  15  33
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.

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)
#> ------------------------------------------------------------------
#> CHECKPOINT: jaccard 1->2 = 0.331   |   jaccard 2->3 = 0.384   |   avg degree = 2.26, 2.32, 2.44
#> ------------------------------------------------------------------

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 0.331 and 0.384. 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 (2.26 → 2.32 → 2.44), 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.

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.

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:

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.

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:

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.

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, so this aside builds a minimal one inline:

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)])
#> Tie changes between subsequent observations:
#>  periods        0 =>  0   0 =>  1   1 =>  0   1 =>  1   Distance Jaccard   Missing
#>   1 ==>   2      2278        59        56        57       115     0.331         0 (0%)
#>   2 ==>   3      2278        56        50        66       106     0.384         0 (0%)

That writes s50_prefit_report.txt into your working directory, 133 lines of densities, degrees, missingness, and the change table shown above. Its Jaccard column is exactly the 0.331 and 0.384 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

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.

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 with an actual number.

How SAOM Estimation Works

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.

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:

Starts with parameter guesses
Simulates network evolution
Calculates the simulated target statistics
Compares them with the observed targets
Adjusts the parameters in a direction that reduces the joint mismatch
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.

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

checkpoint(`waves` = 3, `actors` = nrow(s501),
           `ties w1/w2/w3` = c(sum(s501), sum(s502), sum(s503)),
           label = "END BLOCK A")
#> ------------------------------------------------------------------
#> END BLOCK A: waves = 3   |   actors = 50   |   ties w1/w2/w3 = 113, 116, 122
#> ------------------------------------------------------------------

Block B: Fit it, converge it, read it

B1 · The first fitted model

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

friend  <- sienaDependent(s50_friend_array)
smoke1  <- coCovar(s50s[, 1])       # Smoking at wave 1, treated as fixed
dat1    <- sienaDataCreate(friend, smoke1)
dat1
#> Dependent variables:  friend 
#> Number of observations: 3 
#> 
#> Nodeset                  Actors 
#> Number of nodes              50 
#> 
#> Dependent variable friend          
#> Type               oneMode         
#> Observations       3               
#> Nodeset            Actors          
#> Densities          0.046 0.047 0.05
#> 
#> Constant covariates:  smoke1

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.

eff1 <- getEffects(dat1)
eff1 <- includeEffects(eff1, recip, transTrip, cycle3)
#>   effectNumber effectName          shortName include fix   test  initialValue parm
#> 1 16           reciprocity         recip     TRUE    FALSE FALSE          0   0   
#> 2 22           transitive triplets transTrip TRUE    FALSE FALSE          0   0   
#> 3 42           3-cycles            cycle3    TRUE    FALSE FALSE          0   0
eff1 <- includeEffects(eff1, sameX, interaction1 = "smoke1")
#>   effectNumber effectName  shortName include fix   test  initialValue parm
#> 1 310          same smoke1 sameX     TRUE    FALSE FALSE          0   0

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

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.

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:

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 (455 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.

alg <- sienaAlgorithmCreate(projname = NULL, seed = 6886, n3 = 1000)
#> If you use this algorithm object, siena07 will create/use an output file Siena.txt .
m1 <- fit_or_load("m1_network_only",
  siena07(alg, data = dat1, effects = eff1,
          batch = TRUE, verbose = FALSE, returnDeps = TRUE))

m1
#> Estimates, standard errors and convergence t-ratios
#> 
#>                                    Estimate   Standard   Convergence 
#>                                                 Error      t-ratio   
#> 
#> Rate parameters: 
#>   0.1      Rate parameter period 1  6.5170  ( 1.0881   )             
#>   0.2      Rate parameter period 2  5.2088  ( 0.8519   )             
#> 
#> Other parameters: 
#>   1.  eval outdegree (density)     -2.7780  ( 0.1579   )   -0.0065   
#>   2.  eval reciprocity              2.4753  ( 0.2161   )   -0.0088   
#>   3.  eval transitive triplets      0.6627  ( 0.1512   )   -0.0073   
#>   4.  eval 3-cycles                -0.0937  ( 0.2916   )   -0.0098   
#>   5.  eval same smoke1              0.1401  ( 0.1482   )    0.0233   
#> 
#> Overall maximum convergence ratio:    0.0630 
#> 
#> 
#> Total of 1862 iteration steps.

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 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); this one is the raw console form on purpose, so you have seen it once.

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

The full toolbox: Covariates your own data will have, interactions, and parallel fits (optional self-study)

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.

drink_vc <- varCovar(s50a[, 1:2])   # Waves 1 and 2 only: 2 periods, 2 columns
dim(drink_vc)
#> [1] 50  2

(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:

drinkdist <- coDyadCovar(abs(outer(s50a[, 1], s50a[, 1], "-")))
datB <- sienaDataCreate(friend, drink_vc, drinkdist)
effB <- getEffects(datB)
effB <- includeEffects(effB, recip, transTrip)
#>   effectNumber effectName          shortName include fix   test  initialValue parm
#> 1 16           reciprocity         recip     TRUE    FALSE FALSE          0   0   
#> 2 22           transitive triplets transTrip TRUE    FALSE FALSE          0   0
effB <- includeEffects(effB, egoX, altX, simX, interaction1 = "drink_vc")
#>   effectNumber effectName          shortName include fix   test  initialValue parm
#> 1 268          drink_vc alter      altX      TRUE    FALSE FALSE          0   0   
#> 2 283          drink_vc ego        egoX      TRUE    FALSE FALSE          0   0   
#> 3 329          drink_vc similarity simX      TRUE    FALSE FALSE          0   0
effB <- includeEffects(effB, X, interaction1 = "drinkdist")
#>   effectNumber effectName shortName include fix   test  initialValue parm
#> 1 239          drinkdist  X         TRUE    FALSE FALSE          0   0

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

A genuine interaction: includeInteraction(). B1’s stuck-box 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?”:

effB <- includeInteraction(effB, egoX, recip, interaction1 = c("drink_vc", ""))
#>   effectNumber effectName                 shortName include fix   test  initialValue parm
#> 1  16          reciprocity                recip     TRUE    FALSE FALSE          0   0   
#> 2 283          drink_vc ego               egoX      TRUE    FALSE FALSE          0   0   
#> 3 455          drink_vc ego x reciprocity unspInt   TRUE    FALSE FALSE          0   0   
#>   effect1 effect2
#> 1   0      0     
#> 2   0      0     
#> 3 283     16

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:

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

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.

cat("overall maximum convergence ratio (tconv.max):", round(m1$tconv.max, 4), "\n")
#> overall maximum convergence ratio (tconv.max): 0.063
cat("largest individual t-ratio            :", round(max(abs(m1$tconv)), 4), "\n")
#> largest individual t-ratio            : 0.0233

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

starved <- sienaAlgorithmCreate(projname = NULL, seed = 42, n3 = 200, nsub = 2)
#> If you use this algorithm object, siena07 will create/use an output file Siena.txt .
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")
#> ------------------------------------------------------------------
#> STARVED RUN: tconv.max = 0.343   |   max |t| = 0.173
#> ------------------------------------------------------------------

Both diagnostics fail: tconv.max is 0.343, above the 0.25 publication threshold, and the worst individual t-ratio is 0.173, 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.

full <- sienaAlgorithmCreate(projname = NULL, seed = 42, n3 = 1000)
#> If you use this algorithm object, siena07 will create/use an output file Siena.txt .
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")
#> ------------------------------------------------------------------
#> REPAIR ROUND 1: tconv.max = 0.236   |   max |t| = 0.064
#> ------------------------------------------------------------------

One round fixes it: tconv.max 0.236, worst individual t-ratio 0.064. 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.

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")
#> ------------------------------------------------------------------
#> REPAIR ROUND 2: tconv.max = 0.128   |   max |t| = 0.082
#> ------------------------------------------------------------------

And now it’s comfortable: 0.128. 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.”

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:

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.

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.

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, 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 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

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.

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")
#> objective change, tie to a non-nominator: -2.778
cat("objective change, tie to a nominator    :", round(f_reciprocated, 3), "\n")
#> objective change, tie to a nominator    : -0.303
cat("relative choice weight                  :", round(exp(b_recip), 2), "\n")
#> relative choice weight                  : 11.89

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 2.48. Third, exponentiate only if the relative weight helps: holding the other modeled features constant, that reciprocating move receives about 11.9 times the choice weight of an otherwise comparable move that does not create reciprocity.

Now the same logic for closure:

b_tt <- th[["transitive triplets"]]
cat("one new transitive triplet :", round(b_tt, 3), "\n")
#> one new transitive triplet : 0.663
cat("closing 3 triplets at once :", round(3 * b_tt, 3),
    "  (relative weight", round(exp(3 * b_tt), 2), ")\n")
#> closing 3 triplets at once : 1.988   (relative weight 7.3 )

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 0.66 and the move receives about 1.94 times the relative weight of an otherwise comparable move that closes none. The contribution accumulates: a tie that closes three paths adds $3(0.66)=1.99$ 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

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)
#>      one more outgoing nomination           one returned nomination 
#>                            -2.778                             2.475 
#> two friend-of-friend paths closed 
#>                             1.325
combined_score <- sum(combined_move)
cat("combined change in evaluation score:", round(combined_score, 3), "\n")
#> combined change in evaluation score: 1.023
cat("relative choice weight versus no change:", round(exp(combined_score), 2), "\n")
#> relative choice weight versus no change: 2.78

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 1.02 points in the evaluation score and receives about 2.78 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

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

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")
#> Same-smoking coefficient: 0.14
cat(
  "Conditional probability ratio:",
  round(same_smoking_probability_ratio, 2), "\n"
)
#> Conditional probability ratio: 1.15
cat(
  "Percent difference:",
  round(100 * (same_smoking_probability_ratio - 1)), "%\n"
)
#> Percent difference: 15 %

At a modeled friendship-change opportunity, a nomination to a student with the same smoking status is about 15% 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.

Optional Detail: Removing the Entire Same-Smoking Coefficient

A different question: remove the entire coefficient

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	(_{}=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 (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.

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.

eff1_smoke_on <- updateTheta(eff1, m1)
eff1_smoke_zero <- setEffect(
  eff1_smoke_on,
  sameX,
  interaction1 = "smoke1",
  initialValue = 0
)
#>   effectNumber effectName  shortName include fix   test  initialValue parm
#> 1 310          same smoke1 sameX     TRUE    FALSE FALSE          0   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.

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 (ij) exists.

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

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.

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

Model-based first difference for the same-smoking term
interval	eligible_dyads	smoke_term_zero	smoke_term_fitted	first_difference
Waves 1 to 2	1392	2.23%	2.88%	0.64 pp
Waves 2 to 3	1385	2.07%	2.59%	0.52 pp
Pooled across intervals	2777	2.15%	2.73%	0.58 pp

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 2.15% to 2.73%, an increase of 0.58 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.

More Detail: Why This Is a Model-Based First Difference, Not an Ordinary AME

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:

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?
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$?
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:

Begin from the same observed network and behavior state at wave $t$.
Simulate many paths to wave $t+1$ under the baseline values.
Make one substantively defensible change while holding the other starting information fixed.
Simulate the paths again, ideally using common random-number streams so the paired difference contains less Monte Carlo noise.
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.

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

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.

b_smoke <- th[["same smoke1"]]
exp(______)

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.

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")
#> same-smoking relative choice weight: 1.15  (t = 0.95 )

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

checkpoint(`recip` = b_recip, `transTrip` = b_tt,
           `tconv.max` = m1$tconv.max, label = "END BLOCK B")
#> ------------------------------------------------------------------
#> END BLOCK B: recip = 2.475   |   transTrip = 0.663   |   tconv.max = 0.063
#> ------------------------------------------------------------------

# 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

C1 · Selection and influence, jointly

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.

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
#> Dependent variables:  friend, drink 
#> Number of observations: 3 
#> 
#> Nodeset                  Actors 
#> Number of nodes              50 
#> 
#> Dependent variable friend          
#> Type               oneMode         
#> Observations       3               
#> Nodeset            Actors          
#> Densities          0.046 0.047 0.05
#> 
#> Dependent variable drink   
#> Type               behavior
#> Observations       3       
#> Nodeset            Actors  
#> Range              1 - 5   
#> 
#> Constant covariates:  smoke1

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:

eff2 <- getEffects(dat2)

#: Network equation: structure, then SELECTION on drinking ---
eff2 <- includeEffects(eff2, recip, transTrip, cycle3)
#>   effectNumber effectName          shortName include fix   test  initialValue parm
#> 1 17           reciprocity         recip     TRUE    FALSE FALSE          0   0   
#> 2 23           transitive triplets transTrip TRUE    FALSE FALSE          0   0   
#> 3 43           3-cycles            cycle3    TRUE    FALSE FALSE          0   0
eff2 <- includeEffects(eff2, egoX, altX, simX, interaction1 = "drink")
#>   effectNumber effectName       shortName include fix   test  initialValue parm
#> 1 424          drink alter      altX      TRUE    FALSE FALSE          0   0   
#> 2 439          drink ego        egoX      TRUE    FALSE FALSE          0   0   
#> 3 485          drink similarity simX      TRUE    FALSE FALSE          0   0

#: Behavior equation: shape, then INFLUENCE from friends ---
eff2 <- includeEffects(eff2, name = "drink", avAlt, interaction1 = "friend")
#>   effectNumber effectName          shortName include fix   test  initialValue parm
#> 1 729          drink average alter avAlt     TRUE    FALSE FALSE          0   0

The Complete Co-Evolution 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 below.)

m2 <- fit_or_load("m2_coevolution",
  siena07(alg, data = dat2, effects = eff2,
          batch = TRUE, verbose = FALSE, returnDeps = TRUE))

m2
#> Estimates, standard errors and convergence t-ratios
#> 
#>                                            Estimate   Standard   Convergence 
#>                                                         Error      t-ratio   
#> Network Dynamics 
#>    1. rate constant friend rate (period 1)  6.5274  ( 1.1674   )    0.0085   
#>    2. rate constant friend rate (period 2)  5.1070  ( 0.8356   )    0.0037   
#>    3. eval outdegree (density)             -2.7682  ( 0.1445   )    0.0290   
#>    4. eval reciprocity                      2.3807  ( 0.2177   )    0.0264   
#>    5. eval transitive triplets              0.6555  ( 0.1337   )    0.0462   
#>    6. eval 3-cycles                        -0.0759  ( 0.2647   )    0.0417   
#>    7. eval drink alter                     -0.0617  ( 0.1143   )    0.0239   
#>    8. eval drink ego                        0.0545  ( 0.1128   )    0.0311   
#>    9. eval drink similarity                 1.4396  ( 0.6129   )    0.0718   
#> 
#> Behavior Dynamics
#>   10. rate rate drink (period 1)            1.3249  ( 0.3838   )   -0.0298   
#>   11. rate rate drink (period 2)            1.7988  ( 0.5323   )    0.0508   
#>   12. eval drink linear shape               0.4190  ( 0.2369   )   -0.0163   
#>   13. eval drink quadratic shape           -0.6158  ( 0.3892   )   -0.0222   
#>   14. eval drink average alter              1.3621  ( 0.9458   )    0.0472   
#> 
#> Overall maximum convergence ratio:    0.1631 
#> 
#> 
#> Total of 3189 iteration steps.

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

the two rows this whole session is about, plus structure for context
effect	est	se	t
reciprocity	2.3807	0.2177	10.94
transitive triplets	0.6555	0.1337	4.90
drink similarity	1.4396	0.6129	2.35
drink average alter	1.3621	0.9458	1.44

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

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?	1.44	0.61	A positive selection-consistent pattern
Drinking	Does a student’s drinking move toward the average among current friends?	1.36	0.95	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

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.

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

Start by naming the coefficient and the outcome. The estimate of 1.44 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.

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

Difference in reported drinking categories	Drinking-similarity score
0	1.00
1	0.75
2	0.50
3	0.25
4	0.00

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 1.44, so moving one category closer changes the evaluation score of that possible nomination by

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"
)
#> Change in similarity score: 0.25
cat(
  "Change in nomination's evaluation score:",
  round(score_change_one_category, 2), "\n"
)
#> Change in nomination's evaluation score: 0.36
cat(
  "Relative choice weight:",
  round(relative_choice_weight_one_category, 2), "\n"
)
#> Relative choice weight: 1.43

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 43% 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:

st <- selectionTable(m2, xd = dat2, name = "friend", vname = "drink",
                     levls = 1:5, as.matrix = TRUE, silent = TRUE)
round(st, 3)
#>        1      2      3      4      5
#> 1  0.484  0.062 -0.359 -0.781 -1.202
#> 2  0.179  0.477  0.055 -0.366 -0.788
#> 3 -0.127  0.171  0.470  0.048 -0.374
#> 4 -0.432 -0.134  0.164  0.462  0.041
#> 5 -0.738 -0.439 -0.141  0.157  0.455

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.

it <- influenceTable(m2, xd = dat2, netname = "friend", behname = "drink",
                     levls = 1:5, as.matrix = TRUE, silent = TRUE)
round(it, 3)
#>        1      2      3      4      5
#> 1  2.448  1.975  0.271 -2.665 -6.832
#> 2 -0.431  0.459  0.116 -1.457 -4.262
#> 3 -3.309 -1.058 -0.038 -0.249 -1.693
#> 4 -6.188 -2.574 -0.192  0.958  0.877
#> 5 -9.067 -4.091 -0.347  2.166  3.447

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

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 (_{}=0)	Unchanged
Influence term removed	Unchanged	Same as the complete model except (_{}=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.

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.

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
)
#>   effectNumber effectName       shortName include fix   test  initialValue parm
#> 1 485          drink similarity simX      TRUE    FALSE FALSE          0   0
eff2_no_influence <- setEffect(
  eff2_fitted,
  name = "drink",
  avAlt,
  interaction1 = "friend",
  initialValue = 0
)
#>   effectNumber effectName          shortName include fix   test  initialValue parm
#> 1 729          drink average alter avAlt     TRUE    FALSE FALSE          0   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.

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.

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.

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

What the drinking-similarity selection term changes
Starting drinking gap	Initially absent nominations	Drinking-similarity coefficient set to zero	Fitted model	Average first difference
Same reported drinking	1049	2.04%	3.42%	1.38%
One category apart	1781	2.11%	2.85%	0.75%
Two or more categories apart	1841	1.82%	1.55%	-0.27%

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 3.42% to 2.04%. 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.

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:

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

Average distance from the starting drinking level of current friends
Summary	Student-intervals	Influence term set to zero	Fitted model	Average first difference
Pooled across intervals	76	0.88	0.82	-0.06

With the influence term set to zero, students finish about 0.88 drinking categories from their starting friends’ average. With the fitted influence term, that distance is about 0.82 categories, a reduction of 0.06 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 1.36 with SE 0.95, 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 to ask whether the fitted process reproduces network and behavior features beyond the statistics targeted in estimation.

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

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.

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)

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.

smokeDep <- sienaDependent(array(s50s, dim = c(50, 1, 3)), type = "behavior")
dat3 <- sienaDataCreate(friend, smokeDep)
eff3 <- getEffects(dat3)
eff3 <- includeEffects(eff3, recip, transTrip)
#>   effectNumber effectName          shortName include fix   test  initialValue parm
#> 1 16           reciprocity         recip     TRUE    FALSE FALSE          0   0   
#> 2 22           transitive triplets transTrip TRUE    FALSE FALSE          0   0
eff3 <- includeEffects(eff3, egoX, altX, simX, interaction1 = "smokeDep")
#>   effectNumber effectName          shortName include fix   test  initialValue parm
#> 1 239          smokeDep alter      altX      TRUE    FALSE FALSE          0   0   
#> 2 254          smokeDep ego        egoX      TRUE    FALSE FALSE          0   0   
#> 3 300          smokeDep similarity simX      TRUE    FALSE FALSE          0   0
eff3 <- includeEffects(eff3, name = "smokeDep", avAlt, interaction1 = "friend")
#>   effectNumber effectName             shortName include fix   test  initialValue parm
#> 1 541          smokeDep average alter avAlt     TRUE    FALSE FALSE          0   0

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)

effect	est	se	t
smokeDep similarity	0.773	0.435	1.78
smokeDep average alter	1.895	1.350	1.40

cat("tconv.max:", round(m3$tconv.max, 3), "\n")
#> tconv.max: 0.214

tt3 <- setNames(th3 / se3, m3$effects$effectName)
round(tt3[c("smokeDep similarity", "smokeDep average alter")], 2)
#>    smokeDep similarity smokeDep average alter 
#>                   1.78                   1.40

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 1.78 and smoking influence about 1.4. 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:

eff3s <- getEffects(dat3)
eff3s <- includeEffects(eff3s, recip, transTrip)
#>   effectNumber effectName          shortName include fix   test  initialValue parm
#> 1 16           reciprocity         recip     TRUE    FALSE FALSE          0   0   
#> 2 22           transitive triplets transTrip TRUE    FALSE FALSE          0   0
eff3s <- includeEffects(eff3s, egoX, altX, simX, interaction1 = "smokeDep")
#>   effectNumber effectName          shortName include fix   test  initialValue parm
#> 1 239          smokeDep alter      altX      TRUE    FALSE FALSE          0   0   
#> 2 254          smokeDep ego        egoX      TRUE    FALSE FALSE          0   0   
#> 3 300          smokeDep similarity simX      TRUE    FALSE FALSE          0   0
eff3s <- includeEffects(eff3s, name = "smokeDep", avSim, interaction1 = "friend")
#>   effectNumber effectName                  shortName include fix   test  initialValue
#> 1 498          smokeDep average similarity avSim     TRUE    FALSE FALSE          0  
#>   parm
#> 1 0

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_influence_t avSim_influence_t       avSim_tconv 
#>             1.400             1.650             0.129

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 1.65 under avSim versus 1.4 under avAlt, both short of significance, so at n = 50 these data do not distinguish between the two mechanisms.

C2 · Goodness of fit

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.

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")
#> ------------------------------------------------------------------
#> GOF: indegree p = 0.589   |   behavior p = 0.912   |   triadcensus p = 0.007
#> ------------------------------------------------------------------

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=$ 0.589	Unmodeled differences in how popular students are
Behavior distribution	$p=$ 0.912	A poor description of how drinking changes
Triad census	$p=$ 0.007	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, along with the pointer for writing your own auxiliary.)

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.

eff2b <- includeEffects(eff2, inPopSqrt)
#>   effectNumber effectName                   shortName include fix   test  initialValue
#> 1 86           indegree - popularity (sqrt) inPopSqrt TRUE    FALSE FALSE          0  
#>   parm
#> 1 0
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")
#> ------------------------------------------------------------------
#> AFTER ADDING inPopSqrt: tconv.max = 0.266   |   triad p, before = 0.007   |   triad p, after = 0.004
#> ------------------------------------------------------------------

It didn’t work, and that is the lesson. The triad census p-value went from 0.007 to 0.004, no better, arguably worse - and tconv.max came out at 0.266, 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.

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

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:

eff_sc <- setEffect(eff2, inPopSqrt, fix = TRUE, test = TRUE, initialValue = 0)
#>   effectNumber effectName                   shortName include fix  test initialValue parm
#> 1 86           indegree - popularity (sqrt) inPopSqrt TRUE    TRUE TRUE          0   0
m_sc <- fit_or_load("m2_scoretest",
  siena07(alg, data = dat2, effects = eff_sc,
          batch = TRUE, verbose = FALSE, prevAns = m2))
m_sc
#> Estimates, standard errors and convergence t-ratios
#> 
#>                                            Estimate   Standard   Convergence 
#>                                                         Error      t-ratio   
#> Network Dynamics 
#>    1. rate constant friend rate (period 1)  6.5619  ( 1.1269   )    0.0267   
#>    2. rate constant friend rate (period 2)  5.0962  ( 0.8597   )   -0.0469   
#>    3. eval outdegree (density)             -2.7710  ( 0.1480   )    0.0127   
#>    4. eval reciprocity                      2.3966  ( 0.2073   )    0.0535   
#>    5. eval transitive triplets              0.6448  ( 0.1429   )   -0.0087   
#>    6. eval 3-cycles                        -0.0642  ( 0.2880   )   -0.0013   
#>    7. eval indegree - popularity (sqrt)     0.0000  (     NA   )    0.3189   
#>    8. eval drink alter                     -0.0569  ( 0.1187   )    0.0260   
#>    9. eval drink ego                        0.0505  ( 0.1227   )   -0.0395   
#>   10. eval drink similarity                 1.4389  ( 0.6125   )    0.0607   
#> 
#> Behavior Dynamics
#>   11. rate rate drink (period 1)            1.3262  ( 0.3472   )   -0.0112   
#>   12. rate rate drink (period 2)            1.7765  ( 0.4717   )    0.0133   
#>   13. eval drink linear shape               0.4174  ( 0.2243   )    0.0116   
#>   14. eval drink quadratic shape           -0.6017  ( 0.3455   )   -0.0492   
#>   15. eval drink average alter              1.3189  ( 0.8436   )    0.0226   
#> 
#> Overall maximum convergence ratio:    0.2199 
#> 
#> 
#> Score test for 1 parameter:
#> chi-squared = 5.27, p = 0.0217.
#> 
#> Total of 3266 iteration steps.

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:

st_inpop <- score.Test(m_sc)
st_inpop
#> Tested effects:
#>  friend: indegree - popularity (sqrt) eval 
#> chi-squared = 5.27, d.f. = 1; one-sided Z = -2.29; two-sided p = 0.022.

Same chi-square, same p, but now st_inpop$pvalue is 0.022, a value you can put in a table or an inline sentence instead of retyping it from the console. The result is $\chi^2=5.27$ with $p=0.022. 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 -2.29, 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: 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:

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

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.

tt <- sienaTimeTest(m2)
summary(tt)
#> Joint significance test of time heterogeneity:
#> chi-squared =   9.18, d.f. = 10, p= 0.5148, 
#> where H0: The following parameters are zero:
#>  ( 1) (*)Dummy2:outdegree (density)
#>  ( 2) (*)Dummy2:reciprocity
#>  ( 3) (*)Dummy2:transitive triplets
#>  ( 4) (*)Dummy2:3-cycles
#>  ( 5) (*)Dummy2:drink alter
#>  ( 6) (*)Dummy2:drink ego
#>  ( 7) (*)Dummy2:drink similarity
#>  ( 8) (*)Dummy2:drink linear shape
#>  ( 9) (*)Dummy2:drink quadratic shape
#>  (10) (*)Dummy2:drink average alter
#> 
#> Individual significance tests and one-step estimators:
#>                                 Initial Est. One Step Est. p-Value
#> outdegree (density)             -2.7682      -2.5981        0.0000
#> reciprocity                      2.3807       1.8911        0.0000
#> transitive triplets              0.6555       0.5633        0.0000
#> 3-cycles                        -0.0759       0.0943        0.7740
#> drink alter                     -0.0617      -0.1882        0.5890
#> drink ego                        0.0545       0.0949        0.6290
#> drink similarity                 1.4396       1.8109        0.0190
#> drink linear shape               0.4190       0.4056        0.0770
#> drink quadratic shape           -0.6158      -0.3300        0.1140
#> drink average alter              1.3621       1.1070        0.1500
#> (*)Dummy2:outdegree (density)    0.0000      -0.4098        0.7850
#> (*)Dummy2:reciprocity            0.0000       0.9519        0.3260
#> (*)Dummy2:transitive triplets    0.0000       0.2066        0.5490
#> (*)Dummy2:3-cycles               0.0000      -0.3694        0.5490
#> (*)Dummy2:drink alter            0.0000       0.2313        0.1110
#> (*)Dummy2:drink ego              0.0000       0.0105        0.2370
#> (*)Dummy2:drink similarity       0.0000      -0.4896        0.7100
#> (*)Dummy2:drink linear shape     0.0000       0.0621        0.8050
#> (*)Dummy2:drink quadratic shape  0.0000      -0.4959        0.1610
#> (*)Dummy2:drink average alter    0.0000       0.2846        0.4750
#> 
#> Effect-wise joint significance tests
#> (i.e. each effect across all dummies):
#>                       chi-sq. df p-value
#> outdegree (density)      0.07  1   0.791
#> reciprocity              0.96  1   0.327
#> transitive triplets      0.36  1   0.549
#> 3-cycles                 0.36  1   0.549
#> drink alter              2.54  1   0.111
#> drink ego                1.40  1   0.237
#> drink similarity         0.14  1   0.708
#> drink linear shape       0.06  1   0.806
#> drink quadratic shape    1.97  1   0.160
#> drink average alter      0.51  1   0.475
#> 
#> Period-wise joint significance tests
#> (i.e. each period across all parameters):
#>          chi-sq. df p-value
#> Period 1    9.18 10   0.515
#> Period 2    9.18 10   0.515
#> 
#> Use the following indices for plotting:
#>   (1) outdegree (density)
#>   (2) reciprocity
#>   (3) transitive triplets
#>   (4) 3-cycles
#>   (5) drink alter
#>   (6) drink ego
#>   (7) drink similarity
#>   (8) drink linear shape
#>   (9) drink quadratic shape
#>   (10) drink average alter
#> 
#> If you would like to fit time dummies to your model,
#> use the includeTimeDummy function.
#> Type "?sienaTimeTest" for more information on this output.

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 demonstrates so you have the idiom when your own data needs it.

The full toolbox: Joint tests, the time-dummy remedy, more GOF, and the table your paper needs (optional self-study)

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

Multipar.RSiena(m2, 4, 5, 6)   # H0: recip = transTrip = cycle3 = 0
#> Tested effects:
#>  friend: reciprocity eval 
#>  friend: transitive triplets eval 
#>  friend: 3-cycles eval 
#> chi-squared = 219.42, d.f. = 3;  p < 0.001.
testSame.RSiena(m2, 7, 8)      # H0: drink alter effect == drink ego effect
#> Tested effects:
#>  friend: drink alter eval  ==  friend: drink ego eval 
#> chi-squared = 0.30, d.f. = 1; one-sided Z = -0.55; two-sided p = 0.581.

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 = 0.515), 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:

# 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")
#> ------------------------------------------------------------------
#> PERIOD-2 RECIPROCITY DUMMY: dummy est = 0.28   |   dummy t = 0.912   |   tconv.max = 0.183
#> ------------------------------------------------------------------

The dummy is the period-2 shift in reciprocity, and its t-ratio (0.91) 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 used three auxiliaries; the fourth stock one is OutdegreeDistribution, and descriptives.sienaGOF() prints the simulated envelope the GOF plot draws:

g_out <- sienaGOF(m2, OutdegreeDistribution, varName = "friend",
                  verbose = FALSE, join = TRUE)
round(descriptives.sienaGOF(g_out), 2)
#>                0     1     2     3     4      5      6      7      8
#> max        25.00 57.00 73.00 91.00 98.00 100.00 100.00 100.00 100.00
#> perc.upper 20.00 46.00 68.00 85.00 95.00  99.00 100.00 100.00 100.00
#> mean       12.82 35.75 57.28 75.15 87.99  94.99  98.17  99.39  99.83
#> median     13.00 36.00 57.00 75.00 88.00  95.00  99.00 100.00 100.00
#> perc.lower  7.00 26.00 46.00 64.00 78.00  88.00  94.00  97.00  98.00
#> min         2.00 22.00 37.00 58.00 68.00  77.00  86.00  94.00  95.00
#> sd          3.51  5.20  5.81  5.35  4.39   2.97   1.85   0.95   0.47
#> obs         8.00 30.00 55.00 76.00 93.00 100.00 100.00 100.00 100.00
#> p>          0.90  0.84  0.64  0.42  0.10   0.00   0.00   0.00   0.00
#> p>=         0.95  0.88  0.69  0.50  0.15   0.02   0.27   0.61   0.87

Outdegree fits too (p = 0.516), 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 built its table by hand so you’d know where $theta and $covtheta live. For the manuscript, one line:

siena.table(m2, type = "html", sig = TRUE, d = 3)
#> Results for m2 written to m2.html .

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.

# Optional extension space for a co-evolution question.
# Keep exploratory terms here rather than changing the core specification.

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")
#> ------------------------------------------------------------------
#> END BLOCK C: selection t = 2.349   |   influence t = 1.44   |   triad GOF p = 0.007
#> ------------------------------------------------------------------

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

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:

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

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

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

Did they report the Jaccard index, or otherwise establish the model applies? (A2)
Did they report the overall maximum convergence ratio, not just individual t-ratios, and is it under 0.25? (B2)
Is the influence standard error larger than the estimate? (C1)
Did they run GOF, and on which auxiliary statistics? Did any fail, and did they act on it? (C2)
Do they interpret coefficients as probabilities or marginal effects? (They shouldn’t - B3)
Do they claim selection versus influence from a null result on one of them? (C1)

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

What we did not show

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, 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’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: 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

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.

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

sessionInfo()
#> R version 4.3.3 (2024-02-29)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 24.04.3 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] RSiena_1.5.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] Matrix_1.6-5          ggnewscale_0.5.2      gtable_0.3.6         
#>  [4] jsonlite_2.0.0        dplyr_1.2.1           compiler_4.3.3       
#>  [7] tidyselect_1.2.1      Rcpp_1.1.1-1.1        parallel_4.3.3       
#> [10] scales_1.4.0          yaml_2.3.12           fastmap_1.2.0        
#> [13] coda_0.19-4.1         lattice_0.22-5        ggplot2_4.0.3        
#> [16] R6_2.6.1              generics_0.1.4        igraph_2.2.2         
#> [19] knitr_1.51            htmlwidgets_1.6.4     MASS_7.3-60.0.1      
#> [22] ggrepel_0.9.6         tibble_3.3.1          statnet.common_4.13.0
#> [25] netify_1.5.3          pillar_1.11.1         RColorBrewer_1.1-3   
#> [28] rlang_1.2.0           xfun_0.55             S7_0.2.2             
#> [31] otel_0.2.0            cli_3.6.6             magrittr_2.0.5       
#> [34] network_1.20.0        digest_0.6.39         grid_4.3.3           
#> [37] lifecycle_1.0.5       vctrs_0.7.3           evaluate_1.0.5       
#> [40] glue_1.8.1            farver_2.1.2          rmarkdown_2.30       
#> [43] tools_4.3.3           pkgconfig_2.0.3       htmltools_0.5.9

Check	Result	What a miss would lead us to examine
In-degree distribution	\(p=\) 0.589	Unmodeled differences in how popular students are
Behavior distribution	\(p=\) 0.912	A poor description of how drinking changes
Triad census	\(p=\) 0.007	Missing closure, cycles, or local subgroup structure