Checking early warning signals

early warning signals

time series

resilience

ecology tutorial

Surrogate tests in R, false-positive rates and why early warning signals are neither necessary nor sufficient: a rising signal is a hypothesis, not a forecast.

Author

Tidy Ecology

Published

2026-06-06

Modified

2026-10-01

Corrected 1 October 2026: the post said the fast-approach results came from thirty noise seeds; it now says they come from the seeds that collapse within the record, since the rest do not tip before it ends.

Corrected 27 September 2026: the post called early warning signals necessary but not sufficient evidence of instability; its own fast-approach runs show a fold with no warning, so it now says they are neither necessary nor sufficient.

The earlier tutorials in this series showed indicators that rise before a fold: lag-1 autocorrelation and variance in time, the same on residuals under different analyst choices, and variance and Moran’s I in space. Each ended on the same warning: a rising indicator is a hypothesis, not a forecast. This post does the reckoning. It asks three questions of any early warning claim. Is the rise distinguishable from chance? How often does the test cry wolf when nothing is happening? And can the indicator rise without a transition, or a transition happen without the indicator? The answers are sobering, and they are the reason to treat these tools as a screen rather than a prediction.

The shared machinery

We reuse the collapsing grazing series and the same rolling-window indicators. The lag-1 autocorrelation is computed with the direct order-1 estimator, which matches ar.ols with intercept = FALSE exactly and is quicker inside a Monte Carlo loop.

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

gauss_detrend <- function(y, bw){ tg <- seq_along(y)
  trend <- sapply(tg, function(i){ w <- dnorm(tg, i, bw); sum(w*y)/sum(w) }); y - trend }
ar1_fast <- function(w){ w <- w - mean(w); L <- length(w)
  sum(w[2:L]*w[1:(L-1)]) / sum(w[1:(L-1)]^2) }              # == ar.ols(order 1, no intercept)
roll_ar1 <- function(res, wf = 0.5){ win <- round(wf*length(res)); s <- 1:(length(res)-win+1)
  sapply(s, function(k) ar1_fast(res[k:(k+win-1)])) }
roll_var <- function(res, wf = 0.5){ win <- round(wf*length(res)); s <- 1:(length(res)-win+1)
  sapply(s, function(k) var(res[k:(k+win-1)])) }
kt <- function(v) cor(seq_along(v), v, method = "kendall")

sim_grazing <- function(seed, T = 600, dt = 0.01, sigma = 0.30, c0 = 1.0, c1 = 2.72, x0 = 8.9){
  set.seed(seed); ns <- T/dt; cvec <- seq(c0, c1, length.out = ns); x <- numeric(ns); x[1] <- x0
  for(i in 2:ns){ x[i] <- x[i-1] + f(x[i-1], cvec[i-1])*dt + sigma*sqrt(dt)*rnorm(1)
    if(x[i] < 0.001) x[i] <- 0.001 }
  keep <- seq(1, ns, by = 1/dt); data.frame(t = (keep-1)*dt, x = x[keep], c = cvec[keep]) }

d <- sim_grazing(seed = 4238)
tr_idx <- which(d$x < 3.5)[1]
seg <- d$x[1:(tr_idx - 10)]; n <- length(seg)
res <- gauss_detrend(seg, 0.10*n)
tau_obs <- kt(roll_ar1(res))

The observed autocorrelation trend statistic on this known approach to a fold is 0.69. The question is what that number is worth.

Check one: is the rise significant?

A positive Kendall tau is only evidence against a null of no change. We build that null with surrogates: series that share the record’s short-term correlation structure but carry no trend in it. An autoregressive bootstrap fits an order-1 model to the residuals and resamples its innovations. A phase-randomised surrogate keeps the power spectrum, and so the autocorrelation, but scrambles the phases. For each surrogate we recompute the trend statistic; the observed value’s position in that distribution is a p-value.

arboot_taus <- function(res, B, seed){ set.seed(seed); phi <- ar1_fast(res)
  e <- res[-1] - phi*res[-length(res)]; e <- e - mean(e); nn <- length(res)
  replicate(B, { innov <- sample(e, nn, replace = TRUE); y <- numeric(nn); y[1] <- innov[1]
    for(i in 2:nn) y[i] <- phi*y[i-1] + innov[i]; kt(roll_ar1(y)) }) }
