Reproducing Geographically Optimal Similarity

1 Method Overview & Reproduction Scope

1.1 Core Idea

Most spatial prediction starts from distance. Kriging borrows from nearby samples; a regression borrows from a fitted equation. The geographical similarity principle — the Third Law of Geography — borrows from somewhere else entirely: similar geographical configurations produce similar values, wherever they happen to be. Two sites 60 km apart with the same terrain, vegetation, soil chemistry and distance to mining can inform each other better than two adjacent sites that differ in all of them.

A geographical configuration is just the vector of covariate values at a location. Similarity compares two such vectors directly — there is no regression coefficient anywhere in the model, and no explicit relationship between response and covariates is ever estimated. That is the point: the method describes the comprehensive degree of approximation of a geographical structure instead of an explicit relationship between variables.

The earlier basic configuration similarity (BCS) model applies that principle by averaging over every observation, weighted by similarity. GOS adds one idea and one parameter: at each unknown location, keep only the most similar κ share of observations and discard the rest. Observations with low similarity are not neutral — they are noise dressed up as evidence, and averaging them in drags every prediction towards the study-area mean.

• Research problem: similarity-based prediction uses all samples at every location, so large data sets cost more to compute and predict worse, because unrelated observations dilute the informative ones.
• One-line contribution: a cross-validated threshold λ that selects, at each unknown location, the small share of observations with optimal similarity — plus an uncertainty measure read straight off the retained similarities.
• Why it matters here: on the packaged Zn data, λ = 0.08 — 8 % of 885 observations, about 71 samples per cell — beats using all of them on MAE in 50 of 50 cross-validation splits, and lowers mean prediction uncertainty at ζ = 0.9 by 67 %.

Published reference · GOS Paper Fig. 1 — the four stages this page follows. Stage 1 selects covariates and turns them into a geographical configuration at every location, sampled and unsampled. Stage 2 computes a per-covariate similarity E, composes them into a single sample-scale similarity S with the operator P, and collects one similarity vector S per unknown location. Stage 3 finds the threshold λ that minimises cross-validation error and keeps only the observations above it — the red dots in the bottom right panel. Stage 4 predicts and maps uncertainty. Note that no arrow in this diagram carries a distance. Song (2023), Mathematical Geosciences 55:295–320 (open access).

1.2 Method Logic

GOS runs in four stages, and each stage is a small number of equations. Equation numbers below are the paper's, and every one of them is exercised by the pipeline in §3.

Stage 1 — characterising the geographical configuration (Eq. 1)

Two screens, in order. First keep the covariates significantly correlated with the response; then drop the most collinear one repeatedly until every variance inflation factor is under 4, the paper's conservative threshold. What survives is the configuration:

where e_i is the value of covariate X_i at a given location. The configuration must be computable at the unsampled locations too, which is what rules out covariates you only measured at the samples.

Stage 2 — assessing similarity (Eq. 2–6)

For an unknown location v_β and an observation u_α, each covariate contributes its own Gaussian similarity, and the covariate-scale similarities are composed into one number:

Three details in there decide everything downstream. σ is the standard deviation of the covariate over the observation and prediction locations together, so the scale is set by the whole study area, not by the samples alone. δ is a root-mean-square deviation that makes the Gaussian bandwidth adaptive rather than fixed. And P is the minimum operator — the paper adopts Zhu et al.'s min, so a single badly matched covariate is enough to make two locations dissimilar no matter how well the others agree. Similarity is a conjunction, not an average.

Collecting the similarities of one unknown location against every observation gives a vector (Eq. 5), and the similarity-weighted mean of the observed values is the BCS prediction:

Stage 3 — determining the optimal similarity (Eq. 7–10)

κ is the share of observations retained at each prediction location, defined as one minus the quantile probability τ of the similarity vector:

κ = 1 keeps everything and is exactly BCS. For each candidate κ the observations are split 50/50, the held-out values are predicted with Eq. 6, and the cross-validation error is recorded; repeating the split many times and taking the minimum gives λ:

The similarity vector is then truncated to the observations above the corresponding similarity threshold S_λ (Eq. 10). λ is a property of the data set, not of a location: one number is searched once and used everywhere.

Stage 4 — prediction and uncertainty (Eq. 11–12)

The GOS prediction is the same weighted mean as Eq. 6, taken over the retained observations only:

and the uncertainty is one minus a quantile of the retained similarities:

with ζ ∈ {0.9, 0.95, 0.99, 0.995, 0.999, 1}. This is the model's most elegant move: the uncertainty needs no error model and no second pass, because the quality of the evidence is already in the similarities. At ζ = 1 the expression collapses to 1 − max(similarity) — "how good is the single best match?" — and at that setting GOS and BCS give identical answers, since truncating the vector never removes its maximum. The pipeline reproduces that identity exactly (§3.2, Step 4).

