Reproducing the Optimal Parameters-based Geographical Detector (OPGD)

1 Method Overview & Reproduction Scope

1.1 Core Idea

The geographical detector measures spatial stratified heterogeneity (SSH): the tendency for a response to be more alike within strata of an explanatory variable than between them. Its core statistic, the Q value (q = 1 − SSW/SST), reads as the share of the response's variance that a variable's spatial strata explain — 0 when the strata carry no signal, 1 when they explain it entirely. Unlike a regression coefficient it needs no linearity assumption and works on categorical and stratified continuous variables alike.

The catch is that a continuous variable has no natural strata: it must be discretised first, and the Q value depends heavily on how — which breaking method (equal, natural, quantile, geometric, standard-deviation) and how many intervals. Choose badly and a genuinely important variable looks weak. The optimal parameters-based geographical detector (OPGD) removes that arbitrariness: for every continuous variable it searches the grid of (method × break-number) and keeps the combination that maximises the Q value, then runs the detectors on those optimal strata.

• Research problem: the plain geographical detector's result is an artefact of the analyst's discretisation choice, so two studies of the same data can disagree on which driver matters.
• One-line contribution: OPGD selects the discretisation (and the spatial unit size) that best expresses each variable's spatial stratified heterogeneity, making the four detectors — factor, interaction, risk and ecological — reproducible and comparable.
• Why it matters: a single optimised run answers four questions at once — which variables drive the response, which pairs interact, which zones are high-risk, and whether two drivers differ significantly — with no model tuning and no linearity assumption.

1.2 Method Logic

OPGD wraps the classical geographical detector (Wang et al. 2010) in two parameter optimisations — spatial discretisation and spatial scale — and then exposes the four detectors. It sits between the original detector and the localised methods that later build on it:

The schematic below is the whole method on one page:

Published reference · OPGD

Key concepts and equations

• Q value — the factor detector (Wang et al. 2010):
  q = 1 − SSW / SST = 1 − (1/Nσ²) Σ_h N_hσ_h²
  SST is the total variance of the response; SSW is the variance that remains within the strata of a variable. q ∈ [0,1] is the share the strata explain, with an F/non-central test of significance — Wang et al. (2010) doi.
• Optimal discretisation (the OPGD enhancement): for a continuous variable, compute q for every (method, n-intervals) pair over the methods {equal, natural, quantile, geometric, sd} and a break-number range, and keep the pair with the highest q. The variable is then cut at those breaks before any detector runs — Song et al. (2020) PDF.
• Optimal spatial scale (SESU): re-run the model across a series of spatial unit sizes and adopt the scale where the Q values stabilise, so the finding is not an artefact of the grid resolution.
• The four detectors, all off the same optimal strata: the factor detector ranks variables by q; the interaction detector compares q(x₁∩x₂) against q(x₁) and q(x₂) and types the pair into five classes; the risk detector t-tests the response means between strata to name high/low zones; the ecological detector F-tests whether two variables differ significantly in their impact.

The whole analysis in one call

The GD package's gdm() function does everything above in a single call — optimal discretisation of the continuous variables, then all four detectors, plus the plots. The block below is the vignette's quick-start (GD vignette, §1.3); paste it into R to run the entire NDVI analysis at once. The rest of this page unpacks that one call into transparent, config-driven steps that write one CSV per detector.

## install and library the package
install.packages("GD", dep = TRUE)
library(GD)

## Example 1 — NDVI: ndvi_40
## optional discretisation methods: equal, natural, quantile, geometric, sd, manual
discmethod <- c("equal", "natural", "quantile")
discitv    <- c(4:6)

## "gdm" runs optimal discretisation + all four detectors + plots in one step.
## Here Climatezone and Mining are categorical; Tempchange and GDP are continuous.
ndvigdm <- gdm(NDVIchange ~ Climatezone + Mining + Tempchange + GDP,
               continuous_variable = c("Tempchange", "GDP"),
               data = ndvi_40,
               discmethod = discmethod, discitv = discitv)   # ~3 s
