---
title: "Frailty Relational Event Models with GLMM"
subtitle: "<i>Random intercepts and slopes for actor heterogeneity</i>"
author: "remverse Team"
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{Frailty Relational Event Models with GLMM}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(dpi = 72, collapse = TRUE, comment = "#>")
# Frailty fits need lme4 (declared in Suggests). If it is unavailable - e.g. a
# CRAN check run with _R_CHECK_FORCE_SUGGESTS_=false - the modelling chunks
# render code-only instead of erroring the vignette build.
run_glmm <- requireNamespace("lme4", quietly = TRUE)
```

---

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

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

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

```{r standard}
fit_fixed <- remstimate(reh = reh, stats = stats)
summary(fit_fixed)
```

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

```{r sender-frailty, message=FALSE, warning=FALSE, eval=run_glmm}
fit_sender <- remstimate(
  reh    = reh,
  stats  = stats,
  random = ~ (1 | actor1)
)
summary(fit_sender)
```

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:

```{r receiver-frailty, message=FALSE, warning=FALSE, eval=run_glmm}
fit_receiver <- remstimate(
  reh    = reh,
  stats  = stats,
  random = ~ (1 | actor2)
)
summary(fit_receiver)
```

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

```{r crossed, message=FALSE, warning=FALSE, eval=run_glmm}
fit_crossed <- remfrailty(reh = reh, stats = stats)
summary(fit_crossed)
```

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:

```{r random-slope, message=FALSE, warning=FALSE, eval=run_glmm}
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)
```

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:

```{r comparison, eval=run_glmm}
cat("Fixed-effects only: AIC =", AIC(fit_fixed),    "\n")
cat("Sender frailty:     AIC =", AIC(fit_sender$backend_fit),   "\n")
cat("Receiver frailty:   AIC =", AIC(fit_receiver$backend_fit), "\n")
cat("Crossed frailties:  AIC =", AIC(fit_crossed$backend_fit),  "\n")
```

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:

```{r comparison-bic, eval=run_glmm}
cat("Fixed-effects only: BIC =", BIC(fit_fixed),    "\n")
cat("Sender frailty:     BIC =", BIC(fit_sender$backend_fit),   "\n")
cat("Receiver frailty:   BIC =", BIC(fit_receiver$backend_fit), "\n")
cat("Crossed frailties:  BIC =", BIC(fit_crossed$backend_fit),  "\n")
```
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.

```{r diagnostics, out.width="50%", dev=c("jpeg"), dev.args = list(bg = "white"), eval=run_glmm}
diag_crossed <- diagnostics(fit_crossed, reh, stats)
diag_crossed

plot(diag_crossed)
```

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

```{r diagnostics-fixed, eval=run_glmm}
diag_fixed_only <- diagnostics(fit_crossed, reh, stats, use_ranef = FALSE)
diag_fixed_only
```

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:

```{r ranef, eval=run_glmm}
re <- lme4::ranef(fit_crossed$backend_fit)
re
```

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