Interaction strengths from community time series

R
MAR models
time series
species interactions
ecology tutorial
What a MAR(1) coefficient means, how much a sixty-point series can say about it, and why shared weather blurs the estimated effects of synchronised species.
Author

Tidy Ecology

Published

2026-06-11

Modified

2026-10-14

Corrected 14 October 2026: the post compared the diagonal bias with a signed average of the off-diagonal biases, in which cells of opposite sign cancel, and said the direct interactions do not look stronger than they are; it now compares mean absolute biases and says most direct interactions look slightly stronger; and a section heading no longer says weak links are worse than absent.

Corrected 7 October 2026: the post said the small-sample bias lands on the diagonal whatever the parameters, and now says it does so here but partly moves into an interaction in the shared environment; it said one run in five would call the weak competitor a mutualist, and now adds that almost none of those wrong-signed estimates pass the screen; and the description no longer says a shared environment ruins the off-diagonal estimates.

Corrected 30 September 2026: the post said a shared environment breaks identifiability, that the data could no longer separate the interaction from the shared response, and called this fatal; it now says the coefficients stay identified but are estimated much less precisely.

Fit a MAR(1) model and you get sixteen numbers for a four-species community. Fitting a MAR(1) model to community time series showed the fitting is three lines of base R and that the estimator is consistent. Consistent is a statement about infinite data. This post asks what sixty points buy you.

The answer has three parts, and only the first is the one people expect.

Setup: packages, colours and plot theme 
library(ggplot2)
library(grid)

te_ink    <- "#16241d"
te_forest <- "#275139"
te_gold   <- "#c9b458"
te_rust   <- "#b5534e"
te_sage   <- "#93a87f"
te_line   <- "#dad9ca"
te_muted  <- "#5d6b61"

theme_te <- function(base_size = 12) {
  theme_minimal(base_size = base_size) +
    theme(
      text             = element_text(colour = te_ink),
      plot.title       = element_text(face = "bold", size = rel(1.02), colour = te_ink),
      plot.subtitle    = element_text(colour = te_muted, size = rel(0.88)),
      axis.title       = element_text(colour = te_ink, size = rel(0.92)),
      axis.text        = element_text(colour = te_muted, size = rel(0.82)),
      panel.grid.major = element_line(colour = te_line, linewidth = 0.3),
      panel.grid.minor = element_blank(),
      panel.background = element_rect(fill = "#f5f4ee", colour = NA),
      plot.background  = element_rect(fill = "#f5f4ee", colour = NA),
      legend.position  = "bottom",
      legend.title     = element_text(size = rel(0.82)),
      legend.text      = element_text(size = rel(0.78)),
      legend.key.size  = unit(0.8, "lines")
    )
}

