Reproducing the Local Indicator of Stratified Power (LISP) Model

# =============================================================================
# 00-setup.R — check and install required packages
# Run once: Rscript R/00-setup.R
# =============================================================================

cran_packages <- c(
  "yaml",          # config loading
  "sf",            # spatial data
  "sdsfun",        # spatial weights, discretize_vector
  "geocomplexity", # geocd_vector — GC patterns
  "gdverse",       # gd() — geographical detector q-statistic
  "FNN",           # k-nearest neighbours
  "rpart",         # GOZH decision tree
  "automap",       # autofitVariogram — local range
  "dplyr",
  "parallel",
  # figures
  "ggplot2", "reshape2", "mgcv", "scales", "patchwork", "ggrepel"
)

# localsp provides lisp(), the LISP factor detector.
# If unavailable on CRAN for your R version, install from GitHub:
#   remotes::install_github("stscl/localsp")
other_packages <- c("localsp")

missing <- setdiff(cran_packages, rownames(installed.packages()))
if (length(missing) > 0) {
  cat("Installing:", paste(missing, collapse = ", "), "\n")
  install.packages(missing)
}

for (pkg in other_packages) {
  if (!pkg %in% rownames(installed.packages())) {
    cat(sprintf("Package '%s' missing — trying CRAN, then GitHub.\n", pkg))
    tryCatch(install.packages(pkg), error = function(e) NULL)
    if (!pkg %in% rownames(installed.packages())) {
      if (!"remotes" %in% rownames(installed.packages())) {
        install.packages("remotes")
      }
      remotes::install_github(paste0("stscl/", pkg))
    }
  }
}

cat("\nAll packages present:\n")
for (pkg in c(cran_packages, other_packages)) {
  ok <- requireNamespace(pkg, quietly = TRUE)
  cat(sprintf("  %-14s %s\n", pkg, ifelse(ok, "OK", "MISSING")))
}

# Record the environment for reproducibility.
dir.create("results", showWarnings = FALSE)
writeLines(capture.output(sessionInfo()), "results/session-info.txt")
cat("\nSession info written to results/session-info.txt\n")

# =============================================================================
# LPI project configuration — the single file to edit when porting the
# pipeline to a new domain. Every R script reads this file; none of them
# hard-code variable names, coordinates, or parameters.
#
# Method: LPI (Local Pattern Interaction) = LISP local stratified power
#         + geocomplexity (GC) patterns + local interactions + GOZH total effect
# =============================================================================

project:
  name: "lpi-application-template"
  domain: "example"              # short label for your application domain
  crs_input: 4326                # CRS of the coordinates in the analysis table
  crs_projected: 32650           # metric CRS for distance computation — EDIT
                                 # to a projected CRS suited to your region
                                 # (units must be metres; the demo data sit in
                                 # UTM zone 50N = 32650)

data:
  use_example: true              # true  = 01-prepare-data.R simulates a demo
                                 #         dataset so the pipeline runs end to end
                                 # false = 01-prepare-data.R reads your raw data
                                 #         (edit that script for your source format)
  analysis_table: "data/processed/analysis-table.csv"
  boundary: ""                   # optional study-area boundary for maps, e.g.
                                 # "data/boundary/study-area.gpkg" ("" = no boundary)
  id_col: "uid"
  x_col: "lon"                   # longitude / easting column in analysis table
  y_col: "lat"                   # latitude / northing column in analysis table

response:
  direct_col: "y"                # observed response column in the analysis table
  response_name: "z"             # standard internal name used by the pipeline

predictors:
  x_vars: ["x1", "x2", "x3", "x4", "x5", "x6"]
  # Display labels used in every figure (keep short; GC labels are derived
  # automatically as "GC(<label>)").
  labels:
    x1: "Variable 1"
    x2: "Variable 2"
    x3: "Variable 3"
    x4: "Variable 4"
    x5: "Variable 5"
    x6: "Variable 6"

params:
  gc:                            # Step 03 geocomplexity
    knn_share: 0.02              # k = max(knn_min, knn_share * n)
    knn_min: 30
    method: "spvar"              # geocomplexity::geocd_vector method
  lisp:                          # Steps 04-05 local range + local q
    threshold_multiplier: 2      # local extent = multiplier * variogram range
    discnum: 4                   # discretisation bins inside local windows
    discmethod: "quantile"
  interaction:                   # Step 06 local pairwise interactions
    top_n_vars: 5                # pairs are formed among the top-N variables
                                 # by mean local q (N=5 gives 10 pairs)
    min_local_n: 10              # fallback neighbourhood size if the distance
                                 # circle contains fewer points
  global:                        # Step 08 global benchmark (OPGD-style)
    n_strata: 5
    top_n_interactions: 10
  cores: "auto"                  # "auto" = detectCores() - 1, or an integer

validation:
  holdout:                       # Step 09 subsample stability
    n_repeats: 20
    drop_share: 0.05
    top_n_vars: 10
    seed: 2026
  sensitivity:                   # Step 10 parameter sensitivity
    threshold_multipliers: [1.5, 2, 3]
    discnums: [3, 4, 5]
    top_n_vars: 5

figures:
  # Fixed six-class palette for quantile-binned maps (colour-method mapping
  # stays constant across all figures).
  palette: ["#1F6B4A", "#4F9C93", "#8FC8BE", "#E6A8B6", "#D85A8E", "#8F1E5F"]
  interaction_colors:
    xx: "#1F6B4A"                # X x X pairs
    xgc: "#8FC8BE"               # X x GC pairs
    gcgc: "#8F1E5F"              # GC x GC pairs
  point_size: 0.6
  width_single: 5                # inches, single-panel figures
  width_full: 14                 # inches, four-panel map rows

