Reproducing the Generalized Covariate Field

1 Method Overview & Reproduction Scope

1.1 Core Idea

Most attempts to improve spatial prediction add model: a deeper forest, a cleverer kernel, another layer. GCF starts from the opposite observation — the covariates themselves carry spatial information that point values throw away. A location is not just its covariate vector: it sits inside a neighbourhood of covariate values, with local dependence, heterogeneity, multiscale variation and outlyingness of its own. The generalized covariate field (GCF) turns that surrounding structure into explanatory variables.

GCF is a feature transform, not a learner — there is no fitting inside it. Each covariate is expanded into two families of features at every location, sampled and unsampled: spatial pattern features ψ (eleven local operators over a series of buffer radii) and neighbourhood distribution features Z_x (quantiles of the covariate over multiscale neighbourhoods). The sweeps are collinear by construction, so a functional reduction compresses them into interpretable candidates, and a stability selection keeps the few that matter. Any learner then fits on raw covariates + selected features exactly as it would have fit on raw covariates alone.

• Research problem: prediction accuracy is constrained by spatial heterogeneity and sparse, uneven sampling; existing methods respond by increasing algorithmic complexity while underusing the spatial information embedded in the covariates themselves.
• One-line contribution: expand the covariate space — pattern operators + neighbourhood quantiles + leakage-free reduction + spatial-block stability selection — so that a fixed learner sees the spatial structure it could not extract from point values.
• Why it matters here: in the simulation, GCF lifts the spatial-block CV R² of every one of seven learners — RF from 0.419 to 0.657 (+56.5 %) — and the reproduction of the paper's Table 2 below is byte-identical. In the case study it lifts RF spatial R² from 0.296 to 0.346 (+16.7 %) with a 3.0 % RMSE reduction.

1.2 Method Logic

GCF runs in four steps, followed by the paper's validation design. The notation follows the paper's Section 2; every step below is exercised by the pipeline in §3.

Step 1 — spatial pattern features ψ

For each covariate and each buffer radius b, eleven local operators describe the covariate's spatial behaviour around location v: LISA on the local z-score (spatial dependence), local Geary's c (local contrast), log local variance (heterogeneity), rank-binned quantile entropy (local diversity), geocomplexity (single scale), log scale-variance across the scale series (single scale), the local variogram exponent (roughness of local spatial structure), and positive/negative z-score and MAD outlier strengths (local outlyingness). The simulation uses buffers {2, 4, 6} on the unit grid → 87 ψ columns from 3 covariates; the case study uses {20, …, 100} km → 12 × 83 = 996.

Step 2 — neighbourhood distribution features Z_x

For each location v, covariate x and buffer b, summarise the covariate's distribution over the support points inside the buffer by quantile levels τ:

Sim: 3 covariates × 3 buffers × 11 quantile levels = 99 columns. Case: 12 × 9 × 21 = 2,268. Together with ψ this is what "the field" in the method's name means: every covariate becomes a field of spatial descriptors around every location, computable at unsampled locations too.

Step 3a — functional reduction (leakage-free)

The quantile/buffer sweeps are strongly collinear, so they are reduced to interpretable functionals — a per-row transform computed once, with no response involved. Two scale bands are fixed a priori: fine ({20, 30} km in the case, {2} in the sim) and broad ({90, 100} km / {6}). For each covariate × band, the band-averaged quantile curve yields five D (context) functionals:

and each buffered ψ operator is band-averaged into P (pattern) candidates (the two single-scale operators pass through). With the raw covariates X forced in, the reduced candidate field is 93 columns in the sim (X 3, P 60, D 30) and 372 in the case (X 12, P 240, D 120) — the paper's Table 4, which step 4 rebuilds byte-for-byte.

Step 3b — selection: rf_imp kernel + spatial-block stability + group voting

Selection runs on training rows only (inside each CV fold, so it is leakage-safe). Three moves:

1. Per-subsample selector. On a subsample, fit a random forest (ranger, 200 trees) and keep the top K = 20 derived features by impurity importance: s_b,j = 1{rank(−imp_j) ≤ 20}. The raw covariates are never candidates for removal.
2. Spatial-block stability. The training rows fall in spatial blocks; for b = 1…B (B = 80) draw ⌊0.7 M⌋ of the M blocks, run the selector, and record the selection frequency freq_j = (1/B) Σ_b s_b,j.
3. Medium group voting. Group the derived features by (covariate × category). A group qualifies if any member is kept in at least π = 0.6 of the resamples — fire_g = (1/B) Σ_b 1{Σ_j∈Gg s_b,j > 0} ≥ π — and contributes its highest-frequency member. Voting de-dilutes collinear blocks: sibling features split freq so none clears π alone, while the block as a whole fires reliably.

The output S = covariates ∪ group representatives is typically ~20–25 features in the case study and 8 in the simulation.

The validation design — dual five-fold cross-validation

Two partitions of the same samples: ordinary random five-fold (interpolation), and spatial-block five-fold (transfer): square blocks of side L = 2 × the residual variogram range (132 km in the case; 6 in the sim) are assigned whole to folds by shuffled round-robin (seed 36). Metrics are per-fold mean ± sd:

Seven learners are compared, each with the paper's fixed settings — RF (randomForest, 500 trees), Cubist, GLMNET (ridge), LM, XGBoost, SVM (radial), KNN — first on the raw covariates, then on the per-fold re-selected GCF sets. Because both arms share folds and learners, any difference is attributable to the representation, not to the algorithm.

1.3 Reproduction Scope

Everything tagged Demo run regenerates from this folder with Rscript run-all.R on R 4.6.0 with gcf 0.1.0. The scope follows a deliberate split — simulation: full re-run; case study: lightweight rebuild — and three things are worth stating plainly before any number below is quoted.

How to use this page

• To reproduce the demo: run Rscript run-all.R test, then Rscript run-all.R. Every re-run number and figure below regenerates in about 140 seconds.
• To use GCF on your own data: §2.3 gives the data contract — one table with projected coordinates, a response at the sampled rows and covariates everywhere — and §4.1 the porting steps.
• To check rather than trust: §3.4 lists the 28 verification checks, including the byte-identity of every rebuilt table against data/paper-tables/.

2 Setup, Structure & Data Contract

2.1 Environment & Dependencies

The method is one package written for the paper and validated bit-identical against its analysis engine; the learners are the seven the paper compares.

# The method, and both input data sets (CRAN submission in preparation;
# until acceptance, install gcf from the released tarball)
install.packages("gcf")
# The seven learners of the per-model comparison
install.packages(c("ranger", "randomForest", "glmnet", "Cubist",
                   "xgboost", "e1071", "kknn"))

2.2 Project Structure & Configuration

One entry script runs the six numbered steps, the table builder, or the tests. Each step writes to results/, tables/, figs/ or data/derived/ and reads what the previous step wrote, so any step can be re-run in isolation after a configuration change.

GCF/
├── gcf.html                       # this guide
├── run-all.R                      # entry point: full run, single steps, tests
├── config/project-config.R        # the ONLY file to edit for a new domain
├── R/
│   ├── 00-config.R, 01-helpers.R  # paths, IO, folds, the 7 learners, base-R maps
│   ├── 10-sim-field.R             # sim: gcf_field -> 93 candidates + full-data selection
│   ├── 20-sim-cv.R                # sim: dual five-fold design + 10 per-fold selections
│   ├── 30-sim-permodel.R          # sim: 7 learners x base/GCF  ->  Table 2 (byte-check)
│   ├── 40-case-field.R            # case: gcf_field on the 958 observed cells -> Table 4
│   ├── 50-case-tables.R           # case: Tables 5-7 from case_final.rds (byte-check)
│   ├── 60-case-prediction.R       # case: both prediction surfaces, mapped
│   └── 90-tables.R                # LaTeX tables + session info
├── data/                          # sim-grid.csv + bio-grid.csv (written by steps 1 and 4),
│                                  #   case_final.rds, paper-tables/, provenance.md
├── tests/                         # test-01 method properties; test-02 reproduction
├── env/                           # requirements.md, runtimes.csv, session-info.txt
└── results/  tables/  figs/       # generated output, never hand-edited

Main configuration file

config/project-config.R is the single file a reader edits. Every value in it is the paper's setting; nothing else in the pipeline hard-codes a column name, a buffer or a seed.

