Reproducing Local Pathways of Association

1 Method Overview & Reproduction Scope

1.1 Core Idea

Path analysis answers a question correlation cannot: not whether climate and vegetation move together, but how much of that movement runs through soil. A structural equation model estimates the strength of each arrow in a hypothesised causal diagram, and reports it as a path coefficient λ.

Over a study area the size of the Tibetan Plateau, one λ per arrow is a claim that the same causal chain operates the same way at every location — from the monsoon-fed southeastern gorges to the cold desert of the northwest, 3,000 km away. LPA drops that assumption. It fits the same structural equation model separately inside a neighbourhood around every location, so each arrow becomes a surface rather than a scalar, and the variation of that surface becomes the result.

• Research problem: path coefficients estimated globally absorb regional differences into an average, and the average can be an artefact — the paper's own λ for climate on soil runs from −0.969 to 0.997 across space while the global model reports a single 0.677.
• One-line contribution: LPA determines a local range from the data with Ripley's K, estimates the structural equation model inside that range at every location, and maps the resulting path coefficients and their significance.
• Why it matters: a global coefficient tells you a relationship exists somewhere; a local one tells you where, and where it reverses. On this case, three of the seven pathways are significant at fewer than 71 % of locations, which a single starred coefficient in a conventional SEM table would never reveal.

Published reference · LPA

1.2 Method Logic

Key concepts and equations

• Ripley's K counts neighbours at a distance. For a distance r, K(r) is the expected number of further points within r of a typical point, scaled by the study area A and the point count n:
  K(r) = (A/n²) Σ_i Σ_j≠i I(d_ij ≤ r)
  Under complete spatial randomness K(r) = πr². Anything above that is clustering.
• Besag's L makes it readable, and its peak is the local range. The square-root transform turns the quadratic expectation into a straight line at zero, so the curve can be read directly:
  L(r) = √(K(r)/π) − r  r_opt = arg max_r L(r)
  L(r) > 0 means the points are more clustered than random at that distance, and the distance at which L peaks is the scale at which spatial interaction is strongest. That distance becomes the radius of every local model.
• The local structural equation model. At location u, only the observations within d_a(X_i) of u enter the fit:
  {X_i(v) : ‖u − v‖ ≤ d_a(X_i)}
  Inside that window the usual path equations apply, with every coefficient now carrying the location as an argument:
  X₂ = λ₂₁(u) X₁ + ε₁  X₃ = λ₃₂(u) X₂ + λ₃₁(u) X₁ + ε₂
  Estimation is ordinary maximum likelihood, λ̂ = arg min_λ Σ(y_i − ŷ_i)². Nothing about the estimator is new — the locality lives entirely in which rows are passed to it.
• Across space, λ becomes a function of position. Repeating the fit at u₁ … u_n yields λ_ij(u) = f(u), a surface that can be mapped, summarised and significance-tested like any other spatial variable.

1.3 Reproduction Scope

This page reproduces two different things, and keeping them apart is the point.

How to use this page

Read §2 to get the pipeline running, then §3.2 to confirm it reproduces the published table before trusting anything else it prints. §3.3 walks the method one step at a time, each with the code, the console output it produces and the figure it draws. §4 is for porting the model to another study area.

2 Setup, Structure & Data Contract

2.1 Environment & Dependencies

Two packages carry the method: spatstat.explore for Ripley's K and lavaan for the structural equation model. Neither is called through a wrapper — LPA is an arrangement of standard tools rather than a new estimator, which is why it ports easily to other disciplines.

# The two that do the work
install.packages(c("lavaan", "spatstat.explore", "spatstat.geom"))

# Figures
install.packages(c("ggplot2", "patchwork"))

2.2 Project Structure & Configuration

