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).
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.291This 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 ' ' 1The 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 ' ' 1Based 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 ' ' 1The 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 ' ' 1The 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.504A 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.9Both 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.264361500Actors 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