Checking a coexistence analysis

coexistence

competition

model checking

ecology tutorial

Three diagnostics for an annual plant coexistence result in R: mutual invasibility, whether your design identifies self-limitation, and stability to model form.

Author

Tidy Ecology

Published

2026-07-20

The last three posts fitted competition coefficients, turned them into niche and fitness differences, and found where the coexistence verdict is fragile. This one collects the checks into a short routine you can run on your own analysis. Three questions catch most of the ways a coexistence claim goes wrong: whether a short experiment actually tested coexistence, whether your design could estimate the coefficient the verdict needs, and whether the verdict survives a reasonable change of model form.

library(ggplot2)
theme_te <- theme_minimal(base_size = 12) +
  theme(panel.grid.minor = element_blank(),
        plot.background = element_rect(fill = "white", colour = NA),
        panel.background = element_rect(fill = "white", colour = NA))

lamB <- 30; aAA <- 0.10; aBB <- 0.06; aAB <- 0.03; aBA <- 0.06
gA <- 0.60; gB <- 0.70; sA <- 0.10; sB <- 0.20

Diagnostic 1: did a short experiment really test coexistence?

Coexistence is about the long run: each species must recover when it drifts rare. A common shortcut is to grow the two species together from an even mixture and check that both are still present after a few generations. That shortcut can pass while the pair is heading for exclusion, because a rare-when-common species is not the same as a species growing from rare. Take a pair where A clearly excludes B, and watch what a mixed culture does.

lamA <- 72                                   # A strongly superior: exclusion expected
step <- function(N) {
  FA <- lamA / (1 + aAA * gA * N[1] + aAB * gB * N[2])
  FB <- lamB / (1 + aBA * gA * N[1] + aBB * gB * N[2])
  c((1 - gA) * sA * N[1] + gA * N[1] * FA,
    (1 - gB) * sB * N[2] + gB * N[2] * FB)
}
gens <- 50
M <- matrix(0, gens + 1, 2); M[1, ] <- c(20, 20)
for (t in 1:gens) M[t + 1, ] <- step(M[t, ])
Bshare <- M[, 2] / M[1, 2]                    # B relative to its starting abundance

From an even start, B does not decline politely. It climbs to 924 per cent of its starting abundance by generation 2, and is still above its starting level at generation 10 (123 per cent). Only later does the decline set in: 17 per cent at generation 20, and below one per cent only by generation 35. For an annual plant that is 35 years. A field experiment that ran for the ten or fifteen generations most experiments can manage would show B thriving and report coexistence.

The invasion criterion gives the right answer analytically. Introduce B rare against A at A’s single-species equilibrium and compute its growth rate:

Nstar <- function(lam, g, s, a) (g * lam / (1 - (1 - g) * s) - 1) / (a * g)
NAstar <- Nstar(lamA, gA, sA, aAA)
ldgr_B <- log((1 - gB) * sB + gB * lamB / (1 + aBA * gA * NAstar))

The invasion growth rate of B is -0.191, below zero, so B is excluded despite its early surge in the mixture. The check is not “are both species present after a while” but “does each species grow when rare”: compute the invasion growth rate directly, or run the mixture long enough to see the loser turn over, which for slow exclusion is longer than any real experiment.

df <- rbind(
  data.frame(gen = 0:gens, N = M[, 1], species = "A"),
  data.frame(gen = 0:gens, N = M[, 2], species = "B"))
ggplot(df, aes(gen, N, colour = species)) +
  annotate("rect", xmin = 0, xmax = 15, ymin = 1, ymax = Inf,
           fill = "grey85", alpha = 0.55) +
  annotate("text", x = 7.5, y = 1.4, label = "typical experiment",
           size = 3.3, colour = "grey30") +
  geom_line(linewidth = 0.9) +
  scale_y_log10() +
  scale_colour_manual(values = c("A" = "#275139", "B" = "#b5534e")) +
  labs(x = "generation", y = "seed density (log scale)", colour = "species") +
  theme_te + theme(legend.position = "top")

Log-scale density time series for two species from an equal start. Species A rises to a plateau; species B rises for about a dozen generations then falls steeply toward zero. A shaded vertical band over the first fifteen generations marks the reach of a typical experiment, where both species appear to persist. — Figure 1: A mixed culture of a pair that the invasion criterion rules out. Both species are shown from an even start on a log density scale. Species B rises for the first dozen generations, peaking well above its starting abundance, before collapsing toward zero. The shaded band marks the span a typical field experiment can cover, within which the pair looks like it is coexisting.