LPA/
├── run-all.R                        one command runs everything
├── config/project-config.R          the only file you edit to port the model
├── R/
│   ├── 00-config.R  01-helpers.R    paths, logging, neighbourhoods, the local fit
│   ├── 10-prepare-data.R            read, project, describe
│   ├── 20-optimal-range.R           Ripley's K -> r_opt          (Step 1)
│   ├── 30-local-sem.R               one SEM per location         (Steps 2-3)
│   ├── 40-reproduce-published.R     rebuild Table 2, then test the refit
│   ├── 50-radius-sensitivity.R      how much the radius decides the answer
│   ├── 90-tables.R                  LaTeX tables
│   └── p01..p04                     figures
├── data/
│   ├── lpa-tibetan-plateau.csv      1,801 locations, the published case
│   ├── LSEM_out.xlsx                the same data as the authors distributed it
│   ├── original-violin-plot.R       the authors' own plotting script, unedited
│   └── derived/                     points.rds, r-opt.rds, neighbour counts
├── results/  tables/  figs/         everything the pipeline writes
├── tests/                           31 assertions against the paper, 10 on the method
├── env/                             session info, package versions
└── assets/                          reference figures from the article

Every choice the method makes is a named constant in one file. Porting the model to another study area means editing config/project-config.R and nothing else.

# -- Input ---------------------------------------------------------------
INPUT_FILE <- "data/lpa-tibetan-plateau.csv"
COORD_X    <- "x"          # longitude, decimal degrees
COORD_Y    <- "y"          # latitude,  decimal degrees

# -- The optimal local range ---------------------------------------------
# Edge correction matters more than anything else here: with no correction
# the L curve peaks at roughly half the published distance.
K_CORRECTION <- "border"       # "border", "isotropic" or "none"
K_WINDOW     <- "rectangle"    # bounding box, or "hull"
K_R_MAX      <- 1200           # km, upper end of the search grid
MIN_LOCAL_N  <- 60             # a fit needs enough neighbours to be worth trusting

RADIUS_OVERRIDE    <- NULL     # set a number to force a radius
RADIUS_SENSITIVITY <- c(400, 500, 600, 705.29, 800, 900)

# -- The structural equation model ---------------------------------------
SEM_MODEL <- '
  Plant   =~ EVI + NDVI
  Climate =~ T + LST + P
  Soil    =~ PH + watercontent

  Soil  ~ Climate
  Plant ~ Climate + Soil
'

# -- Reporting conventions -----------------------------------------------
LAMBDA_PLOT_RANGE <- c(-1, 1)  # the filter Table 2 and Fig. 5 apply
SIG_LEVEL         <- 0.05

2.3 Data Contract

One table, one row per location. Columns 1–2 are the coordinates, columns 3–11 the observed variables, and the remaining 14 are the authors' own output — seven path coefficients and their p-values — which is what makes an exact reproduction of Table 2 possible at all.

Every observed variable arrives standardised — mean 0, standard deviation 1 to machine precision — so no scaling happens in the pipeline, and a coefficient outside [−1, 1] is a property of the fit rather than of the units.

Published reference · LPA

2.4 The Model Specification, and Why It Is a Reconstruction

Fig. 5 of the article draws the model: four measurement arrows into three latent variables, and three structural arrows between them. Read literally, it does not identify.

Published reference · LPA

A latent variable measured by a single indicator — Plant by NDVI alone, Soil by water content alone — has its loading fixed at 1 by convention, because there is nothing to estimate it against. The published output contradicts that in two places at once: the NDVI loading varies from 0.004 to 1.000 across locations, and the water-content loading is negative at most of them, averaging −0.763. Fitting the diagram as drawn confirms it — both loadings come back as exactly 1.000 everywhere, and the structural coefficients blow up past ±900 as the optimiser works on an underdetermined problem.

The data file resolves it. It carries three variables no arrow in Fig. 5 uses — EVI, T and PH — and giving each latent a second indicator makes every loading free, negative loadings possible, and the whole model identified:

