---
title: "OD migration with Kronecker covariance"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{OD migration with Kronecker covariance}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  fig.width = 6,
  fig.height = 4,
  fig.align = "center"
)
set.seed(24)
```

## The OD setting

In an **origin-destination (OD)** model, each observation is a flow between
two regions. Migration counts between US states, commuting flows between
zip codes, or trade flows between countries all fit this pattern.

The random-effects design has two indicators per observation: one for the
*origin* region and one for the *destination* region. With $G$ regions,
the random-effects matrix $Z$ is $N \times 2G$, and every row sums to 2.

The natural covariance structure on the resulting $2G$-vector of random
effects is the **Kronecker** form
$$
\Sigma_\alpha \;=\; \underbrace{\Sigma_{2\times 2}}_{\text{origin/dest}}
                  \;\otimes\;
                  \underbrace{\Sigma_{\text{spatial}}}_{\text{between regions}},
$$
where $\Sigma_{2\times 2}$ captures the origin/destination variance plus
their covariance, and $\Sigma_{\text{spatial}}$ ($G \times G$) captures
spatial dependence between regions (often known a priori from geography).

This reduces $2G(2G+1)/2$ free covariance parameters to just **3** for
$\Sigma_{2\times 2}$ — a substantial saving when $G$ is moderate.

## Setup

```{r setup}
library(cevcmm)
```

## Load the bundled example

The package ships a simulated OD migration dataset of $N = 3000$
observations: $G = 10$ regions over 30 years, with all
$10 \times 10$ origin–destination pairs observed annually.

```{r load}
path <- system.file("extdata", "od_migration.csv", package = "cevcmm")
od   <- read.csv(path)
str(od)
head(od)
```

The columns:

- `origin`, `dest`: integer region IDs in $\{1, \dots, 10\}$
- `year`: 1 through 30
- `t`: year normalised to $[0, 1]$
- `wage_diff`: a continuous covariate whose effect on flows is allowed to
  vary smoothly with `t`
- `log_flow`: the response (log migration count)

## Build the random-effects design

```{r build-z}
G <- 10L
N <- nrow(od)

Z <- matrix(0, N, 2L * G)
Z[cbind(seq_len(N), od$origin)]    <- 1   # origin indicators in cols 1..G
Z[cbind(seq_len(N), G + od$dest)]  <- 1   # destination indicators in cols (G+1)..2G

# Verify the OD structure: every row sums to 2 (one origin + one dest)
table(rowSums(Z))
```

## Choose an initial spatial covariance

In a real analysis $\Sigma_{\text{spatial}}$ would come from geographic
distance. Here we use an exponential decay in region-index distance as
a reasonable proxy:

```{r sigma-spatial}
Sigma_spatial <- outer(seq_len(G), seq_len(G),
                       function(i, j) exp(-abs(i - j) / 3))
round(Sigma_spatial[1:4, 1:4], 3)
```

## Fit

A single `vcmm()` call with `re_cov = "kronecker"`. The package detects
the constant row-sum automatically and applies the identifiability shift
described in `vignette("distributed-fitting")`.

```{r fit}
fit <- vcmm(y = od$log_flow,
            X = od$wage_diff,
            Z = Z,
            t = od$t,
            method             = "csl",
            re_cov             = "kronecker",
            n_groups           = G,
            Sigma_spatial_init = Sigma_spatial,
            control            = vcmm_control(
              sigma_eps       = 0.6,
              sigma_alpha     = sqrt(0.5),
              update_variance = TRUE))
fit
```

## Interpret the estimated $\Sigma_{2 \times 2}$

```{r sigma-2x2}
Sigma_2x2_hat <- fit$re_cov_state$Sigma_left
round(Sigma_2x2_hat, 3)

# Correlation between origin and destination effects
corr_OD <- Sigma_2x2_hat[1, 2] /
           sqrt(Sigma_2x2_hat[1, 1] * Sigma_2x2_hat[2, 2])
