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