SEM_MODEL <- '
  Plant   =~ EVI + NDVI          # EVI is the reference indicator
  Climate =~ T + LST + P         # T  is the reference indicator
  Soil    =~ PH + watercontent   # PH is the reference indicator

  Soil  ~ Climate                # lambda_7
  Plant ~ Climate + Soil         # lambda_6, lambda_5
'

The reference indicators are what the published signs imply. With temperature anchoring Climate, land surface temperature loads positively (published mean +0.798) and precipitation negatively (−0.370), which is the pattern of a cold arid plateau. With pH anchoring Soil, water content loads negatively (−0.763). Fitted this way, the model returns coefficients of the right sign and the right magnitude at most locations.

3 Reproduction Pipeline, Results & Validation

3.1 Pipeline Execution

cd LPA                     # this folder
Rscript run-all.R test     # 1. check against the published values (2 s)
Rscript run-all.R          # 2. full pipeline (~11 min)

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

Published reference · LPA

3.2 Reproduce the Published Values

Table 2 — reproduced exactly, once two conventions are right

Table 2 of the article validates LPA by putting the local coefficients beside the global ones. Every figure it prints follows from the distributed output file, but only under two conventions the table itself does not state, and both have to be right at once or nothing matches.

• The summaries cover only λ inside [−1, 1]. This is the filter in the authors' own plotting script. Applied to the raw columns instead, the mean for the climate-to-plant path is +3.57 rather than the published +0.335.
• The share significant is taken over the 1,283 locations with a complete set of p-values, not over all 1,801 and not over the filtered subset. Using the filtered subset instead moves four of the seven percentages off the printed value.

What the global column already tells you

Read the table's last two columns together and the argument for the method appears without any map. Every global coefficient is significant at p < 0.01, yet the local range around it crosses zero on six of the seven pathways. The global −0.528 for soil on plant is the average of a surface running from −0.998 to +0.774, and the pathway is significant at only 70 % of locations. The global estimate is not wrong; it is a summary of something that has no single value.

3.3 Core Analytical Steps

Five stages: purpose → code → output → result → how to read it. Step numbers match the script numbers (Step 1 is R/20-optimal-range.R).

Step 1

Identify the optimal local range

What this step does

Ripley's K counts how many neighbours a typical location has within each distance r, and compares that count to what complete spatial randomness would give. Besag's L rescales the comparison so it can be read against zero, and the distance at which L peaks is taken as the range over which spatial interaction is strongest. That distance becomes the radius of every local model in Steps 2 and 3, which makes this the single most consequential number in the pipeline.

Code

R/01-helpers.R · the estimator Copy

ripley_L <- function(d, correction = K_CORRECTION, window = K_WINDOW,
                     r_max = K_R_MAX, r_step = K_R_STEP) {
  W <- if (identical(window, "hull"))
    spatstat.geom::convexhull.xy(d$px, d$py)
  else
    spatstat.geom::owin(range(d$px), range(d$py))
  pp <- spatstat.geom::ppp(d$px, d$py, window = W, checkdup = FALSE)
  K  <- spatstat.explore::Kest(pp, correction = correction,
                               r = seq(0, r_max, by = r_step))
  col <- setdiff(names(K), c("r", "theo"))[1]
  L   <- sqrt(K[[col]] / pi) - K$r          # Besag's L
  i   <- which.max(L)                       # the optimal local range
  list(r = K$r, K = K[[col]], theo = K$theo, L = L,
       r_opt = K$r[i], L_max = L[i], correction = correction)
}

Example output

== Step 20 — optimal local range ========================= L peaks at r = 705.29 km (correction 'border', rectangle window), L_max = 128.94 L is within 5 of its maximum from 687 to 723 km — the peak is a plateau 36 km wide the border correction stops being reliable beyond 751 km wrote results/ripley-k-curve.csv (1201 rows) wrote results/optimal-range.csv (4 rows) correction window r_opt_km L_max published_km difference_km selected none rectangle 378.51 62.990 707.29 -328.78 FALSE border rectangle 705.29 128.942 707.29 -2.00 TRUE isotropic rectangle 805.99 98.027 707.29 98.70 FALSE border hull 568.03 48.920 707.29 -139.26 FALSE plateau_from_km plateau_to_km NA NA 686.95 723.18 NA NA NA NA local range carried forward: 705.29 km (published 707.29 km, -2.00 km) neighbourhood size: median 596, range 163-909; 0 locations below the minimum of 60

