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

  1. remify() (Arena et al. 2026a) — process a raw edgelist into a standardized event history, specifying the model type (“tie” or “actor”), risk set, timing, and event characteristics (e.g., types, directionality).
  2. remstats() (Meijerink-Bosman, Arena, and Leenders 2024) — compute network statistics (endogenous, exogenous, and interactions), with options for memory decay and case-control sampling.
  3. remstimate() (Arena et al. 2026b) — estimate model parameters via MLE or HMC. Model fit can then be assessed with diagnostics().

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)
#> Loading required package: remify
#> Loading required package: remstats
#> Loading required package: remstimate

1 Data

We use the edgelist0 dataset (115 directed events among 10 actors) and the edgelist0_actors dataset (time-varying actor attributes).

data(edgelist0)
data(edgelist0_actors)

head(edgelist0)
#>   time actor1 actor2 setting
#> 1 24.0      5      3    work
#> 2 25.5      5      1    work
#> 3 31.3      5      1    work
#> 4 31.5      2      4  social
#> 5 33.9      3      1  social
#> 6 35.9      4      5    work
head(edgelist0_actors)
#>   name time job
#> 1    1    0   1
#> 2    2    0   0
#> 3    3    0   1
#> 4    4    0   0
#> 5    5    0   0

Events record a time, a sender (actor1), a receiver (actor2), a setting (event type: "work" or "social"), and a weight. The edgelist0_actors table provides job scores measured from time point 0 for each actor.

2 Quick start: tie-oriented model

The tie-oriented model (model = "tie") models the event rate at the dyadic (tie) level (Butts 2008). The three-step pipeline processes the event history, computes network statistics, and estimates the model.

# Step 1: process the event history
reh <- remify(edgelist = edgelist0, model = "tie", directed = TRUE)
summary(reh)
#> Relational Event Network
#> (processed for tie-oriented modeling):
#>  > events = 1000 (time points = 995)
#>  > actors = 5
#>  > (event) types = 1
#>  > riskset = full
#>      >> included dyads = 20
#>  > directed = TRUE
#>  > ordinal = FALSE
#>  > weighted = FALSE
#>  > time length ~ 4027 
#>  > interevent time 
#>       >> minimum ~ 0.1 
#>       >> maximum ~ 26.3

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

# Step 3: estimate the model
fit <- remstimate(reh = reh, stats = stats, method = "MLE")
summary(fit)
#> Relational Event Model (tie oriented) 
#> 
#> Call:
#> ~inertia() + reciprocity() + indegreeSender() + outdegreeReceiver() +     outdegreeSender()
#> 
#> 
#> Coefficients (MLE with interval likelihood):
#> 
#>                      Estimate    Std. Err     z value Pr(>|z|)    Pr(=0)
#> baseline          -4.4272e+00  5.9052e-02 -7.4971e+01    0e+00 < 2.2e-16
#> inertia            1.1036e-02  6.3043e-04  1.7505e+01    0e+00 < 2.2e-16
#> reciprocity        4.3270e-03  9.1000e-04  4.7550e+00    0e+00  0.000388
#> indegreeSender    -1.6488e-03  5.3436e-04 -3.0856e+00    2e-03  0.212577
#> outdegreeReceiver -1.7421e-03  4.9155e-04 -3.5441e+00    4e-04  0.055760
#> outdegreeSender   -2.5312e-03  5.8484e-04 -4.3281e+00    0e+00  0.002691
#> Null deviance: 10756.99 on 994 degrees of freedom
#> Residual deviance: 9901.172 on 988 degrees of freedom
#> Chi-square: 855.8142 on 6 degrees of freedom, asymptotic p-value 0 
#> AIC: 9913.172 AICC: 9913.257 BIC: 9942.582

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 of dyads to repeat past interactions, reciprocity captures the tendency to reciprocate received events, and the degree effects capture sender and receiver activity, all as function of the volume of past events. These endogenous statistics can be scaled in different manners (Section 4.5) which also affects the interpretation of the corresponding effects. The column Pr(>|z|) is the 2-sided p value for testing whether a parameter equals zero. The column Pr(=0) is the posterior probability that a parameter is zero based on the default Bayes factor of the BFpack package (Mulder et al. 2021).

3 Quick start: actor-oriented model

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?) (Stadtfeld, Hollway, and Block 2017). Currently the actor-oriented model is only supported for directed events.

# Step 1: process the event history
reh_ao <- remify(edgelist = edgelist0, model = "actor", directed = TRUE, riskset = "active")
summary(reh_ao)
#> Relational Event Network
#> (processed for actor-oriented modeling):
#>  > events = 1000 (time points = 995)
#>  > actors = 5
#>  > (event) types = 1
#>  > riskset = active
#>  > sender model riskset: 5 / 5 actors
#>  > receiver model riskset: 4 receivers per sender
#>  > directed = TRUE
#>  > ordinal = FALSE
#>  > weighted = FALSE
#>  > time length ~ 4027 
#>  > interevent time 
#>       >> minimum ~ 0.1 
#>       >> maximum ~ 26.3

