OD migration with Kronecker covariance

In cevcmm: Communication-Efficient Varying Coefficient Mixed-Effects Models

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

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.

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

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:

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

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

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.

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.

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

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:

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

cevcmm documentation built on July 24, 2026, 5:07 p.m.