# Frailty Relational Event Models with GLMM

### Random intercepts and slopes for actor heterogeneity

#### remverse Team

Standard relational event models assume that all actors and dyads share the same set of parameters. In practice, actors differ in their baseline propensity to interact — some actors are inherently more active than others, beyond what is captured by endogenous statistics. A frailty REM addresses this by introducing random intercepts (and optionally random slopes) for actors, borrowing ideas from survival analysis and generalized linear mixed models (GLMMs). Additionally, random effects can be introduced for the other statistics as well across actors and dyads.

In `remverse`, a GLMM is requested simply by passing a `random` formula to `remstimate()`; the convenience wrapper `remfrailty()` builds the standard sender/receiver frailty structure for you. Under the hood the stacked relational event model is fit as a Poisson GLMM (via `glmmTMB` or `lme4` for interval timing, or `coxme` for ordinal timing).

    library(remverse)

# 1 When to use a frailty model

A frailty model is appropriate when:

- Actor heterogeneity is expected but not fully captured by observed covariates. Some actors may simply be more sociable or more central, independent of their network position.
- Overdispersion is present in the event rates — the standard Poisson model underpredicts the variance of event counts.
- Shrinkage is desirable — actors with few observations borrow strength from the population-level estimates rather than producing noisy individual-level estimates.

The random intercept acts as a multiplicative frailty on the event rate: actors with a positive random intercept have a higher-than-average baseline rate, while those with a negative intercept have a lower-than-average rate.

Additionally, random effects can be specified for the other statistics as well via the `random` argument.

# 2 Data

We use the `randomREH3` event history and the `info3` actor attributes provided with `remverse`. Statistics use exponential memory decay (half-life 2000). To keep the GLMM fits fast in this vignette, we compute statistics on the last stretch of the sequence via the `first` argument (a GLMM stacks the full case-control design, which is heavier than a plain MLE fit).

    data(randomREH3)
    data(info3)

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

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

# 3 The standard model (no random effects)

As a baseline, we first estimate a standard fixed-effects-only model.

    fit_fixed <- remstimate(reh = reh, stats = stats)
    summary(fit_fixed)
    #> Relational Event Model (tie oriented)
    #>
    #> Call:
    #> ~inertia(scaling = "std") + reciprocity(scaling = "std")
    #>
    #>
    #> Coefficients (MLE with interval likelihood):
    #>
    #>               Estimate   Std. Err    z value Pr(>|z|)    Pr(=0)
    #> baseline     -7.558156   0.077014 -98.139386        0 < 2.2e-16
    #> inertia       0.418518   0.071551   5.849216        0 6.380e-07
    #> reciprocity   0.401295   0.071797   5.589292        0 2.821e-06
    #> Null deviance: 4833.688 on 294 degrees of freedom
    #> Residual deviance: 4559.24 on 291 degrees of freedom
    #> Chi-square: 274.4476 on 3 degrees of freedom, asymptotic p-value 0
    #> AIC: 4565.24 AICC: 4565.323 BIC: 4576.291

This model assumes all actors share the same baseline rate and the same coefficients. Any unobserved actor-level heterogeneity is absorbed into the residual.

# 4 Random intercepts for actors

## 4.1 Sender frailty

A random intercept for `actor1` (the sender) allows each actor to have a different baseline sending rate:

\[\log \lambda_{ij}(t) = \boldsymbol{\beta}^\top \mathbf{x}_{ij}(t) + u_i, \quad u_i \sim \mathcal{N}(0, \sigma^2_u)\]

where \(u_i\) is the sender-specific random intercept. Random effects are requested via the `random` argument, which uses `glmmTMB` (default) or `lme4`.

    fit_sender <- remstimate(
      reh    = reh,
      stats  = stats,
      random = ~ (1 | actor1)
    )
    summary(fit_sender)
    #>  Family: poisson  ( log )
    #> Formula:          obs ~ -1 + baseline + inertia + reciprocity + (1 | actor1) +
    #>     offset(log_interevent)
    #> Data: df
    #>
    #>       AIC       BIC    logLik -2*log(L)  df.resid
    #>    2384.9    2411.6   -1188.4    2376.9      5876
    #>
    #> Random effects:
    #>
    #> Conditional model:
    #>  Groups Name        Variance Std.Dev.
    #>  actor1 (Intercept) 0.002387 0.04886
    #> Number of obs: 5880, groups:  actor1, 5
    #>
    #> Conditional model:
    #>             Estimate Std. Error z value Pr(>|z|)
    #> baseline    -7.56005    0.08122  -93.08  < 2e-16 ***
    #> inertia      0.41611    0.07377    5.64 1.69e-08 ***
    #> reciprocity  0.40835    0.08597    4.75 2.04e-06 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The output reports fixed-effect coefficients (comparable to the standard model) and the estimated variance of the sender random intercept (\(\hat{\sigma}^2_u\)). A large variance indicates substantial sender-level heterogeneity not captured by the fixed effects. (With only five actors here, the variance component is estimated from only five actors (‘groups’), which is generally too low for accurate estimation of the heterogeneity). Based on the above output, the standard deviation of senders’ activity is (only) `0.04886`.

