cf_dglm: Coarse-to-fine dynamic (space-time) spatial GLMMs (CF-DGLMMs)

In spCF: Coarse-to-Fine Spatial and Spatio-Temporal Modeling

Coarse-to-fine dynamic (space-time) spatial GLMMs (CF-DGLMMs)

Description

Prediction and regression via a separable space-time cascade. Given the scales selected by cf_dglm_hv, the model is refitted on the full sample and predictions (with standard deviations) are produced at sample and, optionally, prediction sites. The link-scale linear predictor is g(\mu_{i,t}) = x_{i,t}'\beta + \sum_k f_k(s_i,t) + offset, where each scale-k field f_k couples a per-knot AR(1) Kalman smoother in time with kernel kriging in space.

Usage

cf_dglm(
  y,
  x = NULL,
  coords,
  time,
  offset = NULL,
  x0 = NULL,
  coords0 = NULL,
  time0 = NULL,
  offset0 = NULL,
  mod_hv,
  robust_se = TRUE,
  se_type = c("prediction", "mean"),
  se_method = c("opt", "classic"),
  keep_scales = TRUE
)

Arguments

y	Vector of response variables (N x 1).
x	Matrix of covariates (N x K).
coords	Matrix of 2-dimensional point coordinates (N x 2). The space-time panel may be unbalanced (observed locations may differ across time points).
time	Vector of time indices (N x 1); must use the same time points as in cf_dglm_hv.
offset	Optional. Offset variable (N x 1), consistent with glm.
x0	Optional. Matrix of covariates at prediction sites (N0 x K).
coords0	Optional. Coordinates at prediction sites (N0 x 2).
time0	Optional. Time points at prediction sites (N0 x 1). May include time points with no observations. They are predicted, after the fit and without changing it, from the smoothed per-knot AR(1) states: a time point between two training time points is bridged between their states, the step being split in proportion to the time differences; a time point after (before) the training period is forecast (backcast), the number of AR(1) steps being the time difference over the median spacing of the training time points. The same applies to the time-varying coefficients (random walk). The fit, and beta_tv, therefore do not depend on time0, and predict.cf_dglm gives the same predictions for new sites and times later.
offset0	Optional. Offset at prediction sites (N0 x 1).
mod_hv	Output object of cf_dglm_hv.
robust_se	Logical; if TRUE (default), the constant-coefficient standard errors (and the coefficient-uncertainty term of the predictive SE) use a spatial-block cluster-robust sandwich that accounts for the cascade field being a correlated random component. The naive model-based covariance treats the field as a known offset and severely understates the SEs; the robust version restores near-nominal coverage. Set FALSE for the naive vcov(glm) SEs.
se_type	Type of predictive uncertainty in pred/pred_q. "prediction" (default) returns the holdout-calibrated OBSERVATION predictive for a new data point (Gaussian: mean uncertainty + residual variance; Poisson: negative-binomial count predictive; binomial: temperature-calibrated probability with pred_sd = sqrt(p(1-p))). Negative binomial (negbin): negative-binomial count predictive with the fitted dispersion. Other families use a moment-matched observation predictive with the holdout-estimated dispersion (Gamma: gamma; inverse.gaussian: inverse Gaussian; quasipoisson: negative binomial; quasibinomial: beta for proportions, as binomial for 0/1 data; otherwise normal), with the mean-uncertainty scale calibrated to 95% holdout coverage. The signal versions are kept in pred_signal/pred_q_signal. "mean" returns the signal (mean) uncertainty only (previous behaviour). See other$calibration.
se_method	Cluster-robust coefficient-SE estimator (used when robust_se = TRUE). "opt" (default) splits the sandwich meat into a field-removed observation-noise part and a field part that adds the calibrated field variance back with a within-block exp(-d/h) correlation (h = median committed bandwidth); this is near-nominal. A refit-free leverage leave-one-out ceiling then caps the field term, preventing over-coverage for count (Poisson) responses while leaving already-calibrated families unchanged. "classic" keeps the realised field inside the working residual (the previous behaviour), which is valid but conservative.
keep_scales	If TRUE (default), the scale-wise processes Z, Z_sd, Z0 and Z0_sd are kept in the output. They hold one column per selected scale for every observation (site x time point), which makes them the largest part of the fitted object for large data, and above all for long panels. With FALSE they are dropped (NULL); predictions, standard errors, sd_summary and the maps of spCFmap for the total prediction are unchanged, but sp_scalewise needs them.

Details