ndvigdm        # printed summary: factor, interaction, risk, ecological
plot(ndvigdm)  # all detector figures

1.3 Reproduction Scope

Everything tagged Demo run on this page regenerates from the scripts in this folder; everything tagged Published reference is a static comparison image from the source article. The demo runs end to end with R 4.6.0 and GD 10.9, and the built-in test reproduces the published NDVI ranking (Song et al. 2020): Precipitation and Climate zone are the dominant drivers of vegetation change.

Published applications of the OPGD model

How to use this page

• To reproduce the demo: this folder ships the R scripts (R/), the configuration (config/project-config.R) and the tests (tests/); the data live inside the GD package. Run Rscript run-all.R test, then Rscript run-all.R; every output on this page regenerates.
• To apply OPGD to your own problem: read §2.3 (data contract), point INPUT_FILE at your CSV and name the response, categorical and continuous columns in the §2.2 config, then run the same scripts (§4.1).
• To write the paper: §4.2 maps every result file to a manuscript element and gives the section skeleton the published applications share.

2 Setup, Structure & Data Contract

2.1 Environment & Dependencies

The pipeline is pure R and wraps the CRAN GD package — the reference implementation of the geographical detector and OPGD. Two packages cover everything; the demo datasets ship inside GD itself, so a fresh machine needs no external download.

# Required — the OPGD reference implementation (ships the demo data)
install.packages("GD")

# Required for the figures
install.packages("ggplot2")

2.2 Project Structure & Configuration

One entry script drives everything: run-all.R runs the seven pipeline steps, the five figure scripts, or the two tests, and every step writes to data/derived/, results/, tables/ or 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.

OPGD/
├── opgd.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, 01-helpers.R  # library layer: paths, IO, GD guard
│   ├── 10-prepare-data.R          # load the analysis table, report roles
│   ├── 20-optimal-discretization.R# the OPGD enhancement (optidisc)
│   ├── 30-factor-detector.R       # relative importance (gd)
│   ├── 40-interaction-detector.R  # pairwise interaction (gdinteract)
│   ├── 50-risk-detector.R         # sub-region means + t-test (gdrisk, riskmean)
│   ├── 60-ecological-detector.R   # F-test between variables (gdeco)
│   ├── 70-sesu-scale.R            # spatial-scale effects (gdm over ndvi_*)
│   ├── 90-tables.R                # LaTeX tables + session info
│   └── p01 … p05                  # figure scripts
├── data/reference/                # provenance of the built-in GD data
├── tests/                         # test-01 reproduces the paper; test-02 properties
├── env/                           # requirements.md, session-info.txt
├── results/  tables/  figs/       # generated output, never hand-edited

Main configuration file

DATASET     <- "ndvi_40"                     # a GD built-in table (or set INPUT_FILE)
INPUT_FILE  <- NULL                          # e.g. "data/raw/your.csv" to use your own data

RESPONSE    <- "NDVIchange"                   # the dependent variable
CATEGORICAL <- c("Climatezone", "Mining")    # already-discrete explanatory vars
CONTINUOUS  <- c("Tempchange", "Precipitation",
                 "GDP", "Popdensity")        # these get optimal discretisation

DISC_METHODS   <- c("equal","natural","quantile","geometric","sd")
DISC_INTERVALS <- 3:7                        # break numbers to search

SESU_DATASETS <- c("ndvi_20","ndvi_30","ndvi_40","ndvi_50")  # spatial-scale series
SESU_SIZES    <- c(20, 30, 40, 50)

Key parameters and where they come from

2.3 Data Contract

The pipeline consumes one table, one row per spatial unit: a numeric response, a set of categorical explanatory variables, and a set of continuous ones. Units can be anything area-referenced — grid cells (this demo), road segments, counties, regions.

Example analysis-ready table

3 Reproduction Pipeline, Results & Validation

3.1 Pipeline Execution

