Reproducing the Degree of Spatial Interpretability (DSI)

A complete walkthrough of the degree of spatial interpretability (DSI) and its local extension, the degree of geocomplexity (DG) — measuring how much of the spatial structure in the data a prediction model actually explains rather than passes through into its residuals — from raw data to manuscript.
Yongze Song

Yongze Song

Associate Professor, Curtin University

Special Issue Editor (SIE), Associate Editor: Int. J. Applied Earth Observation (IF 8.2) · Associate Editor: GIScience & Remote Sensing (IF 6.8)

yongze.song@curtin.edu.au · yongzesong.com

⬇ Download code & example data (zip, 119 KB)
R/ pipeline scripts · config/project-config.R · data/reference/dsi-paper-demo-dataset.csv (the paper's demo data) · tests/ verification · run-all.R — unzip, then run Rscript run-all.R test from the DSI/ folder root.
To cite the DSI and DG models and their R packages and codes in publications, please use:
Liu H, Song Y, Yi W (2026). Degree of spatial interpretability. IJGIS. doi · PDF
Liu H, Song Y, Yi W, Zhang P (2026). Degree of geocomplexity for diagnosing spatial pattern loss beyond prediction accuracy. GIScience & Remote Sensing 63(1):2657087. doi · PDF
Zhang Z, Song Y, Luo P, Wu P (2023). Geocomplexity explains spatial errors. IJGIS 37(7):1449–1469. doi · PDF

1 Method Overview & Reproduction Scope

1.1 Core Idea

Conventional model validation asks how accurate a model is — R², RMSE, MAE. The degree of spatial interpretability (DSI) asks a different question that accuracy metrics cannot answer: how much of the spatial structure present in the data has the model actually explained, rather than passed through into its residuals? A model can lower its errors by smoothing toward the local mean while leaving the spatial arrangement of those errors intact — high accuracy, low spatial interpretability.

The idea rests on a principle from spatial statistics: any spatial structure in the response that a model fails to capture is transferred to its residuals. Measuring a spatial characteristic on the response field and again on the residual field, the relative decrease tells you how much the model absorbed. DSI turns that into a number; its local extension, the degree of geocomplexity (DG), turns it into a map.

Reader anchor

Models in a domain are selected on accuracy → accuracy cannot see whether a model reproduces spatial structure → DSI measures exactly that → applied to a set of models, the two criteria often select different models → therefore selection should weigh both. Everything on this page serves that line.

What you need

R 4.1+ with three required packages (sf, spdep, rpart), the folder from the download button above, and about ten minutes. Optional packages add more models and the DG maps. No internet, no API keys, no paid software; the full demo runs in about 6 seconds.

1.2 Method Logic

DSI is built from three measurements, each applied twice — once to the response field, once to a model's residual field. Two are global (autocorrelation, heterogeneity); one is local (geocomplexity, used by DG). The method composes four earlier building blocks:

Moran's I 1950global spatial autocorrelation of a field
Moran, Biometrikacited
GD 2010geographical detector: Q value of spatial heterogeneity
Wang et al., IJGIScited
GEOC 2023geocomplexity indicator Gi; explains model errors
Zhang et al., IJGIS · PDF
DSI 2026residual-based attenuation η per characteristic, combined into θ
Liu et al., IJGIS · PDF
DG 2026the same attenuation localised: a map of where structure is explained
Liu et al., GISci RS · PDF
Published reference · DSI DSI calculation flowchart from Liu et al. 2026
DSI paper Fig. 1 — the calculation process. Spatial characteristics of the dependent variable (δ₀) and of the modelling residuals (δr) are compared to obtain η for each characteristic, then combined into the θ range. Liu, Song & Yi (2026), IJGIS.

Key concepts and equations

Three equations drive the whole pipeline:

1.3 Reproduction Scope

Everything tagged Demo run on this page regenerates from the scripts and data in this folder; everything tagged Published reference is shown as a static comparison image or table from the two source articles. The demo runs end to end with R 4.6.0, and the built-in test reproduces every published value of the DSI paper's Table 3 (row lm) within 0.005.

ItemDemo templatePublished reference
Data 1,024 synthetic grid points, 6 predictors (data/derived/analysis-data.csv, seed 42) + the paper's own 586-point demo dataset (data/reference/) DSI: 115,083 Australian plant-diversity plots aggregated to 1,958 0.5° grid cells · DG: tree biomass of eastern Australia
Tables results/accuracy-metrics.csv, results/dsi-metrics.csv, results/dg-summary.csv and the sensitivity summaries, regenerated by the scripts DSI Table 3, reproduced value-by-value by run-all.R test and shown beside the demo tables in §3.2–§3.3
Figures fig01–fig11 PDFs and PNGs regenerated by R/p01–p06 (web copies in assets/) DSI Figs 1 and 6 and DG Fig 8, embedded as static reference images
What the template regenerates versus what is shown as published reference.

Published applications of this exact method

ApplicationSpatial unitResponseModelsStatus
Plant-diversity prediction, Australia (Liu, Song & Yi 2026, DSI — reference values on this page) 0.5°×0.5° grid cell (1,958)Species richness 10 regression / machine-learning modelsPublished (IJGIS)
Tree-biomass prediction, eastern Australia (Liu et al. 2026, DG — reference figures on this page) Sample plotAbove-ground biomass Model set scored by local DGPublished (GISci RS)
Synthetic demonstration (this page — all demo outputs) 1,024 grid pointsSimulated response y with planted structure 5 models (9 configured; 4 skipped without their packages)Reproducible demo
Two published applications share the method; the synthetic demo makes it reproducible without external data.

How to use this page

2 Setup, Structure & Data Contract

2.1 Environment & Dependencies

The pipeline is pure R and runs from one entry point. Required packages stop the run if missing; optional packages only disable the model or module that needs them, and the pipeline reports what it skipped — so a run always produces output. Install everything from the R console:

R console
# Required — the pipeline stops without these
install.packages(c("sf", "spdep", "rpart"))

# Optional — the local DG module and the figures
install.packages(c("geocomplexity", "ggplot2", "scales", "patchwork"))

# Optional — the full model set of the demo configuration
install.packages(c("ranger", "xgboost", "Cubist", "FNN",
                   "pls", "earth", "kernlab", "gbm"))
PackageRole in the pipelineRequired?
sf, spdepprojected coordinates, KNN weights, Moran's I (Step 3)yes
rpartregression-tree strata for the Q value (Step 3)yes
geocomplexitylocal geocomplexity Gi for the DG module (Step 4)DG only
ranger, xgboost, Cubist, FNNrf, xgbTree, cubist, knn models (Step 2)per model
pls, earth, kernlab, gbmpls, earth, svmRadial, gbm — the rest of the published model setper model
ggplot2, scales, patchworkfigures (R/p01–p06)figures only
Package roles. Models whose package is missing are skipped with a console message, never a failure — the run below fits 5 of the 9 configured models.
!System libraries for sf

sf needs GDAL, GEOS and PROJ. On macOS: brew install gdal geos proj. On Debian/Ubuntu: apt install libgdal-dev libgeos-dev libproj-dev libudunits2-dev. The complete environment is recorded in env/session-info.txt; package install commands live in env/requirements.md.

2.2 Project Structure & Configuration

One entry script drives everything: run-all.R runs the six pipeline steps, the six figure scripts, or the verification tests, and every step writes to data/derived/, results/, tables/ and figs/. Each step reads what the previous one wrote, so a single step can be re-run after a config change without repeating the whole pipeline. Deleting the four output folders and re-running regenerates everything shown on this page.

folder tree
DSI/
├── dsi.html                       # this guide
├── run-all.R                      # entry point: full run, single steps, figures, test
├── 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 … 05-models.R  # library layer: config, helpers, metrics,
│   │                              #   DSI, DG, model registry
│   ├── 09-make-synthetic-data.R   # built-in demo generator (INPUT_FILE = NULL)
│   ├── 10 … 60                    # pipeline: data, models, DSI, DG,
│   │                              #   sensitivity, tables
│   └── p01 … p06                  # figure scripts
├── data/
│   ├── raw/                       # your input CSV goes here
│   ├── reference/                 # the DSI paper's demo dataset + original code
│   └── derived/                   # generated: analysis data, split, residuals
├── tests/                         # test-01 reproduces the published values;
│                                  #   test-02 checks 17 metric properties
├── env/                           # requirements.md, session-info.txt
├── results/  tables/  figs/       # generated output, never hand-edited

Main configuration file

config/project-config.R (key lines)
INPUT_FILE    <- NULL          # NULL = synthetic demo; else "data/raw/your.csv"
RESPONSE      <- "y"           # name of the dependent-variable column
PREDICTORS    <- NULL          # NULL = every column except lon/lat/response
DATA_CRS      <- 4326          # CRS of the lon/lat columns
PROJECTED_CRS <- 3577          # a projected CRS in METRES (not degrees!)
K_MORAN       <- 8             # KNN neighbours for Moran's I (6-10 dense, 10-20 sparse)
K_GEOC        <- 18            # neighbourhood for geocomplexity (DG)
K_SENSITIVITY <- c(6, 8, 10, 14, 18, 22)   # k values re-run in Step 5
STRATIFY_METHOD <- "tree"      # "tree" (paper) or "kmeans" — never the response
MODELS <- c("rpart","rf","xgbTree","cubist","knn",
            "pls","earth","svmRadial","gbm")
TRAIN_FRACTION <- 0.7          # 70/30 split, as in both papers
SEED           <- 42
Single adaptation point

The config is the only file to edit when porting the pipeline to a new domain — no R script hard-codes variable names, coordinates, or parameters. When writing the manuscript, quote parameter values from this file, never from memory.

Key parameters and where they come from

ParameterDefaultOrigin / rationale
K_MORAN8KNN weights for Moran's I; the DSI paper used k = 8 for its 0.5° Australian grid; sensitivity in Step 5
K_GEOC18geocomplexity neighbourhood of the DG paper; must be large enough that neighbours share neighbours
STRATIFY_METHODtreea regression tree on the predictors defines the Q-value strata, reproducing the DSI paper; "kmeans" is the Step 5 alternative
TRAIN_FRACTION, SEED0.7, 42stratified 70/30 split; residuals are always computed on the held-out 30%
MODELS9 namesthe demo model set; entries without their package are skipped with a message
Free parameters, their defaults, and the published rules they follow.

2.3 Data Contract

The pipeline consumes one table, one row per sampling location. Everything except the coordinates and the response is treated as a predictor — used both to fit models and to define the Q-value strata. Units can be anything point-referenced: grid cells (DSI application), sample plots (DG application), stations, polygon centroids.

ColumnConfig keyTypeMeaningRequired?
lon, latDATA_CRSnumericcoordinates, in the CRS you declare (demo: EPSG:4326)yes
yRESPONSEnumericthe dependent variableyes
predictorsPREDICTORS / DROP_COLSnumericevery remaining column unless listed in DROP_COLSyes
Required schema of the input table. Column names are free — they are declared in the config, not hard-coded.

Example analysis-ready table

lonlaty temperatureprecipitationelevation soil_phdistance_kmnoise_var
114.00−39.0012.51224.956555.82260.166.576120.60−1.825
115.29−39.0011.21823.490593.68303.256.518113.751.420
116.58−39.0010.98723.692541.58366.646.818176.090.724
First rows of data/derived/analysis-data.csv, produced by Step 1. 1,024 rows total; your own data takes exactly this shape — coordinates, a response, and predictor columns.
!Predictors define the strata, so choose them deliberately

The Q statistic partitions space using the predictors. An identifier column left in by accident (a plot ID, a date) will be treated as a covariate and will partition space — put such columns in DROP_COLS. Dropping a predictor for collinearity also changes the strata, so make that decision before the analysis and report it in the paper.

The generator behind the demo (or: how to build a controlled test)

When INPUT_FILE is NULL, the pipeline builds a synthetic field with known structure — a smooth autocorrelated component, a three-region regime term, and a spatially correlated error. Because the truth is known, a model that recovers it should score high DSI and a model that ignores space should not. The same generator is the starting point for the simulation section of a methods paper.

R/09-make-synthetic-data.R (excerpt)
# Smooth spatial field from an exponential covariance (autocorrelation)
smooth_field <- function(range_km, sd = 1) {
  covm <- sd^2 * exp(-d_km / range_km)
  L <- chol(covm + diag(1e-6, n))
  as.numeric(t(L) %*% rnorm(n))
}
# Regime term: three east-west bands (heterogeneity -> picked up by Q)
regime_effect <- c(-4, 0, 5)[cut(grid$lon, breaks = 3, labels = FALSE)]
# Spatially correlated error, so a good aspatial fit still leaves structure
spatial_error <- 2.2 * smooth_field(range_km * 0.5, 1)

y <- 12 + 0.9*(temperature - mean(temperature)) + 1.5*(soil_ph - mean(soil_ph)) +
     regime_effect + spatial_error + rnorm(n, sd = 1)

3 Reproduction Pipeline, Results & Validation

3.1 Pipeline Execution

Six pipeline steps, six figure scripts and two verification tests run from one entry point. Each step writes to data/derived/, results/, tables/ or figs/ and reads what the previous step wrote, so any step can be re-run in isolation after a config change.

Full run
Rscript run-all.R (~6 s on demo data)
Verification
Rscript run-all.R test — do this first
Output folders
results/ tables/ figs/ data/derived/
terminal
cd DSI                    # this folder
Rscript run-all.R test    # 1. verify the implementation reproduces the paper
Rscript run-all.R         # 2. full pipeline on the built-in synthetic data

Rscript run-all.R 30 40   # re-run single steps (10 20 30 40 50 60)
Rscript run-all.R figures # only the figure scripts
Project : dsi-demo Domain : synthetic demonstration == Step 1/6 Prepare data ====================================== INPUT_FILE is NULL; generating the synthetic demonstration dataset analysis set: 1024 points, 6 predictors == Step 2/6 Fit models and extract residuals ================== stratified split: 714 training / 310 test (70% / 30%) fitting 5 model(s): rpart, rf, xgbTree, cubist, knn best R2: rf (0.366) | lowest RMSE: rf (4.122) == Step 3/6 Compute DSI ======================================= Spearman rho between R2 and theta_probable: 0.800 accuracy and interpretability disagree: rf leads on R2, xgbTree on theta_probable == Step 4/6 Compute DG ======================================== locally best model: xgbTree at 32.6% of points (101 of 310) == Step 5/6 Sensitivity of DSI to measurement settings ======== DSI is stable across neighbourhood sizes (all changes below 5%). == Step 6/6 Build manuscript tables =========================== Finished in 6.3 s Outputs: results/ tables/ figs/ data/derived/

3.2 Reproduce the Published Values

Before trusting any number on your own data, confirm the code reproduces the published method. run-all.R test runs the template's metric functions on the demo dataset released with the DSI paper (data/reference/dsi-paper-demo-dataset.csv, 586 test points) and compares nine values against Table 3 (row lm) and Figure 6 of Liu, Song & Yi (2026).

== test-01 Reproduce Liu et al. (2026) Table 3, row lm ================== quantity published computed difference result ---------------------------------------------------------- delta_a_0 0.409 0.410 0.0012 pass delta_h_0 0.513 0.513 0.0001 pass theta_max 0.817 0.817 0.0003 pass ... PASS — the template reproduces every published value within 0.005 == test-02 Metric properties =========================================== pass a residual field with no structure gives eta = 1 pass an amplified residual field gives eta < 0 ... 17 passed, 0 failed
QuantityPublishedComputedΔResult
δa₀ — Moran's I of y0.4090.4100.0012pass
δar — Moran's I of residuals0.1890.1900.0007pass
δh₀ — Q of y0.5130.5130.0001pass
δhr — Q of residuals0.2030.2030.0003pass
ηa0.5380.5380.0005pass
ηh0.6040.6050.0009pass
θmin0.5380.5380.0005pass
θprobable0.6040.6050.0009pass
θmax0.8170.8170.0003pass
Published vs. recomputed — the DSI paper's demo dataset, linear model, 586 test points. All nine values agree within 0.005; the 0.001-level gaps are rounding in the printed table.
!The Q statistic has two forms; only one reproduces the paper

The geographical detector's Q can be written with within-stratum sample variances or as one minus the ratio of within-stratum to total sum of squares. On the demo response these give 0.490 vs 0.513. Only the sum-of-squares form matches the published δh values, so that is what q_statistic() implements. The reproduction test locks it: change the form and the test fails.

R/02-spatial-metrics.R — the Q statistic (sum-of-squares form)
q_statistic <- function(x, strata) {
  strata <- as.factor(strata)
  keep <- is.finite(x) & !is.na(strata)
  x <- x[keep]; strata <- droplevels(strata[keep])
  sst <- sum((x - mean(x))^2)               # total sum of squares
  if (sst <= 0) return(NA_real_)
  ssw <- sum(vapply(split(x, strata),       # within-stratum sum of squares
                    function(v) sum((v - mean(v))^2), numeric(1)))
  1 - ssw / sst
}
Checkpoint

run-all.R test prints PASS with 17 passed, 0 failed. The implementation is faithful to the published method — proceed to your own analysis. If the test fails after you edit the metric code, the edit broke the method.

3.3 Core Analytical Steps

Five stages, each with the same anatomy: purpose → code → output → result → how to read it. Steps are numbered 1–5 on this page while the script files keep their own numbers (Step 1 runs R/10-prepare-data.R, Step 3 runs R/30-compute-dsi.R, and so on) — R/00–05 is the shared library layer.

Step 1

Prepare data — raw input → analysis-ready table

What this step does

Reads the input (or generates the synthetic demo when INPUT_FILE is NULL), drops rows with missing coordinates/response and duplicate coordinates, optionally log-transforms the response, reports predictor pairs with |r| > 0.8 without acting on them, and writes data/derived/analysis-data.csv. Everything downstream reads only the derived file, so the raw data is never modified.

Code

terminal
Rscript run-all.R 10

Example output

== Step 1/6 Prepare data ====================================== INPUT_FILE is NULL; generating the synthetic demonstration dataset 1024 rows, 9 columns response: y 6 predictor(s): temperature, precipitation, elevation, soil_ph, distance_km, noise_var no predictor pair exceeds |r| = 0.8 wrote results/variable-summary.csv (7 rows) wrote data/derived/analysis-data.csv (1024 rows) analysis set: 1024 points, 6 predictors
Demo run Synthetic study area, response field
fig01a — the response y across 1,024 grid points. The east-west colour gradient is the regime term; the local patchiness is the autocorrelated component — both are structure a model should reproduce.

Paper reference

In Liu, Song & Yi (2026) this step aggregates 115,083 Australian plant-diversity plots to 1,958 0.5°×0.5° grid cells (richness as the response, 19 abiotic predictors — PDF). The synthetic demo mirrors that structure at smaller scale so the pipeline runs without downloading the case data.

How to read this result
  • Confirm there is structure to explain. DSI is an attenuation ratio; if the response carries no spatial pattern, every η is noise over noise. A visibly patchy response map (like this one) is the precondition for the whole analysis.
  • Two kinds of structure, deliberately. The gradient (heterogeneity) is what the Q statistic detects; the local clustering (autocorrelation) is what Moran's I detects. A case that shows only one makes only one η interpretable.
  • Collinearity is reported, not fixed. The step flags |r| > 0.8 pairs so you decide; because dropping a predictor reshapes the Q-value strata, that decision belongs in the manuscript, not in a silent default.
Checkpoint

data/derived/analysis-data.csv exists with 1,024 rows and 6 predictors, and the console reports no predictor pair above |r| = 0.8, before you move on.

Step 2

Fit models & accuracy — predictions, residuals, the conventional ranking

What this step does

Splits the data 70/30 by quantiles of the response, fits every available model from the configured set, and computes residuals on the held-out test set. Accuracy (R², RMSE, MAE) is the conventional ranking that DSI will be contrasted against.

!Residuals must come from held-out data

Training residuals of an overfitted model approach zero, which drives η toward a spurious 1. Both the paper and this template evaluate on the test set; the split is stratified so the test set keeps the response distribution and therefore a comparable δ₀.

Code

terminal
Rscript run-all.R 20

Example output

== Step 2/6 Fit models and extract residuals ================== stratified split: 714 training / 310 test (70% / 30%) ! skipped (package not installed): pls [pls], earth [earth], svmRadial [kernlab], gbm [gbm] fitting 5 model(s): rpart, rf, xgbTree, cubist, knn wrote data/derived/test-residuals.csv (310 rows) wrote results/accuracy-metrics.csv (5 rows) best R2: rf (0.366) | lowest RMSE: rf (4.122)
ModelRMSEMAE
rf0.3664.1223.383
xgbTree0.3184.4073.530
knn0.3094.3203.557
cubist0.3094.3553.479
rpart0.2174.7793.879
results/accuracy-metrics.csv — 310 test points, ranked by R². The conventional view; on its own it would send you to rf.
Demo run Observed vs predicted per model
fig02 — observed vs. predicted on the test set, per model.

Paper reference

In the DSI paper's Australian case, rf and cubist achieved the best accuracy (R² ≈ 0.48–0.49), with svmRadial, gbm, gam, xgbTree close behind and earth, lm, knn, pls weakest — a ranking the DSI analysis then partly overturns (PDF).

How to read this result
  • This is the baseline you will challenge, not the answer. The accuracy table establishes the conventional hierarchy; the entire point of DSI is to show that hierarchy is incomplete.
  • Small R² gaps hide large behavioural differences. Four models sit within 0.31–0.37 R² — a reviewer would call them "comparable." Step 3 shows their θmax ranges from 0.47 to 0.74, a difference accuracy cannot see.
  • Deeper thinking: the modest R² (≈ 0.37) is expected when predictors are purely abiotic and a spatially correlated error is present by construction — exactly the regime where spatially structured residuals, and therefore DSI, matter most.
Checkpoint

results/accuracy-metrics.csv holds 5 rows with rf best at R² = 0.366; models with missing packages are listed as skipped, not failed.

Step 3 · the core result

Compute DSI — degree of spatial interpretability per model

What this step does

Under one fixed set-up — the same neighbours, strata, and evaluation points for every model — it measures Moran's I and the Q statistic on the response and on each residual field, forms η for each characteristic (Eq. 1), and combines them into the θ range (Eq. 2).

Code

R/03-dsi.R — the indicator (called by step 30)
dsi_eta <- function(delta_0, delta_r) {
  if (abs(delta_0) < 1e-8) return(NA_real_)     # nothing to explain
  (delta_0 - delta_r) / delta_0
}
dsi_theta <- function(etas) {
  etas <- etas[!is.na(etas)]
  list(theta_min      = min(etas),
       theta_probable = max(etas),
       theta_max      = 1 - prod(1 - etas))
}
terminal
Rscript run-all.R 30

Example output

== Step 3/6 Compute DSI ======================================= Moran's I weights: 8-nearest neighbours, row-standardised, EPSG:3577 Q-value strata: tree method, 20 strata wrote results/dsi-metrics.csv (5 rows) No diagnostic flags raised. wrote results/accuracy-vs-dsi.csv (5 rows) Spearman rho between R2 and theta_probable: 0.800 accuracy and interpretability disagree: rf leads on R2, xgbTree on theta_probable
Modelδarηa δhrηh θminθprobθmax
xgbTree0.4260.3370.2080.6010.3370.6010.735
rpart0.5160.1970.3010.4230.1970.4230.537
knn0.5380.1630.2960.4330.1630.4330.526
rf0.5610.1260.2900.4450.1260.4450.515
cubist0.5800.0970.3050.4160.0970.4160.473
results/dsi-metrics.csv — δa₀ = 0.643, δh₀ = 0.522 (shared by all models), ranked by θmax.
Demo run DSI theta range per model, demo
fig03 — the θ triangle per model, built exactly like the paper's Fig 6: base from θ₁ (θmin) to θ₃ (θmax), apex and dashed line at θ₂ (θprobable), best model on the top row. xgbTree leads.
Published reference · DSI DSI theta distribution, paper Fig 6
DSI paper Fig. 6 — the same figure in Liu et al. (2026), where cubist leads (θ 0.799–0.987) — a direct visual template for the demo panel.

Paper reference — Table 3 (all models)

Modelδarηa δhrηh
cubist0.0270.9340.1030.799
rf0.0390.9050.1270.752
xgbTree0.0480.8830.0950.815
svmRadial0.0590.8560.1570.694
gbm0.0780.8090.1220.762
gam0.0950.7680.1330.741
earth0.1420.6530.1640.680
knn0.1430.6500.1860.637
lm0.1890.5380.2030.604
pls0.2720.3350.2740.466
Liu, Song & Yi (2026) Table 3 — δa₀ = 0.409, δh₀ = 0.513 (shared). The lm row is the one run-all.R test reproduces (PDF).
How to read this result
  • η is an attenuation, not an error. ηa = 0.337 for xgbTree means it removed 34% of the response's autocorrelation from its residuals; ηh = 0.601 means it removed 60% of the heterogeneity. Higher is better.
  • Negative η is a finding, not a bug — watch for it. η < 0 means residuals more spatially clustered than the response: the model manufactures pattern, the typical signature of a global linear or additive fit that cannot bend across a regime boundary. The pipeline flags any negative η in results/dsi-diagnostics.txt (none occurs in this five-model demo).
  • Read the two components, not just θ. A model can reach a respectable θmax from one strong component and one weak one. xgbTree here is far better on heterogeneity (0.60) than autocorrelation (0.34) — the wide θ span (0.34→0.74) says so at a glance.
  • Demo vs paper — same logic, different winner. In the paper's real data every tree/ensemble model reaches η > 0.75 and cubist leads; on the harder synthetic case the absolute η are lower and xgbTree leads. The ranking mechanism reproduces even when the numbers differ, which is what a method port must demonstrate.
  • Deeper thinking: that rf (best accuracy) sits mid-table on θ is the payload — it interpolates well within dense clusters yet leaves neighbourhood-scale residual structure. Accuracy rewarded the interpolation; DSI penalised the leftover structure.
Checkpoint

results/dsi-metrics.csv holds 5 rows; the console reports "No diagnostic flags raised" and ρ(R², θprobable) = 0.800 with rf and xgbTree as the divergent pair.

Step 4 · optional module

Compute DG — point-wise geocomplexity & the locally optimal model

What this step does

Computes local geocomplexity of the response and of each residual field at every test location, forms DG point-by-point (Eq. 3), summarises the distribution, and assigns each location the model with the highest local DG — producing a map that argues no single model wins everywhere. Needs the geocomplexity package; set RUN_DG <- FALSE to skip.

Code

terminal
Rscript run-all.R 40

Example output

== Step 4/6 Compute DG ======================================== geocomplexity neighbourhood: 18 nearest neighbours wrote results/dg-pointwise.csv (1550 rows) wrote results/dg-summary.csv (5 rows) wrote results/dg-optimal-model.csv (310 rows) locally best model: xgbTree at 32.6% of points (101 of 310) 5 model(s) are locally optimal somewhere
ModelMean DGMedian DGShare DG > 0
rpart31.550.51263.2%
xgbTree22.360.50464.8%
cubist29.800.40762.6%
knn16.880.40261.9%
rf28.450.34959.0%
results/dg-summary.csv — ranked by median DG; note the heavy-tailed means.
Demo run Locally optimal model map, demo
fig08 — the model with the highest local DG at each of 310 test points. xgbTree is locally best at 33% of points, but all five models win somewhere — no single model dominates the map.
Published reference · DG Locally optimal model map, DG paper Fig 8
DG paper Fig. 8 — the same construction in Liu et al. (2026) for tree-biomass in eastern Australia: rf best at 23.6% of locations, knn 19.3%, xgbTree 16.6%, each with clear geographic niches.

Paper reference

The DG paper reports mean DG values from −0.57 (lm) to 1.61 (rf) with an unbounded, heavy-tailed distribution — the reason this template reports the median beside the mean (PDF).

!Quote the median, not the mean

DG has local complexity in its denominator, so points where complexity is near zero produce extreme ratios. The demo means (17–32) are dominated by a handful of such points; the medians (0.35–0.51) are the honest summary. The pipeline warns when the mean and median disagree in sign.

How to read this result
  • DG localises what DSI totals. DSI said xgbTree explains the most structure overall; DG shows it does so at most — but not all — locations, and names the second-best model elsewhere.
  • "No single best model" is the practical claim. The map is the argument for local model selection: in a real study you might deploy different models in different regions, an option accuracy metrics never surface.
  • Rank ≠ magnitude. A model can be locally best at many points while explaining little in absolute terms — always read the map together with the summary table before calling a model good.
  • Deeper thinking: with a dense, near-regular synthetic grid the niches are diffuse; in the paper's clustered biomass plots they are sharp and geographic (coastal vs inland). Sparse or clustered sampling makes the map more decisive — and, below a few multiples of K_GEOC, noisier.
Checkpoint

results/dg-summary.csv holds 5 rows with medians between 0.349 and 0.512, and results/dg-optimal-model.csv assigns all 310 test points; the console reports xgbTree locally best at 32.6%.

Step 5

Sensitivity — is the result an artefact of the settings?

What this step does

Two settings are analyst choices, not properties of the data: the neighbourhood size k behind Moran's I, and the rule that defines the Q-value strata. This step recomputes DSI across k ∈ {6, 8, 10, 14, 18, 22} and under an alternative partition (k-means instead of tree), so the manuscript can state how much the ranking moves.

Code

terminal
Rscript run-all.R 50

Example output

== Step 5/6 Sensitivity of DSI to measurement settings ======== wrote results/sensitivity-k.csv (30 rows) k varied over {6, 8, 10, 14, 18, 22}; largest change in theta_probable is 0.00% (cubist at k = 6) DSI is stable across neighbourhood sizes (all changes below 5%). wrote results/sensitivity-strata.csv (5 rows) model ranking agreement across stratification methods: rho = -0.900 wrote results/sensitivity-summary.txt
Demo run theta_probable across neighbourhood sizes
fig10 — θprobable across k ∈ {6,8,10,14,18,22}. The lines are flat — the ranking does not depend on the neighbourhood size (largest change < 5%, the paper's robustness bar).
Demo run eta_heterogeneity under two stratifications
fig11 — ηh under tree vs k-means strata. Points off the diagonal show the heterogeneity component does shift with the partition (rank agreement ρ = −0.90).

Paper reference

Section 5.2 of Liu, Song & Yi (2026) compares standard vs. spatial block cross-validation for rf and finds all DSI components move by less than 5% (Δθprobable = 3.4%, Δθmin = 4.8%, Δθmax = 1.3%) — the basis for its robustness claim (PDF).

How to read this result
  • Neighbourhood size is safe here; report it anyway. θprobable is invariant to k for models where heterogeneity dominates (ηh does not use k), which is why the lines are flat. State the chosen k and its stability so a reviewer cannot suspect cherry-picking.
  • Stratification is not neutral — this is the honest caveat. Tree and k-means strata rank the heterogeneity component differently (ρ = −0.90). Name the stratification rule in the Methods section and, when rankings differ, report both. Hiding this is the kind of fragility reviewers probe.
  • Deeper thinking: the disagreement is larger on the demo than it would be in a strongly heterogeneous real case, because the synthetic regime term is only one of several signals. The more the response is genuinely regime-driven, the more any sensible partition agrees.
Checkpoint

results/sensitivity-summary.txt states both findings: stable across k (all changes below 5%) and partition-dependent across stratification rules (ρ = −0.900, report both rankings).

3.4 Validation & Interpretation

The whole project converges on one figure: interpretability against accuracy. Everything else is scaffolding for the moment where the two criteria are placed on the same axes and disagree.

Demo run DSI vs accuracy, demo
fig05 — θ (vertical span θmin–θmax, point at θprobable) against R², RMSE and MAE. Models high-left are accurate and interpretable; models low-right buy accuracy while leaving spatial structure in their residuals.
Published reference · DSI DSI paper Fig 7 — DSI against R2, RMSE and MAE per model
DSI paper Fig. 7 — the published counterpart of demo fig05: DSI(θ) with its θmin–θmax span against R², RMSE and MAE for the ten Australian plant-diversity models. The annotations carry the paper's headline: highest DSI is cubist, lowest RMSE is rf — the two criteria pick different winners (PDF).
Published reference · DG DG paper Fig 10 — mean DG with 95% CI against accuracy metrics
DG paper Fig. 10 — the same dual-criterion reading at the local scale: mean DG with 95% confidence intervals against R², RMSE and MAE across ten models. Models ranked highly by accuracy do not automatically rank highly on spatial interpretability — the trend is positive but far from tight (PDF).
The demo divergence
rf leads on accuracy (R² 0.366) yet sits fourth of five on θmax (0.515). xgbTree, second on accuracy, leads interpretability (θprobable 0.601). Rank correlation ρ(R², θprobable) = 0.80 — correlated, not identical.
The paper divergence
In Liu, Song & Yi (2026), rf beats cubist on fit (R² 0.487 vs 0.478) yet cubist has the higher DSI; knn out-fits lm yet DSI separates them the other way. Same lesson, real data.
How to interpret the divergence — several angles
  • Mechanism (why it happens). Minimising squared error rewards smoothing toward the local mean, which lowers RMSE while leaving the neighbourhood-scale arrangement of errors intact. Accuracy scores the magnitude of errors; DSI scores their spatial organisation. The two can move apart whenever error magnitude and error structure come apart.
  • Statistical grounding. The divergence is not a quirk of the data — it follows from residual diagnostics: structure the model does not capture is transferred to its residuals (Anselin 1988; Zhang et al. 2023). A low RMSE achieved by averaging still fails this test.
  • Decision angle (what to do about it). If the deliverable is a single aggregate number, accuracy is the right criterion. If the deliverable is a map — a regional allocation, a targeted intervention — then a model that reproduces spatial structure is worth a small accuracy cost, and DSI is how you find it.
  • The negative-η angle. Models below zero on θmin — none in this five-model demo, but the usual suspects are global linear and additive fits — are the strongest evidence that accuracy alone misleads: competitive on R² yet actively distorting spatial pattern. Report them prominently when they occur — they are the paper's best argument.
  • What DSI is not. DSI does not replace accuracy and is generally correlated with it (demo: ρ = 0.80). Its job is to catch the cases where the two part ways — model selection should read both axes, not swap one for the other.

Validation summary

Validation itemFile / figureDecision rule
Published-value reproductionrun-all.R test (test-01)all nine Table-3 values within 0.005 — the code is the published method
Metric propertiestests/test-02-metric-properties.R17 property checks pass (η bounds, θ ordering, Q invariances, Moran significance)
Neighbourhood sensitivityresults/sensitivity-k.csv, fig10θ changes below 5% across k ∈ {6…22} (demo: 0.00%)
Stratification sensitivityresults/sensitivity-strata.csv, fig11rank agreement reported (demo: ρ = −0.90 → state the rule, report both rankings)
Four decision rules that must hold — or be honestly caveated — before the results are written up.
Reading order for your own run

1) results/dsi-diagnostics.txt — what must be caveated (non-significant δ₀, negative η, wide θ). 2) results/accuracy-dsi-disagreement.txt — the rank correlation and the divergent model pair. 3) results/sensitivity-summary.txt — whether the ranking survives the settings. 4) the figures above.

