---
title: "Relational Event Models"
engine: knitr
editor: source
subtitle: "Ancillary self-study · ICPSR Network Analysis: Advanced Topics"
---

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

{{< include _ancillary_header.qmd >}}

::: {.callout-note title="What this extends, and why it is here"}
**Extends:** the temporal ERGM session (Day 11) and the SAOM session (Day 12). Those methods work with networks observed in waves. This document starts from a different kind of record: a time-ordered stream in which the individual interactions are observed.

**The problem it solves:** you have one row per interaction, with a sender, receiver, and timestamp, and the order of those interactions is part of the question. Diplomatic exchanges, protest and repression, legislative speeches and responses, battlefield encounters, calls, and messages can all arrive in this form.

**What you need:** `relevent` and `goldfish`, both available from CRAN. The phone-call data ship with `goldfish`, and the two fitted objects ship in `cache/`.
:::

# 1. Start With the Observed Unit

A panel network records whether a relationship was present at a small number of observation times. A relational event dataset records each interaction as it happens. The outcome is therefore not simply whether $i$ and $j$ ever interact. It is which eligible pair interacts next, when that interaction occurs, or both.

That distinction matters because aggregation can remove the pattern you care about. Ten calls between two students become one tie in a binary snapshot. An immediate return call and a call returned three weeks later both become reciprocity. A relational event model keeps the sequence available to the estimator.

Observing the sequence does not solve every inferential problem. You must still define who could have interacted, decide how the history enters, handle censoring and tied timestamps, and defend any causal language. What it does solve is narrower: the model does not have to invent an unobserved ordering between widely separated panel waves because the ordering is in the data.

::: {.callout-important title="The first question"}
Is the order or timing of the individual interactions part of the substantive quantity you want? If no, a panel or aggregated network may be enough. If yes, preserve the event stream.
:::

# 2. The Ordinal Model in Plain Language

We begin with an **ordinal relational event model**. It conditions on the fact that an event occurred and asks which eligible sender-receiver pair produced it. The waiting time between events is not used.

Immediately before event $m$, the model performs four jobs:

1. Define the risk set $\mathcal{R}_m$, the dyads that could have produced the next event.
2. Calculate history statistics for every eligible dyad, such as recent interaction, receiver popularity, or an immediate return of the preceding event.
3. Give each eligible dyad a score $\eta_{ij,m}=\mathbf{x}_{ij,m}^{\mathsf T}\boldsymbol\beta$.
4. Compare the dyad that actually acted with all the dyads that could have acted.

If event $m$ is the observed ordered pair $(i_m,j_m)$ and $\mathcal{H}_{m-1}$ is the history immediately before it, the model assigns

$$
\Pr\{(i_m,j_m)\mid\mathcal{H}_{m-1}\}
=
\frac{\exp\{\eta_{i_mj_m,m}\}}
{\sum_{(i,j)\in\mathcal{R}_m}\exp\{\eta_{ij,m}\}}.
$$

The fitted coefficients maximize the sum of the log probabilities assigned to the observed sequence:

$$
\ell(\boldsymbol\beta)
=
\sum_m\left[
\eta_{i_mj_m,m}
-
\log\sum_{(i,j)\in\mathcal{R}_m}\exp\{\eta_{ij,m}\}
\right].
$$

The numerator scores the pair that interacted. The denominator scores every pair that could have interacted. The estimator chooses coefficients that give the observed sequence more relative weight.

::: {.callout-note collapse="true" title="Depth: How to read an ordinal coefficient"}
A positive coefficient means that an eligible dyad with a larger value of that history statistic receives more relative weight in the next-event choice, holding the other included statistics fixed.

Exponentiating a one-unit coefficient difference gives a relative choice weight. This resembles the SAOM interpretation from Day 12 because both models compare one realized action with a set of alternatives. It is not an unconditional probability that a tie exists, and it is not a clock-time hazard ratio.

The scale of a statistic matters. A binary indicator for an immediate return and a continuous popularity measure do not have comparable one-unit changes, so ranking effects by raw coefficient magnitude is usually misleading.
:::

# 3. The Phone-Call Stream

The `Social_Evolution` data shipped with `goldfish` contain calls among students living in an MIT residence. We use the calls observed during a six-week portion of the study. One row records one directed call.

```{r}
#| label: event-data
data("Social_Evolution", package = "goldfish")

head(calls, 4)
c(
  events = nrow(calls),
  students = nrow(actors),
  days = round(diff(range(calls$time)) / 86400)
)
```

There are `r nrow(calls)` calls among `r nrow(actors)` students over about `r round(diff(range(calls$time)) / 86400)` days. The first check is not a model check. It is a data check: confirm that the events are ordered, actor identifiers match the actor table, self-events are handled as intended, and the beginning and end of observation are understood.

```{r}
#| label: event-checks
stopifnot(
  all(calls$sender %in% actors$label),
  all(calls$receiver %in% actors$label),
  all(diff(calls$time) >= 0)
)

checkpoint(
  events = nrow(calls),
  actors = nrow(actors),
  tied_times = sum(duplicated(calls$time))
)
```

