The S-map and state-dependent dynamics

EDM

time series

nonlinear dynamics

ecology tutorial

Build the S-map in base R, test whether population dynamics are state dependent, and see when its interaction coefficients are real and when they are noise.

Author

Tidy Ecology

Published

2026-06-10

Modified

2026-10-08

Corrected 8 October 2026: the post put the weak link’s failure down to its size, called the result no artefact of a hard example and said strength and identifiability are the same axis; a rerun with Y driving X as strongly as X drives Y, now in the code, still does not recover that link, so the post now says size alone does not decide. It also said the strong cross term validates beautifully against the truth, and now says it follows the truth from one of its two clusters to the other.

Corrected 1 October 2026: the post said regularised S-maps were introduced because of collinearity, but their authors introduced them because the plain S-map is sensitive to process noise; it also described the strong-link panel as one tight line of points, where it shows a tight cluster and a much looser one, and now says so.

Simplex projection answers a yes or no question: is there state-dependent structure that a linear model misses? The S-map answers the follow-up. It measures how much state dependence there is, on a continuous dial, and it hands you a local linear model at every point on the attractor.

That second part is what made the method popular in ecology. If the local model’s coefficients are the partial derivatives of the dynamics, then the coefficient on species j in species i’s equation is the interaction strength between them, at that moment, and it can change through time. Interaction strengths that fluctuate with the state of the community: that is a genuinely new kind of measurement, and it is also where the method is easiest to over-read.

This post builds the S-map, uses it as a nonlinearity test, and then shows the same estimator recovering one interaction coefficient and fabricating another, in the same system, from the same data, at the same settings.

From a global line to a local one

Fit a linear autoregressive model and you get one set of coefficients for the whole series. The system is assumed to behave the same way at low density as at high density. The S-map keeps the linear form but refits it for every point you want to predict, weighting the other points by how far they are from that point in state space:

\[w_i = \exp\!\left(-\theta \frac{d_i}{\bar{d}}\right)\]

where d_i is the distance from the target point to library point i, and dbar the mean of those distances. Then it solves a weighted least squares problem and uses those coefficients for that one prediction.

The dial is theta. At theta = 0 every point gets weight one, the weighting disappears, and the S-map is the global linear model. As theta grows, distant points fade out and the fit becomes local, so the effective model changes as the system moves around its attractor. Nothing else about the method changes: same form, same data, one number.

That construction is what makes the test work. Comparing theta = 0 against theta > 0 is not comparing two different models with different numbers of parameters; it is the same model at two settings of one dial. If skill improves as the fit localises, the dynamics depend on the state.

The engine

make_block <- function(x, E, tau = 1) {
  n <- length(x)
  tt <- seq((E - 1) * tau + 1, n)
  M <- matrix(NA_real_, nrow = length(tt), ncol = E)
  for (j in seq_len(E)) M[, j] <- x[tt - (j - 1) * tau]
  list(M = M, t = tt)
}

## S-map over an arbitrary state matrix M predicting target y
smap_block <- function(M, y, theta = 0, loo = TRUE) {
  n <- nrow(M)
  D <- as.matrix(dist(scale(M)))          # standardise: coordinates share a scale
  A0 <- cbind(1, M)
  cf <- matrix(NA_real_, n, ncol(M) + 1)
  pred <- rep(NA_real_, n)
  for (i in seq_len(n)) {
    use <- if (loo) seq_len(n) != i else rep(TRUE, n)
    d <- D[i, use]
    w <- exp(-theta * d / mean(d))
    A <- A0[use, , drop = FALSE] * w
    b <- y[use] * w
    sv <- svd(A)                          # SVD solve: weights make A ill-conditioned
    tol <- max(dim(A)) * .Machine$double.eps * max(sv$d)
    co <- as.numeric(sv$v %*% (ifelse(sv$d > tol, 1 / sv$d, 0) * crossprod(sv$u, b)))
    cf[i, ] <- co
    pred[i] <- sum(c(1, M[i, ]) * co)
  }
  list(coef = cf, pred = pred, rho = cor(pred, y))
}

smap <- function(x, E, tau = 1, tp = 1, theta = 0, loo = TRUE) {
  bl <- make_block(x, E, tau)
  ok <- (bl$t + tp) <= length(x)
  smap_block(bl$M[ok, , drop = FALSE], x[bl$t[ok] + tp], theta = theta, loo = loo)
}

Two details earn their place. Distances are computed on standardised coordinates, so a variable measured in thousands does not swamp one measured in units. And the weighted system is solved by SVD rather than solve, because at high theta most weights are near zero and the design matrix becomes ill conditioned; the singular value cutoff keeps that from blowing up silently.

Check the claim before using it

