---
title: "Day 8: Netify, Building the Object and Mapping the Whole Course"
subtitle: "Advanced Network Analysis · ICPSR"
author: "Shahryar Minhas"
date: last-modified
engine: knitr
editor: source
format:
  html:
    theme: [cosmo, _theme.scss]
    toc: true
    toc-depth: 2
    toc-location: left
    toc-title: "Today"
    number-sections: true
    embed-resources: true
    code-copy: true
    code-tools: true
    fig-width: 8
    fig-height: 5
    fig-align: center
    df-print: kable
execute:
  echo: true
  warning: false
  message: false
  freeze: false
knitr:
  opts_chunk:
    comment: "#>"
---

::: {.callout-important title="Run This Setup First"}
Run the setup chunk below before beginning the walkthrough. It loads the packages, finds the Day 8 folder whether you are using the downloaded `day8_netify` folder or the internal course source folder, and defines the local data loader used throughout the document.
:::

```{r}
#| label: setup
#| message: false
#| warning: false
editor_path <- if (requireNamespace("rstudioapi", quietly = TRUE) &&
                   rstudioapi::isAvailable()) {
  tryCatch(rstudioapi::getSourceEditorContext()$path, error = function(e) "")
} else ""
day8_candidates <- unique(c(
  if (nzchar(editor_path)) dirname(editor_path) else character(),
  "day8_netify", "1_netify", "."
))
day8_candidates <- day8_candidates[
  vapply(file.path(day8_candidates, "_setup.R"), file.exists, logical(1))
]
if (!length(day8_candidates)) {
  stop("Cannot find the Day 8 folder. Open the Day 8 QMD from its downloaded folder, ",
       "or set the working directory to the folder containing the Day 8 materials.", call. = FALSE)
}
day8_dir <- normalizePath(day8_candidates[[1]], mustWork = TRUE)
source(file.path(day8_dir, "_setup.R"))
require_pkgs(c("netify", "igraph", "ggplot2", "network", "sandwich",
               "lmtest", "dplyr", "generics", "digest"))
suppressPackageStartupMessages({
  library(netify)
  library(igraph)
  library(ggplot2)
})
theme_set(theme_bw())
set.seed(6886)
```

::: {.callout-note appearance="simple"}
## How to use this document

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

Everything in here runs. There is no separate slide deck hiding the code. This *is* the script, with the reasoning left between the chunks. You can read it from top to bottom afterwards without losing the thread of the analysis.

A few conventions you will see all week:

- **Inline numbers are live.** When a sentence says a network has some density, that number was computed from the object on the line above it, not typed in by me. So if you re-run with different data, the prose updates too.
- **Expensive fits are cached.** Later days fit models that take minutes. Those are wrapped in `cache_fit()`, which stores the result in `cache/` so the first person pays for it once and everyone after loads it in milliseconds. Today has almost nothing expensive in it, so almost everything just runs.
- **No absolute paths.** Every file this document touches lives in its own folder. Unzip `1_netify/`, open this, and it runs.

Today you need three packages, all on CRAN:

```r
install.packages(c("netify", "igraph", "ggplot2", "statnet", "sandwich",
                   "lmtest", "dplyr", "digest"))
```
:::

# Where we are: Day 8 {#sec-where}

Welcome back. The first seven days gave you the foundations of network analysis. Today begins the advanced modeling sequence: getting a coefficient, an interval, a prediction, and a test while taking network dependence seriously. We will get to that problem in the second half of today, and it is the thread that ties Days 8 through 14 together.

Here is where we are going, so the plumbing has a point. Before the first break I will show you a perfectly ordinary-looking regression, the kind in every methods textbook, that screams "significant" **about 50% of the time when the true effect is exactly zero**. Everything that follows is about how you build the object so that you are not the person who published it.

But before you can model a network you have to *have* one, as an object, in the right shape, with the right things attached. That sounds like plumbing, and half of it is. The other half is a series of quiet decisions: Is a tie directed? Is it weighted? Who counts as an actor? What does a missing dyad mean? Every one of those decisions changes the model you get at the end. Get the object wrong and no amount of clever modeling later will save you.

So today has two parts:

- **Part one: the object.** We build networks from raw data with `netify`, justify every argument as we go, attach node- and dyad-level covariates, read what the object is telling us, plot it, and hand it off in the exact shape each modeling package wants. This is the substrate the rest of the course stands on.
- **Part two: the map.** Once you have the object, what do you actually *do* with it? There is a small map of methods, one cell per day for the rest of the week, and it all comes from a single problem. We will draw that map and point at Day 9.

## Let's Start With Your Research Question {#sec-political-question}

Before `netify()` gets an argument, the research question has already made most of the important decisions. Consider three familiar political science projects:

- **Alliance formation:** states are the actors, a defense commitment is the tie, direction depends on how the treaty is coded, and the risk set excludes state-pairs that did not coexist.
- **Legislative collaboration:** legislators are the actors, cosponsorship or joint voting is the tie, and a Congress or legislative session defines the time window.
- **Conflict and cooperation:** states or armed groups are the actors, events are the raw records, and aggregation determines whether the outcome is any interaction, the number of interactions, or their intensity.

Those are not cosmetic differences. They produce different outcomes and different populations. A claim about why states initiate conflict needs a directed network. A claim about whether alliances cluster needs a binary or weighted alliance definition. A claim about change needs comparable actor rosters and time periods.

The construction sequence for every example today is:

1. **Name the substantive event.** What happened between which political actors?
2. **Define the tie.** Decide its direction, scale, and time window.
3. **Define eligibility.** Determine which actor-pairs could have produced an observed tie.
4. **Build the network outcome.** Convert the raw records into an adjacency matrix or longitudinal array.
5. **Attach covariates by actor identity.** Do not rely on row order.
6. **Audit the result.** Check actors, ties, missingness, direction, and dimensions before fitting anything.

The code in this session follows those six steps. The point is not simply to create an R object. It is to make the estimand visible before a model quietly inherits it.

## The Applied Compass for Today {#sec-model-compass}

Even today, before the main modeling sequence begins, every tool has to earn its place in an applied argument.

| Tool or approach | Question we want to understand | What it contributes | What makes it distinct | What the result sounds like in actual units |
|---|---|---|---|---|
| Network construction | Who could relate to whom, what counts as a tie, and over what period? | Turns political records into a defensible outcome and risk set | It defines the estimand before any estimator sees the data | "One cell records how many cooperative statements country A directed toward country B during one calendar year. A value of 12 means twelve coded statements. A zero means both countries existed and no statement was recorded. A missing value means we cannot treat that country-pair-year as an observed opportunity with no cooperation." |
| Descriptive network summaries | Which countries initiated the most actions, which countries were the targets of the most actions, and which pairs acted toward one another? | Makes the relational pattern visible and suggests puzzles worth modeling | It describes the observed graph without claiming a data-generating process | "An out-strength of 240 means a country initiated 240 coded cooperative statements across all partner countries that year. An in-strength of 240 means other countries directed 240 such statements toward it. Those are different political actions, not two labels for the same kind of activity." |
| Naive dyadic regression | How are measured country and country-pair characteristics associated with the number of actions from one country toward another if country-pairs are treated as conditionally independent? | Provides a familiar benchmark and exposes what the independence assumption buys | It uses ordered country-pairs as rows and does not represent shared-country or broader network dependence | "If the outcome is a count of cooperative statements, a coefficient is the expected change in that count for one ordered country-pair when the predictor rises by one unit, holding the listed covariates constant. We should not yet trust its interval because one country contributes a row for every potential partner." |
| Country-clustered uncertainty preview | How much does recognizing that the same country initiates actions toward many partners change uncertainty in that benchmark? | Recalculates uncertainty under one limited dependence pattern | It changes the covariance estimate, not the fitted mean or the substantive dependence model | "The estimated change in the number of cooperative statements stays exactly the same. Its interval changes after we tell the calculation that all actions initiated by the same country may share unmeasured causes. This still does not account for repeated targets, reciprocal action, or dependence involving third countries." |

The question at the end of every code block is therefore not only "Did it run?" It is "What social fact did this object or estimator let us describe, and what did it leave unresolved?"

Let's start abstract, because not every project in this course studies countries.

One more convention: the collapsible boxes titled **"the full toolbox"** are self-study. They contain the rest of the netify surface in a form you can run at home, and nothing in the main path depends on them.

# A tie is a decision, not a cell in a spreadsheet {#sec-abstract}

Strip away the substance and every network in this course is the same object: a set of **actors**, and a set of **ties** running between them. A tie has a **sender** ($i$) and a **receiver** ($j$). That is the whole vocabulary. It does not matter whether $i$ and $j$ are students who named each other as friends, firms that traded, brain regions that co-activated, or countries that signed a treaty: the machinery is identical, and we will say "sender" and "receiver" all week so nobody has to translate out of political science.

The smallest useful example is a classroom. `netify` ships one: thirty students, and a list of who named whom as a friend. Friendship nomination is directed: you can name someone who does not name you back, so this is a directed network, and it is a good place to meet the mechanics before we touch anything with a foreign ministry in it.

```{r}
#| label: classroom-load
classroom_edges <- load_day8_data("classroom_edges", day8_dir)
classroom_nodes <- load_day8_data("classroom_nodes", day8_dir)
head(classroom_edges)
```

This is an **edgelist**: one row per tie, a `from` and a `to`. This is the rawest form network data comes in, and it is what `netify()` most wants to eat. The nodes come with attributes too, including gender, grade, and GPA, which we will attach in a moment.

```{r}
#| label: classroom-build
class_net <- netify(
  classroom_edges,
  actor1    = "from",     # The sender column
  actor2    = "to",       # The receiver column
  symmetric = FALSE       # Nomination is directed: i -> j is not j -> i
)
class_net
```

That print-out is worth reading slowly, because it is `netify` telling you what it thinks you just built. It is **unipartite** (one kind of actor, not students-and-clubs), **asymmetric** (directed, because we said so), **binary** (a tie either exists or it does not; there are no weights yet), and **cross-sectional** (one snapshot, no time). Thirty actors. Those four facts are the identity of the object, and every one of them was a decision you made or accepted a default on. A small Rosetta note, because the labels on screen are not always the words we say: `netify` prints **Asymmetric** (= directed), **Binary Weights** (= unweighted), and **Cross-Sectional** (= one time point).

::: {.callout-warning}
## The message about `missing_to_zero` is a modeling decision in disguise

`netify` mentions that missing dyads were filled with zeros. For *describing* a classroom that is usually fine: if two students are not listed, they are not friends. But "we never observed this pair" and "this pair is definitely not tied" are different claims about the world, and three arguments control which one you are making: `missing_to_zero`, `diag_to_NA`, and `nodelist`. They are cheap to see on a toy, so let's not take them on faith:

```{r}
#| label: missing-diag-toy
toy <- data.frame(
  from = c("ana", "ana", "ben", "cal"),
  to   = c("ben", "cal", "ana", "ana"),
  mins = c(20, 5, 12, 3)            # Minutes talked, say
)

# The defaults: unobserved pairs become 0, the diagonal becomes NA
toy_zero <- netify(toy, actor1 = "from", actor2 = "to",
                   symmetric = FALSE, weight = "mins")
peek(toy_zero, from = 1:3, to = 1:3)

# The preserve-uncertainty version: unobserved pairs stay NA, and
# `nodelist` adds dee, an actor we know exists but who has no rows
toy_na <- netify(toy, actor1 = "from", actor2 = "to",
                 symmetric = FALSE, weight = "mins",
                 missing_to_zero = FALSE,
                 nodelist = c("ana", "ben", "cal", "dee"))
peek(toy_na, from = 1:4, to = 1:4)
```

(`peek()` gets a proper introduction in the next section: it just shows you a corner of the matrix.) Three things moved. Ben→cal flipped from `0` to `NA`: we never saw that pair, and the default zero was a *claim*, not an observation. The diagonal is `NA` in both because `diag_to_NA = TRUE` by default: self-ties are usually meaningless, and you flip it for the networks where they are not. And dee, an isolate with no rows in the edgelist, exists at all only because `nodelist` put her there: without it her non-ties would silently vanish from every density you compute. On a later day the model will believe whichever claim you encoded, so make these three deliberately.
:::

# Reading the object {#sec-read-object}

You now have a network. Before the summary tools, one flat truth to keep the object from ever feeling like a black box: a netify object *is* the adjacency matrix, with bookkeeping attached. `peek()` shows you any corner of it without drowning the console in a thirty-by-thirty print:

```{r}
#| label: classroom-peek
peek(class_net, from = 1:5, to = 1:5)
```

Rows are senders, columns are receivers, a `1` is a nomination, and the diagonal is `NA` because you cannot name yourself. Any time an object confuses you this week, `peek()` first: it is the fastest way to see what you actually built.

Now the two functions that tell you almost everything else, both of which you will use every single day of this course.

`summary()` describes the network as a *whole*:

```{r}
#| label: classroom-summary
summary(class_net)
```

Read a few of those off. The density is `r round(summary(class_net)$density, 3)`. Of all the directed pairs that *could* be friends, that fraction actually are, so this is a sparse network, as classrooms are. Reciprocity here is `r round(summary(class_net)$reciprocity, 3)`, and the sign matters: `netify` reports reciprocity as the *correlation* between the $i\to j$ and $j\to i$ tie values, so it runs from $-1$ to $+1$, not $0$ to $1$. At `r round(summary(class_net)$reciprocity, 2)` it is essentially zero: if anything a hair *below* chance. In fact not a single one of these nominations is returned: `summary(class_net)$mutual` is exactly `r summary(class_net)$mutual`. That these students name friends who do not name them back is already a substantive feature of this network.

`summary()` prints a dozen more columns; ignore them for now: density and reciprocity are all you need today, and the rest arrive when the models that use them do.

`summary_actor()` drops one level down and describes each *actor*:

```{r}
#| label: classroom-actor
sa_class <- summary_actor(class_net)
head(sa_class[, c("actor", "degree_in", "degree_out", "degree_total")])
```

`degree_out` is how many people a student named; `degree_in` is how many named *them*. Those are different questions: being friendly and being popular are not the same, and the fact that a directed network lets you separate them is exactly the kind of structure we will spend the week modeling. The most "popular" student here receives `r max(sa_class$degree_in)` nominations while the busiest nominator sends out `r max(sa_class$degree_out)`.

Now attach the node attributes, so those degrees can be *explained* rather than just ranked:

```{r}
#| label: classroom-attach-nodes
class_net <- add_node_vars(
  class_net,
  classroom_nodes,
  actor     = "student",              # The id column in classroom_nodes
  node_vars = c("gender", "grade", "gpa")
)
class_net
```

The object now reports **nodal features**: gender, grade, and GPA are riding along with each actor. When we get to modeling, "do students befriend others of the same gender?" becomes a coefficient, and it will need those attributes to be attached in exactly this way.

Finally, look at it. A picture of a thirty-node network is not a rigorous object, but it is how you catch the gross mistakes: an isolated component you did not expect, one actor tied to everyone, a direction reversed. And because the attributes are attached, the picture can *carry* them: any nodal feature can drive any aesthetic through the `node_*_by` arguments.

```{r}
#| label: classroom-plot
#| fig-cap: "The friendship network. Nodes are students, colored by gender and sized by GPA; arrows point from the nominator to the nominated. Even by eye you can see a handful of much-named students and several who only send."
plot(class_net, node_color_by = "gender", node_size_by = "gpa")
```

::: {.callout-tip}
## Checkpoint: you can already do the core loop