Published reference · GOS Paper Table 1 and Fig. 3 — the case study's five steps: preprocess the trace-element data, build and select the geological variables, collect the natural and social variables, fit and cross-validate GOS against four other models, then predict and evaluate. The tutorial reproduces steps 1, 3, 4 and 5; step 2, the geological variables, is the part the packaged data cannot supply (§1.3). Table 1 above the diagram is the published summary of the raw Cu / Zn / Pb samples that §3.2 Step 1 rebuilds for Zn. Song (2023), Mathematical Geosciences 55:295–320 (open access).

1.3 Reproduction Scope

Everything tagged Demo run regenerates from this folder with Rscript run-all.R; everything tagged Published reference is a static image from the article. The demo runs end to end on R 4.6.0 with geosimilarity 3.9. Three things are worth stating plainly before any number below is quoted.

What does reproduce is the whole GOS estimation path and the qualitative findings: the correlation-and-VIF screen (which lands on exactly the covariate set the package vignette uses), the κ search of Eqs. 7–9 and its characteristic curve shape, the Eq. 11 prediction over all 13,132 grid cells, the Eq. 12 uncertainty at all six ζ levels including the exact GOS = BCS identity at ζ = 1, the collapse of BCS predictions towards the mean that GOS avoids, and the accuracy ordering GOS < BCS < MLR on both MAE and RMSE.

How to use this page

• To reproduce the demo: run Rscript run-all.R test, then Rscript run-all.R. Every number and figure below regenerates in under two minutes.
• To try GOS without installing anything: the online GOS calculator runs the same CRAN package in your browser, on this same example data or on a CSV of your own.
• To use GOS on your own data: §2.3 gives the data contract — a sample table and a prediction table sharing covariate columns — and §4.1 the porting steps.

2 Setup, Structure & Data Contract

2.1 Environment & Dependencies

This is the shortest dependency list in the series. The method is one CRAN package written by the article's author; everything else is base R.

# The method, and both input data sets
install.packages("geosimilarity")

2.2 Project Structure & Configuration

One entry script runs the five 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.

GOS/
├── gos.html                       # this guide
├── run-all.R                      # entry point: full run, single steps, tests
├── assets/                        # figures shown in this guide (web copies)
├── config/project-config.R        # the ONLY file to edit for a new domain
├── R/
│   ├── 00-config.R, 01-helpers.R  # paths, IO, CV splits, base-R figure helpers
│   ├── 10-prepare-data.R          # load zn + grid, log-transform, outlier screen
│   ├── 20-select-variables.R      # correlation + VIF  ->  the configuration (Eq. 1)
│   ├── 30-best-kappa.R            # the kappa search   ->  lambda (Eqs. 7-9)
│   ├── 40-gos-prediction.R        # grid prediction + uncertainty (Eqs. 11-12)
│   ├── 50-model-comparison.R      # MLR / BCS / GOS on identical CV splits
│   └── 90-tables.R                # LaTeX tables + session info
├── data/                          # written by step 10: verbatim CSV copies of the
│                                  #   packaged zn and grid, plus provenance.md
├── tests/                         # test-01 method properties; test-02 results regenerate
├── 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. It is read by R/00-config.R and used by every step; nothing else in the pipeline hard-codes a column name, a threshold or a palette.

# -- 2. Input data -----------------------------------------------------------
RESPONSE      <- "Zn"        # the raw response column (ppm)
RESPONSE_LOG  <- "logZn"     # name given to the log-transformed response
LOG_TRANSFORM <- TRUE        # paper Sect. 3.2, step 1

CANDIDATES <- c("Elevation", "Slope", "Aspect", "Water",
                "NDVI", "SOC", "pH", "Road", "Mine")

# -- 3. Outlier screening ----------------------------------------------------
OUTLIER_COEF <- 2.5          # removeoutlier(): coef x IQR of the log response

# -- 4. Variable selection (paper Sect. 2.2.1) -------------------------------
COR_ALPHA     <- 0.05        # keep covariates significantly correlated with y
VIF_THRESHOLD <- 4           # the paper's conservative multicollinearity limit

# -- 5. Optimal similarity threshold (Eqs. 7-9) ------------------------------
KAPPA_GRID        <- c(seq(0.01, 0.10, 0.01), seq(0.2, 1.0, 0.1))
BESTKAPPA_NREPEAT <- 10      # cross-validation repeats inside gos_bestkappa()
BESTKAPPA_NSPLIT  <- 0.5     # training share of each split

# -- 6. Model comparison (paper Table 3) -------------------------------------
CV_REPEATS <- 50             # the paper's 50 repeated 50/50 splits
CV_MODELS  <- c("MLR", "BCS", "GOS")

