The GCF workflow: from spatial variables to better predictions

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.

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.

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.

data(sim_grid)
head(sim_grid)
#>            y1         x1           x2          x3 x y
#> 1  0.68156075 -0.7437869 -0.390083041  0.01132319 1 1
#> 2  0.01164551  0.6994206  0.005002826 -0.70919621 2 1
#> 3 -0.26875636 -0.7354008  0.308990469  0.77466509 3 1
#> 4  0.28741097  1.2951746  0.366099153 -0.99436781 4 1
#> 5  1.51411986  1.2656499  2.273028010  0.94935635 5 1
#> 6  2.22486305 -0.2625940  1.066057112 -0.10049406 6 1

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.

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)
#> Generalized covariate field (GCF)
#>   locations:  900
#>   variables:  3 (x1, x2, x3)
#>   candidates: 93  [X (raw) 3 | P (pattern) 60 | D (context) 30]
#> 
#> Candidate variables per input variable and category:
#>     
#>       X  P  D
#>   x1  1 20 10
#>   x2  1 20 10
#>   x3  1 20 10
#> 
#> Reduction: d_mode = functional | fine band {2} | broad band {6}
#> Full sweeps kept in $psi (87 pattern features) and $zx (99 context-quantile features).

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).

head(field$meta, 10)
#>              feature group category operator scale
#> 1                 x1    x1        X     <NA>  <NA>
#> 2      x1_D_fine_med    x1        D     <NA>  fine
#> 3      x1_D_fine_iqr    x1        D     <NA>  fine
#> 4   x1_D_fine_lotail    x1        D     <NA>  fine
#> 5   x1_D_fine_hitail    x1        D     <NA>  fine
#> 6     x1_D_fine_skew    x1        D     <NA>  fine
#> 7     x1_D_broad_med    x1        D     <NA> broad
#> 8     x1_D_broad_iqr    x1        D     <NA> broad
#> 9  x1_D_broad_lotail    x1        D     <NA> broad
#> 10 x1_D_broad_hitail    x1        D     <NA> broad

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).

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
#> GCF variable selection (rf_imp + spatial-block stability + group voting)
#>   candidates: 93 | resamples B = 10 | pi threshold = 0.6
#>   selected:   8 (3 forced raw + 5 derived)
#>   forced:  x1, x2, x3
#>   derived: x1_D_fine_med, x2_D_fine_med, x2_P_gc, x3_D_fine_med, x3_P_gc

The selection frequencies show how consistently each derived variable is kept across the spatial subsamples:

round(head(sel$freq, 8), 2)
#>    x1_D_fine_med x1_D_fine_lotail    x2_D_fine_med x2_D_fine_lotail 
#>                1                1                1                1 
#> x2_D_fine_hitail          x2_P_gc    x3_D_fine_med x3_D_fine_lotail 
#>                1                1                1                1

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.

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.

library(randomForest)
#> randomForest 4.7-1.2
#> Type rfNews() to see new features/changes/bug fixes.

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")]
#>   partition feature_set        R2      RMSE
#> 1    random        base 0.4570685 0.7402535
#> 2    random         GCF 0.6918846 0.5573999
#> 3   spatial        base 0.4176858 0.7419110
#> 4   spatial         GCF 0.6568244 0.5606579

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)
#>         R2_gain_pct RMSE_drop_pct
#> random         51.4          24.7
#> spatial        57.3          24.4

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:

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.)