Reproducing the Spatial Outliers Explaining Heterogeneity

1 Method Overview & Reproduction Scope

1.1 Core Idea

Most spatial heterogeneity models read a covariate at the analysis unit and ask how much of the response it explains. That reads only the first dimension. The SOH model adds a second dimension: the configuration of local anomalies at nearby locations where the response was never observed. Two units with identical air temperature can sit in entirely different neighbourhoods, one beside an anomalously cold pocket and one in a uniform surround, and SOH claims that difference explains variation the first dimension cannot reach.

The model turns those anomalies into multi-scale explanatory variables called spatial outlier patterns (SOPs), stratifies the analysis units with a regression tree, and scores the stratification with the geographical detector power of determinant (PD).

• Research problem: covariate values at the analysis unit ignore the spatial configuration of anomalies in the surroundings, so conventional stratified-heterogeneity models leave second-dimension variation unexplained.
• One-line contribution: SOH converts local anomalies at unsampled grid cells into positive and negative SOP variables across many neighbourhood radii, and shows they out-explain the covariates they were derived from.
• Why it matters: covariates that look weak at the unit (soil carbon in the barley case: PD 0.032) can be among the strongest drivers once their anomaly pattern is read (SOP.SOC: 0.415) — a discovery the first dimension cannot make.

1.2 Method Logic

Three papers form one line, and each supplies a piece the next one uses; the geographical detector supplies the scoring statistic underneath all of them.

Key concepts and equations

Four building blocks drive the whole pipeline, one from each reference:

• Stratified power (PD / q-statistic): q = 1 − SSW/SST over the strata of a partition, with a non-central F significance test; reimplemented natively in code/R/fn_pd.R — Wang et al. (2010) doi
• The second dimension: information at locations outside the samples carries spatial association, generated from a searching range plus an operator over the surrounding field — Song (2022) doi · PDF
• Spatial outlier patterns (Equations 3–4): grid values are standardised locally within each neighbourhood N_r(u), and z-scores beyond ±s local standard deviations are accumulated into SOP⁺ and SOP⁻ — outliers treated as signal, not noise — Ren et al. (2025) doi · PDF
• The SOH model (Equations 5–6): SOP blocks become explanatory variables, a CART tree stratifies the units, PD scores every configuration, and a conditional permutation null tests whether the second dimension adds anything beyond the first — Ren et al. (2026) doi · PDF

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 Ren et al. (2026). The demo runs end to end with R 4.6.0 — all 13 result tables regenerate byte-identically to the copies shipped in expected-results/tables/.

Published application of this exact pipeline

How to use this page

• To reproduce the demo: this folder already contains everything — the R scripts (code/), the configuration (config/), and the synthetic example data (data/demo/). Run Rscript code/99_run_all.R from the folder root; every expected output is shown on this page so you can verify each step.
• To apply SOH to your own problem: read §2.3 (the two-table data contract), edit the §2.2 config and the ADAPT blocks of the preparation script, then run the same scripts unchanged (§4.1).
• To write the paper: §4.2 maps every result table to its manuscript element and settles the application-versus-extension decision first.

2 Setup, Structure & Data Contract

2.1 Environment & Dependencies

The template runs on a bare R installation. Only rpart and ggplot2 are required; the rest add features and degrade gracefully when absent. The power-of-determinant q statistic is implemented natively, so the GD package is optional. The pipeline was run with R 4.6.0; a full run records its own session-info.txt under 03_results/logs/.

# minimum to run the pipeline
install.packages(c("rpart", "ggplot2"))

# optional, each unlocks one feature
install.packages(c(
  "patchwork",  # multi-panel figures (Moran, augmentation)
  "spdep",      # global Moran's I of the response
  "sf",         # read shapefile study areas in your own data prep
  "GD"          # cross-check the native q statistic against the published implementation
))

2.2 Project Structure & Configuration

Everything study-specific lives in one file, config/config.R. No path, variable name, buffer radius, or threshold belongs anywhere else. The numbered scripts run in order and each reads from and writes to data/interim/, so any step can be re-run on its own after a configuration change.