# =============================================================================
# fn-config.R — configuration loader and small shared utilities
# All pipeline scripts source this file first and run from the project root.
# =============================================================================

load_config <- function(path = "config/project-config.yaml") {
  if (!file.exists(path)) {
    stop("Config not found at '", path,
         "'. Run scripts from the project root, e.g. Rscript R/03-geocomplexity.R")
  }
  yaml::read_yaml(path)
}

get_cores <- function(cfg) {
  if (identical(cfg$params$cores, "auto")) {
    max(1L, parallel::detectCores() - 1L)
  } else {
    as.integer(cfg$params$cores)
  }
}

ensure_dir <- function(path) {
  dir.create(path, showWarnings = FALSE, recursive = TRUE)
  invisible(path)
}

read_analysis_table <- function(cfg) {
  path <- cfg$data$analysis_table
  if (!file.exists(path)) {
    stop("Analysis table not found at '", path,
         "'. Run R/01-prepare-data.R first.")
  }
  read.csv(path)
}

# Column-name helpers -----------------------------------------------------

gc_names_of <- function(cfg) paste0("gc_", cfg$predictors$x_vars)

all_vars_of <- function(cfg) c(cfg$predictors$x_vars, gc_names_of(cfg))

# Named label vectors for figures: x labels from config, GC labels derived.
get_labels <- function(cfg) {
  x_vars   <- cfg$predictors$x_vars
  x_labels <- unlist(cfg$predictors$labels)[x_vars]
  gc_labels <- setNames(paste0("GC(", x_labels, ")"), gc_names_of(cfg))
  c(x_labels, gc_labels)
}

# Read the response (with coordinates) produced by 02-response-intake.R
read_response <- function(cfg) {
  path <- "results/step02-response.csv"
  if (!file.exists(path)) {
    stop("Response file not found at '", path,
         "'. Run R/02-response-intake.R first.")
  }
  read.csv(path)
}

# Read the local extent estimated by 04-local-range.R
read_local_range <- function(cfg) {
  path <- "results/step04-local-range.csv"
  if (!file.exists(path)) {
    stop("Local range file not found at '", path,
         "'. Run R/04-local-range.R first.")
  }
  read.csv(path)
}

# Assemble the modelling frame: analysis table + response + geocomplexity.
# Every script from Step 05 onward starts from this call, so any script can
# be re-run in isolation without recomputing upstream steps.
build_model_frame <- function(cfg) {
  d  <- read_analysis_table(cfg)
  rz <- read_response(cfg)
  d[[cfg$response$response_name]] <- rz[[cfg$response$response_name]]

  gc_path <- "results/step03-geocomplexity.csv"
  if (file.exists(gc_path)) {
    gc <- read.csv(gc_path)
    for (v in gc_names_of(cfg)) d[[v]] <- gc[[v]]
  }
  d
}

# =============================================================================
# 01-prepare-data.R — build the analysis table
#
# Output: data/processed/analysis-table.csv with EXACTLY this schema:
#   <id_col>   unique unit id                     (config: data$id_col)
#   <x_col>    longitude or easting               (config: data$x_col)
#   <y_col>    latitude or northing               (config: data$y_col)
#   y          the observed response variable     (config: response$direct_col)
#   x1..xp     the explanatory variables listed in predictors$x_vars
#
# Two modes (config: data$use_example):
#   true  — simulate a synthetic demo dataset so the whole pipeline runs
#           end to end with no external data (recommended first run)
#   false — adapt the marked section below to your own raw data
# =============================================================================

source("R/functions/fn-config.R")
cfg <- load_config()
ensure_dir("data/processed")

if (isTRUE(cfg$data$use_example)) {

  # ---------------------------------------------------------------------------
  # EXAMPLE MODE — synthetic spatial dataset (n = 600 point units)
  # Design: two latent regimes (west/east) with different dominant drivers,
  # so local q values genuinely vary across space; x5 x x6 carry a joint
  # (interaction) signal; y is the directly observed response.
  # ---------------------------------------------------------------------------
  set.seed(42)
  n <- 600

  lon <- runif(n, 116.0, 117.0)   # ~100 km east-west extent
  lat <- runif(n, 36.0, 36.9)

  east <- (lon - min(lon)) / diff(range(lon))   # 0 west -> 1 east
  north <- (lat - min(lat)) / diff(range(lat))

  smooth_field <- function(fx, fy, phase = 0, noise = 0.15) {
    field <- sin(2 * pi * fx * east + phase) * cos(2 * pi * fy * north) +
      0.5 * north
    as.numeric(scale(field + rnorm(n, 0, noise)))
  }

  x1 <- smooth_field(1.0, 0.5)
  x2 <- smooth_field(0.5, 1.0, phase = 1)
  x3 <- smooth_field(1.5, 0.7, phase = 2)
  x4 <- as.numeric(scale(east + rnorm(n, 0, 0.3)))
  x5 <- smooth_field(0.8, 1.2, phase = 3)
  x6 <- smooth_field(1.2, 0.4, phase = 4)

  # Observed response: driver x1 dominates in the west, x2 in the east,
  # x5*x6 adds a nonlinear-enhance interaction, x3 a weak global effect.
  latent <- (1 - east) * 1.2 * x1 + east * 1.2 * x2 +
    0.8 * x5 * x6 + 0.3 * x3 + rnorm(n, 0, 0.4)
  y_obs <- (latent - min(latent)) / diff(range(latent))

  d <- data.frame(
    uid = seq_len(n),
    lon = lon,
    lat = lat,
    y   = y_obs,
    x1 = x1, x2 = x2, x3 = x3, x4 = x4, x5 = x5, x6 = x6
  )

  write.csv(d, cfg$data$analysis_table, row.names = FALSE)
  cat(sprintf("Example analysis table written: %s (%d rows)\n",
              cfg$data$analysis_table, nrow(d)))
  cat("Set data$use_example to false in config/project-config.yaml\n")
  cat("and adapt this script once your own data are ready.\n")

} else {

  # ---------------------------------------------------------------------------
  # REAL-DATA MODE — adapt this section to your domain.
  # Typical steps:
  #   1. read raw files from data/raw/
  #   2. aggregate to your spatial unit (grid cell, county, segment, station)
  #   3. compute or join coordinates for each unit
  #   4. rename columns to match config (id, coords, response, x_vars)
  #   5. drop units with missing coordinates; document any imputation
  # ---------------------------------------------------------------------------
  raw_path <- "data/raw/raw-data.csv"   # EDIT: your raw source file
  if (!file.exists(raw_path)) {
    stop("Real-data mode: expected raw data at '", raw_path,
         "'. Edit R/01-prepare-data.R for your source format.")
  }
  raw <- read.csv(raw_path)

  # EDIT from here: aggregation and renaming for your domain. The example
  # below simply passes columns through and keeps complete rows.
  needed <- c(cfg$data$id_col, cfg$data$x_col, cfg$data$y_col,
              cfg$response$direct_col, cfg$predictors$x_vars)
  missing_cols <- setdiff(needed, names(raw))
  if (length(missing_cols) > 0) {
    stop("Raw data lacks required columns: ",
         paste(missing_cols, collapse = ", "),
         ". Adapt R/01-prepare-data.R or config/project-config.yaml.")
  }

  d <- raw[complete.cases(raw[, needed]), needed]
  write.csv(d, cfg$data$analysis_table, row.names = FALSE)
  cat(sprintf("Analysis table written: %s (%d rows, %d dropped as incomplete)\n",
              cfg$data$analysis_table, nrow(d), nrow(raw) - nrow(d)))
}