Diagnostic 2: can your design estimate self-limitation?

The verdict turns on whether intraspecific competition exceeds interspecific competition, so the intraspecific coefficient has to be estimable from your data. If conspecific density never varied, it is not, and no fitting routine will tell you the niche difference is unreliable unless you look. The test is direct: try to fit the full model on a slice of data where conspecific density is held fixed.

set.seed(207)
gp <- c(0, 2, 4, 8, 16, 32)
d <- expand.grid(Df = gp, Dh = gp); d <- d[rep(seq_len(nrow(d)), each = 6), ]
d$seeds <- rnbinom(nrow(d), mu = 45 / (1 + aAA * d$Df + aAB * d$Dh), size = 10)
additive <- d[d$Df == 0, ]                    # conspecific density fixed at zero
outcome <- tryCatch(
  nls(seeds ~ lam / (1 + aii * Df + aij * Dh), data = additive,
      start = list(lam = 40, aii = 0.05, aij = 0.05)),
  error = function(e) conditionMessage(e))
outcome

[1] "singular gradient matrix at initial parameter estimates"

The singular gradient is the design telling you that self-limitation is not identified: with conspecific density fixed, the intraspecific coefficient never enters the predictions, so the data carry no information about it. A design that varied only heterospecific density gives you an interspecific coefficient and a niche difference resting on a coefficient you never estimated. If this fit fails on your data, the coexistence question is unanswerable with that data, however good the interspecific estimates look.

Diagnostic 3: does the verdict survive a change of model form?

The invasion growth rate is computed at the resident’s equilibrium, far past the fitted densities, so the competition model’s shape out there matters. Refit with a model that fits the data comparably but bends differently, and see whether the conclusions hold. Here we fit Beverton-Holt and Ricker to the same borderline pair and compare not just the verdict but the quantities the study would report.

lamA <- 52                                    # borderline pair from the previous post
simfoc <- function(seed, lam, aii, aij) {
  set.seed(seed)
  g <- expand.grid(Df = gp, Dh = gp); g <- g[rep(seq_len(nrow(g)), each = 6), ]
  g$seeds <- rnbinom(nrow(g), mu = lam / (1 + aii * g$Df + aij * g$Dh), size = 10); g
}
A <- simfoc(111, lamA, aAA, aAB)
bh <- nls(seeds ~ lam / (1 + aii * Df + aij * Dh), data = A,
          start = list(lam = 50, aii = 0.05, aij = 0.05))
rk <- nls(seeds ~ lam * exp(-aii * Df - aij * Dh), data = A,
          start = list(lam = 50, aii = 0.02, aij = 0.02))

eq_density <- function(cf, form) {
  N <- 50
  for (i in 1:2000) {
    F <- if (form == "bh") cf["lam"] / (1 + cf["aii"] * gA * N)
         else cf["lam"] * exp(-cf["aii"] * gA * N)
    N <- (1 - gA) * sA * N + gA * N * F
  }
  gA * N                                        # germinated resident density
}
eq_bh <- eq_density(coef(bh), "bh"); eq_rk <- eq_density(coef(rk), "rk")

The two models fit the data equally well (residual standard deviation 10.6 against 10.5 seeds), yet they disagree completely on how abundant the resident is. Beverton-Holt puts species A’s equilibrium at 355 germinated plants; Ricker puts it at 9, a 39-fold difference on identical data. Both still say the pair coexists, so the qualitative verdict happens to survive. Everything quantitative about it does not: the predicted resident abundance, and with it the invasion margin and any biomass or yield the study reports, swing by more than an order of magnitude with a model choice the data cannot make.

traj <- function(cf, form, tmax = 60) {
  N <- numeric(tmax + 1); N[1] <- 50 / gA
  for (t in 1:tmax) {
    F <- if (form == "bh") cf["lam"] / (1 + cf["aii"] * gA * N[t])
         else cf["lam"] * exp(-cf["aii"] * gA * N[t])
    N[t + 1] <- (1 - gA) * sA * N[t] + gA * N[t] * F
  }
  gA * N
}
tf <- rbind(
  data.frame(gen = 0:60, N = traj(coef(bh), "bh"), model = "Beverton-Holt"),
  data.frame(gen = 0:60, N = traj(coef(rk), "rk"), model = "Ricker"))
ggplot(tf, aes(gen, N, colour = model)) +
  geom_line(linewidth = 1) +
  scale_y_log10() +
  scale_colour_manual(values = c("Beverton-Holt" = "#275139", "Ricker" = "#b5534e")) +
  labs(x = "generation", y = "resident A density (log scale)", colour = NULL) +
  theme_te + theme(legend.position = "top")

