# The Relational Event Modeling Pipeline

### A comprehensive tutorial

####

The `remverse` suite provides a three-step pipeline for relational event modeling:

1. `remify()` — process a raw edgelist into a standardized event history
2. `remstats()` — compute network statistics (endogenous and exogenous)
3. `remstimate()` — estimate model parameters

This vignette walks through the full pipeline for a range of modeling choices: tie-oriented vs. actor-oriented models, interval vs. ordinal timing, different risk set definitions, typed events, and a variety of endogenous and exogenous statistics. Each section builds on the same dataset so the results can be compared directly.

    library(remverse)

# 1 Data

Throughout the `remverse` vignettes we use the simulated dataset `randomREH3`, provided with `remverse`. It is a sequence of 999 directed events among 5 actors, generated from a known relational event model so that the recovered coefficients can be checked against the truth.

    data(randomREH3)   # event history
    data(info3)        # actor attributes

    head(randomREH3)
    #>   time actor1 actor2 end setting duration
    #> 2   18      5      1  24  social        6
    #> 3   57      3      2  70  social       13
    #> 4   64      1      5  85  social       21
    #> 5  127      5      1 159  social       32
    #> 6  228      1      4 239    work       11
    #> 7  274      2      4 320    work       46
    head(info3)
    #>   name  time age sex
    #> 1    1     0  24   0
    #> 2    2     0  28   0
    #> 3    3     0  19   1
    #> 4    4     0  20   1
    #> 5    4 10000  21   1
    #> 6    5     0  28   0

Each row of `randomREH3` records a `time` (event start time), a `sender` (`actor1`), a `receiver` (`actor2`), a `setting` (event type: `"work"` or `"social"`), an `end` time, and a `duration` (`end - time`). The `duration` column can double as an intrinsic event `weight`. The `end`/`duration` columns are only used by the duration model (see the Duration Relational Event Models vignette); the tie- and actor-oriented models in this vignette use `time`, `actor1`, `actor2`, and `setting` only.

`info3` provides actor attributes: `age` and `sex`, indexed by actor `name` and a `time` of change (so age can be time-varying — actors 4 and 5 age during the observation window).

# 2 Tie-oriented models

## 2.1 Basic model (interval timing)

The simplest specification: a directed tie-oriented model with interval timing (exact event times), a full risk set, and a few endogenous statistics.

    # Step 1: process the event history
    reh <- remify(edgelist = randomREH3, model = "tie", directed = TRUE)
    summary(reh)
    #> Relational Event History processed for tie-oriented modeling:
    #>  > events = 999 (time points = 993)
    #>  > actors = 5
    #>  > event types = 1
    #>  > riskset = full
    #>      >> included dyads = 20
    #>  > directed = TRUE
    #>  > ordinal = FALSE
    #>  > weighted = FALSE
    #>  > time length ~ 50676
    #>  > interevent time
    #>       >> minimum ~ 1
    #>       >> median  ~ 36
    #>       >> maximum ~ 448

    # Step 2: define effects and compute statistics
    effects <- ~ inertia(scaling = "std") + reciprocity(scaling = "std") +
                 outdegreeSender(scaling = "std") + isp(scaling = "std")
    stats <- remstats(reh = reh, tie_effects = effects)
    stats
    #> Relational Event Network Statistics
    #> > Model: tie-oriented
    #> > Computation method: per time point
    #> > Dimensions: 992 time points x 20 dyads x 5 statistics
    #> > Statistics:
    #>   >> 1: baseline
    #>   >> 2: inertia
    #>   >> 3: reciprocity
    #>   >> 4: outdegreeSender
    #>   >> 5: isp

    # Step 3: estimate the model
    fit <- remstimate(reh = reh, stats = stats)
    summary(fit)
    #> Relational Event Model (tie oriented)
    #>
    #> Call:
    #> ~inertia(scaling = "std") + reciprocity(scaling = "std") + outdegreeSender(scaling = "std") +     isp(scaling = "std")
    #>
    #>
    #> Coefficients (MLE with interval likelihood):
    #>
    #>                    Estimate    Std. Err     z value Pr(>|z|)    Pr(=0)
    #> baseline          -7.253328    0.043175 -167.999084    0e+00 < 2.2e-16
    #> inertia            0.421529    0.048117    8.760442    0e+00 6.811e-16
    #> reciprocity        0.201510    0.047128    4.275845    0e+00  0.003363
    #> outdegreeSender   -0.121893    0.034797   -3.502931    5e-04  0.063840
    #> isp               -0.363749    0.060352   -6.027102    0e+00 4.075e-07
    #> Null deviance: 15813.97 on 992 degrees of freedom
    #> Residual deviance: 15227.19 on 987 degrees of freedom
    #> Chi-square: 586.7778 on 5 degrees of freedom, asymptotic p-value 0
    #> AIC: 15237.19 AICC: 15237.25 BIC: 15261.69