# -- 7. Compute --------------------------------------------------------------
CORES <- 1                   # one core: reproducible, and matches the web app
SEED  <- 42

2.3 Data Contract

GOS needs two plain tables: the sample table, carrying the response and the covariates at the observation locations, and the prediction table, carrying the same covariates wherever you want a value. Neither needs to be a spatial object.

Units, and why they matter less than you would think

In the demo, distances (Water, Road, Mine) are in kilometres, Elevation in metres, Slope and Aspect in degrees, NDVI and pH dimensionless, and SOC is a mass fraction. Equation 3 divides every deviation by that covariate's own σ, so a covariate's units cancel and rescaling one from metres to kilometres changes nothing. What does matter is that a covariate be measured the same way in both tables — a sample column in metres against a grid column in kilometres is a silent, catastrophic error the model cannot detect.

3 Reproduction Pipeline, Results & Validation

3.1 Pipeline Execution

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

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

Project : gos-zn
Domain  : trace element (Zn) concentration, Leonora mining region, WA
== Step 1/5  Prepare data ======================================
   zn:   894 samples x 12 columns
   grid: 13132 prediction cells x 12 columns
   outlier screen (coef = 2.5): 9 of 894 rows dropped, 885 retained
   Shapiro-Wilk W on log Zn (screened): 0.9860
== Step 2/5  Select variables ==================================
   VIF screen kept 7 of 7 covariates (max VIF 2.608)
   selected: NDVI, SOC, Mine, Water, pH, Slope, Road
   vignette: Slope, Water, NDVI, SOC, pH, Road, Mine
   the two sets are identical
== Step 3/5  Optimal similarity threshold ======================
   kappa grid: 19 values from 0.01 to 1.00, 10 repeats of a 50/50 split
   lambda = 0.08  (cross-validation RMSE 0.6668)
   RMSE at kappa = 1 (BCS) is 0.6729 — GOS lowers it by 0.91%
   only 8% of the 885 observations are used at each prediction location
== Step 4/5  GOS prediction on the 1 km grid ===================
   885 observations -> 13132 grid cells, kappa = 0.08
   GOS predictions: mean 3.353, sd 0.365, range 2.396-4.302 (log Zn)
   BCS predictions: mean 3.403, sd 0.138, range 2.396-4.075 (log Zn)
   BCS spread is 38% of the GOS spread — the clustering the paper reports
   at zeta = 1 the two models agree to 0.00e+00, as the paper notes
== Step 5/5  Model comparison ==================================
   [screened] 885 samples, 50 repeats, 442 training / 443 testing rows each
   MLR  MAE 0.5456  RMSE 0.6865   GOS reduces MAE +5.8%, RMSE +4.2%
   BCS  MAE 0.5220  RMSE 0.6653   GOS reduces MAE +1.5%, RMSE +1.1%
   GOS  MAE 0.5140  RMSE 0.6578
   RMSE ranking, screened  : GOS < BCS < MLR
   RMSE ranking, unscreened: GOS < BCS < MLR
Finished in 107.0 s

3.2 Core Analytical Steps

Six stages: purpose → code → output → result → how to read it. Steps 4 and 5 both come out of R/40-gos-prediction.R; the rest map one-to-one onto the numbered scripts.

Step 1 · R/10-prepare-data.R

Prepare the response — log-transform, then screen

What this step does

Loads the two packaged tables, writes verbatim CSV copies of them, log-transforms the response because trace-element data are strongly right-skewed, and applies removeoutlier() to the transformed values. Both the screened and the unscreened sample sets are kept, because the paper and the package vignette disagree about whether to screen at all, and step 5 reports the comparison both ways.

Code

R — from R/10-prepare-data.R

zn   <- as.data.frame(geosimilarity::zn)     # 894 x 12, no download
grid <- as.data.frame(geosimilarity::grid)   # 13,132 x 12

zn$logZn <- log(zn$Zn)                       # paper Sect. 3.2, step 1
out_idx  <- geosimilarity::removeoutlier(zn$logZn, coef = 2.5)
d        <- zn[-out_idx, ]                   # 885 rows retained

Example output

Series	n	Mean	Min	Median	Max	σ	CV
Zn, ppm — all samples	894	42.55	2	29	353.5	37.11	0.872
log Zn — all samples	894	3.436	0.693	3.367	5.868	0.801	0.233
Zn, ppm — outliers screened	885	41.64	5	29	181	32.95	0.791
log Zn — outliers screened	885	3.438	1.609	3.367	5.198	0.772	0.225

results/response-summary.csv, in the shape of the paper's Table 1. The highlighted row is the sample set every later step uses. The paper reports 966 Zn samples with mean 43.52, median 30, max 353.5, σ 38.2 and CV 0.88 — the same population seen through 72 fewer samples, with an identical maximum.