SOH/
├── soh.html                     # this guide
├── assets/                      # figures shown in this guide (demo + paper)
├── config/config.R              # every study-specific setting, in one file
├── data/
│   ├── demo/                    # synthetic example tables (points.csv, grids.csv)
│   ├── interim/ processed/      # created on first run / your own data
├── code/
│   ├── 00_setup.R               # config, libraries, function loading
│   ├── 01_prepare_data.R        # assemble and validate the two tables
│   ├── 02_generate_sops.R       # build the second dimension
│   ├── 03_pd_individual.R       # Moran's I and single-factor PD
│   ├── 04_pd_interactions.R     # SOP-SOP, covariate-SOP, categories
│   ├── 05_validation.R          # gains, scenarios, permutation nulls
│   ├── 06_sensitivity.R         # scale and threshold sweeps
│   ├── 07_figures.R             # all figures
│   ├── 99_run_all.R             # the whole pipeline
│   ├── R/                       # fn_utils, fn_sop, fn_pd, fn_experiments, fn_plots
│   └── gee/                     # Google Earth Engine covariate extraction
├── expected-results/tables/     # T01–T13 of the run
└── 03_results/                  # tables | figures | logs (auto-generated on run)

Main configuration file

To run the built-in synthetic demo used throughout this page, leave data$mode as "demo"; to run your own study, set it to "user" and give the column names of your two tables.

# config/config.R  (excerpt — the full file documents every setting, with the
# values used by the published barley study recorded beside each one)
cfg <- list(
  study = list(id = "demo", response_label = "Response (unit)"),
  data  = list(
    mode = "demo",                 # "demo" | "user"
    crs_epsg = 3577,               # projected CRS with metric units
    dist_unit = "km"               # unit of x/y AND of the buffers below
  ),
  variables = list(
    codes  = c("V1","V2","V3","V4","V5","V6"),
    category = c(V1="C1", V2="C1", V3="C2", V4="C2", V5="C3", V6="C3")
  ),
  sop = list(
    buffers = seq(20, 200, 20),    # neighbourhood radii, in dist_unit
    sd_threshold = 2,              # outlier threshold, in local SD
    statistic = "zsum"             # accumulate z-scores beyond the threshold
  ),
  pd = list(cart_cp = 0.01)        # rpart complexity parameter
)

2.3 Data Contract

The SOH model needs exactly two tables on a common projected coordinate system with metric units.

points — one row per analysis unit, the support of the response. This is what the model explains. Columns: id, x, y, response, and one column per covariate.

grids — one row per fine grid cell covering the study area and a margin beyond it. These unsampled cells are the second dimension. Same columns as points, minus response.

3 Reproduction Pipeline, Results & Validation

3.1 Pipeline Execution

The whole pipeline runs end to end in about fifteen seconds on the demo. Every step reads its input from data/interim/ and writes its output there, so any single step can be re-run on its own after a configuration change. Step numbers on this page match the script numbers: Step 1 runs code/01_prepare_data.R, and so on.

Rscript code/99_run_all.R

# 01 prepare  -> T01, validation
# 02 sops     -> T02, sopvars.rds
# 03 indiv    -> T03, T04, F03, F05
# 04 interact -> T05-T08, F06
# 05 validate -> T09-T11, F07, F08
# 06 sens     -> T12, T13, F09, F10
# 07 figures  -> all figures written to 03_results/figures/

3.2 Core Analytical Steps

Four analytical stages build the model (Steps 1–4), and two validation stages test it (Steps 5–6, in §3.3). Each stage has the same anatomy: purpose → script → result → how to read it, with the demo output on the left of every figure pair and the published counterpart from Ren et al. (2026) on the right.

Step 1

Data preparation

What this step does

Build the two tables (synthetic in demo mode, your ingestion code in user mode), then run validation checks that catch the failures which silently produce meaningless PD values: unprojected coordinates, a grid coarser than the units, a missing margin, constant covariates.

