Reproducing the Singularity Regression Kriging (SRK) Model

1 Method Overview & Reproduction Scope

1.1 Core Idea

Regression kriging splits a spatial variable into a deterministic trend and a spatially correlated residual, predicts the trend from covariates, and kriges the rest. The weak link is the trend model. Ordinary kriging (OK) assumes there is no trend at all — a constant mean — and regression kriging usually assumes a linear one. Both leave nonlinear, multiscale structure sitting in the residual, where a variogram fitted under an assumption of stationarity cannot represent it.

Ren et al. (2026) change exactly one thing: what the trend model is allowed to see. Alongside each environmental covariate X_k they compute its singularity index α_k(s) — the exponent describing how the local intensity of that covariate scales as the measuring window grows. A value of 2 means the field is locally uniform; below 2 means local enrichment; above 2 means depletion. The random forest is then trained on the covariates and their singularity indices, and ordinary kriging is applied to what the forest still cannot explain.

• Research problem: in heterogeneous terrain the response is nonlinear, multiscale and non-Gaussian, so a kriging model with a simple trend degrades exactly where the interesting anomalies are.
• One-line contribution: SRK builds anomaly descriptors from the covariates at a ladder of spatial scales, hands them to the trend model, and leaves the residual closer to stationary — which is what the kriging step assumes in the first place.
• Why the covariate, not the response, matters: a singularity index built from the response could only be evaluated where the response is known — that is, never at a prediction location, and never inside a held-out block. Building it from covariates is what lets the feature survive spatial block cross-validation instead of leaking through it.
• What it buys: on the paper's simulation, reproduced below, SRK lifts R² from 0.876 to 0.972 on a near-normal field, and from 0.736 to 0.940 on a long-tailed one — the margin widens as the distribution departs from Gaussian. On the real cobalt data it lifts R² from 0.15 (ordinary kriging) to 0.36 under spatial blocking, and stays ahead of the strongest machine-learning benchmark on every one of 20 random-forest seeds (§3.2 Step 5).

Published reference · SRK Fig. 4

Published reference · SRK Fig. 3

1.2 Method Logic (Eq. 1–13)

SRK chains three standard tools in a specific order. Each link is simple; the model is what the middle one is given to work with.

Key concepts and equations

• Local covariate intensity (Eq. 2). At location s and scale r, take the mean absolute covariate value over every reference point inside the square window A(s, r):
  Ĉ_k(A(s, r)) = (1 / |A(s, r)|) Σ_{s′ ∈ A(s,r)} |X_k(s′)|(2)
  A scale carrying fewer than three reference points is discarded, not guessed at.
• The scaling law and the index (Eq. 1, 3–4). Across scales r₁ < … < r_J the theoretical relation C ∝ r^α−2 becomes a straight line in logarithms, and α is its slope plus two:
  log Ĉ_k(A(s,r)) = (α_k(s) − 2) log r + c_k(s) + e_k,r  α_k(s) = β̂_1,k(s) + 2(3–4)
  α = 2 is a spatially uniform field, α < 2 is local enrichment, α > 2 is local depletion. Locations with fewer than two usable scales are set to the neutral 2.
• The augmented feature set (Eq. 5). The q singularity columns that survive the SD filter join the p original covariates:
  F(s_i) = { X₁(s_i), …, X_p(s_i), α_k1(s_i), …, α_kq(s_i) }(5)
• Trend and residual (Eq. 6–7). A 500-tree random forest g(·) fits the trend, and the residual is what it misses:
  µ̂(s_i) = g(F(s_i))  ε(s_i) = Z(s_i) − µ̂(s_i)(6–7)
• Residual kriging and the final prediction (Eq. 8–10). Ordinary kriging interpolates ε under the unbiasedness constraint, and the two parts are added:
  Ẑ(s₀) = µ̂(s₀) + Σ_i λ_i ε(s_i), Σ_i λ_i = 1  δ(s₀) = √σ_K²(s₀)(8–10)
  The kriging standard error δ is the model's uncertainty measure.
• Accuracy (Eq. 11–13). R², RMSE and MAE on the held-out validation data, never on the training fit.

1.3 Reproduction Scope