Demo run fig01 — the response before and after preprocessing: histogram with kernel density (a) and normal Q-Q plot (b) of the 894 raw values in ppm, and the same two views (c, d) of the 885 log values that survive the screen. Look at the Q-Q panels: the raw data bend hard away from the line in the upper tail, the log values sit on it (Shapiro-Wilk W = 0.986, skewness 0.04). Dashed lines mark the mean.

Published reference · GOS Paper Fig. 4 — the published counterpart, for all three trace elements. The middle column is Zn: the same two-humped log density and the same near-linear Q-Q plot the demo reproduces from the packaged subset.

How to read this result

• The transform is doing the work, not the screen. The log takes CV from 0.872 to 0.233; the outlier screen then removes 9 rows and moves CV only from 0.233 to 0.225. Order matters — screening raw ppm values would have thrown away far more.
• The screen is a documented disagreement, not a default. The package vignette applies removeoutlier(coef = 2.5); the paper explicitly does not, because high trace-element values "may indicate the clusters of mineral deposits". Rather than pick a side silently, step 5 runs the whole comparison on both sample sets. The ranking is identical either way.
• Everything downstream is in log units. The predictions, the MAE and the RMSE are all log Zn. The back-transformed pred_zn_ppm column is a plain exp(pred) with no smearing or lognormal bias correction — fine for mapping, not for a mass balance.

# stage 1 — correlation with the log-transformed response
ct  <- stats::cor.test(d$logZn, d[[v]])       # for each candidate v
sig <- cors$variable[cors$p < 0.05]

# stage 2 — iterative VIF pruning at the paper's threshold of 4
pr       <- vif_prune(d, sig, VIF_THRESHOLD)  # VIF_i = 1 / (1 - R2_i)
selected <- pr$keep

Step 3 · R/30-best-kappa.R · the parameter

Find the optimal similarity threshold λ (Eqs. 7–9)

What this step does

Scores 19 candidate κ values by repeated 50/50 cross-validation and takes the minimum. This is the paper's actual contribution, reduced to one function call and one number.

Code

R — from R/30-best-kappa.R

f <- logZn ~ NDVI + SOC + Mine + Water + pH + Slope + Road

bk <- geosimilarity::gos_bestkappa(f, data = d,
        kappa   = c(seq(0.01, 0.10, 0.01), seq(0.2, 1.0, 0.1)),
        nrepeat = 10, nsplit = 0.5, cores = 1)

lambda <- bk$bestkappa      # 0.08
curve  <- bk$cvmean         # one mean RMSE per candidate kappa

Example output

Demo run fig02 — cross-validation RMSE of log Zn against κ, averaged over 10 repeated 50/50 splits: the full candidate range (a) and a zoom to κ ≤ 0.1 where the grid is finest (b). The open circle and dashed line mark λ = 0.08, RMSE = 0.6668. Read the shape, not just the minimum: a steep fall as the first few per cent of observations are added, a shallow trough, then a slow rise to a plateau at κ = 1.

Published reference · GOS Paper Fig. 7 — the published curves for Cu, Zn and Pb. The middle panel is Zn, with λ = 0.04. Same shape, same message, same order of magnitude; the demo's 0.08 differs because the configuration differs (§1.3).

κ	0.01	0.02	0.05	0.07	0.08	0.09	0.10	0.20	0.50	1.00
CV RMSE	0.6944	0.6797	0.6687	0.6670	0.6668	0.6670	0.6673	0.6691	0.6724	0.6729

Ten of the 19 rows of results/kappa-rmse.csv. λ = 0.08 with RMSE 0.6668; κ = 1, which is BCS, gives 0.6729, so choosing a threshold lowers cross-validation RMSE by 0.91 %. Only 8 % of the 885 observations — about 71 samples — contribute at each prediction location.

How to read this result

• The left arm of the curve is the more interesting half. RMSE falls 0.6944 → 0.6668 between κ = 0.01 and 0.08 — nine samples per cell is too few, and the prediction is unstable. The right arm rises only 0.6668 → 0.6729. Under-using the data costs about four times as much as over-using it, which is worth knowing before you pick a coarse κ grid.
• The minimum is shallow, and you should say so. Between κ = 0.05 and 0.10 the RMSE spans 0.6687 to 0.6673 — a 0.2 % range. λ is not sharply identified, and with a coarser grid or a different covariate set it moves easily. That is exactly why the paper's 0.04 and this run's 0.08 are the same finding: only a few per cent of observations are needed.
• λ is deterministic. gos_bestkappa() seeds each repeat internally (1…nrepeat), so every candidate κ is scored on identical splits and the answer does not depend on the outer set.seed(). tests/test-01 asserts this by running the search twice under deliberately different outer seeds.
• κ = 1 is BCS, exactly. The last point on the curve is not an approximation of the older model — it is the older model. That is what makes the comparison in step 6 fair: one implementation, one code path, one parameter changed.

