# Duration Relational Event Models with remverse

### A tutorial

####

A Duration Relational Event Model (DuREM) extends the standard relational event model (REM) to events that have both a start time and an end time. Instead of a single process governing when the next event fires, there are two linked sub-models:

- a start model, which describes which inactive dyad initiates a new event and when;
- an end model, which describes which currently active dyad terminates and when.

The two sub-models are estimated jointly. Statistics for each sub-model can draw on (a) the accumulated event history — duration-weighted versions of standard remstats effects such as `inertia()` or `reciprocity()` — and (b) the current active state of the network, i.e. the set of events that have started but not yet ended. These active-state statistics can only be used for relational events with duration.

The remverse DuREM pipeline follows the same three-step pattern as a standard REM:

    remify(edgelist, duration = TRUE, ...)   →  remify_durem object
      remstats(reh, start_effects = ..., end_effects = ...)  →  remstats_durem object
        remstimate(reh, stats, ...)   →  remstimate_durem object

This vignette walks through each step using the `randomREH3` dataset provided with `remverse`.

# 1 Data

## 1.1 The `randomREH3` dataset

`randomREH3` is a simulated sequence of 999 directed duration events among 5 actors. Each row records one event with an observed start time and end time:

Column | Description
--- | ---
`time` | Start time of the event
`actor1` | Initiating actor (integer ID)
`actor2` | Receiving actor (integer ID)
`end` | End time of the event (`NA` for right-censored events)
`setting` | Context of the interaction: `"work"` or `"social"`
`duration` | Length of the event (`end - time`)

    library(remverse)   # loads remify, remstats, remstimate; remdata is a dependency
    #> Loading required package: remdata

    data(randomREH3)
    data(info3)         # actor attributes (age, sex)
    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

Although the data contains both an `end` column and a `duration`, only one of these columns (with these exact names) is necessary for a DuREM analysis.

The data were generated from a known duration model: starting an event depended on inertia, reciprocity, sender out-degree, incoming shared partners, and sender age; ending an event depended on the dyad’s total active degree and the age difference of the two actors.

    cat(sprintf(
        "%d events: %d complete, %d right-censored\n",
        nrow(randomREH3),
        sum(!is.na(randomREH3$end)),
        sum( is.na(randomREH3$end))
    ))
    #> 999 events: 999 complete, 0 right-censored

This particular dataset has no right-censored events (every event has an observed end), but the pipeline fully supports censoring (`end = NA`); see §5.1.

# 2 Step 1 — Processing with `remify()`

Pass `duration = TRUE` to tell remify that the edgelist contains duration information. DuREM analysis are only possible under a `tie` oriented model. The function returns an object of class `c("remify_durem", "remify")`.

    reh <- remify(
        edgelist  = randomREH3,
        duration  = TRUE,
        model     = "tie",
        directed  = TRUE        # actor1 initiates; the end process is undirected by default
    )
    #> Warning: 2 events start while another event for the same dyad is already
    #> active. Check events: 320, 565.

## 2.1 The `remify_durem` object

    reh
    #> Relational Event Network
    #> (processed for tie-oriented modeling with duration):
    #>  > 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
    #>       >> maximum ~ 448
    #>  > duration
    #>      >> start = directed (actor1 initiates)
    #>      >> end   = undirected (combined dyad-level rate)
    #>      >> event duration (complete events)
    #>           >>> minimum ~ 1
    #>           >>> median  ~ 12
    #>           >>> maximum ~ 326

    summary(reh)
    #> Relational Event Network
    #> (processed for tie-oriented modeling with duration):
    #>  > 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
    #>       >> maximum ~ 448
    #>  > duration
    #>      >> start = directed (actor1 initiates)
    #>      >> end   = undirected (combined dyad-level rate)
    #>      >> event duration (complete events)
    #>           >>> minimum ~ 1
    #>           >>> median  ~ 12
    #>           >>> maximum ~ 326

The `plot()` method for a `remify_durem` object shows the duration-aware descriptives — including the distribution of event durations and how many events are active over time — in addition to the standard relational event descriptives of a `remify` object.

    plot(reh)

Alongside the standard `remify` slots (`$M`, `$N`, `$D`, `$edgelist`, `$meta`, …), a `remify_durem` object carries DuREM-specific elements.

`$durem` — a named list of DuREM metadata:

    reh$durem
    #> $dur_directed_end
    #> [1] FALSE
    #>
    #> $dur_type_exclusive
    #> [1] FALSE
    #>
    #> $has_who_ended
    #> [1] FALSE
    #>
    #> $has_censored
    #> [1] FALSE
    #>
    #> $n_complete
    #> [1] 999
    #>
    #> $n_censored
    #> [1] 0

