---
title: "Dyadic Latent Class Relational Event Models"
subtitle: "<i>Finite mixture models for relational event heterogeneity</i>"
author: ""
date: ""
output:
  rmarkdown::html_document:
    toc: true
    toc_depth: 3
    theme: spacelab
    highlight: pygments
    code_folding: show
header-includes:
  - \usepackage{amsmath,amssymb}
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Dyadic Latent Class Relational Event Models}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(dpi = 72, collapse = TRUE, comment = "#>")
# Mixture fits need flexmix (declared in Suggests). If it is unavailable - e.g.
# a CRAN check run with _R_CHECK_FORCE_SUGGESTS_=false - the modelling chunks
# render code-only instead of erroring the vignette build.
has_flexmix <- requireNamespace("flexmix", quietly = TRUE)
```

---

<script>
window.MathJax = {
  "HTML-CSS": {
    availableFonts: [],
    preferredFont: null,
    webFont: "TeX"
  }
};
</script>


Standard relational event models assume that all dyads follow the same relational dynamics — a single set of parameters governs the entire event process. In many empirical networks this is too restrictive. Some dyads exhibit strong reciprocity while others are driven by inertia; some have a high baseline rate while others interact rarely. A **dyadic latent class REM** (Lakdawala et al., 2026) captures this by assuming that dyads belong to one of $K$ distinct latent classes, each with its own parameter vector.

The model is a finite mixture of relational event models:

$$\mathscr{L}(\boldsymbol{\Theta}) = \prod_{(i,j) \in \mathcal{D}} \sum_{k=1}^{K} \pi_k \prod_{m: e_m \in \mathcal{R}_{ij}} f(e_m \mid \boldsymbol{\beta}_k)$$

where $\pi_k$ is the mixing proportion of class $k$ and $\boldsymbol{\beta}_k$ its class-specific parameters. The key feature is that all events of the same dyad share a latent class — the *dyadic* grouping that distinguishes this model from event-level mixtures.

In `remverse` this is requested through the `mixture` argument of `remstimate()`, or the convenience wrappers `remixture()` (general mixture) and `dlcrem()` (the dyadic-latent-class special case with a `~ (1 | dyad)` grouping). Fitting uses `flexmix`.

```{r load}
library(remverse)
```

---

# 1 Data

We use the `randomREH3` event history provided with `remverse`. Statistics use exponential memory decay (half-life set to 2000). To train the endogenous statistics we start the computation from the 100th event via the `first` argument.

```{r data}
data(randomREH3)

reh <- remify(edgelist = randomREH3, model = "tie", directed = TRUE)

stats <- remstats(
  reh = reh,
  tie_effects = ~ inertia(scaling = "std") + reciprocity(scaling = "std"),
  first = 100, memory = "decay", memory_value = 2000
)
```

---

# 2 Fitting the model

`dlcrem()` fits the dyadic latent class model directly: it groups events by `dyad` and gives every included statistic a class-specific coefficient.

```{r fit-k2, message=FALSE, eval=has_flexmix}
fit_k2 <- dlcrem(reh = reh, stats = stats, k = 2)
fit_k2
```

Equivalently, `remixture(reh, stats, random = ~ (1 | dyad), k = 2)` or
`remstimate(reh, stats, mixture = list(k = 2, random = ~ (1 | dyad)))` fit the same model. All statistics vary across the $K$ classes; the `dyad` grouping ensures that events from the same dyad are assigned to the same class.

## 2.1 Interpreting the output

The output reports mixing proportions ($\pi_k$) and per-component coefficients ($\boldsymbol{\beta}_k$):

```{r coefs, eval=has_flexmix}
coef(fit_k2)
```

Each column is one latent class, sorted by mixing proportion (largest first). Differences between columns reveal *how* the classes differ — e.g. if reciprocity differs, some dyads are reciprocal while others are unidirectional; if the baseline differs, the classes capture high-rate vs. low-rate dyads.

---

# 3 Selecting the number of classes

The number of classes $K$ is usually unknown. Passing a vector to `k` fits all models at once and returns a list; `bic_table()` compares them by BIC.

```{r select-k, message=FALSE, eval=has_flexmix}
fit_multi <- remixture(
  reh    = reh,
  stats  = stats,
  random = ~ (1 | dyad),
  k      = 1:4,
  nrep   = 5
)

