Spatial early warning signals

Spatial variance, Moran’s I and spatial skewness measured on a coupled grazing lattice approaching a fold in R, with base stats only, and their caveats.

Critical slowing down has a spatial face. In a system where neighbouring patches interact, a slow recovery rate means a local perturbation spreads and lingers, so the correlation length grows. Snapshots of such a system, taken as a driver approaches a fold, should show rising spatial variance and rising spatial correlation even without any time series (Dakos et al. 2010; Kefi et al. 2014). A single well-sampled map can, in principle, carry the warning. This post builds a coupled grazing lattice, drives it towards collapse, and measures three proposed spatial indicators with base stats only.

A coupled grazing lattice

Each cell follows the same grazing dynamics as the temporal tutorial, with diffusive coupling to its four rook neighbours on a periodic grid, plus local noise. The grazing pressure c is held fixed within each snapshot and stepped up between snapshots.

r <- 1; K <- 10; h <- 1
f <- function(x, c) r*x*(1 - x/K) - c*x^2/(x^2 + h^2)

sim_spatial <- function(seed, L = 25, D = 0.30, sigma = 0.20, dt = 0.01,
                        c_snaps = c(1.0, 1.6, 2.0, 2.3, 2.5, 2.58), burn = 40, hold = 25){
  set.seed(seed)
  X <- matrix(8.9, L, L)
  neigh_mean <- function(M)                       # rook mean, periodic
    (M[c(L, 1:(L-1)), ] + M[c(2:L, 1), ] + M[, c(L, 1:(L-1))] + M[, c(2:L, 1)]) / 4
  step <- function(M, c){
    lap <- neigh_mean(M) - M
    M2 <- M + (f(M, c) + D*lap)*dt + sigma*sqrt(dt)*matrix(rnorm(L*L), L, L)
    M2[M2 < 0.001] <- 0.001; M2
  }
  snaps <- vector("list", length(c_snaps))
  for(s in seq_along(c_snaps)){
    cc <- c_snaps[s]
    for(i in 1:round(burn/dt)) X <- step(X, cc)   # relax towards a stationary field
    for(i in 1:round(hold/dt)) X <- step(X, cc)
    snaps[[s]] <- X
  }
  list(c = c_snaps, snaps = snaps, L = L)
}
sp <- sim_spatial(seed = 4240)

Three spatial indicators

Spatial variance is the variance across all cells. Spatial skewness is their third standardised moment. Spatial correlation we measure with Moran’s I, coded by hand with row-standardised rook weights on the periodic grid, the same construction used in the Moran’s I tutorial.

moran_grid <- function(M){
  L <- nrow(M); z <- as.vector(M) - mean(M); n <- L*L
  idx <- function(i, j) ((i-1) %% L)*L + ((j-1) %% L) + 1
  num <- 0
  for(i in 1:L) for(j in 1:L){
    nb <- c(z[idx(i-1, j)], z[idx(i+1, j)], z[idx(i, j-1)], z[idx(i, j+1)])
    num <- num + z[idx(i, j)] * sum(nb)/4          # row-standardised weights (1/4 each)
  }
  num / sum(z^2)                                    # (n / sum of weights) = 1 here
}
skew <- function(v){ v <- v - mean(v); mean(v^3) / (mean(v^2))^1.5 }

tab <- data.frame(
  c     = sp$c,
  mean  = sapply(sp$snaps, mean),
  s_var = sapply(sp$snaps, function(M) var(as.vector(M))),
  s_skew = sapply(sp$snaps, function(M) skew(as.vector(M))),
  moran = sapply(sp$snaps, moran_grid))

tau_var  <- cor(tab$c, tab$s_var,  method = "kendall")
tau_skew <- cor(tab$c, tab$s_skew, method = "kendall")
tau_moran <- cor(tab$c, tab$moran, method = "kendall")
var_fold <- tab$s_var[nrow(tab)] / tab$s_var[1]

As the driver climbs from 1.0 to 2.58, the mean biomass falls from 8.89 to 5.23, tracking the sinking upper equilibrium. Spatial variance rises from 0.020 to 0.057, a 2.9-fold increase (Kendall tau = 1.00). Moran’s I moves from 0.036 to 0.230 (Kendall tau = 1.00). These are the spatial analogue of rising variance and rising lag-1 autocorrelation in time.

Figure 1: Biomass fields at three driver values. Far from the fold the field is smooth; approaching it, correlated patches of high and low biomass emerge as the mean sinks.

Figure 2: The three spatial indicators against the driver. Spatial skewness wanders with no reliable trend.

The skewness caveat, and a deeper one

Spatial skewness is supposed to rise in magnitude near a fold, as the field starts to feel the lower state. Here it does not cooperate: its Kendall tau is 0.33, and it wanders between snapshots with no clear direction. A single snapshot is a noisy estimator of a third moment, and one realisation is not enough to pin it down. This is worth stating plainly: not every indicator that has been proposed works on a given system, and reporting only the ones that rose would be the same selection error the checking tutorial warns about.

The deeper caveat applies to spatial variance and Moran’s I as well. Spatial variance and correlation rise for reasons other than an approaching fold. Self-organised vegetation patterns, environmental gradients, dispersal limitation, and measurement grain all raise spatial correlation without any loss of resilience. A patchy map is consistent with an approaching transition, but it is equally consistent with a landscape that is simply patchy. Spatial early warning signals narrow the hypotheses; they do not confirm one. The tools that turn a rising indicator into a tested claim are the subject of the final tutorial in this series.

References

• May 1977 Nature 269(5628):471-477 (10.1038/269471a0)
• Dakos et al. 2010 Theoretical Ecology 3:163-174 (10.1007/s12080-009-0060-6)
• Dai, Korolev and Gore 2013 Nature 496(7445):355-358 (10.1038/nature12071)
• Kefi et al. 2014 PLoS ONE 9(3):e92097 (10.1371/journal.pone.0092097)

Related tutorials

• Early warning signals and critical slowing
• Detrending and bandwidth in early warning signals
• Checking early warning signals
• Spatial autocorrelation and Moran’s I