The baseline intercept is automatically included for interval timing (`ordinal = FALSE`). It captures the average event rate across all dyads. The endogenous effects quantify how the accumulated event history shapes future interaction patterns: `inertia` measures the tendency to repeat past interactions, `reciprocity` captures the tendency to reciprocate received events, `outdegreeSender` captures sender activity, and `isp` (incoming shared partners) captures a form of triadic closure. Standardising the statistics (`scaling = "std"`) puts the coefficients on a comparable scale.

By default `remstimate()` fits the model by maximum likelihood (the `frequentist` approach). The estimation approach and model structure are controlled by the `approach`, `random`, `penalty`, and `mixture` arguments — see the dedicated vignettes on frailty (GLMM), penalized, and mixture models.

### Diagnostics

    diag <- diagnostics(object = fit, reh = reh, stats = stats)
    diag
    #> Diagnostics for a Relational Event Model
    #> ------------------------------------------
    #> Model       : tie
    #> Actors      : 5
    #> Events      : 992
    #>
    #> Tie model
    #>   Statistics : inertia, reciprocity, outdegreeSender, isp
    #>   Recall     : mean rank = 0.724  |  prob ratio = 1.61  |  top 5% = 9.8%  |  lowest 20% = 2.3%

    plot(x = fit, reh = reh, diagnostics = diag)

    #> Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
    #> collapsing to unique 'x' values

The `diagnostics()` function computes recall statistics: for each time point, the observed event is ranked among all competing dyads based on its predicted rate. A well-fitting model assigns high ranks to observed events.

## 2.2 Ordinal timing

When only the order of events matters (or when exact timestamps are unreliable), set `ordinal = TRUE`. The model reduces to a stratified Cox proportional-hazard model — no baseline intercept is estimated.

    reh_ord <- remify(edgelist = randomREH3, model = "tie", ordinal = TRUE)

    stats_ord <- remstats(reh = reh_ord,
                          tie_effects = ~ inertia(scaling = "std") +
                                          reciprocity(scaling = "std"))
    fit_ord <- remstimate(reh = reh_ord, stats = stats_ord)
    summary(fit_ord)
    #> Relational Event Model (tie oriented)
    #>
    #> Call:
    #> ~inertia(scaling = "std") + reciprocity(scaling = "std")
    #>
    #>
    #> Coefficients (MLE with ordinal likelihood):
    #>
    #>              Estimate  Std. Err   z value Pr(>|z|)    Pr(=0)
    #> inertia      0.435997  0.041208 10.580335        0 < 2.2e-16
    #> reciprocity  0.298866  0.042690  7.000838        0  7.17e-10
    #> Null deviance: 5936.275 on 992 degrees of freedom
    #> Residual deviance: 5406.784 on 990 degrees of freedom
    #> Chi-square: 529.4906 on 2 degrees of freedom, asymptotic p-value 0
    #> AIC: 5410.784 AICC: 5410.796 BIC: 5420.584

The coefficients are comparable to the interval model but the intercept and inter-event-time offset are absent.

## 2.3 Risk set variations

The risk set is the collection of dyads considered “at risk” of an event at each time point — the denominator against which each observed event is compared. It is chosen with the `riskset` argument of `remify()`, which takes four values.

The default, `riskset = "full"`, puts every possible directed dyad at risk at every time point: \(D = N(N-1)\) dyads (or \(N(N-1)/2\) undirected), fixed over the whole sequence. This is the natural choice for small, densely connected networks like the one used here. In larger or sparser networks, however, a full risk set is often not realistic: it assumes that any actor could plausibly interact with any other at any moment, and it grows quadratically in the number of actors, which quickly becomes both substantively implausible and computationally expensive. The alternatives below restrict the risk set to something more defensible.

### Active risk set

`riskset = "active"` includes only dyads observed at least once anywhere in the sequence — a common, computationally lighter choice for sparse networks.

    reh_active <- remify(edgelist = randomREH3, model = "tie", riskset = "active")

    cat("Full risk set:  ", reh$D, "dyads\n")
    #> Full risk set:   20 dyads
    cat("Active risk set:", reh_active$activeD, "dyads\n")
    #> Active risk set: 20 dyads

    stats_active <- remstats(reh = reh_active,
                             tie_effects = ~ inertia(scaling = "std") +
                                             reciprocity(scaling = "std"))
    fit_active <- remstimate(reh = reh_active, stats = stats_active)
    summary(fit_active)
    #> Relational Event Model (tie oriented)
    #>
    #> Call:
    #> ~inertia(scaling = "std") + reciprocity(scaling = "std")
    #>
    #>
    #> Coefficients (MLE with interval likelihood):
    #>
    #>                Estimate    Std. Err     z value Pr(>|z|)    Pr(=0)
    #> baseline      -7.182955    0.039269 -182.917348        0 < 2.2e-16
    #> inertia        0.443331    0.044428    9.978525        0 < 2.2e-16
    #> reciprocity    0.291811    0.045432    6.423043        0 3.465e-08
    #> Null deviance: 15813.97 on 992 degrees of freedom
    #> Residual deviance: 15278.25 on 989 degrees of freedom
    #> Chi-square: 535.7245 on 3 degrees of freedom, asymptotic p-value 0
    #> AIC: 15284.25 AICC: 15284.27 BIC: 15298.95

### Active-saturated risk set

