Checking a MAR(1) model

R
MAR models
model diagnostics
time series
ecology tutorial
Observation error drags every MAR(1) stability metric towards calm while the residual diagnostics stay silent. Here is the size of it, and why the residuals cannot help.
Author

Tidy Ecology

Published

2026-08-23

The three previous posts built a MAR(1) analysis: the fit, the interaction coefficients, the stability metrics. Every number in that chain comes out of one matrix, \(B\), and inherits whatever is wrong with it.

This post is about the largest thing that is wrong with it. It has three properties that make it worse than an ordinary nuisance: it is big, it always pushes the same direction, and the standard diagnostics do not see it.

The problem: you never count the community, you sample it

MAR(1) treats \(X_t\) as the community’s log abundance. What you have is a net haul, a quadrat, a transect. The count carries sampling error on top of the real state, and that error is not part of the model as written. Write \(\tau^2\) for its variance.

Add it to the simulation and refit, three hundred times per noise level.

tau_grid <- c(0, 0.02, 0.05, 0.10, 0.20)

att <- do.call(rbind, lapply(tau_grid, function(tau2) {
  Bh <- array(0, c(300, p, p))
  for (m in 1:300) {
    Xm <- sim_mar1(60, Atrue, Btrue, Strue, seed = 7000 + m)
    set.seed(107000 + m)
    Bh[m, , ] <- mar1_cls(Xm + matrix(rnorm(60 * p, 0, sqrt(tau2)), 60, p))$B
  }
  Bm <- apply(Bh, c(2, 3), mean)
  data.frame(tau2 = tau2,
             mean_diag = mean(diag(Bm)),
             pct_of_truth = 100 * mean(diag(Bm)) / mean(diag(Btrue)),
             return_rate = max(Mod(eigen(Bm)$values)),
             det_index = abs(det(Bm))^(2 / p))
}))
round(att, 4)
  tau2 mean_diag pct_of_truth return_rate det_index
1 0.00    0.4678      89.0993      0.6118    0.2974
2 0.02    0.3753      71.4797      0.4907    0.1919
3 0.05    0.2911      55.4414      0.3923    0.1146
4 0.10    0.2123      40.4362      0.3015    0.0599
5 0.20    0.1362      25.9420      0.2069    0.0240

Read the table down, not across.

At \(\tau^2 = 0\), with perfect observation, the mean diagonal is already 89 per cent of the truth. That is the small-sample bias from the previous post, and on its own it is a mild annoyance.

At \(\tau^2 = 0.05\), which is a modest amount of counting error next to a process variance of 0.09, the diagonal is down to 55 per cent. At \(\tau^2 = 0.20\) it is 26 per cent. The interaction matrix has lost three quarters of itself.

Now look at what that does to the stability metrics, which is the part that matters. The true return rate is 0.6555. The estimates run 0.6118, 0.4907, 0.3923, 0.3015, 0.2069 as the observation error grows.

The community looks more stable the worse your data is. Sloppier counting buys you a faster return rate, a lower determinant index, and a more comfortable story about resilience.

That is not a random error you can average away across studies. It is a systematic gradient, and it runs in the direction of the answer people like.

p1 <- ggplot(att, aes(tau2, mean_diag)) +
  geom_hline(yintercept = mean(diag(Btrue)), colour = te_ink, linetype = 2, linewidth = 0.5) +
  annotate("text", x = 0.11, y = mean(diag(Btrue)) + 0.018, label = "true mean b_ii",
           colour = te_ink, size = 3) +
  geom_line(colour = te_rust, linewidth = 0.8) +
  geom_point(colour = te_rust, size = 2.2) +
  ylim(0, 0.58) +
  labs(title = "The interaction matrix shrinks",
       subtitle = "300 refits per level, n = 60",
       x = "Observation error variance", y = "Mean estimated diagonal") +
  theme_te(11)