phase_taus <- function(y, B, seed){ set.seed(seed); nn <- length(y); amp <- Mod(fft(y))
  replicate(B, { ph <- runif((nn - 1) %/% 2, 0, 2*pi)
    phases <- c(0, ph, if(nn %% 2 == 0) 0 else NULL, -rev(ph))
    s <- Re(fft(amp*exp(1i*phases), inverse = TRUE)/nn); kt(roll_ar1(s)) }) }

ta <- arboot_taus(res, 500, 71)
tp <- phase_taus(res, 500, 72)
p_ar <- (1 + sum(ta >= tau_obs)) / 501
p_ph <- (1 + sum(tp >= tau_obs)) / 501
q95_ar <- quantile(ta, 0.95)

Here is the result. The autoregressive-bootstrap null puts its 95th percentile at 0.72, against the observed 0.69. The p-value is 0.068; the phase-randomised surrogates give 0.080. Compare either p-value with the conventional 0.05 before calling the rise a warning, because a series built to approach a fold, with a clean rise in autocorrelation, still has to clear a stationary null. The null distribution of the trend statistic is broad: rolling-window estimates from overlapping windows are strongly autocorrelated, so they drift up or down by chance far more than independent points would.

A histogram of surrogate Kendall tau values spread from about minus 0.9 to 0.9, centred near zero. A solid line marks the observed value at 0.69, sitting left of a dashed line at the 95th percentile near 0.72. — Figure 1: Surrogate null distribution of the autocorrelation trend statistic (autoregressive bootstrap). The solid line is the observed value, the dashed line the 95th percentile of the surrogates, and the surrogates come from a stationary process with the same short-term structure.

Check two: how often does it cry wolf?

If the null distribution is that broad, a naive test will over-declare warnings. We simulate many stationary records, the grazing model held at a constant c = 2.0 on the vegetated branch, with no approaching fold, and apply two tests to each: the naive Kendall test, which treats the window estimates as independent, and the surrogate test from above.

sim_stat <- function(seed, cc = 2.0, n = 200, dt = 0.05, sigma = 0.30, x0 = 7.32){
  set.seed(seed); ns <- n/dt; x <- numeric(ns); x[1] <- x0
  for(i in 2:ns){ x[i] <- x[i-1] + f(x[i-1], cc)*dt + sigma*sqrt(dt)*rnorm(1)
    if(x[i] < 0.001) x[i] <- 0.001 }
  x[seq(1, ns, by = 1/dt)] }

M <- 300; nb <- 200; naive_rej <- 0; surro_rej <- 0; tau_null <- numeric(M)
for(m in 1:M){
  y <- sim_stat(5000 + m); rr <- gauss_detrend(y, 0.10*nb)
  a <- roll_ar1(rr); to <- kt(a); tau_null[m] <- to
  pn <- cor.test(seq_along(a), a, method = "kendall")$p.value
  if(!is.na(pn) && pn < 0.05 && to > 0) naive_rej <- naive_rej + 1
  ps <- (1 + sum(arboot_taus(rr, 150, 9000 + m) >= to)) / 151
  if(ps < 0.05) surro_rej <- surro_rej + 1
}
naive_fp <- naive_rej/M; surro_fp <- surro_rej/M
sd_null <- sd(tau_null); frac_half <- mean(tau_null > 0.5)

On stationary data with nothing approaching, the trend statistic has a standard deviation of 0.47 and exceeds 0.5 in 21 per cent of runs by chance alone. The naive test declares a significant rising warning in 44 per cent of these null records: it cries wolf almost half the time. The surrogate test brings the false-positive rate back to 5 per cent, close to the nominal five. The lesson is blunt: a positive trend statistic assessed with an off-the-shelf significance test is nearly meaningless, and any observed value has to be read against this spread rather than against zero.

Two panels. Left: a histogram of Kendall tau on null data, spread from about minus 0.9 to 0.9 around zero, with the observed 0.69 marked to the right of centre. Right: two bars, the naive test near 0.44 and the surrogate test near 0.05, with a dashed line at 0.05. — Figure 2: Left: the trend statistic on stationary null records, centred at zero but very wide, with the observed value marked. Right: false-positive rate of the naive test versus the surrogate test; the dashed line is the nominal five per cent.