Code

terminal Copy

Rscript code/01_prepare_data.R
# demo mode injects known local anomalies into half the covariates, so the
# whole pipeline runs before you have any data of your own.

Example output

01 PREPARE DATA mode: demo — generating a synthetic dataset with known anomalies demo tables written to data/demo/ validation passed table written: 03_results/tables/T01_data_summary.csv (7 rows) analysis units: 300 grid cells: 3600 covariates: 6 cached: data/interim/data.rds

Demo run

Demo run

◆How to read this result

• (what the demo builds) A smooth autocorrelated covariate field (the first dimension) plus compact injected anomalies (the second dimension), and a response that depends on both. Because the anomaly locations are known, the demo doubles as a correctness check.
• (the knob that matters) cfg$demo$outlier_signal_weight controls how much the response depends on the second dimension. Set it to 0 and the second dimension carries no information — the falsification case the method must be unable to claim (see Step 5).
• (what to check first) The buffer diagnostic printed here reports whether your largest radius sits at the recommended 10–20% of the maximum pairwise distance between units. Too large and SOPs stop being local; too small and neighbourhoods hold too few cells for a standard deviation to mean anything.

✓Checkpoint

T01_data_summary.csv written; the console reports 300 analysis units, 3,600 grid cells, 6 covariates, and "validation passed" before you move on.

Step 2

Build the second dimension (SOP)

What this step does

For each covariate and each neighbourhood radius, standardise the grid values locally and accumulate the anomalies. Positive and negative patterns are kept separate throughout, because a neighbourhood surrounded by unusually high values and one surrounded by unusually low values are different geographic situations.

◆Equations 3–4

z(x_v) = ( x_v − mean_Nr(u) ) / sd_Nr(u)
SOP⁺(X,r,u) = Σ z(x_v) for v in N_r(u) with z(x_v) > s
SOP⁻(X,r,u) = Σ z(x_v) for v in N_r(u) with z(x_v) < −s

Each covariate yields 2 × (number of radii) SOP columns per unit. The standardisation is local to each neighbourhood, not global — that is what makes an SOP a local anomaly rather than a global extreme. s = 2 reproduces the published model.

Code

terminal Copy

Rscript code/02_generate_sops.R

# under the hood (code/R/fn_sop.R):
sopvars <- soh_generate_sop_all(points, grids, cfg)
# -> named list, one data.frame per covariate, each with 2 * length(buffers) columns
#    named p_b20, p_b40, ... (positive) and n_b20, n_b40, ... (negative)

Example output

A diagnostic table (T02) reports, per covariate, the share of empty SOP cells. A column that is entirely zero carries no information, usually because the buffer is too small to hold min_cells grid cells, or the threshold is too strict for a smooth covariate.

variable	n_columns	zero_share	mean_positive	mean_negative
V1	20	0.33	20.3	−44.9
V2	20	0.33	38.1	−26.6
V3	20	0.43	38.3	−39.4
V4	20	0.59	37.9	−16.7

Demo run

Published reference · SOH

Same construct, same behaviour: at small radii almost nothing exceeds two local standard deviations; as the radius grows, coherent positive and negative regions emerge and sharpen. That growth is the signal the scale-sensitivity step in Step 6 quantifies.

Paper reference

Ren et al. (2026) Fig. 6 (PDF); the SOP construction follows Ren et al. (2025) (PDF) and the second-dimension idea of Song (2022) (PDF).

◆How to read this result

• (paper reading) The SOP maps show that anomaly structure is scale-dependent: the same covariate is locally unremarkable at 20 km but strongly patterned at 200 km. Positive and negative regions are spatially organised, not scattered — evidence that they encode real geography rather than noise.
• (mechanism) An SOP is high where a unit sits in a neighbourhood whose covariate field departs sharply from its local average. It reads the shape of the surroundings, which a single covariate value at the unit cannot express.
• (case reasoning) In barley, the positive SOP of soil pH concentrates where alkaline pockets interrupt otherwise acidic cropland; in a health study, expect the positive SOP of PM2.5 to light up around point sources embedded in cleaner surroundings. The map is a hypothesis generator for where second-dimension effects will matter.
• (what to check) If your SOP maps are blank at every radius, the covariate is too smooth to have a second dimension. Lower sd_threshold, widen the smallest buffer, or accept it as a null result rather than tuning until something appears.