Field | Meaning
--- | ---
`n_complete` | Events with an observed end time
`n_censored` | Right-censored events (`end = NA`)
`has_censored` | Whether any events are censored
`dur_directed_end` | Whether the end process is directed
`has_who_ended` | Whether a `who_ended` column is present
`dur_type_exclusive` | Whether being active in one type blocks other types

`$edgelist_dual` — the dual-event expansion: one `"start"` row and one `"end"` row per complete event (censored events contribute only a `"start"` row), sorted by time. The `status` column states whether the observation at that time point is either the `start` of a new event or an `end` of an active event. The DuREM models these observed outcomes jointly by specifying separate sub-models for starting an event and for ending an event.

    head(reh$edgelist_dual, 12)
    #>    time actor1 actor2 status duration .eidx
    #> 1    18      5      1  start        6     1
    #> 2    24      5      1    end        6     1
    #> 3    57      3      2  start       13     2
    #> 4    64      1      5  start       21     3
    #> 5    70      3      2    end       13     2
    #> 6    85      1      5    end       21     3
    #> 7   127      5      1  start       32     4
    #> 8   159      5      1    end       32     4
    #> 9   228      1      4  start       11     5
    #> 10  239      1      4    end       11     5
    #> 11  274      2      4  start       46     6
    #> 12  302      1      5  start       12     7

# 3 Step 2 — Statistics with `remstats()`

For duration models, `remstats()` accepts two formulas — `start_effects` and `end_effects` — one per sub-model. Each formula may mix two kinds of terms:

1. History-weighted effects — the familiar tie-oriented effects (`inertia()`, `reciprocity()`, `totaldegreeDyad()`, `send()`, `difference()`, …), computed on the duration-weighted event history. An overview can be shown when running .
2. Active-state effects — effects that describe the currently active network (events that have started but not yet ended). These are the `active*()` family and are unique to the duration model. An overview can be shown when running .

`remstats()` inspects each term, routes it to the appropriate backend, and combines the results together automatically.

## 3.1 Active-state effect functions

Which active-state effects are available depends on whether the sub-model is directed. The start model is directed whenever `directed = TRUE`; the end model is undirected by default (set `dur_directed_end = TRUE` in `remify()` to make it directed).

Directed active-state effects (start model, or end model when `dur_directed_end = TRUE`):

Effect | Description
--- | ---
`activeTie()` | Is there currently an active event from actor \(i\) to actor \(j\)?
`activeOutdegreeSender()` | Number of active events actor \(i\) (sender) currently has as initiator
`activeIndegreeReceiver()` | Number of active events actor \(j\) (receiver) currently has as target
`activeTotaldegreeSender()` | Total active degree of actor \(i\) (as sender or receiver)
`activeTotaldegreeReceiver()` | Total active degree of actor \(j\)
`activeTotaldegreeDyad()` | Combined active degree of both actors
`activeReciprocalTie()` | Is the reverse dyad \(j \to i\) currently active?
`activeSharedPartners_otp()` / `_itp()` / `_osp()` / `_isp()` | Active shared-partner effects (outgoing/incoming two-path, outgoing/incoming shared partner)

Undirected active-state effects (end model, default):

Effect | Description
--- | ---
`activeTie()` | Is the dyad currently active?
`activeDegreeMin()` | Minimum active degree of the two actors
`activeDegreeMax()` | Maximum active degree of the two actors
`activeDegreeDyad()` | Combined active degree of the two actors
`activeSharedPartners()` | Number of shared active partners

All active-state effects share the optional arguments `scaling = c("none", "std")` and `consider_type = c("ignore", "separate", "interact")` (see §3.3).

> Note. A directed active-state effect such as `activeOutdegreeSender()` cannot be used in the (undirected) end model — it would be rejected as an unknown effect for that sub-model. Put directed active-state effects in the directed start model, and use the undirected active-state effects (or history-weighted effects) in the default end model.

## 3.2 Computing statistics

Here we formulate a DuREM where the start model mixes history-weighted effects with a directed active-state effect, while the (undirected) end model uses history-weighted effects only.

    # Start model (directed): history-weighted effects + active out-degree of sender
    start_fx <- ~ inertia(scaling = "std") +
                  reciprocity(scaling = "std") +
                  activeOutdegreeSender(scaling = "std")

    # End model (undirected): history-weighted effects
    end_fx   <- ~ totaldegreeDyad(scaling = "std") +
                  difference("age", attr_actors = info3, scaling = "std")

    dstats <- remstats(
        reh              = reh,
        start_effects    = start_fx,
        end_effects      = end_fx,
        memory           = "decay",
        memory_value     = 2000,
        display_progress = FALSE
    )

    dstats
    #> Relational Event Network Statistics
    #> > Model: tie-oriented (duration)
    #> > Computation method: per time point
    #> > Start dimensions: 1957 time points x 20 dyads x 4 statistics
    #> > Start statistics:
    #>   >> 1: baseline.start
    #>   >> 2: activeOutdegreeSender.start
    #>   >> 3: inertia.start
    #>   >> 4: reciprocity.start
    #> > End dimensions: 1957 time points x 10 dyads x 3 statistics
    #> > End statistics:
    #>   >> 1: baseline.end
    #>   >> 2: totaldegreeDyad.end
    #>   >> 3: difference_age.end