This pipeline reproduces the paper's simulation experiment and its case study up to spatial block cross-validation. The article maps two trace elements; this tutorial follows cobalt end to end, because every step is written as a loop over the ELEMENTS list and one element is enough to show each step working — putting zinc back is a single config entry (§2.3). Everything tagged Demo run regenerates from the folder behind the download button; everything tagged Published reference is a static image from the article. The published prediction maps, cross-sections and uncertainty maps need the 97,234-cell 500 m prediction grid, which is out of scope here; §3.3 shows those figures from the article and says what it would take to regenerate them.

1.4 Which Implementation This Page Follows

Two bodies of R code carry the SRK name, and they are not the same thing.

This page follows the authors' scripts. R/02-srk-core.R is a port of them into the step-numbered layout the rest of this tutorial series uses, and every place it departs from them is marked DEPARTURE in the source and listed in §2.2.

2 Setup, Code & Data

2.1 Environment & Dependencies

Four packages, all on CRAN. Nothing else is needed: the singularity indices, the variance filter, the spatial block folds, the accuracy metrics and every figure are base R.

install.packages(c("sp", "gstat", "automap", "ranger"))

2.2 The Method in One File

R/02-srk-core.R holds the whole model. Three functions carry Eq. 1–10; everything else in the file is the five benchmarks, the metrics and the block folds, which exist so the cross-validation loop can treat all six models identically.

Eq. 2 — local intensity at every scale, in one pass

The authors' singularity_safe() loops over evaluation points and scales. The port keeps the definition and vectorises the search: one Chebyshev distance matrix per chunk of evaluation points, then one logical mask per scale. That is what makes the sensitivity analysis in step 6 affordable — the intensity matrix is computed once and re-cut for every shorter scale ladder.

srk_intensity <- function(x_ref, y_ref, z_ref, x_eval, y_eval, scales,
                          min_pts_per_scale = 3L, use_abs = TRUE,
                          chunk = 2000L) {
  if (use_abs) z_ref <- abs(z_ref)
  n_eval <- length(x_eval)
  J      <- length(scales)
  C      <- matrix(NA_real_, n_eval, J)
  # Chunked so a large evaluation set never allocates one huge distance matrix.
  for (ii in split(seq_len(n_eval), ceiling(seq_len(n_eval) / chunk))) {
    d <- pmax(abs(outer(x_eval[ii], x_ref, "-")),
              abs(outer(y_eval[ii], y_ref, "-")))    # square window = Chebyshev ball
    for (j in seq_len(J)) {
      inw <- d <= scales[j]
      cnt <- rowSums(inw)
      val <- as.numeric(inw %*% z_ref) / cnt
      val[cnt < min_pts_per_scale] <- NA_real_       # too few neighbours: no estimate
      C[ii, j] <- val
    }
  }
  C
}

Eq. 3–4 — the index is a slope plus two

srk_alpha <- function(C, scales, min_valid_scales = 2L, eps = 1e-9) {
  logr <- log(scales)
  vapply(seq_len(nrow(C)), function(i) {
    cr <- C[i, ]
    ok <- is.finite(cr)
    if (sum(ok) < min_valid_scales) return(NA_real_)   # caller substitutes 2
    yv <- log(ifelse(cr[ok] > 0, cr[ok], eps))         # a window of zeros is possible
    xv <- logr[ok]
    xd <- xv - mean(xv)
    sum(xd * (yv - mean(yv))) / sum(xd * xd) + 2
  }, numeric(1))
}

tests/test-01 transcribes the authors' loop verbatim and checks that these two functions return the same numbers to 1e−10, including which locations come back as NA.

Eq. 5–9 — the model

srk_run() is the whole method: build or look up the singularity columns, drop the flat ones, fit the forest, krige what is left, add the two together. It shares its argument list with the five benchmarks, so SRK_MODELS can be iterated over.