✓Checkpoint

T02_sop_diagnostics.csv shows zero_share 0.33–0.59 on the demo (no covariate above 0.9), and data/interim/sopvars.rds is cached.

Step 3

Individual power of determinant

What this step does

Before scoring the SOPs, establish that the response is spatially structured at all, then compare each covariate against its own SOP block. Each covariate and each SOP block is stratified by a CART tree and scored by PD (Equation 6); the core SOH claim is that the SOP block reaches a higher PD than the covariate it was derived from.

Code

terminal + R Copy

Rscript code/03_pd_individual.R
# global Moran's I of the response with 999 permutations (spdep)

individual <- soh_exp_individual(points, sopvars, cfg)
#   label     type      qv     sig       n_strata
#   SOP.V2    SOP     0.614   3.1e-10       19
#   SOP.V1    SOP     0.556   7.8e-10       16
#   ...
#   V1        variable 0.263  9.5e-10       10   <- strongest raw covariate

Spatial autocorrelation of the response

Demo run

Published reference · SOH

Both responses are strongly, significantly autocorrelated. This is a prerequisite the SOH claim rests on, and also the caveat attached to every PD value below.

Covariate versus its SOP block

Demo run

Published reference · SOH

The rank shift is the headline of this step. Soil variables that are weak on their own (SOC 0.032) become among the strongest once their spatial anomaly pattern is read (SOP.SOC 0.415).

Paper reference

Ren et al. (2026) Figs. 5 and 7 (PDF); the PD statistic follows Wang et al. (2010).

◆How to read this result

• (paper reading) Each SOP block outscoring its covariate is the SOH claim in its simplest form: the configuration of local anomalies explains more heterogeneity than the covariate's own value. In Ren et al. (2026) SOP.pHc (0.467) more than doubles pHc (0.205), and SOP.SOC (0.415) beats SOC (0.032) more than tenfold.
• (mechanism) A high SOP-block PD means units that share a neighbourhood anomaly signature also share response levels — spatial dependence acting through the surroundings, not the site.
• (case reasoning) Rank the SOP-minus-covariate gap per variable. The largest gaps flag drivers that act through spatial context. In barley the soil variables gain most, consistent with patchy soil chemistry that a plot-level average washes out; in an urban port expect amenity and accessibility variables to gain most, because what a neighbourhood is near matters more than what it contains.
• (what to check — and hold your claim) A high SOP-only PD is not yet evidence. SOPs summarise the covariate field around a unit, so they correlate with the covariate at the unit by construction, and there are many of them. Withhold the second-dimension claim until the conditional null in Step 5 clears. Also confirm n_strata is well below the number of units, so the tree is stratifying rather than memorising.

✓Checkpoint

T03_moran.csv reports Moran's I = 0.585 (p = 0.001); T04_individual_pd.csv is led by SOP.V2 = 0.614 against 0.263 for the strongest raw covariate on the demo data.

Step 4

Interactions

What this step does

PD is computed for every pair of factors and classified in geographical-detector terms — whether two factors together explain more than the sum of their parts (enhance, nonlinear), more than the stronger alone (enhance, bivariate), or less (weaken). Interactions run at two levels: individual SOP pairs and thematic category blocks.

Code

terminal + R Copy

Rscript code/04_pd_interactions.R
pairs_ss  <- soh_exp_pairs(points, sopvars, individual, cfg, mode = "sop_sop")
pairs_vs  <- soh_exp_pairs(points, sopvars, individual, cfg, mode = "var_sop")
categories <- soh_exp_categories(points, sopvars, cfg)
#   label            qv     q1     q2     interaction
#   SOP.V1 & SOP.V2  0.695  0.556  0.614  enhance, bivariate