# Step 2: specify statistics for sender and receiver choice models separately
sender_effects <- ~ outdegreeSender() + indegreeSender()
choice_effects <- ~ inertia() + reciprocity()

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

# Step 3: estimate the model
fit_ao <- remstimate(reh = reh_ao, stats = stats_ao, method = "MLE")
summary(fit_ao)
#> Relational Event Model (actor oriented) 
#> 
#> Call rate model **for sender**:
#> 
#>  ~outdegreeSender() + indegreeSender()
#> 
#> 
#> Coefficients rate model (MLE with interval likelihood):
#> 
#>                    Estimate    Std. Err     z value Pr(>|z|)    Pr(=0)
#> baseline        -3.3295e+00  5.3856e-02 -6.1822e+01        0 < 2.2e-16
#> outdegreeSender  4.0300e-03  2.3258e-04  1.7327e+01        0 < 2.2e-16
#> indegreeSender  -1.7079e-03  3.9846e-04 -4.2862e+00        0  0.003221
#> Null deviance: 7987.17 on 994 degrees of freedom
#> Residual deviance: 7737.266 on 991 degrees of freedom
#> Chi-square: 249.9032 on 3 degrees of freedom, asymptotic p-value 0 
#> AIC: 7743.266 AICC: 7743.291 BIC: 7757.972
#> -------------------------------------------------------------------------------- 
#> 
#> Call choice model **for receiver**:
#> 
#>  ~inertia() + reciprocity()
#> 
#> 
#> Coefficients choice model (MLE with interval likelihood):
#> 
#>                Estimate    Std. Err     z value Pr(>|z|) Pr(=0)
#> inertia      0.00698851  0.00051304 13.62178835    0.000 <2e-16
#> reciprocity -0.00042811  0.00110705 -0.38670921    0.699 0.9669
#> Null deviance: 2769.816 on 994 degrees of freedom
#> Residual deviance: 2544.144 on 992 degrees of freedom
#> Chi-square: 225.6718 on 2 degrees of freedom, asymptotic p-value 0 
#> AIC: 2548.144 AICC: 2548.156 BIC: 2557.948

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

4 Modeling choices

The sections below describe the modeling choices available at each step of the pipeline. Most choices apply to both tie-oriented and actor-oriented models; exceptions are noted.

4.1 Timing: interval vs. ordinal

By default, the model uses exact event times (ordinal = FALSE), resulting in a rate model with a baseline intercept. When only the order of events is known or only the order is of interest, set ordinal = TRUE. The model then reduces to a stratified Cox proportional-hazard model — no baseline intercept can be estimated. This choice applies to both tie and actor models.

# Tie model with ordinal timing
reh_ord <- remify(edgelist = edgelist0, model = "tie", ordinal = TRUE)
stats_ord <- remstats(reh = reh_ord, tie_effects = ~ inertia() + reciprocity())
fit_ord <- remstimate(reh = reh_ord, stats = stats_ord, method = "MLE")
summary(fit_ord)
#> Relational Event Model (tie oriented) 
#> 
#> Call:
#> ~inertia() + reciprocity()
#> 
#> 
#> Coefficients (MLE with ordinal likelihood):
#> 
#>               Estimate   Std. Err    z value Pr(>|z|)    Pr(=0)
#> inertia     1.0491e-02 2.1418e-04 4.8980e+01        0 < 2.2e-16
#> reciprocity 2.9302e-03 5.3452e-04 5.4819e+00        0 9.401e-06
#> Null deviance: 5950.486 on 994 degrees of freedom
#> Residual deviance: 5021.903 on 992 degrees of freedom
#> Chi-square: 928.5823 on 2 degrees of freedom, asymptotic p-value 0 
#> AIC: 5025.903 AICC: 5025.915 BIC: 5035.707
# Actor model with ordinal timing
reh_ao_ord <- remify(edgelist = edgelist0, model = "actor", ordinal = TRUE, riskset = "active")
stats_ao_ord <- remstats(
  reh = reh_ao_ord,
  sender_effects = ~ outdegreeSender(),
  receiver_effects = ~ inertia() + reciprocity()
)
fit_ao_ord <- remstimate(reh = reh_ao_ord, stats = stats_ao_ord, method = "MLE")
summary(fit_ao_ord)
#> Relational Event Model (actor oriented) 
#> 
#> Call rate model **for sender**:
#> 
#>  ~outdegreeSender()
#> 
#> 
#> Coefficients rate model (MLE with ordinal likelihood):                  Estimate   Std. Err    z value Pr(>|z|)    Pr(=0)
#> outdegreeSender 6.0909e-03 3.2132e-04 1.8956e+01        0 < 2.2e-16
#> Null deviance: 3199.563 on 994 degrees of freedom
#> Residual deviance: 2840.317 on 993 degrees of freedom
#> Chi-square: 359.2458 on 1 degrees of freedom, asymptotic p-value 0 
#> AIC: 2842.317 AICC: 2842.321 BIC: 2847.219
#> -------------------------------------------------------------------------------- 
#> 
#> Call choice model **for receiver**:
#> 
#>  ~inertia() + reciprocity()
#> 
#> 
#> Coefficients choice model (MLE with ordinal likelihood):               Estimate    Std. Err     z value Pr(>|z|) Pr(=0)
#> inertia      0.00698851  0.00051304 13.62178835    0.000 <2e-16
#> reciprocity -0.00042811  0.00110705 -0.38670921    0.699 0.9669
#> Null deviance: 2755.953 on 994 degrees of freedom
#> Residual deviance: 2544.144 on 992 degrees of freedom
#> Chi-square: 211.8088 on 2 degrees of freedom, asymptotic p-value 0 
#> AIC: 2548.144 AICC: 2548.156 BIC: 2557.948

