---
title: "Duration Relational Event Models with remverse"
subtitle: "<i>A tutorial</i>"
author: ""
date: ""
output:
  rmarkdown::html_document:
    toc: true
    toc_depth: 3
header-includes:
  - \usepackage{amsmath,amssymb}
  - \usepackage{bm}
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Duration Relational Event Models with remverse}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(dpi = 72, collapse = TRUE, comment = "#>")
```

---

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`) |

```{r load-data}
library(remverse)   # loads remify, remstats, remstimate; remdata is a dependency

data(randomREH3)
data(info3)         # actor attributes (age, sex)
head(randomREH3)
```

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.

```{r data-summary}
cat(sprintf(
    "%d events: %d complete, %d right-censored\n",
    nrow(randomREH3),
    sum(!is.na(randomREH3$end)),
    sum( is.na(randomREH3$end))
))
```

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")`.

```{r remify}
reh <- remify(
    edgelist  = randomREH3,
    duration  = TRUE,
    model     = "tie",
    directed  = TRUE        # actor1 initiates; the end process is undirected by default
)
```

## 2.1  The `remify_durem` object

```{r print-reh}
reh
```

```{r summary-reh}
summary(reh)
```

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.

```{r plot-reh, out.width="70%", dev=c("jpeg"), dev.args = list(bg = "white")}
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:

```{r durem-slot}
reh$durem
```

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

```{r edgelist-dual}
head(reh$edgelist_dual, 12)
```

---

# 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 \code{?tie_effects}.
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 \code{?active_effects}.

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

```{r remstats, message=FALSE}
# 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
)
```

```{r print-dstats}
dstats
```

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

```{r dstats-dims}
head(dstats$stacked$remstats_stack)
```

## 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 \code{extend_riskset_by_type = TRUE} in of the `remify` object). In case of `C` event types, this results in `C^2` distinct statistics.|

```{r consider-type-demo, eval=FALSE}
# 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

```{r remstimate}
fit <- remstimate(reh = reh, stats = dstats)
```

```{r print-fit}
fit
```

## 4.2  Coefficients and summary

```{r coef}
coef(fit)
```

```{r summary-fit}
summary(fit)
```

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

```{r durem-diag, out.width="60%", dev=c("jpeg"), dev.args = list(bg = "white")}
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:

```{r model-comparison, message=FALSE}
# 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)
logLik(fit_simple)

AIC(fit)
AIC(fit_simple)
```

---

# 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:

```{r censored}
reh$durem$n_censored
reh$durem$n_complete
```

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

```{r simultaneous}
# Time points in the dual edgelist shared by more than one row
tab <- table(reh$edgelist_dual$time)
head(tab[tab > 1], 6)
```

## 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:

```{r directed-end, eval=FALSE}
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