srk_run <- function(train, test, yvar, xvars, sv = NULL, ...) {

  ## Step 1 — one singularity column per covariate (Eq. 1-4).
  ##   Covariates are observable at prediction locations, so the union of train
  ##   and test may supply the reference values without leaking the response.
  ...
  for (cc in sv_cols) {
    tr[[cc]][!is.finite(tr[[cc]])] <- neutral          # 2 = uniform field
    te[[cc]][!is.finite(te[[cc]])] <- neutral
  }

  ## Step 2 — variance filter, then the forest (Eq. 5-7).
  sv_keep <- srk_select_features(tr, sv_cols, sd_threshold)
  tri     <- srk_trend(tr, te, yvar, c(xvars, sv_keep), ntree = ntree,
                       seed = seed, use_coords = use_coords)
  tr$resid <- tr[[yvar]] - tri$trend_tr

  ## Step 3 — ordinary kriging of the residual, added back (Eq. 8-9).
  kr <- srk_krige_residuals(tr, te, "resid", fitter = fitter)

  list(model = "SRK", pred = tri$trend_te + kr$pred, trend = tri$trend_te,
       kriged = kr$pred, se = kr$se, sv_kept = sv_keep, ...)
}

SRK_MODELS <- list(OK = bench_ok, IDW = bench_idw, LM = bench_lm,
                   RF = bench_rf, RFK = bench_rfk, SRK = srk_run)

2.3 Project Structure & Configuration

SRK/
├── run-all.R                   # one command runs everything
├── config/project-config.R     # the only file you edit to port this elsewhere
├── R/
│   ├── 00-config.R             # paths, output file names
│   ├── 01-helpers.R            # logging, IO, figure devices
│   ├── 02-srk-core.R           # THE METHOD + the five benchmarks
│   ├── 10-simulate.R           # paper Sec. 2.3  -> Table 1, Figs. 1-2
│   ├── 20-case-data.R          # download, tidy, Spearman screen (Sec. 4.1-4.2)
│   ├── 30-singularity.R        # covariate singularity features (Sec. 4.3)
│   ├── 40-block-cv.R           # six models, 5-fold spatial block CV (Table 2)
│   ├── 50-seed-stability.R     # how much of the SRK-RFK gap is forest noise
│   ├── 60-sensitivity.R        # scale x threshold grid (Fig. 7)
│   ├── 90-tables.R             # LaTeX tables
│   └── p01..p06-*.R            # figures
├── data/                       # provenance.md; the CSV lands here on first run
├── results/                    # every number quoted on this page
├── tables/                     # the same numbers as LaTeX
├── figs/                       # png + pdf of every figure
├── env/                        # requirements.md, session-info.txt, runtimes.csv
└── tests/                      # 18 method-property + 29 reproduction checks

Porting SRK to another study area means editing one file. The block below is the part of config/project-config.R that carries the method's parameters; everything on this page is a function of these values.

# One entry per response variable. `xvars` is the covariate set the paper's
# Spearman screen (p < 0.05) retained; R/20-case-data.R re-runs the screen and
# reports whether it reproduces this set.
#
# The paper maps two trace elements; this tutorial follows Co end to end. Every
# step is a loop over ELEMENTS, so adding the second element back is one entry --
# uncomment the Zn block and re-run:
#
#   Zn = list(file = "dt_Zn_linear.csv", yvar = "Zn", unit = "ppm",
#             xvars = c("ss", "ifi", "hm", "elevation", "slope", "aspect")),
ELEMENTS <- list(
  Co = list(file = "dt_Co_linear.csv", yvar = "Co", unit = "ppm",
            xvars = c("ifi", "hm", "elevation"))
)

# -- Singularity features (paper Eq. 1-4) ------------------------------------
SV_SCALES         <- seq(2000, 20000, by = 2000)  # window half-widths, map units
MIN_PTS_PER_SCALE <- 3L      # a scale needs this many neighbours
MIN_VALID_SCALES  <- 2L      # an index needs this many valid scales
SV_NEUTRAL        <- 2.0     # value assigned when estimation fails
SV_SD_THRESHOLD   <- 0.5     # drop sv features flatter than this

# -- Trend model and residual kriging ----------------------------------------
RF_NTREE      <- 500L
RF_SEED       <- 42L         # the authors' scripts leave this unset
RF_USE_COORDS <- TRUE        # their case-study code feeds x, y to the forest
VARIOGRAM_FIT <- "automap"   # as published

# -- Validation ---------------------------------------------------------------
CV_FOLDS      <- 5L
CV_BLOCK_SIZE <- 15000       # spatial block edge, map units
CV_SEED       <- 123L
SEED_REPEATS  <- 20L         # forest seeds used by R/50-seed-stability.R

