Coarse-to-fine dynamic space-time modeling: An application example

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.1154274  1.4256394  1.8781148
#> cos12     -0.5160965 0.1617919 -0.8332087 -0.1989843
#> 
#> ---- 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.83837195
#> 2      validation_RMSE 2.93591766
#> 3       validation_MAE 0.02263026

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. Interior gaps are interpolated and time points beyond the last observed one are forecast, in both cases through the AR(1) predict step of each scale. 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 
#> 4.538025 4.533792 6.872447

The average predictive SD is markedly larger at the forecast month than at the two observed months. Adding a future time point does not alter the fit itself: the accepted bandwidths and the training-time predictions are unchanged.

Multiscale space-time pattern extraction

As in the spatial case, sp_scalewise synthesizes the scale-wise processes whose bandwidths fall inside a given range. For a cf_dglm fit it additionally averages over a time window, given by time_range, and returns one row per location:

mod_l <- sp_scalewise(mod, bw_range = c(150000, Inf))   # large scale (>= 150 km)
mod_s <- sp_scalewise(mod, bw_range = c(0, 150000))     # small scale (< 150 km)
head(mod_l$pred, 3)
#>         px      py     pred    pred_sd n_time
#> 1 538708.5 5947030 2.724960 0.06305079     59
#> 2 545413.6 5930802 2.470467 0.06237242     60
#> 3 551796.2 5991947 2.853334 0.06719856     50

The n_time column records how many time points were averaged at each location. Using time_range, the large-scale process can be compared between seasons; here we average the three winter months of 2005 and the three summer months of the same year:

win <- sp_scalewise(mod, bw_range = c(150000, Inf), time_range = c(49, 51))
smr <- sp_scalewise(mod, bw_range = c(150000, Inf), time_range = c(55, 57))
sea <- st_as_sf(data.frame(win$pred[, c("px", "py")],
                           winter = win$pred$pred,
                           summer = smr$pred$pred),
                coords = c("px", "py"))
plot(sea[, c("winter", "summer")], pch = 20, cex = 1.2,
     axes = TRUE, key.pos = 4, nbreaks = 20)

The large-scale process is far more structured in winter than in summer, which is consistent with the seasonal contrast seen in the predictive maps. Note that time_range must contain at least one time point of each site set present in the fit: because the prediction grid above was requested only at months 49, 55 and 61, a window such as c(1, 12) would select no prediction site and raise an error.

Setting time_range to its default, c(-Inf, Inf), averages over all 60 months and gives the time-averaged multiscale decomposition, which is the direct space-time analogue of the decomposition produced for cf_lm and cf_glm fits.

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.

spCFmap()

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.