# Dyadic Latent Class Relational Event Models

### Finite mixture models for relational event heterogeneity

####

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

    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.

    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.

    fit_k2 <- dlcrem(reh = reh, stats = stats, k = 2)
    fit_k2
    #> REM - tie model - MIXREM [flexmix, k = 2 ]
    #>
    #> Mixing proportions:
    #> [1] 0.5557 0.4443
    #>
    #> Coefficients per component:
    #>             Component.1 Component.2
    #> baseline        -3.6645     -2.8896
    #> inertia          0.6362      0.3322
    #> reciprocity      0.3004      0.2148
    #>
    #> logLik: -3242.52  BIC: 6553.58

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\)):

    coef(fit_k2)
    #>             Component.1 Component.2
    #> baseline     -3.6644640  -2.8895626
    #> inertia       0.6361504   0.3322290
    #> reciprocity   0.3003923   0.2147843

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.

    fit_multi <- remixture(
      reh    = reh,
      stats  = stats,
      random = ~ (1 | dyad),
      k      = 1:4,
      nrep   = 5
    )

    bic_table(fit_multi)
    #>    k      BIC delta_BIC
    #> k2 2 6553.580   0.00000
    #> k1 1 6564.748  11.16761
    #> k3 3 6588.376  34.79576
    #> k4 4 6625.489  71.90851

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.

    plot(fit_multi)
    #> BIC per k:
    #>    k      BIC delta_BIC
    #> k2 2 6553.580   0.00000
    #> k1 1 6564.748  11.16761
    #> k3 3 6588.376  34.79576
    #> k4 4 6625.489  71.90851
    #> Use plot(x[["k2"]], ...) to plot a specific fit.

## 3.1 Inspecting the selected model

    bics   <- bic_table(fit_multi)
    best_k <- bics$k[which.min(bics$BIC)]
    cat("Best K:", best_k, "\n")
    #> Best K: 2

    fit_best <- fit_multi[[paste0("k", best_k)]]
    fit_best
    #> REM - tie model - MIXREM [flexmix, k = 2 ]
    #>
    #> Mixing proportions:
    #> [1] 0.5558 0.4442
    #>
    #> Coefficients per component:
    #>             Component.1 Component.2
    #> baseline        -3.6640     -2.8893
    #> inertia          0.6360      0.3321
    #> reciprocity      0.3005      0.2147
    #>
    #> logLik: -3242.52  BIC: 6553.58

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

    diag_k2 <- diagnostics(fit_k2, reh, stats)
    diag_k2
    #> Diagnostics - Mixture REM (k = 2 )
    #> Actors: 5  Events: 993
    #>
    #> Mixing proportions:
    #> Component.1 Component.2
    #>      0.5557      0.4443
    #>
    #> Recall (posterior-weighted):
    #>   Joint           mean rank = 0.756  |  prob ratio = 1.90  |  top 5% = 16.9%
    #>
    #> Recall by component:
    #>   Component.1     mean rank = 0.753  |  prob ratio = 2.02  |  top 5% = 16.8%
    #>   Component.2     mean rank = 0.754  |  prob ratio = 1.56  |  top 5% = 16.9%

    plot(diag_k2)
    #> Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
    #> collapsing to unique 'x' values

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

    fit_mle <- remstimate(reh = reh, stats = stats)

    cat("MLE logLik:    ", round(as.numeric(logLik(fit_mle)), 2), "\n")
    #> MLE logLik:     -6839.3
    cat("MIXREM logLik: ", round(as.numeric(logLik(fit_k2)),  2), "\n")
    #> MIXREM logLik:  -3242.52

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