# -- Simulation --------------------------------------------------------------
SIM_VARS    <- c("x1", "x2", "x3");  SIM_Y <- "y1";  SIM_COORDS <- c("x", "y")
SIM_BUFFERS <- c(2, 4, 6)          # buffer radii, grid units
SIM_PROBS   <- seq(0, 1, 0.1)      # 11 quantile levels for the Zx sweep
SIM_D_NORM  <- 4                   # normalisation distance (middle buffer)
SIM_FINE    <- 2;  SIM_BROAD <- 6  # reduction bands {2} / {6}
SIM_L       <- 6                   # spatial-block side

# -- Case study ---------------------------------------------------------------
CASE_BUFFERS <- seq(20, 100, 10)   # 9 buffer radii, km
CASE_PROBS   <- seq(0, 1, 0.05)    # 21 quantile levels
CASE_D_NORM  <- 100                # km
CASE_FINE    <- c(20, 30);  CASE_BROAD <- c(90, 100)
CASE_L       <- 132                # block side, km (2 x variogram range 66)

# -- Selection (rf_imp + spatial-block stability + group voting) --------------
SEL_B    <- 80L;   SEL_FRAC  <- 0.7;  SEL_PI   <- 0.6
SEL_KTOP <- 20L;   SEL_TREES <- 200L; SEL_SEED <- 1L

# -- Cross-validation ---------------------------------------------------------
K_FOLDS <- 5L;  SEED_SHUF <- 36L;  SEED_RANDOM <- 1L

# -- Learners -----------------------------------------------------------------
LEARNERS <- c("RF", "Cubist", "GLMNET", "LM", "XGBoost", "SVM", "KNN")
RF_TREES <- 500L                   # trees of the RF learner (not the selector)

# -- Compute ------------------------------------------------------------------
CORES <- 5L   # wall clock only: every selection re-seeds its own RNG

2.3 Data Contract

GCF needs one plain table: projected coordinates, a response at the sampled rows, and covariates at every row where you want features. Unlike a purely tabular method, the coordinates are load-bearing — buffers, blocks and folds are all defined on them.

The shipped result object

3 Reproduction Pipeline, Results & Validation

3.1 Pipeline Execution

cd GCF                     # this folder
Rscript run-all.R test     # 1. 28 checks on the method and on the committed results
Rscript run-all.R          # 2. full pipeline

Rscript run-all.R 30 50    # re-run single steps (10 20 30 40 50 60 90)

Project : gcf-reproduce
Domain  : sim: 30 x 30 synthetic grid | case: SWAFR plant species richness
== Step 10  Simulation: GCF candidate field ==============
   sim_grid: 900 cells, response y1, covariates x1, x2, x3
   gcf_field (psi + Zx + reduction) took 5.3 s
   candidates: 93  [X (raw) 3 | P (pattern) 60 | D (context) 30]
   gcf_select on all 900 cells (B = 80) took 12.3 s
   selected 8 variables: x1, x2, x3, x1_D_fine_med, x2_D_fine_med,
                         x2_P_gc, x3_D_fine_med, x3_P_gc
== Step 20  Simulation: CV folds + per-fold selection ====
   random folds: 180/180/180/180/180
   spatial folds (25 blocks): 180/180/180/180/180
   selected set sizes — spatial: 8/8/8/8/8 | random: 8/8/8/8/8
== Step 30  Simulation: per-model dual five-fold CV ======
   RF 5.8 s · Cubist 1.4 s · GLMNET 0.2 s · LM 0.0 s
   XGBoost 0.9 s · SVM 0.2 s · KNN 0.1 s
   best learner by spatial GCF RMSE: LM
   Table 2 vs paper: IDENTICAL byte-for-byte (max |diff| = 0)