bic_table(fit_multi)
```

The model with the lowest BIC is preferred. If $K = 1$ wins, there is no evidence for dyadic heterogeneity in the included effects. The `nrep` argument sets the number of random initializations per $K$, guarding against local optima. Based on these results, the model with `k = 2` latent classes is preferred by the BIC.

```{r plot-k, out.width="60%", fig.align="center", dev=c("jpeg"), dev.args = list(bg = "white"), eval=has_flexmix}
plot(fit_multi)
```

## 3.1 Inspecting the selected model

```{r best-model, eval=has_flexmix}
bics   <- bic_table(fit_multi)
best_k <- bics$k[which.min(bics$BIC)]
cat("Best K:", best_k, "\n")

fit_best <- fit_multi[[paste0("k", best_k)]]
fit_best
```

---

# 4 Diagnostics

MIXREM diagnostics compute recall using posterior-weighted linear predictors: at each time point, a dyad's predicted rate is a weighted average across classes, weighted by its posterior class probabilities.

$$\hat{\lambda}_{ij}(t) = \sum_{k=1}^{K} P(Z_{ij} = k \mid \text{data}) \cdot \exp(\boldsymbol{\beta}_k^\top \mathbf{x}_{ij}(t))$$

As elsewhere in `remverse`, `diagnostics()` takes the fit together with the `reh` and `stats` it was built from.

```{r diagnostics, out.width="50%", dev=c("jpeg"), dev.args = list(bg = "white"), eval=has_flexmix}
diag_k2 <- diagnostics(fit_k2, reh, stats)
diag_k2

plot(diag_k2)
```

The diagnostics object also carries per-component recall (`recall_by_component`): how well each class's parameters predict that class's events. A component with much higher recall follows a clearer, more predictable pattern.

---

# 5 Comparison with MLE and GLMM

The three approaches address different sources of heterogeneity:

| Approach | Heterogeneity | Parameters | Interpretation |
|--------|--------------|------------|----------------|
| MLE | None assumed | Single $\boldsymbol{\beta}$ | One-size-fits-all |
| GLMM | Continuous (normal) | $\boldsymbol{\beta} + \mathbf{u}$ | Actor-level frailty |
| MIXREM | Discrete (classes) | $\boldsymbol{\beta}_1, \ldots, \boldsymbol{\beta}_K$ | Qualitatively different dyad types |

Use **GLMM** when actor effects vary continuously around a population mean; use **MIXREM** when you expect *distinct types* of relational dynamics (e.g. collaborative vs. hierarchical dyads).

```{r compare-loglik, message=FALSE, eval=has_flexmix}
fit_mle <- remstimate(reh = reh, stats = stats)

cat("MLE logLik:    ", round(as.numeric(logLik(fit_mle)), 2), "\n")
cat("MIXREM logLik: ", round(as.numeric(logLik(fit_k2)),  2), "\n")
```

---

# 6 Practical considerations

**Convergence and local optima.** Finite mixture models are sensitive to starting values. Use `nrep >= 3` (the default); if results vary across runs, increase it.

**Sample size.** Each class needs enough data to estimate its parameters. With $P$ parameters and $K$ classes the model has $K \times P + (K-1)$ free parameters (`P` parameters per class, and `K-1` mixing proportions). The effective sample size per class is the number of dyad-event pairs assigned to it, not the total number of events.

**Identifiability.** If two components converge to the same parameters, the model is over-specified — reduce $K$.

**Grouping variable.** `dyad` grouping (via `dlcrem()`) is natural for tie-oriented models. For a different grouping, use `remixture(reh, stats, random = ~ (1 | actor1), ...)`.

---

#### References

Lakdawala, R., Leenders, R., & Mulder, J. (2026). Not all bonds are created equal: Dyadic latent class models for relational event data. *Social Networks*. https://doi.org/10.1016/j.socnet.2026.06.006

Grün, B., & Leisch, F. (2008). FlexMix version 2: Finite mixtures with concomitant variables and varying and constant parameters. *Journal of Statistical Software*, 28(4), 1–35.