Example output

Demo run

Published reference · SOH

The published matrices are denser (12 covariates, four categories) but tell the same story: interactions lift PD above every individual factor, and the lift is strongest across thematic boundaries.

Interaction (Ren et al. 2026)	PD
SOP.C1 × SOP.C3 (climate × soil)	0.686
SOP.C3 × SOP.C4 (soil × vegetation)	0.666
SOP.C2 × SOP.C3 (topography × soil)	0.663
Variable.C1 × SOP.C3 (climate value × soil anomaly)	0.72

Paper reference

Ren et al. (2026) Fig. 8 (PDF); the interaction typing follows the geographical detector of Wang et al. (2010).

◆How to read this result

• (paper reading) The strongest interactions in Ren et al. (2026) all involve the soil-pH and soil-carbon SOPs crossed with a climate factor. Climate sets the broad envelope of barley productivity; soil anomalies explain the departures within it, so the two together partition the response far better than either alone (SOP.C1 × SOP.C3 = 0.686 against 0.467 for the best single factor).
• (mechanism) An 'enhance, bivariate' or 'enhance, nonlinear' label means the two factors carve the study area along different spatial seams. Stacking them creates finer, more homogeneous strata than either seam alone.
• (case reasoning) Read the interaction table as a map of which driver pairs co-govern the outcome. Cross-category interactions (climate × soil) are more interesting than within-category ones, because they point to a coupled process rather than two facets of one. In a health port, a pollution-SOP × density-SOP interaction would flag neighbourhoods where an emission anomaly and a crowding anomaly compound.
• (what to check) The interaction-type distribution (T08) is a result in itself. A model where every pair merely 'enhances bivariate' behaves differently from one full of 'enhance nonlinear' pairs — the latter signals genuine synergy that no additive model would find.

✓Checkpoint

T05_pairs_sop_sop.csv is led by SOP.V1 & SOP.V2 = 0.695 ('enhance, bivariate') on the demo data; T08_interaction_types.csv holds the type distribution.

3.3 Results, Validation & Interpretation

Three checks in increasing strength close the pipeline: how much the second dimension adds per variable, which of the three overall models wins, and — decisively — whether the gain survives a null that controls for both dimensionality and the SOP–covariate correlation (Step 5); a scale and threshold sweep then shows the conclusions do not hinge on the two free parameters (Step 6).

Step 5

Validation — the headline

What this step does

Quantify the gain from the second dimension per covariate (augmentation), compare the three overall scenarios (A1 SOPs only, A2 both dimensions, A3 covariates only), and run the two permutation nulls. The conditional null is an addition to the source study and the single most important part of the reproduction.

Code

terminal + R Copy

Rscript code/05_validation.R
aug <- soh_exp_augmentation(points, sopvars, cfg)   # per variable and per category

Augmentation: gain from adding the SOP

Demo run

Published reference · SOH

The demo and the paper produce the same shape — a uniform, sizeable positive gain across all variables — which is exactly what a working SOH reproduction should look like.

The three overall scenarios

Demo run

Published reference · SOH

In both, adding SOPs to the covariates raises overall PD, and SOPs on their own are competitive with the entire conventional covariate set.

The two permutation nulls — the decisive test

Two mechanisms inflate the SOP-only PD regardless of whether local anomalies carry any information: sheer dimensionality (hundreds of SOP columns), and the fact that an SOP correlates with the covariate field it summarises. Only a conditional null controls for both.

Test	observed	null mean	empirical p	verdict (demo)
Marginal — SOP-only vs shuffled SOP	0.699	0.400	0.02	clears (necessary, weak)
Conditional — covariates+SOP vs covariates+shuffled SOP	0.700	0.612	0.02	clears (decisive)
reference: covariates only (A3)	0.586	—	—	—

!Why the marginal null is not enough — the falsification test