# =============================================================================
# fn-lpi-core.R — the generalizable LPI method core
# Source: cc005 LPI pipeline (main_analysis.R), generalized and config-driven.
#
# Components (in pipeline order):
#   build_spatial()          projected coordinates + full distance matrix
#   estimate_local_range()   semivariogram range -> local extent (threshold)
#   run_lisp_factor()        LISP local q per variable (localsp::lisp)
#   run_local_interaction()  local pairwise interaction q + 5-class typing
#   run_interactions_all()   parallel loop of run_local_interaction()
#   run_gozh()               GOZH total effect: rpart zones -> gdverse::gd
#   run_local_gozh()         local GOZH over k-nearest neighbourhoods
#   run_global_gd()          OPGD-style global q for one variable
#   run_global_interaction_q()  global q for one variable pair
# =============================================================================

# --- Spatial scaffolding -----------------------------------------------------

# Projected coordinate matrix (metres) and full pairwise distance matrix.
build_spatial <- function(d, cfg) {
  sf_pts  <- sf::st_as_sf(d, coords = c(cfg$data$x_col, cfg$data$y_col),
                          crs = cfg$project$crs_input)
  sf_proj <- sf::st_transform(sf_pts, crs = cfg$project$crs_projected)
  coords_mat <- sf::st_coordinates(sf_proj)
  list(sf_proj    = sf_proj,
       coords_mat = coords_mat,
       distmat    = as.matrix(dist(coords_mat)))
}

# Semivariogram of the response on projected coordinates; the local extent
# (distance threshold) is multiplier * range, following LISP (Hu et al. 2024).
estimate_local_range <- function(d, cfg, spatial = NULL) {
  if (is.null(spatial)) spatial <- build_spatial(d, cfg)
  response <- cfg$response$response_name
  sp_obj <- as(spatial$sf_proj, "Spatial")
  v <- automap::autofitVariogram(
    as.formula(paste(response, "~ 1")), sp_obj
  )
  vg_range  <- v$var_model$range[2]
  threshold <- vg_range * cfg$params$lisp$threshold_multiplier
  list(variogram = v, range = vg_range, threshold = threshold)
}

# --- LISP factor detector ------------------------------------------------

# Local q (power of determinant) of every variable at every location.
# Returns list(q = data.frame, sig = data.frame) with NA imputed:
# q -> 0 (no explanatory power), p -> 1 (not significant).
run_lisp_factor <- function(d, cfg, vars, threshold, distmat,
                            discnum = NULL, cores = NULL) {
  response <- cfg$response$response_name
  if (is.null(discnum)) discnum <- cfg$params$lisp$discnum
  if (is.null(cores))   cores   <- get_cores(cfg)

  d_lisp <- d[, c(response, vars)]
  lisp_result <- localsp::lisp(
    formula    = as.formula(paste(response, "~ .")),
    data       = d_lisp,
    threshold  = threshold,
    distmat    = distmat,
    discnum    = discnum,
    discmethod = cfg$params$lisp$discmethod,
    cores      = cores
  )
  lisp_df <- as.data.frame(lisp_result)

  pd_cols  <- grep("^pd_",  names(lisp_df), value = TRUE)
  sig_cols <- grep("^sig_", names(lisp_df), value = TRUE)

  base_cols <- data.frame(
    LocationID = seq_len(nrow(d)),
    Longitude  = d[[cfg$data$x_col]],
    Latitude   = d[[cfg$data$y_col]]
  )
  local_q   <- cbind(base_cols, lisp_df[, pd_cols,  drop = FALSE])
  local_sig <- cbind(base_cols, lisp_df[, sig_cols, drop = FALSE])
  names(local_q)[-(1:3)]   <- sub("^pd_",  "", pd_cols)
  names(local_sig)[-(1:3)] <- sub("^sig_", "", sig_cols)

  var_cols <- setdiff(names(local_q), c("LocationID", "Longitude", "Latitude"))
  local_q[var_cols]   <- lapply(local_q[var_cols],
                                function(x) { x[is.na(x)] <- 0;   x })
  local_sig[var_cols] <- lapply(local_sig[var_cols],
                                function(x) { x[is.na(x)] <- 1.0; x })

  list(q = local_q, sig = local_sig)
}