4.2 Directionality

When events have no inherent directionality, set directed = FALSE in remify(). The risk set then contains unordered actor pairs rather than ordered dyads, and only statistics defined for undirected networks are available. For example, sp() (shared partners) counts the number of third actors through which the two actors have interacted with in the past. Undirected events are only supported for the tie-oriented model.

reh_undir <- remify(edgelist = edgelist0, model = "tie", directed = FALSE)
cat("Directed dyads:", reh$D, "\n")
#> Directed dyads: 20
cat("Undirected pairs:", reh_undir$D, "\n")
#> Undirected pairs: 10

stats_undir <- remstats(
  reh = reh_undir,
  tie_effects = ~ inertia() + sp()
)

fit_undir <- remstimate(reh = reh_undir, stats = stats_undir, method = "MLE")
summary(fit_undir)
#> Relational Event Model (tie oriented) 
#> 
#> Call:
#> ~inertia() + sp()
#> 
#> 
#> Coefficients (MLE with interval likelihood):
#> 
#>             Estimate    Std. Err     z value Pr(>|z|)    Pr(=0)
#> baseline -3.6133e+00  6.2363e-02 -5.7939e+01        0 < 2.2e-16
#> inertia   8.0259e-03  3.0243e-04  2.6538e+01        0 < 2.2e-16
#> sp       -9.7984e-03  9.6258e-04 -1.0179e+01        0 < 2.2e-16
#> Null deviance: 9372.078 on 994 degrees of freedom
#> Residual deviance: 8679.267 on 991 degrees of freedom
#> Chi-square: 692.8109 on 3 degrees of freedom, asymptotic p-value 0 
#> AIC: 8685.267 AICC: 8685.291 BIC: 8699.972

The sp() statistic also supports scaling: "none" gives the raw shared partner count, while "std" standardizes across all pairs at each time point. Setting unique = TRUE counts only third actors that create a new shared partner (i.e., not previously connected to both actors in the pair).

4.3 Risk set variations

The risk set defines which dyads (or actors) are at risk of producing an event at each time point. This choice applies to both tie and actor models.

Full risk set (default)

All possible dyads are at risk at every time point. For a directed network with \(N\) actors, this gives \(N(N-1)\) dyads.

Active risk set

In certain situations, it may be unrealistic that all dyads in the network are at risk to initiate events (e.g., in the case of very large networks). In this case, the user could approximate the actual risk set by taking the active risk set which includes only dyads that are observed in the event history data. This may severely reduce the computation of endogenous statistics as only statistics are computed that are in this (typically much) smaller risk set.

reh_active <- remify(edgelist = edgelist0, 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() + reciprocity())
fit_active <- remstimate(reh = reh_active, stats = stats_active, method = "MLE")
summary(fit_active)
#> Relational Event Model (tie oriented) 
#> 
#> Call:
#> ~inertia() + reciprocity()
#> 
#> 
#> Coefficients (MLE with interval likelihood):
#> 
#>                Estimate    Std. Err     z value Pr(>|z|) Pr(=0)
#> baseline    -4.8094e+00  4.1773e-02 -1.1513e+02     0.00 <2e-16
#> inertia      8.3389e-03  2.3809e-04  3.5024e+01     0.00 <2e-16
#> reciprocity  5.0533e-04  6.1339e-04  8.2383e-01     0.41 0.9574
#> Null deviance: 10756.99 on 994 degrees of freedom
#> Residual deviance: 9973.755 on 991 degrees of freedom
#> Chi-square: 783.2304 on 3 degrees of freedom, asymptotic p-value 0 
#> AIC: 9979.755 AICC: 9979.78 BIC: 9994.461

In addition to the active risk set, the active_saturated risk set extends the observed dyads by including their reverse direction (e.g., if only A \(\rightarrow\) B is observed and B \(\rightarrow\) A is not observed, B \(\rightarrow\) A is also included) and, for typed events, all event types for each observed actor pair.

Manual risk set

A manual risk set allows you to restrict the set of dyads at risk to a predefined set. Observed dyads not included in the manual specification are added automatically.

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

reh_manual <- remify(
  edgelist       = edgelist0,
  model          = "tie",
  riskset        = "manual",
  manual.riskset = my_riskset
)
#> Warning in remify(edgelist = edgelist0, 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: 22 dyads