Main finding narrative

The demo recovers the published mechanism on fresh data: accuracy and interpretability correlate (ρ = 0.80) but pick different winners — rf on R², xgbTree on θprobable — and the best-fitting model sits only fourth of five on θmax. DG then shows the winner is local, not global: xgbTree is best at only 33% of locations and all five models win somewhere. The validation quartet shows none of this is an artefact of the sample or the settings — except the stratification rule, whose effect is quantified (ρ = −0.90) and must be reported.

4 Adaptation, Writing & Reproducibility

4.1 Port to Your Domain

The pipeline was designed to move between domains without touching the analytical scripts: the config declares the data, the preparation step builds the analysis table, and everything downstream is generic. Two published applications illustrate the span:

DomainResponseSpatial unitKey adaptation
Plant-diversity prediction (Liu, Song & Yi 2026) Species richness0.5° grid cell 19 abiotic predictors; Australian Albers CRS (EPSG:3577); k = 8
Tree-biomass prediction (Liu et al. 2026, DG) Above-ground biomassSample plot log-transformed response; geocomplexity neighbourhood k = 18
Any point-referenced unit with one numeric response and a set of continuous predictors fits the data contract of §2.3.

What to edit

  1. config/project-config.R — point INPUT_FILE at your CSV, name the RESPONSE column, set PROJECTED_CRS to a metric CRS for your region, and revisit K_MORAN (6–10 dense samples, 10–20 sparse). The single adaptation point of §2.2.
  2. Optionally list DROP_COLS (IDs, dates) and set LOG_TRANSFORM_RESPONSE <- TRUE for right-skewed responses (biomass, income).
  3. Nothing else: run Rscript run-all.R test once (the published-value check does not depend on your data), then Rscript run-all.R, and read results/dsi-diagnostics.txt before drafting anything.