# --- Local pairwise interaction detector ----------------------------------

classify_interaction <- function(qv1, qv2, qv12) {
  dplyr::case_when(
    qv12 < min(qv1, qv2)                          ~ "Weaken, nonlinear",
    qv12 >= min(qv1, qv2) & qv12 < max(qv1, qv2)  ~ "Weaken, uni-",
    qv12 == qv1 + qv2                             ~ "Independent",
    qv12 > max(qv1, qv2) & qv12 < qv1 + qv2       ~ "Enhance, bi-",
    qv12 > qv1 + qv2                              ~ "Enhance, nonlinear",
    TRUE                                          ~ NA_character_
  )
}

# One location: local window by distance threshold (fallback: min_local_n
# nearest points), quantile discretization of each pair member, single and
# combined-zone q via gdverse::gd, then 5-class interaction typing.
run_local_interaction <- function(i, data, distmat, threshold,
                                  response, var_pairs, discnum = 4,
                                  min_local_n = 10) {
  dists     <- distmat[i, ]
  local_idx <- which(dists <= threshold)
  if (length(local_idx) < min_local_n) {
    local_idx <- order(dists)[1:min(min_local_n, nrow(data))]
  }
  local_data <- data[local_idx, , drop = FALSE]

  empty_row <- function(v1, v2) {
    data.frame(var1 = v1, var2 = v2,
               qv1 = NA, qv2 = NA, qv12 = NA, pv12 = NA,
               interaction = NA, LocationID = i, stringsAsFactors = FALSE)
  }

  pair_results <- lapply(var_pairs, function(pair) {
    v1 <- pair[1]; v2 <- pair[2]
    if (length(unique(local_data[[v1]])) < 2 ||
        length(unique(local_data[[v2]])) < 2) {
      return(empty_row(v1, v2))
    }

    tryCatch({
      disc_v1 <- tryCatch(
        sdsfun::discretize_vector(local_data[[v1]], n = discnum,
                                  method = "quantile"),
        error = function(e) NULL)
      disc_v2 <- tryCatch(
        sdsfun::discretize_vector(local_data[[v2]], n = discnum,
                                  method = "quantile"),
        error = function(e) NULL)
      if (is.null(disc_v1) || is.null(disc_v2)) stop("discretization failed")

      tmp <- local_data
      tmp[["disc_v1"]]  <- disc_v1
      tmp[["disc_v2"]]  <- disc_v2
      tmp[["disc_v12"]] <- paste(disc_v1, disc_v2, sep = "_")

      y <- tmp[[response]]
      if (sum((y - mean(y))^2) == 0) stop("zero total SS")

      gd_v1  <- gdverse::gd(as.formula(paste(response, "~ disc_v1")),  data = tmp)
      gd_v2  <- gdverse::gd(as.formula(paste(response, "~ disc_v2")),  data = tmp)
      gd_v12 <- gdverse::gd(as.formula(paste(response, "~ disc_v12")), data = tmp)

      qv1  <- as.numeric(gd_v1$factor[["Q-statistic"]])
      qv2  <- as.numeric(gd_v2$factor[["Q-statistic"]])
      qv12 <- as.numeric(gd_v12$factor[["Q-statistic"]])
      pv12 <- as.numeric(gd_v12$factor[["P-value"]])

      data.frame(var1 = v1, var2 = v2,
                 qv1 = qv1, qv2 = qv2, qv12 = qv12, pv12 = pv12,
                 interaction = classify_interaction(qv1, qv2, qv12),
                 LocationID = i, stringsAsFactors = FALSE)
    }, error = function(e) empty_row(v1, v2))
  })

  do.call(rbind, pair_results)
}

# Parallel driver over all locations. Returns one long data.frame.
run_interactions_all <- function(d, cfg, var_pairs, distmat, threshold) {
  response    <- cfg$response$response_name
  discnum     <- cfg$params$lisp$discnum
  min_local_n <- cfg$params$interaction$min_local_n
  n_cores     <- get_cores(cfg)

  worker <- function(i) {
    run_local_interaction(
      i = i, data = d, distmat = distmat, threshold = threshold,
      response = response, var_pairs = var_pairs,
      discnum = discnum, min_local_n = min_local_n
    )
  }

  if (n_cores > 1) {
    cl <- parallel::makeCluster(n_cores)
    on.exit(parallel::stopCluster(cl), add = TRUE)
    parallel::clusterExport(
      cl,
      varlist = c("d", "distmat", "threshold", "response", "var_pairs",
                  "discnum", "min_local_n",
                  "run_local_interaction", "classify_interaction"),
      envir = environment()
    )
    parallel::clusterEvalQ(cl, {
      library(gdverse); library(sdsfun); library(dplyr)
    })
    inter_list <- parallel::parLapply(cl, seq_len(nrow(d)), worker)
  } else {
    inter_list <- lapply(seq_len(nrow(d)), function(i) {
      if (i %% 200 == 0) cat(sprintf("  interaction: location %d / %d\n",
                                     i, nrow(d)))
      worker(i)
    })
  }

  do.call(rbind, Filter(Negate(is.null), inter_list))
}

