<figure><img src="image_1.png" /></figure>

HAL
open science

# Metacommunity dynamics and the detection of species associations in co-occurrence analyses: Why patch disturbance matters

Vincent Calcagno, Nik J. Cunniffe, Frédéric Marie Hamelin

► To cite this version:

Vincent Calcagno, Nik J. Cunniffe, Frédéric Marie Hamelin. Metacommunity dynamics and the detection of species associations in co-occurrence analyses: Why patch disturbance matters. Functional Ecology, 2022, 36 (6), pp.1483-1499. 10.1111/1365-2435.14047. hal-03713124

HAL Id: hal-03713124

https://hal.science/hal-03713124

Submitted on 23 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Metacommunity dynamics and the detection of species associations in co-occurrence analyses: Why patch disturbance matters

Vincent Calcagno<fnref n="1" /> ID

Nik J. Cunniffe<fnref n="2" /> ID

Frédéric M. Hamelin<fnref n="3" /> ID

**Correspondence**

Nik J. Cunniffe

Email: njc1001@cam.ac.uk

**Handling Editor:** Raul Costa-Pereira

## Abstract

1. Many statistical methods attempt to detect species associations—and so infer interspecific interactions—from species co-occurrence patterns. Habitat heterogeneity and out-of-equilibrium colonization histories are well recognized as potentially causing species associations, even when interactions are absent. The potential for patch disturbance, a classical component of metacommunity dynamics, to also drive spurious species associations has however been overlooked.

2. Using a new general metacommunity model, we derive mathematical predictions regarding how patch disturbance would affect the patterns of species associations detected in 'null' co-occurrence matrices. We also conduct numerical simulations to test our predictions and to compare the performance of several widespread statistical methods, including direct tests of pairwise independence, matrix permutation approaches and joint species distribution modelling.

3. We show how classical metacommunity dynamics can produce statistical associations, both positive and negative, even when species do not interact, when there is no habitat heterogeneity, and at equilibrium. This occurs as soon as there is some rate of patch disturbance (i.e. simultaneous extinction of several species in a patch) and/or a finite life span of patches, a common feature of a broad range of plant, animal or microbial systems.

4. Patch disturbance can compromise species co-occurrence analyses and cause the artefactual detection of species associations if not taken into account. Including patch age (i.e. the time since the last patch disturbance event) as a covariate in a joint species distribution model can resolve the artefact. However, this requires additional data that often are not available in practice. We argue that the consequences of patch disturbance should not be underestimated when analysing species distribution patterns in metacommunity-like systems.

## KEYWORDS

Hmsc, independence, interactions, JSDM, metapopulation, permutations, Sim2, Sim9

<footnote n="1">¹Institut Sophia Agrobiotech, Université Côte d'Azur, INRAE, CNRS, Sophia-Antipolis, France</footnote>

<footnote n="2">²Department of Plant Sciences, University of Cambridge, Cambridge, United Kingdom</footnote>

<footnote n="3">³Institut Agro, Univ Rennes, INRAE, IGEPP, Rennes, France</footnote>
# 1 | INTRODUCTION

Community ecology has long been concerned with inferring interactions between species (e.g. competition or facilitation) from co-occurrence data (Caswell, 1976; Cohen, 1970; Connor & Simberloff, 1979; Diamond, 1975; Forbes, 1907). Recently, metagenomics approaches have renewed this use of co-occurrence patterns as a means to infer species interaction networks, in environments ranging from soils to the human microbiome (Barberán et al., 2012; Faust & Raes, 2012). If species do not interact, intuition suggests the proportion of patches (or hosts) where species co-occur should be the product of the proportions of patches occupied by each species. Testing for pairwise independence, that is, looking for statistical patterns of association, consequently underpins most methods for inferring species interactions from presence-absence data (Gotelli, 2000; Gotelli & Ulrich, 2012; Ovaskainen & Abrego, 2020). Positive associations (excess of co-occurrences) may indicate facilitation, whereas negative associations (deficit of co-occurrences) may indicate competition.

However, it is well known that several pitfalls can compromise this approach, since correlation is not equivalent to interaction (Barner et al., 2018; Blanchet et al., 2020; Molina & Stone, 2020). As a simple example, patch heterogeneity and differential habitat preferences among species (i.e. environmental filtering) can cause species associations, even in the absence of any interaction between the pairs of species involved. For instance, if some patch types host more species than others because of their larger size or more favourable conditions, the occurrence of a species in a patch increases the odds of occurrence of other species, simply because it makes it more likely the patch is of higher quality. This would generate positive association signals if not controlled for. Succession or non-equilibrium dynamics, that is, species expanding through space, concomitantly or differentially, are also expected to cause similar patterns (D'Amen et al., 2018). More recent statistical and logical frameworks incorporate and test for alternative mechanisms to species interactions (Blois et al., 2014; D'Amen et al., 2018; Ulrich et al., 2017). However, it is commonly expected that in the absence of any of the above confounding factors and others (such as dispersal limitation), or after having appropriately controlled for them, classical extinction/recolonization dynamics (metacommunity dynamics; Leibold et al., 2004) would not introduce particular species associations, unless ecological interactions are actually occurring (Opedal, von Numers, et al., 2020).

Here we show that a specific component of metacommunity dynamics, patch disturbance, has consequences that have been overlooked so far. Patch disturbance here means the stochastic occurrence of patch extinctions (Hastings, 1980, 2003). Patch disturbance events denote extreme but possibly frequent events, such as fires, droughts, floods and others, causing the local extinction of all, or perhaps some subset of all, species from a patch (Leibold et al., 2004; Sousa, 1984). This can also correspond to the actual destruction and disappearance of a patch, for instance in ecosystems managed by humans, such as agricultural fields or forest stands (Ovaskainen et al., 2017), or in systems for which 'patches' are

actually hosts with a finite life span, such as parasite or microbiote communities (Friedman & Alm, 2012; Hamelin et al., 2019). Patch disturbances are ubiquitous in natural and anthropogenic systems (Jentsch & White, 2019; Pickett & White, 2013). Although naturally regarded as a part of metapopulation and metacommunity frameworks (e.g. Calcagno et al., 2011; Hastings, 1980, 2003; Leibold et al., 2004; Levin & Paine, 1974), patch disturbance is often omitted (e.g. Slatkin, 1974) or confounded with species-specific extinctions (e.g. Thompson et al., 2020). Although this may seem a benign mathematical simplification, we will show here that this omission can have non-trivial implications.

We first introduce a general metacommunity model of patch dynamics, that allows for patch disturbance to occur at some, potentially time-varying, rate, and in which patches can have a finite life span. The model generalizes many existing particular models, and is 'null' in the sense that it lacks the ingredients regarded as causing species statistical associations: species do not interact, the habitat is completely homogeneous and all species disperse uniformly with no patch preferences. We use this model to generate null sample co-occurrence matrices, drawn from the steady-state occupancies predicted from the model. We then apply different approaches as often are used to detect species associations on these null synthetic data. The methods considered include direct statistical tests of independence (e.g. Veech, 2013), different permutation-based matrix analyses (e.g. Gotelli, 2000) and joint species distribution modelling (e.g. Ovaskainen & Abrego, 2020). In each case, we first derive mathematical predictions concerning the expected performance of the method in question, and then test these predictions with numerical simulations.

We show that the occurrence of patch disturbance will yield characteristic patterns of spurious species associations, for all methods that do not explicitly model the effect of patch age (time since last disturbance event) as a covariate. We show that a quantity we term species 'fastness', that captures how quickly species occurrence probability recovers after disturbance under the effect of recolonization/extinction dynamics, is a key determinant of the expected patterns of spurious associations. We show that the use of joint species distribution modelling with patch age as an environmental covariate (e.g. Ovaskainen et al., 2017) is the only way to achieve satisfactory performance (absence of spurious signals). Since patch age would often be unknown, we explore whether patch richness (defined as the number of species in a patch) could suffice as a proxy for patch age, but we find it does not. We suggest the implications of patch disturbance as a constituent of metacommunity dynamics should not be underestimated, and that explicit data on patch age would generally be required for accurate species co-occurrence analysis in metacommunity-like systems.

# 2 | MATERIALS AND METHODS

This study is based entirely upon theoretical and numerical investigation of a mathematical model and did not generate or use
experimental data. Therefore, neither ethical approval nor site per-mission were required.

## 2.1 | A general metacommunity model

We consider a metacommunity model describing the extinction-recolonization dynamics of s non-interacting species over a large number of identical patches (Caswell, 1976). Each species (i) has its own colonization rate $c_i$ per occupied patch, as well as a constant rate of immigration from outside the metacommunity, $m_i$ (Hastings, 1987). Any species can go extinct in a patch, at some species-dependent rate $e_i$. In addition, patches can undergo catastrophic disturbances, after which all species in the patch immediately go extinct, irrespective of species composition (Calcagno et al., 2011; Hastings, 1980; Leibold et al., 2004). Such catastrophes occur at rate $\mu_x$, where x is the ‘age’ of a patch, that is, the time since its appearance or since it experienced the last catastrophic disturbance event (Hastings, 2003; Levin & Paine, 1974). Finally, there can be some maximal patch age X after which a catastrophic disturbance systematically occurs (Olivieri et al., 1995). The model is illustrated in Figure 1.