round(corr_OD, 3)
```

The diagonal entries give the origin and destination variance components;
the off-diagonal entry summarises whether a region that tends to *send*
many migrants also tends to *receive* many. A positive correlation says
yes — high-traffic regions are high-traffic in both directions.

The true simulation values were $\Sigma_{2\times2} = \begin{pmatrix}
0.60 & 0.25 \\ 0.25 & 0.50 \end{pmatrix}$ (correlation 0.46).

**On small $G$, the estimated $\hat\Sigma_{2\times 2}$ will differ from
truth.** The 3 free parameters of $\Sigma_{2\times 2}$ are estimated from
effectively a single $G \times 2$ realisation of $M$, plus an EM
correction for the posterior uncertainty in $\alpha$ given the data. For
$G = 10$, sampling variability in the empirical $\text{cor}(M)$ has
standard error $\approx 1/\sqrt{G-2} \approx 0.35$, and the EM correction
inflates the diagonals further to account for posterior uncertainty. The
qualitative pattern — positive variance components, positive OD
correlation, "high-traffic" regions that send *and* receive heavily —
is the right takeaway at this $G$; the exact numerical values rely on a
larger network or repeated realisations.

## Recover the per-region random effects

The internal column-stacking convention is $\alpha = \mathrm{vec}_{\text{col}}(M)$
where $M$ is $G \times 2$ — the first column holds origin effects, the
second holds destination effects.

```{r ranef}
alpha_hat <- fit$alpha
M_hat <- matrix(alpha_hat, nrow = G, ncol = 2L)
colnames(M_hat) <- c("origin", "dest")
rownames(M_hat) <- paste0("region_", seq_len(G))
round(M_hat, 3)
```

A quick visual of the two effect series. Notice how `origin` and `dest`
track each other closely across regions — that's the visual signature of
the high empirical OD correlation discussed above.

```{r ranef-plot, fig.cap = "Estimated origin and destination effects per region."}
matplot(seq_len(G), M_hat, type = "b", pch = 19, lty = 1, lwd = 2,
        col = c("steelblue", "darkorange"),
        xlab = "Region", ylab = "Random effect")
abline(h = 0, lty = 2, col = "grey60")
legend("topright", c("origin", "dest"),
       col = c("steelblue", "darkorange"), pch = 19, lwd = 2, bty = "n")
```

## Estimated varying coefficient

```{r vcoef, fig.cap = "Estimated time-varying effect of wage_diff on log flow."}
t_grid <- seq(0, 1, length.out = 100L)
vc     <- varying_coef(fit, t_new = t_grid, k = 1L, se.fit = TRUE)

plot(t_grid, vc$fit, type = "l", lwd = 2, col = "steelblue",
     xlab = "t (year, normalised)",
     ylab = expression(hat(beta)[1](t)),
     ylim = range(vc$fit - 2 * vc$se.fit, vc$fit + 2 * vc$se.fit,
                  1.5 * sin(2 * pi * t_grid)))
polygon(c(t_grid, rev(t_grid)),
        c(vc$fit + 2 * vc$se.fit, rev(vc$fit - 2 * vc$se.fit)),
        col = adjustcolor("steelblue", alpha.f = 0.25), border = NA)
lines(t_grid, 1.5 * sin(2 * pi * t_grid),
      col = "red", lty = 2, lwd = 2)
legend("topright", c("estimate", "truth"),
       col = c("steelblue", "red"), lty = c(1, 2), lwd = 2, bty = "n")
```

The wage-difference effect oscillates with time: in the simulation it
follows $1.5\sin(2\pi t)$. Unlike $\Sigma_{2\times 2}$, the
varying-coefficient curve is supported by all $N = 3000$ observations and
tracks the truth essentially perfectly across the 30-year window.

## Distributing this fit

The same OD problem can be fit in the distributed setting. Each node
computes a `node_summary()` on its slice and ships the result; the
central node aggregates and calls `fit_from_summaries()` with
**`rowsum_constant = 2`** so the identifiability shift matches what
`vcmm()` applies automatically:

```{r distributed-od, eval = FALSE}
# Split N obs across 3 nodes
node_id <- sample.int(3L, N, replace = TRUE)
splits  <- split(seq_len(N), node_id)

design <- build_vcmm_design(X = od$wage_diff, t = od$t)
X_d    <- design$X_design

summaries <- lapply(splits, function(idx)
  node_summary(od$log_flow[idx],
               X_d[idx, , drop = FALSE],
               Z[idx, , drop = FALSE]))

fit_dist <- fit_from_summaries(
  summaries,
  penalty            = design$penalty,
  control            = vcmm_control(sigma_eps = 0.6,
                                    sigma_alpha = sqrt(0.5),
                                    update_variance = TRUE),
  method             = "csl",
  re_cov             = "kronecker",
  n_groups           = G,
  Sigma_spatial_init = Sigma_spatial,
  rowsum_constant    = 2
)

# Same answer as the pooled fit above, up to BLAS noise
max(abs(fit$beta  - fit_dist$beta))
max(abs(fit$alpha - fit_dist$alpha))
```

See `vignette("distributed-fitting", package = "cevcmm")` for the full
explanation of the distributed API.

## Where to go next

- **Basic usage** with diagonal random effects:
  `vignette("getting-started", package = "cevcmm")`.
- **Distributed fitting** in detail:
  `vignette("distributed-fitting", package = "cevcmm")`.

## Reference

Jalili, L. and Lin, L.-H. (2025). *Scalable and Communication-Efficient Varying Coefficient Mixed Effect Models: Methodology, Theory, and Applications.* arXiv:2511.12732; under review at *Journal of the American Statistical Association*.
