# 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 via `glmnet`, with cross-validated selection of the penalty strength.
- Bayesian (`approach = "Bayesian"`): shrinkage priors (default horseshoe) via `shrinkem`. 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).

    library(remverse)

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

    data(randomREH3)
    data(info3)

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

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

With 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.156

The 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.00000000

Equivalently, `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.000

The `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 15

The 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: 10

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

    plot(diag_bayes)

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

    plot(diag_glmnet)

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

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