4.4 Typed events

When events carry a type label (here, "work" vs. "social" in the column setting of the edgelist), different choices can be made how this can be modeled. When a column of the edgelist contains a column type, the event types are automatically taken into account when processing (remifying) the data. Alternatively, one can use the argument event_type to specify which column of the edgelist contains the event types. The risk set can optionally be extended by type (extend_riskset_by_type = TRUE; default is FALSE), so that each dyad-type combination is a separate unit at risk. This is only possible in the tie-oriented model.

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

When computing the endogenous statistics, the consider_type argument controls how past event types enter the computation.

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. This allows one to model separate effects for the separate event types.

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"

With consider_type = "interact" (requires extend_riskset_by_type = TRUE), each base statistic is expanded into \(C^2\) statistics, one for each combination of past event type and riskset element type. For example, with types work and social, inertia(consider_type = "interact") produces four statistics: inertia.work.work, inertia.work.social, inertia.social.work, and inertia.social.social. This allows the model to estimate a separate effect for each past-type \(\times\) future-type combination.

reh_typed_ext <- remify(
  edgelist   = edgelist0,
  model      = "tie",
  event_type = "setting",
  extend_riskset_by_type = TRUE
)
stats_int <- remstats(
  reh = reh_typed_ext,
  tie_effects = ~ inertia(consider_type = "interact")
)
dimnames(stats_int)[[3]]
#> [1] "baseline"              "inertia.social.social" "inertia.social.work"  
#> [4] "inertia.work.social"   "inertia.work.work"

Estimation with typed events

Different consider_type settings can be freely mixed within the same model.

stats_typed <- remstats(
  reh = reh_typed_ext,
  tie_effects = ~ inertia(consider_type = "interact") +
                   reciprocity(consider_type = "ignore")
)
fit_typed <- remstimate(reh = reh_typed_ext, stats = stats_typed, method = "MLE")
summary(fit_typed)
#> Relational Event Model (tie oriented) 
#> 
#> Call:
#> ~inertia(consider_type = "interact") + reciprocity(consider_type = "ignore")
#> 
#> 
#> Coefficients (MLE with interval likelihood):
#> 
#>                          Estimate    Std. Err     z value Pr(>|z|)    Pr(=0)
#> baseline              -5.4926e+00  4.1924e-02 -1.3101e+02   0.0000 < 2.2e-16
#> inertia.social.social  5.4267e-02  8.5555e-03  6.3430e+00   0.0000 5.784e-08
#> inertia.social.work    5.9420e-03  7.3931e-03  8.0372e-01   0.4216   0.95803
#> inertia.work.social   -1.7874e-02  4.5098e-03 -3.9633e+00   0.0001   0.01209
#> inertia.work.work      1.1122e-02  3.9017e-03  2.8506e+00   0.0044   0.35159
#> reciprocity           -1.3507e-04  6.2320e-04 -2.1674e-01   0.8284   0.96855
#> Null deviance: 12141.89 on 994 degrees of freedom
#> Residual deviance: 11292.97 on 988 degrees of freedom
#> Chi-square: 848.9196 on 6 degrees of freedom, asymptotic p-value 0 
#> AIC: 11304.97 AICC: 11305.06 BIC: 11334.38

The type-separated coefficients of inertia tell us whether inertia operates differently across event types. For instance, a strong negative effect of inertia.work.social suggests that past work interactions tend to result in fewer social events.

Typed events and consider_type work the same way in the actor-oriented model:

reh_ao_typed <- remify(
  edgelist   = edgelist0,
  model      = "actor",
  event_type = "setting"
)
stats_ao_typed <- remstats(
  reh = reh_ao_typed,
  sender_effects = ~ outdegreeSender(),
  receiver_effects = ~ inertia(consider_type = "separate") + reciprocity()
)
fit_ao_typed <- remstimate(reh = reh_ao_typed, stats = stats_ao_typed, method = "MLE")
summary(fit_ao_typed)
#> Relational Event Model (actor oriented) 
#> 
#> Call rate model **for sender**:
#> 
#>  ~outdegreeSender()
#> 
#> 
#> Coefficients rate model (MLE with interval likelihood):
#> 
#>                    Estimate    Std. Err     z value Pr(>|z|)    Pr(=0)
#> baseline        -3.4596e+00  4.6876e-02 -7.3805e+01        0 < 2.2e-16
#> outdegreeSender  3.7031e-03  2.2214e-04  1.6670e+01        0 < 2.2e-16
#> Null deviance: 7987.17 on 994 degrees of freedom
#> Residual deviance: 7758.331 on 992 degrees of freedom
#> Chi-square: 228.8389 on 2 degrees of freedom, asymptotic p-value 0 
#> AIC: 7762.331 AICC: 7762.343 BIC: 7772.134
#> -------------------------------------------------------------------------------- 
#> 
#> Call choice model **for receiver**:
#> 
#>  ~inertia(consider_type = "separate") + reciprocity()
#> 
#> 
#> Coefficients choice model (MLE with interval likelihood):
#> 
#>                   Estimate    Std. Err     z value Pr(>|z|) Pr(=0)
#> inertia.social  0.00578609  0.00757659  0.76367959   0.4451 0.9593
#> inertia.work    0.00760241  0.00389432  1.95217856   0.0509 0.8242
#> reciprocity    -0.00039993  0.00112302 -0.35612416   0.7217 0.9673
#> Null deviance: 2769.816 on 994 degrees of freedom
#> Residual deviance: 2544.119 on 991 degrees of freedom
#> Chi-square: 225.697 on 3 degrees of freedom, asymptotic p-value 0 
#> AIC: 2550.119 AICC: 2550.143 BIC: 2564.824