Demo run

Published reference · LPA

How to read this result

The reproduction lands 2 km from the published value, a difference of 0.3 %, and that residual is a projection choice the article does not state: scaling longitude by cos(latitude) gives 705.29 km, treating both axes at 111.32 km per degree gives 710.11 km, and the published number sits between them.

Agreeing to 2 km is less impressive than it sounds, and the shaded band in the figure is why. L stays within 5 units of its maximum from 687 to 723 km, so the curve has a plateau rather than a spike, and any radius in that 36 km window is an equally good reading of the same data. The honest statement is that the published range and this one fall on the same plateau — not that the third decimal was recovered.

◆Why the curve is cut short

The figure stops at 748 km although the search grid runs to 1,200. The border correction only counts points whose whole neighbourhood of radius r lies inside the window, and by 750 km almost none do: K falls from 2.28 million to 1.54 million in a single step and L drops to −50. That last stretch is the estimator running out of data, not the plateau ending, so trim_collapse() removes the collapsing tail before plotting. The peak is unaffected — it sits 43 km before the cut — but a figure drawn without the trim shows one cliff and nothing else.

!The edge correction is the result

Without a correction, pairs that fall outside the study window are simply never counted, K is biased downwards, and the peak lands at 378 km — a local model built on that radius would use a quarter of the data and answer a different question. Ripley's isotropic correction gives 806 km. The convex hull as window gives 568 km. Only the reduced-sample border correction on the bounding box reproduces the article. Report which one you used; the number means little without it.

Note also where the estimator stops: the border correction cannot estimate K beyond 750 km for this window, so the peak at 705 km sits 44 km inside the limit of what is observable. tests/test-02 asserts that margin, because a peak resting on the last estimable distance would be an artefact of the estimator rather than a feature of the pattern.

Demo run

# Every location within radius r of location i, itself included.
# This one line is what makes a path coefficient local.
neighbours_within <- function(d, i, r) {
  which(sqrt((d$px - d$px[i])^2 + (d$py - d$py[i])^2) <= r)
}

fit_local_sem <- function(dat, model = SEM_MODEL, paths = PATHS) {
  fit <- try(suppressWarnings(lavaan::sem(model, data = dat)), silent = TRUE)
  if (inherits(fit, "try-error")) return(...)            # keep NA, record why
  if (!lavaan::lavInspect(fit, "converged")) return(...)  # do not invent a number
  est <- lavaan::standardizedSolution(fit)
  # pull the seven reported pathways by their lavaan names
  mapply(function(l, o, r) est$est.std[est$lhs == l & est$op == o & est$rhs == r],
         paths$lhs, paths$op, paths$rhs)
}

Step 3

Local pathway analysis across space

What this step does

Repeat Step 2 at every location. Each of the seven path coefficients stops being a scalar and becomes a surface, λ_ij(u), which can be mapped and summarised like any other spatial variable.

Code

R/30-local-sem.R · the loop Copy

for (k in seq_along(idx)) {
  i   <- idx[k]
  sub <- d[neighbours_within(d, i, R_OPT), SEM_VARS, drop = FALSE]
  sub <- sub[stats::complete.cases(sub), , drop = FALSE]
  if (nrow(sub) < MIN_LOCAL_N) {                 # too few neighbours to trust
    diag_n[k] <- nrow(sub); diag_note[k] <- "too few neighbours"; next
  }
  f <- fit_local_sem(sub)
  L[k, ] <- f$lambda; P[k, ] <- f$p              # seven coefficients, seven p-values
  diag_n[k] <- f$n; diag_ok[k] <- f$converged
}