Try the κ search in your browser

The online GOS calculator runs this exact search — the same CRAN geosimilarity package, compiled to WebAssembly — on the same example data or on a CSV of your own, and draws this same curve with λ marked. No R installation, nothing uploaded to a server. It is the fastest way to see how the curve reacts to a different covariate set before committing to a pipeline run.

Step 4 · R/40-gos-prediction.R

Predict the grid, and map the uncertainty (Eqs. 11–12)

What this step does

One gos() call predicts log Zn at all 13,132 grid cells at κ = λ = 0.08 and returns the Eq. 12 uncertainty at all six ζ levels in the same tibble. A second call at κ = 1 produces the BCS counterpart used in step 5. Together they take about ten seconds.

Code

R — from R/40-gos-prediction.R

g <- geosimilarity::gos(f, data = d, newdata = grid,
                        kappa = lambda, cores = 1)

# g holds: pred, uncertainty90, uncertainty95, uncertainty99,
#          uncertainty99.5, uncertainty99.9, uncertainty100

b <- geosimilarity::gos(f, data = d, newdata = grid,   # BCS, for step 5
                        kappa = 1, cores = 1)

Example output

Demo run fig03 — the GOS prediction of log Zn at κ = λ = 0.08 (a) and its uncertainty at ζ = 0.99 (b). Predictions run 2.40–4.30 in log Zn (11.0–73.9 ppm back-transformed), and the map has sharp local contrast: a dark north–south high through the centre of the study area against pale flanks. The uncertainty map is nearly blank by contrast (mean 0.033) with dark filaments where no observation shares the local configuration — those, not the prediction map, are the places to distrust. The uncertainty scale is clipped at its 99th percentile, marked ≥, because a handful of isolated cells reach 1.0 and a raw range would flatten the whole panel to one shade.

Demo run fig04 — the same uncertainty at all six ζ levels on one shared scale. Uncertainty falls monotonically as ζ rises, from a mean of 0.072 at ζ = 0.9 to 0.026 at ζ = 1, while the pattern barely moves: the poorly matched places stay in the same places. ζ is a strictness dial on the evidence, not a different map.

Published reference · GOS Paper Fig. 9 — published predictions for Cu (a), Zn (b) and Pb (c) from all five models. Compare the large GOS panel with the small BCS panel below it in each column: the same loss of contrast that fig05 quantifies here.

ζ	0.9	0.95	0.99	0.995	0.999	1
GOS mean Θ	0.0724	0.0545	0.0332	0.0295	0.0265	0.0258
BCS mean Θ	0.2203	0.1656	0.0801	0.0602	0.0351	0.0258
Reduction by GOS	67.1 %	67.1 %	58.5 %	51.0 %	24.5 %	0 %

Mean Eq. 12 uncertainty over the 13,132 cells, computed from results/grid-prediction.csv and grid-prediction-bcs.csv. The paper reports a 28.5–74.1 % reduction over ζ = 0.9 to 0.999; this run gives 24.5–67.1 % — same magnitude, same direction.

How to read this result

• ζ = 1 is an exact identity, and it holds to machine precision. At ζ = 1, Eq. 12 becomes 1 − max(similarity), and truncating a vector never removes its maximum — so GOS and BCS must agree exactly even though their predictions differ everywhere. The largest difference across all 13,132 cells is 0.00e+00. That is a genuine test of the implementation, not a coincidence, and tests/test-01 checks it on synthetic data too.
• The gap between GOS and BCS shrinks as ζ rises, for a mechanical reason. At ζ = 0.9, BCS's quantile falls in the long tail of poor matches it is obliged to keep, so its uncertainty is triple. By ζ = 0.999 both models are reading their few best matches, and the reduction collapses to 24.5 %. Report the ζ you used, or the number means nothing.
• Uncertainty here means "was there anything to learn from?", not a variance. A cell can have a confident prediction with high Θ if it is simply unlike everything observed. This is an extrapolation warning rather than a confidence interval, and it is available at no extra cost because it falls out of the similarities already computed.

You cannot predict the grid in chunks

σ in Eq. 3 is computed over the observation and prediction locations together, and δ in Eq. 4 divides by the number of prediction locations. Splitting newdata into blocks therefore changes both, and the numbers drift. This is a property of the method, not a bug: tests/test-02 re-predicts all 13,132 cells from scratch rather than checking a sample of them, and thinning the grid would shift the values slightly too.