met <- rbind(data.frame(tau2 = att$tau2, v = att$return_rate, w = "return rate",
                        truth = max(Mod(eigen(Btrue)$values))),
             data.frame(tau2 = att$tau2, v = att$det_index, w = "|det B|^(2/p)",
                        truth = abs(det(Btrue))^(2 / p)))

p2 <- ggplot(met, aes(tau2, v, colour = w)) +
  geom_hline(aes(yintercept = truth, colour = w), linetype = 2, linewidth = 0.5) +
  geom_line(linewidth = 0.8) +
  geom_point(size = 2.2) +
  scale_colour_manual(values = c(`return rate` = te_forest, `|det B|^(2/p)` = te_gold)) +
  labs(title = "So does every metric built on it",
       subtitle = "Dashes are the true values",
       x = "Observation error variance", y = "Metric", colour = NULL) +
  theme_te(11)

two_panel(p1, p2)
Two panels. Left, a line falling steeply from just under the true mean diagonal towards zero as observation error grows. Right, two lines for return rate and the determinant index both falling steeply away from their true values.
Figure 1: Left: the mean estimated diagonal against observation error variance, with the truth marked. Right: the two headline stability metrics falling with the same noise. Every point is 300 refits of the identical community.

The diagnostics are silent

The natural defence is the residual check. A MAR(1) fit has residuals; if the model is wrong, the residuals should say so. Ljung-Box at lag 5 on each species, four hundred replicates, five per cent level.

lb <- do.call(rbind, lapply(c(0, 0.05, 0.20), function(tau2) {
  rej <- matrix(FALSE, 400, p)
  for (m in 1:400) {
    Xt <- sim_mar1(60, Atrue, Btrue, Strue, seed = 300000 + m)
    set.seed(600000 + m)
    E <- mar1_cls(Xt + matrix(rnorm(60 * p, 0, sqrt(tau2)), 60, p))$E
    rej[m, ] <- sapply(1:p, function(i) Box.test(E[, i], lag = 5, type = "Ljung-Box")$p.value < 0.05)
  }
  data.frame(tau2 = tau2, sp1 = mean(rej[, 1]), sp2 = mean(rej[, 2]),
             sp3 = mean(rej[, 3]), sp4 = mean(rej[, 4]), any = mean(apply(rej, 1, any)))
}))
round(lb, 3)
  tau2   sp1   sp2   sp3   sp4   any
1 0.00 0.030 0.025 0.022 0.030 0.108
2 0.05 0.022 0.032 0.048 0.032 0.130
3 0.20 0.030 0.030 0.022 0.020 0.098

At \(\tau^2 = 0\) the test fires on 0.107 of replicates, which is roughly what four tests at five per cent should do. At \(\tau^2 = 0.20\), where the interaction matrix has collapsed by 74 per cent, the test fires on 0.098. If anything it fires slightly less.

The residuals are white. Perfectly, honestly white, at every noise level. There is no diagnostic signal at all, because there is nothing wrong with the residuals: attenuation does not leave serial correlation behind, it just rescales the coefficient. A model can be 74 per cent wrong about the quantity you care about and pass its residual check every time.

lbl <- do.call(rbind, lapply(1:nrow(lb), function(i) {
  data.frame(tau2 = factor(sprintf("%.2f", lb$tau2[i])),
             species = paste("Species", 1:4),
             rate = as.numeric(lb[i, c("sp1", "sp2", "sp3", "sp4")]))
}))

p1 <- ggplot(lbl, aes(tau2, rate, fill = species)) +
  geom_col(position = position_dodge(0.8), width = 0.7) +
  geom_hline(yintercept = 0.05, colour = te_rust, linetype = 2, linewidth = 0.7) +
  annotate("text", x = 0.7, y = 0.062, label = "nominal", colour = te_rust, size = 3) +
  scale_fill_manual(values = c(te_forest, te_sage, te_gold, te_muted)) +
  ylim(0, 0.1) +
  labs(title = "Ljung-Box, lag 5", subtitle = "400 replicates per level",
       x = "Observation error variance", y = "Rejection rate", fill = NULL) +
  theme_te(11)