Tied timestamps deserve a substantive decision. If the measurement system only records time to the nearest minute, an apparent ordering inside that minute may be artificial. Do not let row order silently decide a sequence the data did not observe.

# 4. Prepare the Stream for `relevent`

`relevent::rem.dyad()` expects integer actor identifiers and a time-ordered matrix with time, sender, and receiver columns.

```{r}
#| label: relevent-prep
ids <- setNames(seq_len(nrow(actors)), actors$label)
el <- cbind(
  time = calls$time - min(calls$time) + 1,
  snd = ids[calls$sender],
  rec = ids[calls$receiver]
)
el <- el[base::order(el[, "time"]), , drop = FALSE]

head(el, 4)
```

The risk set in this example is every ordered pair of distinct students. That is an assumption, not a property of the software. In another application, institutional rules, geography, office holding, or entry and exit may make some dyads ineligible at particular times.

::: {.callout-warning title="An incorrect risk set changes the question"}
If an actor could not have produced an event, that actor should not be in the denominator. Treating impossible dyads as available alternatives makes the estimator compare the observed event with events that could not occur.
:::

# 5. Fit an Ordinal Model

We include four features of the history:

- `RRecSnd`: recency of the reverse-direction event $j\rightarrow i$
- `RSndSnd`: recency of the same-direction event $i\rightarrow j$
- `NTDegRec`: the receiver's accumulated degree
- `PSAB-BA`: an immediate reversal, where the event after $a\rightarrow b$ is $b\rightarrow a$

```{r}
#| label: relevent-fit
rem_calls <- fit_or_load(
  "rem_calls",
  {
    invisible(capture.output(
      fit <- relevent::rem.dyad(
        el,
        n = nrow(actors),
        effects = c("RRecSnd", "RSndSnd", "NTDegRec", "PSAB-BA"),
        ordinal = TRUE,
        hessian = TRUE
      )
    ))
    fit
  }
)

summary(rem_calls)
```

The fit is cached because the point of this document is interpretation, not watching an optimizer. Set `options(session.refit = TRUE)` before running the chunk if you want to estimate it from the event stream yourself.

```{r}
#| label: relevent-table
rem_coef <- rem_calls$coef
rem_se <- setNames(
  sqrt(diag(rem_calls$cov)),
  names(rem_coef)
)

data.frame(
  statistic = names(rem_coef),
  estimate = round(unname(rem_coef), 3),
  standard_error = round(unname(rem_se[names(rem_coef)]), 3),
  row.names = NULL
)
```

# 6. Read What Happened, Not Only the Coefficients

Start with one concrete sequence. The code below finds calls that immediately reverse the preceding call and calculates how quickly those reversals occurred.

```{r}
#| label: immediate-returns
ordered_calls <- calls[base::order(calls$time), ]
immediate_return <- which(
  ordered_calls$sender[-1] ==
    ordered_calls$receiver[-nrow(ordered_calls)] &
  ordered_calls$receiver[-1] ==
    ordered_calls$sender[-nrow(ordered_calls)]
) + 1
return_gaps <- ordered_calls$time[immediate_return] -
  ordered_calls$time[immediate_return - 1]

c(
  immediate_returns = length(immediate_return),
  median_gap_seconds = median(return_gaps)
)
```

The strongest result is turn-taking. When the immediately preceding event was $a\rightarrow b$, the $b\rightarrow a$ alternative receives substantially more relative weight than an otherwise comparable eligible dyad. In this stream, `r length(immediate_return)` calls immediately reverse the event before them, and the median gap is `r median(return_gaps)` seconds.

Same-direction recency and receiver popularity also help distinguish the pair that calls next. Generic reverse-direction recency is much weaker once the model gives immediate returns their own statistic.

The applied conclusion is not “these students are reciprocal” in a broad, timeless sense. It is more specific: **return calls are concentrated immediately after the incoming call**. Aggregating the calls into a snapshot would preserve a coarse reciprocity count while discarding this sequence-specific pattern.

::: {.callout-tip title="A reusable interpretation"}
Holding the other history features fixed, a candidate event receives more or less relative weight when it has the named feature. Always name the candidate event, the comparison alternatives, and the time scale represented by the statistic.
:::

# 7. The Same Stream Through `goldfish`

`goldfish` fits relational event models with a different set of history statistics and also supports actor-oriented event models. We first fit an ordinary dyad-oriented model with inertia, reciprocity, and transitivity.

```{r}
#| label: goldfish-fit
gf_rem <- fit_or_load(
  "gf_rem",
  {
    call_network <- goldfish::defineNetwork(nodes = actors, directed = TRUE)
    call_network <- goldfish::linkEvents(
      x = call_network,
      changeEvent = calls,
      nodes = actors
    )
    calls_dep <- goldfish::defineDependentEvents(
      events = calls,
      nodes = actors,
      defaultNetwork = call_network
    )
    invisible(capture.output(
      fit <- goldfish::estimate(
        calls_dep ~ inertia + recip + trans,
        model = "REM"
      )
    ))
    fit
  }
)

data.frame(
  effect = rownames(gf_rem$names),
  estimate = round(gf_rem$parameters, 3),
  standard_error = round(gf_rem$standardErrors, 3),
  row.names = NULL
)
```

