---
title: "Regularized Relational Event Models"
subtitle: "<i>Variable selection and shrinkage for REMs</i>"
author: ""
date: ""
output:
  rmarkdown::html_document:
    toc: true
    toc_depth: 3
    theme: spacelab
    highlight: pygments
    code_folding: show
header-includes:
  - \usepackage{amsmath,amssymb}
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Regularized Relational Event Models}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(dpi = 72, collapse = TRUE, comment = "#>")
# Penalized fits need glmnet (frequentist) / shrinkem (Bayesian), both declared
# in Suggests. If a backend is unavailable - e.g. a CRAN check run with
# _R_CHECK_FORCE_SUGGESTS_=false - the relevant chunks render code-only instead
# of erroring the vignette build.
has_glmnet   <- requireNamespace("glmnet",   quietly = TRUE)
has_shrinkem <- requireNamespace("shrinkem", quietly = TRUE)
```

---

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).

```{r load}
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`.

```{r data}
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.

```{r kitchen-sink}
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]]
```

That gives us `r dim(stats_rich)[3]` 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

```{r mle-rich}
fit_mle <- remstimate(reh = reh, stats = stats_rich)
summary(fit_mle)
```

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.

```{r shrinkem, message=FALSE, eval=has_shrinkem}
fit_bayes <- rempenalty(reh = reh, stats = stats_rich, approach = "Bayesian")
summary(fit_bayes)
round(coef(fit_bayes), 3)
```

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

```{r glmnet-default, message=FALSE, eval=has_glmnet}
fit_glmnet <- rempenalty(reh = reh, stats = stats_rich, approach = "frequentist")
summary(fit_glmnet)
coef(fit_glmnet)
```

Equivalently, `remstimate(reh, stats_rich, penalty = list(alpha = 1))` fits the same model.

## 4.1 Comparing coefficients

```{r compare-coefs, eval=has_glmnet}
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
```

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

```{r alpha, message=FALSE, eval=has_glmnet}
# 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")
cat("Ridge non-zero:", sum(coef(fit_ridge) != 0), "of", length(coef(fit_ridge)), "\n")
```

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$:

```{r lambda, message=FALSE, eval=has_glmnet}
fit_min <- rempenalty(reh, stats_rich, approach = "frequentist",
                      lambda_select = "min")

cat("lambda.1se non-zero:", sum(coef(fit_glmnet) != 0), "\n")
cat("lambda.min non-zero:", sum(coef(fit_min)    != 0), "\n")
```

---

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

```{r diagnostics, out.width="50%", dev=c("jpeg"), dev.args = list(bg = "white"), eval=has_glmnet}
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")
cat("BAYES recall: ", round(diag_bayes$recall$summary$mean_rel_rank, 3), "\n")
cat("GLMNET recall: ", round(diag_glmnet$recall$summary$mean_rel_rank, 3), "\n")

plot(diag_mle)
plot(diag_bayes)
plot(diag_glmnet)
```

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`:

```{r show-comparison, eval=has_glmnet}
comparison
```

---

# 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