cmp <- rbind(
  data.frame(tau2 = lb$tau2, v = lb$any, w = "diagnostic fires (any species)"),
  data.frame(tau2 = att$tau2[c(1, 3, 5)], v = att$pct_of_truth[c(1, 3, 5)] / 100,
             w = "interaction strength retained")
)

p2 <- ggplot(cmp, aes(tau2, v, colour = w)) +
  geom_line(linewidth = 0.8) +
  geom_point(size = 2.4) +
  scale_colour_manual(values = c(`diagnostic fires (any species)` = te_muted,
                                 `interaction strength retained` = te_rust)) +
  ylim(0, 1) +
  labs(title = "Flat where it hurts", subtitle = "Both on the same noise axis",
       x = "Observation error variance", y = "Proportion", colour = NULL) +
  theme_te(11) +
  guides(colour = guide_legend(ncol = 1))

two_panel(p1, p2)
Two panels. Left, bars of Ljung-Box rejection rate at three noise levels, all close to the nominal five per cent line. Right, two lines against noise level: a flat one for the diagnostic and a steeply falling one for the retained interaction strength.
Figure 2: Left: Ljung-Box rejection rates sit at the nominal level regardless of observation error. Right: the same information as a contrast, with the collapse of the interaction matrix on the same axis. The diagnostic is flat where the damage is steepest.

The mechanism is exact, and it is an old one

None of this is mysterious. Observation error in a MAR(1) model is measurement error in a predictor, because the lagged state is the predictor. That is regression dilution, and it has a closed form: the estimate converges to the truth times the reliability ratio.

\[\hat b \;\to\; b \cdot \frac{V}{V + \tau^2}\]

Check it at \(p = 1\), where the algebra is clean, with long series so that small-sample bias is out of the way.

b <- 0.7; s2 <- 0.25
V1 <- s2 / (1 - b^2)

rel <- do.call(rbind, lapply(c(0, 0.05, 0.25, 1.00), function(tau2) {
  bh <- numeric(400)
  for (m in 1:400) {
    set.seed(50000 + m)
    x <- numeric(3000)
    for (t in 2:3000) x[t] <- b * x[t - 1] + rnorm(1, 0, sqrt(s2))
    set.seed(150000 + m)
    xo <- x + rnorm(3000, 0, sqrt(tau2))
    bh[m] <- coef(lm(xo[-1] ~ xo[-3000]))[2]
  }
  data.frame(tau2 = tau2, empirical = mean(bh),
             predicted = b * V1 / (V1 + tau2),
             gap = mean(bh) - b * V1 / (V1 + tau2),
             reliability = V1 / (V1 + tau2))
}))
round(rel, 4)
  tau2 empirical predicted    gap reliability
1 0.00    0.6996    0.7000 -4e-04      1.0000
2 0.05    0.6346    0.6352 -6e-04      0.9074
3 0.25    0.4628    0.4636 -7e-04      0.6623
4 1.00    0.2296    0.2303 -7e-04      0.3289

The stationary variance is 0.4902, and the formula predicts every row to within 0.0007. This is not an approximation that happens to work. It is the same reliability ratio as Measurement error and regression dilution, arriving in a time series because the lag makes the response into next step’s predictor.

ggplot(rel, aes(predicted, empirical)) +
  geom_abline(slope = 1, intercept = 0, colour = te_line, linewidth = 0.8) +
  geom_point(aes(colour = factor(tau2)), size = 3.4) +
  scale_colour_manual(values = c(te_ink, te_forest, te_gold, te_rust)) +
  labs(title = "b times V/(V + tau squared)",
       subtitle = "Predicted against empirical, 400 replicates per level",
       x = "Predicted by the dilution formula", y = "Empirical mean estimate",
       colour = "Observation error variance") +
  theme_te()