# --- GOZH total effect -----------------------------------------------------

# rpart regression tree on all predictors -> leaf zones -> combined-zone q.
run_gozh <- function(data, response, predictors) {
  formula_tree <- as.formula(
    paste(response, "~", paste(predictors, collapse = " + "))
  )
  tree <- tryCatch(rpart::rpart(formula_tree, data = data, method = "anova"),
                   error = function(e) NULL)
  if (is.null(tree)) return(c(q = NA_real_, p = NA_real_))

  data$.zone <- as.character(as.numeric(tree$where))
  result <- tryCatch(
    gdverse::gd(as.formula(paste(response, "~ .zone")), data = data),
    error = function(e) NULL)
  if (is.null(result)) return(c(q = NA_real_, p = NA_real_))

  qv <- tryCatch(as.numeric(result$factor[["Q-statistic"]]),
                 error = function(e) NA_real_)
  pv <- tryCatch(as.numeric(result$factor[["P-value"]]),
                 error = function(e) NA_real_)
  if (length(qv) == 0) qv <- NA_real_
  if (length(pv) == 0) pv <- NA_real_
  c(q = qv, p = pv)
}

# Local GOZH: run_gozh() inside each location's k-nearest neighbourhood.
run_local_gozh <- function(d, cfg, predictors, coords_mat, k_local) {
  response <- cfg$response$response_name
  n <- nrow(d)
  local_q <- rep(NA_real_, n)
  local_p <- rep(NA_real_, n)
  for (i in seq_len(n)) {
    if (i %% 200 == 0) cat(sprintf("  local GOZH: location %d / %d\n", i, n))
    nn_idx <- FNN::get.knnx(coords_mat, coords_mat[i, , drop = FALSE],
                            k = k_local)$nn.index
    res <- run_gozh(d[nn_idx, ], response, predictors)
    local_q[i] <- res["q"]
    local_p[i] <- res["p"]
  }
  data.frame(local_gozh = local_q, local_gozh_p = local_p)
}

# --- Global benchmarks (OPGD-style) ----------------------------------------

run_global_gd <- function(data, response, predictor, n_strata = 5) {
  tryCatch({
    df <- data.frame(y = data[[response]], x = data[[predictor]])
    df <- df[complete.cases(df), ]
    if (nrow(df) < 30) return(c(q = NA_real_, p = NA_real_))
    df$x_d <- cut(df$x, breaks = n_strata, labels = FALSE,
                  include.lowest = TRUE)
    result <- gdverse::gd(y ~ x_d, data = df)
    qv <- tryCatch(as.numeric(result$factor[["Q-statistic"]]),
                   error = function(e) NA_real_)
    pv <- tryCatch(as.numeric(result$factor[["P-value"]]),
                   error = function(e) NA_real_)
    if (length(qv) == 0) qv <- NA_real_
    if (length(pv) == 0) pv <- NA_real_
    c(q = qv, p = pv)
  }, error = function(e) c(q = NA_real_, p = NA_real_))
}

run_global_interaction_q <- function(data, response, v1, v2, n_strata = 5) {
  tryCatch({
    df <- data.frame(y = data[[response]], x1 = data[[v1]], x2 = data[[v2]])
    df <- df[complete.cases(df), ]
    if (nrow(df) < 30) return(c(q = NA_real_, p = NA_real_))
    df$x1_d <- cut(df$x1, breaks = n_strata, labels = FALSE,
                   include.lowest = TRUE)
    df$x2_d <- cut(df$x2, breaks = n_strata, labels = FALSE,
                   include.lowest = TRUE)
    df$x12 <- interaction(df$x1_d, df$x2_d, drop = TRUE)
    result <- gdverse::gd(y ~ x12, data = df)
    qv <- tryCatch(as.numeric(result$factor[["Q-statistic"]]),
                   error = function(e) NA_real_)
    pv <- tryCatch(as.numeric(result$factor[["P-value"]]),
                   error = function(e) NA_real_)
    if (length(qv) == 0) qv <- NA_real_
    if (length(pv) == 0) pv <- NA_real_
    c(q = qv, p = pv)
  }, error = function(e) c(q = NA_real_, p = NA_real_))
}

# =============================================================================
# 11-figures.R — the eight core figures of an LPI application paper
# Reads results/step*.csv only; no recomputation. Figure plan and the
# argumentation task of each figure: figures/figure-plan.md
#
#   fig01 response index maps            fig05 interaction strength bars
#   fig02 X vs GC(X) trend curves        fig06 interaction type map + counts
#   fig03 local q histograms (X, GC)     fig07 local GOZH map + histogram
#   fig04 local q maps (top X + top GC)  fig08 variogram + local vs global q
# =============================================================================

source("R/functions/fn-config.R")
source("R/functions/fn-plot-helpers.R")

library(ggplot2)
library(reshape2)
library(sf)
library(dplyr)
library(mgcv)
library(patchwork)

cfg <- load_config()
ensure_dir("figures")

x_vars   <- cfg$predictors$x_vars
gc_names <- gc_names_of(cfg)
all_vars <- all_vars_of(cfg)
labels   <- get_labels(cfg)
resp     <- cfg$response$response_name
w_full   <- cfg$figures$width_full
w_single <- cfg$figures$width_single

d_z    <- read.csv("results/step02-response.csv")
d_gc   <- read.csv("results/step03-geocomplexity.csv")
d_q    <- read.csv("results/step05-local-q-values.csv")
d_iall <- read.csv("results/step06-all-pairwise-interactions.csv")
d_ibst <- read.csv("results/step06-local-best-interaction.csv")
d_gozh <- read.csv("results/step07-local-gozh.csv")
d_vt   <- read.csv("results/step08-validation-table.csv")
d_vg   <- read.csv("results/step04-variogram-points.csv")
d_lr   <- read.csv("results/step04-local-range.csv")