4.5 Endogenous effects, exogenous effects, and scaling

Endogenous effects

The full list of available statistics is obtained with tie_effects() for the tie model or actor_effects() for the actor model. Many endogenous statistics accept a scaling argument that controls how raw counts are transformed:

  • "none" (default): raw counts are used as-is. For inertia(), this is the number of past events on the dyad.
  • "prop": raw counts are divided by a normalizing quantity. For inertia(), the count is divided by the sender’s outdegree at time \(t\), yielding the proportion of the sender’s past events directed at this receiver.
  • "std": raw counts are standardized (zero mean, unit variance) across all dyads at each time point.
reh <- remify(edgelist = edgelist0, model = "tie", directed = TRUE)

stats_none <- remstats(reh = reh, tie_effects = ~ inertia(scaling = "none"))
stats_prop <- remstats(reh = reh, tie_effects = ~ inertia(scaling = "prop"))
stats_std  <- remstats(reh = reh, tie_effects = ~ inertia(scaling = "std"))

# Compare values at the last time point for the first dyad
cat("Raw count:      ", stats_none[114, 1, "inertia"], "\n")
#> Raw count:       2
cat("Proportional:   ", stats_prop[114, 1, "inertia"], "\n")
#> Proportional:    0.25
cat("Standardized:   ", stats_std[114, 1, "inertia"], "\n")
#> Standardized:    -0.5079315

remst_none <- remstimate(reh, stats = stats_none)
coef(remst_none)
#>     baseline      inertia 
#> -4.796843530  0.008384646
remst_prop <- remstimate(reh, stats = stats_prop)
coef(remst_prop)
#>  baseline   inertia 
#> -5.759835  3.920526
remst_std <- remstimate(reh, stats = stats_std)
coef(remst_std)
#>   baseline    inertia 
#> -4.7560614  0.6087506

The different types of scaling change the interpretation of the corresponding effects. Depending on the application, different types of scaling may be preferred. Caution is advised when using raw counts (scaling = "none") because the statistics will then be unbounded, potentially resulting in process explosion.

Exogenous effects

Exogenous actor attributes and dyadic covariates can be included in the model alongside endogenous effects. The edgelist0_actors dataset provides time-varying attributes.

exo_effects <- ~ inertia(scaling = "std") +
                   reciprocity(scaling = "std") +
                   send("job", edgelist0_actors) +
                   same("job", edgelist0_actors)
stats_exo <- remstats(
  reh = reh,
  tie_effects = exo_effects
)
fit_exo <- remstimate(reh = reh, stats = stats_exo, method = "MLE")
summary(fit_exo)
#> Relational Event Model (tie oriented) 
#> 
#> Call:
#> ~inertia(scaling = "std") + reciprocity(scaling = "std") + send("job",     edgelist0_actors) + same("job", edgelist0_actors)
#> 
#> 
#> Coefficients (MLE with interval likelihood):
#> 
#>               Estimate   Std. Err    z value Pr(>|z|)    Pr(=0)
#> baseline     -4.887699   0.071910 -67.969816   0.0000 < 2.2e-16
#> inertia       0.537229   0.019781  27.159057   0.0000 < 2.2e-16
#> reciprocity   0.074080   0.033341   2.221905   0.0263  0.727592
#> send_job     -0.210875   0.093617  -2.252518   0.0243  0.713808
#> same_job      0.415294   0.092020   4.513086   0.0000  0.001189
#> Null deviance: 10756.99 on 994 degrees of freedom
#> Residual deviance: 9483.757 on 989 degrees of freedom
#> Chi-square: 1273.228 on 5 degrees of freedom, asymptotic p-value 0 
#> AIC: 9493.757 AICC: 9493.818 BIC: 9518.266
  • send("job"): quantifies whether actors who have a job tend to initiate more events.
  • same("job"): quantifies whether actors having the same job status tend to send more events to each other.

Exogenous effects are also available in the actor model, where they can enter the sender rate model, the receiver choice model, or both:

rate_effects_exo <- ~ outdegreeSender(scaling = "std") + send("job", edgelist0_actors)
choice_effects_exo <- ~ inertia(scaling = "std") + reciprocity(scaling = "std") + same("job", edgelist0_actors)

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, method = "MLE")
summary(fit_ao_exo)
#> Relational Event Model (actor oriented) 
#> 
#> Call rate model **for sender**:
#> 
#>  ~outdegreeSender(scaling = "std") + send("job", edgelist0_actors)
#> 
#> 
#> Coefficients rate model (MLE with interval likelihood):
#> 
#>                   Estimate   Std. Err    z value Pr(>|z|) Pr(=0)
#> baseline         -3.087328   0.049494 -62.378099   0.0000 <2e-16
#> outdegreeSender   0.578656   0.036976  15.649602   0.0000 <2e-16
#> send             -0.318084   0.098880  -3.216858   0.0013 0.1515
#> Null deviance: 7987.17 on 994 degrees of freedom
#> Residual deviance: 7525.149 on 991 degrees of freedom
#> Chi-square: 462.0205 on 3 degrees of freedom, asymptotic p-value 0 
#> AIC: 7531.149 AICC: 7531.173 BIC: 7545.854
#> -------------------------------------------------------------------------------- 
#> 
#> Call choice model **for receiver**:
#> 
#>  ~inertia(scaling = "std") + reciprocity(scaling = "std") + same("job",     edgelist0_actors)
#> 
#> 
#> Coefficients choice model (MLE with interval likelihood):
#> 
#>              Estimate  Std. Err   z value Pr(>|z|) Pr(=0)
#> inertia      0.590418  0.045437 12.994315   0.0000 <2e-16
#> reciprocity  0.090911  0.049404  1.840135   0.0657 0.8529
#> same         0.239312  0.092492  2.587389   0.0097 0.5259
#> Null deviance: 2769.816 on 994 degrees of freedom
#> Residual deviance: 2448.297 on 991 degrees of freedom
#> Chi-square: 321.5196 on 3 degrees of freedom, asymptotic p-value 0 
#> AIC: 2454.297 AICC: 2454.321 BIC: 2469.002

Interactions between endogenous and exogenous effects

Interaction effects can be used to examine how endogenous tendencies (e.g., inertia) depend on actors’ attributes:

stats_ix <- remstats(
  reh = reh,
  tie_effects = ~ inertia() +
                   send("job", edgelist0_actors) +
                   inertia():send("job", edgelist0_actors)
)

fit_ix <- remstimate(reh = reh, stats = stats_ix, method = "MLE")
summary(fit_ix)
#> Relational Event Model (tie oriented) 
#> 
#> Call:
#> ~inertia() + send("job", edgelist0_actors) + inertia():send("job",     edgelist0_actors)
#> 
#> 
#> Coefficients (MLE with interval likelihood):
#> 
#>                     Estimate    Std. Err     z value Pr(>|z|)    Pr(=0)
#> baseline         -4.5679e+00  4.5087e-02 -1.0131e+02   0.0000 < 2.2e-16
#> inertia           7.6668e-03  2.4565e-04  3.1210e+01   0.0000 < 2.2e-16
#> send_job         -5.9210e-01  1.3868e-01 -4.2694e+00   0.0000  0.003461
#> inertia:send_job -6.2953e-03  8.8552e-03 -7.1092e-01   0.4771  0.960765
#> Null deviance: 10756.99 on 994 degrees of freedom
#> Residual deviance: 9905.106 on 990 degrees of freedom
#> Chi-square: 851.88 on 4 degrees of freedom, asymptotic p-value 0 
#> AIC: 9913.106 AICC: 9913.146 BIC: 9932.713

The interaction between inertia and send_job quantifies whether the tendency to repeat past interactions differs depending on the sender’s job.

4.6 Memory

The memory argument of remstats controls how the transpired time of past events contributes to endogenous statistics at each time point (Arena, Mulder, and Leenders 2023). This applies to both tie and actor models.

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

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

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

# Compare the inertia statistic at the last time point for a single dyad
cat("Full memory:  ", stats_full[900, 2, "inertia"], "\n")
#> Full memory:   19
cat("Window (500): ", stats_win[900, 2, "inertia"], "\n")
#> Window (500):  3
cat("Decay (200):  ", stats_dec[900, 2, "inertia"], "\n")
#> Decay (200):   1.453096

full memory accumulates all past events equally; window memory discards events older than the threshold (specified in memory_value); decay memory down-weights older events exponentially via a half-life parameter (specified in memory_value).

4.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, substantially reducing computation (Lerner and Lomi 2020). This can be done by setting sampling = TRUE and setting samp_num to be the total number of dyads in the risk set per event (which includes the observed event(s)). Case-control sampling is currently only available for the tie-oriented model.