Step 5 · R/40-gos-prediction.R · the clearest reproduced finding

Why the threshold matters — BCS collapses towards the mean, GOS does not

What this step does

Puts the κ = 1 prediction on the same colour scale as the κ = λ one, and draws the distribution of all 13,132 predictions from both. Two views of a single fact.

Example output

Demo run fig05 — the BCS prediction (a), drawn on the identical colour scale as the GOS map in fig03, and the density of all 13,132 predictions from both models (b). The map has lost its structure: what was a strong central high in fig03 is now a mottled mid-tone. The density panel says why — BCS is a single narrow spike on the observed mean (dashed line), GOS a broad multi-modal spread.

Prediction set	Mean	σ	Range (log Zn)	Range (ppm)	Within ±0.1 of the observed mean
GOS (κ = 0.08)	3.353	0.365	2.396 – 4.302	11.0 – 73.9	12.8 %
BCS (κ = 1)	3.403	0.138	2.396 – 4.075	11.0 – 58.8	70.0 %

Computed from results/grid-prediction.csv and grid-prediction-bcs.csv. BCS's spread is 38 % of GOS's, and 70 % of its cells sit within ±0.1 of the observed mean of 3.438 against 12.8 % of GOS cells — a 5.5-fold concentration.

Published reference · GOS Paper Fig. 10 — the published counterpart: density distributions of Cu, Zn and Pb predictions from all five models, with the vertical lines at the mean. The orange BCS curve spikes on the mean in every panel while the red GOS curve stays broad. That is the same fact, on the full data set, for all three elements.

How to read this result

• This is the paper's Fig. 10 finding, reproduced cleanly and with no caveats. Of everything on this page it is the most robust: it needs no cross-validation, no accuracy metric and no held-out sample — it is visible in the predictions themselves, on 13,132 cells.
• Averaging over everything is a shrinkage estimator you did not ask for. Weight 885 observations by similarities that are mostly small and mostly similar to each other, and the weighted mean converges on the unweighted one. The informative observations are outvoted. Discarding 92 % of them is what restores the local signal.
• Both models share a lower bound. Both ranges start at exactly 2.396 — some locations are dominated by the same single best-matching observation under either rule. The divergence is at the top: BCS reaches 4.075 where GOS reaches 4.302, so the high-Zn areas are precisely where BCS's dilution costs the most, and those are the areas a mining study cares about.
• Spread is not accuracy. A model could scatter its predictions wildly and be worse. Step 6 is what closes that loop: the extra contrast comes with lower error, not higher.

Step 6 · R/50-model-comparison.R

Accuracy comparison — MLR, BCS and GOS on identical splits

What this step does

Scores three predictors of log Zn over 50 repeated 50/50 splits, the paper's protocol. Splits are generated once and shared, so every model sees exactly the same training and testing rows in every repeat and the differences between models are paired rather than merely averaged. The whole comparison is then run a second time on all 894 unscreened samples, with λ re-derived for that sample set so GOS is not handicapped by a threshold tuned elsewhere.

Code

R — from R/50-model-comparison.R

one_repeat <- function(d, idx, lam) {
  tr <- d[idx, ]; te <- d[-idx, ]

  p_mlr <- stats::predict(stats::lm(f, data = tr), newdata = te)
  p_bcs <- geosimilarity::gos(f, tr, te, kappa = 1,   cores = 1)$pred
  p_gos <- geosimilarity::gos(f, tr, te, kappa = lam, cores = 1)$pred
  ...
}

splits <- cv_splits(nrow(d), CV_REPEATS, CV_SPLIT, SEED)  # shared by all three

Example output

Model	MAE	sd	RMSE	sd	MAE reduction by GOS	RMSE reduction by GOS
MLR — multivariate linear regression	0.5456	0.0124	0.6865	0.0154	5.81 %	4.18 %
BCS — basic configuration similarity (κ = 1)	0.5220	0.0125	0.6653	0.0174	1.54 %	1.12 %
GOS — geographically optimal similarity (κ = 0.08)	0.5140	0.0119	0.6578	0.0183	—	—

results/model-comparison.csv — 885 outlier-screened samples, λ = 0.08, 50 repeated 50/50 splits, all errors in log Zn units. RMSE ranking GOS < BCS < MLR. The paper reports, for Zn, GOS reducing MAE by 2.6 % and RMSE by 0.9 % against BCS, and 6.0 % / 1.2 % against MLR: same direction, same order of magnitude, with the MAE-against-MLR figure almost exact.

Demo run fig06 — MAE (a) and RMSE (b) over the 50 splits. White diamonds are the means reported in the table, the dashed line the GOS mean. Note how much the boxes overlap: split-to-split spread (sd ≈ 0.012 for MAE) is larger than the gap between BCS and GOS (0.008), which is exactly why the models must be scored on identical splits and compared pairwise.