Seven pipeline steps, five figure scripts and two verification tests run from one entry point. 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 config change.

cd OPGD                    # this folder
Rscript run-all.R test     # 1. verify the code reproduces the published ranking
Rscript run-all.R          # 2. full pipeline on the built-in ndvi_40 data

Rscript run-all.R 20 30    # re-run single steps (10 20 30 40 50 60 70 90)
Rscript run-all.R figures  # only the figure scripts

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 OPGD optimisation and factor detector on the paper's own ndvi_40 data and checks the published ranking — Precipitation and Climate zone are the dominant drivers of NDVI change (Song et al. 2020).

3.3 Core Analytical Steps

Seven stages, each with the same anatomy: purpose → code → output → result → how to read it. Steps are numbered 1–7 on this page and the script files share those numbers (Step 2 runs R/20-optimal-discretization.R, and so on).

Rscript run-all.R 10

Step 2 · the OPGD enhancement

Optimal discretisation — the best strata for each continuous variable

What this step does

For every continuous variable it computes the factor-detector Q across the grid of five methods × break numbers 3–7, keeps the combination with the highest Q, and cuts the variable at those breaks. This is the step that makes the whole analysis reproducible — the discretised table it writes feeds every detector.

Code

R/20-optimal-discretization.R — the search (via GD::optidisc) Copy

for (v in CONTINUOUS) {
  od <- optidisc(as.formula(paste(RESPONSE, "~", v)), data = d,
                 discmethod = DISC_METHODS, discitv = DISC_INTERVALS)[[1]]
  bestq <- max(od$qv.matrix, na.rm = TRUE)   # highest Q over method x n-intervals
  d[[v]] <- cut(d[[v]], unique(od$itv), include.lowest = TRUE)  # optimal strata
}

terminal Copy

Rscript run-all.R 20

Example output

Variable	Best method	Intervals	Q value
Precipitation	equal	7	0.869
Tempchange	sd	7	0.327
Popdensity	quantile	7	0.196
GDP	quantile	7	0.099

Demo run

Published reference · OPGD

!The choice is not neutral — report it

Precipitation's Q ranges from 0.37 (geometric, 3 intervals) to 0.87 (equal, 7) across the grid — a factor of more than two, from the binning alone. State the search grid and the chosen method/breaks in the Methods section; that is what makes a geographical-detector result comparable across studies.

◆How to read this result

• The optimal method differs per variable. Precipitation prefers equal breaks, temperature standard-deviation breaks, the skewed GDP and population quantile breaks — one global rule would have mis-scored at least one of them.
• Higher break numbers usually help, but not always. Most curves rise toward 7 intervals; where a curve dips (temperature under equal breaks at n=4) the search simply avoids it. Cap the range by sample size — too many strata leave too few observations each.
• A flat, low curve is a real finding. GDP never exceeds Q≈0.10 under any binning: it is genuinely a weak driver of NDVI change here, not a victim of bad discretisation.

✓Checkpoint

results/optimal-discretization.csv holds one row per continuous variable, and data/derived/discretized-data.csv (713 rows) carries the optimal strata for Steps 3–6.

Step 3 · the core result

Factor detector — relative importance of every variable

What this step does

Runs the factor detector on the optimally discretised table: the Q value of every variable, with a significance test, ranked. This is the headline of most geographical-detector papers — which drivers explain the spatial pattern of the response.

Code

R/30-factor-detector.R — via GD::gd Copy

dcz <- read.csv("data/derived/discretized-data.csv")   # optimal strata from Step 2
g   <- gd(NDVIchange ~ ., data = dcz)                   # factor detector
g$Factor                                                # variable, qv, sig

terminal Copy

Rscript run-all.R 30

Example output

Demo run

Variable	Q value	Significance (p)	Rank
Precipitation	0.869	2.6e−10	1
Climatezone	0.822	7.3e−10	2
Tempchange	0.327	3.1e−10	3
Popdensity	0.196	2.2e−10	4
Mining	0.141	6.7e−10	5
GDP	0.099	9.6e−10	6