reh_samp <- remify(edgelist = edgelist0, model = "tie", directed = TRUE, event_type = "setting")
stats_sampled <- remstats(
  reh = reh_samp,
  tie_effects = exo_effects,
  sampling = TRUE,
  samp_num = 20,
  seed = 42
)
fit_sampled <- remstimate(reh = reh_samp, stats = stats_sampled, method = "MLE")
summary(fit_sampled)
#> Relational Event Model (tie oriented) 
#> 
#> Call:
#> ~inertia(scaling = "std") + reciprocity(scaling = "std") + send("job",     edgelist0_actors) + same("job", edgelist0_actors)
#> 
#> 
#> Coefficients (MLE with interval likelihood):
#> 
#>               Estimate   Std. Err    z value Pr(>|z|)    Pr(=0)
#> baseline     -4.887819   0.071908 -67.973202   0.0000 < 2.2e-16
#> inertia       0.523009   0.019288  27.116047   0.0000 < 2.2e-16
#> reciprocity   0.072035   0.032504   2.216209   0.0267  0.730089
#> send_job     -0.211061   0.093615  -2.254562   0.0242  0.712866
#> same_job      0.415577   0.092018   4.516260   0.0000  0.001173
#> Null deviance: 10748.22 on 994 degrees of freedom
#> Residual deviance: 9479.315 on 989 degrees of freedom
#> Chi-square: 1268.902 on 5 degrees of freedom, asymptotic p-value 0 
#> AIC: 9489.315 AICC: 9489.375 BIC: 9513.823

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

5 Estimation

5.1 Maximum likelihood estimation

MLE is the default estimation method (method = "MLE") and has been used throughout this vignette. For interval timing, the model uses the full likelihood; for ordinal timing, a stratified partial likelihood (equivalent to a Cox proportional-hazard model). With case-control sampling, the case-control likelihood is used.

5.2 Bayesian estimation with HMC

Both tie-oriented and actor-oriented models can be estimated with Hamiltonian Monte Carlo via method = "HMC". By default, independent normal priors with mean 0 and variance 100 are placed on all coefficients. Custom priors can be specified via the prior argument. When using va:

P <- length(dimnames(stats_exo)[[3]])
my_prior <- list(
  mean = rep(0, P),
  vcov = diag(c(100, rep(1, P - 1)))  # wide prior on intercept, N(0,1) on effects
)

fit_exo_hmc <- remstimate(
  reh    = reh,
  stats  = stats_exo,
  method = "HMC",
  prior  = my_prior,
  burnin = 200,
  nsim   = 500,
  thin   = 2,
  L      = 100,
  seed   = 42
)

summary(fit_exo_hmc)
#> Relational Event Model (tie oriented) 
#> 
#> Call:
#> ~inertia(scaling = "std") + reciprocity(scaling = "std") + send("job",     edgelist0_actors) + same("job", edgelist0_actors)
#> 
#> 
#> Posterior Modes (HMC with interval likelihood):
#> 
#>             Post.Mode   Post.SD     Q2.5%      Q50%  Q97.5%  Pr(=0|x)
#> baseline    -4.889010  0.037965 -4.968681 -4.882815 -4.8182 < 2.2e-16
#> inertia      0.536315  0.011427  0.513806  0.537613  0.5581 < 2.2e-16
#> reciprocity  0.067568  0.017594  0.037355  0.077854  0.1051 0.0193956
#> send_job    -0.218039  0.047728 -0.303379 -0.206877 -0.1274 0.0009255
#> same_job     0.411246  0.052798  0.310082  0.406858  0.5259 2.111e-12
#> Log posterior: -4742.312

The HMC output includes posterior means, standard deviations, and trace plots for convergence assessment. Note that the algorithm takes considerably longer than MLE.

6 Diagnostics

The diagnostics() function computes several goodness-of-fit measures. The plot() function produces three diagnostic plots by default (which = 1:3).

6.1 Tie model diagnostics

diag <- diagnostics(object = fit_exo, reh = reh, stats = stats_exo)
diag
#> Diagnostics for a Relational Event Model
#> ------------------------------------------
#> Model       : tie
#> Actors      : 5
#> Events      : 994
#> 
#> Tie model
#>   Statistics : inertia, reciprocity, send_job, same_job 
#>   Recall     : mean rank = 0.717  |  top 5% = 41.1%
plot(x = fit_exo, reh = reh, diagnostics = diag)

Waiting-time diagnostics (plot 1)

Under a correctly specified model, the rescaled inter-event times \(\Delta t_m \cdot \sum_{(i,j) \in \mathcal{R}} \hat{\lambda}_{ij}(t_m)\) should follow an Exp(1) distribution. The Q-Q plot (left) compares the empirical quantiles against the theoretical Exp(1) quantiles: points falling on the diagonal indicate good fit. Systematic departures reveal specific types of misfit — a convex curve (points above the line) suggests the model underestimates waiting times for large gaps, while a concave curve suggests the opposite. The density histogram (right) overlays the empirical distribution with the Exp(1) density (dashed line) for a visual comparison.

These plots are only produced for interval timing models (ordinal = FALSE).

Schoenfeld residuals (plot 2)

For each effect included in the model, a standardized Schoenfeld residual is computed at every time point (Schoenfeld 1982). The residual measures the difference between the observed covariate value and its expected value under the model at that time point. Under a correctly specified model, the residuals should fluctuate randomly around zero with no systematic trend.

A smoothing spline (red line) is fitted through the residuals. If the spline substantially deviates from the zero line (dashed), the effect of that statistic may be changing over time — violating the proportional-hazards assumption that the coefficient is constant throughout the observation period.