## 4.2 Receiver frailty

Similarly, a random intercept for `actor2` captures receiver-level heterogeneity — some actors are more attractive as receivers regardless of the network history:

    fit_receiver <- remstimate(
      reh    = reh,
      stats  = stats,
      random = ~ (1 | actor2)
    )
    summary(fit_receiver)
    #>  Family: poisson  ( log )
    #> Formula:          obs ~ -1 + baseline + inertia + reciprocity + (1 | actor2) +
    #>     offset(log_interevent)
    #> Data: df
    #>
    #>       AIC       BIC    logLik -2*log(L)  df.resid
    #>    2379.5    2406.2   -1185.8    2371.5      5876
    #>
    #> Random effects:
    #>
    #> Conditional model:
    #>  Groups Name        Variance Std.Dev.
    #>  actor2 (Intercept) 0.04376  0.2092
    #> Number of obs: 5880, groups:  actor2, 5
    #>
    #> Conditional model:
    #>             Estimate Std. Error z value Pr(>|z|)
    #> baseline    -7.59807    0.12368  -61.43  < 2e-16 ***
    #> inertia      0.41981    0.07734    5.43 5.70e-08 ***
    #> reciprocity  0.46264    0.07842    5.90 3.65e-09 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Based on the above output, the standard deviation of receivers’ popularity is `0.2092`, which is larger than the heterogeneity of the senders’ activity.

## 4.3 Crossed sender and receiver frailties

Both can be included simultaneously as crossed random effects. This is the standard frailty-REM specification, and `remfrailty()` builds it for you: for a directed tie model it uses `(1 | actor1) + (1 | actor2)`.

    fit_crossed <- remfrailty(reh = reh, stats = stats)
    summary(fit_crossed)
    #>  Family: poisson  ( log )
    #> Formula:          obs ~ -1 + baseline + inertia + reciprocity + (1 | actor1) +
    #>     (1 | actor2) + offset(log_interevent)
    #> Data: df
    #>
    #>       AIC       BIC    logLik -2*log(L)  df.resid
    #>    2381.5    2414.9   -1185.8    2371.5      5875
    #>
    #> Random effects:
    #>
    #> Conditional model:
    #>  Groups Name        Variance  Std.Dev.
    #>  actor1 (Intercept) 5.921e-09 7.695e-05
    #>  actor2 (Intercept) 4.376e-02 2.092e-01
    #> Number of obs: 5880, groups:  actor1, 5; actor2, 5
    #>
    #> Conditional model:
    #>             Estimate Std. Error z value Pr(>|z|)
    #> baseline    -7.59808    0.12368  -61.43  < 2e-16 ***
    #> inertia      0.41981    0.07734    5.43 5.69e-08 ***
    #> reciprocity  0.46264    0.07842    5.90 3.65e-09 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The model estimates separate variance components for sender and receiver heterogeneity, and the fixed-effect coefficients are adjusted for both sources of unobserved variation. (An explicit `remstimate(..., random = ~ (1 | actor1) + (1 | actor2))` gives the same fit.) Again we see a considerable larger heterogeneity in the receivers’ (`actor1`) popularity than in the heterogeneity of the senders’ (`actor2`) activity.

# 5 Random slopes

Beyond intercepts, we can allow the effect of an endogenous statistic to vary across actors. For example, some actors may be more habitual than others — their inertia effect is stronger — modeled with a random slope for inertia on the sender:

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

    fit_slope <- remstimate(
      reh    = reh,
      stats  = stats_inertia,
      random = ~ (1 + inertia | actor1)
    )
    summary(fit_slope)
    #>
    #> ** fit did NOT converge - estimates unreliable **
    #>  Family: poisson  ( log )
    #> Formula:          obs ~ -1 + baseline + inertia + reciprocity + (1 + inertia |
    #>     actor1) + offset(log_interevent)
    #> Data: df
    #>
    #>       AIC       BIC    logLik -2*log(L)  df.resid
    #>    2387.1    2427.2   -1187.6    2375.1      5874
    #>
    #> Random effects:
    #>
    #> Conditional model:
    #>  Groups Name        Variance  Std.Dev. Corr
    #>  actor1 (Intercept) 0.0005941 0.02437
    #>         inertia     0.0119960 0.10953  1.00
    #> Number of obs: 5880, groups:  actor1, 5
    #>
    #> Conditional model:
    #>             Estimate Std. Error z value Pr(>|z|)
    #> baseline    -7.56611    0.07902  -95.75  < 2e-16 ***
    #> inertia      0.46292    0.09943    4.66 3.23e-06 ***
    #> reciprocity  0.42552    0.08008    5.31 1.07e-07 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The random effects now include a variance for the intercept, a variance for the inertia slope, and their correlation. A positive correlation would mean that actors with a high baseline sending rate also have a stronger inertia effect. Random slopes substantially increase model complexity; with small datasets the model may struggle to converge, which the `summary()` output would indicate.