Set cfg$demo$outlier_signal_weight = 0 so the response depends only on the smooth covariates, and re-run. The SOP-only PD still reaches 0.745 and still clears the marginal null — yet covariates alone score 0.825, and the conditional null fails at p = 1.0. A study that reported that 0.745 as evidence for the second dimension would be reporting an artefact. Always report the conditional null.

Paper reference

Ren et al. (2026) Figs. 9–10 (PDF); the permutation-null layer is an addition of this template.

◆How to read this result

• (paper reading) Ren et al. (2026) report the combined model as best (0.716) and interpret the SOP gain as evidence that outlier patterns explain heterogeneity beyond the covariates. The augmentation figure supports this variable by variable.
• (the deeper point) A high SOP-only PD proves nothing on its own, because SOPs are built from the covariate field and there are hundreds of them. The claim lives or dies on the conditional null: real SOPs must beat shuffled SOPs once the covariates are already in the model.
• (case reasoning) When you port SOH, run the falsification test first on your own data — build a response you know depends only on covariates and confirm the conditional null fails. If it passes on a known-null response, your pipeline has a leak (often an edge-effect or an id misalignment between points and grids).
• (what to report) Put the conditional null in the body of the paper, not the supplement. A reviewer who has to hunt for it will assume it was buried, and it is the one number that separates a real second-dimension effect from SOPs quietly proxying the covariates.

✓Checkpoint

T11_permutation_null.csv shows the conditional null clearing at p = 0.02 (observed 0.700 vs null mean 0.612) with the signal on — and failing at p = 1.0 in the signal-off falsification case.

Step 6

Scale and threshold sensitivity

What this step does

Sweep the neighbourhood radius and the outlier threshold. The headline methodological finding of the reference study came from the scale sweep: PD of the SOP-only and combined models rose with radius and stabilised near 200 km, while the covariate-only model was flat. That plateau is the scale at which local anomaly structure stops adding context.

Code

terminal + R Copy

Rscript code/06_sensitivity.R
scale_df <- soh_exp_scale(points, sopvars, cfg)        # PD vs largest radius
thr_df   <- soh_exp_threshold(points, grids, cfg)      # PD vs 1.5 / 2 / 2.5 SD

Example output

Demo run

Published reference · SOH

Identical qualitative behaviour. The plateau radius is the scale finding to report — the distance beyond which widening the neighbourhood adds no explanatory context.

Paper reference

Ren et al. (2026) Fig. 10b (PDF).

◆How to read this result

• (paper reading) The ~200 km plateau means barley's second dimension is a broad-scale phenomenon: local anomalies matter up to a regional extent, beyond which the neighbourhood becomes so wide that it converges on the global average and SOPs stop being local.
• (mechanism) Read the plateau on the SOP-only curve, not the combined one. The combined curve is anchored by the covariate floor and would report a plateau it inherited rather than earned.
• (case reasoning) The plateau radius is a domain fingerprint. A commuting-scale outcome (housing, crashes) will plateau at a few kilometres; a climate-driven outcome (crops, fire) at hundreds. If your PD is still rising at the largest radius in the sweep, extend cfg$scale_sweep before claiming any threshold — the finding is the location of the plateau, and you cannot report one you have not reached.
• (what to check) Because CART is greedy, the scale curve need not rise monotonically and the combined model can dip below SOP-only. That is expected behaviour, not a bug. Interpret the plateau, not each step; the threshold sweep confirms the conclusions are stable to the one free parameter that a reviewer will question.

✓Checkpoint

T12_scale_sensitivity.csv shows the plateau inside the sweep range; T13_threshold_sensitivity.csv keeps the combined PD within a 0.032 spread across thresholds 1.5 / 2 / 2.5 SD on the demo data.

code/07_figures.R regenerates all ten figures (F01–F10) from the cached interim data alone — no recomputation — into 03_results/figures/; the web copies shown on this page live in assets/demo_figures/.

Validation summary

4 Adaptation, Writing & Reproducibility

4.1 Port to Your Domain

A candidate problem needs four things. The binding one is covariates at a resolution finer than your analysis units — without that, there is no second dimension to extract.