If you can read the four lines of that print-out (Unipartite / Asymmetric / Binary Weights / Cross-Sectional: `netify`'s own words for one-mode / directed / unweighted / one time point), run `summary()` and `summary_actor()`, and attach node attributes, you have the whole descriptive loop. Everything from here is the same loop on messier, more interesting data.
:::

::: {.callout-note collapse="true"}
## The full toolbox: know your object (self-study)

The debugging vocabulary for when an object you built is not what you thought it was. None of this is analysis; all of it saves you at 11pm.

```{r}
#| label: toolbox-know-object
# What am I holding? Seven predicates ...
c(netify       = is_netify(class_net),
  binary       = is_binary(class_net),
  longitudinal = is_longitudinal(class_net),
  multilayer   = is_multilayer(class_net),
  bipartite    = is_bipartite_netify(class_net),
  directed     = is_directed_netify(class_net),
  symmetric    = is_symmetric_netify(class_net))

# ... and three dimension counters (these work on ANY netify object,
# and they are the documented way to ask, not length())
c(actors = n_actors(class_net), periods = n_periods(class_net),
  layers = n_layers(class_net))

# The attached node attributes, back as a data.frame
head(nodal_data(class_net), 3)

# A full coherence check, run it when something feels off; the usual
# culprit it catches is an actor-name mismatch between your tables
all(unlist(validate_netify(class_net, verbose = FALSE)))

# And the escape hatch: it really is just a matrix
m <- as.matrix(class_net)
m[1:3, 1:3]
```

`get_raw()` is the same escape hatch for any shape: on a longitudinal object it hands you the full actor × actor × time array.
:::

# A two-mode interlude: bipartite networks {#sec-bipartite}

Everything so far has one kind of actor. Plenty of you have data with two kinds, such as students and clubs, people and events, firms and boards, or words and documents, where ties run only *between* the modes and never within them. That is a **bipartite** network, and `netify` builds it the moment you say so:

```{r}
#| label: bipartite-build
affil <- data.frame(
  student = c("sam", "sam", "lena", "lena", "omar", "omar", "noor"),
  club    = c("chess", "debate", "debate", "robotics", "chess",
              "robotics", "debate")
)
bip_net <- netify(affil, actor1 = "student", actor2 = "club",
                  mode = "bipartite")
bip_net
```

The print-out now says **bipartite**, and rows and columns are different kinds of thing:

```{r}
#| label: bipartite-peek
peek(bip_net, from = 1:4, to = 1:3)
```

The matrix is a rectangle, not a square: rows are students, and columns are clubs. The familiar "who shares a club with whom" projection is a *choice you make later*, not something the object does behind your back. A rectangle also reads naturally as tiles:

```{r}
#| label: bipartite-heatmap
#| fig-width: 5
#| fig-height: 3.5
#| fig-cap: "The affiliation matrix as a heatmap. For two-mode data this is often the more direct picture: no layout algorithm, no projection, just the matrix."
plot(bip_net, style = "heatmap")
```

We stay unipartite for the rest of the week, but if your data are two-mode you are not locked out of the modeling half of this course: Day 10's AME machinery has a bipartite variant, and this object is exactly what it eats.

# From events to a network {#sec-events}

Classrooms come pre-cleaned. Observed data almost never does. A common raw form in the social sciences is **event data**: a long log of things that happened between pairs of actors, one row per event, and it is *not* a network until you decide how to turn events into ties.

`netify` ships an example from ICEWS, a dataset of interactions between countries coded from news reports, with one row per directed country-pair per year. "Country" is just the substantive name for "actor" here: if you study firms or people, read `i` and `j` as your senders and receivers.

```{r}
#| label: icews-load
icews <- load_day8_data("icews", day8_dir)
dim(icews)
head(icews[, c("i", "j", "year", "verbCoop", "matlConf", "i_polity2")])
```

Nearly three hundred thousand rows. Each has a sender `i`, a receiver `j`, a `year`, and several *counts*: `verbCoop` is how many cooperative statements $i$ made toward $j$ that year, `matlConf` how many materially conflictual acts, and so on. There are also attributes measured on the sender: `i_polity2` (a regime-type score), `i_log_gdp`, `i_log_pop`.

We are not going to model a hundred and fifty countries today: big networks plot badly and distract from the mechanics. So pick a coherent set of major players. This is the first substantive **actor-handling** decision, and it is yours to make explicitly:

```{r}
#| label: icews-subset
keep <- c(
  "United States", "China", "Russian Federation",
  "United Kingdom", "France", "Germany", "Japan", "India",
  "Brazil", "Israel", "Iran, Islamic Republic Of", "Iraq",
  "Egypt", "Saudi Arabia", "Pakistan", "Afghanistan",
  "Korea, Republic Of", "Syrian Arab Republic"
)
ic <- icews[icews$i %in% keep & icews$j %in% keep, ]
c(rows = nrow(ic), actors = length(unique(ic$i)),
  years = length(unique(ic$year)))
```

::: {.callout-note collapse="true"}
## Note the country names are exact strings, and that bites everyone once

`"Iran, Islamic Republic Of"`, not `"Iran"`. `"Korea, Republic Of"`, not `"South Korea"`. If you type the friendly name, `netify` will not error: it will simply build a network that silently omits Iran, and you may not notice until much later. Actor identity is *string identity*, and reconciling names across sources (ISO codes, COW codes, UN names) is an important and under-taught part of this work. When in doubt, `sort(unique(c(ICEWS$I, ICEWS$j)))` and read the list.
:::

Now the two decisions that define the tie. **Is it directed, and is it weighted?**

```{r}
#| label: icews-cross-section
ic14 <- ic[ic$year == 2014, ]

net14 <- netify(
  ic14,
  actor1    = "i",
  actor2    = "j",
  symmetric = FALSE,       # Cooperation i->j need not equal j->i
  weight    = "verbCoop"   # The tie carries a count, not just 0/1
)
net14
```

Two arguments carried all the meaning there:

- **`symmetric = FALSE`.** The United States can direct far more cooperative signals at a small state than it gets back. If you set `symmetric = TRUE` on directed data, `netify` mirrors one triangle of the matrix onto the other and *throws half your information away*: quietly. So the rule is: symmetry is a claim about the world, not a convenience. Set it to match how the data were generated.
- **`weight = "verbCoop"`.** With no weight you get a binary tie: did $i$ ever cooperate with $j$? With a weight, the tie carries *how much*. Which you want depends on your question, and it changes the model family downstream (a count is not a coin flip). Here we keep the count.

The density is now `r round(summary(net14)$density, 3)`: among these major powers essentially everyone signals everyone, which is exactly why a *binary* network of great-power cooperation would be boring and a *weighted* one is not.

# `sum_dyads`: when a pair shows up more than once {#sec-sumdyads}

Here is the argument people skip and then get wrong. We just built a single year. But `netify` needs one tie per pair, so what happens when the raw data has the *same pair on multiple rows*? That is not a corner case; it is what all event data looks like before you aggregate it.

Watch what happens if we hand `netify` all thirteen years at once but *do not tell it about time*. Each country-pair now appears up to thirteen times:

```{r}
#| label: sumdyads-count
pair <- paste(ic$i, ic$j)
table(table(pair))   # How many rows does each directed pair get?
```

Every pair shows up thirteen times. `netify` cannot put thirteen numbers in one cell of the matrix, so it makes you decide with **`sum_dyads`**:

```{r}
#| label: sumdyads-build
net_pooled <- netify(
  ic,
  actor1    = "i",
  actor2    = "j",
  symmetric = FALSE,
  weight    = "verbCoop",
  sum_dyads = TRUE     # Collapse the repeated rows by summing the weight
)
net_pooled
```

`sum_dyads = TRUE` says "these repeated rows are pieces of one tie: add them up." Now the weight on the tie from the United States to China is *total* verbal cooperation from 2002 to 2014. That is a legitimate object: a pooled, cross-sectional summary of a decade. The mean tie weight jumps to `r round(summary(net_pooled)$mean_edge_weight, 0)`, roughly thirteen times the single-year figure, exactly because we stacked thirteen years into it.

::: {.callout-important}
## `sum_dyads` is the difference between "a tie" and "an average tie"

Whenever more than one row maps to the same pair, you *must* decide what the tie means: the sum of the events (`sum_dyads = TRUE`) or, if you aggregate by hand first, their average, maximum, or most recent value. `netify` will not guess for you, and the default (`FALSE`) assumes each pair is already unique. In `netify` 1.5.3, repeated weighted dyads trigger a warning and an automatic switch to `sum_dyads = TRUE`, so the rows are summed rather than dropped. Do not rely on that fallback: make the choice explicit so your code says which world you are in.
:::

::: {.callout-note collapse="true"}
## The full toolbox: event data before it becomes a network (self-study)

`sum_dyads = TRUE` aggregates *inside* `netify()`. Sometimes you want to aggregate first and inspect the result before you commit to a network. `aggregate_dyad()` is the standalone version of the same move. `netify` ships UCDP records of organized violence in Mexico for practicing it: one row per recorded incident, with a `best` estimate of deaths.

```{r}
#| label: toolbox-aggregate-dyad
mexico <- load_day8_data("mexico", day8_dir)
mex_agg <- aggregate_dyad(mexico, actor1 = "side_a", actor2 = "side_b",
                          time = "year", weight = "best",
                          symmetric = FALSE)
head(mex_agg, 3)
mex_net <- netify(mex_agg, actor1 = "side_a", actor2 = "side_b",
                  time = "year", symmetric = FALSE, weight = "best")
mex_net
```

`r format(nrow(mexico), big.mark = ",")` raw incident rows collapse to `r nrow(mex_agg)` dyad-years *before* the network exists, and you got to read the aggregation on the way through, which is exactly what you want when the aggregation rule itself is a modeling decision. Should you use the sum, the maximum, or the most recent value? One companion utility: when the aggregation is *undirected* (`symmetric = TRUE` treats i→j and j→i as the same pair), `gen_symm_id()` builds the order-free pair id that makes that work: both orderings of a pair get the same id:

```{r}
#| label: toolbox-gen-symm-id
ex <- data.frame(i = c("usa", "china"), j = c("china", "usa"))
gen_symm_id(ex, actor1 = "i", actor2 = "j")   # Same id, both rows
```

(`myanmar`, a second UCDP event dataset, also ships with the package if you want a homework-sized playground.)
:::

But pooling a decade into one number also *throws away the time*. Usually the time is the interesting part. So instead of collapsing the years, let's keep them.

::: {.brk}
Break. [You have built cross-sectional networks, learned what `symmetric`, `weight`, and `sum_dyads` each claim about the world, and met the caching pattern. After the break: time, covariates, and handing the object to a model.]{.sub}
:::

# Going longitudinal {#sec-longitudinal}

Give `netify` a `time` column and it stops making one network and starts making one network *per period*, lined up so the same actor is the same row in every year. This is the object that days 11 and 12 need.

```{r}
#| label: longitudinal-build
net_long <- netify(
  ic,
  actor1    = "i",
  actor2    = "j",
  time      = "year",       # <- The only new argument
  symmetric = FALSE,
  weight    = "verbCoop"
)
net_long
```

No `sum_dyads` needed now: once time is a dimension, each pair *is* unique within a year. The object reports **`r n_periods(net_long)` periods** (`n_periods()` is how you ask; `n_actors()` and `n_layers()` are its siblings), and `summary()` on a longitudinal object gives you one row per period. But great-power *cooperation* is nearly saturated every single year: its density sits near one and barely budges, so it is the wrong place to watch a network move. Read the same countries on material *conflict* instead:

```{r}
#| label: longitudinal-summary
net_long_conf <- netify(
  ic, actor1 = "i", actor2 = "j", time = "year",
  symmetric = FALSE, weight = "matlConf"   # Conflict, not cooperation
)
sl <- summary(net_long_conf)
sl[, c("net", "num_actors", "density", "mean_edge_weight")]
```

Now the panel *remembers the news*. Density climbs from `r round(sl$density[sl$net == 2002], 2)` in 2002 to `r round(sl$density[sl$net == 2014], 2)` by 2014, and the mean edge weight, the average amount of conflict carried on an active tie, jumps from `r round(sl$mean_edge_weight[sl$net == 2002], 0)` in 2002 to `r round(sl$mean_edge_weight[sl$net == 2003], 0)` in 2003. Watch 2003: that spike is the Iraq invasion landing in the object. A single cross-section would have hidden it; the longitudinal object lets a historical event show up as a number.

And because the object knows it is a panel, the tools you already have turn into figures for free. `plot()` facets by period, filtered here to four years so each panel remains legible, and `plot()` on the summary table turns its columns into time series:

```{r}
#| label: longitudinal-facets
#| fig-height: 6
#| fig-cap: "Material conflict among the same eighteen countries in four selected years: same actors, one facet per period. The 2003 panel is the Iraq invasion, visible as structure rather than as a row in a table."
plot(net_long_conf,
     time_filter   = as.character(c(2002, 2003, 2009, 2014)),
     mutate_weight = log1p)   # Event counts are skewed; log for readability
```

```{r}
#| label: longitudinal-stats-plot
#| fig-height: 6
#| fig-cap: "The summary() table as time series. The table above and this figure are the same object."
plot(sl)
```

One more question a panel invites: *how much does the network change year to year?* `compare_networks()` can return several pairwise comparisons at once. Here we request `method = "all"`, so it computes the ordinary edge-pattern correlations we display below and also runs a quadratic assignment procedure (QAP) test that permutes actor labels jointly across rows and columns:

```{r}
#| label: longitudinal-compare
cn <- compare_networks(net_long_conf, method = "all",
                       n_permutations = 500, seed = 6886,
                       return_details = TRUE)
cn
```

```{r}
#| label: longitudinal-compare-detail
cor_mat <- cn$details$correlation_matrix
c(`2002 vs 2003` = round(cor_mat["2002", "2003"], 2),
  `2002 vs 2014` = round(cor_mat["2002", "2014"], 2))
```

The 2002-2003 pair correlates at `r round(cn$details$correlation_matrix["2002", "2003"], 2)`; across the full twelve-year gap that drops to `r round(cn$details$correlation_matrix["2002", "2014"], 2)`. The conflict network has memory, although one adjacent pair is not enough to establish a smooth rate of decay, which is exactly the structure Day 12 (SAOM) turns into a model instead of a description.

::: {.callout-note collapse="true"}
## A reference for comparing networks

`compare_networks()` can compare four different things, and the choice should follow the question. `what = "edges"` compares the dyad cells themselves; `what = "structure"` compares graph-level summaries such as density and reciprocity; `what = "nodes"` checks who entered or left; and `what = "attributes"` compares the distributions of attached nodal attributes. For a longitudinal object it compares **all pairs of waves**, not only adjacent waves.

The edge metrics are not interchangeable:

- **Correlation** asks whether the two sets of edge weights rise and fall together. Shared zeros count, so a sparse pair of networks can correlate highly even when the positive ties overlap only modestly.
- **Jaccard** thresholds the weights and asks what fraction of ties present in at least one network appear in both. For a weighted network the default `edge_threshold = 0` means any positive weight is a tie.
- **Hamming** reports the fraction of eligible dyads whose thresholded tie status differs.
- **QAP** compares the observed Pearson edge correlation with a reference distribution produced by relabeling actors jointly across rows and columns. It is a test against that relabeling null, not a correction for omitted variables and not a causal design.

The same object can answer the structure and composition questions without being rebuilt:

```{r}
#| label: toolbox-compare-what
structure_over_time <- compare_networks(net_long_conf, what = "structure")
actors_over_time <- compare_networks(net_long_conf, what = "nodes")
structure_over_time$summary
actors_over_time$summary
```

Before interpreting any comparison, check that the networks use the same actor identities, direction, observation window, weight scale, and missing-data rule. An actor who did not exist is not the same as an actor with zero ties, and a one-month event count is not directly comparable to a one-year count.
:::

::: {.callout-note collapse="true"}
## `actor_time_uniform`: who is allowed to appear and when

We got a clean panel because all eighteen of these countries appear in all thirteen years. Empirical panels are often ragged: states enter, firms merge, and people drop out of a survey wave. `netify` handles this with `actor_time_uniform`. Left at its default it assumes your actor set is the same every period; set it `FALSE` and composition is inferred from when each actor appears in the data, treating an absent actor as *absent* rather than as present-with-no-ties. Those are different claims, and a toy makes the difference visible: dee only shows up in year 2:

```{r}
#| label: actor-entry-exit
toy_l <- data.frame(
  i  = c("ana", "ben", "ana",   "ben", "dee", "ana"),
  j  = c("ben", "ana", "cal",   "cal", "ana", "dee"),
  yr = c(1, 1, 1,   2, 2, 2))

uniform <- netify(toy_l, actor1 = "i", actor2 = "j", time = "yr",
                  symmetric = FALSE)                    # The default
ragged  <- netify(toy_l, actor1 = "i", actor2 = "j", time = "yr",
                  symmetric = FALSE, actor_time_uniform = FALSE)

summary(uniform)[, c("net", "num_actors", "density")]
summary(ragged)[,  c("net", "num_actors", "density")]

peek(uniform, from = 1:4, to = 1:4, time = 1)  # dee: a row of zeros
peek(ragged,  from = 1:4, to = 1:4, time = 1)  # dee: not there at all
```

Same data, different year-1 density, because "Dee sent no ties" and "Dee did not exist yet" are different denominators. If the automatic inference is not the composition you want, `actor_pds` lets you hand `netify` an explicit actor-by-period record instead. The model downstream will believe whichever claim you encoded; Day 12 (SAOM) leans on this directly when actors join and leave between waves.
:::

Now the caching pattern, once, on the one object that takes more than a blink. Building the network for *all* one-hundred-fifty-odd countries across all thirteen years requires appreciable computation:

```{r}
#| label: cache-demo
net_world <- cache_fit("net_world", netify(
  icews,
  actor1    = "i",
  actor2    = "j",
  time      = "year",
  symmetric = FALSE,
  weight    = "verbCoop",
  nodal_vars = c("i_polity2", "i_log_gdp", "i_log_pop")
), dir = file.path(day8_dir, "cache"))
c(actors = n_actors(net_world), periods = n_periods(net_world))
```

The *call* stays right there on the screen: there is no hidden second version that ran. `cache_fit()` hashes that code and saves the result to the Day 8 `cache/net_world.rda` file, then reloads it in milliseconds until you change the code. At that point it notices and rebuilds, which is how every slow model fit later in the week ships without stalling a live class.

Once the validated world-sized object exists, you rarely rebuild from raw data again: you *carve*. `subset_netify()` takes any combination of actors, periods, and layers and returns a new netify object with the covariates still attached:

```{r}
#| label: subset-netify
net_recent <- subset_netify(net_world,
                            actors = keep,                      # Our 18 majors
                            time   = as.character(2010:2014))   # Last 5 years
net_recent
```

This is the pattern every later day's data prep starts with: build wide once, then subset down to the estimation window you can defend. (Actor names are string-exact here too: misspell `"Iran, Islamic Republic Of"` and it is silently not in your subset.)

# Attaching covariates {#sec-covariates}

A network you only ever describe can stay bare. A network you want to *model* needs its covariates attached, and there are two straightforward ways to do it.

**Inline, at build time**: name the columns and let `netify` pull them:

```{r}
#| label: attach-inline
net_cov <- netify(
  ic14,
  actor1     = "i",
  actor2     = "j",
  symmetric  = FALSE,
  weight     = "verbCoop",
  nodal_vars = c("i_polity2", "i_log_gdp"),   # Measured on the actor
  dyad_vars  = "matlConf",                    # Measured on the pair
  dyad_vars_symmetric = FALSE                 # Conflict i->j != j->i
)
net_cov
```

Two kinds of covariate, and the distinction is the whole ballgame later today:

- **Nodal** covariates attach to an *actor*: regime type, GDP, population. One value per country. These get copied into every dyad that country is in, and that copying is precisely what will cause the trouble in part two.
- **Dyadic** covariates attach to a *pair*: material conflict from $i$ to $j$, distance between them, whether they share an alliance. One value per ordered pair, and `dyad_vars_symmetric` tells `netify` whether the $i\to j$ value equals the $j\to i$ value (distance: yes; conflict: no).

**Or after the fact**, when your covariates live in a separate table, which they usually do:

```{r}
#| label: attach-after
node_tab <- unique(ic14[, c("i", "i_polity2", "i_log_gdp")])
names(node_tab)[1] <- "country"

net_bare <- netify(ic14, actor1 = "i", actor2 = "j",
                   symmetric = FALSE, weight = "verbCoop")
net_bare <- add_node_vars(
  net_bare, node_tab, actor = "country",
  node_vars = c("i_polity2", "i_log_gdp")
)
net_bare <- add_dyad_vars(
  net_bare, ic14, actor1 = "i", actor2 = "j",
  dyad_vars = "matlConf", dyad_vars_symmetric = FALSE
)
net_bare
```

Same object, either road. `add_node_vars()` and `add_dyad_vars()` matter because they do the *alignment* for you: they match on actor identity, not on row order, so a covariate can never silently end up attached to the wrong actor. If you have ever built a design matrix by `cbind()`-ing two tables that turned out to be sorted differently, you know why that guarantee is worth having.

::: {.callout-note collapse="true"}
## The full toolbox: reshaping and rescuing objects (self-study)

Five utilities you will not need today but will want the week you meet messy data, plus the doors in from other formats:

```{r}
#| label: toolbox-manipulation
# 1. Drop actors with missing covariates (returns the cleaned object)
drop_na_actors(net_cov)

# 2. Combine networks built separately, `align_actors` decides how
#    the two actor sets get reconciled ("union" keeps everyone)
eu   <- subset_netify(net14, actors = c("United Kingdom", "France", "Germany"))
asia <- subset_netify(net14, actors = c("China", "Japan", "India"))
bind_netifies(eu, asia, align_actors = "union")

# 3. Swap a dyadic covariate in as THE network, demoting the old
#    network to a covariate, here conflict becomes the outcome
pivot_dyad_to_network(net_cov, dyad_var = "matlConf")

# 4. Decompose back into plain data.frames (edges + nodes)
dec <- decompose_netify(net_cov)
head(dec$edge_data, 3)

# 5. Arrived with an igraph or network object instead of an edgelist?
#    to_netify() is the door in (round-tripped here for the demo)
to_netify(to_igraph(net_cov), weight = "weight")
```

(In the `bind_netifies()` output, the cross-bloc cells are *missing*, not zero: the two source networks never observed each other's dyads, which is `missing_to_zero` logic doing the right thing.) And if your collaborator hands you a Gephi or Pajek file, there are readers for those too:

```{r}
#| label: toolbox-readers
#| eval: false
read_graphml("your_file.graphml")
read_gml("your_file.gml")
read_pajek("your_file.net")
```
:::

# Actor Statistics and Plotting with ICEWS {#sec-plotting}

Back to `summary_actor()`, now on something with more going on than a classroom. On a directed, weighted network it gives you the full menu: in- and out-degree, strength (degree weighted by tie value), betweenness, closeness, hub and authority scores:

```{r}
#| label: icews-actor-stats
sa14 <- summary_actor(net14)
head(sa14[order(-sa14$strength_sum_out),
          c("actor", "degree_total", "strength_sum_out",
            "strength_sum_in")])
```

Sort by out-strength and you get the states that initiated the largest total number of cooperative statements toward other states. Sort by in-strength and you get the states that were the targets of the largest total number of cooperative statements. **`r sa14$actor[which.max(sa14$strength_sum_out)]`** initiated the most in this network. These rankings can name different states because initiating diplomatic cooperation and being the target of it are different political actions. That distinction is the seed of the entire modeling week. Remember it for part two.

Now the picture. Two small moves come first because raw event counts are so right-skewed that a direct plot is one blinding U.S.-China edge and seventeen ghosts. `mutate_weights()` transforms the tie weights (the original column survives as a dyad feature, so nothing is lost), and `add_node_vars()` will attach *anything*, including the actor statistics you just computed, so the figure can be driven by the statistics:

```{r}
#| label: icews-transform
net14_log <- mutate_weights(net14, transform_fn = log1p,
                            new_name = "log_verbCoop")
net14_log <- add_node_vars(net14_log, sa14, actor = "actor",
                           node_vars = "strength_sum_out")
```

```{r}
#| label: icews-plot
#| fig-cap: "Great-power cooperation in 2014. Node size is out-strength, the summary_actor() column we just attached, and edge transparency carries the logged tie weight. The point is not to read individual ties but to see that a few actors sit at the centre of a dense core."
plot(net14_log, node_size_by = "strength_sum_out")
```

That is the whole aesthetic grammar: `node_size_by`, `node_color_by`, `edge_alpha_by`, `remove_isolates`, and a `layout` argument for the algorithm. This is the same mapping you used to color the classroom by gender, now fed by computed statistics instead of survey columns. For a weighted matrix there is also a more complete view than a node-link diagram: the heatmap shows every cell, including the empty ones, with no layout algorithm in between.

```{r}
#| label: icews-heatmap
#| fig-height: 6
#| fig-cap: "The same network as tiles. Weighted networks often read better this way: the asymmetries, who sends more than they receive, are visible cell by cell."
plot(net14_log, style = "heatmap")
```

The actor statistics themselves deserve a figure too: ranking tables do not survive a slide:

```{r}
#| label: actor-stats-plot
#| fig-height: 4
#| fig-cap: "Out- and in-strength by actor: the ranking table above as a picture, and the first look at the sender/receiver asymmetry that Day 9 turns into parameters."
plot_actor_stats(sa14, across_actor = TRUE,
                 specific_stats = c("strength_sum_out", "strength_sum_in"))
```

And `summary_actor()` is extensible: hand any function of the adjacency matrix to `other_stats` and it becomes a column beside the built-ins. Here is the pattern:

```{r}
#| label: actor-stats-custom
sa_custom <- summary_actor(net14, other_stats = list(
  n_partners_out = function(mat) rowSums(mat > 0, na.rm = TRUE)))
head(sa_custom[, c("actor", "degree_total", "n_partners_out")], 3)
```

Remember it when the centrality you need is not on the menu.

::: {.callout-note collapse="true"}
## Depth: plotting is its own craft, and `netify` gives you the seams

`plot()` returns a `ggplot` object, so everything you know about `ggplot2` applies: you can add layers, change scales, facet a longitudinal network by year. Under the hood `netify` computes a layout and hands it back through `get_node_layout()` / `get_edge_layout()` if you want to freeze coordinates across time so a longitudinal network does not jump around frame to frame. For a static teaching plot the defaults are fine; for a paper figure you will want that control.

When the high-level plot method almost does what you want, `net_plot_data()` exposes the two data frames behind it. The edge frame contains `x1`, `y1`, `x2`, and `y2`; the node frame contains `x` and `y`. From there it is ordinary `ggplot2`:

```{r}
#| label: toolbox-manual-plot
#| fig-height: 5
#| fig-cap: "The 2014 cooperation network assembled from the edge and node data frames that sit underneath plot.netify()."
plot_parts <- net_plot_data(net14_log)$net_dfs

ggplot() +
  geom_segment(
    data = plot_parts$edge_data,
    aes(x = x1, y = y1, xend = x2, yend = y2,
        alpha = log_verbCoop),
    colour = "grey45", linewidth = 0.35
  ) +
  geom_point(
    data = plot_parts$nodal_data,
    aes(x = x, y = y, size = strength_sum_out),
    colour = "#18453B"
  ) +
  scale_alpha_continuous(range = c(0.05, 0.65), guide = "none") +
  scale_size_continuous(range = c(2, 9), name = "out-strength") +
  coord_equal() +
  theme_void()
```

For a longitudinal figure, compute one layout and reuse it across waves, or use `static_actor_positions = TRUE` in the high-level plot method. Allowing a new layout in every facet can make stable actors look as though they moved. The `manual_plotting` vignette in the netify documentation develops this seam layer further.
:::

::: {.callout-note collapse="true"}
## The full toolbox: styles, themes, ego networks, layers (self-study)

**Styles.** Every plot above used the defaults. `netify` ships preset styles, each supplied as a function that you pass to `style =`, as well as named palettes:

```{r}
#| label: toolbox-styles
head(list_network_styles(), 4)
head(list_palettes(), 4)
```

```{r}
#| label: toolbox-style-tufte
#| fig-cap: "The classroom again, in the Tufte preset."
plot(class_net, style = style_tufte)
```

`plot_with_style()` and the `theme_netify()` / `theme_publication_netify()` family cover the rest of the aesthetics, and `animate_netify()` renders a longitudinal object as an animation. It needs `gganimate` and additional render time, so it is shown rather than run:

```{r}
#| label: toolbox-animate
#| eval: false
animate_netify(net_long_conf, static_actor_positions = TRUE)
```

**Ego networks.** `ego_netify()` carves out one actor's neighborhood, and the plot method has layouts built for exactly that shape:

```{r}
#| label: toolbox-ego
#| fig-cap: "The United States' 2014 cooperation neighborhood, radial layout: ego at the centre, alters arranged around it."
usa_ego <- ego_netify(net14, ego = "United States")
plot(usa_ego, layout = "radial")
```

An ego network contains the focal actor, its directly connected alters, and the alter-to-alter ties observed inside that boundary. Every statistic is therefore local to the boundary. A high betweenness score here means central inside this extracted neighborhood, not globally influential in the parent network.

The boundary is part of the analysis, not a plotting preference. In a weighted network it changes when the tie threshold or neighborhood direction changes:

```{r}
#| label: toolbox-ego-sensitivity
w14 <- as.matrix(net14)
strong_cut <- quantile(w14[w14 > 0], 0.90, na.rm = TRUE)
usa_out <- ego_netify(net14, ego = "United States",
                      ngbd_direction = "out")
usa_out_strong <- ego_netify(net14, ego = "United States",
                             ngbd_direction = "out",
                             threshold = strong_cut)
c(any_positive_outgoing = n_actors(usa_out),
  strongest_outgoing = n_actors(usa_out_strong))
```

Report the direction, threshold, time window, and parent network whenever you compare egos. `remove_ego_edges()` is useful when you want a picture that reveals alter-to-alter ties, but removing those edges changes the graph and should not silently change the object used for summaries. One last distinction matters: carving an ego network out of a complete network is not the same data design as collecting an egocentric survey. The latter may not observe all alter-to-alter ties, so full-network measures do not become valid just because the result can be plotted as a network.

**Multilayer.** Two relations on the same actors can live in one object, with cooperation and conflict side by side and `n_layers()` now returning 2:

```{r}
#| label: toolbox-multilayer
net14_conf <- netify(ic14, actor1 = "i", actor2 = "j",
                     symmetric = FALSE, weight = "matlConf")
ml <- layer_netify(list(net14, net14_conf),
                   layer_labels = c("coop", "conflict"))
ml
```

The `ego_networks` and `multilayer_networks` vignettes take both to full depth; no later day fits models to these shapes, which is why they live in this self-study section.
:::

# Handing the object to a model {#sec-export}

The benefit of doing this through one object rather than ad hoc `reshape` calls is portability. Every modeling package in this course wants the network in *its own* shape: `statnet` and `ergm` want a `network` object, `igraph` wants an `igraph`, and the latent-variable models in `lame` want a stack of matrices. You built the network once with the covariates aligned. When you export it, that alignment comes along.

One prerequisite move first. Two of the destinations, `ergm` on Day 11 and the binary family of models on Day 10, want a 0/1 network, and ours carries counts. `binarize()` makes the weighted-to-binary decision explicit instead of letting it hide in an export default:

```{r}
#| label: export-binarize
w <- as.matrix(net14)
net14_bin <- binarize(net14, threshold = median(w[w > 0], na.rm = TRUE))
c(density_weighted = summary(net14)$density,
  density_binary   = summary(net14_bin)$density)
```

Thresholding at the median positive weight turns "essentially everyone signals everyone" (density `r round(summary(net14)$density, 2)`) into a network of *above-typical* cooperation (density `r round(summary(net14_bin)$density, 2)`) that a binary model can actually learn from. The threshold is *your* claim: zero, the median, a substantive cutoff. Saying it out loud is the point.

Now the doors. Each export is one line, and the covariates ride along in whatever form the destination expects:

```{r}
#| label: export-statnet
net_stat <- to_statnet(net_cov)
net_stat
network::list.vertex.attributes(net_stat)   # The nodal covariates rode along
```

```{r}
#| label: export-igraph
net_ig <- to_igraph(net_cov)
# ... and igraph's entire toolkit is immediately usable. The three most
# eigen-central states, straight off the exported object:
tail(sort(round(igraph::eigen_centrality(net_ig)$vector, 2)), 3)
```

`to_statnet()` attaches nodal covariates as vertex attributes and dyadic ones as network-level matrices ready for `ergm`'s `edgecov()`; `to_igraph()` hands you every centrality and community routine igraph owns, no re-wrangling. The third door is `to_lame()`, and Day 10 lives behind it, so today we only read the sign on it: the export splits the covariates into the sender / receiver / dyad arrays `ame()` needs and even writes the model call you will run:

```{r}
#| label: export-lame-call
to_lame(net_cov)$ame_call
```

And nothing is trapped in the object. `unnetify()` closes the loop back to the dyadic data.frame we started from, covariates and all: the shape you would want for a merge, a regression table, or a colleague who has never heard of any of this:

```{r}
#| label: export-unnetify
head(unnetify(net_cov), 3)
```

**The network object is the single source of truth, and the modeling packages are just different views of it.** This is the entire reason `netify` exists, and it is why Day 8 is about the object: get it right once, and days 9 through 14 are downstream of it.

::: {.callout-note collapse="true"}
## The full toolbox: three lines for the tidyverse (self-study)

The broom trio works on netify objects too (this needs `dplyr` installed): `as_tibble()` and `tidy()` give the tidy edgelist, `glance()` the one-row network summary: handy inside a `purrr` pipeline over many networks.

```{r}
#| label: toolbox-broom
dplyr::as_tibble(net14)[1:3, ]
generics::glance(net14)[, 1:6]
```
:::

::: {.brk}
Break. [The object is built. After the break we ask the only question that matters: now that you have a network, what do you *do* with it? That question is the rest of the course.]{.sub}
:::

# Now what do you DO with it? {#sec-map}

You have spent the descriptive week that precedes this course, and most of today, *describing* networks: degree, closeness, betweenness, communities, and just now `summary_actor()` on the trade data we are about to load. Description is useful and it is where everyone starts. But it has a ceiling, and the rest of this course lives above that ceiling. The next example shows where that ceiling is, followed by a map of the models that take us beyond it.

We switch to a cross-sectional directed trade network of thirty countries because it gives us a clean export example with meaningful source-target direction and dyadic covariates. Day 9 uses different applications, but the object-shaping work is the same.

```{r}
#| label: trade-load
# The .rda ships in this day's data/ folder alongside several day-9
# objects; load it into its own environment and take only the one
# data.frame we need, so nothing else lands in your workspace unannounced
e <- new.env()
load(file.path(day8_dir, "data", "trade_netlet.rda"), envir = e)
trade <- e$trade
trade_net <- netify(
  trade, actor1 = "Var1", actor2 = "Var2",
  symmetric  = FALSE, weight = "trade",
  nodal_vars = c("pop1", "gdp1", "polity1"),
  dyad_vars  = c("conflicts", "distance", "shared_igos"),
  dyad_vars_symmetric = c(FALSE, TRUE, TRUE)
)
sa_trade <- summary_actor(trade_net)
head(sa_trade[order(-sa_trade$strength_sum_out),
              c("actor", "strength_sum_out", "strength_sum_in")])
```

## The ceiling of description {#sec-ceiling}

Description tells you *who* trades the most, who sits at the centre, and who clusters with whom. It is necessary, and you should always do it. But the moment you want to say something *causal or explanatory*, such as "Do democracies trade more with each other?" or "Does sharing an alliance raise trade?", three questions arrive that description simply cannot answer:

- **Why?** Is the pattern homophily, reciprocity, transitivity, or a covariate? A centrality ranking cannot tell you.
- **How sure?** Description hands you a number with no uncertainty around it. No interval, no test.
- **What next?** Description cannot predict a tie it has not seen.

`netify` itself hands you the first rung above that ceiling, and it is worth climbing it right now so you can feel exactly where it ends. The question we just posed, whether democracies trade more with one another, is a **homophily** question, and `homophily()` answers its *descriptive* half with a simple dyad-permutation benchmark. It holds the observed ties fixed, computes attribute similarity for every eligible dyad, shuffles those similarity values across dyads, and asks whether connected pairs look more similar than that shuffled reference would suggest. This treats dyads as exchangeable, so it is an exploratory benchmark, not a node-label permutation, a correction for network dependence, or a causal test.

```{r}
#| label: ceiling-homophily
set.seed(6886)   # homophily()'s permutations draw from the RNG
h <- homophily(trade_net, attribute = "polity1")
h[, c("homophily_correlation", "similarity_difference", "p_value")]
```

The observed regime-similarity correlation is a modest `r round(h$homophily_correlation, 3)`, with a permutation p-value of `r round(h$p_value, 2)`: under this crude dyad-exchangeability reference, we do not detect regime homophily. That is a descriptive result under one simple reference distribution, not a covariate-adjusted answer.

The same permutation logic also answers "is this network's structure surprising at all?": `compare_to_null()` simulates networks from a null model at the observed density and asks where your statistic falls among them:

```{r}
#| label: ceiling-null
null_recip <- compare_to_null(trade_net,
  fn    = function(net) c(reciprocity = summary(net)$reciprocity),
  n_sim = 200, model = "erdos_renyi", seed = 6886, verbose = FALSE)
null_recip[, c("metric", "observed", "null_mean", "p_value", "extreme")]
```

Observed reciprocity of `r round(null_recip$observed, 2)` against a null distribution centred at `r round(null_recip$null_mean, 2)`: p = `r round(null_recip$p_value, 3)`, about as extreme as two hundred draws can call anything. Trade flows are highly reciprocal under this reference distribution. Day 9's social relations model gives reciprocity a parameter on a different application, while Days 11 and 13 reuse the broader logic of comparing what we observed with an explicit reference world.

But look hard at what both of these tests have in common: **one statistic at a time, nothing held constant, no prediction.** A permutation p on homophily cannot tell you whether democracies trade because they are democracies or because they are rich and near each other: it has no way to hold anything else fixed. That is the ceiling. The instant you want covariates adjusted for one another, you reach for a regression. And that is where the one problem of the whole course walks in.

::: {.callout-note collapse="true"}
## The full toolbox: attribute analysis and null-model primitives (self-study)

**Attribute analysis.** `homophily()` has siblings. `mixing_matrix()` is the categorical version: who ties with whom, across groups (back to the classroom, where gender is categorical):

```{r}
#| label: toolbox-mixing
#| fig-width: 5
#| fig-height: 4
mm <- mixing_matrix(class_net, attribute = "gender")
mm$summary_stats[, c("attribute", "assortativity", "diagonal_proportion")]
plot_mixing_matrix(mm)
```

`dyad_correlation()` screens dyadic covariates against the outcome network, previewing Day 10's `Xdyad` logic. `plot_homophily()` draws the distribution behind the permutation test you ran above:

```{r}
#| label: toolbox-dyadcorr
dyad_correlation(trade_net, dyad_vars = "distance")[
  , c("dyad_var", "correlation", "p_value")]
```

```{r}
#| label: toolbox-plot-homophily
#| fig-height: 4
#| fig-cap: "Polity similarity for connected versus unconnected pairs, the picture behind the permutation test, and visibly why it came back null."
plot_homophily(h, netlet = trade_net, type = "distribution",
               attribute = "polity1")
```

`attribute_report()` runs the whole battery, including homophily, mixing, dyadic and centrality correlations, and attribute summaries, in one call. Its output is long, so here is just its one-line verdict; **remember this function when Homework 3 asks for attribute EDA**:

```{r}
#| label: toolbox-attribute-report
ar <- attribute_report(trade_net, node_vars = "polity1",
                       dyad_vars = "distance")
ar$overall_summary$summary_points
```

**How to read this battery without over-reading it.** These functions are useful screens, but none of them is a jointly adjusted network model:

- For a continuous attribute, `homophily(method = "correlation")` uses the negative absolute difference between the two actors as its similarity score; the reported homophily correlation relates that score to tie presence. For a categorical attribute, similarity is simply same category versus different category.
- A mixing matrix shows where observed ties fall across groups. Raw counts are partly driven by how many dyads each group pairing makes possible, while row normalization answers the directional question "where do ties from this group go?" A large diagonal is not automatically a mechanism or a causal effect.
- `dyad_correlation()` is a bivariate screen. Its usual p-value comes from an ordinary correlation test and does not adjust for dyads sharing actors. It is especially misleading if the proposed predictor is the outcome itself, or a component used to construct the outcome.
- `attribute_report()` organizes these diagnostics in one place. It does not make them conditional on one another, and it does not turn centrality into influence.

For longitudinal or multilayer objects, the functions repeat the descriptive calculation by wave or layer. That is useful for finding patterns to investigate, but it is not a model of change over time or dependence across layers.

**Null-model primitives.** `compare_to_null()` is a convenience wrapper, and the parts are yours to use directly. `simulate()` draws whole networks from a null family: `"erdos_renyi"` (density only), `"configuration"` (preserves the degree sequence), or `"dyad_permutation"` (relabels actors and preserves everything else). `bootstrap_netlet()` puts an interval on any statistic of the observed network:

```{r}
#| label: toolbox-simulate
sims <- simulate(net14, nsim = 3, model = "configuration", seed = 6886)
summary(sims[[1]])$density   # Same degrees as net14, shuffled structure
```

```{r}
#| label: toolbox-bootstrap
bootstrap_netlet(trade_net,
  fn = function(net) c(reciprocity = summary(net)$reciprocity),
  n_boot = 200, seed = 6886, verbose = FALSE)
```

Choosing *which* null to simulate from is itself the modeling decision. "Surprising given density" and "surprising given the degree distribution" are different questions, and Day 11 makes that choice into a career.
:::

## The one problem {#sec-problem}

The moment you want to *explain* a tie with a covariate, you write a regression with one row per dyad. And a regression assumes the rows are independent. **In a network they are not**, because the same actor sits in a whole band of rows at once. Let me make that concrete before I prove it hurts.

```{r}
#| label: coupling
n <- 30                        # Actors in the trade network
dyad_rows      <- n * (n - 1)  # Rows in a directed dyadic regression
rows_per_actor <- 2 * (n - 1)  # Rows that contain any one given actor
c(dyad_rows = dyad_rows,
  rows_per_actor = rows_per_actor,
  share = round(rows_per_actor / dyad_rows, 3))
```

`r dyad_rows` rows *look* like `r dyad_rows` independent facts. But nudge one country, perhaps because Japan trades more than its size predicts, and `r rows_per_actor` of those rows, about `r round(100 * rows_per_actor / dyad_rows)`% of the entire regression, move together, because Japan is the sender or the receiver in every one of them. The rows are not independent draws. They are coupled in bands, one band per actor.

Does that coupling actually break anything? Let's not argue about it: let's simulate a world where we *know* the truth and watch the regression lie.

```{r}
#| label: p-value-lies
set.seed(6886)
n <- 25

# One experiment: a predictor x with NO true effect on y.
# The common actor-level case, x is an ACTOR attribute (gdp, regime,
# ideal point) copied across every dyad that actor is in.
sim_once <- function(coupled_x = TRUE) {
  a <- rnorm(n)      # Actor "sociability" -> structures the ERRORS
  z <- rnorm(n)      # An actor attribute with ZERO true effect
  d <- expand.grid(i = 1:n, j = 1:n)
  d <- d[d$i != d$j, ]
  d$x <- if (coupled_x) z[d$i] else rnorm(nrow(d))  # z[d$i]: actor i's single value STAMPED onto every dyad i sends, this is the coupling
  d$y <- a[d$i] + a[d$j] + rnorm(nrow(d))  # y depends on WHO is in it
  coef(summary(lm(y ~ x, data = d)))["x", "Pr(>|t|)"]
}

p_actor <- replicate(500, sim_once(coupled_x = TRUE))
p_dyad  <- replicate(500, sim_once(coupled_x = FALSE))

c(actor_attribute = mean(p_actor < 0.05),
  dyad_specific   = mean(p_dyad  < 0.05))
```

Read that carefully, because it is the single most important number today. We built a world where `x` has **no** effect on `y` whatsoever. A well-calibrated test should call it "significant" about 5% of the time. When the predictor is an actor attribute, which GDP, population, democracy, and ideal point all are, naive OLS calls it significant **`r round(100 * mean(p_actor < 0.05))`% of the time**. That is a coin flip wearing a lab coat. When the predictor varies dyad by dyad, the false positive rate is a well-behaved `r round(100 * mean(p_dyad < 0.05))`%. The damage comes from actor-level coupling, the same coupling the arithmetic above described. One line of the simulation creates it: `d$x <- z[d$i]` takes an actor-length vector `z` and indexes it by the dyad's sender column, stamping each actor's single value across every dyad it appears in. That is the "copied into every dyad that country is in" move you made by hand back in @sec-covariates, now shown to be the source of the false result.

::: {.callout-important}
## The Checkpoint to Carry Forward

Before you model anything of your own, ask: **Is your predictor an actor attribute or a dyadic attribute?** If it is an actor attribute, as most of the interesting ones are, your *effective* sample size is closer to your number of actors than your number of dyads (for the 30-country trade network you built today, that is about 30, not its 870 directed dyads or the 25-actor toy the simulation above used), and a naive p-value is not to be trusted. This is re-runnable and it stands alone; it is the one thing to remember if you remember nothing else from today.
:::

::: {.callout-note collapse="true"}
## A first taste of one fix, previewing Day 13

You do not have to change the model to make the p-value stop lying: one option is to change only how the standard errors are computed. Cluster them on the sender:

```{r}
#| label: cluster-preview
suppressMessages({ library(sandwich); library(lmtest) })
set.seed(6886); n <- 25
fp_clustered <- mean(replicate(500, {
  a <- rnorm(n); z <- rnorm(n)
  d <- expand.grid(i = 1:n, j = 1:n); d <- d[d$i != d$j, ]
  d$x <- z[d$i]; d$y <- a[d$i] + a[d$j] + rnorm(nrow(d))
  m <- lm(y ~ x, data = d)
  coeftest(m, vcov = vcovCL(m, cluster = ~ i))["x", "Pr(>|t|)"]
}) < 0.05)
fp_clustered
```

Clustering on the sender pulls the false-positive rate from ~50% down to about `r round(100 * fp_clustered)`%: *without changing a single coefficient*, just by telling the standard errors that rows sharing a sender are not independent. It still overshoots 5% because with only ~25 actors there are too few clusters for the robust standard errors to be exactly calibrated. Day 13 treats that small-network limitation directly and develops the full dyadic cluster-robust calculation.
:::

## The map {#sec-mapfig}

The dependence will not go away because **the dependence *is* the network**: it is the thing you are trying to study, not noise to be scrubbed out. You do not get to ignore it. You get to *choose where to put it*. There are only a handful of places, and each one is a method, and each method is a day.

The cleanest way to see it is one equation. Write $\eta_{ij}$ for the conditional linear predictor of a tie, so $\mathbb E(Y_{ij}\mid\text{covariates, latent terms})=g^{-1}(\eta_{ij})$. Every day for the rest of the week adds, replaces, or refuses one piece:

$$
\eta_{ij} \;=\; \underbrace{\beta^{\top} x_{ij}}_{\text{observed covariates}}
\;+\; \underbrace{\,a_i + b_j\,}_{\substack{\text{additive}\\ \text{SRM · Day 9}}}
\;+\; \underbrace{\;\alpha_{z_i z_j}\;}_{\substack{\text{discrete roles}\\ \text{blockmodel · Day 9}}}
\;+\; \underbrace{-\lVert z_i - z_j\rVert}_{\substack{\text{shared position}\\ \text{distance · Day 10}}}
\;+\; \underbrace{\,u_i^{\top}\Lambda\, v_j\,}_{\substack{\text{sender/receiver factors}\\ \text{factor/AME · Day 10}}}
$$

This is a map, not a recommendation to put every latent term in one kitchen-sink model. The block, distance, and factor terms are usually alternative ways to represent the same residual structure. For binary ties, $g^{-1}$ is commonly the inverse logit or Probit; for Gaussian outcomes, the residual error is added on the outcome scale.

The two geometry terms are **not** the same move, and Day 10 is built around exactly that contrast. The distance idea gives every actor a *single* latent position $z_i$ and a tie forms when two actors sit close: a symmetric distance, $-\lVert z_i - z_j\rVert$, so $i\to j$ and $j\to i$ get the same pull. The factor model instead gives each actor a *separate* sender vector $u_i$ and receiver vector $v_j$, so the bilinear term $u_i^{\top}\Lambda v_j$ can carry directed, asymmetric structure that one shared position cannot. Day 10 opens with the distance intuition, the cleanest picture of latent geometry there is, and then shows, on data where we know the truth, exactly where it breaks and why the factor model contains it and goes further.

Those four added terms are **one philosophy**: represent dependence with latent actor structure in the conditional mean or link. Conditional on those latent effects, distinct dyads may be independent; after we average over them, dyads that share actors are dependent. This can reduce distortion in the covariate estimates, *provided* the latent effects are not themselves correlated with the covariates. When an actor trait like GDP is tangled up with how much that actor sends, the random-effect route can bias the very coefficient you care about, and that is where the robust-inference camp pushes back. Three other methods take a different route:

- **ERGM (Day 11)** does not decompose the variance. It writes the probability of the *whole graph* directly, $p(Y\mid\theta)\propto\exp\{\theta^{\top} g(Y)\}$, where $g(Y)$ counts edges, triangles, shared partners. Dependence goes **into the model of the structure itself**. Day 11 ends by carrying that same machinery across panel waves: the temporal ERGM.
- **SAOM (Day 12)** treats the graph as the moving part and models $P(Y^{(t+1)}\mid Y^{(t)})$ as actors add and drop ties. Dependence goes **across time**.
- **Dyadic cluster-robust inference (Day 13)** keeps the naive $y_{ij}=\beta^{\top}x_{ij}+\varepsilon_{ij}$, with no extra term at all, and changes only how $\mathrm{Var}(\hat\beta)$ is computed. Dependence goes **into the standard errors** (you just saw the four-line preview). Day 13 builds this out into a full session on causal inference with network data: what randomization buys you, what spillovers do to your estimand, and which claims an observational design can support.

And then **Day 14** refuses to pick a cell at all. It puts the models head to head on one dataset, offers a decision framework for *which question each one actually answers*, and tours what is out there beyond this course (relational event models, multilayer networks).

The whole map fits in one table. One column, the *measurement role*, gets explained immediately afterward:

```{r}
#| label: the-map
#| echo: false
map <- data.frame(
  `where the dependence goes` = c(
    "describe it, don't model",
    "into LATENT EFFECTS (additive)",
    "into LATENT EFFECTS (discrete)",
    "into LATENT EFFECTS (distance)",
    "into LATENT EFFECTS (factor)",
    "into the MODEL (structure)",
    "across TIME",
    "into the STANDARD ERRORS",
    "all of the above, head-to-head"),
  method = c(
    "centrality, communities",
    "social relations model",
    "stochastic blockmodel",
    "latent distance model",
    "latent factor / AME",
    "ERGM",
    "SAOM (RSiena)",
    "dyadic cluster-robust SEs + causal design",
    "which model when? + the frontier"),
  `the one sentence` = c(
    "who is central, who clusters, but not why, no uncertainty, no prediction",
    "some actors just send/receive more; reciprocity, soak it into random effects",
    "actors fall into blocks; ties depend only on block membership",
    "actors near each other in a latent space tie, homophily as distance",
    "actors with similar latent roles tie, captures triangles description missed",
    "a coefficient on the triangle asks whether there is closure beyond nodal traits",
    "model how ties are added and dropped wave to wave",
    "keep the naive model, fix the inference, and ask what causal claim survives",
    "one dataset, every model, and a framework for choosing among them"),
  `the measurement role` = c(
    "not applicable",
    "who over-sends, who over-receives",
    "recovered roles / blocs",
    "the map itself: latent positions",
    "latent relations: who leans toward whom",
    "not applicable",
    "not applicable",
    "not applicable",
    "not applicable"),
  day = c("Days 1-8", "9", "9", "10", "10", "11", "12", "13", "14"),
  check.names = FALSE)
knitr::kable(map)
```

Every later day fills in one or two rows of this table and points back at it. Day 9 takes the additive and discrete latent-effect rows, Day 10 takes the two continuous latent-effect rows, Days 11 through 13 take one row each, and Day 14 runs the whole table on a single dataset. This sequence continues the foundation you built in Days 1 through 7 rather than starting a separate course.

## The latent terms: nuisance, or the thing itself {#sec-measurement}

One column needs more explanation: **the measurement role**. Everything above treats the latent terms as a *nuisance*: the $a_i$'s, the blocks, the $z_i$'s and $u_i$'s give the modeled dependence somewhere to live while we estimate the covariate associations. Days 9 and 10 often use them that way. There is also a second use that I find especially interesting: sometimes the latent structure **is the estimand**. You did not come for $\beta$ at all; you came for the map.

A useful analogy comes from test theory. Psychometricians call it *item response theory*, or IRT, and the acronym will come back in a moment. No single exam question measures ability, but forty noisy questions together help locate a student on an unobserved trait. Ideal-point models use the same logic: nobody directly observes ideology, yet hundreds of roll-call votes help place each legislator on a recovered dimension. Networks have a close analogue. Every tie is one noisy observation about a pair, and the positions the model recovers can summarize relational features we do not directly measure, such as affinity or alignment. *Who has moved closer to whom* in trade, votes, or conflict can be the finding itself, not merely a step toward a coefficient.

So why is this not just IRT with more rows? The inference is harder for three reasons:

- **The map is not identified on its own.** Rotate, reflect, or slide every position at once and the likelihood does not move. In Day 10's factor variant, even the overall scale can be traded against the loadings. Only *relative* geometry means anything, so positions are interpretable only after Procrustes alignment, and the number of dimensions is chosen a priori, not estimated: Day 10 teaches both fixes.
- **One graph is one realization.** IRT gets many students crossed with many items; you get one network. $n(n-1)$ cells, but coupled in the bands from @sec-problem, with no repeated items per pair: information about any one position accrues slowly, and "large-$n$" theory on a single graph remains nonstandard. Day 10's uncertainty statements will therefore be more careful than you may be used to.
- **The space absorbs whatever the mean model omits.** Leave a relevant covariate out and it does not disappear: it can reappear *as latent structure*, so your "map of affinity" partly reflects the specification you chose. The same problem runs in the other direction. If the latent terms correlate with covariates you *did* include, the exogeneity argument for $\beta$ fails too, with bias on the order of $\mathrm{Cov}(x,\text{latent})/ \mathrm{Var}(x)$. Day 13 takes this up as a design problem.

The first reason fits in four lines, so let's not take it on faith:

```{r}
#| label: rotation-invariance
set.seed(6886)
z  <- matrix(rnorm(10 * 2), 10, 2)   # Ten actors in a 2-d latent space
th <- pi / 3                         # Rotate the whole map by 60 degrees
Rm <- matrix(c(cos(th), -sin(th), sin(th), cos(th)), 2, 2)
rot_gap <- max(abs(dist(z) - dist(z %*% Rm)))
rot_gap                              # Distances are all the model sees
```

`r format(signif(rot_gap, 2))`: machine zero. Two maps that look completely different on paper are *identical* to the likelihood, and a reflection does the same trick exactly. Any substantive story you tell about "the x-axis" of a latent space is a story about an arbitrary rotation until you have pinned it down.

This map also has a rival cartography that often appears in econometrics. Econometricians hit the same dependence and made the opposite bet: put it in **fixed effects**, not random ones. Condition the actor-level heterogeneity out entirely (Graham 2017, *Econometrica*: the tetrad logit for network formation, which treats each actor's degree heterogeneity as incidental parameters and conditions on a sufficient statistic for them), bias-correct the incidental-parameter contamination instead (Fernández-Val & Weidner 2018, *Annual Review of Economics*), or absorb it in two-way exporter/importer terms and fit the gravity model by PPML: Santos Silva and Tenreyro (2006) showed the logged version breaks under heteroskedasticity, and Jochmans (2017) showed that fixed-effect bias in these dyadic models is the same order as the standard error, so naive intervals under-cover. Then use dyadic-robust variance estimates (Cameron and Miller 2014; Tabord-Meehan 2019 supplies the published asymptotics), which is the estimator introduced on Day 13.

The sharpest point of contact is Bai (2009, *Econometrica*): his *interactive fixed effects* term $\lambda_i^{\top} f_t$ is the same bilinear object as AME's $u_i^{\top}\Lambda v_j$, estimated as fixed parameters instead of given a prior. The fixed-effects approach never has to assume the heterogeneity is uncorrelated with your covariates, a practical advantage that we will take seriously, at the price that the latent relations are conditioned away rather than measured. So the two camps split cleanly: fixed effects plus robust standard errors is robust to *correlated* heterogeneity and silent about the relations themselves; random or latent structure is efficient and produces measures of the latent relations, but owes an exogeneity assumption in return. Neither dominates: the choice tracks exactly the question this section just posed: is the latent structure your nuisance, or your estimand? Day 14's reading list opens the doorway into this literature (de Paula 2017 is the survey to start with).

Keep the one-line version: when the latent structure is your *nuisance* you mostly need it to stay uncorrelated with your covariates; when it is your *estimand* you additionally need it identified and its uncertainty quantified. Those are different demands, and Day 10 is where we start paying both bills.

## The Model-Choice Tension {#sec-tension}

Two cells of that map, "dependence in the model" (ERGM) and "dependence in the variance" (latent space / AME), are not just different tools; they answer *different questions*, and the field debates how to choose between them. I wrote the AME paper, so I want to state the distinction directly rather than sell you my side: **ERGM answers a question about mechanism that AME cannot; AME answers a question about covariate inference and prediction after modeling latent dependence. Pick the one that matches your question, not the one that "wins."** The full version, including ERGM's old degeneracy reputation and what modern `ergm` 4.x does about it, is in the box. Day 14 stages the comparison on one dataset.

::: {.callout-note collapse="true"}
## ERGM vs. Latent Space: The Substantive Disagreement

Two cells of that map, "dependence in the model" (ERGM) and "dependence in the variance" (latent space / AME), are not just different tools. They answer *different questions*, and the field debates how to choose between them. I wrote the AME paper, so I want to state the tradeoff clearly rather than sell you my side:

- **ERGM makes dependence the substance.** If your *theory is about structure*, such as "conflict spreads by closing triangles," "friendship reciprocates," "popularity begets popularity": ERGM hands you a coefficient on exactly that mechanism and a hypothesis test for it. A latent-variable model *cannot* give you a "transitivity effect" to interpret; it absorbs the triangle-generating structure into geometry, and you can fit it beautifully while never learning the mechanism was closure. That is a genuine loss, and it is ERGM's home turf.
- **Latent space / AME makes dependence a nuisance to represent**, so the covariate effects are estimated after accounting for the modeled latent structure and the model can predict held-out ties. It is rarely degenerate and scales comfortably to large graphs. But it will not name the structural mechanism for you.

The short version is: **ERGM answers a question about mechanism that AME cannot; AME answers a question about covariate inference and prediction after modeling latent dependence. They optimize different things. Pick the one that matches your question, not the one that "wins."** Naive triangle-based ERGMs can be prone to degeneracy, while geometrically weighted terms such as `gwesp` and `gwdsp` and modern `ergm` estimation make many specifications much more stable, though they still require care. ERGMs also remain harder to fit at scale than latent-variable models. Days 10 and 11 let you work through that tradeoff on the same data, and Day 14 puts the two side by side.
:::

## Tomorrow: the social relations model, then blocks {#sec-tomorrow}

The dependence is part of the network, not noise, and you have to decide where to put it. The cheapest place to start, and the one whose bookkeeping every later method reuses, is the variance. We say some actors just send more and some receive more, some pairs reciprocate, and we let *random effects* carry all of that so the covariates can return to their conditional role, as long as those random effects are not themselves entangled with the covariates, the caveat we just met. That is the social relations model, and it is the first half of tomorrow. The second half takes the same variance philosophy in a discrete direction: actors sorted into *blocks*, with ties depending only on which block they occupy. That is the stochastic blockmodel. One day, two ways of putting the dependence in latent actor effects: continuous actor effects before the break, discrete roles after it.

And here is the object we start from: the exact trade network you built two sections ago, handed to the model with one line:

```{r}
#| label: handoff
lame_data <- to_lame(trade_net)   # Day 10 uses the same matrix-list contract
sapply(lame_data[c("Y", "Xrow", "Xcol", "Xdyad")],
       function(x) paste(dim(x), collapse = " x "))
```

Tomorrow's very first question is the one this half of today just made unavoidable: *What exactly goes wrong if I just run the regression?* By the end of Day 9 one coefficient will have flipped sign, with the standard errors left untouched. That sign flip is the p-value lie you saw today, caught and corrected by putting the dependence where it belongs.

## After the Code Runs: Say What Network You Built {#sec-after-fit}

Today does not end with a model coefficient. It ends with a defensible measurement statement. Before interpreting density, centrality, reciprocity, or any other network summary, write down four things in ordinary political science language:

1. **Who can appear?** Name the actors and the eligibility rule. A state-year network, for example, should not quietly treat states that do not yet exist as eligible isolates.
2. **What makes a tie?** State the event, threshold, direction, weight, and time window. "Conflict" is too vague if the code actually records at least one material-conflict event initiated during a calendar year.
3. **What do zero and missing mean?** A zero is an observed eligible non-tie. Missing means the tie was not observed. A structural zero means the tie was impossible by design.
4. **What descriptive pattern deserves explanation?** Name the politically interesting concentration, asymmetry, clustering, or change over time, but do not turn that pattern into a mechanism before fitting a model that represents the relevant alternatives.

::: {.callout-tip title="A human version of the result"}
"We construct an annual directed network of material conflict among states. A tie records at least one material-conflict event initiated by one state against another during the year. State pairs that do not coexist are excluded rather than coded as peaceful dyads. The resulting network is sparse, and a small number of states initiate a large share of the recorded conflict. That concentration motivates separate model terms for how often each state initiates conflict and how often each state is targeted."
:::

That paragraph does useful substantive work. It identifies the political object, tells the reader what an absence means, reports the main descriptive pattern, and explains why the next model needs actor-level heterogeneity. It does not claim that a high-centrality state is powerful, that reciprocity proves retaliation, or that a descriptive association is causal. Those interpretations require additional theory, comparison, and evidence.

# Your turn {#sec-yourturn}

::: {.callout-tip icon=false}
## Exercises

The first two are the core workflow and everyone should do them. The third is a substantial stretch that supports part two's argument. The fourth is on your own data. That is the one that makes the workflow portable.

**1. (Required: the build loop.)** Take the `icews` subset `ic` from @sec-events, but build a network on `matlConf` (material *conflict*) instead of `verbCoop`, for the year 2010. Is it more or less dense than the cooperation network? Run `summary_actor()` and name the state that initiated the largest number of conflict events. In one sentence, explain why the state that initiates the most conflict need not be the state targeted by the most conflict.

**2. (Required: covariates and export.)** Attach `i_polity2` as a nodal covariate and `matlConf` as a dyadic covariate to your network from (1), then export it three ways (`to_statnet`, `to_igraph`, `to_lame`) and confirm the covariate survived each trip. Why does `to_lame()` split one nodal covariate into a sender column *and* a receiver column?

**3. (Stretch: make the p-value lie, then stop it.)** In the simulation of @sec-problem, change `n` from 25 to 60 and re-run. Does the actor-attribute false-positive rate improve? (It does not: this is not a small-sample problem, and understanding *why* is the point.) Then modify the data-generating process so the null predictor varies on the *dyad* rather than the actor, confirm the false-positive rate returns to ~5%, and write one sentence on why the coupling, rather than the predictor or sample size, was the culprit.

**4. (Your own data.)** Build a `netify` object from a network you actually work with. Decide `symmetric` and `weight` deliberately and say why in a comment. Then answer the checkpoint from @sec-problem: is your key predictor an actor attribute or a dyadic attribute? If it is an actor attribute, write down what your *effective* sample size really is.

```{r}
#| label: solutions
#| code-fold: true
#| code-summary: "Solutions, open after you've tried them"
#| eval: false

## 1. A conflict network
ic10 <- ic[ic$year == 2010, ]
net_conf <- netify(ic10, actor1 = "i", actor2 = "j",
                   symmetric = FALSE, weight = "matlConf")
summary(net_conf)$density        # Far sparser than cooperation:
                                 # most pairs simply have no conflict
sa_conf <- summary_actor(net_conf)
head(sa_conf[order(-sa_conf$strength_sum_out),
             c("actor", "strength_sum_out", "strength_sum_in")])
# Sender vs receiver differ for conflict because initiating conflict
# and being targeted by it are different roles, a state can be a
# frequent target (high in-strength) while rarely initiating (low
# out-strength). Cooperation is far more reciprocal, so its two
# rankings look more alike. That asymmetry is exactly what the
# sender/receiver effects of Day 9 are built to capture.

## 2. Covariates survive the export
net_conf <- add_node_vars(net_conf,
  unique(ic10[, c("i", "i_polity2")]), actor = "i",
  node_vars = "i_polity2")
net_conf <- add_dyad_vars(net_conf, ic10, actor1 = "i", actor2 = "j",
  dyad_vars = "matlConf", dyad_vars_symmetric = FALSE)
to_statnet(net_conf)   # i_polity2 is a vertex attribute
to_igraph(net_conf)    # i_polity2 is a vertex attribute
str(to_lame(net_conf)[c("Xrow", "Xcol")])
# to_lame splits a nodal covariate into Xrow (its value for the
# SENDER of the dyad) and Xcol (its value for the RECEIVER) because
# the model asks two separate questions of the same trait: do
# high-polity states SEND more ties, and do they RECEIVE more? Those
# get two different coefficients, so the covariate has to appear in
# two aligned places. This is the sender/receiver split again.

## 3. n does not save you
sim_once <- function(coupled_x = TRUE, n = 60) {
  a <- rnorm(n); z <- rnorm(n)
  d <- expand.grid(i = 1:n, j = 1:n); d <- d[d$i != d$j, ]
  d$x <- if (coupled_x) z[d$i] else rnorm(nrow(d))
  d$y <- a[d$i] + a[d$j] + rnorm(nrow(d))
  coef(summary(lm(y ~ x, data = d)))["x", "Pr(>|t|)"]
}
set.seed(6886)
mean(replicate(500, sim_once(TRUE,  60)) < 0.05)  # ~0.63
mean(replicate(500, sim_once(FALSE, 60)) < 0.05)  # ~0.05
# Raising n from 25 to 60 more than doubles the number of DYADS but
# the coupling is unchanged: every actor still sits in a whole band
# of rows, so the effective information about an actor-level x still
# grows only like the number of ACTORS. More dyads of the same few
# actors is not more independent evidence. The predictor that varies
# by dyad has independent rows under this simulation, so its test is calibrated at
# any n. The culprit was never the sample size, it was where the
# variation lived.

## 4. Your own data
# my_net <- netify(my_edges, actor1 = "...", actor2 = "...",
#                  symmetric = ..., weight = "...")
# summary(my_net); summary_actor(my_net)
# Effective n for an actor-level covariate ~ number of actors,
# NOT number of dyads.
```
:::

# Reading {#sec-reading}

- Minhas, Dorff, Gallop, Foster, Liu, Tellez & Ward (2022), "Taking Dyads Seriously," *Political Science Research and Methods*, `minhas_etal_2022.pdf`. The paper behind the whole spine of this course. Sections 1 and 2 are exactly the "why can't I just run the regression?" argument you saw in part two; the rest is where we are going.
- Ward, Stovel & Sacks (2011), "Network Analysis and Political Science," *Annual Review of Political Science*: `ward_etal_2011.pdf`. A gentle, discipline-spanning tour of why networks need their own methods. Read this one first if today was your first serious contact with network modeling.

::: {.callout-note collapse="true"}
## Appendix: what we did not show, and where it lives

Today's main path plus the toolbox asides cover the netify surface a working modeler needs. The deliberate omissions, named so you know they are omissions and not oversights:

- **Manual plot assembly**: `net_plot_data()`, `assemble_netify_plot()`, the `netify_node()` / `netify_edge()` / `netify_label()` geom layers, and `get_node_layout()` / `get_edge_layout()` for freezing coordinates across panels; the `manual_plotting` vignette in the netify docs is the walkthrough.
- **Package internals**: `new_netify()` and the `get_adjacency*()` family are the low-level constructors that `netify()` wraps; the `internals` vignette if you are curious.
- **`to_lame()` / `to_amen()` in anger**: Day 10 exercises the whole pipeline; today you only read `$ame_call` off the export.
- **`to_dbn()` and `from_lame_fit()`**: the dynamic-bilinear-network and fitted-model bridges; days 10 and 14 own those.
- **`netify_measurements()`, `get_actor_time_info()`, `melt()`, and the long-form `netify_to_*` aliases**: superseded for our purposes by `n_actors()` / `n_periods()` / `n_layers()`, `subset_netify()`, and the short `to_*` names you learned today. Same story for `add_vertex_attributes()` / `add_edge_attributes()` (straight aliases of the `add_node_vars()` / `add_dyad_vars()` you already know) and the `decompose_igraph()` / `decompose_statnet()` variants of `decompose_netify()` for foreign objects.
- **Modeling ego or multilayer objects**: the objects themselves appeared in the plotting toolbox aside; the `ego_networks` and `multilayer_networks` vignettes go deeper, and no later day fits models to them.

Two habits to leave with. First: whenever you set `symmetric`, `weight`, `sum_dyads`, or `missing_to_zero`, say out loud what claim about the world you are making: the defaults are conveniences, not truths. Second: open `?netify_workflows` once tonight. It is the cheat sheet for this entire course: every pipeline recipe (ERGM, lame, temporal, large networks) in one help page, and it will still be there in August when you cannot remember which day we did what.
:::

```{r}
#| label: session
#| code-fold: true
#| code-summary: "Session info"
sessionInfo()
```