If theta = 0 really is the global linear model, then the S-map at that setting, run on the full library with no leave-one-out, must reproduce lm exactly. Worth checking rather than believing.

gen_ricker <- function(n, burn = 500, r = 3.0, x0 = 0.6) {
  N <- n + burn; x <- numeric(N); x[1] <- x0
  for (t in 1:(N - 1)) x[t + 1] <- x[t] * exp(r * (1 - x[t]))
  x[(burn + 1):N]
}

set.seed(4265)
n <- 300
z <- gen_ricker(n) * exp(rnorm(n, 0, 0.05))    # 5% observation error

s0 <- smap(z, E = 2, tp = 1, theta = 0, loo = FALSE)
X <- embed(z, 3)
m_lm <- lm(X[, 1] ~ X[, 2] + X[, 3])
coef_gap <- max(abs(sweep(s0$coef, 2, coef(m_lm))))
fit_gap <- max(abs(s0$pred - fitted(m_lm)))

The S-map’s coefficients match lm’s, and its fitted values match lm’s fitted values, to within floating-point rounding. The S-map at theta = 0 is not like a linear autoregressive model, it is one.

The nonlinearity test

Now turn the dial. The chaotic Ricker series from the previous post, and its linear counterpart: an autoregressive model fitted to it and simulated from, so the two share their autocorrelation structure and differ only in whether a deterministic map is underneath.

ar_fit <- ar(z, order.max = 10, method = "mle")
set.seed(9265)
z_lin <- as.numeric(arima.sim(n = n, model = list(ar = ar_fit$ar),
                              sd = sqrt(ar_fit$var.pred)))

thetas <- c(0, 0.1, 0.3, 0.5, 1, 2, 3, 4, 6, 8)
rho_chaos <- sapply(thetas, function(th) smap(z, E = 2, tp = 1, theta = th)$rho)
rho_lin <- sapply(thetas, function(th) smap(z_lin, E = 2, tp = 1, theta = th)$rho)

gain_chaos <- max(rho_chaos) - rho_chaos[1]
gain_lin <- max(rho_lin) - rho_lin[1]

Two line plots of forecast correlation against theta from zero to eight. The left panel for the chaotic series rises from about 0.80 to 0.99. The right panel for the linear series stays flat near 0.80 and declines slightly after theta of one. — Figure 1: Forecast skill against the S-map localisation parameter. The chaotic series improves steeply as the fit becomes local; the linear series with matching autocorrelation gains nothing and then degrades.

The chaotic series goes from rho = 0.796 at theta = 0 to 0.993 at its best, a gain of +0.197. The linear series starts at 0.796 and its best is 0.796, a gain of $3.0 \times 10^{-4}$, which is nothing next to the chaotic gain. Push theta further and the linear series gets actively worse (0.747 at theta = 8), because localising a fit that has no local structure to find just throws away data.

This is the cleanest test in the toolkit. Same model, same data, one dial, and the answer is a comparison against a value the method itself produces at theta = 0. A dashed line marks that baseline in both panels.

Interaction coefficients that move

The dial test uses one species. The coefficients need more. If you have measured several species, use their contemporaneous values as the state instead of lagged values of one, and predict one species forward. The local coefficients are then estimates of

\[\frac{\partial\, x_{i,t+1}}{\partial\, x_{j,t}}\]

at each point in time: how much species j currently affects species i.

The test system is two competitors with asymmetric coupling, and we know its Jacobian exactly. Species X drives Y with strength 0.10; Y drives X with strength 0.02. Both are chaotic, and their abundances are almost uncorrelated.

gen_coupled <- function(n, burn = 500, rx = 3.8, ry = 3.5,
                        bxy = 0.02, byx = 0.1, x0 = 0.4, y0 = 0.2) {
  N <- n + burn; x <- numeric(N); y <- numeric(N); x[1] <- x0; y[1] <- y0
  for (t in 1:(N - 1)) {
    x[t + 1] <- x[t] * (rx - rx * x[t] - bxy * y[t])
    y[t + 1] <- y[t] * (ry - ry * y[t] - byx * x[t])
  }
  data.frame(x = x[(burn + 1):N], y = y[(burn + 1):N])
}

set.seed(4265)
nc <- 400
d <- gen_coupled(nc)
M <- cbind(x = d$x[-nc], y = d$y[-nc])

## true Jacobian entries, from differentiating the model by hand
true_YX <- -0.10 * M[, "y"]                          # d y[t+1] / d x[t] : strong link
true_XY <- -0.02 * M[, "x"]                          # d x[t+1] / d y[t] : weak link
true_XX <- 3.8 - 2 * 3.8 * M[, "x"] - 0.02 * M[, "y"]