`riskset = "active_saturated"` extends the active risk set by adding, for every observed actor pair, the reverse direction (if \(A \to B\) is observed, \(B \to A\) is also placed at risk) and — when events are typed — all event types for that pair. This encodes the assumption that observing any interaction between two actors implies that both directions (and all types) are possible for them, without opening up the full risk set to actor pairs that never interact at all.

    reh_sat <- remify(edgelist = randomREH3, model = "tie",
                      riskset = "active_saturated")
    cat("Active-saturated risk set:", reh_sat$activeD, "dyads\n")
    #> Active-saturated risk set: 20 dyads

### Manual risk set

`riskset = "manual"` restricts the risk set to a user-supplied set of dyads, passed via `manual_riskset`. Observed dyads not included in the manual specification are added automatically.

    # Define a manual risk set: a subset of directed dyads among actors 1-5
    my_riskset <- data.frame(
      actor1 = c(1, 2, 3, 4, 5),
      actor2 = c(2, 3, 4, 5, 1)
    )

    reh_manual <- remify(
      edgelist       = randomREH3,
      model          = "tie",
      riskset        = "manual",
      manual_riskset = my_riskset
    )
    #> Warning in remify(edgelist = randomREH3, model = "tie", riskset = "manual", :
    #> 15 observed dyads were added to the manual risk set.
    cat("Manual risk set:", reh_manual$activeD, "dyads\n")
    #> Manual risk set: 20 dyads

### Extending the risk set by event type

Orthogonal to the four risk-set types above is the `extend_riskset_by_type` argument, relevant only when events carry a type (via `event_type`). With `extend_riskset_by_type = TRUE`, each actor pair is duplicated for every event type, so a dyad-type combination — not just a dyad — is the unit at risk, and the risk set has size \(D = N(N-1) \times C\) (directed, \(C\) types). With `FALSE` (the default), the event type is treated as a mark on events only and does not expand the risk set. This choice is revisited in the typed-events section below.

## 2.4 Typed events

When events carry a type label (here, `"work"` vs. `"social"` in the `setting` column), the `consider_type` argument controls how statistics account for event types. The risk set can optionally be extended by type (`extend_riskset_by_type = TRUE`), so that each dyad-type combination is a separate unit at risk.

    reh_typed <- remify(
      edgelist   = randomREH3,
      model      = "tie",
      event_type = "setting"
    )

### `consider_type = "ignore"`

Event types are not distinguished — statistics aggregate over all types.

    stats_ign <- remstats(
      reh = reh_typed,
      tie_effects = ~ inertia(consider_type = "ignore")
    )
    dimnames(stats_ign)[[3]]
    #> [1] "baseline" "inertia"

### `consider_type = "separate"`

One statistic per type: separate counts of past `"work"` events and past `"social"` events.

    stats_sep <- remstats(
      reh = reh_typed,
      tie_effects = ~ inertia(consider_type = "separate")
    )
    dimnames(stats_sep)[[3]]
    #> [1] "baseline"       "inertia.social" "inertia.work"

### `consider_type = "interact"`

The statistic conditions on the type of the event at the current time point, allowing the effect of event history to depend on the triggering type.

    stats_int <- remstats(
      reh = reh_typed,
      tie_effects = ~ inertia(consider_type = "interact")
    )
    #> Warning in prepare_tomstats(effects = effects, reh = reh, attr_actors =
    #> attr_actors, : "interact" requires extend_riskset_by_type = TRUE; coercing to
    #> "separate".
    dimnames(stats_int)[[3]]
    #> [1] "baseline"       "inertia.social" "inertia.work"

### Estimation with typed events

    stats_typed <- remstats(
      reh = reh_typed,
      tie_effects = ~ inertia(consider_type = "separate") +
                       reciprocity(consider_type = "separate")
    )

    fit_typed <- remstimate(reh = reh_typed, stats = stats_typed)
    summary(fit_typed)
    #> Relational Event Model (tie oriented)
    #>
    #> Call:
    #> ~inertia(consider_type = "separate") + reciprocity(consider_type = "separate")
    #>
    #>
    #> Coefficients (MLE with interval likelihood):
    #>
    #>                       Estimate    Std. Err     z value Pr(>|z|)    Pr(=0)
    #> baseline           -7.2518e+00  5.0306e-02 -1.4415e+02   0.0000 < 2.2e-16
    #> inertia.social     -4.9709e-03  4.2046e-03 -1.1823e+00   0.2371 0.9399697
    #> inertia.work        5.4111e-02  1.3646e-02  3.9653e+00   0.0001 0.0119841
    #> reciprocity.social -1.7126e-02  4.2649e-03 -4.0155e+00   0.0001 0.0098328
    #> reciprocity.work    6.3544e-02  1.3626e-02  4.6633e+00   0.0000 0.0005967
    #> Null deviance: 15813.97 on 992 degrees of freedom
    #> Residual deviance: 15527.04 on 987 degrees of freedom
    #> Chi-square: 286.9303 on 5 degrees of freedom, asymptotic p-value 0
    #> AIC: 15537.04 AICC: 15537.1 BIC: 15561.54