two_panel <- function(p1, p2) {
  grid.newpage()
  pushViewport(viewport(layout = grid.layout(1, 2)))
  print(p1, vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
  print(p2, vp = viewport(layout.pos.row = 1, layout.pos.col = 2))
}

Btrue <- matrix(c( 0.50, -0.10, -0.35,  0.00,
                  -0.15,  0.60,  0.00,  0.00,
                   0.30,  0.00,  0.45, -0.30,
                   0.00,  0.00,  0.25,  0.55), 4, 4, byrow = TRUE)
Atrue <- c(1.20, 0.90, 0.40, 0.30)
Strue <- matrix(c(0.09, 0.03, 0, 0, 0.03, 0.09, 0, 0,
                  0, 0, 0.04, 0.01, 0, 0, 0.01, 0.04), 4, 4, byrow = TRUE)
Scol  <- Strue; Scol[1, 2] <- Scol[2, 1] <- 0.087
p <- 4

sim_mar1 <- function(n, A, B, S, seed, burn = 200) {
  set.seed(seed); p <- length(A); L <- chol(S)
  x <- solve(diag(p) - B, A); X <- matrix(0, n, p)
  for (t in seq_len(burn + n)) {
    x <- A + B %*% x + as.vector(rnorm(p) %*% L)
    if (t > burn) X[t - burn, ] <- x
  }
  X
}
mar1_cls <- function(X) {
  n <- nrow(X); Y <- X[-1, , drop = FALSE]; D <- cbind(1, X[-n, , drop = FALSE])
  cf <- solve(crossprod(D), crossprod(D, Y)); E <- Y - D %*% cf
  list(A = cf[1, ], B = t(cf[-1, , drop = FALSE]), Sigma = crossprod(E) / nrow(Y), E = E)
}
mar1_se <- function(X) {
  n <- nrow(X); pp <- ncol(X); Y <- X[-1, , drop = FALSE]; D <- cbind(1, X[-n, , drop = FALSE])
  XtXi <- solve(crossprod(D)); q <- nrow(Y); k <- ncol(D); SE <- matrix(0, pp, pp)
  for (i in 1:pp) {
    r <- Y[, i] - D %*% (XtXi %*% crossprod(D, Y[, i]))
    SE[i, ] <- sqrt(sum(r^2) / (q - k) * diag(XtXi))[-1]
  }
  SE
}

What the coefficients mean

\(b_{ij}\) for \(i \neq j\) is the effect of species \(j\) on species \(i\)’s per-capita growth rate, per unit of log abundance. Negative means suppression, positive means facilitation, zero means no direct link. The sign convention is the ecologist’s, so a competitor sits at a negative off-diagonal.

The diagonal is a different animal. \(b_{ii}\) measures density dependence, and the useful way to read it is as persistence:

  • \(b_{ii} = 1\): no density dependence at all. A shock never decays and the species does a random walk.
  • \(b_{ii} = 0.9\): strong persistence, slow return, long memory of the last bad winter.
  • \(b_{ii} = 0\): total density dependence. Whatever happened last time step is gone; the species snaps straight back to its mean.
  • \(b_{ii} < 0\): overcompensation, the boom-and-bust of an overshooting population.

So the diagonal and the off-diagonal are not two flavours of the same number. One is a memory length, the other is a direct effect.

The bias is asymmetric

Simulate the same community 600 times at \(n = 60\) and average the estimates.

mc <- function(n, S, M, seed0) {
  Bh <- array(0, c(M, p, p)); Sg <- array(0, c(M, p, p))
  for (m in 1:M) {
    Xm <- sim_mar1(n, Atrue, Btrue, S, seed = seed0 + m)
    f <- mar1_cls(Xm)
    Bh[m, , ] <- f$B
    Sg[m, , ] <- abs(f$B / mar1_se(Xm)) > 2   # crude 2-SE screen
  }
  list(B = Bh, sig = Sg)
}

r60  <- mc(60, Strue, 600, 60000)
bias <- apply(r60$B, c(2, 3), mean) - Btrue
round(bias, 3)
       [,1]   [,2]   [,3]   [,4]
[1,] -0.057 -0.006 -0.022 -0.018
[2,]  0.002 -0.079 -0.011  0.009
[3,]  0.009  0.000 -0.036 -0.014
[4,] -0.007  0.005  0.013 -0.054
c(diagonal     = mean(diag(bias)),
  off_diagonal = mean(bias[row(bias) != col(bias)]))
    diagonal off_diagonal 
 -0.05642935  -0.00342599

The diagonal is biased by -0.0564. The off-diagonal cells average -0.0034, but that average lets cells of opposite sign cancel: their mean absolute bias is 0.0097, about 6 times smaller than the diagonal’s, and in these replicates the average estimate of 5 of the 6 true links sits slightly further from zero than the truth.

Read that in ecological units. The estimated diagonal is systematically too low, and a low \(b_{ii}\) means fast return to the mean. Density dependence looks stronger than it is, and most direct interactions look only slightly stronger than they are. This is the small-sample bias of an autoregressive coefficient, known since Marriott and Pope (1954), and here it lands mainly on the diagonal, where the lagged self-terms sit. The asymmetry is a property of this parameter set rather than a law: in the shared environment below, part of the bias moves off the diagonal.

It is an \(O(1/n)\) effect, so it disappears with length. Watch how slowly.

sweep <- do.call(rbind, lapply(c(25, 40, 60, 120, 400), function(n) {
  r <- mc(n, Strue, 250, n * 1000)
  s <- apply(r$sig, c(2, 3), mean)
  z0 <- which(Btrue == 0)
  data.frame(n = n,
             bias  = mean(diag(apply(r$B, c(2, 3), mean) - Btrue)),
             mc_sd = mean(apply(r$B, c(2, 3), sd)),
             false_pos = mean(s[z0]),
             detection = mean(s[-z0]))
}))
round(sweep, 4)
    n    bias  mc_sd false_pos detection
1  25 -0.1344 0.2173    0.0793    0.4116
2  40 -0.0808 0.1578    0.0567    0.5940
3  60 -0.0606 0.1218    0.0613    0.7356
4 120 -0.0252 0.0804    0.0547    0.8524
5 400 -0.0090 0.0416    0.0500    0.9644
cells <- data.frame(
  err  = as.vector(bias),
  kind = ifelse(as.vector(diag(4)) == 1, "diagonal", "off-diagonal"),
  cell = paste0("b", as.vector(row(bias)), as.vector(col(bias)))
)

p1 <- ggplot(cells, aes(err, reorder(cell, err), colour = kind)) +
  geom_vline(xintercept = 0, colour = te_ink, linewidth = 0.4) +
  geom_point(size = 2.2) +
  scale_colour_manual(values = c(diagonal = te_rust, `off-diagonal` = te_sage)) +
  labs(title = "Where the bias lands", subtitle = "600 replicates, n = 60",
       x = "Mean estimate minus truth", y = NULL, colour = NULL) +
  theme_te(11)

p2 <- ggplot(sweep, aes(n, bias)) +
  geom_hline(yintercept = 0, colour = te_ink, linewidth = 0.4) +
  geom_line(colour = te_forest, linewidth = 0.7) +
  geom_point(colour = te_forest, size = 2.2) +
  geom_vline(xintercept = 60, colour = te_gold, linetype = 2, linewidth = 0.5) +
  annotate("text", x = 72, y = min(sweep$bias) / 2, label = "a real series",
           hjust = 0, size = 3, colour = te_muted) +
  scale_x_log10(breaks = sweep$n, labels = sweep$n) +
  labs(title = "It goes away like 1/n", subtitle = "Mean bias of the diagonal",
       x = "Series length", y = "Diagonal bias") +
  theme_te(11)

two_panel(p1, p2)
Two panels. Left, a dot plot of estimation error per matrix cell showing all four diagonal cells clearly below zero and twelve off-diagonal cells scattered around zero. Right, a line showing diagonal bias rising towards zero as series length increases from twenty-five to four hundred.
Figure 1: Left: mean estimation error for each of the sixteen cells at n = 60, split by position. Every diagonal cell sits below zero; the off-diagonal cells straddle it. Right: the diagonal bias shrinks with series length, but sixty points still leaves close to half of it.

At 25 points the diagonal is off by -0.1344, which is a quarter of a typical coefficient. At 400 it is -0.0090 and gone. Sixty points, a respectable monitoring series, still carries -0.0606.

The screen is calibrated, and that is not enough

This community has six true zeros. Screen every cell at two standard errors and count how often a zero cell fires.

sg <- apply(r60$sig, c(2, 3), mean)
z0 <- which(Btrue == 0)
c(false_positive_on_6_true_zeros = mean(sg[z0]),
  detection_on_10_true_effects   = mean(sg[-z0]))
false_positive_on_6_true_zeros   detection_on_10_true_effects 
                    0.05888889                     0.73450000

0.059 on the true zeros. The screen is calibrated: it is doing about what a five per cent test should do, and nothing is broken.

Now multiply. Six zero cells at 0.059 is 0.35 false interactions per fitted model, on average. A ten-species community has ninety off-diagonal cells; if all of them were true zeros, a single honest, correctly calibrated analysis would still hand you around 5 interactions that do not exist. They will not be flagged. They will look exactly like the real ones, and the web you draw will carry them as arrows.

This is the multiple testing problem wearing an ecological hat, and the fix is the usual one: control the false discovery rate across the matrix rather than testing each cell as if it were the only cell. The false discovery rate and Benjamini-Hochberg has the mechanics.

A shared environment inflates the variance

Now the third part, and the one that reaches furthest. Keep \(B\) identical and correlate the environmental shocks of species 1 and 2 at 0.97, a deliberately extreme version of a shared weather signal, chosen so the mechanism shows through 300 replicates.

Xb <- sim_mar1(60, Atrue, Btrue, Strue, seed = 4269)
Xc <- sim_mar1(60, Atrue, Btrue, Scol,  seed = 4269)

vf <- function(M) {
  Ml <- M[-nrow(M), ]
  sapply(1:p, function(j) 1 / (1 - summary(lm(Ml[, j] ~ Ml[, -j]))$r.squared))
}
rbind(baseline = vf(Xb), shared = vf(Xc))
             [,1]     [,2]     [,3]     [,4]
baseline 1.078509 1.063237 1.151190 1.132686
shared   4.126688 4.387299 1.130292 1.366421

The predictors in a MAR(1) design are the species, so environmental synchrony is collinearity in the design matrix. The variance inflation factor for species 1 goes from 1.08 to 4.13 and for species 2 from 1.06 to 4.39. Species 3 and 4, whose environments were left alone, barely move.

The cost lands on the coefficients that involve the synchronised pair.

sb <- apply(Bb, c(2, 3), sd); sc <- apply(Bc, c(2, 3), sd)
tab <- t(sapply(list(c(1, 1), c(1, 2), c(2, 1)), function(ij) {
  c(baseline = sb[ij[1], ij[2]], shared = sc[ij[1], ij[2]],
    ratio = sc[ij[1], ij[2]] / sb[ij[1], ij[2]])
}))
rownames(tab) <- c("b[1,1]", "b[1,2]", "b[2,1]")
round(tab, 3)
       baseline shared ratio
b[1,1]    0.120  0.261 2.170
b[1,2]    0.125  0.286 2.289
b[2,1]    0.130  0.258 1.988

The sampling standard deviation of \(b_{11}\) goes from 0.120 to 0.261, a factor of 2.17. For \(b_{12}\) it is 2.29. Nothing about the ecology changed. Two species started responding to the same weather, and the data became much less able to tell their interaction apart from their shared response: the coefficients are still identified, but their sampling spread more than doubled. The small-sample bias moves as well: across the same shared-environment replicates the mean error of \(b_{11}\) is +0.013 and that of the weak competitor \(b_{12}\) is -0.093, so the weak competitor now looks stronger than it is while the self-term of species 1 is close to unbiased.

Which is the ordinary lesson of collinearity, arriving through a door ecologists do not usually watch. Collinearity and VIF is the static version of exactly this; here the correlated predictors are the community itself.

The honest limit

Three separate things degrade a MAR(1) interaction estimate, and they do not degrade it in the same way.

The small-sample bias is systematic, sits mainly on the diagonal under the baseline shocks, and shrinks predictably with length. You can reason about it. The screening error rate is calibrated per cell and catastrophic per matrix, and the fix is a multiple-comparisons fix. Neither of those is fatal.

The third comes closer. Environmental synchrony inflates the variance of exactly the coefficients you most want, and no amount of care in the fitting recovers what the design does not contain. A long series helps; a synchronised community needs something other than the abundance series to break the tie.

So the reportable output of a sixty-point MAR(1) fit is a sign structure with a false-discovery control on it, not a table of interaction strengths to two decimal places. Community stability from a MAR(1) model is the next question: what happens when you feed these coefficients into a stability metric.

References

Newsletter

Get updates by email

An occasional email when tutorials are added or substantially corrected. No spam; unsubscribe anytime.

By subscribing you agree to receive these emails and confirm your address once. See the privacy policy.