The full-sample fit is a SINGLE coarse-to-fine cascade sweep, mirroring the relationship between cf_glm and cf_glm_hv: it reuses the same single greedy sweep that cf_dglm_hv performs for scale selection, plus prediction. Within the sweep, for each band (coarse to fine) the GLM working response/weights are refreshed (IRLS folded into the sweep, as cf_glm's per-band glm() does), the scale is fit and accumulated, and the constant and time-varying coefficients are backfit. (The earlier outer-IRLS implementation is archived as cf_dglm_iter under misc/.)

The field variance of the mean is assembled scale by scale so that it grows smoothly with the distance to the data and reaches the marginal variance of the fitted total field (the sill, var(sum_r z_r) on the link scale) far from it, as a stationary process reverts to its marginal variance. For scale r, r_r = V^d_r / P_{0,r} is the fraction of the scale's prior variance left after the data, from a distance-aware gPoE variance (each knot informs a site through its kernel correlation w, conditional variance w^2 P + (1 - w^2) P_0); the calibrated variance is c_r \tau r_r / (1 + r_r(\tau - 1)), where the caps c_r are proportional to the variance of each fitted scale and sum to the sill, and \tau (mod_hv$other$tau_stage) scales the information of the data and is solved from the holdout moment equation in cf_dglm_hv. The sill is floored at a direct estimate of the field variance (working residual variance of the GLM minus a nearest-neighbour nugget), and when no scale is accepted that estimate is added as unmodeled field variance. Point predictions and coefficient estimates do not depend on these bounds. For binomial responses the field variance is left uncapped.

Value

A list (class "cf_dglm") mirroring cf_glm: beta, sd_summary, e_summary, pred, pred0, pred_q, pred0_q, bands, Z, Z_sd, Z0, Z0_sd, other, call, plus

beta_tv, beta_tv_sd

Time-varying coefficients and their standard deviations, one row per training time point and one column per covariate named in tvc (plus a time column). NULL when tvc was not used in cf_dglm_hv.

pred_signal, pred_q_signal

The signal (mean) predictive kept alongside the observation predictive when se_type = "prediction".

As in cf_glm, the quantile tables pred_q, pred0_q and pred_q_signal are not stored but computed on access (at the 15 levels 0.005, 0.025, 0.05, 0.1, ..., 0.9, 0.95, 0.975, 0.995; predict.cf_dglm gives other levels), and Z, Z_sd, Z0, Z0_sd are NULL when keep_scales = FALSE. The temporal parameters of the fitted cascade are in other$rho (AR(1) autocorrelation), other$Q (innovation variance) and other$tau (holdout-calibrated field-variance factor); the first two are shown by print.

Author(s)

Daisuke Murakami

References

Murakami, D. (2026). Fast covariance-free spatiotemporal modeling via coarse-to-fine learning. *ArXiv preprint*, 2608.03449.

See Also

cf_dglm_hv, cf_glm

Examples

### Monthly PM10 at 63 German background stations, 2001-2005 (the data set
### behind the "Demo (air, space-time)" entry of spCFmap(); see the
### spCF_dglm vignette for a fuller walk-through).
require(sf)
air    <- read.csv(system.file("shiny", "spCFmap",
                               "example_spacetime_air.csv", package = "spCF"))
pts    <- st_as_sf(air, coords = c("lon", "lat"), crs = 4326)
coords <- st_coordinates(st_transform(pts, 25832))  # UTM 32N, in metres

### The annual cycle is a fixed effect; the space-time process takes the rest
x      <- data.frame(sin12 = sin(2 * pi * air$month / 12),
                     cos12 = cos(2 * pi * air$month / 12))

### Holdout validation optimizing the number of spatial scales
mod_hv <- cf_dglm_hv(y = air$pm10, x = x, coords = coords, time = air$time)

### Prediction sites: a regular 25 km grid covering the convex hull of the
### network (as in the spCF_dglm vignette), at 10 time points every six months
### over the observed period: June 2001 (time = 6), December 2001 (12), ...,
### December 2005 (60)
uni     <- unique(as.data.frame(coords))
hull    <- st_convex_hull(st_union(st_as_sf(uni, coords = c("X", "Y"))))
gcen    <- st_make_grid(hull, cellsize = 25000, what = "centers")
gcen    <- gcen[st_intersects(gcen, hull, sparse = FALSE)[, 1]]
gxy     <- st_coordinates(gcen)
ng      <- nrow(gxy)
tp      <- seq(6, 60, by = 6)
month0  <- rep(c(6, 12), length.out = length(tp))  # calendar month (June, December)
coords0 <- do.call(rbind, replicate(length(tp), gxy, simplify = FALSE))
time0   <- rep(tp, each = ng)
x0      <- data.frame(sin12 = sin(2 * pi * rep(month0, each = ng) / 12),
                      cos12 = cos(2 * pi * rep(month0, each = ng) / 12))

### Space-time modeling and prediction
mod    <- cf_dglm(y = air$pm10, x = x, coords = coords, time = air$time,
                  x0 = x0, coords0 = coords0, time0 = time0, mod_hv = mod_hv)
mod

round(mod$bands / 1000, 1)              # accepted bandwidths, in km
round(c(rho = mod$other$rho, Q = mod$other$Q), 3)  # AR(1) parameters

### Mapping the predictions for June 2005 and December 2005
grid_sf <- st_sf(Jun2005 = mod$pred0$pred[time0 == 54],
                 Dec2005 = mod$pred0$pred[time0 == 60],
                 geometry = gcen, crs = 25832)
plot(grid_sf, pch = 15, cex = 1.5, axes = TRUE, key.pos = 4, nbreaks = 20)

### Multiscale extraction, averaged over the observed months
mod_s1 <- sp_scalewise(mod, bw_range = c(150000, Inf))  # large scale
mod_s2 <- sp_scalewise(mod, bw_range = c(0, 150000))    # small scale

### The same fit, explored interactively over a basemap
# spCFmap(mod, crs = 25832)

### Prediction with predict(): the model can be fitted WITHOUT prediction
### sites and times (no x0, coords0, time0) and used to predict at any sites
### and time points later, without the training data: e.g. the grid in
### December 2005 (as above) and in March 2006 (time = 63, three months beyond
### the data), with a 90 percent prediction interval
mod_f  <- cf_dglm(y = air$pm10, x = x, coords = coords, time = air$time,
                  mod_hv = mod_hv)                  # no x0, coords0, time0
p60    <- predict(mod_f, x0 = x0[time0 == 60, ], coords0 = gxy,
                  time0 = rep(60, ng))
all.equal(p60$pred, mod$pred0$pred[time0 == 60])   # same as cf_dglm(..., time0)
p      <- predict(mod_f, x0 = data.frame(sin12 = rep(sin(2 * pi * 3 / 12), ng),
                                        cos12 = rep(cos(2 * pi * 3 / 12), ng)),
                  coords0 = gxy, time0 = rep(63, ng), probs = c(0.05, 0.95))
head(p)

spCF documentation built on Oct. 5, 2026, 5:07 p.m.