Co-occurrence is not interaction

JSDM

community ecology

ecology tutorial

ggplot2

Why a residual correlation between two species cannot identify competition in R: an observational equivalence, an unidentified arrow, and a missing gradient.

Author

Tidy Ecology

Published

2026-07-28

A joint model hands you a matrix of residual correlations, and the temptation is immediate. Two species correlate negatively after the environment is accounted for, so they must be competing. Two correlate positively, so one facilitates the other, or they share a mutualist.

This post takes that inference apart. Not with a simulation showing it is hard, but with an algebraic identity showing it is impossible from this data alone, followed by a measurement of how badly a single unrecorded gradient distorts the matrix.

Two species, two mechanisms, one distribution

Set up two species that share a hidden environmental gradient and do not interact at all. Each responds to the gradient with its own loading and has its own independent noise.

library(ggplot2)

theme_te <- function(base_size = 11) {
  theme_minimal(base_size = base_size) +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_line(colour = "#dad9ca", linewidth = 0.3),
          axis.title = element_text(colour = "#46604a"),
          axis.text = element_text(colour = "#5d6b61"),
          plot.title = element_text(colour = "#16241d", face = "bold"),
          plot.subtitle = element_text(colour = "#5d6b61"),
          legend.title = element_text(colour = "#46604a"),
          legend.text = element_text(colour = "#5d6b61"),
          plot.background = element_rect(fill = "white", colour = NA),
          panel.background = element_rect(fill = "white", colour = NA))
}

lam_a <- 0.9; lam_b <- 0.7
sd_a  <- 0.6; sd_b  <- 0.5

Sigma_driver <- matrix(c(lam_a^2 + sd_a^2, lam_a * lam_b,
                         lam_a * lam_b,    lam_b^2 + sd_b^2), 2, 2)
r_driver <- cov2cor(Sigma_driver)[1, 2]

Now build a second world with no hidden gradient at all. Species A varies on its own; species B responds directly to the local abundance of species A with coefficient beta, plus its own noise. This is a biotic interaction in the plainest sense: change A and B changes.

var_a <- lam_a^2 + sd_a^2
beta  <- lam_a * lam_b / var_a
var_b <- (lam_b^2 + sd_b^2) - beta^2 * var_a

Sigma_inter <- matrix(c(var_a, beta * var_a,
                        beta * var_a, beta^2 * var_a + var_b), 2, 2)

equiv_gap <- max(abs(Sigma_driver - Sigma_inter))
r_inter <- cov2cor(Sigma_inter)[1, 2]

The two covariance matrices are equal to within 1.1e-16. Both give a residual correlation of 0.6771. Since both worlds are multivariate normal with the same mean and the same covariance, they are the same distribution: every data set the first can produce, the second produces with identical probability. The likelihoods are identical, so no estimator, no prior, no amount of extra plots can separate them.

The interaction coefficient that reproduces the shared-gradient world exactly is 0.5385. An ecologist who fitted the interaction model to data from the shared-gradient world would recover that number, report it with a tight confidence interval, and be describing a mechanism that does not exist.

The arrow does not point anywhere

The equivalence is worse than a two-way ambiguity, because the interaction story itself can be told in either direction.

beta_ab <- Sigma_driver[1, 2] / Sigma_driver[1, 1]
beta_ba <- Sigma_driver[1, 2] / Sigma_driver[2, 2]

rebuild <- function(b, v_from, v_noise) {
  matrix(c(v_from, b * v_from, b * v_from, b^2 * v_from + v_noise), 2, 2)
}
gap_ab <- max(abs(rebuild(beta_ab, Sigma_driver[1, 1],
                          Sigma_driver[2, 2] - beta_ab^2 * Sigma_driver[1, 1]) -
                  Sigma_driver))
gap_ba <- max(abs(rebuild(beta_ba, Sigma_driver[2, 2],
                          Sigma_driver[1, 1] - beta_ba^2 * Sigma_driver[2, 2])[2:1, 2:1] -
                  Sigma_driver))