Check three: neither necessary nor sufficient

Even a correct, significant rise is not proof of an approaching fold, and a fold can arrive with no rise at all. Two simulations make the point. First a false alarm: hold c constant, so there is no bifurcation anywhere, but let the process noise grow through the record. Nothing is destabilising, yet the variance indicator climbs.

sim_risingnoise <- function(seed, cc = 2.0, n = 300, dt = 0.05, s0 = 0.12, s1 = 0.55, x0 = 7.32){
  set.seed(seed); ns <- n/dt; sv <- seq(s0, s1, length.out = ns); x <- numeric(ns); x[1] <- x0
  for(i in 2:ns){ x[i] <- x[i-1] + f(x[i-1], cc)*dt + sv[i-1]*sqrt(dt)*rnorm(1)
    if(x[i] < 0.001) x[i] <- 0.001 }
  x[seq(1, ns, by = 1/dt)] }
y_fa <- sim_risingnoise(seed = 4241); r_fa <- gauss_detrend(y_fa, 0.10*length(y_fa))
tau_var_fa <- kt(roll_var(r_fa)); tau_ar_fa <- kt(roll_ar1(r_fa))

With no bifurcation in sight, the variance trend statistic is +0.98. It rises because these indicators respond to any change in the fluctuation structure, and a noisier environment is not an approaching fold.

Now the opposite failure: a real fold that arrives too fast to announce itself. Critical slowing needs a slow, quasi-static approach; when the driver is pushed quickly, the system tips through the fold before the indicators develop. We ramp the same driver over a much shorter record, run thirty noise seeds, and keep the runs that collapse before the record ends.

sim_fast <- function(seed, T = 90, dt = 0.01, sigma = 0.30, c0 = 1.0, c1 = 2.95, x0 = 8.9){
  set.seed(seed); ns <- T/dt; cvec <- seq(c0, c1, length.out = ns); x <- numeric(ns); x[1] <- x0
  for(i in 2:ns){ x[i] <- x[i-1] + f(x[i-1], cvec[i-1])*dt + sigma*sqrt(dt)*rnorm(1)
    if(x[i] < 0.001) x[i] <- 0.001 }
  keep <- seq(1, ns, by = 1/dt); data.frame(t = (keep-1)*dt, x = x[keep], c = cvec[keep]) }

tau_fast <- c()
for(sd in 1:30){ df <- sim_fast(sd); tp <- which(df$x < 3.5)[1]
  if(is.na(tp) || tp < 40) next; sg <- df$x[1:(tp-4)]; if(length(sg) < 50) next
  tau_fast <- c(tau_fast, kt(roll_ar1(gauss_detrend(sg, 0.10*length(sg))))) }
med_fast <- median(tau_fast); lo_fast <- min(tau_fast); hi_fast <- max(tau_fast)

Over the 24 seeds that collapse within the record, the fast approach gives a median trend statistic of 0.48, below the slow-approach 0.69, and across those seeds it ranges from -0.36 to +0.86: some runs give a negative statistic, no warning at all.

Two panels. Left: a biomass series with constant mean but visibly widening fluctuations over time. Right: a strip of points showing fast-approach Kendall tau scattered from below zero to high positive values, with a horizontal line marking the slow-approach value. — Figure 3: Left: constant driver with growing noise, no bifurcation, yet rising variance. Right: warning strength across noise seeds for a fast approach to a real fold, with the slow-approach value marked.

Using early warning signals honestly

None of this retires the method. Critical slowing is a real phenomenon, and on long, slowly-forced systems the indicators do carry information. The point is what to claim from them. Three habits follow from the checks above. Assess significance against surrogates, never against a naive test, because the naive false-positive rate here was nearly half. Report the sensitivity to detrending and window choice from the previous tutorial, since a warning that survives only one setting is a choice, not a signal. And remember the base rate. Boettiger and Hastings (2012) call the common error a prosecutor’s fallacy: studying only the systems that collapsed inflates the apparent skill, because the probability of a warning given a transition is not the probability of a transition given a warning. Early warning signals are a screening tool that narrows hypotheses. They are evidence consistent with instability, neither required for a transition nor proof of one, and they are best reported as a tested hypothesis rather than a forecast.