Let $p_{i,x,t}$ be the fraction of patches that have age x and are occupied by species i at time t (Hastings, 1991). The fraction of patches occupied by species i is $p_{i,*}$, and the fraction of patches with age x at time t is $p_{i,x,t}$. The fraction of patches with age x that

are not yet occupied by species i is therefore $p_{i,x,t} - p_{i,x,t}$. We refer to Appendix S1 (Section S1-1.1) for more technical definitions. The general metacommunity model can be written as, for all i = 1, 2, ..., s,

$$ \frac{\partial p_{i,x,t}}{\partial x} + \frac{\partial p_{i,x,t}}{\partial t} = -(\mu_x + e_i)p_{i,x,t} + (c_i p_{i,*} + m_i)(p_{i,x,t} - p_{i,x,t}). \quad (1) $$

This framework generalizes classical metacommunity models (e.g. Bell, 2001; Cohen, 1970; Hanski, 1982; Hanski & Gyllenberg, 1997; Hastings, 1980, 1987; Slatkin, 1974; Tilman, 1994; Tilman et al., 1994), deriving from MacArthur and Wilson (1963)'s 'mainland-island' model ($c_i = 0$), and from Levins (1969)'s metapopulation model ($m_i = 0$). Such metacommunity models, when they consider patch catastrophic extinctions at all, assume a constant disturbance rate $\mu_x = \mu$, but this need not be the case. The rate of disturbance might be increasing with patch age, producing an accelerating failure time model, for instance if a patch is a host that ages and suffers higher mortality with aging; this is termed as Type 1 survivorship curve in Begon and Townsend (2020). Conversely, if patches have some variability in their risk of disturbance, younger patches would go extinct at a relatively high rate, whereas older patches are relatively more disturbance resistant, causing a decline of $\mu_x$ with patch age (Type 3 survivorship curve). Another common situation is when patches have a finite lifetime and systematically get destroyed after some

<figure><img src="image_4.png" /><figcaption>FIGURE 1 Graphical summary of our model and approach. Top: null metacommunity model, its assumptions and the processes it describes. Note that patch disturbances may occur at a certain rate, possibly dependent on the time since the last disturbance. Alternatively, they may occur after some prescribed amount of time (introducing a maximum life span of patches). Bottom: from the metacommunity model, using specific parameter values, we can generate sample co-occurrence matrices of arbitrary size. We can then subject these matrices to standard analytical methods to search for signals of species association. Considering that our model assumes no species interaction, no habitat heterogeneity or differences in habitat selection among species, we expect no significant association for any pair of species. Any such association would be a spurious signal (type I error, i.e. false positive)</figcaption></figure>
prescribed duration X (Pickett & White, 2013), as may happen for instance in agricultural settings where crops are harvested after some fixed time. Our framework encompasses all such situations. As special cases, the classical ‘mainland-island’ and ‘Levins’ models correspond to $c_i = 0$ and $m_i = 0$, respectively (Gotelli, 1991; Hanski & Gyllenberg, 1993), together with $X \to \infty$ and $\mu_x = \mu$, a constant (see Appendix S1, Section S1-2).

## 2.2 | Steady-state occupancies and co-occurrence patterns

At steady state, we can drop the t subscripts. The overall occupancy of species i, $p_{i,*}$ can be derived explicitly in the special cases of the mainland-island or Levins models, for which simple expressions exist (see Appendix S1, Section S1-2.1). For instance, $p_{i,*} = 1 - (e_i + \mu)/c_i$ in the Levins model. In general, no such solution exists when the rate of patch disturbance depends on age, but $p_{i,*}$ can be computed numerically, for example, using a recursive algorithm outlined in Appendix S1, Section S1-1.4.

It would be possible to derive a separate equation for the dynamics of every possible patch state, that is, the fraction of patches that are occupied by any given combination of species, and solve these (as done in Appendix S1, Section S1-2.2, for two species). However, this is computationally inefficient. To obtain random co-occurrence matrices drawn from the above model, we instead use a Monte-Carlo procedure that generates the desired number of patches and draws their species composition. We provide a comprehensive set of R functions that make it easy to compute the steady-state occupancy of each species and generate random species co-occurrence matrices from the general model (1), for any number of species, patches and any shape of the disturbance function $\mu_x$ (see Appendix S2). We will use this method to generate null co-occurrence matrices that involve no species interactions, just metacommunity dynamics (Figure 1).

## 2.3 | Introducing relative distribution profiles and ‘fastness’

As illustrated in Figure 2a, at steady state, there exists a stable distribution of patch ages, and species occupancies vary as a function of patch age. We here derive a quantity $\pi_{i/x}$ representing how species i is distributed over patch ages. We term this quantity the relative distribution profile of a species (Figure 2b). Values smaller (larger) than one imply the species i is rarer (more frequent) in patches of age x, relative to its overall occupancy $p_{i,*}$.

We first define the fraction of patches of age x that are occupied by species i as $p_{i|x} = p_{i,x}/p_{*,x}$, for which we obtain an explicit expression from model 1 at steady state (Appendix S1, Section S1-1.3). The relative distribution profile of species i is then derived as the fraction of patches of age x that are occupied by species i, relative to its overall probability of occupancy $p_{i,*}$.

$$ \pi_{i/x} = \frac{p_{i|x}}{p_{i,*}} = \frac{1}{p_{i,*}} \frac{c_i p_{i,*} + m_i}{c_i p_{i,*} + m_i + e_i} (1 - \exp(-(c_i p_{i,*} + m_i + e_i)x)). \quad (2) $$

As shown in Figure 2b, the mean value of the profile, for any species, is equal to one by definition (Appendix S1, Section S1-1.5).

To facilitate interpretation, we will discriminate species based on their ‘fastness’: fast species are those which have very flat relative distribution profiles that quickly jump from zero and saturate at a value close to one (Figure 2b). Slow species, on the contrary, are those which have relative distribution profiles that slowly increase from zero and reach much larger values as the patch age gets larger.

Formally, we will quantify the fastness of a relative distribution profile as its variance over patch ages, $Var(\pi_{i/x})$. A maximally fast species has $\pi_{i/x} = 1$ for all x, that is, zero variance (Figure 2). Conversely, a very slow species has a relative distribution that gradually climbs to very large values, implying a huge variance. It can be shown that the variance of a species is also strictly decreasing with the initial slope of its relative distribution profile (Appendix S1, Section S1-1.7). The variance of a species profile is thus a metric of ‘slowness’, and its inverse is a metric of fastness. Importantly, the fastness of a species is not directly related to its overall occupancy in the metacommunity.

## 2.4 | Methods to test for species associations in co-occurrence matrices

Methods to infer species associations from co-occurrence patterns are all based on the probability of co-occurrence, that is, the fraction of patches in which both of the two species are found. Various methods exist to determine whether the observed probability deviates from statistical independence, under some appropriate ‘null’ model. The probability of co-occurrence of two species i and j, that is, the overall fraction of patches in which the two species are found, will be denoted $q_{i,j,*}$ (at steady state).

A straightforward null expectation is that of species independence: if species are not interacting, they should be distributed independently across patches, and thus the value of $q_{i,j,*}$ should be compared to the null value $p_{i,*}p_{j,*}$, using a standard test for contingency tables (Veech, 2013). In this article, we will use pairwise Fisher tests for independence as an exemplar of this type of method.

The prediction of species independence is straightforward and intuitive, but is known to suffer from several limitations. To circumvent these, more sophisticated null models are often preferred. A widespread approach consists of using permutation schemes intended to break species associations in co-occurrence matrices while retaining important species and patch differences (Gotelli, 2000). In this context, the C-score (Stone & Roberts, 1990) is the metric most commonly used to quantify the tendency of species to be segregated or aggregated in co-occurrence matrices. The partial C-score between two species, $C_{i,j}$, can be expressed as

$$ C_{i,j} = N^2(p_{i,*} - q_{i,j,*})(p_{j,*} - q_{i,j,*}), \quad (3) $$

with N the total number of sites (patches) in the matrix.
<figure><img src="image_5.png" /><figcaption>FIGURE 2 Example steady-state solutions of Equation (1), obtained by numerical simulation. Three species differing in fastness are used for illustration. (a) In this example, the stable patch-age distribution was exponential with mean 10 units of time (top curve). For any patch age x, each species has some patch occupancy $p_{i|x}$ (bottom curves), leaving a certain fraction of patches empty. Note that in this example the fastest species are more abundant overall, but this need not be the case. (b) Relative distribution profiles ($\pi_{i|x}$). The faster species has steeper initial slope and attains lower asymptotic values. As a species gets faster and faster, its profile converges to the horizontal dashed line (i.e. a value of 1 for all patch ages). As shown in the inset, faster species have a smaller variance in $\pi_{i|x}$ values, and as a species gets faster and faster, the distribution converges to a Dirac delta function at value 1 (vertical dashed line). Parameters: $m_i = 0.01$, $\mu = 0.1$, $X = \infty$ and $e_i = (0.05, 0.1, 0.2)$ and $c_i = (0.2, 0.3, 1)$ for slow, intermediate and fast species, respectively</figcaption></figure>