th_grid <- c(0, 0.5, 1, 2, 3, 4, 6, 8)
rho_X <- sapply(th_grid, function(th) smap_block(M, d$x[-1], theta = th)$rho)
rho_Y <- sapply(th_grid, function(th) smap_block(M, d$y[-1], theta = th)$rho)
th_best <- th_grid[which.max(rho_Y)]

fitX <- smap_block(M, d$x[-1], theta = th_best)      # x[t+1] ~ (x[t], y[t])
fitY <- smap_block(M, d$y[-1], theta = th_best)      # y[t+1] ~ (x[t], y[t])

est_YX <- fitY$coef[, 2]     # coefficient on x in the y equation : strong link
est_XY <- fitX$coef[, 3]     # coefficient on y in the x equation : weak link
est_XX <- fitX$coef[, 2]     # self effect : strong

cor_strong <- cor(est_YX, true_YX)
cor_weak <- cor(est_XY, true_XY)
cor_self <- cor(est_XX, true_XX)

## the same model with Y driving X as strongly as X drives Y, theta chosen the same way
d_eq <- gen_coupled(nc, bxy = 0.1)
M_eq <- cbind(x = d_eq$x[-nc], y = d_eq$y[-nc])
th_eq <- th_grid[which.max(sapply(th_grid, function(th) smap_block(M_eq, d_eq$y[-1], theta = th)$rho))]
cor_eq   <- cor(smap_block(M_eq, d_eq$x[-1], theta = th_eq)$coef[, 3], -0.1 * M_eq[, "x"])  # Y on X, now strong
cor_eq_y <- cor(smap_block(M_eq, d_eq$y[-1], theta = th_eq)$coef[, 2], -0.1 * M_eq[, "y"])  # X on Y, as before

Cross validation picks theta = 8, the top of the grid searched, and rho is still inching up at that edge; the forecasts are near perfect either way (rho = 0.99997 for the Y equation). Now compare each estimated coefficient against the truth we can compute analytically.

Two scatter plots of estimated coefficient against true coefficient. The left panel for the strong link shows two clusters along the one-to-one line: a tight one at the most negative true values and a looser cloud, spread well above and below the line, at the less negative ones. The right panel for the weak link shows a diffuse vertical scatter with no relationship to the truth. — Figure 2: Estimated S-map interaction coefficients against the true Jacobian entries, for the strong link and the weak link in the same system at the same settings.

The strong link is recovered: the estimated coefficients track the true ones at a correlation of 0.909, and their mean, -0.0631, sits on the true mean of -0.0637. The self effect comes back too, at 0.997. Given that the two species’ abundances correlate at only 0.11, recovering this time variation from the raw series is a real result, although most of it is Y alternating between high and low values from one step to the next; within each of the two clusters in the left panel the estimates track the variation in the truth only weakly.

The weak link is not recovered. The estimates correlate with the truth at -0.054: no relationship at all.

The part to be careful about

Look at the right panel again. Those estimates are not flat, not obviously broken, and not implausible. Their mean is -0.0160 against a true mean of -0.0127, which is close enough to pass a sanity check. Plot them against time and you get a wiggly line that looks exactly like a fluctuating interaction strength. Every diagnostic available within the analysis is green: the forecasts are near perfect, theta was chosen by cross validation, the self terms validate closely against the truth in the same fit, and the strong cross term follows the truth from one of its two clusters to the other.

And the time course is noise.

The obvious suspect is size. The weak coefficient is about 5 times smaller than the strong one and roughly two orders of magnitude smaller than the self term, which ranges over several units, and a local linear fit has to separate that contribution from everything else happening in a neighbourhood of the attractor. Size is not the whole story, though. The code block that defines gen_coupled() ends by rerunning the same model with Y driving X as strongly as X drives Y, with theta again chosen by cross validation: the coefficient on Y in the X equation is still not recovered, correlating with its truth at -0.008, while the one on X in the Y equation comes back at 0.954. In this model a coefficient as large as the recovered one can fail too, so size alone does not decide which ones a local linear fit can recover. Nothing warns you, because the quantity that fails is not the one the diagnostics measure: prediction is dominated by the large terms, so the model predicts brilliantly while the small coefficient it reports is arbitrary.

This is an honest limit of S-map interaction strengths, shown here on one simulated system, and it matters in community ecology, where most pairwise effects are weak and a few are strong. The practical reading:

The average of a coefficient can be informative even when its fluctuations are not. Here the weak link’s mean was roughly right while its time course was noise. Reporting “the interaction varied through time from A to B” is a much stronger claim than reporting its average, and it needs much more support.
Size alone does not tell you which coefficients to trust: in the rerun above a link as strong as the recovered one failed in the other equation, so a large coefficient needs the same check as a small one.
Validate on something. A simulation with a known Jacobian, a perturbation experiment, an independent estimate of the interaction. If nothing anchors the coefficient, its wiggles are decoration.
At theta = 0 the coefficients are useless even for the strong link, because a single global line through a curved map has no local derivative to report. Localisation is what makes the coefficients mean anything, and it is also what makes them variable.