Example output

lambda path kind n_compared correlation sign_agreement lambda1 Precipitation =~ Climate Measurement 1612 0.779 93.1 lambda2 LST =~ Climate Measurement 306 0.826 93.5 lambda3 Water content =~ Soil Measurement 1362 0.310 98.7 lambda4 NDVI =~ Plant Measurement 380 0.714 100.0 lambda5 Soil ~ Plant Structural 731 0.374 76.7 lambda6 Climate ~ Plant Structural 1061 0.489 82.2 lambda7 Climate ~ Soil Structural 1023 -0.128 72.4

Demo run

Demo run

How to read this result

The refit reproduces the measurement half of the model well and the structural half poorly. The four measurement pathways agree in sign at 93–100 % of locations and correlate 0.31 to 0.83; the three structural pathways agree in sign at 72–82 % and one of them, climate on soil, does not correlate at all. Weighted across every comparable location, sign agreement is 87.8 %.

!What that split means

Loadings are pinned by the correlation between an indicator and its latent, so a measurement path is robust to how the rest of the model is written. A structural path between two latents depends on how both latents were constructed, so any difference in the specification propagates straight into it. The pattern in the table is what you would expect if the reconstruction in §2.4 has the right variables and not quite the right model — which is exactly the claim being made for it. Read the structural surfaces on this page as the published ones (they are, in fig03); read the refit as evidence that the method behaves as described, not as an independent confirmation of the estimates.

Demo run

Step 4

Pathway assessment and significance

What this step does

Three questions of the seven surfaces: what does the distribution of each coefficient look like, where is it significant, and how much of the study area does it actually describe.

Code

R/01-helpers.R · the reporting conventions Copy

# Table 2's summary: mean [min, max] over the plotted range, and the share of
# locations whose p-value clears SIG_LEVEL — over a different denominator.
summarise_lambda <- function(lambda, p, range = LAMBDA_PLOT_RANGE,
                             alpha = SIG_LEVEL) {
  v  <- lambda[is.finite(lambda) & lambda >= range[1] & lambda <= range[2]]
  pv <- p[is.finite(p)]
  data.frame(n_in_range = length(v), mean = mean(v), min = min(v), max = max(v),
             n_p = length(pv), pct_significant = 100 * mean(pv < alpha))
}

Example output

lambda path n n_kept pct_kept min max lambda1 Precipitation =~ Climate 1732 1732 100.00 -9.665e-01 0.7466 lambda2 LST =~ Climate 1732 1212 69.98 -6.052e-01 21.8025 lambda3 Water content =~ Soil 1475 1472 99.80 -1.364e+00 0.4304 lambda4 NDVI =~ Plant 1801 1015 56.36 3.525e-03 10.3234 lambda5 Soil ~ Plant 1475 906 61.42 -5.439e+02 1564.6905 lambda6 Climate ~ Plant 1732 1164 67.21 -6.477e+02 3910.1002 lambda7 Climate ~ Soil 1473 1145 77.73 -1.353e+02 38.0763

Demo run

Demo run

Demo run

Demo run

How to read this result

Two of these panels carry a warning rather than a finding. The maps in fig03 and fig04 look complete, but λ₄ is drawn at 56 % of locations and λ₅ at 50 % — the rest are either missing or outside the plotted range. A reader of the published figure has no way to tell the difference between a location that was estimated as unremarkable and one that was dropped, so results/map-coverage.csv reports the coverage of every panel explicitly.

RADIUS_SENSITIVITY <- c(400, 500, 600, 705.29, 800, 900)

# Fitting every location at six radii would take the best part of an hour,
# so the check runs on a fixed random subsample.
for (r in RADIUS_SENSITIVITY) {
  for (k in seq_along(idx)) {
    sub <- d[neighbours_within(d, idx[k], r), SEM_VARS, drop = FALSE]
    f <- fit_local_sem(sub[stats::complete.cases(sub), , drop = FALSE])
    L[k, ] <- f$lambda
  }
}

