Relational Event Models

Ancillary self-study · ICPSR Network Analysis: Advanced Topics

Author

Shahryar Minhas

Published

July 29, 2026

ImportantOptional: not covered in the 2026 course

This document is not part of any taught session. Nobody is going to lecture it at you, nobody is going to grade you on it, and you can ignore it completely without falling behind.

It exists because it is a thing you will plausibly need next, and because it is easier to hand you something that works than to hand you a citation and wish you luck. Everything here runs end to end on the same machine setup you used in class. Every number and every figure on this page was produced by a chunk in this file: there are no screenshots of R output anywhere in the ancillary package.

How to work through it alone. Get the source (click the </> Code button at the top right, then View Source, and copy it into a .qmd file: or open the .qmd you were given directly), run the chunks in order, and read the prose between them: the prose is where the argument lives. Slow fits are identified in the prose, and precomputed fits load a shipped result while still showing you the code that produced it. Folded Depth boxes are optional extras; folded Stuck? boxes are hints. Open whichever you need and ignore the other. Exercises have a Solution tab: you are working alone, so use it.

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

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

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.

data("Social_Evolution", package = "goldfish")

head(calls, 4)
#>         time   sender receiver increment
#> 1 1220733470 Actor 72 Actor 50         1
#> 2 1221102974 Actor 43 Actor 51         1
#> 3 1221784293 Actor 43 Actor 51         1
#> 4 1221785882 Actor 43 Actor 22         1
c(
  events = nrow(calls),
  students = nrow(actors),
  days = round(diff(range(calls$time)) / 86400)
)
#>   events students     days 
#>      439       84       42

There are 439 calls among 84 students over about 42 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.

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))
)
#> ------------------------------------------------------------------
#> CHECKPOINT: events = 439   |   actors = 84   |   tied_times = 0
#> ------------------------------------------------------------------

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.

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)
#>             time snd rec
#> Actor 72       1  72  50
#> Actor 43  369505  43  51
#> Actor 43 1050824  43  51
#> Actor 43 1052413  43  22

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.

WarningAn 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\)
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)
#> Relational Event Model (Ordinal Likelihood)
#> 
#>          Estimate Std.Err Z value Pr(>|z|)    
#> NTDegRec  7.25357 0.64195 11.2993   <2e-16 ***
#> RRecSnd   0.12253 0.16318  0.7509   0.4527    
#> RSndSnd   6.24911 0.14072 44.4091   <2e-16 ***
#> PSAB-BA   3.58481 0.17631 20.3323   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Null deviance: 7769.999 on 439 degrees of freedom
#> Residual deviance: 3525.881 on 435 degrees of freedom
#>  Chi-square: 4244.118 on 4 degrees of freedom, asymptotic p-value 0 
#> AIC: 3533.881 AICC: 3533.973 BIC: 3550.219

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.

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
)
#>   statistic estimate standard_error
#> 1  NTDegRec    7.254          0.642
#> 2   RRecSnd    0.123          0.163
#> 3   RSndSnd    6.249          0.141
#> 4   PSAB-BA    3.585          0.176

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.

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)
)
#>  immediate_returns median_gap_seconds 
#>                153                 20

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, 153 calls immediately reverse the event before them, and the median gap is 20 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.

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

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
)
#>    effect estimate standard_error
#> 1 inertia    6.198          0.151
#> 2   recip    1.434          0.115
#> 3   trans    0.223          0.157

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.

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

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.

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.

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
data.frame(
  package = c("R", "relevent", "goldfish"),
  version = c(
    as.character(getRversion()),
    as.character(packageVersion("relevent")),
    as.character(packageVersion("goldfish"))
  )
)
#>    package version
#> 1        R   4.3.3
#> 2 relevent   1.2.1
#> 3 goldfish  1.6.12