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.