The type-separated coefficients tell us whether inertia and reciprocity operate differently across event types — for instance, whether `"work"` interactions are more habitual than `"social"` ones.

### Extending the risk set by type

The `consider_type` argument (above) changes how statistics are computed, but by default the risk set still contains one entry per dyad. Setting `extend_riskset_by_type = TRUE` in `remify()` goes a step further and makes each dyad-type combination its own unit at risk, so the model treats “a work event from \(i\) to \(j\)” and “a social event from \(i\) to \(j\)” as competing, distinct events. The risk set then has size \(D = N(N-1) \times C\).

    reh_typed_ext <- remify(
      edgelist               = randomREH3,
      model                  = "tie",
      event_type             = "setting",
      extend_riskset_by_type = TRUE
    )
    cat("Dyad-only risk set:      ", reh_typed$D, "\n")
    #> Dyad-only risk set:       20
    cat("Type-extended risk set:  ", reh_typed_ext$D, "\n")
    #> Type-extended risk set:   40

Use `extend_riskset_by_type = TRUE` when the choice of type is itself part of the event process you want to model; keep the default `FALSE` when the type is merely an attribute recorded alongside each event.

## 2.5 Exogenous statistics

Exogenous actor attributes and dyadic covariates can be included alongside endogenous effects. The `info3` dataset provides the actor covariate `age`.

### Actor-level effects

    stats_exo <- remstats(
      reh = reh,
      tie_effects = ~ inertia(scaling = "std") +
                       reciprocity(scaling = "std") +
                       send("age", attr_actors = info3, scaling = "std") +
                       receive("age", attr_actors = info3, scaling = "std") +
                       difference("age", attr_actors = info3, scaling = "std")
    )

    fit_exo <- remstimate(reh = reh, stats = stats_exo)
    summary(fit_exo)
    #> Relational Event Model (tie oriented)
    #>
    #> Call:
    #> ~inertia(scaling = "std") + reciprocity(scaling = "std") + send("age",     attr_actors = info3, scaling = "std") + receive("age", attr_actors = info3,     scaling = "std") + difference("age", attr_actors = info3,     scaling = "std")
    #>
    #>
    #> Coefficients (MLE with interval likelihood):
    #>
    #>                   Estimate    Std. Err     z value Pr(>|z|)    Pr(=0)
    #> baseline         -7.238625    0.041354 -175.041072   0.0000 < 2.2e-16
    #> inertia           0.402356    0.054313    7.408075   0.0000 3.813e-11
    #> reciprocity       0.240859    0.055845    4.312992   0.0000  0.002869
    #> send_age         -0.018312    0.036131   -0.506837   0.6123  0.965156
    #> receive_age      -0.061821    0.035888   -1.722593   0.0850  0.877202
    #> difference_age   -0.322162    0.038607   -8.344572   0.0000 2.387e-14
    #> Null deviance: 15813.97 on 992 degrees of freedom
    #> Residual deviance: 15200.25 on 986 degrees of freedom
    #> Chi-square: 613.7194 on 6 degrees of freedom, asymptotic p-value 0
    #> AIC: 15212.25 AICC: 15212.34 BIC: 15241.65

- `send("age")`: the sender’s age — do older actors initiate more events?
- `receive("age")`: the receiver’s age — are older actors more often chosen as receivers?
- `difference("age")`: the absolute age difference within the dyad — do events flow between similarly aged actors?

### Interactions between endogenous and exogenous effects

Effects can be interacted to test whether endogenous tendencies (e.g., inertia) depend on actor attributes:

    stats_ix <- remstats(
      reh = reh,
      tie_effects = ~ inertia(scaling = "std") +
                       send("age", attr_actors = info3, scaling = "std") +
                       inertia(scaling = "std"):send("age", attr_actors = info3, scaling = "std")
    )

    fit_ix <- remstimate(reh = reh, stats = stats_ix)
    summary(fit_ix)
    #> Relational Event Model (tie oriented)
    #>
    #> Call:
    #> ~inertia(scaling = "std") + send("age", attr_actors = info3,     scaling = "std") + inertia(scaling = "std"):send("age", attr_actors = info3,     scaling = "std")
    #>
    #>
    #> Coefficients (MLE with interval likelihood):
    #>
    #>                     Estimate    Std. Err     z value Pr(>|z|) Pr(=0)
    #> baseline           -7.153090    0.038310 -186.717881   0.0000 <2e-16
    #> inertia             0.674666    0.029403   22.945394   0.0000 <2e-16
    #> send_age           -0.064890    0.039264   -1.652652   0.0984 0.8894
    #> inertia:send_age   -0.069468    0.030007   -2.315038   0.0206 0.6836
    #> Null deviance: 15813.97 on 992 degrees of freedom
    #> Residual deviance: 15301.09 on 988 degrees of freedom
    #> Chi-square: 512.8841 on 4 degrees of freedom, asymptotic p-value 0
    #> AIC: 15309.09 AICC: 15309.13 BIC: 15328.69

A significant interaction would indicate that the tendency to repeat past interactions depends on the sender’s age.

## 2.6 Memory types

