--- title: "The GCF workflow: from spatial variables to better predictions" author: "Yongze Song" output: rmarkdown::html_vignette: toc: true vignette: > %\VignetteIndexEntry{The GCF workflow: from spatial variables to better predictions} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4 ) ``` The generalized covariate field (GCF) method expands each spatial covariate into two complementary families of features -- spatial-pattern features and neighbourhood-distribution features -- and selects a stable subset of them for geospatial prediction. GCF is prediction-oriented feature construction: it enriches the covariate space rather than the model structure, and the selected variables feed any downstream regression learner. This vignette walks through the full workflow on the paper's simulation data: 1. generate the GCF variables with `gcf_field()`; 2. select a stable subset with `gcf_select()`; 3. compare a machine-learning model fitted on the raw covariates against the same model fitted on the GCF variables, under both random and spatial-block cross-validation, using the external **randomForest** package (the learner adopted in the paper). To keep the build time of this vignette short, the selection uses B = 10 stability resamples and reduced tree counts; the paper's full settings are B = 80 resamples with 200-tree importance kernels and 500-tree final forests, and give the same qualitative conclusions. ```{r setup} library(gcf) has_rf <- requireNamespace("randomForest", quietly = TRUE) ``` ## The simulation data `sim_grid` is the paper's simulation dataset: a 30 x 30 regular grid (900 cells, unit spacing) with a response `y1` and three spatially structured covariates `x1`, `x2`, `x3`. ```{r data} data(sim_grid) head(sim_grid) ``` ## Step 1-3: generate the GCF variables `gcf_field()` maps the covariates to their GCF variables in three steps: spatial-pattern features (`gcf_psi()`, 11 operators over buffer radii), neighbourhood-distribution features (`gcf_zx()`, buffer-wise quantiles), and a functional reduction (`gcf_reduce()`) that collapses the collinear buffer/quantile sweeps into interpretable functionals per variable and scale band. The paper's simulation settings are buffers 2, 4, 6 (with the LISA normalization radius `d_norm = 4`), quantile levels 0, 0.1, ..., 1, a fine scale band {2} and a broad scale band {6}. No response is used at any point of the generation. ```{r field} field <- gcf_field(sim_grid, coords = c("x", "y"), vars = c("x1", "x2", "x3"), buffers = c(2, 4, 6), probs = seq(0, 1, 0.1), d_norm = 4, fine_band = 2, broad_band = 6) summary(field) ``` The candidate field holds the 3 raw covariates (category `X`) plus their GCF variables: band-averaged pattern operators (category `P`) and functionals of the neighbourhood quantile curve (category `D` -- median, IQR, low tail, high tail, skew per scale band). ```{r meta} head(field$meta, 10) ``` ## Step 3b: select stable GCF variables `gcf_select()` screens the candidates with the rf_imp kernel: on each of `B` spatial-block subsamples it fits a random forest and keeps the top derived variables by impurity importance; a (variable x category) group qualifies when it fires in at least `pi_thr` of the subsamples, and contributes its most frequently kept member. The raw covariates are always kept. The spatial blocks come from `gcf_blocks()`; the paper's simulation uses blocks of side 6 (twice the residual variogram range, as in the case study). ```{r select} blocks <- gcf_blocks(sim_grid[, c("x", "y")], size = 6) sel <- gcf_select(field, y = sim_grid$y1, blocks = blocks, B = 10, seed = 1) sel ``` The selection frequencies show how consistently each derived variable is kept across the spatial subsamples: ```{r freq} round(head(sel$freq, 8), 2) ``` ## Does GCF improve prediction? A cross-validated comparison We now compare two feature sets with the paper's adopted learner, a random forest from the external **randomForest** package: - **base**: the raw covariates `x1`, `x2`, `x3`; - **GCF**: the selected GCF variable set. Following the paper, the comparison uses five-fold cross-validation under two partitions: ordinary random folds, and spatial-block folds that assign whole blocks to folds (shuffled round-robin), which tests spatial transferability. To stay leakage-free, the GCF variable selection is re-run inside each training fold. ```{r folds} folds_random <- function(n, K = 5, seed = 1) { set.seed(seed) sample(rep_len(seq_len(K), n)) } folds_blockwise <- function(block_id, K = 5, seed = 36) { set.seed(seed) ub <- sample(unique(block_id)) as.integer(stats::setNames(rep_len(seq_len(K), length(ub)), ub)[block_id]) } fold_rd <- folds_random(nrow(sim_grid)) fold_sp <- folds_blockwise(blocks) ``` The comparison needs the **randomForest** package (a suggested, not required, dependency of gcf); the two chunks below are evaluated only when it is installed. ```{r cv, eval = has_rf} library(randomForest) y <- sim_grid$y1 X <- field$candidates x_cols <- field$meta$feature[field$meta$category == "X"] cv_rf <- function(fold, cols_by_fold, ntree = 300) { ks <- sort(unique(fold)) r2 <- rmse <- numeric(length(ks)) for (i in seq_along(ks)) { te <- which(fold == ks[i]); tr <- which(fold != ks[i]) cols <- cols_by_fold[[i]] set.seed(1) fit <- randomForest(X[tr, cols, drop = FALSE], y[tr], ntree = ntree) p <- as.numeric(predict(fit, X[te, cols, drop = FALSE])) r2[i] <- 1 - sum((y[te] - p)^2) / sum((y[te] - mean(y[te]))^2) rmse[i] <- sqrt(mean((y[te] - p)^2)) } c(R2 = mean(r2), RMSE = mean(rmse)) } # per-fold GCF selection on the training part of each fold (leakage-free) select_by_fold <- function(fold) { lapply(sort(unique(fold)), function(k) { gcf_select(field, y = y, blocks = blocks, train = which(fold != k), B = 10, seed = 1)$selected }) } results <- do.call(rbind, lapply( list(random = fold_rd, spatial = fold_sp), function(fold) { S <- select_by_fold(fold) base <- cv_rf(fold, rep(list(x_cols), 5)) gcfv <- cv_rf(fold, S) data.frame(feature_set = c("base", "GCF"), R2 = c(base["R2"], gcfv["R2"]), RMSE = c(base["RMSE"], gcfv["RMSE"])) })) results$partition <- rep(c("random", "spatial"), each = 2) rownames(results) <- NULL results[, c("partition", "feature_set", "R2", "RMSE")] ``` ```{r improvement, eval = has_rf} imp <- do.call(rbind, lapply(split(results, results$partition), function(d) { data.frame(partition = d$partition[1], R2_gain_pct = 100 * (d$R2[2] - d$R2[1]) / d$R2[1], RMSE_drop_pct = 100 * (d$RMSE[1] - d$RMSE[2]) / d$RMSE[1]) })) round(imp[, -1], 1) ``` The GCF variables raise the cross-validated R-squared and lower the RMSE under both partitions, with the larger gain under the spatial-block partition -- the setting that matters for predicting into unsampled areas. This mirrors the paper's simulation result (its Table 2, computed with B = 80, seven learners, and 500-tree forests). ## The case study data The package also ships the paper's case study, `bio_grid`: vascular plant species richness over the Southwest Australian Floristic Region on a 10-km grid (6229 cells, 958 of them observed) with twelve environmental covariates. The paper generates the GCF variables on the projected kilometre coordinates with buffers 20--100 km, 21 quantile levels, scale bands {20, 30} and {90, 100} km, and selects with blocks of side 132 km: ```{r case, eval = FALSE} data(bio_grid) obs <- bio_grid[bio_grid$observed, ] covs <- c("Elevation", "Slope", "Precipitation", "Radiation", "DistWater", "DistBuilt", "SoilN", "SoilC", "SoilClay", "SoilDepth", "SoilpH", "SoilBD") field <- gcf_field(obs, coords = c("xkm", "ykm"), vars = covs, buffers = seq(20, 100, 10), probs = seq(0, 1, 0.05), d_norm = 100, fine_band = c(20, 30), broad_band = c(90, 100)) blocks <- gcf_blocks(obs[, c("xkm", "ykm")], size = 132) sel <- gcf_select(field, y = obs$richness, blocks = blocks, B = 80) ``` (Not run here: the full case-study selection takes several minutes.) ## Reference Song, Y. (2026). Generalized covariate field (GCF): spatial-pattern and neighbourhood-distribution feature expansion improves geospatial prediction. *International Journal of Geographical Information Science*, 40, 1--29. The package source is available from GitHub: .