A scatter of empirical against predicted autoregressive coefficient for four noise levels, all four points lying on the one-to-one line running from about 0.23 up to 0.70.
Figure 3: The dilution formula against 400 simulations per noise level, at p = 1 with 3000-step series. The points sit on the line; the largest gap across the whole range is under 0.001.

So the attenuation is not a bug to be diagnosed. It is the correct behaviour of least squares given a noisy predictor, and asking the residuals about it is asking the wrong question. The residuals describe the fit; the fit is fine. It is the target that moved.

The remedy has to come from outside

The reliability ratio needs \(\tau^2\), and \(\tau^2\) is not in the data. One series of counts cannot separate the wobble that is the community from the wobble that is your net. Something external has to supply it:

  • Repeated sampling within a time step. Two hauls on the same day differ only by observation error, which estimates \(\tau^2\) directly. This is the cleanest answer and it is a field-design decision, not an analysis decision.
  • A state-space model. Put the latent true state and the observation layer in explicitly and let the Kalman filter separate them. The Gompertz state-space model is the one-species version of exactly this, and MARSS (Holmes, Ward and Wills 2012) does the multivariate case.

The second option is the popular one, and it comes with its own honest limit that is worth knowing before you reach for it. The two variances are only weakly identified from a single series: Knape (2008) showed that for realistic parameter values the density dependence and the two error variances trade off against each other in a way one series cannot resolve, and Auger-Methe et al. (2016) found that even simple linear Gaussian state-space models hit parameter and state estimation problems more often than their popularity suggests. Dennis et al. (2006) is the careful treatment of the Gompertz case.

So the state-space model is not a machine that removes observation error. It is a machine that lets you state your assumption about observation error and propagate it honestly. If you have replicate samples, it works well. If you are hoping it will find \(\tau^2\) for you out of sixty points, it may hand you a confident number that a slightly different starting value would have changed.

That is the same shape as every honest limit in this cluster: the naive AR(1) shortcut needs \(\Sigma\) from outside, the interaction estimates need external information to break environmental synchrony, and the correction here needs \(\tau^2\) from outside. The data does not contain its own error rate.

A checklist, and what it cannot cover

Before reporting a MAR(1) result:

  1. Series length. Report it. Below about 100 points the diagonal carries visible small-sample bias, and the fix is arithmetic, not judgement.
  2. Observation error. State what you assume \(\tau^2\) is and where the assumption came from. “We assumed none” is a defensible answer only if you say it out loud, because it sets every stability metric too low.
  3. Environmental synchrony. Check the VIF of the lagged design. Synchronised species have unidentifiable interactions and no fitting method fixes that.
  4. Multiplicity. A \(p\)-species fit tests \(p^2\) cells. Control the false discovery rate before drawing the web.
  5. Residuals. Run Ljung-Box, and understand that passing it rules out one narrow failure mode and nothing else.

The last item is the honest core of this post. Every check listed above is a falsification test: it can reject a model, and it can never confirm one. And the assumptions that matter most in a MAR(1) analysis are not on the list at all, because no diagnostic reaches them:

  • Linearity on the log scale. MAR(1) is exactly right for Gompertz dynamics and a local approximation for anything else. Certain, Barraquand and Gardmark (2018) found the coefficients recover signs and ranks under nonlinear dynamics but not values, and the residuals do not complain.
  • No missing species. An unmeasured species that drives two of your taxa shows up as an interaction between them. Nothing in the fit distinguishes that from the real thing.
  • Stationarity. A community that is on its way somewhere has no stationary distribution for the return rate to be about.

A MAR(1) analysis that reports its series length, its observation-error assumption and its multiplicity correction, and that claims signs rather than strengths, is doing what the method supports. One that quotes a return rate to four decimals from sixty uncorrected counts is quoting the sampling protocol.

References

Newsletter

Get new tutorials by email

New R and QGIS tutorials for ecologists, straight to your inbox. No spam; unsubscribe anytime.

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