The `memory` argument controls which past events contribute to endogenous statistics at each time point.

    # Full memory (default): entire event history
    stats_full <- remstats(reh = reh, tie_effects = ~ inertia(), memory = "full")

    # Window memory: only events within the last 2000 time units
    stats_win <- remstats(reh = reh, tie_effects = ~ inertia(),
                          memory = "window", memory_value = 2000)

    # Decay memory: exponential decay with half-life 1000
    stats_dec <- remstats(reh = reh, tie_effects = ~ inertia(),
                          memory = "decay", memory_value = 1000)

    # Compare the inertia statistic at the last time point for the first dyad
    M <- dim(stats_full)[1]
    cat("Full memory:   ", stats_full[M, 1, "inertia"], "\n")
    #> Full memory:    4
    cat("Window (2000): ", stats_win[M, 1, "inertia"], "\n")
    #> Window (2000):  0
    cat("Decay (1000):  ", stats_dec[M, 1, "inertia"], "\n")
    #> Decay (1000):   3.348464e-07

Full memory accumulates all past events equally; window memory discards events older than the threshold; decay memory down-weights older events exponentially.

Full memory is the default, but it is often not the most realistic choice: it implies that an interaction from the distant past influences the current rate exactly as strongly as one that happened moments ago. In most social processes, influence fades. The decay memory captures this with a single half-life parameter (`memory_value`): the weight of a past event halves every `memory_value` time units. Rather than fixing the half-life arbitrarily, it can be tuned like any other modeling choice — by fitting the model across a grid of half-lives and comparing information criteria (AIC/BIC) or predictive recall, and selecting the value that optimizes them (see the model-assessment section below). For the `randomREH3` data a half-life around 2000 is close to optimal, and we adopt `memory = "decay"`, `memory_value = 2000` as the default in the more advanced vignettes.

It is instructive to see where the three numbers above come from. The first dyad (the directed pair 1 → 2) has four past events in `randomREH3`, at times 3951, 5190, 9987 and 29149, and the last time point at which the statistic is evaluated is around 50676. Under full memory each of the four events counts once, giving 4. Under window memory with a width of 2000, only events in the trailing interval — here roughly the last 2000 time units before the final time point — are counted; since the most recent 1 → 2 event (at 29149) is far outside that window, the count is 0. Under decay memory with a half-life of 1000, each past event is down-weighted by one half for every 1000 time units that have elapsed since it occurred (measured, in the default `"pt"` scheme, from the previous time point). The three oldest events are effectively forgotten, and the most recent one — about 21500 time units, or roughly 21.5 half-lives, in the past — contributes a weight of about `0.5^21.5`, so the whole statistic collapses to roughly 3.3e-07. The window value of 0 and the tiny decay value therefore tell the same story from two angles: this dyad has simply not been active recently.

## 2.7 Case-control sampling

For large networks, computing statistics over the full risk set at every time point is expensive. Case-control sampling draws a random subset of non-event dyads as the comparison set.

    stats_sampled <- remstats(
      reh = reh,
      tie_effects = ~ inertia(scaling = "std") + reciprocity(scaling = "std"),
      sampling = TRUE,
      samp_num = 5,
      seed = 42
    )

    fit_sampled <- remstimate(reh = reh, stats = stats_sampled)
    summary(fit_sampled)
    #> Relational Event Model (tie oriented)
    #>
    #> Call:
    #> ~inertia(scaling = "std") + reciprocity(scaling = "std")
    #>
    #>
    #> Coefficients (MLE with interval likelihood):
    #>
    #>                Estimate    Std. Err     z value Pr(>|z|)    Pr(=0)
    #> baseline      -7.145876    0.038694 -184.675932        0 < 2.2e-16
    #> inertia        0.465151    0.056679    8.206732        0 7.470e-14
    #> reciprocity    0.387469    0.056400    6.870020        0 1.776e-09
    #> Null deviance: 15813.97 on 992 degrees of freedom
    #> Residual deviance: 15146.15 on 989 degrees of freedom
    #> Chi-square: 667.8199 on 3 degrees of freedom, asymptotic p-value 0
    #> AIC: 15152.15 AICC: 15152.18 BIC: 15166.85

The estimates are consistent as the sample size grows but will differ slightly from the full-risk-set estimates due to sampling variability.

## 2.8 Model assessment and variable selection

Having seen the individual modeling choices, how do we decide which effects to keep? Different criteria are available directly from a fitted model: the AIC and BIC (in `fit$AIC` / `fit$BIC`, trading fit against the number of parameters, with BIC penalizing complexity more heavily), and the median recall from `diagnostics()` (a prediction-based measure — higher is better). We use a `decay` memory with half-life 2000 from here on. Additionally, p-values, `Pr(>|z|)`, or posterior probability of a null value, `Pr(=0)`, (based on the default Bayes factor methodology of R packages `BFpack` and `bain`) can be used to inform user about including or excluding certain effects.

