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.