# 6 Model comparison

We can compare models using information criteria:

    cat("Fixed-effects only: AIC =", AIC(fit_fixed),    "\n")
    #> Fixed-effects only: AIC = 4565.24
    cat("Sender frailty:     AIC =", AIC(fit_sender$backend_fit),   "\n")
    #> Sender frailty:     AIC = 2384.862
    cat("Receiver frailty:   AIC =", AIC(fit_receiver$backend_fit), "\n")
    #> Receiver frailty:   AIC = 2379.504
    cat("Crossed frailties:  AIC =", AIC(fit_crossed$backend_fit),  "\n")
    #> Crossed frailties:  AIC = 2381.504

A lower AIC indicates a better trade-off between fit and complexity. If a frailty model improves AIC substantially, actor-level heterogeneity is an important feature of the data.

Alternatively, the BIC, which has a stronger penalty for more complex models, thus, preferring simpler models, could also be used:

    cat("Fixed-effects only: BIC =", BIC(fit_fixed),    "\n")
    #> Fixed-effects only: BIC = 4576.291
    cat("Sender frailty:     BIC =", BIC(fit_sender$backend_fit),   "\n")
    #> Sender frailty:     BIC = 2411.579
    cat("Receiver frailty:   BIC =", BIC(fit_receiver$backend_fit), "\n")
    #> Receiver frailty:   BIC = 2406.221
    cat("Crossed frailties:  BIC =", BIC(fit_crossed$backend_fit),  "\n")
    #> Crossed frailties:  BIC = 2414.9

## Both models have a preference towards the `Receiver frailty` model.

# 7 Diagnostics

The GLMM diagnostics compute recall using the full linear predictor, including the best linear unbiased predictors (BLUPs) of the random effects. As throughout `remverse`, `diagnostics()` takes the fit together with the `reh` and `stats` objects it was built from.

    diag_crossed <- diagnostics(fit_crossed, reh, stats)
    diag_crossed
    #> Diagnostics for a Relational Event Model
    #> ------------------------------------------
    #> Model       : tie
    #> Actors      : 5
    #> Events      : 294
    #>
    #> Tie model
    #>   Statistics : inertia, reciprocity
    #>   Recall     : mean rank = 0.793  |  prob ratio = 2.34  |  top 5% = 17.6%  |  lowest 20% = 2.4%

    plot(diag_crossed)

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

For comparison, the recall can also be computed using fixed effects only (without the random-effect predictions), by setting `use_ranef = FALSE`:

    diag_fixed_only <- diagnostics(fit_crossed, reh, stats, use_ranef = FALSE)
    diag_fixed_only
    #> Diagnostics for a Relational Event Model
    #> ------------------------------------------
    #> Model       : tie
    #> Actors      : 5
    #> Events      : 294
    #>
    #> Tie model
    #>   Statistics : inertia, reciprocity
    #>   Recall     : mean rank = 0.793  |  prob ratio = 2.34  |  top 5% = 17.6%  |  lowest 20% = 2.4%

The difference between the two recall values quantifies how much of the model’s predictive power comes from the random effects vs. the fixed effects.

# 8 Interpreting random effects

The individual random effects (BLUPs) can be extracted from the fitted model to identify actors with unusually high or low activity:

    re <- lme4::ranef(fit_crossed$backend_fit)
    re
    #> $actor1
    #>     (Intercept)
    #> 1  4.891003e-09
    #> 2  5.867691e-08
    #> 3  2.502861e-10
    #> 4 -2.865893e-08
    #> 5 -3.216297e-08
    #>
    #> $actor2
    #>    (Intercept)
    #> 1  0.005184616
    #> 2  0.112582916
    #> 3 -0.245931796
    #> 4 -0.114052301
    #> 5  0.264361500

Actors with large positive sender random effects are more active senders than average; actors with large positive receiver effects are more popular as receivers.

# 9 Summary

Specification | Call | Interpretation
--- | --- | ---
Sender frailty | `remstimate(random = ~ (1 \| actor1))` | Sender-level heterogeneity in baseline rate
Receiver frailty | `remstimate(random = ~ (1 \| actor2))` | Receiver-level heterogeneity (attractiveness)
Crossed frailties | `remfrailty(reh, stats)` | Both sender and receiver heterogeneity
Sender random slope | `remstimate(random = ~ (1 + inertia \| actor1))` | Actor-varying inertia effect

The GLMM approach integrates naturally into the `remverse` pipeline. The same `remify()` and `remstats()` calls are used; only the `remstimate()` step gains a `random` argument (or is replaced by `remfrailty()`). Diagnostics and plotting work identically.

#### References

Juozaitiene, R., & Wit, E. C. (2024). Nodal heterogeneity can induce ghost triadic effects in relational event models. Psychometrika, 89(1), 151-171. https://doi.org/10.1007/s11336-024-09952-x

Mulder, J., & Hoff, P. D. (2024). A latent variable approach for modeling relational data with multiple receivers. Annals of Applied Statistics. https://doi.org/10.1214/24-AOAS1885