!Common pitfalls

PROJECTED_CRS must be in metres — degrees silently distort every KNN neighbourhood and Moran's I with it. Q-value strata come from the predictors, never the response: quantile bins of the response drive its own Q towards 1 and make η meaningless. Residuals are computed on held-out data — training residuals of an overfitted model produce a spurious η near 1. The Q statistic uses the sum-of-squares form; the reproduction test locks it, so never edit the metric code without re-running run-all.R test. And for DG, quote the median — the mean is heavy-tailed by construction.

4.2 Write the Paper

The manuscript mirrors the pipeline: Methods sections describe the indicator and the experiment design 1:1 with the scripts, every Results sentence carries a number pulled from tables/*.tex or results/ so nothing is retyped, and the Discussion is built on the accuracy–interpretability divergence and the sensitivity sweeps. Section 2 (the indicator) is largely written already, because its definition does not change between domains — the writing work concentrates in the case study, the Results, and the mechanism of the divergence.

Manuscript sectionTemplate outputWriting job
Methodsconfig/project-config.Rdefine the three spatial characteristics of the response, their residual attenuation ratios (η), and the θ aggregation, with the model list and every parameter quoted from the config, never from memory
Resultstables/table-accuracy.tex, table-dsi.tex, table-theta.tex, table-dg.tex + figs/fig02–09order accuracy → DSI → divergence → local (DG) → sensitivity, mirroring the pipeline steps 1:1; every finding sentence carries a number; report δ₀ before any η and caveat non-significant characteristics
Discussionresults/sensitivity-*.csv + figs/fig10–11explain why accuracy and interpretability diverge (error minimisation rewards smoothing toward the local mean while leaving neighbourhood-scale structure intact); contribution, limitations
Where each manuscript section draws its material from. Fill the abstract from the tables first — accuracy ranking from table-accuracy, the θ range and leader from table-theta, the divergence case from table-accuracy-vs-dsi — because the abstract fixes the argument every other section must serve.

Manuscript skeleton

manuscript outline
1  Introduction            — 6 paragraphs, narrowing: background → gap → aim   (~1,200 words)
2  The indicator           — 2.1 Concept (no formulas)                          (~1,400 words;
                             2.2 Calculation process (ALL formulas live here):   largely
                                 2.2.1 Spatial characteristics of the response   pre-written)
                                       (autocorrelation, heterogeneity,
                                       geocomplexity; significance δ0)
                                 2.2.2 Residual attenuation of each
                                       characteristic (η)
                                 2.2.3 θ aggregation (θ_min, θ_probable, θ_max)
                                 2.2.4 Local degree of interpretability (DG)
3  Case study              — study area, data, experiment design — no results  (~1,300 words)
4  Results                 — accuracy → DSI → divergence → local (DG) →        (~1,600 words)
                             sensitivity, mirroring the pipeline steps 1:1;
                             every finding sentence carries a number;
                             one interpretive takeaway per subsection
5  Discussion              — why the divergence happens, contribution,         (~1,200 words)
                             limitations; no subsection headings
6  Conclusions             — 6-7 sentences, each mapping to a contribution     (~300 words)

Results-to-manuscript map

Manuscript elementSource fileShown in
Study-area / response figurefigs/fig01a–c§3.1 of paper
Accuracy table + scattertables/table-accuracy.tex, figs/fig02§4.1
Main table (DSI per model)tables/table-dsi.tex, tables/table-theta.tex§4.2 (≈ published Table 3)
θ-range figurefigs/fig03, figs/fig04§4.2
Divergence figurefigs/fig05, tables/table-accuracy-vs-dsi.tex§4.3
DG maps + optimal-model mapfigs/fig06–09, tables/table-dg.tex§4.4
Robustness numbersresults/sensitivity-*.csv, figs/fig10–11§4.5 / response to reviewers
Every manuscript element traces to one regenerable file.

4.3 Final Reproducibility Package

Everything an editor, reviewer, or reader needs sits in this folder; the three blocks below are what a submission's data-and-code availability statement should point to.

Code
numbered scripts in R/ + run-all.R + config/project-config.R + tests/
Data
data/reference/dsi-paper-demo-dataset.csv (published demo) — or your raw data in data/raw/
Results
figs/fig01–fig11 (PDF+PNG) + tables/*.tex as the manuscript tables + env/session-info.txt
Before claiming reproducible
  • Rscript run-all.R test passes — the code still reproduces the published values (change the metric code and this is the tripwire).
  • Delete results/, tables/, figs/ and data/derived/, re-run Rscript run-all.R — every number and figure on this page regenerates (8.1 s).
  • CRS, k, stratification rule, split fraction and seed are recorded in the manuscript and match config/project-config.R.
  • Every finding sentence traces to a file in results/ or tables/ — none retyped; δ₀ is reported before any η.
  • env/session-info.txt is included in the deposit.
Anticipate the reviewer

Have the answer pre-computed: Why this k? → design justification + sensitivity-summary.txt. Would a different stratification change the ranking?sensitivity-strata.csv, both rankings reported. Is DG's mean trustworthy given the heavy tail? → the median is already in dg-summary.csv. Is the test set spatially contaminated? → the DSI paper's block-CV check moved all components by less than 5%; add block CV for your case or state the limitation.

References & credits

Demo outputs were generated by the scripts in this folder. Reference figures and table values are reproduced from the cited articles for comparison and are © their authors and publishers. The upstream reference implementation and demo dataset are credited to the DSI authors.

  1. Moran PAP (1950). Notes on continuous stochastic phenomena. Biometrika 37(1–2):17–23. doi:10.1093/biomet/37.1-2.17
  2. Wang J-F, Li X-H, Christakos G, et al. (2010). Geographical detectors-based health risk assessment… IJGIS 24(1):107–127. doi:10.1080/13658810802443457
  3. Anselin L (1988). Spatial Econometrics: Methods and Models. Springer. doi:10.1007/978-94-015-7799-1
  4. Zhang Z, Song Y, Luo P, Wu P (2023). Geocomplexity explains spatial errors. IJGIS 37(7):1449–1469. doi:10.1080/13658816.2023.2203212 · PDF
  5. Liu H, Song Y, Yi W (2026). Degree of spatial interpretability. IJGIS. doi:10.1080/13658816.2026.2614335 · PDF (source of Table 3 and the Fig 1/Fig 6 reference images on this page)
  6. Liu H, Song Y, Yi W, Zhang P (2026). Degree of geocomplexity for diagnosing spatial pattern loss beyond prediction accuracy. GIScience & Remote Sensing 63(1):2657087. doi:10.1080/15481603.2026.2657087 · PDF (source of the DG Fig 8 reference image on this page)
  7. Reference implementation and demo dataset: github.com/KevinOceanLiu/DSI