References

Scheffer et al. 2009 Nature 461(7260):53-59 (10.1038/nature08227)
Dakos et al. 2012 PLoS ONE 7(7):e41010 (10.1371/journal.pone.0041010)
Boettiger and Hastings 2012 Proceedings of the Royal Society B 279(1748):4734-4739 (10.1098/rspb.2012.2085)
Boettiger, Ross and Hastings 2013 Theoretical Ecology 6(3):255-264 (10.1007/s12080-013-0192-6)
Dakos et al. 2015 Philosophical Transactions of the Royal Society B 370(1659):20130263 (10.1098/rstb.2013.0263)

--- title: "Checking early warning signals" description: "Surrogate tests in R, false-positive rates and why early warning signals are neither necessary nor sufficient: a rising signal is a hypothesis, not a forecast." date: "2026-06-06 14:00" categories: [early warning signals, time series, R, resilience, ecology tutorial] image: thumbnail.png image-alt: "A histogram of a surrogate null distribution for a trend statistic, counts against the Kendall tau of the surrogates, centred near zero, with a solid line marking the observed value inside the range of the surrogates near their upper end and a few surrogate bars beyond it, and a dashed line marking the 95th percentile just to its right. Titled Observed trend against the surrogate null, with the autoregressive-bootstrap p-value as a subtitle." date-modified: "2026-10-01" --- *Corrected 1 October 2026: the post said the fast-approach results came from thirty noise seeds; it now says they come from the seeds that collapse within the record, since the rest do not tip before it ends.* *Corrected 27 September 2026: the post called early warning signals necessary but not sufficient evidence of instability; its own fast-approach runs show a fold with no warning, so it now says they are neither necessary nor sufficient.* The earlier tutorials in this series showed indicators that rise before a fold: [lag-1 autocorrelation and variance in time](../early-warning-signals-critical-slowing/), [the same on residuals under different analyst choices](../detrending-and-bandwidth-choice-ews/), and [variance and Moran's I in space](../spatial-early-warning-signals/). Each ended on the same warning: a rising indicator is a hypothesis, not a forecast. This post does the reckoning. It asks three questions of any early warning claim. Is the rise distinguishable from chance? How often does the test cry wolf when nothing is happening? And can the indicator rise without a transition, or a transition happen without the indicator? The answers are sobering, and they are the reason to treat these tools as a screen rather than a prediction. ## The shared machinery We reuse the collapsing grazing series and the same rolling-window indicators. The lag-1 autocorrelation is computed with the direct order-1 estimator, which matches `ar.ols` with `intercept = FALSE` exactly and is quicker inside a Monte Carlo loop. ```{r} #| label: machinery r <- 1; K <- 10; h <- 1 f <- function(x, c) r*x*(1 - x/K) - c*x^2/(x^2 + h^2) gauss_detrend <- function(y, bw){ tg <- seq_along(y) trend <- sapply(tg, function(i){ w <- dnorm(tg, i, bw); sum(w*y)/sum(w) }); y - trend } ar1_fast <- function(w){ w <- w - mean(w); L <- length(w) sum(w[2:L]*w[1:(L-1)]) / sum(w[1:(L-1)]^2) } # == ar.ols(order 1, no intercept) roll_ar1 <- function(res, wf = 0.5){ win <- round(wf*length(res)); s <- 1:(length(res)-win+1) sapply(s, function(k) ar1_fast(res[k:(k+win-1)])) } roll_var <- function(res, wf = 0.5){ win <- round(wf*length(res)); s <- 1:(length(res)-win+1) sapply(s, function(k) var(res[k:(k+win-1)])) } kt <- function(v) cor(seq_along(v), v, method = "kendall") sim_grazing <- function(seed, T = 600, dt = 0.01, sigma = 0.30, c0 = 1.0, c1 = 2.72, x0 = 8.9){ set.seed(seed); ns <- T/dt; cvec <- seq(c0, c1, length.out = ns); x <- numeric(ns); x[1] <- x0 for(i in 2:ns){ x[i] <- x[i-1] + f(x[i-1], cvec[i-1])*dt + sigma*sqrt(dt)*rnorm(1) if(x[i] < 0.001) x[i] <- 0.001 } keep <- seq(1, ns, by = 1/dt); data.frame(t = (keep-1)*dt, x = x[keep], c = cvec[keep]) } d <- sim_grazing(seed = 4238) tr_idx <- which(d$x < 3.5)[1] seg <- d$x[1:(tr_idx - 10)]; n <- length(seg) res <- gauss_detrend(seg, 0.10*n) tau_obs <- kt(roll_ar1(res)) ``` The observed autocorrelation trend statistic on this known approach to a fold is `r sprintf("%.2f", tau_obs)`. The question is what that number is worth. ## Check one: is the rise significant? A positive Kendall tau is only evidence against a null of no change. We build that null with surrogates: series that share the record's short-term correlation structure but carry no trend in it. An autoregressive bootstrap fits an order-1 model to the residuals and resamples its innovations. A phase-randomised surrogate keeps the power spectrum, and so the autocorrelation, but scrambles the phases. For each surrogate we recompute the trend statistic; the observed value's position in that distribution is a p-value. ```{r} #| label: surrogates arboot_taus <- function(res, B, seed){ set.seed(seed); phi <- ar1_fast(res) e <- res[-1] - phi*res[-length(res)]; e <- e - mean(e); nn <- length(res) replicate(B, { innov <- sample(e, nn, replace = TRUE); y <- numeric(nn); y[1] <- innov[1] for(i in 2:nn) y[i] <- phi*y[i-1] + innov[i]; kt(roll_ar1(y)) }) } phase_taus <- function(y, B, seed){ set.seed(seed); nn <- length(y); amp <- Mod(fft(y)) replicate(B, { ph <- runif((nn - 1) %/% 2, 0, 2*pi) phases <- c(0, ph, if(nn %% 2 == 0) 0 else NULL, -rev(ph)) s <- Re(fft(amp*exp(1i*phases), inverse = TRUE)/nn); kt(roll_ar1(s)) }) } ta <- arboot_taus(res, 500, 71) tp <- phase_taus(res, 500, 72) p_ar <- (1 + sum(ta >= tau_obs)) / 501 p_ph <- (1 + sum(tp >= tau_obs)) / 501 q95_ar <- quantile(ta, 0.95) ``` Here is the result. The autoregressive-bootstrap null puts its 95th percentile at `r sprintf("%.2f", q95_ar)`, against the observed `r sprintf("%.2f", tau_obs)`. The p-value is `r sprintf("%.3f", p_ar)`; the phase-randomised surrogates give `r sprintf("%.3f", p_ph)`. Compare either p-value with the conventional 0.05 before calling the rise a warning, because a series built to approach a fold, with a clean rise in autocorrelation, still has to clear a stationary null. The null distribution of the trend statistic is broad: rolling-window estimates from overlapping windows are strongly autocorrelated, so they drift up or down by chance far more than independent points would. ```{r} #| label: setup-theme #| echo: false #| message: false #| warning: false library(ggplot2) ink <- "#16241d"; body <- "#2c3a31"; forest <- "#275139"; label <- "#46604a" sage <- "#93a87f"; paper <- "#f5f4ee"; line <- "#dad9ca"; faint <- "#5d6b61" amber <- "#cda23f"; brick <- "#b5534e" theme_te <- function(){ theme_minimal(base_size = 11) + theme(plot.background = element_rect(fill = "#f5f4ee", colour = NA), panel.background = element_rect(fill = "#f5f4ee", colour = NA), panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = line, linewidth = 0.3), axis.title = element_text(colour = body), axis.text = element_text(colour = faint), plot.title = element_text(colour = ink, face = "bold", size = 12), plot.subtitle = element_text(colour = label, size = 10), legend.position = "none") } ``` ```{r} #| label: fig-surrogates #| echo: false #| fig-width: 7.4 #| fig-height: 3.8 #| fig-cap: "Surrogate null distribution of the autocorrelation trend statistic (autoregressive bootstrap). The solid line is the observed value, the dashed line the 95th percentile of the surrogates, and the surrogates come from a stationary process with the same short-term structure." #| fig-alt: "A histogram of surrogate Kendall tau values spread from about minus 0.9 to 0.9, centred near zero. A solid line marks the observed value at 0.69, sitting left of a dashed line at the 95th percentile near 0.72." ds <- data.frame(tau = ta) p_surr <- ggplot(ds, aes(tau)) + geom_histogram(bins = 40, fill = sage, colour = "white", alpha = 0.8) + geom_vline(xintercept = tau_obs, colour = forest, linewidth = 0.9) + geom_vline(xintercept = q95_ar, colour = brick, linetype = "dashed", linewidth = 0.6) + annotate("text", x = tau_obs, y = Inf, label = "observed ", colour = forest, hjust = 1, vjust = 2, size = 3.4) + annotate("text", x = q95_ar, y = Inf, label = " 95th pct", colour = brick, hjust = 0, vjust = 4, size = 3.2) + labs(title = "Observed trend against the surrogate null", subtitle = paste0("autoregressive-bootstrap p = ", sprintf("%.3f", p_ar)), x = "Kendall tau of surrogate", y = "count") + theme_te() p_surr ``` ```{r} #| label: thumb #| include: false ggsave("thumbnail.png", p_surr, width = 7.5, height = 5.4, dpi = 100) ``` ## Check two: how often does it cry wolf? If the null distribution is that broad, a naive test will over-declare warnings. We simulate many stationary records, the grazing model held at a constant `c` = 2.0 on the vegetated branch, with no approaching fold, and apply two tests to each: the naive Kendall test, which treats the window estimates as independent, and the surrogate test from above. ```{r} #| label: falsepos sim_stat <- function(seed, cc = 2.0, n = 200, dt = 0.05, sigma = 0.30, x0 = 7.32){ set.seed(seed); ns <- n/dt; x <- numeric(ns); x[1] <- x0 for(i in 2:ns){ x[i] <- x[i-1] + f(x[i-1], cc)*dt + sigma*sqrt(dt)*rnorm(1) if(x[i] < 0.001) x[i] <- 0.001 } x[seq(1, ns, by = 1/dt)] } M <- 300; nb <- 200; naive_rej <- 0; surro_rej <- 0; tau_null <- numeric(M) for(m in 1:M){ y <- sim_stat(5000 + m); rr <- gauss_detrend(y, 0.10*nb) a <- roll_ar1(rr); to <- kt(a); tau_null[m] <- to pn <- cor.test(seq_along(a), a, method = "kendall")$p.value if(!is.na(pn) && pn < 0.05 && to > 0) naive_rej <- naive_rej + 1 ps <- (1 + sum(arboot_taus(rr, 150, 9000 + m) >= to)) / 151 if(ps < 0.05) surro_rej <- surro_rej + 1 } naive_fp <- naive_rej/M; surro_fp <- surro_rej/M sd_null <- sd(tau_null); frac_half <- mean(tau_null > 0.5) ``` On stationary data with nothing approaching, the trend statistic has a standard deviation of `r sprintf("%.2f", sd_null)` and exceeds 0.5 in `r sprintf("%.0f", 100*frac_half)` per cent of runs by chance alone. The naive test declares a significant rising warning in `r sprintf("%.0f", 100*naive_fp)` per cent of these null records: it cries wolf almost half the time. The surrogate test brings the false-positive rate back to `r sprintf("%.0f", 100*surro_fp)` per cent, close to the nominal five. The lesson is blunt: a positive trend statistic assessed with an off-the-shelf significance test is nearly meaningless, and any observed value has to be read against this spread rather than against zero. ```{r} #| label: fig-falsepos #| echo: false #| fig-width: 8 #| fig-height: 3.6 #| fig-cap: "Left: the trend statistic on stationary null records, centred at zero but very wide, with the observed value marked. Right: false-positive rate of the naive test versus the surrogate test; the dashed line is the nominal five per cent." #| fig-alt: "Two panels. Left: a histogram of Kendall tau on null data, spread from about minus 0.9 to 0.9 around zero, with the observed 0.69 marked to the right of centre. Right: two bars, the naive test near 0.44 and the surrogate test near 0.05, with a dashed line at 0.05." p_hist <- ggplot(data.frame(tau = tau_null), aes(tau)) + geom_histogram(bins = 40, fill = sage, colour = "white", alpha = 0.8) + geom_vline(xintercept = 0.5, colour = faint, linetype = "dotted", linewidth = 0.4) + geom_vline(xintercept = tau_obs, colour = forest, linewidth = 0.9) + annotate("text", x = tau_obs, y = Inf, label = " observed", colour = forest, hjust = 0, vjust = 2, size = 3.2) + scale_x_continuous(expand = expansion(mult = c(0.05, 0.20))) + labs(title = "Trend statistic under a stationary null", x = "Kendall tau", y = "count") + theme_te() d_fp <- data.frame(test = factor(c("naive", "surrogate"), levels = c("naive", "surrogate")), fp = c(naive_fp, surro_fp)) p_bar <- ggplot(d_fp, aes(test, fp, fill = test)) + geom_col(width = 0.6) + geom_hline(yintercept = 0.05, colour = brick, linetype = "dashed", linewidth = 0.5) + scale_fill_manual(values = c("naive" = amber, "surrogate" = forest)) + labs(title = "False-positive rate", x = NULL, y = "P(declare warning | no transition)") + theme_te() grid::grid.newpage() grid::pushViewport(grid::viewport(layout = grid::grid.layout(1, 2))) print(p_hist, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 1)) print(p_bar, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 2)) ``` ## Check three: neither necessary nor sufficient Even a correct, significant rise is not proof of an approaching fold, and a fold can arrive with no rise at all. Two simulations make the point. First a false alarm: hold `c` constant, so there is no bifurcation anywhere, but let the process noise grow through the record. Nothing is destabilising, yet the variance indicator climbs. ```{r} #| label: falsealarm sim_risingnoise <- function(seed, cc = 2.0, n = 300, dt = 0.05, s0 = 0.12, s1 = 0.55, x0 = 7.32){ set.seed(seed); ns <- n/dt; sv <- seq(s0, s1, length.out = ns); x <- numeric(ns); x[1] <- x0 for(i in 2:ns){ x[i] <- x[i-1] + f(x[i-1], cc)*dt + sv[i-1]*sqrt(dt)*rnorm(1) if(x[i] < 0.001) x[i] <- 0.001 } x[seq(1, ns, by = 1/dt)] } y_fa <- sim_risingnoise(seed = 4241); r_fa <- gauss_detrend(y_fa, 0.10*length(y_fa)) tau_var_fa <- kt(roll_var(r_fa)); tau_ar_fa <- kt(roll_ar1(r_fa)) ``` With no bifurcation in sight, the variance trend statistic is `r sprintf("%+.2f", tau_var_fa)`. It rises because these indicators respond to any change in the fluctuation structure, and a noisier environment is not an approaching fold. Now the opposite failure: a real fold that arrives too fast to announce itself. Critical slowing needs a slow, quasi-static approach; when the driver is pushed quickly, the system tips through the fold before the indicators develop. We ramp the same driver over a much shorter record, run thirty noise seeds, and keep the runs that collapse before the record ends. ```{r} #| label: fastapproach sim_fast <- function(seed, T = 90, dt = 0.01, sigma = 0.30, c0 = 1.0, c1 = 2.95, x0 = 8.9){ set.seed(seed); ns <- T/dt; cvec <- seq(c0, c1, length.out = ns); x <- numeric(ns); x[1] <- x0 for(i in 2:ns){ x[i] <- x[i-1] + f(x[i-1], cvec[i-1])*dt + sigma*sqrt(dt)*rnorm(1) if(x[i] < 0.001) x[i] <- 0.001 } keep <- seq(1, ns, by = 1/dt); data.frame(t = (keep-1)*dt, x = x[keep], c = cvec[keep]) } tau_fast <- c() for(sd in 1:30){ df <- sim_fast(sd); tp <- which(df$x < 3.5)[1] if(is.na(tp) || tp < 40) next; sg <- df$x[1:(tp-4)]; if(length(sg) < 50) next tau_fast <- c(tau_fast, kt(roll_ar1(gauss_detrend(sg, 0.10*length(sg))))) } med_fast <- median(tau_fast); lo_fast <- min(tau_fast); hi_fast <- max(tau_fast) ``` Over the `r sprintf("%d", length(tau_fast))` seeds that collapse within the record, the fast approach gives a median trend statistic of `r sprintf("%.2f", med_fast)`, below the slow-approach `r sprintf("%.2f", tau_obs)`, and across those seeds it ranges from `r sprintf("%+.2f", lo_fast)` to `r sprintf("%+.2f", hi_fast)`: some runs give a negative statistic, no warning at all. ```{r} #| label: fig-necessary #| echo: false #| fig-width: 8 #| fig-height: 3.6 #| fig-cap: "Left: constant driver with growing noise, no bifurcation, yet rising variance. Right: warning strength across noise seeds for a fast approach to a real fold, with the slow-approach value marked." #| fig-alt: "Two panels. Left: a biomass series with constant mean but visibly widening fluctuations over time. Right: a strip of points showing fast-approach Kendall tau scattered from below zero to high positive values, with a horizontal line marking the slow-approach value." d_fa <- data.frame(t = seq_along(y_fa), x = y_fa) p_fa <- ggplot(d_fa, aes(t, x)) + geom_line(colour = amber, linewidth = 0.4) + labs(title = "Constant driver, growing noise", subtitle = paste0("variance tau = ", sprintf("%+.2f", tau_var_fa)), x = "time", y = "biomass x") + theme_te() d_ft <- data.frame(tau = tau_fast) p_ft <- ggplot(d_ft, aes(x = 1, y = tau)) + geom_hline(yintercept = tau_obs, colour = forest, linewidth = 0.6) + geom_hline(yintercept = 0, colour = faint, linetype = "dotted", linewidth = 0.4) + geom_jitter(width = 0.12, height = 0, colour = brick, size = 2, alpha = 0.8) + annotate("text", x = Inf, y = tau_obs, label = "slow approach", colour = forest, size = 3, hjust = 1.05, vjust = -0.5) + labs(title = "Fast approach across noise seeds", subtitle = paste0("median tau = ", sprintf("%.2f", med_fast)), x = NULL, y = "AR(1) Kendall tau") + theme_te() + theme(axis.text.x = element_blank()) grid::grid.newpage() grid::pushViewport(grid::viewport(layout = grid::grid.layout(1, 2))) print(p_fa, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 1)) print(p_ft, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 2)) ``` ## Using early warning signals honestly None of this retires the method. Critical slowing is a real phenomenon, and on long, slowly-forced systems the indicators do carry information. The point is what to claim from them. Three habits follow from the checks above. Assess significance against surrogates, never against a naive test, because the naive false-positive rate here was nearly half. Report the sensitivity to detrending and window choice from the [previous tutorial](../detrending-and-bandwidth-choice-ews/), since a warning that survives only one setting is a choice, not a signal. And remember the base rate. Boettiger and Hastings (2012) call the common error a prosecutor's fallacy: studying only the systems that collapsed inflates the apparent skill, because the probability of a warning given a transition is not the probability of a transition given a warning. Early warning signals are a screening tool that narrows hypotheses. They are evidence consistent with instability, neither required for a transition nor proof of one, and they are best reported as a tested hypothesis rather than a forecast. ## References - Scheffer et al. 2009 Nature 461(7260):53-59 (10.1038/nature08227) - Dakos et al. 2012 PLoS ONE 7(7):e41010 (10.1371/journal.pone.0041010) - Boettiger and Hastings 2012 Proceedings of the Royal Society B 279(1748):4734-4739 (10.1098/rspb.2012.2085) - Boettiger, Ross and Hastings 2013 Theoretical Ecology 6(3):255-264 (10.1007/s12080-013-0192-6) - Dakos et al. 2015 Philosophical Transactions of the Royal Society B 370(1659):20130263 (10.1098/rstb.2013.0263) ## Related tutorials - [Early warning signals and critical slowing](../early-warning-signals-critical-slowing/) - [Detrending and bandwidth in early warning signals](../detrending-and-bandwidth-choice-ews/) - [Spatial early warning signals](../spatial-early-warning-signals/) - [Interpolated cores and false slowing down](../interpolated-cores-and-false-slowing-down/)