The returned object is of class `"remstats_durem"`. The statistics element can be viewed as a stacked remstats object:

    head(dstats$stacked$remstats_stack)
    #>   obs log_interevent baseline.start activeOutdegreeSender.start inertia.start
    #> 1   1       1.791759              0                   0.0000000             0
    #> 2   0       1.791759              1                  -0.4873397             0
    #> 3   0       1.791759              1                  -0.4873397             0
    #> 4   0       1.791759              1                  -0.4873397             0
    #> 5   0       1.791759              1                  -0.4873397             0
    #> 6   0       1.791759              1                  -0.4873397             0
    #>   reciprocity.start baseline.end totaldegreeDyad.end difference_age.end
    #> 1                 0            1                   0         -0.3721042
    #> 2                 0            0                   0          0.0000000
    #> 3                 0            0                   0          0.0000000
    #> 4                 0            0                   0          0.0000000
    #> 5                 0            0                   0          0.0000000
    #> 6                 0            0                   0          0.0000000
    #>   time_index dyad process actor1 actor2
    #> 1          2    4     end      1      5
    #> 2          2    1   start      1      2
    #> 3          2    2   start      1      3
    #> 4          2    3   start      1      4
    #> 5          2    4   start      1      5
    #> 6          2    5   start      2      1

## 3.3 Event types and `consider_type`

Similar as history-based statistics, the active-state statistics accept a `consider_type` argument that controls how the `setting` event type is handled when the model is built on a typed risk set:

Value | Behaviour
--- | ---
`"ignore"` (default) | Counts active events of any type
`"separate"` | One statistic per type; only type-\(c\) events contribute to the type-\(c\) stat
`"interact"` | One statistic for every combination of past event types and the event type in the riskset (only possible when setting in of the `remify` object). In case of `C` event types, this results in `C^2` distinct statistics.

    # Rename `setting` to `type` and build a typed duration history
    history_typed <- randomREH3
    names(history_typed)[names(history_typed) == "setting"] <- "type"
    reh_typed <- remify(history_typed, duration = TRUE, model = "tie", directed = TRUE)

    # "separate": one activeTie statistic per setting
    dstats_sep <- remstats(
        reh_typed,
        start_effects = ~ activeTie(consider_type = "separate")
    )
    dimnames(dstats_sep$start_stats)[[3]]

# 4 Step 3 — Estimation with `remstimate()`

`remstimate()` dispatches on the `"remstats_durem"` class and fits the joint start/end model. The start and end sub-models are stacked into a single long-format design matrix; start statistics are zero for end rows and vice versa.

## 4.1 Fitting the model

    fit <- remstimate(reh = reh, stats = dstats)

    fit
    #> Relational Event Model (tie oriented, duration)
    #> Estimation: MLE [ glm ]
    #>
    #> Coefficients:
    #>
    #>              baseline.start activeOutdegreeSender.start
    #>                  -7.2005633                   0.0461310
    #>               inertia.start           reciprocity.start
    #>                   0.3852725                   0.4169325
    #>                baseline.end         totaldegreeDyad.end
    #>                  -2.9004091                   0.0501352
    #>          difference_age.end
    #>                   0.3787515
    #>
    #> Null deviance:     28873.87
    #> Residual deviance: 23182.67
    #> AIC: 23196.67  AICC: 23196.73  BIC: 23235.73

## 4.2 Coefficients and summary

    coef(fit)
    #>              baseline.start activeOutdegreeSender.start
    #>                  -7.2005633                   0.0461310
    #>               inertia.start           reciprocity.start
    #>                   0.3852725                   0.4169325
    #>                baseline.end         totaldegreeDyad.end
    #>                  -2.9004091                   0.0501352
    #>          difference_age.end
    #>                   0.3787515

    summary(fit)
    #> Relational Event Model (tie oriented, duration)
    #> Estimation: MLE [ glm ]
    #> --------------------------------------------------
    #>
    #> Coefficients:
    #>                              Estimate  Std. Err   z value Pr(>|z|)
    #> baseline.start              -7.200563  0.039916 -180.3928  < 2e-16 ***
    #> activeOutdegreeSender.start  0.046131  0.064670    0.7133  0.47564
    #> inertia.start                0.385273  0.036346   10.6001  < 2e-16 ***
    #> reciprocity.start            0.416932  0.036146   11.5348  < 2e-16 ***
    #> baseline.end                -2.900409  0.035847  -80.9113  < 2e-16 ***
    #> totaldegreeDyad.end          0.050135  0.030040    1.6689  0.09513 .
    #> difference_age.end           0.378751  0.034933   10.8422  < 2e-16 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Null deviance:     28873.87  on 1957 degrees of freedom
    #> Residual deviance: 23182.67  on 1950 degrees of freedom
    #> Chi-square: 5691.2 on 7 degrees of freedom, p-value 0
    #> AIC: 23196.67  AICC: 23196.73  BIC: 23235.73