“A drives B” with coefficient 0.538 reproduces the covariance to within 1.1e-16. “B drives A” with coefficient 0.851 reproduces it to within 0.0e+00. Both are exact. Nothing in the joint distribution of two simultaneously measured species carries directional information, which should not surprise anyone who has met the confounding problem in a regression context: what is missing here is a design, not a better model.

The same holds for a whole community. Any positive definite covariance matrix factorises as a lower triangular matrix times its transpose, and that factorisation is a recursive system of species affecting species, in whatever order you happen to list them.

S <- 6L
L6 <- cbind(c(0.9, 0.8, -0.7, 0.2, 0.1, -0.3),
            c(0.2, -0.4, 0.5, 0.8, -0.7, 0.6))
Sig6 <- L6 %*% t(L6) + diag(rep(0.4, S))

Ch1 <- t(chol(Sig6))
gap1 <- max(abs(Ch1 %*% t(Ch1) - Sig6))

ord <- c(4, 1, 6, 2, 5, 3)
Ch2 <- t(chol(Sig6[ord, ord]))
back <- matrix(0, S, S)
back[ord, ord] <- Ch2 %*% t(Ch2)
gap2 <- max(abs(back - Sig6))

Listing the species in their original order gives a cascade that rebuilds the covariance to 2.2e-16. Listing them in a shuffled order gives a completely different cascade, with different coefficients and different arrows, that rebuilds the same covariance to 2.2e-16. There are 720 such orderings for six species and every one of them fits perfectly.

What one unrecorded gradient does

The algebra says the mechanisms cannot be told apart. The simulation below says something more practical: the most likely source of residual correlation in a real survey is not biotic at all, and its effect is large.

Ten species respond to three environmental gradients. Two are measured. The third is real, ordinary and unrecorded: soil depth, say, or grazing history. There is no interaction anywhere in the code.

set.seed(385)
n <- 300L
Sp <- 10L
x1 <- runif(n, -2, 2)
x2 <- runif(n, -2, 2)
x3 <- runif(n, -2, 2)

b1 <- rnorm(Sp, 0, 0.8)
b2 <- rnorm(Sp, 0, 0.6)
b3 <- c(1.2, 1.1, 0.9, 0.8, 0.7, -0.9, -1.0, -1.1, 0.1, -0.2)
a0 <- rnorm(Sp, 2, 0.5)

Y <- outer(rep(1, n), a0) + outer(x1, b1) + outer(x2, b2) + outer(x3, b3) +
     matrix(rnorm(n * Sp, 0, 0.7), n, Sp)

lower <- lower.tri(diag(Sp))
R_short <- cor(residuals(lm(Y ~ x1 + x2)))
R_full  <- cor(residuals(lm(Y ~ x1 + x2 + x3)))

z_test <- function(r, nn) 2 * pnorm(-abs(atanh(r) * sqrt(nn - 4)))
flag_short <- sum(z_test(R_short[lower], n) < 0.05)
flag_full  <- sum(z_test(R_full[lower], n) < 0.05)
npairs <- sum(lower)

Fitting the two measured gradients leaves an average absolute residual correlation of 0.494, with the strongest pair at 0.796. Tested pair by pair at the five per cent level, 39 of the 45 pairs come out significant. Add the third gradient and the average falls to 0.055, the strongest pair to 0.140, and the count of significant pairs to 4.

Two overlapping density curves of pairwise residual correlations. The curve for the model missing a gradient is broad and spans minus one to plus one; the curve for the full model is narrow and centred on zero. — Figure 1: The 45 pairwise residual correlations from the same ten species, with and without one unrecorded environmental gradient in the model. The community is identical and contains no biotic interaction.

A proxy is not a substitute

The obvious response is that nobody measures everything, so use whatever you have. The question is how good a stand-in has to be before it helps.