The reciprocity estimate here is strongly positive even though `relevent`'s general reverse-recency term was weak. That is not evidence that one package found reciprocity and the other did not. The statistics assign the sequence to different places.

`goldfish` uses a cumulative reciprocity statistic in this specification and does not include the immediate $ab\text{-}ba$ participation-shift indicator. The immediate returns therefore contribute to its general reciprocity term. The `relevent` specification lets the immediate-return term compete for that pattern. The lesson is the same one you saw with ERGM terms: **the included statistics decide which coefficient receives credit for an observed feature**.

::: {.callout-note collapse="true" title="Depth: Dyad-oriented and actor-oriented event models"}
A dyad-oriented relational event model asks which eligible pair produces the next event.

The dynamic network actor model in `goldfish`, fitted with `model = "DyNAM"`, splits that decision into two stages: which actor acts, then which receiver that actor chooses. This can be useful when the theory itself has those two stages, such as a legislator deciding first to initiate a cosponsorship request and then choosing a colleague.

Neither formulation is universally better. The dyad-oriented model matches a theory about pair-specific event opportunities. The actor-oriented model matches a theory about actors initiating actions and selecting receivers. Choose the one whose decision process matches the question.
:::

# 8. Order-Only and Clock-Time Questions

The two examples above use an ordinal likelihood. They learn from event order and ignore the length of the waiting times. This is often sensible when timestamps are rounded, observation pauses are poorly recorded, or the question is specifically about who acts next.

A clock-time relational event model also describes how quickly the next event arrives. Reach for it when waiting is part of the outcome, such as the speed of diplomatic retaliation, the delay before a protest response, or the time until a legislator answers a request.

Do not call an ordinal coefficient a clock-time hazard ratio. An ordinal model may correctly predict which pair acts next while saying nothing about whether the system waits thirty seconds or thirty days.

# 9. What to Check Before Believing the Result

A converged optimizer only says that the fitting routine found a stable solution to the stated objective. It does not establish that the risk set, history statistics, or substantive interpretation are right.

Before reporting a relational event model, check:

1. **Eligibility:** Who could send to whom at each event?
2. **Observation:** Are starts, stops, missing periods, and actor entry or exit recorded?
3. **Timestamp resolution:** Are tied events genuinely simultaneous, or merely rounded?
4. **History scale:** Does “recent” mean the preceding event, elapsed time, or a decaying count?
5. **Competing statistics:** Would immediate return, cumulative reciprocity, and general recency divide the same pattern differently?
6. **Actor heterogeneity:** Are a few frequent senders or receivers driving what looks like dyadic history?
7. **Prediction and simulation:** Does the model reproduce substantively important sequences, and does it improve held-out event prediction over a simple frequency baseline?
8. **Causal language:** Could common shocks, strategic anticipation, or the observation process explain the sequence?

The model describes conditional event patterns under its assumptions. A positive history coefficient does not by itself show that the preceding event caused the next one.

# 10. Exercises

::: {.panel-tabset}

## Core

Explain why `relevent`'s immediate-return coefficient and `goldfish`'s cumulative reciprocity coefficient can both be positive without estimating the same quantity. Your answer should name the event sequence each statistic credits.

## Application

Choose one timestamped process you might study. In four sentences, define:

1. One event
2. Its sender and receiver
3. The eligible risk set
4. One history statistic that answers a meaningful question

Then say whether event order is enough or whether waiting time is part of the outcome.

## Solution

The immediate-return statistic credits the specific sequence $a\rightarrow b$ followed directly by $b\rightarrow a$. A cumulative reciprocity statistic can credit a reverse-direction event even when other events occur between the two. Both can be positive because immediate returns are one part of the broader reciprocal history. When both are included, they compete to explain that shared part of the sequence.

A strong application answer makes ineligible events explicit. For example, in a stream of legislative speeches and named replies, one event is legislator $i$ replying to legislator $j$. The risk set contains legislators present and permitted to speak at that moment. A recent-reply statistic can ask whether a direct reply is likely to be answered in turn. Clock time matters if the argument concerns speed of response; order alone is enough if the argument concerns who responds next.

:::

# 11. Reading and Versions

- Butts (2008), "A Relational Event Framework for Social Action," *Sociological Methodology*
- Stadtfeld and Block (2017), "Interactions, Actors, and Time: Dynamic Network Actor Models for Relational Events," *Sociological Science*
- Butts, Lomi, Snijders, and Stadtfeld (2023), "Relational Event Models in Network Science," *Network Science*
- Bianchi, Filippi-Mazzola, Lomi, and Wit (2024), "Relational Event Modeling," *Annual Review of Statistics and Its Application*

```{r}
#| label: versions
#| cache: false
data.frame(
  package = c("R", "relevent", "goldfish"),
  version = c(
    as.character(getRversion()),
    as.character(packageVersion("relevent")),
    as.character(packageVersion("goldfish"))
  )
)
```