Two population trajectories on a log density scale converging to very different equilibria. The Beverton-Holt trajectory settles high; the Ricker trajectory settles far lower, more than a tenfold gap, from the same starting density. — Figure 2: Resident population of species A under two competition models fitted to the same data, iterated to equilibrium from the same start on a log density scale. Beverton-Holt settles at a high equilibrium; Ricker, whose fecundity falls faster at high density, settles more than an order of magnitude lower. The observed data cannot separate the two models, so this gap is unconstrained.

The routine

Three checks, in order. First, do not trust a short mixed culture: compute the invasion growth rate, because a pair heading for exclusion can look like coexistence for longer than an experiment runs. Second, confirm your design varied conspecific density, or the intraspecific coefficient the verdict needs is not in your data. Third, refit with a second model form and check whether the quantities you report, not just the yes-or-no verdict, hold up in the extrapolation where the invasion criterion lives. A coexistence result that clears all three is worth believing; one that clears none of them is a statement about the model, not the plants.

References

Grainger TN, Levine JM, Gilbert B 2019. Trends in Ecology and Evolution 34(10):925-935 (10.1016/j.tree.2019.05.007).

Hart SP, Freckleton RP, Levine JM 2018. Journal of Ecology 106(5):1902-1909 (10.1111/1365-2745.12954).

Levine JM, HilleRisLambers J 2009. Nature 461(7261):254-257 (10.1038/nature08251).