xc <- cfg$data$x_col
yc <- cfg$data$y_col

# =============================================================================
# FIGURE 1: Spatial distribution of the response index (and sub-indices)
# =============================================================================

resp_cols <- intersect(
  c(grep(paste0("^", resp, "_"), names(d_z), value = TRUE), resp),
  names(d_z)
)
fig1_panels <- lapply(resp_cols, function(v) {
  make_value_map(d_z[[v]], d_z[[xc]], d_z[[yc]], cfg,
                 title = v, legend_title = "Index")
})
p_fig1 <- wrap_plots(fig1_panels, ncol = min(4, length(fig1_panels)))
ggsave("figures/fig01-response-maps.pdf", p_fig1,
       width = min(w_full, 3.5 * length(fig1_panels) + 1), height = 4.5,
       dpi = 300)
cat("Figure 1 saved.\n")

# =============================================================================
# FIGURE 2: X value vs GC(X), 5%-quantile binned means + GAM smoother
# =============================================================================

# Raw variable values live in the analysis table; GC values in step03.
d_raw  <- read.csv(cfg$data$analysis_table)
d_main <- merge(d_raw, d_gc[, c(cfg$data$id_col, gc_names)],
                by = cfg$data$id_col)
gc_summary <- do.call(rbind, lapply(x_vars, function(v) {
  df <- data.frame(x = d_main[[v]], y = d_main[[paste0("gc_", v)]])
  q_breaks <- quantile(df$x, probs = seq(0, 1, 0.05), na.rm = TRUE)
  df$group <- cut(df$x, breaks = unique(q_breaks),
                  labels = FALSE, include.lowest = TRUE)
  smry <- df %>%
    group_by(group) %>%
    summarise(mean_x = mean(x, na.rm = TRUE),
              mean_y = mean(y, na.rm = TRUE),
              se_y   = sd(y, na.rm = TRUE) / sqrt(sum(!is.na(y))),
              .groups = "drop")
  smry$variable <- labels[[v]]
  smry
}))
gc_summary$variable <- factor(gc_summary$variable, levels = labels[x_vars])

p_fig2 <- ggplot(gc_summary, aes(x = mean_x, y = mean_y)) +
  geom_errorbar(aes(ymin = mean_y - se_y, ymax = mean_y + se_y),
                width = 0, linewidth = 0.4, colour = "grey50") +
  geom_point(size = 1.2, colour = "steelblue") +
  geom_smooth(colour = "red", linewidth = 0.8, se = TRUE,
              fill = "pink", alpha = 0.3) +
  facet_wrap(~ variable, scales = "free", ncol = 3) +
  labs(x = "Variable value", y = "Geocomplexity (GC)") +
  facet_theme()
ggsave("figures/fig02-gc-trends.pdf", p_fig2, width = 8, height = 6, dpi = 300)
cat("Figure 2 saved.\n")

# =============================================================================
# FIGURE 3: Local q distributions — X variables (3a) and GC patterns (3b)
# =============================================================================

plot_q_histograms(d_q, x_vars,   labels, 3, "figures/fig03a-q-hist-x.pdf")
plot_q_histograms(d_q, gc_names, labels, 3, "figures/fig03b-q-hist-gc.pdf")
cat("Figure 3 saved.\n")

# =============================================================================
# FIGURE 4: Spatial maps of local q — top X and top GC variables
# =============================================================================

n_top    <- min(4, length(x_vars))
mean_q_x  <- sort(colMeans(d_q[, x_vars],   na.rm = TRUE), decreasing = TRUE)
mean_q_gc <- sort(colMeans(d_q[, gc_names], na.rm = TRUE), decreasing = TRUE)
top_map_vars <- c(names(mean_q_x)[1:n_top], names(mean_q_gc)[1:n_top])

fig4_panels <- lapply(top_map_vars, function(v) {
  make_value_map(d_q[[v]], d_q$Longitude, d_q$Latitude, cfg,
                 title = labels[[v]], legend_title = "Local q")
})
p_fig4 <- wrap_plots(fig4_panels, ncol = n_top)
ggsave("figures/fig04-local-q-maps.pdf", p_fig4,
       width = w_full, height = 7, dpi = 300)
cat("Figure 4 saved.\n")

# =============================================================================
# FIGURE 5: Interaction strength — top pairs by mean qv12, by pair category
# =============================================================================

pair_label_of <- function(v) {
  ifelse(v %in% names(labels), labels[v], v)
}

top_pairs <- d_iall %>%
  group_by(var1, var2) %>%
  summarise(mean_qv12 = mean(qv12, na.rm = TRUE),
            sd_qv12   = sd(qv12,   na.rm = TRUE),
            .groups = "drop") %>%
  mutate(vartype = vartype_of(var1, var2, x_vars, gc_names)) %>%
  arrange(desc(mean_qv12)) %>%
  slice_head(n = 15) %>%
  mutate(pair_label = paste(pair_label_of(var1), "×", pair_label_of(var2)))
top_pairs$pair_label <- factor(top_pairs$pair_label,
                               levels = rev(top_pairs$pair_label))

p_fig5 <- ggplot(top_pairs,
                 aes(x = pair_label, y = mean_qv12, fill = vartype)) +
  geom_bar(stat = "identity", alpha = 0.7) +
  geom_errorbar(aes(ymin = mean_qv12 - sd_qv12, ymax = mean_qv12 + sd_qv12),
                width = 0.4, colour = "orange", alpha = 0.9, linewidth = 1.0) +
  scale_fill_manual(values = interaction_palette(cfg),
                    name = "Interaction type") +
  coord_flip() +
  labs(x = NULL, y = "Mean interaction q-value (qv12) ± SD") +
  theme_bw() +
  theme(axis.text = element_text(size = 8),
        legend.position = "bottom",
        panel.grid.minor = element_blank())
