This vignette demonstrates space-time prediction and multiscale
analysis using the dynamic extension of the coarse-to-fine spatial
modeling (CFSM) framework. As in the purely spatial case (Murakami et al. 2026b, 2026a), the space-time
process is a sum of scale-wise processes learned sequentially from
coarser to finer scales, and the number of scales is selected by holdout
validation. See Murakami (2026) for the
dynamic extension used here. Each scale additionally couples a per-knot
AR(1) process in time with kernel kriging in space, so the same model
yields maps at individual time points, temporally averaged patterns, and
forecasts at time points that were never observed.
Data and setup
Let us load the required packages
library(spCF)
library(sf)
#> Linking to GEOS 3.14.1, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
The example data are the monthly air quality panel shipped with the
package, which is also the Demo (air, space-time) data set of
the spCFmap() application. It contains monthly mean PM10
concentrations (µg/m³) recorded at 63 background monitoring stations in
Germany from January 2001 to December 2005:
air <- read.csv(system.file("shiny", "spCFmap", "example_spacetime_air.csv",
package = "spCF"))
str(air)
#> 'data.frame': 2769 obs. of 7 variables:
#> $ station: chr "DESH001" "DESH001" "DESH001" "DESH001" ...
#> $ lon : num 9.59 9.59 9.59 9.59 9.59 ...
#> $ lat : num 53.7 53.7 53.7 53.7 53.7 ...
#> $ year : int 2001 2001 2001 2001 2001 2001 2001 2001 2001 2001 ...
#> $ month : int 1 2 3 4 5 6 7 8 9 10 ...
#> $ time : int 1 2 3 4 5 6 7 8 9 10 ...
#> $ pm10 : num 32 21.1 27.4 18.9 20.8 ...
The time column indexes the 60 months (1 = January 2001,
60 = December 2005). The panel is unbalanced: stations enter and leave
the network, so only 2,769 of the 63 x 60 = 3,780 possible station-month
combinations are observed. CFSM does not require a balanced panel, nor
does it require the same locations to be observed at every time
point.
range(table(air$station)) # observations per station
#> [1] 1 60
The station coordinates are given as longitude/latitude. Because the
kernels used by CFSM are distance based, we project them to UTM zone 32N
(EPSG:25832), whose unit is the metre, so that the selected bandwidths
are readable as distances:
pts <- st_as_sf(air, coords = c("lon", "lat"), crs = 4326)
coords <- st_coordinates(st_transform(pts, 25832))
colnames(coords) <- c("px", "py")
PM10 has a pronounced annual cycle, which is a fixed effect rather
than a spatial one. We therefore describe it with a pair of seasonal
harmonics and let the space-time process capture the remaining
structure:
y <- air$pm10 # response
time <- air$time # time index
x <- data.frame(sin12 = sin(2 * pi * air$month / 12),
cos12 = cos(2 * pi * air$month / 12))
The station means over the five years are plotted as follows:
loc <- air[!duplicated(air$station), c("station", "lon", "lat")]
loc$pm10_ave <- tapply(air$pm10, air$station, mean)[loc$station]
loc_sf <- st_as_sf(loc, coords = c("lon", "lat"), crs = 4326)
plot(loc_sf[, "pm10_ave"], pch = 20, cex = 1.3, axes = TRUE,
key.pos = 4, nbreaks = 20)
Coarse-to-fine dynamic spatial GLMM (CF-DGLM)
The cf_dglm_hv function performs the holdout validation
that selects the number of spatial scales. It is used exactly like
cf_lm_hv and cf_glm_hv, with the addition of
the time argument:
mod_hv <- cf_dglm_hv(y = y, x = x, coords = coords, time = time)
#> --- Validation deviance: Basic GLM ---
#> 37460.02
#> --- Validation deviance: Learning multi-scale space-time process ---
#> 24308.55 (Scale 1)
#> 22943.82 (Scale 2)
#> 21737.48 (Scale 3)
#> 20565.34 (Scale 4)
#> 19587.77 (Scale 5)
#> 18719.19 (Scale 6)
#> 17972.16 (Scale 7)
#> 17365.7 (Scale 8)
#> 16885.19 (Scale 9)
#> 16515.12 (Scale 10)
#> 16234.53 (Scale 11)
#> 16016.64 (Scale 12)
#> 15841.46 (Scale 13)
#> 15695.02 (Scale 14)
#> 15582.85 (Scale 15)
#> 15518.05 (Scale 16)
#> 15518.03 (Scale 17)
#> 15518.03 (Scale 18) no improvement
#> 15518.03 (Scale 19) no improvement
#> 15518.03 (Scale 20) no improvement
#> 15518.03 (Scale 21) no improvement
#> 15518.03 (Scale 22) no improvement
As the output shows, the validation deviance decreases as learning
proceeds from the coarsest scale to finer ones and the loop terminates
once further scales stop improving it. The temporal parameters (the
AR(1) coefficient and its innovation variance) are estimated inside this
step; rho and Q may also be fixed by the
user.
Next we define the sites at which predictions are required. Here we
use a regular 25 km grid covering the convex hull of the monitoring
network:
uni <- unique(as.data.frame(coords))
hull <- st_convex_hull(st_union(st_as_sf(uni, coords = c("px", "py"))))
gcen <- st_make_grid(hull, cellsize = 25000, what = "centers")
gcen <- gcen[st_intersects(gcen, hull, sparse = FALSE)[, 1]]
gxy <- st_coordinates(gcen)
nrow(gxy)
#> [1] 585
Every prediction site needs a time point as well as coordinates. We
ask for three of them: January 2005 (time = 49) and July
2005 (time = 55), which are both observed, and January 2006
(time = 61), which lies beyond the end of the data and is
therefore a forecast. The seasonal covariates are evaluated at the
corresponding months:
tp <- c(49, 55, 61) # Jan-2005, Jul-2005, Jan-2006
month0 <- c(1, 7, 1)
ng <- nrow(gxy)
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))
The full model is then trained with cf_dglm:
mod <- cf_dglm(y = y, x = x, coords = coords, time = time,
x0 = x0, coords0 = coords0, time0 = time0, mod_hv = mod_hv)
The estimated regression coefficients, the standard deviations of the
model elements, and the error statistics are displayed as follows:
mod
#> Call:
#> cf_dglm(y = y, x = x, coords = coords, time = time, x0 = x0,
#> coords0 = coords0, time0 = time0, mod_hv = mod_hv)
#>
#> ---- Coefficients -------------------------------------
#> coef coef_se lower_95CI upper_95CI
#> Intercept 18.7435597 0.2964365 18.1625441 19.3245753
#> sin12 1.6518771 0.1133397 1.4297313 1.8740229
#> cos12 -0.5160965 0.1612323 -0.8321119 -0.2000811
#>
#> ---- Space-time parameters ----------------------------
#> parameter value
#> Temporal autocorrelation, AR(1) (rho) 0.2315
#> Temporal innovation variance (Q) 18.6345
#>
#> ---- Standard deviations (model elements) -------------
#> elements standard_deviation
#> 1 xb 1.2261484
#> 2 spatial_scale1 4.3604285
#> 3 spatial_scale2 0.3406890
#> 4 spatial_scale3 0.3233307
#> 5 spatial_scale4 0.3632474
#> 6 spatial_scale5 0.3541019
#> 7 spatial_scale6 0.3450485
#> 8 spatial_scale7 0.3198243
#> 9 spatial_scale8 0.2734035
#> 10 spatial_scale9 0.2584766
#> 11 spatial_scale10 0.2534986
#> 12 spatial_scale11 0.2578029
#> 13 spatial_scale12 0.2725394
#> 14 spatial_scale13 0.2973643
#> 15 spatial_scale14 0.3287485
#> 16 spatial_scale15 0.3597883
#> 17 spatial_scale16 0.3826629
#> 18 spatial_scale17 0.3899164
#>
#> ---- Error statistics ---------------------------------
#> stat value
#> 1 validation_Pseudo-R2 0.8383719
#> 2 validation_RMSE 2.9359177
#> 3 validation_MAE 2.0160181
The seasonal coefficients confirm the expected annual cycle, and the
scale-wise standard deviations show that most of the space-time
variation is carried by the coarsest scales. The space-time parameters
printed above are specific to cf_dglm: rho is
the AR(1) coefficient of the scale-wise processes at the monthly step
(values near one indicate a persistent field, values near zero a field
that is re-drawn every month once the seasonal fixed effect has been
removed), Q is its innovation variance, and
tau is the holdout-calibrated multiplier applied to the
field variance. They are also available individually as
mod$other$rho, mod$other$Q and
mod$other$tau.
The selected bandwidths are stored in the fitted object:
round(mod$bands / 1000, 1) # bandwidth of each accepted scale, in km
#> [1] 297.7 268.0 241.2 217.1 195.3 175.8 158.2 142.4 128.2 115.4 103.8 93.4
#> [13] 84.1 75.7 68.1 61.3 55.2
Space-time prediction
Predictive values at observed time points
The predictive means at the grid cells are extracted by subsetting
pred0 with time0. January and July of 2005 are
mapped side by side:
grid_sf <- st_as_sf(as.data.frame(gxy), coords = c("X", "Y"))
grid_sf$Jan2005 <- mod$pred0$pred[time0 == 49]
grid_sf$Jul2005 <- mod$pred0$pred[time0 == 55]
plot(grid_sf[, c("Jan2005", "Jul2005")], pch = 15, cex = 1.9,
axes = TRUE, key.pos = 4, nbreaks = 20)
The January map is far more structured than the July one, with a
pronounced maximum in the north-west of the network, whereas July is
close to flat. The annual cycle itself is carried by the seasonal
covariates, so what the space-time process adds on top of it is a
winter-specific spatial gradient rather than a uniform seasonal shift.
The predictive standard deviations are mapped in the same way:
grid_sf$Jan2005_sd <- mod$pred0$pred_sd[time0 == 49]
plot(grid_sf[, "Jan2005_sd"], pch = 15, cex = 1.9, axes = TRUE, key.pos = 4,
pal = function(n) hcl.colors(n, "Viridis"))
Uncertainty is smallest near the monitoring stations and grows
towards the edges of the network, as expected.
Forecasting an unobserved time point
time0 may contain time points that carry no
observations. They are predicted from the smoothed AR(1) states of the
knots of each scale, without changing the fit: a time point between two
observed ones is bridged between their states, and a time point after
the last observed one is forecast, by as many AR(1) steps as months have
passed. January 2006 was requested above, so its map is already
available:
grid_sf$Jan2006 <- mod$pred0$pred[time0 == 61]
plot(grid_sf[, c("Jan2005", "Jan2006")], pch = 15, cex = 1.9,
axes = TRUE, key.pos = 4, nbreaks = 20)
The forecast keeps the seasonal fixed effect but the spatial field is
pulled towards its mean, because no data are available to update it. The
predictive standard deviation reflects this honestly:
tapply(mod$pred0$pred_sd, time0, mean)
#> 49 55 61
#> 5.108730 5.104234 6.772798
The average predictive SD is markedly larger at the forecast month
than at the two observed months. Requesting such time points does not
alter the fit itself: the accepted bandwidths and the training-time
predictions are unchanged.
Prediction at new sites and times with predict()
The prediction sites and times need not be given when the model is
fitted. The model can be fitted without x0,
coords0 and time0, and predict()
then predicts at any sites and time points later from the fitted object
alone, without the training data. The result is identical to giving the
same sites and times to cf_dglm(). For example, the
forecast for March 2006 (time = 63) at the grid cells, with
a 90% prediction interval:
mod_f <- cf_dglm(y = y, x = x, coords = coords, time = time,
mod_hv = mod_hv) # no x0, coords0, time0
p63 <- predict(mod_f, coords0 = gxy, time0 = rep(63, ng),
x0 = data.frame(sin12 = rep(sin(2 * pi * 3 / 12), ng),
cos12 = rep(cos(2 * pi * 3 / 12), ng)),
probs = c(0.05, 0.95))
predict() can be applied to mod above as
well, which avoids refitting here:
p63 <- predict(mod, coords0 = gxy, time0 = rep(63, ng),
x0 = data.frame(sin12 = rep(sin(2 * pi * 3 / 12), ng),
cos12 = rep(cos(2 * pi * 3 / 12), ng)),
probs = c(0.05, 0.95))
head(p63)
#> pred pred_sd q0.05 q0.95
#> 1 18.68621 6.869918 7.386202 29.98622
#> 2 18.70710 6.869978 7.406992 30.00721
#> 3 18.72815 6.870021 7.427966 30.02832
#> 4 18.74675 6.870054 7.446514 30.04698
#> 5 18.76139 6.870079 7.461116 30.06166
#> 6 18.77242 6.870100 7.472111 30.07273
Predictive quantiles and the size of the fitted object
mod$pred_q and mod$pred0_q give the
predictive quantiles at the levels 0.005, 0.025, 0.05, 0.1, …, 0.9,
0.95, 0.975, 0.995 (other levels with
predict(..., probs = )); for example a 90% interval at the
forecast month:
head(mod$pred0_q[time0 == 61, c("q0.05", "q0.95")])
#> q0.05 q0.95
#> 1171 3.593651 25.78368
#> 1172 3.948868 26.20798
#> 1173 4.316480 26.62570
#> 1174 4.644991 26.99130
#> 1175 4.903617 27.27908
#> 1176 5.097292 27.49697
A cf_dglm fit holds its predictions for every site and
time point, so for long panels the fitted object can become large. The
quantile tables are therefore not stored but computed when accessed, and
the scale-wise processes Z, Z_sd,
Z0 and Z0_sd, which take one column per
selected scale for every observation, can be dropped with
keep_scales = FALSE. Predictions, standard deviations and
sd_summary are unchanged; only sp_scalewise()
below needs the scale-wise processes:
mod_small <- cf_dglm(y = y, x = x, coords = coords, time = time,
x0 = x0, coords0 = coords0, time0 = time0,
mod_hv = mod_hv, keep_scales = FALSE)
Interactive mapping
The maps above are static. The same results can be explored over a
basemap with spCFmap(), which has two modes.
Called without arguments, it opens the full application, where this
data set is available as Demo (air, space-time). The response
(pm10), the coordinates (lon,
lat) and the time index (time) are
pre-selected, so the model can be fitted, mapped and exported as CSV or
GeoJSON without writing code. The seasonal harmonics used above are not
columns of the file, so the app fits the space-time process without
covariates; user data can of course be uploaded with any covariates
already present as columns.
Passing a fitted model maps that object directly, with the layer
(predictive mean, predictive SD, covariate effect, or a scale-wise
component), the colour scaling and — for a cf_dglm fit —
the time range and the bandwidth range chosen interactively:
spCFmap(mod, crs = 25832)
The crs argument states the reference system of the
coordinates that were passed to the model, so that the results can be
placed on the longitude/latitude basemap. An EPSG code is enough as long
as the coordinates are in the unit of that CRS, which is why the
projected coordinates were kept in metres above; a proj/WKT string may
be given instead, for instance if the coordinates have been
rescaled.
Reference
Murakami, Daisuke. 2026. Fast Covariance-Free Spatiotemporal
Modeling via Coarse-to-Fine Learning. ArXiv.
Murakami, Daisuke, Alexis Comber, Takahiro Yoshida, Narumasa Tsutsumida,
Chris Brunsdon, and Tomoki Nakaya. 2026a. Coarse-to-Fine Spatial
GLMM for Scalable Prediction and Multiscale Analysis. ArXiv.
Murakami, Daisuke, Alexis Comber, Takahiro Yoshida, Narumasa Tsutsumida,
Chris Brunsdon, and Tomoki Nakaya. 2026b. “Coarse-to-Fine Spatial
Modeling: A Scalable, Machine-Learning-Compatible Framework.”
Geographical Analysis 58 (2): e70034.