--- title: "Checking a coexistence analysis" description: "Three diagnostics for an annual plant coexistence result in R: mutual invasibility, whether your design identifies self-limitation, and stability to model form." date: "2026-07-20 15:00" categories: [R, coexistence, competition, model checking, ecology tutorial] image: thumbnail.png image-alt: "Two annual plant species in a mixed culture where the weaker species rises for the first dozen generations before collapsing, illustrating that a short experiment can mimic coexistence that the invasion criterion rules out." --- The last three posts fitted competition coefficients, turned them into niche and fitness differences, and found where the coexistence verdict is fragile. This one collects the checks into a short routine you can run on your own analysis. Three questions catch most of the ways a coexistence claim goes wrong: whether a short experiment actually tested coexistence, whether your design could estimate the coefficient the verdict needs, and whether the verdict survives a reasonable change of model form. ```{r setup} #| message: false #| warning: false library(ggplot2) theme_te <- theme_minimal(base_size = 12) + theme(panel.grid.minor = element_blank(), plot.background = element_rect(fill = "white", colour = NA), panel.background = element_rect(fill = "white", colour = NA)) lamB <- 30; aAA <- 0.10; aBB <- 0.06; aAB <- 0.03; aBA <- 0.06 gA <- 0.60; gB <- 0.70; sA <- 0.10; sB <- 0.20 ``` ## Diagnostic 1: did a short experiment really test coexistence? Coexistence is about the long run: each species must recover when it drifts rare. A common shortcut is to grow the two species together from an even mixture and check that both are still present after a few generations. That shortcut can pass while the pair is heading for exclusion, because a rare-when-common species is not the same as a species growing from rare. Take a pair where A clearly excludes B, and watch what a mixed culture does. ```{r excl-pair} lamA <- 72 # A strongly superior: exclusion expected step <- function(N) { FA <- lamA / (1 + aAA * gA * N[1] + aAB * gB * N[2]) FB <- lamB / (1 + aBA * gA * N[1] + aBB * gB * N[2]) c((1 - gA) * sA * N[1] + gA * N[1] * FA, (1 - gB) * sB * N[2] + gB * N[2] * FB) } gens <- 50 M <- matrix(0, gens + 1, 2); M[1, ] <- c(20, 20) for (t in 1:gens) M[t + 1, ] <- step(M[t, ]) Bshare <- M[, 2] / M[1, 2] # B relative to its starting abundance ``` ```{r excl-num} #| echo: false b10 <- 100 * Bshare[11]; b20 <- 100 * Bshare[21] bpeak <- 100 * max(Bshare); gpeak <- which.max(Bshare) - 1 gunder <- which(Bshare < 0.01)[1] - 1 ``` From an even start, B does not decline politely. It climbs to `r sprintf("%.0f", bpeak)` per cent of its starting abundance by generation `r gpeak`, and is still above its starting level at generation 10 (`r sprintf("%.0f", b10)` per cent). Only later does the decline set in: `r sprintf("%.0f", b20)` per cent at generation 20, and below one per cent only by generation `r gunder`. For an annual plant that is `r gunder` years. A field experiment that ran for the ten or fifteen generations most experiments can manage would show B thriving and report coexistence. The invasion criterion gives the right answer analytically. Introduce B rare against A at A's single-species equilibrium and compute its growth rate: ```{r excl-invasion} Nstar <- function(lam, g, s, a) (g * lam / (1 - (1 - g) * s) - 1) / (a * g) NAstar <- Nstar(lamA, gA, sA, aAA) ldgr_B <- log((1 - gB) * sB + gB * lamB / (1 + aBA * gA * NAstar)) ``` ```{r excl-inv-num} #| echo: false lb_v <- ldgr_B ``` The invasion growth rate of B is `r sprintf("%.3f", lb_v)`, below zero, so B is excluded despite its early surge in the mixture. The check is not "are both species present after a while" but "does each species grow when rare": compute the invasion growth rate directly, or run the mixture long enough to see the loser turn over, which for slow exclusion is longer than any real experiment. ```{r fig-trap} #| fig-cap: "A mixed culture of a pair that the invasion criterion rules out. Both species are shown from an even start on a log density scale. Species B rises for the first dozen generations, peaking well above its starting abundance, before collapsing toward zero. The shaded band marks the span a typical field experiment can cover, within which the pair looks like it is coexisting." #| fig-alt: "Log-scale density time series for two species from an equal start. Species A rises to a plateau; species B rises for about a dozen generations then falls steeply toward zero. A shaded vertical band over the first fifteen generations marks the reach of a typical experiment, where both species appear to persist." df <- rbind( data.frame(gen = 0:gens, N = M[, 1], species = "A"), data.frame(gen = 0:gens, N = M[, 2], species = "B")) ggplot(df, aes(gen, N, colour = species)) + annotate("rect", xmin = 0, xmax = 15, ymin = 1, ymax = Inf, fill = "grey85", alpha = 0.55) + annotate("text", x = 7.5, y = 1.4, label = "typical experiment", size = 3.3, colour = "grey30") + geom_line(linewidth = 0.9) + scale_y_log10() + scale_colour_manual(values = c("A" = "#275139", "B" = "#b5534e")) + labs(x = "generation", y = "seed density (log scale)", colour = "species") + theme_te + theme(legend.position = "top") ``` ## Diagnostic 2: can your design estimate self-limitation? The verdict turns on whether intraspecific competition exceeds interspecific competition, so the intraspecific coefficient has to be estimable from your data. If conspecific density never varied, it is not, and no fitting routine will tell you the niche difference is unreliable unless you look. The test is direct: try to fit the full model on a slice of data where conspecific density is held fixed. ```{r design-check} set.seed(207) gp <- c(0, 2, 4, 8, 16, 32) d <- expand.grid(Df = gp, Dh = gp); d <- d[rep(seq_len(nrow(d)), each = 6), ] d$seeds <- rnbinom(nrow(d), mu = 45 / (1 + aAA * d$Df + aAB * d$Dh), size = 10) additive <- d[d$Df == 0, ] # conspecific density fixed at zero outcome <- tryCatch( nls(seeds ~ lam / (1 + aii * Df + aij * Dh), data = additive, start = list(lam = 40, aii = 0.05, aij = 0.05)), error = function(e) conditionMessage(e)) outcome ``` The singular gradient is the design telling you that self-limitation is not identified: with conspecific density fixed, the intraspecific coefficient never enters the predictions, so the data carry no information about it. A design that varied only heterospecific density gives you an interspecific coefficient and a niche difference resting on a coefficient you never estimated. If this fit fails on your data, the coexistence question is unanswerable with that data, however good the interspecific estimates look. ## Diagnostic 3: does the verdict survive a change of model form? The invasion growth rate is computed at the resident's equilibrium, far past the fitted densities, so the competition model's shape out there matters. Refit with a model that fits the data comparably but bends differently, and see whether the conclusions hold. Here we fit Beverton-Holt and Ricker to the same borderline pair and compare not just the verdict but the quantities the study would report. ```{r modelform} lamA <- 52 # borderline pair from the previous post simfoc <- function(seed, lam, aii, aij) { set.seed(seed) g <- expand.grid(Df = gp, Dh = gp); g <- g[rep(seq_len(nrow(g)), each = 6), ] g$seeds <- rnbinom(nrow(g), mu = lam / (1 + aii * g$Df + aij * g$Dh), size = 10); g } A <- simfoc(111, lamA, aAA, aAB) bh <- nls(seeds ~ lam / (1 + aii * Df + aij * Dh), data = A, start = list(lam = 50, aii = 0.05, aij = 0.05)) rk <- nls(seeds ~ lam * exp(-aii * Df - aij * Dh), data = A, start = list(lam = 50, aii = 0.02, aij = 0.02)) eq_density <- function(cf, form) { N <- 50 for (i in 1:2000) { F <- if (form == "bh") cf["lam"] / (1 + cf["aii"] * gA * N) else cf["lam"] * exp(-cf["aii"] * gA * N) N <- (1 - gA) * sA * N + gA * N * F } gA * N # germinated resident density } eq_bh <- eq_density(coef(bh), "bh"); eq_rk <- eq_density(coef(rk), "rk") ``` ```{r mf-num} #| echo: false eqbh_v <- eq_bh; eqrk_v <- eq_rk; ratio_v <- eq_bh / eq_rk sdbh_v <- sd(residuals(bh)); sdrk_v <- sd(residuals(rk)) ``` The two models fit the data equally well (residual standard deviation `r sprintf("%.1f", sdbh_v)` against `r sprintf("%.1f", sdrk_v)` seeds), yet they disagree completely on how abundant the resident is. Beverton-Holt puts species A's equilibrium at `r round(eqbh_v)` germinated plants; Ricker puts it at `r round(eqrk_v)`, a `r sprintf("%.0f", ratio_v)`-fold difference on identical data. Both still say the pair coexists, so the qualitative verdict happens to survive. Everything quantitative about it does not: the predicted resident abundance, and with it the invasion margin and any biomass or yield the study reports, swing by more than an order of magnitude with a model choice the data cannot make. ```{r fig-modelform} #| fig-cap: "Resident population of species A under two competition models fitted to the same data, iterated to equilibrium from the same start on a log density scale. Beverton-Holt settles at a high equilibrium; Ricker, whose fecundity falls faster at high density, settles more than an order of magnitude lower. The observed data cannot separate the two models, so this gap is unconstrained." #| fig-alt: "Two population trajectories on a log density scale converging to very different equilibria. The Beverton-Holt trajectory settles high; the Ricker trajectory settles far lower, more than a tenfold gap, from the same starting density." traj <- function(cf, form, tmax = 60) { N <- numeric(tmax + 1); N[1] <- 50 / gA for (t in 1:tmax) { F <- if (form == "bh") cf["lam"] / (1 + cf["aii"] * gA * N[t]) else cf["lam"] * exp(-cf["aii"] * gA * N[t]) N[t + 1] <- (1 - gA) * sA * N[t] + gA * N[t] * F } gA * N } tf <- rbind( data.frame(gen = 0:60, N = traj(coef(bh), "bh"), model = "Beverton-Holt"), data.frame(gen = 0:60, N = traj(coef(rk), "rk"), model = "Ricker")) ggplot(tf, aes(gen, N, colour = model)) + geom_line(linewidth = 1) + scale_y_log10() + scale_colour_manual(values = c("Beverton-Holt" = "#275139", "Ricker" = "#b5534e")) + labs(x = "generation", y = "resident A density (log scale)", colour = NULL) + theme_te + theme(legend.position = "top") ``` ## The routine Three checks, in order. First, do not trust a short mixed culture: compute the invasion growth rate, because a pair heading for exclusion can look like coexistence for longer than an experiment runs. Second, confirm your design varied conspecific density, or the intraspecific coefficient the verdict needs is not in your data. Third, refit with a second model form and check whether the quantities you report, not just the yes-or-no verdict, hold up in the extrapolation where the invasion criterion lives. A coexistence result that clears all three is worth believing; one that clears none of them is a statement about the model, not the plants. ## Related tutorials - [Assumptions behind coexistence](../assumptions-behind-coexistence/) - [Niche and fitness differences](../niche-and-fitness-differences/) - [Fitting annual plant competition models](../annual-plant-competition-models/) - [Checking a nonlinear model](../checking-a-nonlinear-model/) ## References Grainger TN, Levine JM, Gilbert B 2019. Trends in Ecology and Evolution 34(10):925-935 (10.1016/j.tree.2019.05.007). Hart SP, Freckleton RP, Levine JM 2018. Journal of Ecology 106(5):1902-1909 (10.1111/1365-2745.12954). Levine JM, HilleRisLambers J 2009. Nature 461(7261):254-257 (10.1038/nature08251).