Getting started with cevcmm

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(1)

What this package fits

cevcmm fits Varying Coefficient Mixed-Effects Models (VCMMs):

$$ y_i = \beta_0(t_i) + \sum_{k=1}^{K} x_{ik}\,\beta_k(t_i) + z_i^\top \alpha + \varepsilon_i, $$

where each $\beta_k(t)$ is a smooth function of $t$ (cubic B-splines with a second-order-difference penalty), $\alpha \sim N(0, \Sigma_\alpha)$ are random effects with one of three covariance structures (diag, kronecker, separable), and $\varepsilon_i \sim N(0, \sigma_\varepsilon^2)$.

The package implements two estimators from Jalili and Lin (2025):

• SS (sufficient-statistics): iterates to convergence using only the aggregated summary $(a,\,b,\,C,\,Z^\top Z,\,Z^\top y,\,X^\top Z)$. Raw data is not needed once the summaries exist.
• CSL (one-step communication-efficient): a single Newton refinement from a pilot estimator. Asymptotically equivalent to SS at a fraction of the cost when the pilot is good enough.

This vignette walks through the simplest case: a single varying coefficient, diagonal random effects, all on one machine.

Setup

library(cevcmm)

Simulate a small dataset

$N = 500$ observations, $K = 1$ varying covariate, $q = 3$ independent random effects.

set.seed(1)
N <- 500L
q <- 3L

t <- runif(N)
x <- runif(N)
Z <- matrix(rnorm(N * q), N, q)

beta_0     <- 2
beta_1_fun <- function(u) sin(2 * pi * u)
alpha_true <- rnorm(q, sd = 0.4)
sigma_eps  <- 0.5

y <- beta_0 + beta_1_fun(t) * x +
     as.vector(Z %*% alpha_true) +
     rnorm(N, sd = sigma_eps)

Fit

A single vcmm() call. With method = "auto" (default), CSL is picked when $N \cdot q > 10^5$ or $q > 50$; otherwise SS. This small example picks SS.

fit <- vcmm(y, X = x, Z = Z, t = t,
            method  = "auto",
            re_cov  = "diag",
            control = vcmm_control(sigma_eps       = 0.5,
                                   sigma_alpha     = 0.4,
                                   update_variance = TRUE))
fit

Inspect the fit

All standard S3 methods work on a vcmm_fit.

# Fixed-effects coefficient vector (intercept + spline basis coefs)
head(coef(fit))

# Same vector, reshaped: intercept + (basis x covariate) matrix
fx <- fixef(fit)
fx$intercept
dim(fx$varying)

# Random effects
ranef(fit)

# Sample size and residual SD
nobs(fit)
fit$sigma_eps

# Log-likelihood, AIC, BIC
logLik(fit)
AIC(fit); BIC(fit)

A full summary() includes a coefficient table with z-tests:

summary(fit)

Recover the varying coefficient

varying_coef() evaluates $\hat\beta_k(t)$ at any grid with optional pointwise standard errors.

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", ylab = expression(hat(beta)[1](t)),
     ylim = range(vc$fit - 2 * vc$se.fit,
                  vc$fit + 2 * vc$se.fit,
                  beta_1_fun(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, beta_1_fun(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 estimate tracks the true $\sin(2\pi t)$ across the unit interval; the band is the pointwise 95% Wald interval.

Predict on new data

new_idx <- sample.int(N, 5L)
newdata <- list(t = t[new_idx],
                X = x[new_idx],
                Z = Z[new_idx, , drop = FALSE])

predict(fit, newdata = newdata)
y[new_idx]

With se.fit = TRUE:

pred <- predict(fit, newdata = newdata, se.fit = TRUE)
cbind(fit = pred$fit, se = pred$se.fit)

Built-in diagnostic plots

plot(fit, which = 1) shows the varying coefficient with its confidence band; which = 2 requires the training data and shows residual diagnostics; which = 3 shows random-effect diagnostics. See ?plot.vcmm_fit.

plot(fit, which = 1)

Where to go next

• Distributed fitting — when data lives on multiple nodes or won't fit in memory: vignette("distributed-fitting", package = "cevcmm").
• Origin-destination migration — Kronecker covariance for OD-flow data: vignette("od-migration", 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.