2.4 Data Contract

SRK needs one table per response: point observations in a projected CRS with one column per covariate. That is the entire input.

The case-study table is not redistributed with this tutorial. R/20-case-data.R downloads it (about 130 KB) from the authors' repository on first run and writes data/provenance.md recording where it came from: GSWA Geochemistry (DMIRS-047) for the concentrations, Geoscience Australia's Surface Geology and DEM for the covariates. Set DATA_SOURCE <- "local" and drop your own CSVs in data/ to run everything on your own study area.

3 Reproduction Pipeline, Results & Validation

3.1 Pipeline Execution

One command runs the seven steps and the six figure scripts, in order, from a clean folder. Step 20 fetches the sample table the first time; nothing else touches the network.

Rscript run-all.R          # full pipeline, about 150 s
Rscript run-all.R test     # 18 + 24 verification checks, about 8 s
Rscript run-all.R 10 40    # only the named steps
Rscript run-all.R figures  # only the figure scripts

3.2 Core Analytical Steps

Step 1 · R/10-simulate.R · 0.5 s

Three distributions, one covariate, and where SRK's margin comes from

The paper's simulation is a controlled test of one claim: SRK's advantage over ordinary kriging should grow as the response departs from Gaussian. A spherical field on a 20 × 20 grid is the response; a second field, mixed 85/15 with the first, is the covariate; both are then pushed through a cube and an exponential to make a skewed and a long-tailed version. Seed 42, a 70/30 split, and the two models are scored on the held-out 120 cells.

R/10-simulate.R Copy

set.seed(SEED)                                     # 42, as in the released script
xy <- expand.grid(x = seq_len(SIM_SIDE), y = seq_len(SIM_SIDE))
sp::gridded(xy) <- ~x + y

sim_field <- function(range) {
  g <- gstat::gstat(formula = z ~ 1, dummy = TRUE, beta = 0,
                    model = gstat::vgm(psill = 1, model = "Sph", range = range),
                    nmax = 20)
  stats::predict(g, newdata = xy, nsim = 1)@data$sim1
}
y_raw     <- sim_field(SIM_RANGE_Y)                # range 10
noise_raw <- sim_field(SIM_RANGE_X)                # range  8

x_raw  <- SIM_MIX * y_raw + (1 - SIM_MIX) * noise_raw   # cor(X, Y) = 0.98
y_base <- y_raw + SIM_SHIFT                             # keep the transforms positive
x_base <- x_raw + SIM_SHIFT

scenarios <- list(
  normal      = list(Y = y_base,             X = x_base),
  skewed      = list(Y = y_base^3 / 100,     X = x_base^3 / 100),
  `long-tail` = list(Y = exp(y_base) / 1000, X = exp(x_base) / 1000))

The scales here are grid cells, not metres — SV_SCALES_SIM is 1 … 10 — and the SD filter is switched off, because a single covariate cannot be compared against anything. That matches the authors' simulation code, which keeps its one singularity column unconditionally.

Demo run

Published reference · SRK Fig. 1

Demo run

Published reference · SRK Fig. 2

Simulated field	R2			RMSE			MAE
	OK	SRK	paper SRK	OK	SRK	paper SRK	OK	SRK	paper SRK
Normal	0.876	0.972	0.972	0.361	0.171	0.172	0.283	0.137	0.137
Skewed	0.833	0.968	0.966	0.361	0.158	0.162	0.265	0.117	0.120
Long-tail	0.736	0.940	0.940	0.178	0.085	0.085	0.107	0.056	0.056

◆What this step establishes

The gain over ordinary kriging is 0.096 in R² on the near-normal field, 0.135 on the skewed one and 0.205 on the long-tailed one — monotonically increasing, which is the paper's central simulation claim and is asserted as a test in tests/test-02. Note also what the simulation does not establish: the covariate here correlates 0.98 with the response. Real covariates do not, and the case study below is where the margin gets realistic.

Step 2 · R/20-case-data.R · 0.5 s

Screen the covariates against the response

The lithology polygons have already been reduced to proximity variables by Eq. 14 — 1 inside the unit, falling linearly to 0 across a 10 km buffer. What remains is the paper's Sec. 3.2.1 screen: keep only covariates with a significant monotonic relation to the element.