== Step 40  Case study: GCF feature generation ===========
   bio_grid: 6229 cells, 958 observed (richness > 0)
   gcf_field on 958 cells (9 buffers, 21 quantiles) took 59.4 s
   candidates: 372  [X (raw) 12 | P (pattern) 240 | D (context) 120]
   Table 4 vs paper: IDENTICAL byte-for-byte
   demo selection: 21 variables (19 also in the paper's 21-variable set)
== Step 50  Case study: rebuild Tables 5-7 ===============
   table05-case-permodel vs paper: IDENTICAL byte-for-byte
   table06-case-featuresets vs paper: IDENTICAL byte-for-byte
   table07-case-shapley vs paper: IDENTICAL byte-for-byte
== Step 60  Case study: prediction surfaces ==============
   cross-section: 93 cells along lat -33.01 (ykm = -3654)
Finished in 132.3 s

3.2 Core Analytical Steps

Six steps: purpose → code → output → result → how to read it. Steps 1–3 are the simulation (fully re-run); steps 4–6 are the case study (feature demo + rebuild). Each maps one-to-one onto a numbered script.

Step 1 · R/10-sim-field.R

Simulation — build the GCF candidate field

What this step does

Loads the packaged 30 × 30 simulation grid (response y1 with smooth spatial structure; three covariates that are nearly unstructured at point level), writes a verbatim CSV copy, and runs the whole feature transform in one call: gcf_field() computes the eleven ψ operators over buffers {2, 4, 6}, the Z_x quantiles at 11 levels, and the functional reduction with bands {2}/{6} — 93 candidates from 3 covariates. It then runs one full-data B = 80 selection as a first look at what GCF keeps.

Code

R — from R/10-sim-field.R

field <- gcf::gcf_field(sim, coords = c("x", "y"), vars = c("x1", "x2", "x3"),
                        buffers = c(2, 4, 6), probs = seq(0, 1, 0.1),
                        d_norm = 4, include_vario_exp = TRUE,
                        vario_buffers = c(2, 4, 6),
                        fine_band = 2, broad_band = 6)
# 900 x 93 candidates: field$candidates, field$meta, field$psi, field$zx

bid <- gcf::gcf_blocks(sim[, c("x", "y")], 6)      # 25 blocks of 36 cells
sel <- gcf::gcf_select(field, y = sim$y1, blocks = bid, B = 80)

Example output

Selected variable	Category	Forced	Selection frequency
x1, x2, x3	X (raw)	yes	—
x1_D_fine_med	D (context)	no	1.000
x2_D_fine_med	D (context)	no	1.000
x2_P_gc	P (pattern)	no	1.000
x3_D_fine_med	D (context)	no	1.000
x3_P_gc	P (pattern)	no	1.000

results/sim-selected-features.csv — the full-data selection keeps 8 of 93 candidates: the three covariates (forced) plus five derived features, every one kept in all 80 stability resamples. The fine-band context median of each covariate is selected — exactly the feature family panel (e) of fig02 shows recovering the response's spatial structure.

Demo run fig01 — the simulation design: the response (a) carries smooth spatial structure; the covariates (b–d) look close to noise at point level. Point values of x explain y only partly, and the missing part is spatial — the designed challenge GCF is built for.

Demo run fig02 — x1 traced through the method: raw (a), ψ LISA and log local variance at the finest buffer (b, c), the buffer-6 neighbourhood median (d), and two reduced candidates (e, f). Compare (d)–(e) with fig01(a): the response's structure emerges from the covariate's neighbourhood distribution.

How to read this result

• Feature generation is deterministic arithmetic. No RNG anywhere in ψ/Z_x/reduction (asserted in test-01), which is half of why this page can promise byte-identical tables.
• The reduction is the honesty step. 3 covariates → 87 ψ + 99 Z_x = 186 raw features would invite overfitting and collinearity; band-averaging to 93 interpretable candidates happens before any response is seen, so nothing about it can leak.
• The selection is unanimous here. Frequency 1.0 on all five derived features means every one of 80 block-subsamples ranked them top-20 — the simulation's signal is strong and stable. The case study (step 4) shows the marginal-group behaviour instead.

Step 2 · R/20-sim-cv.R

Simulation — the dual five-fold design and per-fold selection

What this step does

Builds both partitions of the 900 cells — random five-fold (seed 1) and spatial-block five-fold (25 blocks of side 6, whole blocks round-robined to folds after a seed-36 shuffle) — then re-runs the B = 80 selection inside every training fold of both partitions: ten selections in all. The GCF arm of step 3 uses these per-fold sets, so no validation cell ever influences the features it is predicted with.

Code

R — from R/20-sim-cv.R and R/01-helpers.R

fold_sp <- gcf_folds_rr(bid, K = 5L, seed = 36L)   # whole blocks -> folds
fold_rd <- gcf_folds_random(nrow(sim), K = 5L, seed = 1L)

# leakage-safe: one selection per training fold (x2 partitions = 10)
S_sp <- lapply(1:5, function(k)
  gcf::gcf_select(field, y = sim$y1, blocks = bid,
                  train = which(fold_sp != k), B = 80)$selected)

Example output

Demo run fig03 — the two partitions. Random assignment (a) measures interpolation: every validation cell has training neighbours. Spatial-block assignment (b) measures transfer: whole 6 × 6 blocks are held out, so a validation cell's entire neighbourhood is unseen. The paper reports both, and the gap between them is itself informative.

How to read this result

• All ten selections return 8 variables — the same three covariates plus five derived features per fold (results/sim-selected-by-fold.csv). Stability under refolding is the selection design doing its job.
• Per-fold selection is the leakage discipline. Selecting once on all data and then cross-validating would let validation rows vote on the features. The cost is honest: ten B = 80 selections are ~80 % of this step's 21 seconds.
• Blocks respect geography, folds respect blocks. test-01 asserts every block lands entirely inside one fold — the property that makes spatial-block CV a transfer test rather than a leaky variant of random CV.

Step 3 · R/30-sim-permodel.R

Simulation — seven learners, base vs GCF (paper Table 2)

What this step does

Runs the paper's per-model comparison: every learner fits twice per fold — once on the three raw covariates, once on that fold's re-selected GCF set — under both partitions. Learner settings are ported verbatim from the paper's engine, and every fit starts from set.seed(1), so each fit is a pure function of its inputs. The finished wide table is then compared byte-for-byte against the published Table 2 shipped in data/paper-tables/.

Code

R — from R/30-sim-permodel.R

X <- as.matrix(field$candidates)     # the paper's engine used a matrix
x_cols <- field$meta$feature[field$meta$category == "X"]

for (ln in LEARNERS) for (pt in c("random", "spatial")) {
  fold <- fold_of[[pt]]; Sb <- S_of[[pt]]
  b <- gcf_cv(X, y, fold, function(k) x_cols, ln)                 # baseline
  g <- gcf_cv(X, y, fold, function(k) Sb[[as.character(k)]], ln)  # GCF
  # per-fold mean/sd of R2 and RMSE + relative improvements -> p1
}
t2 <- gcf_permodel_wide(p1, best)     # the paper's Table-2 shape
compare_with_paper(F_TABLE02, "data/paper-tables/table02-sim-permodel.csv")

Result — Table 2, reproduced byte-for-byte

Learner	Best	Random CV — R²			Random CV — RMSE			Spatial CV — R²			Spatial CV — RMSE
		base	GCF	imp.%	base	GCF	imp.%	base	GCF	imp.%	base	GCF	imp.%
LM	*	0.514	0.700	36.0	0.700	0.550	21.4	0.469	0.668	42.4	0.706	0.551	22.0
GLMNET		0.514	0.700	36.2	0.701	0.550	21.4	0.470	0.668	42.1	0.707	0.552	22.0
RF		0.457	0.697	52.6	0.740	0.552	25.4	0.419	0.657	56.5	0.741	0.561	24.3
Cubist		0.476	0.678	42.3	0.727	0.570	21.6	0.438	0.637	45.6	0.728	0.575	21.0
KNN		0.467	0.684	46.5	0.734	0.565	23.0	0.429	0.637	48.5	0.736	0.578	21.5
XGBoost		0.445	0.677	52.1	0.748	0.571	23.7	0.399	0.619	55.2	0.754	0.590	21.8
SVM		0.472	0.645	36.7	0.730	0.598	18.0	0.433	0.591	36.6	0.734	0.624	15.0

results/table02-sim-permodel.csv — byte-identical to the published Table 2 (results/sim-permodel-check.csv lists all 84 numeric cells side by side; every difference is 0). Learners sorted by spatial GCF RMSE; the asterisk marks the adopted learner. Improvements are relative percentages.

Demo run fig04 — Table 2 as a picture: RMSE of the covariate baseline (grey) and GCF (green) per learner, spatial-block (a) and random (b) CV. Every learner moves left; the linear models move furthest left of all.

Why the match is exact, not approximate

Four links in a chain: (1) the gcf package is a verbatim port of the paper's engine, validated to max|diff| = 0 on ψ, Z_x, candidates and selections; (2) folds use the paper's seeds (36 / 1); (3) every selection re-seeds seed = 1; (4) every learner fit starts from set.seed(1) with the paper's exact settings. Same inputs, same code, same package versions — same bytes. If a future learner-package release moves a third decimal, results/sim-permodel-check.csv will name the exact cells.

How to read this result

• All seven learners improve, under both partitions. That pattern — not any single number — is the paper's claim: the gain is attributable to the representation, not to a lucky algorithm.
• GCF-based LM wins. Spatial RMSE 0.551 against 0.706 for baseline LM. The features encode the spatial structure before fitting, so a simple linear map on top suffices — and the baseline's own ordering (LM 0.470 spatial R² vs RF 0.419) shows there was little residual nonlinearity for trees to exploit.
• The spatial gain exceeds the random gain (RF +56.5 % vs +52.6 % R²): neighbourhood features transfer to unseen blocks better than point values do — the property that matters for mapping.

Step 4 · R/40-case-field.R

Case study — real feature generation on the 958 observed cells

What this step does

Runs the identical gcf_field() call the paper uses — 9 buffers from 20 to 100 km, 21 quantile levels, bands {20,30}/{90,100} km, d_norm = 100 — on the 958 observed cells of the SWAFR grid (support = the observed cells; the paper's production run uses the full 6,229-cell grid as support). It rebuilds the paper's Table 4 from the resulting field and then runs one B = 80 selection with 132 km blocks, comparing the outcome with the paper's full-support selection.

Code

R — from R/40-case-field.R

obs <- bio_grid[bio_grid$observed, ]           # 958 sampled cells
field <- gcf::gcf_field(obs, coords = c("xkm", "ykm"), vars = CASE_VARS,
                        buffers = seq(20, 100, 10), probs = seq(0, 1, 0.05),
                        d_norm = 100, include_vario_exp = TRUE,
                        vario_buffers = seq(20, 100, 10),
                        fine_band = c(20, 30), broad_band = c(90, 100))
# 958 x 372 candidates; psi 958 x 996; Zx 958 x 2268

bid <- gcf::gcf_blocks(obs[, c("xkm", "ykm")], 132)   # 38 occupied blocks
sel <- gcf::gcf_select(field, y = obs$richness, blocks = bid, B = 80)

Result — Table 4, rebuilt byte-for-byte

Category	Symbol	Raw count	Reduced count	Construction
Covariates (x)	x	12	12	raw covariates (forced)
Spatial pattern (ψ)	psi	996	240	11 operators, band-averaged
Neighbourhood distribution (Zx)	Zx	2,268	120	5 quantile functionals × 2 scale bands
Total		3,276	372

results/table04-case-predictor-categories.csv — byte-identical to the published Table 4. The counts depend only on the settings (12 covariates × 83 ψ columns; 12 × 9 buffers × 21 quantiles; 12 × (20 P + 10 D) + 12), so they rebuild exactly even on the observed-cell support.

The demo selection against the paper's full-support selection

Agreement	Variables
Selected in both (19)	the 12 covariates (forced), plus Slope_P_lvar_broad, Precipitation_D_broad_hitail, DistBuilt_D_broad_med, SoilC_D_broad_iqr, SoilClay_D_broad_lotail, SoilDepth_D_broad_hitail, SoilpH_D_broad_med
Demo only (2)	Slope_D_broad_skew, Precipitation_P_lvar_fine
Paper only (2)	Precipitation_P_lvar_broad, DistBuilt_P_poutlier_z_broad

results/case-selected-observed-demo.csv — 21 variables selected on the observed-cell field, 19 shared with the paper's 21-variable full-support set. Note what the two disagreements are: the Precipitation local-variance group fires in both runs but contributes a different band's representative, and two genuinely marginal groups swap. The stable core — the seven broad-band context features — survives the support change untouched.

Demo run fig05 — the demo field over the observed cells (10 km cells drawn at 14 km for visibility; axes are EPSG:3577 Albers km): observed richness (a), the strongest covariate precipitation (b), its ψ log local variance at the 20 km buffer (c), the fine-band pattern candidate built from it (d), and two context candidates the paper's selection also keeps — the fine-band median of distance to built areas (e) and the fine-band 10th quantile of soil clay (f).

The support caveat, quantified

Feature values here are not the paper's: Z_x at an observed cell summarises only the observed neighbours, not the full grid. That is why this step's selection is a demonstration, not a reproduction — and the 19/21 overlap is a fair measure of how much the support choice moves the outcome on these data. For production use, build the field on every location you will ever predict (§2.3, rule 1); the numbers the paper reports all come from the full-support field inside case_final.rds.

Step 6 · R/60-case-prediction.R

Case study — the two prediction surfaces

What this step does

Maps the paper's Part-4 result: RF trained on all 958 observed cells — once on the covariates, once on the 21-variable full-support GCF set — predicting all 6,229 grid cells. Both finished surfaces live in case_final$prediction; this step draws them, the difference, and the paper's west–east cross-section, without refitting.

Example output

Demo run fig06 — predicted species richness from GCF-RF (a) and covariate-RF (b) on a shared scale; their difference (c; blue = GCF higher, scale clipped at the 98th percentile of |difference|); and the cross-section at y = −3654 km (d), the paper's transect. Axes are EPSG:3577 Albers km (the paper's own maps use a Web-Mercator display projection; the values are identical).

Surface	Mean	sd	Min	Max	Note
GCF prediction	33.211	15.362	5.164	90.002
covariate (RF) prediction	33.566	14.603	6.689	88.508	r(GCF, base) = 0.953
difference (GCF − base)	−0.355	4.649	−21.714	16.607	GCF higher in 44.4 % of cells

results/case-prediction-summary.csv. The two surfaces agree on the broad gradient (r = 0.953) and disagree exactly where neighbourhood information matters — the difference map is spatially organised, not noise.

How to read this result

• Follow the transect. Along y = −3654 km the GCF surface keeps the high peaks (up to ~80) and the deep troughs that the covariate-only surface dampens — richer features let the learner commit to local extremes instead of hedging towards the regional mean.
• The GCF surface has slightly larger spread (sd 15.36 vs 14.60) with the same mean level — consistent with preserved contrasts rather than added noise, given its lower cross-validated error in Table 5.
• The difference map is the sampling map's shadow. The largest disagreements sit in the sparsely observed east and interior — where a validation block has few nearby samples and the neighbourhood features carry the most extra information.

3.3 Reading the Results Together

3.4 Validation

Rscript run-all.R test runs two scripts and 28 checks in about 16 seconds. They test different things, and it is worth knowing which is which.

What test-01 pins down

• Feature generation is deterministic — gcf_field() twice under a deliberately hostile ambient seed returns identical candidates, and the X columns are the raw covariates unchanged.
• The D functionals behave like quantiles — q_0.10 ≤ median ≤ q_0.90 at every location and the interquartile range is non-negative; and the band defaults (min/max buffer) equal the explicit sim bands.
• The spatial design is exact — the 30 × 30 grid at L = 6 gives 25 blocks of 36 cells; round-robin folds keep whole blocks together, are deterministic for a fixed seed, and random folds are balanced to within one cell.
• Selection invariants — every raw covariate is kept, the result is reproducible for a fixed seed, and every selected variable exists in the candidate field.
• The CV helper is honest — a fold's RMSE from gcf_cv() equals a hand-fitted LM on the same split to 10⁻⁹.

What test-02 pins down

• The sim field rebuilds from scratch to the committed structure (93 candidates; ψ 87, Z_x 99) and matches the stored pipeline field exactly.
• The fold-1 spatial selection re-derives to the committed set.
• Table 2 is byte-identical to the published CSV, and Tables 5, 6, 7 rebuilt from case_final.rds are byte-identical too (checked as data frames to 10⁻¹² and as raw bytes).
• Table 4's counts (12 / 996 / 2,268 / 3,276 raw; 12 / 240 / 120 / 372 reduced) hold, and the shipped prediction grid is intact (6,229 cells, both surfaces finite).

4 Adaptation, Writing & Reproducibility

4.1 Port to Your Domain

GCF applies wherever you have a response measured at some locations and covariates measured over a grid or a dense point set with projected coordinates: soil properties, air quality, species richness, urban indicators, mineral prospectivity. The whole method is two calls:

library(gcf)
field <- gcf_field(mydata, coords = c("xkm", "ykm"),
                   vars = my_covariates,
                   buffers = my_buffers,          # e.g. 1-5 cell diameters
                   fine_band = head(my_buffers, 2),
                   broad_band = tail(my_buffers, 2))
blocks <- gcf_blocks(mydata[, c("xkm", "ykm")], size = 2 * variogram_range)
sel <- gcf_select(field, y = mydata$y[train], blocks = blocks, train = train)
X <- as.matrix(field$candidates)[, sel$selected]   # -> any learner

What to edit

1. config/project-config.R — coordinates, response, covariates (CASE_* block), then the spatial settings below. Nothing else.
2. *_BUFFERS — a series from roughly one cell diameter up to at or beyond the response's variogram range. Integer values (§2.3, rule 3).
3. *_FINE / *_BROAD — the two ends of your buffer series; the defaults (min / max buffer) are a sensible start.
4. *_L — spatial-block side ≈ 2 × the residual variogram range, so validation blocks sit beyond short-range dependence. Fit a variogram to the covariate-detrended response to get the range, as the paper does.
5. Selection settings (SEL_*) — the paper's B = 80, 0.7, π = 0.6, K = 20 travel well; lower B while iterating, restore it for the run you report.
6. Nothing else. Rscript run-all.R test, then Rscript run-all.R.

Choosing what to expand

• Every covariate must exist at the prediction locations — features are computed there too. Remote-sensing, terrain, climate and distance-to-feature layers are ideal; sample-only lab measurements cannot enter.
• Give the field the full support once. Build features on observed + prediction locations together, then subset rows for training. Features built per subset are different variables (step 4 measured the drift), and chunked builds truncate neighbourhoods at chunk boundaries.
• Let the selection do the pruning. The candidate field is deliberately generous (31 candidates per covariate); the stability selection with group voting is the part designed to cope with that. Pre-screening candidates by hand mostly re-implements it, worse.
• Watch scale mismatches. If your buffers are all far below the process scale, the Z_x features approach the point values and add little; if far above, they approach global constants. The buffer series should straddle the variogram range.

4.2 Write the Paper

The published study has a clean architecture worth copying: define the feature families, prove the mechanism on a simulation, then validate on the real case with a per-model table, a feature-set decomposition, and prediction surfaces.

Manuscript skeleton

1  Introduction        — spatial prediction constrained by heterogeneity and
                          sparse sampling; methods grow model complexity while
                          underusing the covariates' own spatial structure (aims)
2  Generalized covariate field
                        — 2.1 spatial pattern features (11 operators x buffers)
                          2.2 neighbourhood distribution features (quantiles)
                          2.3 functional reduction (bands; D and P candidates)
                          2.4 selection: rf_imp + block stability + group voting
                          2.5 validation: dual five-fold CV, L = 2 x range
3  Simulation          — design; per-model table; why linear learners win
4  Case study          — data; the field (Table-4 shape); selection;
                          per-model table; feature sets + Shapley; prediction
5  Discussion          — representation vs model complexity; transfer gains;
                          support and scale choices; limits
6  Conclusions         — 6-8 sentences answering the aims

4.3 Final Reproducibility Package

What is in the download, and what is not

References & credits

Demo outputs were generated by the scripts in this folder using the gcf package and its bundled data; the case-study validation tables and surfaces come from the paper's computed result object, shipped here unchanged.

1. 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. doi:10.1080/13658816.2026.2729719 · PDF (source of the method and of every published number cross-checked on this page)
2. Song Y. gcf: Generalized Covariate Field. R package, version 0.1.0 CRAN submission in preparation (the reference implementation — a verbatim, validated port of the paper's analysis engine — and the source of both data sets)
3. GCF code, data and examples — github.com/yongzesong/gcf (the standalone R implementation gcf.R, the simulation data and a worked example, as published with the paper)
4. Mokany K, McCarthy JK, Falster DS, Gallagher RV, Harwood TD, Kooyman R, Westoby M (2022). Patterns and drivers of plant diversity across Australia. Ecography 2022(11):e06426. doi:10.1111/ecog.06426 (the vascular-plant richness observations behind the case study)
5. Online GCF calculator — yongzesong.com/app/gcf/ (the same gcf feature engine compiled to WebAssembly, running entirely in the browser)