The helper below fits a tie model and returns these quantities, so several candidate specifications can be compared on a common footing. We use `first = 200` so every model is scored on the same set of events.

    assess <- function(formula, memory_value = 2000) {
      s <- remstats(reh, tie_effects = formula, first = 200,
                    memory = "decay", memory_value = memory_value)
      m <- remstimate(reh, s)
      d <- diagnostics(m, reh, s)
      data.frame(
        npar          = length(coef(m)),
        AIC           = round(m$AIC, 1),
        BIC           = round(m$BIC, 1),
        median_recall = round(d$recall$summary$median_rel_rank, 3)
      )
    }

We compare four specifications: an empty (baseline-only) model, an endogenous-only model (`no_age`), that model plus the sender’s `age` (`wage`), and a richer model with extra degree and shared-partner effects.

    f_empty <- ~ 1
    f_noage <- ~ inertia(scaling = "std") + reciprocity(scaling = "std") +
                 outdegreeSender(scaling = "std") + isp(scaling = "std")
    f_wage  <- ~ inertia(scaling = "std") + reciprocity(scaling = "std") +
                 outdegreeSender(scaling = "std") +
                 send("age", attr_actors = info3, scaling = "std") +
                 isp(scaling = "std")
    f_rich  <- ~ inertia(scaling = "std") + reciprocity(scaling = "std") +
                 outdegreeSender(scaling = "std") + indegreeReceiver(scaling = "std") +
                 send("age", attr_actors = info3, scaling = "std") +
                 receive("age", attr_actors = info3, scaling = "std") +
                 isp(scaling = "std") + osp(scaling = "std")

    comparison <- rbind(
      empty  = assess(f_empty),
      no_age = assess(f_noage),
      wage   = assess(f_wage),
      rich   = assess(f_rich)
    )
    comparison
    #>        npar     AIC     BIC median_recall
    #> empty     1 12776.3 12781.0         0.500
    #> no_age    5 12086.1 12109.5         0.842
    #> wage      6 12088.0 12116.1         0.842
    #> rich      9 12086.2 12128.3         0.842

The empty model anchors the scale: its median recall sits at 0.5 (no better than chance) and its information criteria are worse. Adding the endogenous effects improves everything sharply. Among the three non-empty models the median recall is identical (0.842) — with only five actors the risk set has just 20 dyads, so the median relative rank takes few distinct values and is too coarse to separate them here; the information criteria do the discriminating. Both AIC and BIC are lowest for the parsimonious `no_age` model: adding the sender’s `age` (`wage`) or the extra structural effects (`rich`) leaves the fit essentially unchanged while spending more parameters, so both criteria go slightly up. The verdict is that `age` is not supported by these data, and the endogenous-only model is preferred.

The per-effect output of `summary()` gives a complementary, one-model-at-a-time view that corroborates this. For each coefficient it reports a frequentist p-value in the `Pr(>|z|)` column and a posterior probability that the effect equals zero in the `Pr(=0)` column. In the model that includes `age`, the `send.age` row carries a large p-value and a high posterior probability of a zero effect — the same signal the model comparison gave, namely that this effect can be dropped.

    stats_wage <- remstats(reh, tie_effects = f_wage, first = 200,
                           memory = "decay", memory_value = 2000)
    summary(remstimate(reh, stats_wage))
    #> Relational Event Model (tie oriented)
    #>
    #> Call:
    #> ~inertia(scaling = "std") + reciprocity(scaling = "std") + outdegreeSender(scaling = "std") +     send("age", attr_actors = info3, scaling = "std") + isp(scaling = "std")
    #>
    #>
    #> Coefficients (MLE with interval likelihood):
    #>
    #>                    Estimate    Std. Err     z value Pr(>|z|)    Pr(=0)
    #> baseline          -7.523964    0.056829 -132.397108   0.0000 < 2.2e-16
    #> inertia            0.478865    0.045637   10.492805   0.0000 < 2.2e-16
    #> reciprocity        0.220127    0.042121    5.226112   0.0000 3.305e-05
    #> outdegreeSender   -0.126343    0.042465   -2.975208   0.0029    0.2521
    #> send_age           0.010830    0.037216    0.291015   0.7710    0.9643
    #> isp               -0.612400    0.066844   -9.161629   0.0000 < 2.2e-16
    #> Null deviance: 12774.29 on 794 degrees of freedom
    #> Residual deviance: 12076.04 on 788 degrees of freedom
    #> Chi-square: 698.2479 on 6 degrees of freedom, asymptotic p-value 0
    #> AIC: 12088.04 AICC: 12088.15 BIC: 12116.1

### Tuning the memory half-life

The recall and information criteria can also be used to tune the memory value (half-life for a `decay` memory). Refitting the preferred (`no_age`) specification across a grid of half-lives and tracking BIC and recall points to a best value — for these data, BIC bottoms out near a half-life of 2000 and the relative recall is maximal, which is why we use it

    half_lives <- c(500, 1000, 2000, 3000, 4000)
    mem_path <- do.call(rbind, lapply(half_lives, function(h)
      cbind(half_life = h, assess(f_noage, memory_value = h))))
    mem_path[, c("half_life", "BIC", "median_recall")]
    #>   half_life     BIC median_recall
    #> 1       500 12198.0         0.842
    #> 2      1000 12125.5         0.842
    #> 3      2000 12109.5         0.842
    #> 4      3000 12123.7         0.789
    #> 5      4000 12144.3         0.789

    plot(mem_path$half_life, mem_path$BIC, type = "b", pch = 19,
         xlab = "decay half-life", ylab = "BIC", main = "BIC vs. memory half-life")