R/20-case-data.R Copy

sc <- do.call(rbind, lapply(CANDIDATE_XVARS, function(xv) {
  ct <- suppressWarnings(stats::cor.test(d[[yvar]], d[[xv]], method = "spearman"))
  data.frame(element = elem, covariate = xv,
             rho = unname(ct$estimate), p = ct$p.value,
             significant = ct$p.value < 0.05,
             used_by_paper = xv %in% cfg$xvars)
}))

Covariate	Spearman ρ	p	Screen	Used by the paper
hm haematite proximity	+0.431	2e−46	kept	yes
elevation	+0.360	7e−32	kept	yes
ifi felsic intrusive proximity	−0.173	4e−08	kept	yes
imi mafic intrusive proximity	+0.031	0.330	dropped	no
slope	−0.021	0.508	dropped	no
ss sedimentary proximity	+0.005	0.868	dropped	no
aspect	−0.001	0.972	dropped	no

Published reference · SRK Fig. 5

Element	n	mean	median	max	sd	CV	skewness
Co (ppm)	998	13.30	8.94	112.2	13.21	0.994	2.58

Step 3 · R/30-singularity.R · 0.3 s

Build the singularity features, then throw one away

For each retained covariate, the intensity is measured in ten square windows of half-width 2 km to 20 km, and α is the slope of log intensity on log scale plus two. The reference points are the sample locations; the evaluation points are the same locations here, but the same call would evaluate the index on a prediction grid without changing anything.

R/30-singularity.R Copy

for (xv in cfg$xvars) {
  C <- srk_intensity(d$x, d$y, d[[xv]], d$x, d$y, SV_SCALES,
                     min_pts_per_scale = MIN_PTS_PER_SCALE)      # Eq. 2
  a <- srk_alpha(C, SV_SCALES, min_valid_scales = MIN_VALID_SCALES)  # Eq. 3-4
  a[!is.finite(a)] <- SV_NEUTRAL                                 # uniform field
  Cs[[xv]] <- C                    # cached: step 60 re-cuts the ladder from this
  sv[[paste0("sv_", xv)]] <- a
}

Feature	mean α	SD	range	% α < 2	r with Co	SD filter
sv(ifi)	7.02	6.77	0.77 – 21.97	28.1	+0.219	retained
sv(hm)	5.43	6.18	1.08 – 21.95	43.6	−0.227	retained
sv(elevation)	2.00	0.046	1.86 – 2.16	44.0	−0.172	dropped

◆Why terrain has no singularity to speak of