Published reference · GOS Paper Fig. 8 — the published counterpart, with five models rather than three: ordinary kriging and regression kriging are the two columns this tutorial does not reproduce (§1.3). The MLR, BCS and GOS bars are the ones being compared here.

Robustness — the paper's own preprocessing

The vignette screens outliers; the paper explicitly does not. Re-running everything on all 894 samples, with λ re-derived on that sample set (it comes out at 0.08 again), gives:

Model	MAE	sd	RMSE	sd	MAE reduction by GOS	RMSE reduction by GOS
MLR	0.5640	0.0130	0.7184	0.0168	5.66 %	3.69 %
BCS	0.5401	0.0150	0.6994	0.0185	1.48 %	1.08 %
GOS	0.5321	0.0136	0.6919	0.0185	—	—

results/model-comparison-no-outlier-screening.csv — 894 samples, no screening, λ re-derived = 0.08, 50 repeats. Every error rises by roughly 0.02–0.03, because the nine retained extremes are genuinely hard to predict, but the ordering GOS < BCS < MLR and the size of GOS's advantage are unchanged. The preprocessing choice does not carry the result.

Paired over the shared splits

Averages can hide a model that wins narrowly on most splits and loses badly on a few. These do not. On the 885 screened samples, GOS has the lower error than MLR on 50/50 splits by MAE and 50/50 by RMSE, and than BCS on 50/50 by MAE and 47/50 by RMSE. On the 894 unscreened samples: 50/50 and 50/50 against MLR, 49/50 and 47/50 against BCS. Counted from results/model-comparison-repeats.csv, which holds all 300 per-repeat rows.

How to read this result

• The GOS-vs-BCS gap is small, consistent and the point of the paper. 1.5 % MAE is not dramatic, but it is the same 1.5 % on essentially every split, and it comes with a 92 % cut in the observations used per prediction and the restored map contrast of step 5. Cheaper, sharper and slightly more accurate.
• The GOS-vs-MLR gap is what justifies the similarity framing. Both models see exactly the same seven covariates. MLR insists on an explicit linear relationship; GOS never estimates one, and still predicts better by 5.8 % MAE on all 50 splits. The comparative structure of a configuration carries information that a fitted coefficient does not.
• MAE gaps exceed RMSE gaps throughout. 5.81 vs 4.18 % against MLR, 1.54 vs 1.12 % against BCS. GOS's advantage is in typical accuracy rather than in the worst cases — expected, since the hardest cells are the ones with no similar observation anywhere, which is precisely what its uncertainty map flags.
• What this comparison does not say. It says nothing about kriging, which is not run here (§1.3). On these covariates GOS beats a linear model and beats averaging over everything; whether it beats a geostatistical model on this data set is a question this tutorial does not answer.

3.3 Reading the Results Together

3.4 Validation

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

What test-01 pins down

• Eq. 6 is transcribed and matched exactly. The test writes out Eqs. 2–6 by hand — the Gaussian E_i, the σ over observation and prediction locations together, the δ normalisation as the package computes it, the min operator, the weighted mean — and compares it with gos(kappa = 1). Maximum difference: 0.0e+00. That single check verifies the minimum-operator composition and the Eq. 6 weighted mean at once.
• Eq. 12 behaves like a probability. All six uncertainty columns present, all values in [0, 1], non-increasing in ζ, and ζ = 1 the smallest of the six — with GOS and BCS identical there.
• Eq. 9 is deterministic. gos_bestkappa() is run twice under deliberately different outer seeds and must return the same λ and the same curve to 10⁻¹².
• Scale equivariance and range. Shifting and scaling the response shifts and scales every prediction identically, and predictions never leave the observed range — both necessary consequences of a weighted mean with positive weights.
• The helpers. The base-R VIF matches 1/(1 − R²), and the split generator is reproducible, disjoint, correctly sized, and refuses a degenerate split.

What test-02 pins down