ggsave("figures/fig05-interaction-strength.pdf", p_fig5,
       width = 8, height = 6, dpi = 300)
cat("Figure 5 saved.\n")

# =============================================================================
# FIGURE 6: Dominant interaction category — spatial map + frequency bars
# =============================================================================

d_ibst$vartype <- vartype_of(d_ibst$var1, d_ibst$var2, x_vars, gc_names)

p_fig6a <- ggplot() +
  boundary_layer(cfg) +
  geom_point(data = d_ibst,
             aes(x = Longitude, y = Latitude, colour = vartype),
             size = cfg$figures$point_size, alpha = 0.8) +
  scale_colour_manual(values = interaction_palette(cfg),
                      name = "Dominant interaction") +
  labs(x = "Longitude", y = "Latitude") +
  map_theme() +
  guides(colour = guide_legend(override.aes = list(size = 3)))

type_count <- as.data.frame(table(d_ibst$vartype))
colnames(type_count) <- c("Type", "Count")
type_count$Pct <- round(type_count$Count / sum(type_count$Count) * 100, 1)

p_fig6b <- ggplot(type_count, aes(x = Type, y = Count, fill = Type)) +
  geom_bar(stat = "identity", colour = "black", linewidth = 0.3) +
  geom_text(aes(label = paste0(Pct, "%")), vjust = -0.4, size = 4) +
  scale_fill_manual(values = interaction_palette(cfg)) +
  labs(x = "Interaction type", y = "Number of locations") +
  theme_bw() +
  theme(legend.position = "none", panel.grid.minor = element_blank())

p_fig6 <- p_fig6a + p_fig6b + plot_layout(widths = c(1.4, 1))
ggsave("figures/fig06-interaction-type.pdf", p_fig6,
       width = 10, height = 4.5, dpi = 300)
cat("Figure 6 saved.\n")

# =============================================================================
# FIGURE 7: Local GOZH total effect — map + distribution
# =============================================================================

p_fig7a <- make_value_map(d_gozh$local_gozh, d_gozh[[xc]], d_gozh[[yc]], cfg,
                          title = NULL, legend_title = "GOZH q (all vars)")

mean_gozh <- mean(d_gozh$local_gozh, na.rm = TRUE)
p_fig7b <- ggplot(d_gozh, aes(x = local_gozh)) +
  geom_histogram(aes(y = after_stat(density)), bins = 40,
                 fill = "lightblue", colour = "black", alpha = 0.7) +
  geom_density(colour = "red", linewidth = 0.8) +
  geom_vline(xintercept = mean_gozh, colour = "blue",
             linetype = "dashed", linewidth = 0.8) +
  annotate("text", x = Inf, y = Inf,
           label = paste("Mean:", round(mean_gozh, 3)),
           hjust = 1.15, vjust = 1.4, size = 4) +
  labs(x = "Local GOZH q-value (all variables)", y = "Density") +
  theme_bw() +
  theme(panel.grid.minor = element_blank())

p_fig7 <- p_fig7a + p_fig7b + plot_layout(widths = c(1.2, 1))
ggsave("figures/fig07-gozh-total-effect.pdf", p_fig7,
       width = 10, height = 4.5, dpi = 300)
cat("Figure 7 saved.\n")

# =============================================================================
# FIGURE 8: Method diagnostics — variogram (8a) + local vs global q (8b)
# =============================================================================

p_fig8a <- ggplot(d_vg, aes(x = dist_m / 1000, y = gamma)) +
  geom_point(size = 1.5, colour = "steelblue") +
  geom_vline(xintercept = d_lr$range_m[1] / 1000,
             colour = "red", linetype = "dashed") +
  geom_vline(xintercept = d_lr$threshold_m[1] / 1000,
             colour = "blue", linetype = "dashed") +
  annotate("text", x = d_lr$range_m[1] / 1000, y = Inf,
           label = "range", hjust = -0.1, vjust = 1.5,
           size = 3, colour = "red") +
  annotate("text", x = d_lr$threshold_m[1] / 1000, y = Inf,
           label = "local extent", hjust = 1.1, vjust = 1.5,
           size = 3, colour = "blue") +
  labs(x = "Distance (km)", y = "Semivariance") +
  theme_bw() +
  theme(panel.grid.minor = element_blank())

d_cmp <- d_vt %>%
  filter(Variable_type %in% c("X", "GC")) %>%
  filter(!is.na(Global_q))
d_cmp$label <- ifelse(d_cmp$Variable %in% names(labels),
                      labels[d_cmp$Variable], d_cmp$Variable)

p_fig8b <- ggplot(d_cmp, aes(x = Global_q, y = LPI_mean,
                             colour = Variable_type)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed",
              colour = "grey60") +
  geom_errorbar(aes(ymin = LPI_min, ymax = LPI_max),
                width = 0, linewidth = 0.4, alpha = 0.6) +
  geom_point(size = 2) +
  ggrepel::geom_text_repel(aes(label = label), size = 2.5,
                           show.legend = FALSE) +
  scale_colour_manual(values = c("X" = unlist(cfg$figures$palette)[1],
                                 "GC" = unlist(cfg$figures$palette)[6]),
                      name = "Variable type") +
  labs(x = "Global q (OPGD benchmark)",
       y = "Local q, mean with [min, max]") +
  theme_bw() +
  theme(panel.grid.minor = element_blank())

