Regularized Relational Event Models
Variable selection and shrinkage for REMs
When building relational event models, researchers often face a choice between including many network statistics to capture complex dynamics and keeping the model parsimonious to avoid overfitting. Regularized estimation provides a principled middle ground: fit all candidate statistics simultaneously and let the data determine which effects matter through penalization.
remverse supports two flavours of penalized estimation,
both requested through the penalty argument of
remstimate() or the convenience wrapper
rempenalty():
- Frequentist (
approach = "frequentist", the default): an elastic-net penalty viaglmnet, with cross-validated selection of the penalty strength. - Bayesian (
approach = "Bayesian"): shrinkage priors (default horseshoe) viashrinkem. This tends to produce parsimonious set of non-zero effects (via the posterior mode), and has been validated in REMs (Karimova et al., 2023, 2025).
1 Data
We use the randomREH3 event history and
info3 actor attributes provided with remverse.
Statistics use exponential memory decay (half-life 2000; see the
pipeline vignette on why full memory is often unrealistic), and to keep
cross-validation quick they are computed on the last stretch of the
sequence via first.
2 The kitchen-sink model
We start by specifying a rich model with many endogenous and exogenous effects — the kind of model a researcher might consider when exploring which network mechanisms drive the event process.
effects_rich <- ~ inertia(scaling = "std") +
reciprocity(scaling = "std") +
indegreeSender(scaling = "std") +
indegreeReceiver(scaling = "std") +
outdegreeSender(scaling = "std") +
outdegreeReceiver(scaling = "std") +
totaldegreeDyad(scaling = "std") +
isp(scaling = "std") +
osp(scaling = "std") +
itp(scaling = "std") +
otp(scaling = "std") +
send("age", attr_actors = info3, scaling = "std") +
receive("age", attr_actors = info3, scaling = "std") +
difference("age", attr_actors = info3, scaling = "std")
stats_rich <- remstats(reh = reh, tie_effects = effects_rich, first = 500,
memory = "decay", memory_value = 2000)
dimnames(stats_rich)[[3]]
#> [1] "baseline" "inertia" "reciprocity"
#> [4] "indegreeSender" "indegreeReceiver" "outdegreeSender"
#> [7] "outdegreeReceiver" "totaldegreeDyad" "isp"
#> [10] "osp" "itp" "otp"
#> [13] "send_age" "receive_age" "difference_age"That gives us 15 statistics (including the baseline) — several of them (the degree and shared-partner families) strongly correlated with one another. This is exactly the setting where regularization helps.
2.1 Unregularized MLE
fit_mle <- remstimate(reh = reh, stats = stats_rich)
summary(fit_mle)
#> Relational Event Model (tie oriented)
#>
#> Call:
#> ~inertia(scaling = "std") + reciprocity(scaling = "std") + indegreeSender(scaling = "std") + indegreeReceiver(scaling = "std") + outdegreeSender(scaling = "std") + outdegreeReceiver(scaling = "std") + totaldegreeDyad(scaling = "std") + isp(scaling = "std") + osp(scaling = "std") + itp(scaling = "std") + otp(scaling = "std") + send("age", attr_actors = info3, scaling = "std") + receive("age", attr_actors = info3, scaling = "std") + difference("age", attr_actors = info3, scaling = "std")
#>
#>
#> Coefficients (MLE with interval likelihood):
#>
#> Estimate Std. Err z value Pr(>|z|) Pr(=0)
#> baseline -7.681876 0.075971 -101.115568 0.0000 < 2e-16
#> inertia 0.366962 0.089738 4.089253 0.0000 0.00517
#> reciprocity 0.185452 0.088426 2.097252 0.0360 0.71137
#> indegreeSender 0.096183 0.167279 0.574987 0.5653 0.94960
#> indegreeReceiver 0.086965 0.168375 0.516498 0.6055 0.95110
#> outdegreeSender 0.012251 0.138604 0.088391 0.9296 0.95678
#> outdegreeReceiver 0.102047 0.137427 0.742550 0.4578 0.94404
#> totaldegreeDyad -0.198520 0.296214 -0.670191 0.5027 0.94668
#> isp -0.279742 0.171403 -1.632075 0.1027 0.85439
#> osp -0.051311 0.199216 -0.257566 0.7967 0.95556
#> itp -0.279825 0.196007 -1.427627 0.1534 0.88916
#> otp -0.144164 0.191927 -0.751140 0.4526 0.94370
#> send_age -0.077836 0.056903 -1.367887 0.1713 0.89713
#> receive_age -0.010558 0.056332 -0.187429 0.8513 0.95622
#> difference_age -0.166995 0.080896 -2.064303 0.0390 0.72524
#> Null deviance: 8069.974 on 494 degrees of freedom
#> Residual deviance: 7594.033 on 479 degrees of freedom
#> Chi-square: 475.9415 on 15 degrees of freedom, asymptotic p-value 0
#> AIC: 7624.033 AICC: 7625.037 BIC: 7687.071With many correlated parameters, some coefficients have large
standard errors and the model is prone to overfitting — good in-sample
fit, poor prediction on held-out events. Note here that the p-values and
posterior probabilities of zero effects, in columns
Pr(>|z|) and Pr(=0), respectively, already
indicate that there is not a lot of evidence for including many effects.
In this multiple testing problem, however, the question is what
statistically sound thresholds should be for including/excluding
effects. This is achieved via statistical regularization techniques.
3 Bayesian regularization
The Bayesian route applies a shrinkage prior (horseshoe by default) that pulls small effects strongly toward zero while leaving genuinely large effects essentially unpenalized. It is fast and typically yields a clean, parsimonious solution.
fit_bayes <- rempenalty(reh = reh, stats = stats_rich, approach = "Bayesian")
summary(fit_bayes)
#> REM - tie model - SHRINKEM [ horseshoe ] regularization
#>
#> input.est shrunk.mean shrunk.median shrunk.mode shrunk.lower
#> baseline -7.682 -7.637 -7.637 -7.638 -7.787
#> inertia 0.367 0.364 0.363 0.363 0.207
#> reciprocity 0.185 0.206 0.208 0.210 0.038
#> indegreeSender 0.096 0.002 0.000 0.000 -0.104
#> indegreeReceiver 0.087 0.006 0.002 0.000 -0.099
#> outdegreeSender 0.012 -0.038 -0.026 0.000 -0.171
#> outdegreeReceiver 0.102 0.013 0.006 0.000 -0.085
#> totaldegreeDyad -0.199 -0.019 -0.007 0.000 -0.188
#> isp -0.280 -0.306 -0.318 -0.358 -0.562
#> osp -0.051 -0.079 -0.057 0.000 -0.326
#> itp -0.280 -0.205 -0.198 -0.001 -0.515
#> otp -0.144 -0.091 -0.060 0.000 -0.369
#> send_age -0.078 -0.037 -0.029 0.000 -0.141
#> receive_age -0.011 -0.005 -0.002 0.000 -0.088
#> difference_age -0.167 -0.139 -0.139 -0.156 -0.306
#> shrunk.upper nonzero
#> baseline -7.488 TRUE
#> inertia 0.522 TRUE
#> reciprocity 0.362 TRUE
#> indegreeSender 0.115 FALSE
#> indegreeReceiver 0.120 FALSE
#> outdegreeSender 0.057 FALSE
#> outdegreeReceiver 0.132 FALSE
#> totaldegreeDyad 0.115 FALSE
#> isp -0.003 TRUE
#> osp 0.112 FALSE
#> itp 0.027 FALSE
#> otp 0.089 FALSE
#> send_age 0.037 FALSE
#> receive_age 0.072 FALSE
#> difference_age 0.006 FALSE
round(coef(fit_bayes), 3)
#> baseline inertia reciprocity indegreeSender
#> -7.638 0.363 0.210 0.000
#> indegreeReceiver outdegreeSender outdegreeReceiver totaldegreeDyad
#> 0.000 0.000 0.000 0.000
#> isp osp itp otp
#> -0.358 0.000 -0.001 0.000
#> send_age receive_age difference_age
#> 0.000 0.000 -0.156The effects that survive shrinkage identify the network mechanisms
the data actually support; the rest are shrunk to (near) zero. The
baseline is never penalized. We recommend to use the posterior mode as
point estimate shrunk.mode of the effects, resulting in 4
nonzero effects (besides the intercept). (Because the posterior mode is
estimated numerically from the posterior density, some effects may be
reported as marginally nonzero—for example,
itp = -0.001—but can be treated as zero.).
4 Frequentist regularization (elastic net)
The frequentist route uses glmnet. Cross-validation
selects the penalty strength \(\lambda\), and the elastic-net mixing
parameter alpha interpolates between lasso
(alpha = 1, the default) and ridge
(alpha = 0).
fit_glmnet <- rempenalty(reh = reh, stats = stats_rich, approach = "frequentist")
summary(fit_glmnet)
#>
#> Call: glmnet::cv.glmnet(x = Xfit, y = mm$y, weights = mm$weights, offset = mm$offset, nfolds = nfolds, family = if (!ordinal) "poisson" else "binomial", alpha = alpha, penalty.factor = pen_factor, intercept = !is.null(intercept_col))
#>
#> Measure: Poisson Deviance
#>
#> Lambda Index Measure SE Nonzero
#> min 0.000901 44 0.3049 0.008031 9
#> 1se 0.013375 15 0.3129 0.008566 3
coef(fit_glmnet)
#> baseline inertia reciprocity indegreeSender
#> -7.29691386 0.36851135 0.21433623 0.00000000
#> indegreeReceiver outdegreeSender outdegreeReceiver totaldegreeDyad
#> 0.00000000 0.00000000 0.00000000 0.00000000
#> isp osp itp otp
#> -0.06932612 0.00000000 0.00000000 0.00000000
#> send_age receive_age difference_age
#> 0.00000000 0.00000000 0.00000000Equivalently,
remstimate(reh, stats_rich, penalty = list(alpha = 1)) fits
the same model.
4.1 Comparing coefficients
coefs_mle <- coef(fit_mle)
coefs_bayes <- coef(fit_bayes)
coefs_glmnet <- coef(fit_glmnet)
shared <- intersect(intersect(names(coefs_mle), names(coefs_glmnet)), names(coefs_bayes))
comparison <- data.frame(
statistic = shared,
MLE = round(coefs_mle[shared], 3),
BAYES = round(coefs_bayes[shared], 3),
GLMNET = round(coefs_glmnet[shared], 3),
row.names = NULL
)
comparison
#> statistic MLE BAYES GLMNET
#> 1 baseline -7.682 -7.638 -7.297
#> 2 inertia 0.367 0.363 0.369
#> 3 reciprocity 0.185 0.210 0.214
#> 4 indegreeSender 0.096 0.000 0.000
#> 5 indegreeReceiver 0.087 0.000 0.000
#> 6 outdegreeSender 0.012 0.000 0.000
#> 7 outdegreeReceiver 0.102 0.000 0.000
#> 8 totaldegreeDyad -0.199 0.000 0.000
#> 9 isp -0.280 -0.358 -0.069
#> 10 osp -0.051 0.000 0.000
#> 11 itp -0.280 -0.001 0.000
#> 12 otp -0.144 0.000 0.000
#> 13 send_age -0.078 0.000 0.000
#> 14 receive_age -0.011 0.000 0.000
#> 15 difference_age -0.167 -0.156 0.000The selected column shows which effects survived
regularization. Effects shrunk to zero are considered non-informative
given the other effects in the model.
4.2 Alpha: lasso vs. ridge
# Pure lasso (default): maximum sparsity
fit_lasso <- rempenalty(reh, stats_rich, approach = "frequentist", alpha = 1)
# Pure ridge: shrinkage without variable selection
fit_ridge <- rempenalty(reh, stats_rich, approach = "frequentist", alpha = 0)
cat("Lasso non-zero:", sum(coef(fit_lasso) != 0), "of", length(coef(fit_lasso)), "\n")
#> Lasso non-zero: 3 of 15
cat("Ridge non-zero:", sum(coef(fit_ridge) != 0), "of", length(coef(fit_ridge)), "\n")
#> Ridge non-zero: 15 of 15The lasso selects a small subset of effects; ridge retains all but
shrinks them. Elastic net (0 < alpha < 1) sits in
between.
4.3 Lambda selection
By default rempenalty() uses lambda.1se
(the most parsimonious model within one standard error of the
cross-validation minimum). Use lambda_select = "min" for
the deviance-minimizing \(\lambda\):
fit_min <- rempenalty(reh, stats_rich, approach = "frequentist",
lambda_select = "min")
cat("lambda.1se non-zero:", sum(coef(fit_glmnet) != 0), "\n")
#> lambda.1se non-zero: 4
cat("lambda.min non-zero:", sum(coef(fit_min) != 0), "\n")
#> lambda.min non-zero: 105 Diagnostics
Recall diagnostics for a penalized model use its point estimates to
compute linear predictors. As elsewhere in remverse,
diagnostics() takes the fit plus the reh and
stats it was built from.
diag_mle <- diagnostics(fit_mle, reh, stats_rich)
diag_bayes <- diagnostics(fit_bayes, reh, stats_rich)
diag_glmnet <- diagnostics(fit_glmnet, reh, stats_rich)
cat("MLE recall: ", round(diag_mle$recall$summary$mean_rel_rank, 3), "\n")
#> MLE recall: 0.777
cat("BAYES recall: ", round(diag_bayes$recall$summary$mean_rel_rank, 3), "\n")
#> BAYES recall: 0.777
cat("GLMNET recall: ", round(diag_glmnet$recall$summary$mean_rel_rank, 3), "\n")
#> GLMNET recall: 0.776
plot(diag_mle)
#> Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
#> collapsing to unique 'x' values
#> Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
#> collapsing to unique 'x' values
#> Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
#> collapsing to unique 'x' values
In-sample recall may be slightly lower for the regularized model — it
trades a little in-sample fit for better generalization. In this case,
the recall is equal approximately across models while the regularized
solutions coefs_bayes and coefs_glmnet and
much more parsimonious than the unregularized solution
coefs_mle:
comparison
#> statistic MLE BAYES GLMNET
#> 1 baseline -7.682 -7.638 -7.297
#> 2 inertia 0.367 0.363 0.369
#> 3 reciprocity 0.185 0.210 0.214
#> 4 indegreeSender 0.096 0.000 0.000
#> 5 indegreeReceiver 0.087 0.000 0.000
#> 6 outdegreeSender 0.012 0.000 0.000
#> 7 outdegreeReceiver 0.102 0.000 0.000
#> 8 totaldegreeDyad -0.199 0.000 0.000
#> 9 isp -0.280 -0.358 -0.069
#> 10 osp -0.051 0.000 0.000
#> 11 itp -0.280 -0.001 0.000
#> 12 otp -0.144 0.000 0.000
#> 13 send_age -0.078 0.000 0.000
#> 14 receive_age -0.011 0.000 0.000
#> 15 difference_age -0.167 -0.156 0.0006 Summary
| Aspect | MLE | Penalized (frequentist / Bayesian) |
|---|---|---|
| Estimates | Unpenalized | Shrunk toward zero |
| Variable selection | No | Yes |
| Multicollinearity | Problematic | Handled |
| Overfitting risk | Higher with many effects | Lower |
| Uncertainty | Standard errors | glmnet: point estimates; shrinkem: posterior |
Penalized estimation integrates seamlessly into the
remverse pipeline: the same remify() and
remstats() calls are used, and only the
remstimate() step changes (via penalty = ...
or rempenalty()). A natural workflow is to use penalization
to screen a large set of candidate effects, then re-estimate the
selected model with MLE or a GLMM for full inference.
References
Karimova, D., Leenders, R., Meijerink-Bosman, M., & Mulder, J. (2023). Separating the wheat from the chaff: Bayesian regularization in dynamic social networks. Social Networks, 74, 139-155. https://doi.org/10.1016/j.socnet.2023.02.006
Karimova, D., van Erp, S., Leenders, R., & Mulder, J. (2025). Honey, I shrunk the irrelevant effects! Simple and flexible approximate Bayesian regularization. Journal of Mathematical Psychology, 126, 102925. https://doi.org/10.1016/j.jmp.2025.102925
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Methodological), 58(1), 267-288. https://doi.org/10.1111/j.2517-6161.1996.tb02080.x