• The packaged zn still holds 894 rows and the screen still retains 885 — so a package update that changed the data would fail here rather than silently rewrite this page.
• All nine correlations, the seven selected covariates (asserted equal to the vignette's set), and their VIFs.
• λ reproduces as 0.08 and the whole 19-point κ–RMSE curve matches to 10⁻⁶.
• All 13,132 grid predictions match to 10⁻⁶ and the ζ = 0.99 uncertainties to 10⁻⁷ — re-predicting the entire grid, because a subset would legitimately differ.
• Three cross-validation repeats reproduce for all three models to 10⁻⁸, the summary table averages its own per-repeat file, and GOS beats BCS and MLR on both metrics in both preprocessing variants.

4 Adaptation, Writing & Reproducibility

4.1 Port to Your Domain

GOS needs no geometry, no neighbourhood and no variogram, so it travels further than most spatial methods: anywhere you have a response measured at some locations, and covariates measured everywhere. Soil properties, air quality, species abundance, house prices, disease incidence — the contract in §2.3 is the whole requirement.

What to edit

1. config/project-config.R — RESPONSE, and LOG_TRANSFORM if your response is right-skewed and strictly positive (check fig01 for your own data before deciding).
2. CANDIDATES — every covariate available at both the samples and the prediction locations. Step 2 whittles the list down; you do not have to pre-select.
3. KAPPA_GRID — keep it fine below 0.1, where the minimum usually sits, and coarse above. If your first run puts λ at the smallest grid value, extend the grid downwards before believing it.
4. CV_REPEATS — 50 is the paper's protocol and costs about a minute here. Lower it while iterating, restore it for the run you report.
5. Nothing else. Rscript run-all.R test, then Rscript run-all.R.

Choosing covariates

The configuration is the model. Three rules that follow from Eqs. 2–4 rather than from taste:

• Availability at the prediction locations is the binding constraint. A covariate you measured only at your samples is useless to GOS, however predictive. This is why remote sensing, terrain and distance-to-feature layers dominate every published configuration.
• Because P is the minimum, more is not better. Every added covariate can only lower each pairwise similarity, never raise it. A weak covariate does not dilute the signal the way it would in a regression — it actively suppresses matches. Let the VIF and correlation screens do their job, and prefer a compact configuration.
• Prefer covariates that carry the spatial structure you cannot otherwise use. The model never sees coordinates, so anything smooth in space that drives your response — a geological unit, a distance to a river, an elevation band — has to enter through a covariate or not at all. In the paper, the two strongest Zn covariates were distances to related lithology and faults, which are exactly such fields; their absence from the packaged data is the main reason this reproduction differs from the published one.
• Watch for circular and categorical variables. Aspect is dropped by the screen here, and rightly: a Gaussian on degrees treats 359° and 1° as far apart. Decompose it into sin/cos, or leave it out. Unordered categories have no meaningful squared difference at all.

What κ means for your sample size

κ is a share, so the number of observations actually used is κ × m — here 0.08 × 885 ≈ 71 per prediction. That has two consequences worth planning for.

4.2 Write the Paper

The published study has a clean architecture, and it is worth copying: define the four stages, characterise the configuration, show the threshold search, then predict, quantify uncertainty and benchmark.

Manuscript skeleton

1  Introduction        — spatial prediction from similarity, not distance;
                          BCS uses every sample; the gap: which samples,
                          and how uncertain (3 aims)
2  Geographically optimal similarity
                        — 2.1 Characterising configurations (Eq. 1)
                          2.2 Assessing similarity (Eqs. 2-6; P = min)
                          2.3 The optimal threshold (Eqs. 7-10)
                          2.4 Prediction and uncertainty (Eqs. 11-12)
3  Case & data         — study area, response and transform, covariates,
                          availability at the prediction locations
4  Results             — 4.1 the selected configuration (correlation, VIF)
                          4.2 the kappa curve and lambda
                          4.3 prediction map
                          4.4 uncertainty across zeta
                          4.5 GOS against BCS: spread and accuracy
                          4.6 accuracy comparison, paired over splits
5  Discussion          — why fewer samples predict better; what the
                          uncertainty does and does not mean; limits of a
                          covariate-only configuration
6  Conclusions         — 6-8 sentences answering the aims

Results-to-manuscript map

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 CRAN geosimilarity package and its bundled data. Reference figures are reproduced from the article, which is open access, for comparison.

1. Song Y (2023). Geographically Optimal Similarity. Mathematical Geosciences 55:295–320. doi:10.1007/s11004-022-10036-8 · PDFopen access (source of the method, of Figs. 1, 3, 4, 7, 8, 9, 10 and 11 reference images, and of the case reproduced here)
2. Song Y, Lyu W. geosimilarity: Geographically Optimal Similarity. R package, version 3.9. CRAN · vignette (the reference implementation, and the source of both data sets)
3. Zhu A-X, Lu G, Liu J, Qin C-Z, Zhou C (2018). Spatial prediction based on Third Law of Geography. Annals of GIS 24(4):225–240. doi:10.1080/19475683.2018.1534890 (the geographical similarity principle GOS builds on)
4. Zhu A-X, Liu J, Du F, et al. (2015). Predictive soil mapping with limited sample data. European Journal of Soil Science 66(3):535–547. doi:10.1111/ejss.12244 (the BCS / iPSM model that κ = 1 reproduces exactly)
5. Online GOS calculator — yongzesong.com/app/gos/ (the same package compiled to WebAssembly, running entirely in the browser)