Each coefficient quantifies the effect of the corresponding statistic on the start or end event rate. For example, `activeOutdegreeSender.start` captures engagement: a sender who is already in many active events may be less (or more) likely to initiate another, corresponding to a negative (or positive) coefficient, while `totaldegreeDyad.end` describes how a dyad’s accumulated activity in past events relates to how quickly its events terminate.

## 4.3 Diagnostics

    diag <- diagnostics(fit, reh, dstats)
    plot(diag)

## 4.4 Model comparison

Because `remstimate_durem` returns a standard `logLik` object, information criteria compare nested models in the usual way:

    # A simpler model: baseline dynamics only
    dstats_simple <- remstats(
        reh,
        start_effects = ~ inertia(scaling = "std"),
        end_effects   = ~ totaldegreeDyad(scaling = "std"),
        memory        = "decay",
        memory_value  = 2000
    )
    fit_simple <- remstimate(reh = reh, stats = dstats_simple)

    logLik(fit)
    #> 'log Lik.' -11591.34 (df=7)
    logLik(fit_simple)
    #> 'log Lik.' -11708.23 (df=4)

    AIC(fit)
    #> [1] 23196.67
    AIC(fit_simple)
    #> [1] 23424.45

# 5 Additional topics

## 5.1 Duration based event weights

When computing the history-weighted (endogenous) effects, `inertia()`, `reciprocity()`, etc., the duration of past events are used as event weights according to the formulas:

- `(duration + 1)^psi_start` for the `start_effects`
- `(duration + 1)^psi_end` for the `end_effects`.

Thus, the `psi_start` and `psi_end` arguments of `remstats` control the duration-weighting exponent applied to the history-weighted terms (and are ignored for active-state terms). Defaults are `psi_start = 1` and `psi_end = 1`, implying that events with longer duration result in a larger weight in the endogenous statistics. On the other hand, if these arguments would be set to a negative value, then events with a longer duration have a smaller weight in the endogenous statistics. Thus, the `psi_start` and `psi_end` have a comparable role as the `memory_value` argument which affects the weight of events as function of the transpired time. Moreover, if events also have a weight by its (e.g., for quantifying the intensity of an event), specified via `event_weight`, the weight of the past event is computed by multiplying its `weight` with the duration-based weight.

## 5.2 Right-censored events

Right-censored events (`end = NA` in the edgelist) are fully supported. They contribute a start row to the dual edgelist and appear in the active-state counts until the end of the observation window, but never produce an observed end event. The `$durem` slot records how many events are censored:

    reh$durem$n_censored
    #> [1] 0
    reh$durem$n_complete
    #> [1] 999

## 5.3 Simultaneous events and boundary coincidences

The dual edgelist may contain multiple rows at the same time point — two events starting simultaneously, or one event ending exactly when another begins. `remify_durem` handles these automatically, computing the active state consistently just before each time point.

    # Time points in the dual edgelist shared by more than one row
    tab <- table(reh$edgelist_dual$time)
    head(tab[tab > 1], 6)
    #>
    #>  523  945 1240 1954 4329 4944
    #>    2    2    2    2    2    2

## 5.3 Directed end process

By default the end process is undirected — either actor can terminate an event. To model a directed end process (e.g. when a specific actor ends the event), set `dur_directed_end = TRUE`. This makes the directed active-state effects available in `end_effects` as well:

    reh_de <- remify(randomREH3, duration = TRUE, model = "tie",
                     directed = TRUE, dur_directed_end = TRUE)

    dstats_de <- remstats(
        reh_de,
        start_effects = ~ inertia(scaling = "std") + activeOutdegreeSender(scaling = "std"),
        end_effects   = ~ activeOutdegreeSender(scaling = "std")
    )
    fit_de <- remstimate(reh = reh_de, stats = dstats_de)
    coef(fit_de)

#### References

Butts, C. 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

Hoffman, M., Block, P., Elmer, T., & Stadtfeld, C. (2020). A model for the dynamics of face-to-face interactions in social groups. Network Science, 8(S1), S4–S25. https://doi.org/10.1017/nws.2020.3

Lakdawala, R., Leenders, R., Ejbye-Ernst, P., & Mulder, J. (2026). Modelling interaction duration in relational event models. arXiv preprint. https://doi.org/10.48550/arXiv.2602.21000