set.seed(3851)
rhos <- c(0, 0.3, 0.6, 0.9, 0.99)
proxy_mean <- sapply(rhos, function(rho) {
  prox <- rho * x3 + sqrt(1 - rho^2) * rnorm(n)
  Rp <- cor(residuals(lm(Y ~ x1 + x2 + prox)))
  mean(abs(Rp[lower]))
})
base_mean <- mean(abs(R_short[lower]))
removed <- 1 - proxy_mean / base_mean

Line chart with proxy quality on the horizontal axis and mean absolute residual correlation on the vertical axis. The line is nearly flat at first and falls steeply only above a proxy correlation of 0.9. — Figure 2: Mean absolute residual correlation as a stand-in variable is made more and more similar to the unrecorded gradient. The spurious structure survives until the proxy is nearly perfect.

A proxy correlating 0.6 with the truth, which would count as a strong covariate in most surveys, removes only 25 per cent of the spurious structure. At 0.9 it removes 63 per cent. The relationship is steeply nonlinear because what survives is the component of the gradient the proxy misses, and that component still drives every species.

What survives

The residual correlation matrix does carry information, and it is worth being precise about what kind.

It is evidence that the fitted environmental model is incomplete. That is a real and useful finding: a strongly structured residual correlation matrix says your covariates are not enough, in the same way a patterned residual plot does in an ordinary regression.

It is a good basis for prediction, as the first post in this series measured. Conditional prediction does not require knowing why species covary, only that they do.

It is not evidence of interaction, and adding species traits, phylogeny or spatial random effects does not change that. Those additions can explain the correlation better; none of them can rule out an unrecorded gradient, because the unrecorded gradient produces exactly the same distribution.

The honest reporting standard follows from that. Call the entries associations, not interactions. State what environment was measured and what was not. If the question really is whether species A suppresses species B, the answer comes from a removal experiment or a designed gradient, not from a bigger observational matrix.

References

Blanchet, Cazelles, Gravel 2020 Ecology Letters 23(7):1050-1063 (10.1111/ele.13525)

Dormann, Bobrowski, Dehling, Harris, Hartig, Lischke, Moretti, Pagel, Pinkert, Schleuning, Schmidt, Sheppard, Steinbauer, Zeuss, Kraan 2018 Global Ecology and Biogeography 27(9):1004-1016 (10.1111/geb.12759)

Poggiato, Munkemuller, Bystrova, Arbel, Clark, Thuiller 2021 Trends in Ecology and Evolution 36(5):391-401 (10.1016/j.tree.2021.01.002)

Zurell, Pollock, Thuiller 2018 Ecography 41(11):1812-1819 (10.1111/ecog.03315)