p_fig8 <- p_fig8a + p_fig8b + plot_layout(widths = c(1, 1.3))
ggsave("figures/fig08-diagnostics-benchmark.pdf", p_fig8,
       width = 10, height = 4.5, dpi = 300)
cat("Figure 8 saved.\n")

cat("\n=== All figures saved to figures/ ===\n")

# =============================================================================
# fn-plot-helpers.R — shared plotting components
# Source: cc005 LPI pipeline (plot results.R), generalized.
# Colour-method mapping is bound once from config and reused in every figure
# (writing style guide Rule 26: fixed colour mapping across all figures).
# =============================================================================

map_theme <- function() {
  ggplot2::theme_bw() +
    ggplot2::theme(
      axis.text    = ggplot2::element_text(size = 8),
      legend.title = ggplot2::element_text(size = 9),
      legend.text  = ggplot2::element_text(size = 8),
      panel.grid   = ggplot2::element_blank()
    )
}

facet_theme <- function() {
  ggplot2::theme_bw() +
    ggplot2::theme(
      strip.text       = ggplot2::element_text(size = 7.5),
      axis.text        = ggplot2::element_text(size = 7),
      axis.title       = ggplot2::element_text(size = 9),
      panel.grid.minor = ggplot2::element_blank()
    )
}

# Optional study-area boundary layer ("" in config -> empty layer list).
boundary_layer <- function(cfg) {
  if (is.null(cfg$data$boundary) || cfg$data$boundary == "") return(list())
  boundary <- sf::st_read(cfg$data$boundary, quiet = TRUE)
  list(ggplot2::geom_sf(data = boundary, fill = "grey95",
                        colour = "grey60", linewidth = 0.3))
}

# One map panel with its own 6-quantile legend (equal-frequency bins).
make_value_map <- function(values, x_coords, y_coords, cfg,
                           title = NULL, legend_title = "Value") {
  palette  <- unlist(cfg$figures$palette)
  q_breaks <- unique(quantile(values, probs = seq(0, 1, length.out = 7),
                              na.rm = TRUE))
  n_bins <- length(q_breaks) - 1
  bin_labels <- paste0(round(q_breaks[-(n_bins + 1)], 3), "-",
                       round(q_breaks[-1], 3))

  df_map <- data.frame(
    x_plot  = x_coords,
    y_plot  = y_coords,
    value_q = cut(values, breaks = q_breaks, labels = bin_labels,
                  include.lowest = TRUE)
  )

  ggplot2::ggplot() +
    boundary_layer(cfg) +
    ggplot2::geom_point(
      data = df_map,
      ggplot2::aes(x = x_plot, y = y_plot, colour = value_q),
      size = cfg$figures$point_size, alpha = 0.8
    ) +
    ggplot2::scale_colour_manual(
      values = setNames(palette[1:n_bins], bin_labels),
      name   = legend_title,
      guide  = ggplot2::guide_legend(override.aes = list(size = 2.5))
    ) +
    ggplot2::labs(title = title, x = "Longitude", y = "Latitude") +
    map_theme() +
    ggplot2::theme(
      plot.title      = ggplot2::element_text(size = 9, hjust = 0.5),
      legend.key.size = ggplot2::unit(0.4, "cm"),
      legend.text     = ggplot2::element_text(size = 7)
    )
}

# Faceted histogram + density + mean line of local q values.
plot_q_histograms <- function(q_df, var_list, label_map, ncol_val, fname,
                              width = 8, height = 6) {
  q_long <- reshape2::melt(q_df[, c("LocationID", var_list)],
                           id.vars = "LocationID")
  q_long$variable <- factor(label_map[as.character(q_long$variable)],
                            levels = label_map[var_list])
  mean_vals <- aggregate(value ~ variable, data = q_long,
                         FUN = mean, na.rm = TRUE)
  mean_vals$value <- round(mean_vals$value, 3)

  p <- ggplot2::ggplot(q_long, ggplot2::aes(x = value)) +
    ggplot2::geom_histogram(
      ggplot2::aes(y = ggplot2::after_stat(density)), bins = 40,
      fill = "lightblue", colour = "black", alpha = 0.7) +
    ggplot2::geom_density(colour = "red", linewidth = 0.8) +
    ggplot2::geom_vline(
      data = mean_vals, ggplot2::aes(xintercept = value),
      colour = "blue", linetype = "dashed", linewidth = 0.8) +
    ggplot2::geom_text(
      data = mean_vals,
      ggplot2::aes(x = Inf, y = Inf, label = paste("Mean:", value)),
      hjust = 1.15, vjust = 1.4, size = 3, colour = "black",
      inherit.aes = FALSE) +
    ggplot2::facet_wrap(~ variable, scales = "free", ncol = ncol_val) +
    ggplot2::labs(x = "Local q-value", y = "Density") +
    facet_theme()

  ggplot2::ggsave(fname, p, width = width, height = height, dpi = 300)
  cat("Saved:", fname, "\n")
  invisible(p)
}

# Interaction pair category: "X x X", "X x GC", "GC x GC".
vartype_of <- function(var1, var2, x_vars, gc_names) {
  ifelse(var1 %in% x_vars & var2 %in% x_vars, "X × X",
         ifelse(var1 %in% gc_names & var2 %in% gc_names,
                "GC × GC", "X × GC"))
}

interaction_palette <- function(cfg) {
  setNames(
    c(cfg$figures$interaction_colors$xx,
      cfg$figures$interaction_colors$xgc,
      cfg$figures$interaction_colors$gcgc),
    c("X × X", "X × GC", "GC × GC")
  )
}