• A response on discrete analysis units — polygons or points. Continuous fields must be aggregated to units first (a modifiable areal unit problem to acknowledge).
• Covariates finer than the units — remote sensing, gridded climate, high-resolution soil, dense sensors. Anything published only at unit level carries no second dimension.
• Genuine local anomalies — hotspots, pockets, scars, clusters. A perfectly smooth covariate has an empty second dimension.
• Room around the units — covariate coverage must extend past the study boundary by at least the largest buffer.

What to edit

1. Prepare data/processed/points.csv and grids.csv (use code/gee/ if your covariates come from Earth Engine).
2. In config/config.R: set data$mode = "user", the column names, a projected crs_epsg, the variable codes and their category, and the buffers.
3. Write your ingestion code in the two blocks of code/01_prepare_data.R marked ADAPT, then run code/99_run_all.R and work through the §4.3 checklist.

4.2 Write the Paper

The published application shares one manuscript architecture: Methods sections mirror the pipeline steps 1:1, every Results sentence carries a number from a 03_results/tables/ CSV, and the Discussion is built on the conditional null and the scale plateau of Steps 5–6. Decide first whether the paper is an application (SOH used unchanged; the contribution is the domain finding, and the target is the domain's journal — most ports) or a methodological extension (the model itself changed; Section 2 carries the full formalisation plus a simulation study, and the target is IJGIS or Spatial Statistics) — application papers that pose as methods papers get rejected for thin novelty.

Manuscript skeleton (shared with the published application)

1  Introduction            — 6 paragraphs; P4 reviews GD/SDA/SDO; P5 gap sentence
                             "However, there are still challenges in {ISSUE} due to
                             {GAP_1}, {GAP_2}, {GAP_3}."; P6 aim + 3 contributions
2  {MODEL NAME}            — 2.1 Concept (no formulas)
                             2.2 Calculation process (ALL formulas live here):
                                 2.2.1 Neighbourhood buffers and local standardisation
                                 2.2.2 Spatial outlier patterns SOP+ / SOP-   (Eqs 3-4)
                                 2.2.3 CART stratification                    (Eq 5)
                                 2.2.4 Power of determinant + significance    (Eq 6)
3  Case                    — 3.1 Study area & data · 3.2 Pre-processing
                             3.3 Model-based analysis (1:1 with §2.2 steps)
                             3.4 Validation: augmentation · scenarios · permutation
                                 nulls · scale and threshold sensitivity
4  Results                 — subsections mirror 3.3 steps 1:1; every finding sentence
                             carries a number; one interpretive takeaway per subsection
5  Discussion              — 4 paragraphs, 2 tables (capability comparison;
                             reproduction guidance), no subsection headings
6  Conclusions             — 8-10 sentences answering: accomplished? contribution?
                             meaning for readers? limitations & future work?

Results-to-manuscript map

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.

References & credits

Demo outputs were generated by the scripts in this folder. Reference figures are reproduced from Ren et al. (2026) for comparison; cite the method papers, not this walkthrough.

1. Ren K, Song Y, Yang X, Wang X, Chen M, Yu Q (2026). Spatial outliers as a pattern determinant for explaining heterogeneity. International Journal of Geographical Information Science. doi:10.1080/13658816.2026.2682957 · PDF (source of all reference figures on this page)
2. Ren K, Song Y, Yu Q (2025). Second-dimension outliers for spatial prediction. IJGIS 40(6):1915–1942. doi:10.1080/13658816.2025.2580414 · PDF
3. Song Y (2022). The second dimension of spatial association. International Journal of Applied Earth Observation and Geoinformation 111:102834. doi:10.1016/j.jag.2022.102834 · PDF
4. Wang J-F, et al. (2010). Geographical detectors-based health risk assessment. IJGIS 24(1):107–127. doi:10.1080/13658810802443457

Authors' code and data: github.com/renkaigis/Second-dimension_Outlier-driven_Heterogeneity (MIT) and Figshare doi:10.6084/m9.figshare.31169701.