Values of $C_{ij}$ larger than expected indicate segregation (e.g. competition) between the two species, whereas values smaller than expected indicate preferential association (e.g. facilitative interactions). Several null models have been proposed to determine the expected $C_{ij}$ value and test for deviations from it. The most classical methods are the so-called fixed-equiprobable (or sim2) and fixed-fixed (or sim9) permutation schemes that reshuffle matrix values while keeping the row sums fixed, or both the row and column sums fixed, respectively (Gotelli & Ulrich, 2012; Münkemüller et al., 2020). We will consider both approaches.

More recently, flexible methods relying on the joint statistical modelling of species occurrences have gained a lot of popularity.

These models are sufficiently flexible to allow several fixed or random predictors to be incorporated in modelling species occurrences, allowing different aspects of habitat variability or species trait differences to be controlled against. In such a framework, one can infer species associations from their residual covariances, that is, the correlation between their residual probabilities of occurrence. One commonly used example of this approach are hierarchical models of species composition, as implemented in the HMSC R package (Ovaskainen & Abrego, 2020; Tikhonov et al., 2020). This is the method we will use in this article.

In this framework, we will use four different model specifications. First, a null model (Hmsc M0) that includes no latent factor or covariate, and therefore is equivalent to the above simplest methods (Fisher's tests or fixed-equiprobable permutation schemes). Second, differences across patches will be accounted for by including patch richness (i.e. the number of species in a patch) as a latent factor (Hmsc M1). Then, patch age (x) will be explicitly included as an environmental covariate in the model (Hmsc M2). As a variant of M2, patch richness will be also used as a proxy for patch age (Hmsc M2'). In other words, the difference between models Hmsc M1 and Hmsc M2' is that patch richness is a latent factor in the former, but an environmental covariate in the latter.

## 2.5 | Parameter values and numerical simulations

In the following, we will analyse how null co-occurrence matrices generated from our general metacommunity model behave, in terms of producing significant species associations, when subjected to the different methods described above (Figure 1). We will first derive mathematical predictions and test them with numerical simulations.

For all numerical simulations, we will use a synthetic community of 31 species, constructed to exhibit variation in their occupancies and extinction/colonization parameters (Figure 4a; see Appendix S2 for details). Although the parameterization does not correspond to a particular set of real species, the range of variation they present in colonization rates is similar to that found in natural metacommunities, for instance in Caribbean pond snails (Dubart et al., 2019) and grassland plant communities (Tilman, 1994).

We assumed the most classical metacommunity setup, that is, no external immigration ($m_i = 0$) and a constant rate of patch disturbance. The value of the disturbance rate was set to $\mu(x) = 0.1$ or $\mu(x) = 0.2$ per unit of time. The maximum life span of patches was set to $X = 20$ time units. These values correspond to low or intermediate levels of patch disturbance, considering the average species-specific extinction rate was 0.5 and the average colonization rate was 1.5. This means that extinctions caused by patch disturbance were always several times less frequent than other extinction/recolonization events. We chose these values in a conservative mindset, knowing that in metacommunities more dominated by patch disturbance, and/or with more variable species traits, the patterns we report could be more pronounced.
We varied the sample size (number of patches in the co-occurrence matrices) between 300 and 1,500, with intermediate values 500 and 1,000 (for recent studies using similar sample sizes, see Dubart et al., 2019; Façon et al., 2021; Opedal, Ovaskainen, et al., 2020). Finally, we also considered the possibility that not all species might be affected by patch disturbance: some species might be immune to disturbance ('immune species'). To address this possibility, we used the same set of 31 species, but rendered six of them, evenly spaced along trait values, invulnerable to patch disturbance (μ(x) = 0 for them), for each parameter combination.

This corresponded to a total of 2 × 4 × 2 = 16 parameter combinations. For each combination, we generated random co-occurrence matrices from the null metacommunity model (Equation 1). Each time, we applied the different methods introduced above to test for significant species associations, with standard settings recommended in the literature (Gotelli, 2000; Ovaskainen & Abrego, 2020). We computed, for each parameter combination and statistical method, the percent of spurious positive associations (fraction of species pairs declared positively associated) and the percent of spurious negative associations (fraction of species pairs declared negatively associated).

All simulations and figures reported can be reproduced using the accompanying R markdown (Appendix S2).

# 3 | RESULTS

## 3.1 | Mathematical predictions

In this section, we analyse mathematically the action that different methods have if applied on community matrices generated from our general model (1) and what conclusions would follow in terms of the species associations detected. We derive predictions under the limit of large community matrices (many patches and many species). Readers who are not mathematically inclined might decide to skip this section, and focus on the predictions we obtained, as illustrated in Figure 3c-f.

### 3.1.1 | Patch disturbance generates spurious positive associations with methods that do not account for patch age

#### Direct independence tests (e.g. Fisher's tests)

In the absence of any interaction, species are independently distributed among patches within any particular age class (Appendix S1, Section S1-3). This implies that the fraction of patches occupied by species i and j for all patch ages ($q_{i,j}$) will be determined by the steady-state probabilities of patches of age x co-occupied by species i ($p_{i|x}$) and species j ($p_{j|x}$):

$$q_{i,j,\cdot} = \int_{0}^{x} p_{i|x} p_{j|x} p_{\cdot,x} dx, \qquad (4)$$

which can be equivalently expressed as:

$$q_{i,j,\cdot} = p_{i,\cdot} p_{j,\cdot} \int_{0}^{x} p_{\cdot,x} \pi_{i/x} \pi_{j/x} dx = p_{i,\cdot} p_{j,\cdot} (1 + \text{Cov}(\pi_{i/x}, \pi_{j/x})), \quad (5)$$

given that the mean of the profile is one (Appendix S1, Section S1-1.5). Since $\pi_{i/x}$ and $\pi_{j/x}$ are increasing functions of x (Equation 2), from Harris' inequality, we have $\text{Cov}(\pi_{i/x}, \pi_{j/x}) \ge 0$ (Appendix S1, Section S1-1.5). From (5), we therefore have

$$q_{i,j,\cdot} \ge p_{i,\cdot} p_{j,\cdot} \qquad (6)$$

Equality occurs only if at least one of the two species is distributed uniformly over all patch ages, that is, $\pi_{k/x} = 1$ for all x so that $\text{Cov}(\pi_{i/x}, \pi_{j/x}) = 0$. From Equation 2, this is not the case in general. This only occurs if there is no patch disturbance at all ($\mu_x = 0$ and $X \to \infty$), in which case all patches are effectively infinitely old so that $\pi_{i/x}$ and $\pi_{j/x}$ both tend to one, and $\text{Cov}(\pi_{i/x}, \pi_{j/x})$ tends to zero. But in all other cases, species are not distributed independently: they are positively associated, that is, they co-occur more often than expected from the 'null' product of their respective overall occupancies (see Equation 6).

As a consequence, methods based on the test of pairwise independence of species, such as direct Fisher's tests or more recently proposed methods (e.g. Veech, 2013) would consistently detect spurious positive associations (aggregation) among species, even in the absence of any interaction between species. Specifically, positive associations will be strongest, and thus more likely to be statistically detected, for species pairs with a large covariance of their relative profiles (Equation 5), that is, for pairs of species that are both slow. In contrast, species pairs including at least one fast species are more likely to be declared independent.

#### Fixed-equiprobable permutations (Sim2)

The same conclusion holds for more elaborate methods that do not control for patch age, such as the 'fixed-equiprobable' permutation algorithm, also known as 'Sim2' (Gotelli, 2000). To show this, let us consider the commonly used partial C-score $C_{ij}$ as our metric for species association (Equation 3). We will denote with an asterisk the value of quantities after the permutation algorithm has been applied to a co-occurrence matrix.

A fixed-equiprobable matrix permutation algorithm preserves the number of occurrences per row (i.e. overall occupancies of species) but reshuffles species across all patches equiprobably. Therefore, after permutations, species overall occupancies are unchanged ($p_{i,\cdot}^* = p_{i,\cdot}$ for all i), but all species become equiprobable over the different age classes, on average. Therefore, their relative distribution profiles become flat: $\pi_{k/x}^* = 1$ for all x. From Equation 5, it follows

$$q_{i,j,\cdot}^* = p_{i,\cdot}^* p_{j,\cdot}^* = p_{i,\cdot} p_{j,\cdot} \qquad (7)$$