# 3 Actor-oriented models

Actor-oriented models decompose the event process into two steps: a sender rate model (which actor initiates next?) and a receiver choice model (whom does the sender choose?). This requires `directed = TRUE`.

## 3.1 Sender + receiver model

    reh_ao <- remify(edgelist = randomREH3, model = "actor", directed = TRUE)

    rate_effects   <- ~ outdegreeSender(scaling = "std") + indegreeSender(scaling = "std")
    choice_effects <- ~ inertia(scaling = "std") + reciprocity(scaling = "std")

    stats_ao <- remstats(
      reh = reh_ao,
      sender_effects   = rate_effects,
      receiver_effects = choice_effects
    )
    stats_ao
    #> Relational Event Network Statistics
    #> > Model: actor-oriented
    #> > Computation method: per time point
    #> > Sender model:
    #>   >> Dimensions: 992 time points x 5 actors x 3 statistics
    #>   >> Statistics:
    #>       >>> 1: baseline
    #>       >>> 2: outdegreeSender
    #>       >>> 3: indegreeSender
    #> > Receiver model:
    #>   >> Dimensions: 992 events x 5 actors x 2 statistics
    #>   >> Statistics:
    #>       >>> 1: inertia
    #>       >>> 2: reciprocity

    fit_ao <- remstimate(reh = reh_ao, stats = stats_ao)
    summary(fit_ao)
    #> Relational Event Model (actor oriented)
    #>
    #> Call rate model **for sender**:
    #>
    #>  ~outdegreeSender(scaling = "std") + indegreeSender(scaling = "std")
    #>
    #>
    #> Coefficients rate model (MLE with interval likelihood):
    #>
    #>                    Estimate    Std. Err     z value Pr(>|z|) Pr(=0)
    #> baseline        -5.5381e+00  3.1703e-02 -1.7469e+02   0.0000 <2e-16
    #> outdegreeSender  4.9827e-03  3.6671e-02  1.3588e-01   0.8919 0.9690
    #> indegreeSender   6.0544e-02  3.6147e-02  1.6749e+00   0.0940 0.8857
    #> Null deviance: 13046.93 on 992 degrees of freedom
    #> Residual deviance: 13043.82 on 989 degrees of freedom
    #> Chi-square: 3.10366 on 3 degrees of freedom, asymptotic p-value 0.3759173
    #> AIC: 13049.82 AICC: 13049.85 BIC: 13064.52
    #> --------------------------------------------------------------------------------
    #>
    #> Call choice model **for receiver**:
    #>
    #>  ~inertia(scaling = "std") + reciprocity(scaling = "std")
    #>
    #>
    #> Coefficients choice model (MLE with interval likelihood):
    #>
    #>             Estimate Std. Err z value Pr(>|z|) Pr(=0)
    #> inertia      0.16805  0.10140 1.65730   0.0975 0.8886
    #> reciprocity  0.10580  0.10164 1.04091   0.2979 0.9482
    #> Null deviance: 2767.044 on 992 degrees of freedom
    #> Residual deviance: 2722.256 on 990 degrees of freedom
    #> Chi-square: 44.78752 on 2 degrees of freedom, asymptotic p-value 1.88154e-10
    #> AIC: 2726.256 AICC: 2726.268 BIC: 2736.055

The sender model estimates a baseline rate and the effects of out-degree and in-degree on sender activity. The receiver model describes how inertia and reciprocity shape the choice of receiver, conditional on the sender.

## 3.2 Actor-oriented model with exogenous effects

    rate_effects_exo   <- ~ outdegreeSender(scaling = "std") +
                            send("age", attr_actors = info3, scaling = "std")
    choice_effects_exo <- ~ inertia(scaling = "std") + reciprocity(scaling = "std") +
                            receive("age", attr_actors = info3, scaling = "std")

    stats_ao_exo <- remstats(
      reh = reh_ao,
      sender_effects   = rate_effects_exo,
      receiver_effects = choice_effects_exo
    )

    fit_ao_exo <- remstimate(reh = reh_ao, stats = stats_ao_exo)
    summary(fit_ao_exo)
    #> Relational Event Model (actor oriented)
    #>
    #> Call rate model **for sender**:
    #>
    #>  ~outdegreeSender(scaling = "std") + send("age", attr_actors = info3,     scaling = "std")
    #>
    #>
    #> Coefficients rate model (MLE with interval likelihood):
    #>
    #>                    Estimate    Std. Err     z value Pr(>|z|) Pr(=0)
    #> baseline          -5.536886    0.031666 -174.855347   0.0000 <2e-16
    #> outdegreeSender    0.042033    0.050524    0.831932   0.4054  0.957
    #> send              -0.031186    0.050342   -0.619494   0.5356  0.963
    #> Null deviance: 13046.93 on 992 degrees of freedom
    #> Residual deviance: 13046.23 on 989 degrees of freedom
    #> Chi-square: 0.695269 on 3 degrees of freedom, asymptotic p-value 0.8743162
    #> AIC: 13052.23 AICC: 13052.26 BIC: 13066.93
    #> --------------------------------------------------------------------------------
    #>
    #> Call choice model **for receiver**:
    #>
    #>  ~inertia(scaling = "std") + reciprocity(scaling = "std") + receive("age",     attr_actors = info3, scaling = "std")
    #>
    #>
    #> Coefficients choice model (MLE with interval likelihood):
    #>
    #>              Estimate  Std. Err   z value Pr(>|z|) Pr(=0)
    #> inertia      0.067447  0.108193  0.623391   0.5330 0.9629
    #> reciprocity  0.204916  0.109116  1.877969   0.0604 0.8438
    #> receive     -0.117852  0.043621 -2.701750   0.0069 0.4502
    #> Null deviance: 2767.044 on 992 degrees of freedom
    #> Residual deviance: 2714.937 on 989 degrees of freedom
    #> Chi-square: 52.10606 on 3 degrees of freedom, asymptotic p-value 2.843181e-11
    #> AIC: 2720.937 AICC: 2720.962 BIC: 2735.637