Paper reference

Song et al. (2020) report Precipitation and Climate zone as the leading drivers of NDVI change in this region, with the socioeconomic variables (GDP, population) weak — the ranking this demo reproduces (PDF).

◆How to read this result

• Q is a share of variance explained, not a correlation. Precipitation's 0.869 means its (optimal) strata account for 87% of the spatial stratified heterogeneity of NDVI change — no sign, no linearity, just explained structure.
• Significance and magnitude are different questions. Every variable here is statistically significant, yet GDP's Q is only 0.10. A significant but small Q means "real but minor"; read the number, not just the star.
• Climate leads, economy trails — the physical reading. That precipitation and climate zone tower over GDP and population says vegetation change in this arid region is climate-driven; the socioeconomic drivers matter only locally (they resurface in the interaction detector, Step 4).

✓Checkpoint

results/factor-detector.csv holds 6 ranked rows; the console names Precipitation as dominant (Q = 0.869) and GDP as weakest (0.099).

Step 4

Interaction detector — which pairs of variables amplify each other

What this step does

For every pair it computes the joint Q of the two variables overlaid and compares it with the two single Q values, classifying the pair into one of five types (bi-enhance, nonlinear-enhance, uni-weaken, nonlinear-weaken, independent). Interactions routinely exceed the strongest single driver.

Code

terminal Copy

Rscript run-all.R 40

Example output

Demo run

Published reference · OPGD

Published reference · road case

◆How to read this result

• Every pair enhances — the general rule of this method. The joint Q of two variables is almost always higher than either alone, so a pair explains more than a sum of parts would suggest. The question is how, not whether.
• Bi-enhance vs nonlinear-enhance is the key distinction. 13 pairs are bi-enhance (joint Q above both singles but below their sum); two — Mining ∩ GDP (0.248) and Mining ∩ Population (0.359) — are nonlinear-enhance (joint above the sum), meaning mining's effect on NDVI only emerges alongside the socioeconomic variables.
• Weak singles can matter jointly. GDP alone is the weakest driver (0.099), yet Precipitation ∩ GDP reaches 0.900 — the socioeconomic signal is real but conditional on climate, exactly the kind of finding a single-variable ranking hides.

✓Checkpoint

results/interaction-detector.csv holds all 15 pairs; the strongest is Temperature ∩ Precipitation (Q = 0.916), and 2 of 15 are nonlinear-enhance.

Step 5

Risk detector — which zones are high vs low, and is the difference real

What this step does

riskmean() gives the mean response within each stratum of every variable; gdrisk() t-tests every pair of strata and flags the significant differences. Together they turn a Q value into an actionable map of where the response is high.

Code

terminal Copy

Rscript run-all.R 50

Example output

Demo run

== Step 5/7 Risk detector ===================================== wrote results/risk-means.csv (38 rows) wrote results/risk-detector.csv (104 rows) 104 stratum pairs tested; 82 significantly different (t-test)

◆How to read this result

• Read the direction, not just the magnitude. Climate zone's Q says it matters; the risk means say how — the humid Dwc zone has the highest NDVI gain (0.445), the arid Bwk the lowest (0.005). That ordering is the actionable output.
• The plateau is a signal. NDVI gain rises with precipitation up to ~500 mm then stops climbing — beyond that, water stops being the limiting factor. A monotone Q value cannot show this; the risk means can.
• Significance guards against reading noise. 82 of 104 stratum pairs differ significantly; the 22 that do not (mostly among the socioeconomic strata) are differences you should not report as real.

✓Checkpoint

results/risk-means.csv (38 stratum means) and results/risk-detector.csv (104 t-tests) are written; 82 pairs are significant.

Rscript run-all.R 60

Step 7

Spatial scale (SESU) — is the result robust to the unit size

What this step does

Re-runs the model on the same case aggregated to a series of spatial unit sizes (20–50 km, the built-in ndvi_* series) and reads where the Q values stabilise. The scale where the curves flatten is the recommended analysis unit.