--- title: "Co-occurrence is not interaction" description: "Why a residual correlation between two species cannot identify competition in R: an observational equivalence, an unidentified arrow, and a missing gradient." date: "2026-07-28 10:00" categories: [R, JSDM, community ecology, ecology tutorial, ggplot2] image: thumbnail.png image-alt: "Two overlapping density curves of pairwise residual correlations, one wide and spread out, the other narrow and centred on zero." --- A joint model hands you a matrix of residual correlations, and the temptation is immediate. Two species correlate negatively after the environment is accounted for, so they must be competing. Two correlate positively, so one facilitates the other, or they share a mutualist. This post takes that inference apart. Not with a simulation showing it is hard, but with an algebraic identity showing it is impossible from this data alone, followed by a measurement of how badly a single unrecorded gradient distorts the matrix. ## Two species, two mechanisms, one distribution Set up two species that share a hidden environmental gradient and do not interact at all. Each responds to the gradient with its own loading and has its own independent noise. ```{r setup} #| message: false #| warning: false library(ggplot2) theme_te <- function(base_size = 11) { theme_minimal(base_size = base_size) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#dad9ca", linewidth = 0.3), axis.title = element_text(colour = "#46604a"), axis.text = element_text(colour = "#5d6b61"), plot.title = element_text(colour = "#16241d", face = "bold"), plot.subtitle = element_text(colour = "#5d6b61"), legend.title = element_text(colour = "#46604a"), legend.text = element_text(colour = "#5d6b61"), plot.background = element_rect(fill = "white", colour = NA), panel.background = element_rect(fill = "white", colour = NA)) } lam_a <- 0.9; lam_b <- 0.7 sd_a <- 0.6; sd_b <- 0.5 Sigma_driver <- matrix(c(lam_a^2 + sd_a^2, lam_a * lam_b, lam_a * lam_b, lam_b^2 + sd_b^2), 2, 2) r_driver <- cov2cor(Sigma_driver)[1, 2] ``` Now build a second world with no hidden gradient at all. Species A varies on its own; species B responds directly to the local abundance of species A with coefficient `beta`, plus its own noise. This is a biotic interaction in the plainest sense: change A and B changes. ```{r interaction} var_a <- lam_a^2 + sd_a^2 beta <- lam_a * lam_b / var_a var_b <- (lam_b^2 + sd_b^2) - beta^2 * var_a Sigma_inter <- matrix(c(var_a, beta * var_a, beta * var_a, beta^2 * var_a + var_b), 2, 2) equiv_gap <- max(abs(Sigma_driver - Sigma_inter)) r_inter <- cov2cor(Sigma_inter)[1, 2] ``` The two covariance matrices are equal to within `r sprintf("%.1e", equiv_gap)`. Both give a residual correlation of `r sprintf("%.4f", r_driver)`. Since both worlds are multivariate normal with the same mean and the same covariance, they are the *same distribution*: every data set the first can produce, the second produces with identical probability. The likelihoods are identical, so no estimator, no prior, no amount of extra plots can separate them. The interaction coefficient that reproduces the shared-gradient world exactly is `r sprintf("%.4f", beta)`. An ecologist who fitted the interaction model to data from the shared-gradient world would recover that number, report it with a tight confidence interval, and be describing a mechanism that does not exist. ## The arrow does not point anywhere The equivalence is worse than a two-way ambiguity, because the interaction story itself can be told in either direction. ```{r direction} beta_ab <- Sigma_driver[1, 2] / Sigma_driver[1, 1] beta_ba <- Sigma_driver[1, 2] / Sigma_driver[2, 2] rebuild <- function(b, v_from, v_noise) { matrix(c(v_from, b * v_from, b * v_from, b^2 * v_from + v_noise), 2, 2) } gap_ab <- max(abs(rebuild(beta_ab, Sigma_driver[1, 1], Sigma_driver[2, 2] - beta_ab^2 * Sigma_driver[1, 1]) - Sigma_driver)) gap_ba <- max(abs(rebuild(beta_ba, Sigma_driver[2, 2], Sigma_driver[1, 1] - beta_ba^2 * Sigma_driver[2, 2])[2:1, 2:1] - Sigma_driver)) ``` "A drives B" with coefficient `r sprintf("%.3f", beta_ab)` reproduces the covariance to within `r sprintf("%.1e", gap_ab)`. "B drives A" with coefficient `r sprintf("%.3f", beta_ba)` reproduces it to within `r sprintf("%.1e", gap_ba)`. Both are exact. Nothing in the joint distribution of two simultaneously measured species carries directional information, which should not surprise anyone who has met the confounding problem in a regression context: what is missing here is a design, not a better model. The same holds for a whole community. Any positive definite covariance matrix factorises as a lower triangular matrix times its transpose, and that factorisation *is* a recursive system of species affecting species, in whatever order you happen to list them. ```{r cascade} S <- 6L L6 <- cbind(c(0.9, 0.8, -0.7, 0.2, 0.1, -0.3), c(0.2, -0.4, 0.5, 0.8, -0.7, 0.6)) Sig6 <- L6 %*% t(L6) + diag(rep(0.4, S)) Ch1 <- t(chol(Sig6)) gap1 <- max(abs(Ch1 %*% t(Ch1) - Sig6)) ord <- c(4, 1, 6, 2, 5, 3) Ch2 <- t(chol(Sig6[ord, ord])) back <- matrix(0, S, S) back[ord, ord] <- Ch2 %*% t(Ch2) gap2 <- max(abs(back - Sig6)) ``` Listing the species in their original order gives a cascade that rebuilds the covariance to `r sprintf("%.1e", gap1)`. Listing them in a shuffled order gives a completely different cascade, with different coefficients and different arrows, that rebuilds the same covariance to `r sprintf("%.1e", gap2)`. There are `r factorial(S)` such orderings for six species and every one of them fits perfectly. ## What one unrecorded gradient does The algebra says the mechanisms cannot be told apart. The simulation below says something more practical: the most likely source of residual correlation in a real survey is not biotic at all, and its effect is large. Ten species respond to three environmental gradients. Two are measured. The third is real, ordinary and unrecorded: soil depth, say, or grazing history. There is no interaction anywhere in the code. ```{r missing} set.seed(385) n <- 300L Sp <- 10L x1 <- runif(n, -2, 2) x2 <- runif(n, -2, 2) x3 <- runif(n, -2, 2) b1 <- rnorm(Sp, 0, 0.8) b2 <- rnorm(Sp, 0, 0.6) b3 <- c(1.2, 1.1, 0.9, 0.8, 0.7, -0.9, -1.0, -1.1, 0.1, -0.2) a0 <- rnorm(Sp, 2, 0.5) Y <- outer(rep(1, n), a0) + outer(x1, b1) + outer(x2, b2) + outer(x3, b3) + matrix(rnorm(n * Sp, 0, 0.7), n, Sp) lower <- lower.tri(diag(Sp)) R_short <- cor(residuals(lm(Y ~ x1 + x2))) R_full <- cor(residuals(lm(Y ~ x1 + x2 + x3))) z_test <- function(r, nn) 2 * pnorm(-abs(atanh(r) * sqrt(nn - 4))) flag_short <- sum(z_test(R_short[lower], n) < 0.05) flag_full <- sum(z_test(R_full[lower], n) < 0.05) npairs <- sum(lower) ``` Fitting the two measured gradients leaves an average absolute residual correlation of `r sprintf("%.3f", mean(abs(R_short[lower])))`, with the strongest pair at `r sprintf("%.3f", max(abs(R_short[lower])))`. Tested pair by pair at the five per cent level, `r flag_short` of the `r npairs` pairs come out significant. Add the third gradient and the average falls to `r sprintf("%.3f", mean(abs(R_full[lower])))`, the strongest pair to `r sprintf("%.3f", max(abs(R_full[lower])))`, and the count of significant pairs to `r flag_full`. ```{r fig-shift} #| echo: false #| fig-cap: "The 45 pairwise residual correlations from the same ten species, with and without one unrecorded environmental gradient in the model. The community is identical and contains no biotic interaction." #| fig-alt: "Two overlapping density curves of pairwise residual correlations. The curve for the model missing a gradient is broad and spans minus one to plus one; the curve for the full model is narrow and centred on zero." dfd <- data.frame( r = c(R_short[lower], R_full[lower]), model = rep(c("two gradients measured", "all three gradients measured"), each = npairs)) dfd$model <- factor(dfd$model, levels = c("two gradients measured", "all three gradients measured")) ggplot(dfd, aes(r, fill = model, colour = model)) + geom_density(alpha = 0.45, linewidth = 0.7, bw = 0.08) + geom_vline(xintercept = 0, colour = "#5d6b61", linetype = "dashed") + scale_fill_manual(values = c("#c9b458", "#275139"), name = NULL) + scale_colour_manual(values = c("#b5534e", "#1d5b4e"), name = NULL) + labs(x = "pairwise residual correlation", y = "density", title = "One missing gradient, forty-five spurious associations", subtitle = "Ten species, 300 plots, no biotic interaction in the simulation") + theme_te() + theme(legend.position = "bottom") ``` ## A proxy is not a substitute The obvious response is that nobody measures everything, so use whatever you have. The question is how good a stand-in has to be before it helps. ```{r proxy} set.seed(3851) rhos <- c(0, 0.3, 0.6, 0.9, 0.99) proxy_mean <- sapply(rhos, function(rho) { prox <- rho * x3 + sqrt(1 - rho^2) * rnorm(n) Rp <- cor(residuals(lm(Y ~ x1 + x2 + prox))) mean(abs(Rp[lower])) }) base_mean <- mean(abs(R_short[lower])) removed <- 1 - proxy_mean / base_mean ``` ```{r fig-proxy} #| echo: false #| fig-cap: "Mean absolute residual correlation as a stand-in variable is made more and more similar to the unrecorded gradient. The spurious structure survives until the proxy is nearly perfect." #| fig-alt: "Line chart with proxy quality on the horizontal axis and mean absolute residual correlation on the vertical axis. The line is nearly flat at first and falls steeply only above a proxy correlation of 0.9." dfx <- data.frame(rho = rhos, m = proxy_mean) ggplot(dfx, aes(rho, m)) + geom_line(colour = "#275139", linewidth = 0.9) + geom_point(colour = "#275139", size = 2.8) + scale_x_continuous(breaks = rhos) + labs(x = "correlation between the proxy and the true gradient", y = "mean absolute residual correlation", title = "A good proxy is not good enough", subtitle = "Same ten species, same 300 plots") + theme_te() ``` A proxy correlating `r sprintf("%.1f", rhos[3])` with the truth, which would count as a strong covariate in most surveys, removes only `r sprintf("%.0f", 100 * removed[3])` per cent of the spurious structure. At `r sprintf("%.1f", rhos[4])` it removes `r sprintf("%.0f", 100 * removed[4])` per cent. The relationship is steeply nonlinear because what survives is the component of the gradient the proxy misses, and that component still drives every species. ## What survives The residual correlation matrix does carry information, and it is worth being precise about what kind. It is evidence that the fitted environmental model is incomplete. That is a real and useful finding: a strongly structured residual correlation matrix says your covariates are not enough, in the same way a patterned residual plot does in an ordinary regression. It is a good basis for prediction, as the first post in this series measured. Conditional prediction does not require knowing why species covary, only that they do. It is not evidence of interaction, and adding species traits, phylogeny or spatial random effects does not change that. Those additions can explain the correlation better; none of them can rule out an unrecorded gradient, because the unrecorded gradient produces exactly the same distribution. The honest reporting standard follows from that. Call the entries associations, not interactions. State what environment was measured and what was not. If the question really is whether species A suppresses species B, the answer comes from a removal experiment or a designed gradient, not from a bigger observational matrix. ## Related tutorials - [Joint species distribution models in R](../joint-species-distribution-models/) - where these correlations come from. - [Checking a joint species distribution model](../checking-a-jsdm/) - the sampling noise on top of the identification problem. - [Confounding and backdoor adjustment](../confounding-and-backdoor-adjustment/) - the same missing-variable problem in a regression setting. - [Two-species occupancy and interactions](../two-species-occupancy-interactions/) - the detection-model version of the same question. ## References Blanchet, Cazelles, Gravel 2020 Ecology Letters 23(7):1050-1063 (10.1111/ele.13525) Dormann, Bobrowski, Dehling, Harris, Hartig, Lischke, Moretti, Pagel, Pinkert, Schleuning, Schmidt, Sheppard, Steinbauer, Zeuss, Kraan 2018 Global Ecology and Biogeography 27(9):1004-1016 (10.1111/geb.12759) Poggiato, Munkemuller, Bystrova, Arbel, Clark, Thuiller 2021 Trends in Ecology and Evolution 36(5):391-401 (10.1016/j.tree.2021.01.002) Zurell, Pollock, Thuiller 2018 Ecography 41(11):1812-1819 (10.1111/ecog.03315)