## 3.3 Actor-oriented diagnostics

    diag_ao <- diagnostics(object = fit_ao, reh = reh_ao, stats = stats_ao)
    plot(x = fit_ao, reh = reh_ao, diagnostics = diag_ao)

    #> Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
    #> collapsing to unique 'x' values

    #> Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
    #> collapsing to unique 'x' values

# 4 Bayesian estimation with HMC

Basic tie-oriented models can also be estimated by Hamiltonian Monte Carlo (HMC) via `approach = "Bayesian"`. Sampler settings are passed through the `bayes` list. Bayesian HMC is available for the basic tie (and actor) model only; for penalized models the Bayesian approach uses shrinkage priors instead (see the penalized-REM vignette).

    fit_hmc <- remstimate(
      reh      = reh,
      stats    = stats,
      approach = "Bayesian",
      bayes    = list(
        nsim    = 300L,
        nchains = 2L,
        burnin  = 300L,
        L       = 100L,
        epsilon = 0.001,
        thin    = 2L
      ),
      seed = 42
    )

    summary(fit_hmc)
    #> Relational Event Model (tie oriented)
    #>
    #> Call:
    #> ~inertia(scaling = "std") + reciprocity(scaling = "std") + outdegreeSender(scaling = "std") +     isp(scaling = "std")
    #>
    #>
    #> Posterior Modes (HMC with interval likelihood):
    #>
    #>                 Post.Mode   Post.SD     Q2.5%      Q50%  Q97.5%  Pr(=0|x)
    #> baseline        -7.251976  0.026655 -7.304215 -7.251976 -7.1974 < 2.2e-16
    #> inertia          0.422994  0.025578  0.376685  0.421789  0.4697 < 2.2e-16
    #> reciprocity      0.202351  0.023501  0.154995  0.201360  0.2447 2.512e-15
    #> outdegreeSender -0.117026  0.023227 -0.159971 -0.119142 -0.0829 9.678e-05
    #> isp             -0.358685  0.031968 -0.422702 -0.353837 -0.2972 < 2.2e-16
    #> Log posterior: -7613.877

    diag_hmc <- diagnostics(object = fit_hmc, reh = reh, stats = stats)
    plot(x = fit_hmc, reh = reh, diagnostics = diag_hmc)

    #> Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
    #> collapsing to unique 'x' values

The HMC output includes posterior means, standard deviations, and trace plots for convergence assessment. For most applied work the maximum-likelihood fit is sufficient; HMC is useful mainly when full posterior uncertainty is required.

# 5 Putting it all together

The table below summarizes the key modeling choices and how they map to function arguments:

Choice | Argument | Options
--- | --- | ---
Model type | `remify(model = ...)` | `"tie"`, `"actor"`
Timing | `remify(ordinal = ...)` | `FALSE` (interval), `TRUE` (ordinal)
Risk set | `remify(riskset = ...)` | `"full"`, `"active"`, `"active_saturated"`, `"manual"`
Event types | `remify(event_type = ...)` | column name or `NULL`
Type-extended risk set | `remify(extend_riskset_by_type = ...)` | `FALSE` (default), `TRUE`
Type handling | `inertia(consider_type = ...)` | `"ignore"`, `"separate"`, `"interact"`
Memory | `remstats(memory = ...)` | `"full"`, `"window"`, `"interval"`, `"decay"`
Sampling | `remstats(sampling = ...)` | `TRUE`/`FALSE`
Estimation approach | `remstimate(approach = ...)` | `"frequentist"` (default), `"Bayesian"`
Model structure | `remstimate(random=, penalty=, mixture=)` | GLMM / penalized / mixture

These choices can be freely combined — for example, an ordinal actor-oriented model with an active risk set and typed events using `consider_type = "separate"`. The `random`, `penalty`, and `mixture` arguments — covered in the frailty, penalized, and mixture vignettes — extend the basic model without changing the `remify()`/`remstats()` steps.