There is a further wrinkle worth knowing about: with several species the local design matrix is often close to collinear, and the coefficient split between correlated species is then poorly determined even when their sum is not. Regularised versions of the S-map exist (Cenci and colleagues), introduced because the plain version is sensitive to process noise; a ridge-type penalty on the local coefficients also steadies their split when the design is near collinear.

Where to go next

Two of these coefficients said the same thing about the true structure, in opposite directions: X affects Y appreciably, Y affects X hardly at all. The S-map got the strong one right and the weak one wrong, but it needed both species measured to say anything at all.

The next question is what to do when you cannot measure both, or when you can but want to know which one is driving which. Convergent cross mapping attacks that directly, and it needs only the reconstruction from the previous post. It also inherits every fragility on display here, plus some of its own, which is what the checking post is for.

References

Ushio, Hsieh, Masuda, Deyle, Ye, Chang, Sugihara & Kondoh 2018 Nature 554(7692):360-363 (10.1038/nature25504)

Sugihara 1994 Philosophical Transactions of the Royal Society A 348(1688):477-495 (10.1098/rsta.1994.0106)

Deyle, May, Munch & Sugihara 2016 Proceedings of the Royal Society B 283(1822):20152258 (10.1098/rspb.2015.2258)

Cenci, Sugihara & Saavedra 2019 Methods in Ecology and Evolution 10(5):650-660 (10.1111/2041-210X.13150)

Hsieh, Glaser, Lucas & Sugihara 2005 Nature 435(7040):336-340 (10.1038/nature03553)

Sugihara, May, Ye, Hsieh, Deyle, Fogarty & Munch 2012 Science 338(6106):496-500 (10.1126/science.1227079)