Elevation is smooth at 2–20 km. Its mean absolute value barely changes as the window grows, so the log-log slope is almost zero and α sits on 2 — the definition's own value for a uniform field. Lithology proximity is the opposite: it is 1 inside a unit and 0 more than 10 km outside, so the window mean changes sharply with scale and α ranges over an order of magnitude. The SD filter is really a smoothness detector, and the paper's 0.5 is doing exactly the job its Sec. 2.2.1 describes. (In the article's Zn model, where slope and aspect are also covariates, they behave the same way and are dropped for the same reason — SD 0.32 and 0.21.)

Demo run

Published reference · SRK Fig. 6

Feature importance is worth quoting directly, because it is the clearest evidence that the singularity columns earn their place. Fitting SRK on all 998 Co samples gives, in order: hm (23 % of total importance), elevation (17 %), sv_hm (17 %), x (14 %), y (14 %), sv_ifi (11 %), ifi (4 %). Both singularity columns outrank the raw ifi variable they were built from, and together they carry more of the forest than the two coordinates do.

Step 4 · R/40-block-cv.R · 7 s

Six models, five spatial blocks, no leakage

The samples cluster heavily, so a random split would leave almost every test point ringed by its own training neighbours and would flatter every kriging-based model. The paper instead cuts the study area into 15 km squares, shuffles the squares and deals them round-robin to five folds: 120 blocks, 5 folds, whole blocks held out together.

R/02-srk-core.R Copy

srk_block_folds <- function(data, k = 5L, block_size = 15000, seed = 123L) {
  df <- as.data.frame(data)
  bx <- floor((df$x - min(df$x, na.rm = TRUE)) / block_size)
  by <- floor((df$y - min(df$y, na.rm = TRUE)) / block_size)
  block <- paste(bx, by, sep = "_")
  set.seed(seed)
  ub <- sample(unique(block))                       # shuffle whole blocks
  fold <- as.integer(stats::setNames(((seq_along(ub) - 1L) %% k) + 1L, ub)[block])
  attr(fold, "n_blocks") <- length(ub)
  fold
}

R/40-block-cv.R Copy

fold <- srk_block_folds(d, CV_FOLDS, CV_BLOCK_SIZE, CV_SEED)
for (k in seq_len(CV_FOLDS)) {
  tr <- d[fold != k, ]; te <- d[fold == k, ]
  for (mname in names(SRK_MODELS)) {                # OK IDW LM RF RFK SRK
    out <- SRK_MODELS[[mname]](tr, te, yvar, cfg$xvars, sv = sv)
    met <- srk_metrics(te[[yvar]], out$pred)        # Eq. 11-13
    ...
  }
}

Model	This run			Published Table 2
	R2	RMSE	MAE	R2	RMSE	MAE
SRK	0.359	10.39	6.42	0.353	10.45	6.47
RFK	0.340	10.53	6.45	0.345	10.48	6.41
RF	0.337	10.55	6.51	0.342	10.50	6.46
IDW	0.266	11.15	7.42	0.266	11.15	7.42
LM	0.237	11.33	7.48	0.237	11.33	7.48
OK	0.152	12.02	7.99	0.148	12.04	8.03

Demo run

Published reference · SRK Fig. 8

◆The result that matters, and the one that needs an error bar

Matters. Ordinary kriging explains 15 % of the cobalt variance under spatial blocking; adding covariates through a forest takes that to 36 %, and reduces RMSE by 14 % and MAE by 20 %. That is a large, robust, unambiguous effect, and it is the same conclusion the article draws.

Needs an error bar. SRK against RFK is 0.020 in R² here, against 0.008 in the article — and the authors' own released cv_summary_Co.csv has the two the other way round, with RFK ahead by 0.001. Three artefacts, three answers. A gap that size cannot be read off a single run of a model containing an unseeded random forest. That is what step 5 is for.

Step 5 · R/50-seed-stability.R · 85 s · not in the paper

Is the SRK lead real? Put an error bar on it

Three of the six models contain a random forest, and the released scripts call ranger() without a seed. So RF, RFK and SRK produce a different number every time the analysis is run — and the differences the article reports between them are small. There is direct evidence of how small: the authors' own repository ships cv_summary_Co.csv with SRK at R² = 0.344, behind RFK at 0.345, while the article's Table 2 has SRK first at 0.353. Same code, same data, different run.

This step holds the folds and the data fixed and repeats the entire cross-validation over 20 forest seeds, then asks a paired question: on how many of those seeds does SRK actually finish ahead of RFK?

R/50-seed-stability.R Copy

fold <- srk_block_folds(d, CV_FOLDS, CV_BLOCK_SIZE, CV_SEED)   # fixed once

for (s in seq_len(SEED_REPEATS)) {
  for (mname in c("RF", "RFK", "SRK")) {                       # the seeded three
    fm <- do.call(rbind, lapply(seq_len(CV_FOLDS), function(k) {
      out <- SRK_MODELS[[mname]](d[fold != k, ], d[fold == k, ],
                                 yvar, cfg$xvars, sv = sv, seed = s)
      srk_metrics(d[[yvar]][fold == k], out$pred)
    }))
    ...   # mean over the five blocks, one row per (seed, model)
  }
}

Demo run

Model	mean R2	SD	min	max	published	authors' CSV
SRK	0.3580	0.0030	0.3522	0.3618	0.353	0.344
RFK	0.3401	0.0022	0.3349	0.3461	0.345	0.345
RF	0.3369	0.0022	0.3307	0.3419	0.342	0.342

◆For cobalt, the SRK advantage survives the error bar

SRK finishes ahead of RFK on 20 of 20 seeds, by +0.018 in R² on average, with the smallest margin (+0.010) still more than three times the seed-to-seed SD. That is a stronger result than the article claims for cobalt (+0.008), and it is the kind of statement a paired design licenses and a single run does not.

The comparison that never needed an error bar is SRK against the geostatistical benchmarks: +0.21 in R² over ordinary kriging is seventy times the seed noise. Reporting the spread separates the claim that is bulletproof from the one that had to be earned.

!Not all of the offset is the seed

The authors' released SRK value of 0.344 is about five seed-SDs below the mean this pipeline produces, so the forest seed alone does not explain it. SRK's third step is an autoKrige call on the residuals, and the same automap version difference that shifts ordinary kriging on two of the five folds (§3.5) shifts SRK too. Both effects are real, both are small compared with SRK's margin over OK, IDW and LM, and both are worth naming rather than averaging away.

✓Carry this into your own work

Any comparison between two models that both contain a random forest needs this. It costs one loop and it changes what you are entitled to write. If the gap between your method and the nearest benchmark is smaller than the spread across seeds, report the spread — a reviewer who re-runs your code will find it anyway.

Step 6 · R/60-sensitivity.R · 56 s

How much do the two parameters matter?

SRK has exactly two tuning choices: how far up the scale ladder the regression runs, and how flat a singularity feature may be before it is dropped. The paper varies both and reports that performance stays within roughly ±5 % of the baseline. This step re-runs the full cross-validation over a 6 × 5 grid of the two — 30 cells, 150 SRK fits.

Truncating the ladder needs no new neighbourhood search. Step 3 cached the per-scale intensity matrix, so a shorter ladder is a column subset and a re-fitted slope:

R/60-sensitivity.R Copy

Cs <- readRDS(file.path(DERIVED, sprintf("intensity-%s.rds", elem)))

for (mx in SENS_MAX_SCALES) {              # 10, 12, 14, 16, 18, 20 km
  keep_scale <- SV_SCALES <= mx
  for (xv in cfg$xvars) {                  # re-fit Eq. 3-4 on the shorter ladder
    a <- srk_alpha(Cs[[xv]][, keep_scale, drop = FALSE], SV_SCALES[keep_scale],
                   min_valid_scales = MIN_VALID_SCALES)
    a[!is.finite(a)] <- SV_NEUTRAL
    sv[[paste0("sv_", xv)]] <- a
  }
  for (thr in SENS_THRESHOLDS) { ... }     # 0.3, 0.4, 0.5, 0.6, 0.7
}

Demo run

Published reference · SRK Fig. 7

baseline R2	best	at	ΔR2 range	ΔRMSE range	cells within ±5 %
0.359	0.365	12 km / any threshold	−3.3 % to +1.7 %	−0.8 % to +0.8 %	30 / 30

◆The SD threshold is inert here, and it is worth saying why

Read the heatmaps down a column rather than across: every threshold from 0.3 to 0.7 gives the identical number. That follows from step 3 — the retained features have SD 6.18 and 6.77 and the dropped one 0.046, so nothing in the 0.3–0.7 range changes which features are kept. The threshold is not tuning anything on this data; it is a switch sitting far from its own decision boundary. That is a comfortable place to be, and quite different from the CRAN package's simulated data, where the same 0.5 lands on the wrong side of the only feature there is (§1.4).

The scale ladder does move the answer, mildly and non-monotonically: cobalt peaks at 12 km (0.365) against 0.359 at the published 20 km. That is a 1.7 % gain — the same order as the seed noise in step 5, so the location of the optimum is not sharply identified either. The article's own sensitivity analysis reaches the same verdict from the other direction: nothing in the grid is far from anything else.

3.3 What the Article Shows Beyond Cross-Validation

The figure at the top of this page and the two below are the part of the article this pipeline deliberately does not regenerate. They all need the 97,234-cell, 500 m prediction grid — an 8.1 MB covariate table on which the singularity indices, the forest trend and the kriging system must all be evaluated. Adding it is not conceptually different from what step 4 already does: srk_run() takes a pred_data-shaped test frame, and R/30-singularity.R already evaluates α at arbitrary locations. It is a change of scale, not of method — and the point of a tutorial that stops at cross-validation is that every number on this page can be checked in two and a half minutes on a laptop.

Published reference · SRK Fig. 10

Published reference · SRK Fig. 11

3.4 Reading the Results Together

Six steps produce one argument, and it is worth stating in the order a reader will want it.

1. The mechanism is real and it is visible. Singularity indices built from lithology proximity vary over an order of magnitude, correlate with cobalt at r = +0.22 and −0.23, and the forest ranks sv(hm) third of seven features — above both coordinates, and above the ifi covariate that sv(ifi) was derived from.
2. Under controlled conditions the gain is large and behaves as predicted. On the simulation, +0.096, +0.135 and +0.205 in R² over ordinary kriging as the response goes from near-normal to long-tailed. Monotone, as the paper claims.
3. On real data the gain over kriging is large; over a forest it is small but consistent. SRK beats OK by +0.21 in R² and cuts RMSE by 14 %. Against RFK the margin is +0.018 — small, but positive on every one of 20 forest seeds.
4. The residual really does get easier to krige. On the full-sample fit, residual SD falls from 4.90 to 4.72 ppm when the two singularity columns are added — which is the mechanism the article's Fig. 11 turns into lower kriging variance.
5. Nothing here is delicate. Thirty parameter combinations stay within ±5 % of the baseline, and the SD threshold is nowhere near a decision boundary. The quantity that does need care is the SRK–RFK margin, and step 5 measures it rather than asserting it.

3.5 Validation

Rscript run-all.R test runs two suites in about eight seconds. test-01 checks the implementation against the definitions; test-02 checks the outputs against the published article and against the authors' own released fold-level numbers.

4 Adaptation, Writing & Reproducibility

4.1 Port to Your Domain

SRK needs a point response, covariates observable everywhere, and projected coordinates. Nothing about it is specific to geochemistry — soil properties, air pollution, groundwater chemistry, biomass and disease incidence all have the shape it wants. Four edits to config/project-config.R and the pipeline runs on your data.

Design decisions worth making deliberately

• Coordinates in the forest or not. RF_USE_COORDS = TRUE reproduces the published case study; FALSE follows Eq. 5. Leaving them in lets the forest learn location directly, which competes with the job the kriging step is there to do — for cobalt the two coordinates take 28 % of the forest's total importance between them. If your interest is in why the response varies rather than only in prediction accuracy, run it both ways and report both.
• Reference set for the index. Here the reference points are the sample locations, so the index inherits the sampling design. If your covariate exists as a raster, evaluating the intensity on the raster instead removes that dependence entirely and is closer to the theory in Eq. 1–2. That is a real modelling choice, not a detail.
• Anisotropy. Both the article and this pipeline krige isotropically. In a structurally controlled setting — a fault trend, a prevailing wind, a river network — fitting an anisotropic variogram to the residual is the obvious next thing to try; automap will not do it for you.
• Back-transformation. The response is modelled on its raw ppm scale, not logged, even though it is strongly right-skewed. If you log a skewed response, remember that exponentiating the prediction is biased low, and say what correction you applied.

4.2 Write the Paper

The article's own structure is the template, and this pipeline produces a table or a figure for every part of it.

4.3 Final Reproducibility Package

The data are not in the zip by design. They belong to the GSWA geochemistry database and Geoscience Australia, they are already published by the article's authors, and a tutorial should point at the source rather than fork it. R/20-case-data.R fetches the one table it needs in under a second and records where it came from.

unzip srk-code-and-data.zip && cd SRK
Rscript -e 'install.packages(c("sp","gstat","automap","ranger"))'
Rscript run-all.R          # about 150 s
Rscript run-all.R test     # 18 + 24 checks

Sources

• Ren K, Song Y, Chen M, Yu Q (2026). A singularity regression kriging for spatial prediction. GIScience & Remote Sensing 63(1):2690341. doi:10.1080/15481603.2026.2690341 · local PDF
• Authors' code and data: github.com/renkaigis/Singularity_Regression_Kriging · Figshare
• Cheng Q (2017). Singularity analysis of geo-anomalies — the source of the scaling relation in Eq. 1–4.
• Breiman L (2001). Random forests. Machine Learning 45:5–32.
• Hiemstra P et al. (2009). automap; Pebesma E (2004). gstat; Wright M & Ziegler A (2017). ranger.
• GSWA Geochemistry database (DMIRS-047); Geoscience Australia, Surface Geology of Australia (2012) and Digital Elevation Model (2015).