This indicates that the algorithm effectively breaks the positive association generated by patch disturbance. From Equation 6, this implies $q_{i,j,\cdot}^* \le q_{i,j,\cdot}$. Permutated matrices will on average have a deficit of co-occurrences, for all species pairs, compared to the original matrix.
<figure><img src="image_6.png" /><figcaption>FIGURE 3 Mathematical predictions. (a) The set of synthetic species. Occupancy as a function of patch age x ($p_{ijx}$) is shown for each of the 31 species. Species traits were drawn randomly to have variable overall occupancies and fastnesses. (b) Pairwise covariances of relative profiles for all species pairs. Species are ranked according to fastness: species 1 is the fastest (lowest variance), and species 31 is the slowest (highest variance). The variance of the average profile in the community ($\hat{\pi}_x$) was subtracted from covariances to help distinguish lower/higher than average values (colour bar). See Equations (4) and (7). (c) Mathematical predictions associated with methods testing for pairwise species independence (Fisher's test, fixed-equiprobable -Sim2-permutations and the simplest joint species distribution model Hmsc M0). (d-e) Same as (c) based on methods that control for patch richness, that is, the fixed-fixed permutation (Sim9) algorithm (d), and joint species distribution model Hmsc M1 (e). (f) Same as (c) based on joint species distribution model Hmsc M2, a method that controls for patch age. Species are ordered as in (b). The bullets indicate the predictions relative to an example species (species 27) being immune to patch disturbance</figcaption></figure>

Since $p_{i,}$ and $p_{j,}$ are unaffected by permutation, we can directly conclude that $C_{ij}^* > C_{ij}$ as soon as patch disturbances occur. Using the fixed-equiprobable permutation algorithm will therefore consistently yield spurious facilitation signals between virtually all pairs of species, exactly as direct tests for species independence do (Figure 3c).

### Joint species distribution modelling (Hmsc M0)

Joint species distribution models that do not account for patch differences (such as our Hmsc M0), by assuming species to be

independently distributed across patches, will suffer from the same caveat.

Hmsc models predict the probability of presence of some species i in a patch of age x, $p_{i|x}^*$, where the asterisk stands for 'predicted value'.

The model predicts the probability that species i is present in a patch regardless of the patch age, that is,

$$p_{i|x}^* = p_{i,}.$$
It follows that the residual association (covariance) of two species *i*, *j* is

$$\int_{0}^{x} p_{*,x}(p_{i|x} - p_{i,*})(p_{j|x} - p_{j,*})dx$$

which is equal to, dividing by $p_{i,*}p_{j,*}$,

$$p_{i,*}p_{j,*} \int_{0}^{x} p_{*,x}(\pi_{i/x} - 1)(\pi_{j/x} - 1)dx. \quad (8)$$

Since all species have relative distribution profiles smaller than one for small *x*, and greater than one for large *x*, it follows that the integrand is positive for most *x* values, and thus the two species will have positive residual association.

Note that if at least one of the two species is infinitely fast, then the corresponding parenthesis is null for all *x*: the two species will have no residual association (i.e. will be independently distributed).

Hmsc M0 models will thus report positive residual associations for all species pairs. The ability to detect those associations will however depend on statistical power, especially for pairs of fast species, as in previous methods (independence tests and fixed-equiprobable permutations).

### 3.1.2 | Patch disturbance yields a characteristic signature of spurious associations even with methods that control for patch richness

We now turn to methods that can control for systematic differences among patches, in this case emerging from heterogeneity in the time since the last disturbance occurred (patch age). We here distinguish two groups of approaches. We first treat the popular fixed-fixed permutation algorithm, also known as Sim9 or swap (Gotelli, 2000). We then consider joint species distribution models that include patch richness as a factor (as our Hmsc M1).

#### Fixed-fixed permutation approaches (Sim9)

By fixing the row totals, this method effectively maintains the number of species per patch unchanged, and can thus deal with the positive association of species that causes issues with the previous set of approaches.

Similar to the permutation algorithm discussed in the previous section (Sim2), preservation of row totals implies $p_{i,*}^* = p_{i,*}$. Therefore, we only need to study changes in $q_{i,j,*}$ to understand changes in C-scores. Preservation of column totals further means that the number of occurrences will remain the same in every age class, which can be stated as $\sum_{k=1}^{s} p_{k,x}^* = \sum_{k=1}^{s} p_{k,x}$ for all *x*. Species are otherwise reshuffled indifferently across all patch ages. Let $w_i$ be the relative occupancy of species *i* in the matrix, that remains unchanged after permutations:

$$w_i = \frac{p_{i,*}}{\sum_{k=1}^{s} p_{k,*}}.$$

A fixed-fixed permutation algorithm yields:

$$p_{i,x}^* = \left( \sum_{k=1}^{s} p_{k,x} \right) w_i,$$

leading to, for all $i = 1, 2, \dots, s$,

$$\pi_{i/x}^* = \frac{p_{i,x}^*}{p_{*,x} p_{i,*}} = \frac{\sum_{k=1}^{s} p_{k,x}}{\sum_{k=1}^{s} p_{k,*}} = \frac{\sum_{k=1}^{s} p_{k,*} \pi_{k/x}}{\sum_{k=1}^{s} p_{k,*}} = \sum_{k=1}^{s} w_k \pi_{k/x} = \hat{\pi}_{x}.$$

This indicates that species relative distribution profiles are replaced with their weighted-average over all species, with weights equal to species relative overall occupancies. Let us denote this weighted-average by $\hat{\pi}_x$. We note that $\hat{\pi}_x$ can also be expressed as

$$\hat{\pi}_x = \frac{\sum_{k=1}^{s} p_{k|x}}{\sum_{k=1}^{s} p_{k,*}}, \quad (9)$$

that is, $\hat{\pi}_x$ is the mean species richness in patches age of *x* relative to the mean richness over all patch ages.

We thus obtain:

$$q_{i,j,*}^* = p_{i,*}p_{j,*} \int_{0}^{x} p_{*,x} \pi_{i/x}^* \pi_{j/x}^* dx = p_{i,*}p_{j,*} \int_{0}^{x} p_{*,x} \hat{\pi}_x^2 dx = p_{i,*}p_{j,*}(1 + \text{Var}(\hat{\pi})). \quad (10)$$

Comparing with Equation 5, we see that the term correcting for non-independence, $\text{Cov}(\pi_{i/x}, \pi_{j/x})$, is here replaced with $\text{Var}(\hat{\pi})$. The difference $q_{i,j,*} - q_{i,j,*}^*$ has the same sign as $\text{Cov}(\pi_{i/x}, \pi_{j/x}) - \text{Var}(\hat{\pi})$.

The permutation algorithm is thus biased, except if $\text{Cov}(\pi_{i/x}, \pi_{j/x}) = \text{Var}(\hat{\pi})$ for all species pairs. This would occur if all species had identical relative distribution profiles, that is, when all species are 'similar', in the sense $p_{i,x} = \kappa p_{j,x}$ for all *x*, with some positive $\kappa$. This trivially occurs if there is no patch disturbance at all, or if all species have identical extinction and colonization parameters. Beyond that, species may differ in their parameters and overall occupancies and still be similar in that sense, but this requires quite constrained parametric conditions (Appendix S1, Section S1-2.3).

Except for these restricted scenarios, the permutation algorithm will yield spurious associations in both directions. Species pairs whose relative distribution profiles have covariance smaller than the variance of the average community profile ($\text{Cov}(\pi_{i/x}, \pi_{j/x}) < \text{Var}(\hat{\pi})$) will appear to have a spurious negative association. Conversely, species pairs for which $\text{Cov}(\pi_{i/x}, \pi_{j/x}) > \text{Var}(\hat{\pi})$ will appear to have a spurious positive association. Only species pairs for which $\text{Cov}(\pi_{i/x}, \pi_{j/x}) = \text{Var}(\hat{\pi})$ would appear as independently distributed.

In practical terms, species pairs containing at least one species that is 'faster' than average (i.e. whose relative profile is flatter) will have low relative covariance, and thus will be declared negatively associated. Indeed, if the fast species (say *i*) is fast enough (regardless of how slow the other is), then $\pi_{i/x} \approx 1$ for all *x* thus and $\text{Cov}(\pi_{i/x}, \pi_{j/x}) \approx 0 = \text{Var}(\pi_{i/x})$, which is less than the variance of the average. In contrast, species pairs that are both 'slower' than average (i.e. whose relative profiles are steeper) will have large relative
covariance, and will be declared positively associated. Pairs of species that are both close to average fastness may yield no particular signal.

If species are ordered according to fastness (i.e. variance), this produces a characteristic signature of spurious association signals. The pairwise covariance of all species pairs used in reference set of species is shown in Figure 3b. The corresponding mathematical predictions for fixed-fixed permutation algorithms are illustrated in Figure 3d (note the vertical and horizontal blue lines for species 27 are because that species is assumed to be immune to patch disturbance as an illustration of the potential effect of that phenomenon; see below).

### Joint species distribution modelling (Hmsc M1)

In the context of joint species distribution modelling, it has been suggested that an equivalent to the fixed-fixed permutation scheme discussed above would be to include patch richness as a latent factor in the model of species occurrence (our Hmsc M1). This was hypothesised in Ovaskainen and Abrego (2020, p. 335; see Appendix S2 for details).

We can similarly formulate mathematical predictions regarding the behaviour such an approach would have under our null metacommunity model. In the limit of large co-occurrence matrices, the number of species in a patch ($\sum_k p_{k|x}$) tends to become tightly associated with patch age $x$. In other words, one may neglect the variance in patch age for a given species richness. A statistical model that describes the probability of occurrence, while including patch richness as a latent factor, would therefore predict species $i$ to occur in some particular patch of age $x$ with probability $p_{i|x}^* = p_{i,\hat{\pi}_x}$, since $\hat{\pi}_x$ is the relative richness of a patch of age $x$ (Equation 9).

The expected residual for this species in this sort of patch would be $p_{i,}(\pi_{i,x} - \hat{\pi}_x)$ and the residual covariance between two species $i$ and $j$, over all patches, would be, from Equation (8),

$$p_{i,}p_{j,} = \int_{0}^{x} p_{*,}(\pi_{i,x} - \hat{\pi}_{x})(\pi_{j,x} - \hat{\pi}_{x})dx. \quad (11)$$

From the above equation, we can immediately conclude that if at least one of the two species has average fastness, in the sense defined above, then $\pi_{i,x} \approx \hat{\pi}_x$ for all $x$, yielding a residual association of zero (independence). Otherwise, if the two species are faster (slower) than average, then both parentheses in the integral will be first positive (negative), and then negative (positive), as $x$ increases. It follows that the overall residual association will be found positive. By the same argument, if one species is faster than average and the other is slower, the two parentheses will have opposite signs, and the overall residual association will be negative.

As for the above permutation scheme, we expect a mixture of positive and negative associations, unless all species have the same fastness. The average fastness again plays a critical role in determining where the switch from negative to positive associations, or vice versa, occurs. However, the resulting pattern of species

associations is markedly different. We predict a block-like pattern, with faster species positively associated together, slower species also positively associated, and the two groups of species mutually negatively associated. We thus predict a four-block pattern, reminiscent of differential specialization on two habitats, as illustrated in Figure 3e.

In any case, controlling for differences across patches is therefore insufficient to handle the metacommunity dynamical consequences of patch disturbance, even though it brings some improvement. First, correct results are obtained if all species are 'similar' (i.e. have identical fastness). Second, the overall number of spurious associations obtained will likely be reduced, since some species pairs will likely be correctly declared as independently distributed. However, this comes at the price of generating spurious associations of both kinds (positive and negative), not only positive associations.

### 3.1.3 | Patch disturbance requires species-specific modelling of patch age

Efficiently dealing with patch disturbance in co-occurrence analyses therefore requires both patch age and species traits to be taken into account, to model appropriately the differential effect that patch age has on species with different fastnesses. In other words, there must be an interaction between the effect of patch age and the identity of species.

This cannot be done with direct independent tests or simple null model approaches, but can be in the context of joint species distribution models, such as our model Hmsc M2.

The model predicts the probability that species $i$ is present in a patch as an increasing function of patch age, that is, it will yield $p_{i|x}^* \approx p_{i|x}$ where the approximation denotes the fact that the functional form of the regression curve may not exactly match the shape of the $p_{i|x}$ function.

It follows that residuals will be close to zero, and so will be the residual associations, for all species pairs. This is illustrated in Figure 3f.

### 3.1.4 | What if some species are immune to patch disturbance?

So far, we assumed all species in the metacommunity were affected by patch disturbances. In practice, in large sets of species, it may be that only a subset of species are susceptible to disturbance events, while some others are not ('immune' species). Based on our mathematical results, this possibility is straightforward to deal with: the immune species would behave as a species with infinite fastness ($\text{Var}(\pi_{i/x}) = 0$), regardless of its actual colonization/extinction rates (Appendix S1, Section S1-4).

As a consequence, any species pair including at least one immune species is expected to return no overall association, if using pairwise independence tests or equivalent methods treated in Section 2.1;
see Figure 3c). Therefore, fewer spurious associations would be de-
tected with such methods.