3.4 Reading the Results Together

Demo run

The three findings above point the same way. LPA earns its keep not by producing a better number than a global SEM but by showing that the number was never single — and the same evidence that makes that case also limits how far any one local coefficient should be pushed. A coefficient estimated from 596 overlapping neighbours, significant at 62 % of locations, drawn on 64 % of the study area and moving by 0.16 when the radius changes is good evidence of spatial heterogeneity in the climate-to-soil pathway. It is weak evidence about any particular valley.

4 Adaptation, Writing & Reproducibility

4.1 Port to Your Domain

LPA needs point locations, a set of observed variables and a causal diagram worth testing. It does not care whether the variables are ecological — the same structure applies to house prices and amenities, or to exposure, deprivation and health.

What to edit

4.2 Write the Paper

Manuscript skeleton

Manuscript outline

1  Introduction       — associations between variables are assumed constant over
                         space; path analysis is global; the contribution (3 aims)
2  LPA model          — 2.1 The optimal local range (Ripley's K, Besag's L)
                         2.2 The local structural equation model
                         2.3 Estimation at every location; the surfaces
                         2.4 Assessment: significance, convergence, coverage
3  Case & data        — study area, the observed variables, the causal diagram
4  Results            — 4.1 the local range and what fixes it
                         4.2 the seven pathway surfaces and their significance
                         4.3 local against global coefficients
                         4.4 sensitivity to the radius
5  Discussion         — what varies in space and what does not; the identification
                         limits of the reconstructed specification
6  Conclusions        — the local pathways, and the radius they depend on

Questions a reviewer will ask

Why that radius? → results/optimal-range.csv shows the estimate and the three alternatives. Would a different radius change the conclusion? → tables/table-radius-sensitivity.tex; on this case the signs hold and two magnitudes move. How many local models failed? → 134 of 1,801 here, in results/local-fit-diagnostics.csv. Is the map showing every location? → results/map-coverage.csv, and no. Is this just a global model in disguise? → tests/test-02 checks that 59 % of locations sit more than 0.1 away from the global estimate.

4.3 Final Reproducibility Package

References & credits

Demo outputs were generated by the scripts in this folder with R 4.6.0, lavaan 0.7.2 and spatstat.explore 3.8.1. Reference figures are reproduced from the article, which is open access under CC BY, for comparison. The data file and the violin-plot script were provided by the authors.

1. Hu J, Qu R, Song Y, Wu P (2025). Local pathways of association. International Journal of Applied Earth Observation and Geoinformation 139:104531. doi:10.1016/j.jag.2025.104531 · PDF (source of the model, the case study, the Fig. 1, Fig. 4 and Fig. 5 reference images, and of every value locked in tests/test-01)
2. Rosseel Y (2012). lavaan: an R package for structural equation modeling. Journal of Statistical Software 48(2):1–36. doi:10.18637/jss.v048.i02 · CRAN (the estimator behind every local fit)
3. Baddeley A, Rubak E, Turner R (2015). Spatial Point Patterns: Methodology and Applications with R. Chapman and Hall/CRC. spatstat (Ripley's K, and the edge corrections that decide the local range)
4. Ripley BD (1977). Modelling spatial patterns. Journal of the Royal Statistical Society B 39(2):172–192. doi:10.1111/j.2517-6161.1977.tb01615.x (the K function)
5. Besag J (1977). Contribution to the discussion of Dr Ripley's paper. Journal of the Royal Statistical Society B 39(2):193–195. (the L transform whose peak defines the local range)
6. Anselin L (1995). Local indicators of spatial association — LISA. Geographical Analysis 27(2):93–115. doi:10.1111/j.1538-4632.1995.tb00338.x (the local-statistics tradition LPA extends to path analysis)