--- title: "The S-map and state-dependent dynamics" description: "Build the S-map in base R, test whether population dynamics are state dependent, and see when its interaction coefficients are real and when they are noise." date: "2026-06-10 12:00" date-modified: "2026-10-08" categories: [R, EDM, time series, nonlinear dynamics, ecology tutorial] image: thumbnail.png image-alt: "Forecast skill rising sharply with the S-map tuning parameter for a chaotic series and staying flat for a linear series." --- *Corrected 8 October 2026: the post put the weak link's failure down to its size, called the result no artefact of a hard example and said strength and identifiability are the same axis; a rerun with Y driving X as strongly as X drives Y, now in the code, still does not recover that link, so the post now says size alone does not decide. It also said the strong cross term validates beautifully against the truth, and now says it follows the truth from one of its two clusters to the other.* *Corrected 1 October 2026: the post said regularised S-maps were introduced because of collinearity, but their authors introduced them because the plain S-map is sensitive to process noise; it also described the strong-link panel as one tight line of points, where it shows a tight cluster and a much looser one, and now says so.* [Simplex projection](../simplex-projection-and-embedding/) answers a yes or no question: is there state-dependent structure that a linear model misses? The S-map answers the follow-up. It measures *how much* state dependence there is, on a continuous dial, and it hands you a local linear model at every point on the attractor. That second part is what made the method popular in ecology. If the local model's coefficients are the partial derivatives of the dynamics, then the coefficient on species `j` in species `i`'s equation is the interaction strength between them, at that moment, and it can change through time. Interaction strengths that fluctuate with the state of the community: that is a genuinely new kind of measurement, and it is also where the method is easiest to over-read. This post builds the S-map, uses it as a nonlinearity test, and then shows the same estimator recovering one interaction coefficient and fabricating another, in the same system, from the same data, at the same settings. ## From a global line to a local one Fit a linear autoregressive model and you get one set of coefficients for the whole series. The system is assumed to behave the same way at low density as at high density. The S-map keeps the linear form but refits it for every point you want to predict, weighting the other points by how far they are from that point in state space: $$w_i = \exp\!\left(-\theta \frac{d_i}{\bar{d}}\right)$$ where `d_i` is the distance from the target point to library point `i`, and `dbar` the mean of those distances. Then it solves a weighted least squares problem and uses those coefficients for that one prediction. The dial is `theta`. At `theta = 0` every point gets weight one, the weighting disappears, and the S-map *is* the global linear model. As `theta` grows, distant points fade out and the fit becomes local, so the effective model changes as the system moves around its attractor. Nothing else about the method changes: same form, same data, one number. That construction is what makes the test work. Comparing `theta = 0` against `theta > 0` is not comparing two different models with different numbers of parameters; it is the same model at two settings of one dial. If skill improves as the fit localises, the dynamics depend on the state. ## The engine ```{r} #| label: smap-core make_block <- function(x, E, tau = 1) { n <- length(x) tt <- seq((E - 1) * tau + 1, n) M <- matrix(NA_real_, nrow = length(tt), ncol = E) for (j in seq_len(E)) M[, j] <- x[tt - (j - 1) * tau] list(M = M, t = tt) } ## S-map over an arbitrary state matrix M predicting target y smap_block <- function(M, y, theta = 0, loo = TRUE) { n <- nrow(M) D <- as.matrix(dist(scale(M))) # standardise: coordinates share a scale A0 <- cbind(1, M) cf <- matrix(NA_real_, n, ncol(M) + 1) pred <- rep(NA_real_, n) for (i in seq_len(n)) { use <- if (loo) seq_len(n) != i else rep(TRUE, n) d <- D[i, use] w <- exp(-theta * d / mean(d)) A <- A0[use, , drop = FALSE] * w b <- y[use] * w sv <- svd(A) # SVD solve: weights make A ill-conditioned tol <- max(dim(A)) * .Machine$double.eps * max(sv$d) co <- as.numeric(sv$v %*% (ifelse(sv$d > tol, 1 / sv$d, 0) * crossprod(sv$u, b))) cf[i, ] <- co pred[i] <- sum(c(1, M[i, ]) * co) } list(coef = cf, pred = pred, rho = cor(pred, y)) } smap <- function(x, E, tau = 1, tp = 1, theta = 0, loo = TRUE) { bl <- make_block(x, E, tau) ok <- (bl$t + tp) <= length(x) smap_block(bl$M[ok, , drop = FALSE], x[bl$t[ok] + tp], theta = theta, loo = loo) } ``` Two details earn their place. Distances are computed on standardised coordinates, so a variable measured in thousands does not swamp one measured in units. And the weighted system is solved by SVD rather than `solve`, because at high `theta` most weights are near zero and the design matrix becomes ill conditioned; the singular value cutoff keeps that from blowing up silently. ```{r} #| label: setup #| echo: false library(ggplot2) library(grid) te_ink <- "#16241d"; te_forest <- "#275139"; te_gold <- "#c9b458" te_rust <- "#b5534e"; te_paper <- "#f5f4ee"; te_line <- "#dad9ca" te_muted <- "#5d6b61"; te_green <- "#2f8f63" theme_te <- function(base_size = 11) { theme_minimal(base_size = base_size) + theme( plot.title = element_text(colour = te_ink, face = "bold", size = rel(1.0)), plot.subtitle = element_text(colour = te_muted, size = rel(0.85)), axis.title = element_text(colour = te_ink, size = rel(0.9)), axis.text = element_text(colour = te_muted, size = rel(0.8)), panel.grid.major = element_line(colour = te_line, linewidth = 0.3), panel.grid.minor = element_blank(), panel.background = element_rect(fill = te_paper, colour = NA), plot.background = element_rect(fill = te_paper, colour = NA), legend.position = "bottom", legend.title = element_text(colour = te_ink, size = rel(0.8)), legend.text = element_text(colour = te_muted, size = rel(0.75)), strip.text = element_text(colour = te_ink, size = rel(0.85)) ) } two_panel <- function(pA, pB) { grid.newpage() pushViewport(viewport(layout = grid.layout(1, 2))) print(pA, vp = viewport(layout.pos.row = 1, layout.pos.col = 1)) print(pB, vp = viewport(layout.pos.row = 1, layout.pos.col = 2)) } ``` ```{r} #| label: sci-tex-helper #| include: false sci_tex <- function(x, digits = 2) { s <- formatC(x, format = "e", digits = digits) paste0("$", sub("e.*", "", s), " \\times 10^{", as.integer(sub(".*e", "", s)), "}$") } ``` ## Check the claim before using it If `theta = 0` really is the global linear model, then the S-map at that setting, run on the full library with no leave-one-out, must reproduce `lm` exactly. Worth checking rather than believing. ```{r} #| label: identity gen_ricker <- function(n, burn = 500, r = 3.0, x0 = 0.6) { N <- n + burn; x <- numeric(N); x[1] <- x0 for (t in 1:(N - 1)) x[t + 1] <- x[t] * exp(r * (1 - x[t])) x[(burn + 1):N] } set.seed(4265) n <- 300 z <- gen_ricker(n) * exp(rnorm(n, 0, 0.05)) # 5% observation error s0 <- smap(z, E = 2, tp = 1, theta = 0, loo = FALSE) X <- embed(z, 3) m_lm <- lm(X[, 1] ~ X[, 2] + X[, 3]) coef_gap <- max(abs(sweep(s0$coef, 2, coef(m_lm)))) fit_gap <- max(abs(s0$pred - fitted(m_lm))) ``` The S-map's coefficients match `lm`'s, and its fitted values match `lm`'s fitted values, to within floating-point rounding. The S-map at `theta = 0` is not *like* a linear autoregressive model, it is one. ```{r} #| label: fp-guard-1 #| include: false stopifnot(coef_gap / max(abs(coef(m_lm))) < 1e-10, fit_gap / max(abs(fitted(m_lm))) < 1e-10) ``` ## The nonlinearity test Now turn the dial. The chaotic Ricker series from the previous post, and its linear counterpart: an autoregressive model fitted to it and simulated from, so the two share their autocorrelation structure and differ only in whether a deterministic map is underneath. ```{r} #| label: theta-sweep ar_fit <- ar(z, order.max = 10, method = "mle") set.seed(9265) z_lin <- as.numeric(arima.sim(n = n, model = list(ar = ar_fit$ar), sd = sqrt(ar_fit$var.pred))) thetas <- c(0, 0.1, 0.3, 0.5, 1, 2, 3, 4, 6, 8) rho_chaos <- sapply(thetas, function(th) smap(z, E = 2, tp = 1, theta = th)$rho) rho_lin <- sapply(thetas, function(th) smap(z_lin, E = 2, tp = 1, theta = th)$rho) gain_chaos <- max(rho_chaos) - rho_chaos[1] gain_lin <- max(rho_lin) - rho_lin[1] ``` ```{r} #| label: fig-theta #| fig-width: 8 #| fig-height: 4 #| fig-cap: "Forecast skill against the S-map localisation parameter. The chaotic series improves steeply as the fit becomes local; the linear series with matching autocorrelation gains nothing and then degrades." #| fig-alt: "Two line plots of forecast correlation against theta from zero to eight. The left panel for the chaotic series rises from about 0.80 to 0.99. The right panel for the linear series stays flat near 0.80 and declines slightly after theta of one." #| echo: false tA <- data.frame(theta = thetas, rho = rho_chaos) tB <- data.frame(theta = thetas, rho = rho_lin) pA <- ggplot(tA, aes(theta, rho)) + geom_line(colour = te_forest, linewidth = 0.8) + geom_point(colour = te_forest, size = 2) + geom_hline(yintercept = rho_chaos[1], colour = te_muted, linetype = 2, linewidth = 0.4) + labs(title = "Chaotic population", subtitle = "skill climbs as the fit localises", x = expression(theta), y = expression(rho)) + theme_te() pB <- ggplot(tB, aes(theta, rho)) + geom_line(colour = te_rust, linewidth = 0.8) + geom_point(colour = te_rust, size = 2) + geom_hline(yintercept = rho_lin[1], colour = te_muted, linetype = 2, linewidth = 0.4) + labs(title = "Linear series, same ACF", subtitle = "nothing to localise", x = expression(theta), y = expression(rho)) + theme_te() two_panel(pA, pB) ``` The chaotic series goes from `rho = `r sprintf("%.3f", rho_chaos[1])`` at `theta = 0` to `r sprintf("%.3f", max(rho_chaos))` at its best, a gain of `r sprintf("%+.3f", gain_chaos)`. The linear series starts at `r sprintf("%.3f", rho_lin[1])` and its best is `r sprintf("%.3f", max(rho_lin))`, a gain of `r sci_tex(gain_lin, 1)`, which is nothing next to the chaotic gain. Push `theta` further and the linear series gets actively worse (`r sprintf("%.3f", rho_lin[length(thetas)])` at `theta = 8`), because localising a fit that has no local structure to find just throws away data. This is the cleanest test in the toolkit. Same model, same data, one dial, and the answer is a comparison against a value the method itself produces at `theta = 0`. A dashed line marks that baseline in both panels. ## Interaction coefficients that move The dial test uses one species. The coefficients need more. If you have measured several species, use their contemporaneous values as the state instead of lagged values of one, and predict one species forward. The local coefficients are then estimates of $$\frac{\partial\, x_{i,t+1}}{\partial\, x_{j,t}}$$ at each point in time: how much species `j` currently affects species `i`. The test system is two competitors with asymmetric coupling, and we know its Jacobian exactly. Species X drives Y with strength 0.10; Y drives X with strength 0.02. Both are chaotic, and their abundances are almost uncorrelated. ```{r} #| label: coupled gen_coupled <- function(n, burn = 500, rx = 3.8, ry = 3.5, bxy = 0.02, byx = 0.1, x0 = 0.4, y0 = 0.2) { N <- n + burn; x <- numeric(N); y <- numeric(N); x[1] <- x0; y[1] <- y0 for (t in 1:(N - 1)) { x[t + 1] <- x[t] * (rx - rx * x[t] - bxy * y[t]) y[t + 1] <- y[t] * (ry - ry * y[t] - byx * x[t]) } data.frame(x = x[(burn + 1):N], y = y[(burn + 1):N]) } set.seed(4265) nc <- 400 d <- gen_coupled(nc) M <- cbind(x = d$x[-nc], y = d$y[-nc]) ## true Jacobian entries, from differentiating the model by hand true_YX <- -0.10 * M[, "y"] # d y[t+1] / d x[t] : strong link true_XY <- -0.02 * M[, "x"] # d x[t+1] / d y[t] : weak link true_XX <- 3.8 - 2 * 3.8 * M[, "x"] - 0.02 * M[, "y"] th_grid <- c(0, 0.5, 1, 2, 3, 4, 6, 8) rho_X <- sapply(th_grid, function(th) smap_block(M, d$x[-1], theta = th)$rho) rho_Y <- sapply(th_grid, function(th) smap_block(M, d$y[-1], theta = th)$rho) th_best <- th_grid[which.max(rho_Y)] fitX <- smap_block(M, d$x[-1], theta = th_best) # x[t+1] ~ (x[t], y[t]) fitY <- smap_block(M, d$y[-1], theta = th_best) # y[t+1] ~ (x[t], y[t]) est_YX <- fitY$coef[, 2] # coefficient on x in the y equation : strong link est_XY <- fitX$coef[, 3] # coefficient on y in the x equation : weak link est_XX <- fitX$coef[, 2] # self effect : strong cor_strong <- cor(est_YX, true_YX) cor_weak <- cor(est_XY, true_XY) cor_self <- cor(est_XX, true_XX) ## the same model with Y driving X as strongly as X drives Y, theta chosen the same way d_eq <- gen_coupled(nc, bxy = 0.1) M_eq <- cbind(x = d_eq$x[-nc], y = d_eq$y[-nc]) th_eq <- th_grid[which.max(sapply(th_grid, function(th) smap_block(M_eq, d_eq$y[-1], theta = th)$rho))] cor_eq <- cor(smap_block(M_eq, d_eq$x[-1], theta = th_eq)$coef[, 3], -0.1 * M_eq[, "x"]) # Y on X, now strong cor_eq_y <- cor(smap_block(M_eq, d_eq$y[-1], theta = th_eq)$coef[, 2], -0.1 * M_eq[, "y"]) # X on Y, as before ``` Cross validation picks `theta = `r th_best``, the top of the grid searched, and rho is still inching up at that edge; the forecasts are near perfect either way (`rho = `r sprintf("%.5f", max(rho_Y))`` for the Y equation). Now compare each estimated coefficient against the truth we can compute analytically. ```{r} #| label: fig-coefs #| fig-width: 8 #| fig-height: 4 #| fig-cap: "Estimated S-map interaction coefficients against the true Jacobian entries, for the strong link and the weak link in the same system at the same settings." #| fig-alt: "Two scatter plots of estimated coefficient against true coefficient. The left panel for the strong link shows two clusters along the one-to-one line: a tight one at the most negative true values and a looser cloud, spread well above and below the line, at the less negative ones. The right panel for the weak link shows a diffuse vertical scatter with no relationship to the truth." #| echo: false cA <- data.frame(truth = true_YX, est = est_YX) cB <- data.frame(truth = true_XY, est = est_XY) pA <- ggplot(cA, aes(truth, est)) + geom_abline(slope = 1, intercept = 0, colour = te_gold, linewidth = 0.8) + geom_point(colour = te_forest, alpha = 0.5, size = 1.1) + labs(title = "Strong link: X on Y", subtitle = paste0("r = ", sprintf("%.3f", cor_strong)), x = "true coefficient", y = "S-map estimate") + theme_te() pB <- ggplot(cB, aes(truth, est)) + geom_abline(slope = 1, intercept = 0, colour = te_gold, linewidth = 0.8) + geom_point(colour = te_rust, alpha = 0.5, size = 1.1) + labs(title = "Weak link: Y on X", subtitle = paste0("r = ", sprintf("%.3f", cor_weak)), x = "true coefficient", y = "S-map estimate") + theme_te() two_panel(pA, pB) ``` The strong link is recovered: the estimated coefficients track the true ones at a correlation of `r sprintf("%.3f", cor_strong)`, and their mean, `r sprintf("%.4f", mean(est_YX))`, sits on the true mean of `r sprintf("%.4f", mean(true_YX))`. The self effect comes back too, at `r sprintf("%.3f", cor_self)`. Given that the two species' abundances correlate at only `r sprintf("%.2f", cor(d$x, d$y))`, recovering this time variation from the raw series is a real result, although most of it is Y alternating between high and low values from one step to the next; within each of the two clusters in the left panel the estimates track the variation in the truth only weakly. The weak link is not recovered. The estimates correlate with the truth at `r sprintf("%.3f", cor_weak)`: no relationship at all. ## The part to be careful about Look at the right panel again. Those estimates are not flat, not obviously broken, and not implausible. Their mean is `r sprintf("%.4f", mean(est_XY))` against a true mean of `r sprintf("%.4f", mean(true_XY))`, which is close enough to pass a sanity check. Plot them against time and you get a wiggly line that looks exactly like a fluctuating interaction strength. Every diagnostic available *within the analysis* is green: the forecasts are near perfect, `theta` was chosen by cross validation, the self terms validate closely against the truth in the same fit, and the strong cross term follows the truth from one of its two clusters to the other. And the time course is noise. The obvious suspect is size. The weak coefficient is about `r sprintf("%.0f", 0.10/0.02)` times smaller than the strong one and roughly two orders of magnitude smaller than the self term, which ranges over several units, and a local linear fit has to separate that contribution from everything else happening in a neighbourhood of the attractor. Size is not the whole story, though. The code block that defines `gen_coupled()` ends by rerunning the same model with Y driving X as strongly as X drives Y, with `theta` again chosen by cross validation: the coefficient on Y in the X equation is still not recovered, correlating with its truth at `r sprintf("%.3f", cor_eq)`, while the one on X in the Y equation comes back at `r sprintf("%.3f", cor_eq_y)`. In this model a coefficient as large as the recovered one can fail too, so size alone does not decide which ones a local linear fit can recover. Nothing warns you, because the quantity that fails is not the one the diagnostics measure: prediction is dominated by the large terms, so the model predicts brilliantly while the small coefficient it reports is arbitrary. This is an honest limit of S-map interaction strengths, shown here on one simulated system, and it matters in community ecology, where most pairwise effects are weak and a few are strong. The practical reading: - The **average** of a coefficient can be informative even when its **fluctuations** are not. Here the weak link's mean was roughly right while its time course was noise. Reporting "the interaction varied through time from A to B" is a much stronger claim than reporting its average, and it needs much more support. - Size alone does not tell you which coefficients to trust: in the rerun above a link as strong as the recovered one failed in the other equation, so a large coefficient needs the same check as a small one. - Validate on something. A simulation with a known Jacobian, a perturbation experiment, an independent estimate of the interaction. If nothing anchors the coefficient, its wiggles are decoration. - At `theta = 0` the coefficients are useless even for the strong link, because a single global line through a curved map has no local derivative to report. Localisation is what makes the coefficients mean anything, and it is also what makes them variable. There is a further wrinkle worth knowing about: with several species the local design matrix is often close to collinear, and the coefficient split between correlated species is then poorly determined even when their sum is not. Regularised versions of the S-map exist (Cenci and colleagues), introduced because the plain version is sensitive to process noise; a ridge-type penalty on the local coefficients also steadies their split when the design is near collinear. ## Where to go next Two of these coefficients said the same thing about the true structure, in opposite directions: X affects Y appreciably, Y affects X hardly at all. The S-map got the strong one right and the weak one wrong, but it needed *both species measured* to say anything at all. The next question is what to do when you cannot measure both, or when you can but want to know which one is driving which. [Convergent cross mapping](../convergent-cross-mapping/) attacks that directly, and it needs only the reconstruction from the previous post. It also inherits every fragility on display here, plus some of its own, which is what the [checking post](../checking-an-empirical-dynamic-model/) is for. ## References Ushio, Hsieh, Masuda, Deyle, Ye, Chang, Sugihara & Kondoh 2018 Nature 554(7692):360-363 (10.1038/nature25504) Sugihara 1994 Philosophical Transactions of the Royal Society A 348(1688):477-495 (10.1098/rsta.1994.0106) Deyle, May, Munch & Sugihara 2016 Proceedings of the Royal Society B 283(1822):20152258 (10.1098/rspb.2015.2258) Cenci, Sugihara & Saavedra 2019 Methods in Ecology and Evolution 10(5):650-660 (10.1111/2041-210X.13150) Hsieh, Glaser, Lucas & Sugihara 2005 Nature 435(7040):336-340 (10.1038/nature03553) Sugihara, May, Ye, Hsieh, Deyle, Fogarty & Munch 2012 Science 338(6106):496-500 (10.1126/science.1227079) ## Related tutorials - [Simplex projection and delay embedding](../simplex-projection-and-embedding/) - [Convergent cross mapping in ecology](../convergent-cross-mapping/) - [Checking an empirical dynamic model](../checking-an-empirical-dynamic-model/) - [Checking a penalised regression](../checking-a-penalised-regression/)