On the contrary, by increasing the heterogeneity of fastness
among species in the community, the presence of immune spe-
cies would degrade the performance of methods that control for
patch richness. More specifically, when using a fixed-fixed per-
mutation method, any species pair including one immune species
will systematically appear as negatively associated with the oth-
ers, since Cov(πᵢ/ₓ, πⱼ/ₓ) = 0. For the same reason, if using a joint
species distribution model with patch richness as a latent fac-
tor, immune species will tend to cluster with the fastest species,
and to be declared negatively associated with the slower species
(Figure 3d,e).

All the mathematical predictions derived in this section are sum-
marized in Figure 3. We will now proceed to test them with numer-
ical simulations.

3.2 | Numerical simulations

3.2.1 | Patch disturbance generates spurious positive associations with methods that do not account for patch age

Consistent with our mathematical predictions, species were posi-
tively associated overall in the simulated co-occurrence matrices,
even under our 'null' model (Equation 1). Using the reference pa-
rameter set (μ = 0.1, N = 1,000 and no immune species), the odds
of species co-occurrence were on average 75% higher than ex-
pected under independence (Table 1). About 10% of all species
pairs had odds of co-occurrence more than 2.5 times as large as
expected.

Results for methods not controlling for patch richness or patch
age (pairwise Fisher's tests, fixed-equiprobable permutations and
species distribution model Hmsc M0) are shown in Figure 4. As
predicted, all methods detected significant positive associations
for the vast majority of species pairs, except for some of the spe-
cies pairs involving the fastest species. Direct pairwise Fisher's
tests and fixed-fixed permutations yielded very similar results,
with about 75% of species pairs declared as positively associated
(Figure 4a,c). Hmsc M0, which uses a quite different statistical ap-
proach, happened to be more sensitive, and thus more error prone,
declaring more than 90% of species pairs positively associated
(Figure 4e).

Again confirming predictions, the inclusion of species immune
to patch disturbance reduced the number of positive associations
detected, as species pairs including at least one immune species
were reported as independently distributed, with a few exceptions
(Figure 4b,d,f). The main difference between the three methods
was that Hmsc M0 could occasionally go as far as declaring one im-
mune species as negatively associated with the non-immune species
(Figure 4f).

3.2.2 | Patch disturbance yields a characteristic signature of spurious associations even with methods that control for patch richness

Simulation results applying fixed-fixed permutation schemes (Sim9) and the species distribution model Hmsc M1 are presented in Figure 5. As expected from our mathematical analyses, we observe the detection of spurious associations of both kinds (positive and negative). The patterns of species associations closely resemble those of mathematical predictions, be it for fixed-fixed permutations (see Figure 3d) or Hmsc M0 (see Figure 3e). The main difference from predictions is the absence of some predicted associations, due to lack of statistical power. This is particularly the case for species pairs with covariance close to average, for a fixed-fixed permutation algorithm (Figure 5a).

The total number of species pairs for which a spurious associa-
tion is reported is slightly smaller than it was without controlling for
patch richness. It is 56% with a fixed-fixed permutation scheme ver-
sus 73% with a fixed-equiprobable scheme. Similarly, it is 12% with
Hmsc M1 versus more than 90% with Hmsc M0.

However, both positive and negative spurious associations are
detected. There is a slight majority of negative associations with
fixed-fixed permutations, but, in contrast, a majority of positive as-
sociations with Hmsc M1. Unlike what was observed for the previous
set of methods, Hmsc modelling is here less sensitive, and thus less
error prone, than the comparable permutation algorithm. Therefore,
fixed-fixed permutations and Hmsc M1 models differ not only in the
qualitative type of characteristic signature they produce, but also in
their quantitative performance (power and ratio of negative versus
positive associations).

Introducing immune species also yielded results conforming to
mathematical predictions (Figure 5b,d). Interestingly, the presence
of immune species degraded the overall performance of the two
methods, but not to the same extent: the percentage of spurious as-
sociations detected increased from 12% to 28% in the case of Hmsc
M1 and, more modestly, from 57% to 59% in the case of fixed-fixed
permutations.

3.2.3 | Patch disturbance requires species-specific modelling of patch age

Species distribution models can describe how the probability of
occurrence increases, for each species, as a function of patch age.
This can be done rather straightforwardly by providing patch age
as an environmental covariate, over which species occupancies are
regressed, as in our model Hmsc M2.