Regarding static covariates (e.g., sender attributes such as send_job and same_job), the Schoenfeld residuals do not reveal a trend in the effect. However, these residuals are still informative: if residuals for a static covariate are systematically shifted away from zero, this suggests the estimated coefficient may be unreliable — possibly due to confounding with other effects or because the variable does not contribute meaningfully to the model. In contrast, residuals centered around zero indicate that the model’s predictions are consistent with the observed covariate values across events. In the current example, the residuals for same_job are centered around zero, while those for send_job are systematically shifted — suggesting the latter may be redundant or confounded with other effects in the model.

Recall (plot 3)

For each time point, the observed event is ranked among all competing dyads based on the estimated rate from the fitted model. The relative rank is plotted over time, where 1 = top-ranked (perfectly predicted) and 0 = bottom-ranked. A well-fitting model concentrates points near the top.

Three reference lines aid interpretation:

  • Null model (dotted grey at 0.5): the expected rank under an intercept-only model where all dyads are equally likely. Points above this line indicate the model does better than chance.
  • Median rank (dashed black): summarizes overall model performance. A median rank well above 0.5 indicates the model captures meaningful structure in the event process.
  • Top % threshold (solid blue at 0.95): points above this line are events the model ranked in the top 5% of all competing dyads.

The title reports the median rank and the proportion of events in the top 5%.

In this example, we see that the median recall is 0.85 implying that half of the observed events belong to the 15% most likely events according to the fitted model, suggesting a fairly good in-sample predictive performance.

6.2 Actor model diagnostics

The same diagnostic tools are available for actor-oriented models:

diag_ao <- diagnostics(object = fit_ao_exo, reh = reh_ao, stats = stats_ao_exo)
plot(x = fit_ao_exo, reh = reh_ao, diagnostics = diag_ao)

6.3 HMC diagnostics

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

7 Summary

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)
Directionality remify(directed = ...) TRUE, FALSE (tie model only)
Risk set remify(riskset = ...) "full", "active", "active_saturated", "manual"
Event types remify(event_type = ...) column name or NULL
Type handling inertia(consider_type = ...) "ignore", "separate", "interact"
Scaling inertia(scaling = ...) "none", "prop", "std"
Memory remstats(memory = ...) "full", "window", "interval", "decay"
Sampling remstats(sampling = ...) TRUE/FALSE (tie model only)
Estimation remstimate(method = ...) "MLE", "HMC"

These choices can be freely combined — for example, an ordinal actor-oriented model with an active risk set, typed events using consider_type = "separate", and Bayesian estimation via HMC.

Arena, Giuseppe, Rumana Lakdawala, Marlyne Meijerink-Bosman, Diana Karimova, Fabio Generoso Vieira, Mahdi Shafiee Kamalabad, Roger Th. A. J. Leenders, and Joris Mulder. 2026a. Remify: Processing and Transforming Relational Event History Data. https://CRAN.R-project.org/package=remify.
———. 2026b. Remstimate: Remstimate: Optimization Frameworks for Tie-Oriented and Actor-Oriented Relational Event Models. https://CRAN.R-project.org/package=remstimate.
Arena, Giuseppe, Joris Mulder, and Roger Th AJ Leenders. 2023. “How Fast Do We Forget Our Past Social Interactions? Understanding Memory Retention with Parametric Decays in Relational Event Models.” Network Science 11 (2): 267–94. https://doi.org/10.1017/nws.2023.5.
Butts, Carter T. 2008. “A Relational Event Framework for Social Action.” Sociological Methodology 38 (1): 155–200. https://doi.org/10.1111/j.1467-9531.2008.00203.x.
Lerner, Jürgen, and Alessandro Lomi. 2020. “Reliability of Relational Event Model Estimates Under Sampling: How to Fit a Relational Event Model to 360 Million Dyadic Events.” Network Science 8 (1): 97–135. https://doi.org/10.1017/nws.2019.57.
Meijerink-Bosman, Marlyne, Joris Arena Giuseppe and"; Mulder, and Roger Th. A. J. Leenders. 2024. “Discovering Trends of Social Interaction Behavior over Time: An Introduction to Relational Event Modeling.” Behavior Research Methods 56: 997–1023. https://doi.org/10.3758/s13428-022-01821-8.
Mulder, Joris, Donald R Williams, Xin Gu, Andrew Tomarken, Florian Böing-Messing, Anton Olsson-Collentine, Marlyne Meijerink-Bosman, et al. 2021. “BFpack: Flexible Bayes Factor Testing of Scientific Theories in r.” Journal of Statistical Software 100: 1–63. https://doi.org/10.18637/jss.v100.i18.
Schoenfeld, David. 1982. “Partial Residuals for the Proportional Hazards Regression Model.” Biometrika 69 (1): 239–41. https://doi.org/10.1093/biomet/69.1.239.
Stadtfeld, Christoph, James Hollway, and Per Block. 2017. “Dynamic Network Actor Models: Investigating Coordination Ties Through Time.” Sociological Methodology 47 (1): 1–40. https://doi.org/10.1177/0081175017709295.