Code

terminal Copy

Rscript run-all.R 70

Example output

Demo run

◆How to read this result

• Stability of the ranking is the robustness claim. The order of variables does not change from 20 to 50 km, so the factor-detector conclusion is not an artefact of one grid resolution — the reassurance a reviewer asks for.
• Slightly rising Q with coarser units is the modifiable areal unit effect. Aggregation smooths within-unit variance, nudging Q up; report the scale you chose and why, not just the single-scale number.

✓Checkpoint

results/sesu-scale.csv holds Q for 6 variables × 4 unit sizes; the mean Q settles near 0.40 from 40 km.

3.4 Reading the Results Together

The four detectors are one story told from four angles: who matters (factor), with whom (interaction), where (risk), and how confidently (ecological). Read alone, each can mislead; read together, they are hard to argue with.

4 Adaptation, Writing & Reproducibility

4.1 Port to Your Domain

The geographical detector moves between domains without touching the analytical scripts: the config declares the table and the variable roles, and every detector is generic. The published applications span the range:

Published reference · road case

What to edit

1. config/project-config.R — point INPUT_FILE at your CSV (or name a built-in DATASET), set RESPONSE, and split the drivers into CATEGORICAL and CONTINUOUS. Choose the DISC_METHODS and DISC_INTERVALS search grid. The single adaptation point of §2.2.
2. For spatial-scale optimisation, either supply your own multi-scale tables or set RUN_SESU <- FALSE; SESU needs the aggregated series, which for a custom INPUT_FILE you build yourself.
3. Nothing else: run Rscript run-all.R test once (a sanity check on the built-in case), then Rscript run-all.R, and read results/factor-detector.csv before drafting.

4.2 Write the Paper

The published OPGD applications share one manuscript architecture: a Methods section that states the discretisation search grid and the four detectors, a Results section that mirrors the pipeline steps, and every number pulled from a results/ file so nothing is retyped. The download ships the code-and-config layer; the tables regenerate from R/90-tables.R.

Manuscript skeleton (shared by the published applications)

1  Introduction        — spatial stratified heterogeneity; the discretisation
                          problem; the OPGD contribution (3 aims)
2  OPGD model          — 2.1 Q value (factor detector; no case data)
                          2.2 Parameter optimisation (discretisation grid; SESU)
                          2.3 Interaction, risk and ecological detectors
3  Case & data         — study area, response, explanatory variables, unit
4  Results             — 4.1 optimal discretisation (Fig: Q vs intervals)
                          4.2 factor detector (ranking + significance)
                          4.3 interaction detector (bubble matrix)
                          4.4 risk detector (zone means)
                          4.5 ecological detector + scale robustness
5  Discussion          — co-dominance, conditional drivers, mechanism, limits
6  Conclusions         — 6-8 sentences mapping to the aims

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 using the CRAN GD package. Reference figures are reproduced from the cited articles for comparison and are © their authors and publishers.

1. Wang J-F, Li X-H, Christakos G, et al. (2010). Geographical detectors-based health risk assessment and its application in the neural tube defects study of the Heshun region, China. IJGIS 24(1):107–127. doi:10.1080/13658810802443457
2. Song Y, Wright G, Wu P, et al. (2018). Segment-based spatial analysis for assessing road infrastructure performance using monitoring observations and remote sensing data. Remote Sensing 10(11):1696. doi:10.3390/rs10111696 (source of the GD road case and the study-area figure on this page)
3. Song Y, Wang J, Ge Y, Xu C (2020). An optimal parameters-based geographical detector model enhances geographic characteristics of explanatory variables for spatial heterogeneity analysis. GIScience & Remote Sensing 57(5):593–610. doi:10.1080/15481603.2020.1760434 · PDF (source of Figs 2, 5, 8 and the NDVI ranking reproduced on this page)
4. The GD package (CRAN): cran.r-project.org/package=GD · vignette (the reference implementation this template wraps)