As predicted, such an approach can restore satisfying re-
sults, that is, no association was detected. Examples are shown in
Figure 6a,b. Using the same set of species and co-occurrence matrix
as before, the use of model Hmsc M2 suppressed almost all spurious
association signals, with or without immune species, as desired. The
<table><caption>TABLE 1 Performance (% of errors) for each of the five methods tested (Fisher, Fixed-fixed, Hmsc MO, Hmsc M1, Hmsc M2), under the 16 parametric scenarios (sample size N × disturbance rate μ × No Immune species/With immune species). Note that the 'fixed-fixed' algorithms is also known as 'Sim9'. Grey shading corresponds to the default parameterisation used throughout the paper</caption><thead><tr><th rowspan="2">Scenario</th><th colspan="2">Fisher</th><th colspan="2">Fixed-fixed</th><th colspan="2">Hmsc MO</th><th colspan="2">Hmsc M1</th><th colspan="2">Hmsc M2</th></tr><tr><th>%+</th><th>%-</th><th>%tot</th><th>%+</th><th>%-</th><th>%tot</th><th>%+</th><th>%-</th><th>%tot</th><th>%+</th><th>%-</th><th>%tot</th></tr></thead><tbody><tr><td colspan="13"><strong>No immune species</strong></td></tr><tr><td>μ = 0.2</td><td>N = 300</td><td>38.3</td><td>0</td><td>38.3</td><td>11.8</td><td>8.8</td><td>20.6</td><td>75</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td></td><td>N = 500</td><td>58.3</td><td>0</td><td>58.3</td><td>20.2</td><td>19.8</td><td>40</td><td>78.3</td><td>0</td><td>4.5</td><td>0</td><td>0.4</td></tr><tr><td><sup>a</sup></td><td>N = 1,000</td><td>73.8</td><td>0</td><td>73.8</td><td>29.7</td><td>27.3</td><td>57</td><td>98.9</td><td>0</td><td>12.1</td><td>0.9</td><td>2.8</td><td>3.7</td></tr><tr><td>μ = 0.1</td><td>N = 1,500</td><td>79.6</td><td>0</td><td>79.6</td><td>35</td><td>28.6</td><td>63.6</td><td>98.5</td><td>0</td><td>13.8</td><td>20.6</td><td>1.3</td><td>1.9</td><td>3.4</td></tr><tr><td></td><td>N = 300</td><td>19.6</td><td>0</td><td>19.6</td><td>5</td><td>3.6</td><td>8.6</td><td>69.9</td><td>0</td><td>0.2</td><td>0.2</td><td>0</td><td>0</td><td>0</td></tr><tr><td></td><td>N = 500</td><td>45.2</td><td>0</td><td>45.2</td><td>11.8</td><td>10.5</td><td>22.4</td><td>70.1</td><td>0</td><td>1.1</td><td>1.5</td><td>0</td><td>0.6</td><td>0.6</td></tr><tr><td></td><td>N = 1,000</td><td>69.4</td><td>0</td><td>69.4</td><td>19.1</td><td>18.5</td><td>37.6</td><td>98.5</td><td>0</td><td>2.5</td><td>3.4</td><td>0</td><td>0</td><td>0</td></tr><tr><td></td><td>N = 1,500</td><td>75</td><td>0</td><td>75</td><td>25.6</td><td>22.6</td><td>48.2</td><td>87.3</td><td>0</td><td>6.4</td><td>8.4</td><td>0</td><td>0.2</td><td>0.2</td></tr><tr><td colspan="13"><strong>With immune species</strong></td></tr><tr><td>μ = 0.2</td><td>N = 300</td><td>22.6</td><td>0</td><td>22.6</td><td>4.7</td><td>6</td><td>10.7</td><td>50</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td></td><td>N = 500</td><td>37.8</td><td>0</td><td>37.8</td><td>20.4</td><td>17.8</td><td>38.3</td><td>45.2</td><td>0</td><td>7.3</td><td>8.4</td><td>0</td><td>0</td><td>0</td></tr><tr><td></td><td>N = 1,000</td><td>48.4</td><td>0</td><td>48.4</td><td>35.3</td><td>23.4</td><td>58.7</td><td>64.1</td><td>4.7</td><td>18.5</td><td>28.4</td><td>0</td><td>0.4</td><td>0.4</td></tr><tr><td></td><td>N = 1,500</td><td>51.8</td><td>0</td><td>51.8</td><td>41.3</td><td>27</td><td>68.4</td><td>59.3</td><td>0</td><td>20</td><td>25.6</td><td>0</td><td>0</td><td>0</td></tr><tr><td>μ = 0.1</td><td>N = 300</td><td>7.5</td><td>0</td><td>7.5</td><td>1.9</td><td>3.6</td><td>5.6</td><td>45.2</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td></td><td>N = 500</td><td>28.2</td><td>0</td><td>28.2</td><td>10.1</td><td>10.3</td><td>20.4</td><td>45.2</td><td>0</td><td>3</td><td>5.8</td><td>8.8</td><td>0</td><td>0</td></tr><tr><td></td><td>N = 1,000</td><td>44.7</td><td>0</td><td>44.7</td><td>20.4</td><td>19.8</td><td>40.2</td><td>64.5</td><td>3.9</td><td>3.2</td><td>7.3</td><td>10.5</td><td>0</td><td>0</td></tr><tr><td></td><td>N = 1,500</td><td>49.2</td><td>0</td><td>49.2</td><td>27.7</td><td>24.9</td><td>52.7</td><td>64.5</td><td>3.9</td><td>3.2</td><td>7.3</td><td>10.5</td><td>0</td><td>0</td></tr><tr><td>Grand average</td><td></td><td>46.8</td><td>0</td><td>46.8</td><td>20</td><td>17</td><td>37</td><td>69.6</td><td>0.78</td><td>2.5</td><td>6.3</td><td>8.9</td><td>0.17</td><td>0.23</td><td>0.33</td></tr></tbody></table>

<footnote>a This is the scenario used in all figures.</footnote>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 4 Numerical simulations: methods testing for pairwise species independence. (a) Species pairs significantly associated according to Fisher's pairwise tests. (b) Same as (a), with a subset of species immune to patch disturbance (shown as bullets). (c) Species pairs significantly associated according to the fixed-equiprobable (Sim2) permutation algorithm. (d) Same as (c), with a subset of immune species. (e) Species pairs significantly associated according to the simplest joint species distribution model (Hmsc M0). (f) Same as (e), with a subset of immune species. Species are ordered as in Figure 3b. The sample co-occurrence matrix was generated from Equation (1) with N = 1,000 sites (patches) and the illustrative parameter set</figcaption></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

percentage of spurious associations detected was not strictly speaking zero, but remained under a tolerable error rate of 5%.

By regressing the probability of occurrence on patch age, for every species individually, such a model appropriately captures both the overall increase in occupancy with patch age (Figure 6c) and the differential response of the different species, depending on their fastness (Figure 6d).

Note however that this approach requires additional data, namely the age of patches, on top of the species co-occurrence matrix. This would not always be possible. Considering the fair association between patch age and patch richness (Figure 6c), one could consider

using patch richness as a surrogate for patch age. As patch richness is already encoded in the species co-occurrence matrix, this approach, like all previous sets of methods, would not require additional data.

However, if using patch richness as a proxy of patch age (model Hmsc M2'), the method fails to provide correct results (see Figure S1). The association presented in Figure 6c appears not to be strong enough for the model to appropriately capture age-related patch differences. Arguably, in datasets containing even more species (e.g. 50-100), the correlation might well be stronger, but in many situations, patch richness would also vary because of several factors beyond patch age. Therefore, its capacity to serve as a proxy of the
<figure><img src="image_13.png" /><figcaption>**FIGURE 5** Numerical simulations: methods controlling for patch richness. (a) Species pairs significantly associated according to the fixed-fixed (Sim9) permutation algorithm. (b) Same as (a), with a subset of species immune to patch disturbance (shown as bullets). (c) Species pairs significantly associated according to a joint species distribution model (Hmsc M1) that models patch richness as a single latent factor. (d) Same as (c), with a subset of immune species. Same species and sample co-occurrence matrix as in previous figures</figcaption></figure>

latter would be compromised anyway. Consequently, using explicit independent data on patch age is the only generally recommendable strategy.

### 3.2.4 | Effect of changing parameter values and sample size

The performances of the main methods, under the different parameter scenarios, are summarized in Table 1. Confirming the previous results, only Hmsc M2 yielded satisfying results in all cases. All other methods generally detected spurious associations, often at quite high rates, sometimes exceeding 70%. Obviously, reducing sample size down to *N* = 300 reduced the general statistical power, and therefore the rate of spurious detections. In the case of species distribution model Hmsc M1, this sometimes made the percentage of errors fall to acceptable levels. However, this is just indicative of a lack of power, and would also come at the expense of higher type II errors in actual datasets containing genuine species associations.

Decreasing the rate of patch disturbance *μ*, as expected, also decreases the overall amount of spurious associations detected (Table 1). This is simply because this weakens the metacommunity-dynamics-driven patterns, in the absence of which all methods compared here perform correctly.

Interestingly, the effect of introducing some species immune to patch disturbance systematically reduces the rate of spurious detections for methods that do not control for patch richness (Fisher's tests or Hmsc M0) but can have variable consequences for methods that do (fixed-fixed permutations and Hmsc M1). The consequences depend strongly on sample size. At low sample sizes, the main consequence is to reduce the average rate of perturbation, and thus to decrease the number of spurious detections, as statistical power is limiting. At large sizes, on the contrary, the presence of immune species increases the number of spurious detections, as it increases the heterogeneity in fastness among species, to which these methods are sensitive.

Our main overall conclusions are synthesized in Table 2.

## 4 | DISCUSSION

Metacommunity dynamics matter. We demonstrated here how simple metacommunity dynamics can cause the detection of spurious species associations, even in null communities with no interactions. This occurs as soon as patch disturbances exist, that is, catastrophic events that cause all species (or at least a subset of them) to be locally wiped out. The latter is an arguably common element in natural and exploited ecosystems (Pickett & White, 2013), and one often,
<figure><img src="image_14.png" /><figcaption>FIGURE 6 Numerical simulations: modelling the effect of patch age. (a) A joint species distribution model that incorporates patch age as an explicit environmental covariate (Hmsc M2) yields correct results, even with patch disturbance. The few spurious species associations detected are within the type I error rate. Note that log patch age was used as the covariate, since patch age is distributed exponentially. (b) Same as (a), with a subset of species immune to patch disturbance (shown as bullets). (c) Predicted species richness as a function of patch age. Grey dots represent the actual data in the sample co-occurrence matrix. (d) Predicted probability of occurrence as a function of patch age, for the fastest species and the slowest species. Same species and sample co-occurrence matrix as in previous figures</figcaption></figure>

<table><caption>TABLE 2 Are detection methods robust to patch disturbance? Summary of main findings. For each type of approach, we only mention methods tested in this study, but conclusions apply to other similar methods as well. Note that 'Sim2' and 'Sim9' denote the 'fixed-equiprobable' and 'fixed-fixed' algorithms, respectively</caption><thead><tr><td>Type of approach</td><td>Example methods</td><td>Robust?</td><td>Type of bias</td></tr></thead><tbody><tr><td>Direct independence tests</td><td>Fisher's tests; Sim2; Hmsc M0</td><td>No<sup>a</sup></td><td>Positive</td></tr><tr><td>Controlling for patch richness</td><td>Sim9; Hmsc M1; Hmsc M2'</td><td>No<sup>b</sup></td><td>Positive/negative</td></tr><tr><td>Modelling effect of patch age</td><td>Hmsc M2</td><td>Yes</td><td>None</td></tr></tbody></table>

though not always, included in metapopulation and metacommunity models (Leibold et al., 2004).

It is classically assumed that in the absence of any ecological interactions, habitat heterogeneity or species differences in dispersal patterns, species would be distributed independently across patches. This ‘null’ assumption underlies many interpretations of species co-occurrence patterns. Our analysis shows that this intuition is true for any given patch age (time since last disturbance), but it is not true overall. Patch age acts as a hidden structuring variable, the implications of which for co-occurrence patterns have been overlooked so far. Indeed, among the list of processes typically invoked as drivers which might cause departures from species independence, patch age does not normally appear (Gotelli, 2000; Ovaskainen & Abrego, 2020).

This result may seem counter-intuitive at first glance. Indeed, if species are initially independently distributed in a metacommunity, patch disturbances do not per se cause any association to appear: as they occur regardless of species composition, they preserve species independence, unlike interaction-driven extinctions would. Rather, it is the outcome of extinction/recolonization dynamics (Equation 2) that inevitably causes, at steady state, species to be more frequent in older patches. Since all species are more likely to occupy older patches, all appear positively associated overall. We found that this signal was detectable even with realistic sample sizes (number of patches), under relatively low levels of patch disturbance, and for variations among species fastness similar to what has been observed in actual metacommunity systems. Obviously, the amount

<footnote><sup>a</sup>Works only if no patch disturbance.</footnote>

<footnote><sup>b</sup>Works only if no patch disturbance, or all species similar (same fastness).</footnote>
of spurious species associations detected generally declined with sample size (Table 1), but so would the power to detect genuine species association signals. Moreover, the range of sample sizes we used may quickly become conservative, as more and more ecological studies now have datasets with several hundreds or thousands of patches sampled (Dubart et al., 2019; Facon et al., 2021; Opedal, Ovaskainen, et al., 2020). The development of metagenomics and microbiome studies will undoubtedly contribute to this upward trend in sample sizes, and so the phenomenon reported here could become more and more prevalent.

The consequences of patch disturbance cannot simply be dealt with by controlling for patch differences in richness. If doing so, spurious associations, both negative and positive, also appear, as soon as species have different relative distribution profiles (fastnesses). These spurious associations have a different cause: they stem from the fact that slow and fast species are concentrated in different parts of the patch-age spectrum, which can be thought of as a form of niche differentiation. This differential distribution with reference to patch age causes spurious negative association signals to appear. But the signature we report is more subtle than that (see Figure 5). The fastest species (or species immune to patch disturbance, if any) also appear to be segregated from each other, as well as from all the other species. Reciprocally, the slowest species, being more concentrated than others in the older parts of the patch age spectrum, appear positively associated (and segregated from the others).

The generation of spurious associations (Type I error) is in itself problematic. But beyond that, in more realistic communities with actual species interactions, these spurious patterns would interfere with and confound the actual signals. For instance, actual signals of species interactions, for example, an increase in species extinction rates with species number, that is, a 'neutral' form of competition (Bell, 2001; Hastings, 1987), could be hidden in the presence of patch disturbance: the negative association patterns would be cancelled, at least among pairs of slow species, by the positive signal caused by metacommunity dynamics alone. In such situations, false negatives (loss of power) will result. Furthermore, the power of detection is not homogeneously distributed across all species pairs: signals of competition would be weakened for slow species pairs, but reinforced for fast species pairs. Therefore, a diffuse type of competitive interaction with no discernible pattern (all species competitively identical) could turn into a structured pattern, with nonzero nestedness.

A possible solution to these issues consists of using species distribution models that explicitly incorporate patch age as a covariate (Figure 6). However, this requires having independent data on patch ages, which is far from being universal in co-occurrence analyses. While feasible for some types of systems (Ovaskainen et al., 2017), typically parasite or microbiome communities for which patches are hosts of known age, in many cases such data would be difficult or impossible to obtain. It may sometimes be possible to approximate patch age from measurable environmental variables, when patch disturbance causes environmental modifications and patches follow a succession-like cycle after them. Depending on the type of biological system considered, indicator species known to reflect

successional stage, or the dosage of chemical compounds indicative of disturbance (e.g. fire), may correlate with patch age. In the case of parasites, host size or other host attributes may also be used to infer host age. However, this approach would not feasible for all systems, particularly since a sufficiently strong correlation is required. Indeed, in our examples, using patch richness as a proxy of patch age failed to remove the spurious associations, even though the two variables were reasonably well correlated (Figure 6). In general, effectively controlling for patch age would almost certainly require genuine temporal data on patch life cycles (Ovaskainen et al., 2017, Figure 5c).

In general, our work here underscores once more how inferring species interactions from snapshots of co-occurrence patterns can be as hazardous as it is appealing, and requires careful thought about null expectations (Blanchet et al., 2020; Molina & Stone, 2020; Münkemüller et al., 2020). We join the recent burst of cautionary articles in highlighting these concerns, and in any case, we encourage investigators to evaluate whether detected species association patterns could match those reported here, and so could stem from pure metacommunity dynamical effects.

## ACKNOWLEDGEMENTS

The authors thank Maxime Dubart for guidance on using Hmsc and for sharing data on Caribbean snail metacommunity parameter estimates, and the anonymous peer reviewers for helpful suggestions. Colin Guéteme and Héloïse Villesèche participated in preliminary research on the topic.

## CONFLICT OF INTEREST

The authors declare no conflict of interest.

## AUTHORS' CONTRIBUTIONS

Authors contributed similarly to the manuscript preparation; V.C. and F.M.H. contributed most to the mathematical derivations; N.J.C. contributed most to the stochastic simulations and R coding.

## DATA AVAILABILITY STATEMENT

All the codes used to produce the results in this article are available via the doi https://doi.org/10.5281/zenodo.6400509 (Cunniffe, 2002).

## ORCID

Vincent Calcagno https://orcid.org/0000-0002-5781-967X

Nik J. Cunniffe https://orcid.org/0000-0002-3533-8672

Frédéric M. Hamelin https://orcid.org/0000-0003-2653-699X

## REFERENCES

Barberán, A., Bates, S. T., Casamayor, E. O., & Fierer, N. (2012). Using network analysis to explore co-occurrence patterns in soil microbial communities. *The ISME Journal*, 6(2), 343–351.

Barner, A. K., Coblentz, K. E., Hacker, S. D., & Menge, B. A. (2018). Fundamental contradictions among observational and experimental estimates of non-trophic species interactions. *Ecology*, 99(3), 557–566.
Begon, M., & Townsend, C. R. (2020). *Ecology: From individuals to ecosystems*. John Wiley & Sons.

Bell, G. (2001). Neutral macroecology. *Science*, 293(5539), 2413-2418.

Blanchet, F. G., Cazelles, K., & Gravel, D. (2020). Co-occurrence is not evidence of ecological interactions. *Ecology Letters*, 23(7), 1050-1063.
https://doi.org/10.1111/ele.13525

Blois, J. L., Gotelli, N. J., Behrensmeyer, A. K., Faith, J. T., Lyons, S. K., Williams, J. W., Amatangelo, K. L., Bercovici, A., Du, A., Eronen, J. T., Graves, G. R., Jud, N., Labandeira, C., Looy, C. V., McGill, B., Patterson, D., Potts, R., Riddle, B., Terry, R., ... Wing, S. (2014). A framework for evaluating the influence of climate, dispersal limitation, and biotic interactions using fossil pollen associations across the late quaternary. *Ecography*, 37(11), 1095-1108.

Calcagno, V., Massol, F., Mouquet, N., Jarne, P., & David, P. (2011). Constraints on food chain length arising from regional metacommunity dynamics. *Proceedings of the Royal Society B: Biological Sciences*, 278(1721), 3042-3049.

Caswell, H. (1976). Community structure: A neutral model analysis. *Ecological Monographs*, 46(3), 327-354.

Cohen, J. E. (1970). A markov contingency-table model for replicated lotka-volterra systems near equilibrium. *The American Naturalist*, 104(940), 547-560.

Connor, E. F., & Simberloff, D. (1979). The assembly of species communities: Chance or competition? *Ecology*, 60(6), 1132-1140. https://doi.org/10.2307/1936961

Cunniffe, N. J. (2002). Code to support "Metacommunity dynamics and the detection of species associations in co-occurrence analyses: Why patch disturbance matters". *Zenodo*, https://doi.org/10.5281/zenodo.6400509

D'Amen, M., Mod, H. K., Gotelli, N. J., & Guisan, A. (2018). Disentangling biotic interactions, environmental filters, and dispersal limitation as drivers of species co-occurrence. *Ecography*, 41(8), 1233-1244.

Diamond, J. M. (1975). Assembly of species communities. In J. M. Diamond, & M. L. Cody (Eds.), *Ecology and evolution of communities* (pp. 342-444). Harvard University Press.

Dubart, M., Pantel, J. H., Pointier, J.-P., Jarne, P., & David, P. (2019). Modeling competition, niche, and coexistence between an invasive and a native species in a two-species metapopulation. *Ecology*, 100(6), e02700.

Facon, B., Hafsi, A., Charlery de la Masselière, M., Robin, S., Massol, F., Dubart, M., Chiquet, J., Frago, E., Chiroleu, F., Duyck, P.-F., & Ravigné, V. (2021). Joint species distributions reveal the combined effects of host plants, abiotic factors and species competition as drivers of species abundances in fruit flies. *Ecology Letters*, 24(9), 1905-1916. https://doi.org/10.1111/ele.13825

Faust, K., & Raes, J. (2012). Microbial interactions: From networks to models. *Nature Reviews Microbiology*, 10(8), 538-550.

Forbes, S. A. (1907). *On the local distribution of certain Illinois fishes: An essay in statistical ecology* (Vol. 7). Illinois State Laboratory of Natural History.

Friedman, J., & Alm, E. J. (2012). Inferring correlation networks from genetic survey data. *PLoS Computational Biology*, 8(9), e1002687.

Gotelli, N. J. (1991). Metapopulation models: The rescue effect, the propagule rain, and the core-satellite hypothesis. *The American Naturalist*, 138(3), 768-776.

Gotelli, N. J. (2000). Null model analysis of species co-occurrence patterns. *Ecology*, 81(9), 2606-2621.

Gotelli, N. J., & Ulrich, W. (2012). Statistical challenges in null model analysis. *Oikos*, 121(2), 171-180.

Hamelin, F. M., Allen, L. J., Bokil, V. A., Gross, L. J., Hilker, F. M., Jeger, M. J., Manore, C. A., Power, A. G., Rúa, M. A., & Cunniffe, N. J. (2019). Coinfections by noninteracting pathogens are not independent and require new tests of interaction. *PLoS Biology*, 17(12), e3000551.

Hanski, I. (1982). Dynamics of regional distribution: The core and satellite species hypothesis. *Oikos*, 38(2), 210-221. https://doi.org/10.2307/3544021

Hanski, I., & Gyllenberg, M. (1993). Two general metapopulation models and the core-satellite species hypothesis. *The American Naturalist*, 142(1), 17-41.

Hanski, I., & Gyllenberg, M. (1997). Uniting two general patterns in the distribution of species. *Science*, 275(5298), 397-400.

Hastings, A. (1980). Disturbance, coexistence, history, and competition for space. *Theoretical Population Biology*, 18(3), 363-373.

Hastings, A. (1987). Can competition be detected using species co-occurrence data? *Ecology*, 68(1), 117-123. https://doi.org/10.2307/1938811

Hastings, A. (1991). Structured models of metapopulation dynamics. *Biological Journal of the Linnean Society*, 42(1-2), 57-71.

Hastings, A. (2003). Metapopulation persistence with age-dependent disturbance or succession. *Science*, 301(5639), 1525-1526.

Jentsch, A., & White, P. (2019). A theory of pulse dynamics and disturbance in ecology. *Ecology*, 100(7), e02734.

Leibold, M. A., Holyoak, M., Mouquet, N., Amarasekare, P., Chase, J. M., Hoopes, M. F., Holt, R. D., Shurin, J. B., Law, R., Tilman, D., Loreau, M., & Gonzalez, A. (2004). The metacommunity concept: A framework for multi-scale community ecology. *Ecology Letters*, 7(7), 601-613.

Levin, S. A., & Paine, R. T. (1974). Disturbance, patch formation, and community structure. *Proceedings of the National Academy of Sciences of the United States of America*, 71(7), 2744-2747.

Levins, R. (1969). Some demographic and genetic consequences of environmental heterogeneity for biological control. *American Entomologist*, 15(3), 237-240.

MacArthur, R. H., & Wilson, E. O. (1963). An equilibrium theory of insular zoogeography. *Evolution*, 17, 373-387.

Molina, C., & Stone, L. (2020). Difficulties in benchmarking ecological null models: An assessment of current methods. *Ecology*, 101(3), e02945.

Münkemüller, T., Gallien, L., Pollock, L. J., Barros, C., Carboni, M., Chalmandrier, L., Mazel, F., Mokany, K., Roquet, C., Smyčka, J., Talluto, M. V., & Thuiller, W. (2020). Dos and don'ts when inferring assembly rules from diversity patterns. *Global Ecology and Biogeography*, 29(7), 1212-1229.

Olivieri, I., Michalakis, Y., & Gouyon, P.-H. (1995). Metapopulation genetics and the evolution of dispersal. *The American Naturalist*, 146(2), 202-228.

Opedal, Ø. H., Ovaskainen, O., Saastamoinen, M., Laine, A.-L., & van Nouhuys, S. (2020). Host-plant availability drives the spatiotemporal dynamics of interacting metapopulations across a fragmented landscape. *Ecology*, 101(12), e03186. https://doi.org/10.1002/ecy.3186

Opedal, Ø. H., von Numers, M., Tikhonov, G., & Ovaskainen, O. (2020). Refining predictions of metacommunity dynamics by modeling species non-independence. *Ecology*, 101(8), e03067. https://doi.org/10.1002/ecy.3067

Ovaskainen, O., & Abrego, N. (2020). *Joint species distribution modelling: With applications in R*. Cambridge University Press.

Ovaskainen, O., Tikhonov, G., Norberg, A., Guillaume Blanchet, F., Duan, L., Dunson, D., Roslin, T., & Abrego, N. (2017). How to make more out of community data? A conceptual framework and its implementation as models and software. *Ecology Letters*, 20(5), 561-576.

Pickett, S. T., & White, P. S. (2013). *The ecology of natural disturbance and patch dynamics*. Elsevier.

Slatkin, M. (1974). Competition and regional coexistence. *Ecology*, 55(1), 128-134.

Sousa, W. P. (1984). The role of disturbance in natural communities. *Annual Review of Ecology and Systematics*, 15(1), 353-391.

Stone, L., & Roberts, A. (1990). The checkerboard score and species distributions. *Oecologia*, 85(1), 74-79.
Thompson, P. L., Guzman, L. M., De Meester, L., Horváth, Z., Ptacnik, R., Vanschoenwinkel, B., Viana, D. S., & Chase, J. M. (2020). A process-based metacommunity framework linking local and regional scale community ecology. *Ecology Letters*, 23(9), 1314-1329.

Tikhonov, G., Opedal, Ø. H., Abrego, N., Lehikoinen, A., de Jonge, M. M., Oksanen, J., & Ovaskainen, O. (2020). Joint species distribution modelling with the R-package Hmsc. *Methods in Ecology and Evolution*, 11(3), 442-447.

Tilman, D. (1994). Competition and biodiversity in spatially structured habitats. *Ecology*, 75(1), 2-16.

Tilman, D., May, R. M., Lehman, C. L., & Nowak, M. A. (1994). Habitat destruction and the extinction debt. *Nature*, 371(6492), 65-66.

Ulrich, W., Kryszewski, W., Sewerniak, P., Puchałka, R., Strona, G., & Gotelli, N. J. (2017). A comprehensive framework for the study of species co-occurrences, nestedness and turnover. *Oikos*, 126(11), 1607-1616.

Veech, J. A. (2013). A probabilistic model for analysing species co-occurrence. *Global Ecology and Biogeography*, 22(2), 252-260.

## SUPPORTING INFORMATION

Additional supporting information may be found in the online version of the article at the publisher's website.

**How to cite this article:** Calcagno, V., Cunniffe, N. J., & Hamelin, F. M. (2022). Metacommunity dynamics and the detection of species associations in co-occurrence analyses: Why patch disturbance matters. *Functional Ecology*, 36, 1483-1499. https://doi.org/10.1111/1365-2435.14047