General statistical scaling laws for stability in ecological systems

Ecological stability refers to a family of concepts used to describe how systems of interacting species vary through time and respond to disturbances. Because observed ecological stability depends on sampling scales and environmental context, it is notoriously difficult to compare measurements across sites and systems. Here, we apply stochastic dynamical systems theory to derive general statistical scaling relationships across time, space, and ecological level of organisation for three fundamental stability aspects: resilience, resistance, and invariance. These relationships can be calibrated using random or


INTRODUCTION
Ecological stability refers to a range of concepts that describe how interacting systems of species and their environments vary over time (Donohue et al., 2013;Grimm & Wissel, 1997;Holling, 1973;Ives & Carpenter, 2007;Lewontin, 1969;May, 1973;Pimm, 1984;Pimm et al., 2019).Stability metrics are primary tools used to study how systems withstand and recover from disturbances, and are therefore vital for predicting how anthropogenic pressures, such as land use change, global warming, and species loss are likely to influence natural systems (Carpenter et al., 2001;Donohue et al., 2016;Scheffer et al., 2009;Zelnik et al., 2018).However, measures of stability are highly context dependent, and vary with the spatiotemporal scale and ecological level of organisation at which measurements are made (Clark et al., 2019;Domínguez-García et al., 2019;Kéfi et al., 2019;Levin, 1992;Stommel, 1963).Measurements taken across different scales are therefore not directly comparable, which limits opportunities for cross-system comparison and synthesis (Clark et al., 2019;Csillag et al., 2000;Levin, 1992;Wang et al., 2019).This scale dependence is also a major challenge for conservation and management, because the scales that are most relevant for decision making often differ from those at which ecological systems are measured (Carpenter et al., 2001;Chesson, 2000;Isbell et al., 2018;Leibold & Chase, 2018;Levin, 1992).
To address this challenge, our goal here is to identify general statistical scaling laws for stability in ecological systems.Specifically, we seek to identify relationships that can: (1) be fitted using empirical data observed at one set of scales; (2) be extrapolated to accurately describe conditions across other scales; and (3) that are valid regardless of the underlying processes governing dynamics.
Throughout, we focus on three common stability metrics-resilience, resistance, and invariance-and explore how these change across levels of temporal, spatial, and ecological organisation (see Stability and scale definitions below) (Donohue et al., 2016;Grimm & Wissel, 1997).These metrics are especially relevant to ecological studies because they describe how systems respond to regularly occurring, externally imposed pulse perturbations, as might be expected from, for example, extreme weather events or anthropogenic disturbances such as logging, mowing, or fishing (Holling, 1973;Isbell et al., 2018;Ives & Carpenter, 2007;Zelnik et al., 2018).Moreover recent advances in stochastic dynamical systems theory have identified links between these stability metrics and basic statistical properties that can be estimated from empirical data (Arnoldi et al., 2018;Haegeman et al., 2016;Lee et al., 2020).Because these statistical properties arise directly from the mathematical definitions of variance and covariance, the resulting scaling relationships are exceedingly general (Arnoldi et al., 2019;Wang et al., 2017;Zelnik et al., 2018).
We proceed in four parts.First, we present our working definitions for resilience, resistance, and invariance, and for temporal, spatial, and ecological scale.Second, we derive the scaling relationships, and explain how they arise from the basic mathematical definitions of variance and covariance.Third, we introduce three models that simulate species abundance dynamics across a range of scales, and use these to demonstrate the scaling functions.Finally, we introduce methods for fitting the scaling relationships to observational data.These methods can be implemented via our accompanying ecostatscale R package, available via CRAN and archived on GitHub (github.com/adamtclark/ecostatscale) and Zenodo (https://doi.org/10.5281/zenodo.4626672).

M ET HOD S
The statistical scaling functions developed here are rooted in stochastic dynamical systems theory, which extends the general concept of dynamical systems modelling (e.g.ordinary differential equations, such as in logistic growth or Lotka-Volterra interactions) to include effects of stochastic perturbations on dynamic variables (Black & McKane, 2012;Gillespie, 1992;Wilkinson, 2011Wilkinson, , 2018)).In particular, we leverage the fact that when variability is grouped across different subsets of data, it adheres to rules that follow directly from the mathematical definitions of variance and covariance.These rules hold for any set of variables regardless of their distributions, so long as their mean, variance, and covariance exist and can be measured.
Editor: Tim Coulson a wide range of contexts.Moreover deviations between observed vs. extrapolated scaling relationships can reveal information about unobserved heterogeneity across time, space, or species.We anticipate that these methods will be useful for cross-study synthesis of stability data, extrapolating measurements to unobserved scales, and identifying underlying causes and consequences of heterogeneity.community, disturbance, diversity, invariability, invariance, population, resilience, resistance, spatial, temporal As an example, imagine a set of n timeseries {y 1 , y 2 , … y n }, each representing a variable with its own mean and variance.For our purposes, each y i might represent abundance dynamics for different species.By the definition of covariance, the variance of the summed abundance across all species is where <var(y i )> and <cov(y i , y j )> indicate, respectively, average species-level variance and between-species covariance.Following Wang et al. (2017), this relationship can be written as where <ρ(y i , y j )> = <cov(y i , y j )>/√(<var(y i )><var(y j )>), which is roughly analogous to the Pearson correlation coefficient (see Appendix A for more details).Critically, the relationships in Equations 1a-1b hold regardless of the underlying structure of the groups, or the processes by which data were generated.

Stability and scale definitions
Our next step is to define stability and scale in terms that can be related to the statistical rules in Equations 1a-1b.Specifically, we assume that species abundances, and the effects of disturbances on those abundances, can both be defined in terms of statistical distributions with a measurable mean and variance.In practice, this assumption implies that: (1) species abundances fluctuate around a fixed average value; (2) disturbances can be abstracted into discrete pulses that lead to near-instantaneous changes in species abundances; and (3) disturbances do not alter the processes underlying abundance dynamics (e.g.reproductive rates).We will eventually relax some of these assumptions, though even in their strict form, they apply broadly across abundance dynamics that lack strong monotonic trends, and to disturbances that act quickly relative to species growth rates.
We choose to focus on resilience, resistance, and invariance both because they are widely used in ecology (Donohue et al., 2013;Grimm & Wissel, 1997;Ives & Carpenter, 2007;Pimm, 1984;Pimm et al., 2019), and because they are particularly well suited to statistical abstraction (see examples in Figure 1; Arnoldi et al., 2018;Lande 1993;Turelli 1986).Resilience (r) describes the per-capita rate at which abundances recover following a disturbance -high resilience thus implies fast recovery.Given N(t) as abundance at time t, we can therefore calculate r = log(N(t + τ)/N(t))/τ for a sufficiently small unit of time τ.Resistance (σ −1 ) is inversely related to the effect that disturbances have on abundances-thus, high resistance implies either that disturbances are weak, or that they weakly impact species.To quantify disturbance impacts, we define σ = sd(N(t i ) − N(t i − τ)), where t i is the moment of time at which a disturbance occurred, t i − τ is the time immediately before that disturbance, and sd is the standard deviation.Invariance (var(x) −1 ) is inversely related to how widely abundances fluctuate over time-high invariance thus implies that abundances remain relatively constant, and var(x) is simply temporal variance calculated across timeseries x.To avoid confusion, we use the term "invariance" rather than the more common "invariability," as the latter often signifies a specific metric standardised by mean abundance.
We focus on three types of scaling, chosen both because they are widely used in empirical studies, and because they Working definitions of different aspects of stability for this paper.Black line shows dynamics for a basic linear process with equilibrium at x = 0 and subjected to repeated disturbance events, as introduced in the single-patch, single-species model in Equation 3. Disturbance impacts on x are drawn from a normal distribution with mean μ = 0 and standard deviation σ, and disturbances occur with average frequency λ.In this framework, resilience (r) describes rate of return towards equilibrium following a disturbance; resistance (σ −1 ) describes ability to avoid displacement by disturbances; and invariance (var(x) −1 ) describes the tendency to stay near the equilibrium value (blue dashed lines show standard deviation of the abundance fluctuations).Here, and in subsequent figures, we use resistance −1 and invariance −1 to signify these inverse relationships relate to grouping data into different subsets, which lends itself to the definitions in Equations 1a-1b (see examples in Figure 2).Temporal scale describes the rate at which observations are made over time (e.g.measurements per year).Thus, increasing temporal scale leads to larger time gaps between observations.Spatial scale describes the grain at which spatial information is aggregated (e.g.plot size, sample volume).Thus, increasing spatial scale involves grouping information across a larger number of replicates.Ecological scale describes the number of species across which measurements are aggregated (e.g.species vs. functional groups).Thus, increasing ecological scale involves grouping information across a larger number of species.

Statistical scaling relationships
Because they are primarily related to how data are grouped across observations, spatial and ecological scales have similar relationships with stability, which follow directly from Equations 1a-1b, and are described below.Conversely, because our definitions of stability are themselves functions of time, temporal scale has somewhat different effects, primarily related to measurement biases.These are described in the Measuring stability section, below.
The simplest scaling relationship to arise from Equations 1a-1b is that for invariance.Imagine that we have a series of observations of invariance measured at scale b, and wish to extrapolate these to some larger aggregate scale B. For example, B might represent a landscape (e.g. Figure 2f), and b could represent individual plots within that landscape (e.g. Figure 2d).Substituting var(y) with plot-level variance in Equation 1b yields where <var(x b )> and <ρ(x b )> represent, respectively, the average variance and average pairwise correlation measured at scale b. (2a) Working definitions of temporal, spatial, and ecological scales for this paper.Black shaded regions and lines highlight individual measurements at each scale.(a-c) Temporal scale describes the sampling interval (e.g.daily vs. annual sampling).(d-f) Spatial scale describes the total area across which available information is aggregated (e.g.plot size).(g-i) Ecological scale describes the number of species across which information is combined (e.g.species vs. functional groups).Note that although the scale of measurements varies among panels, the total simulated extent is constant within each type of scaling scenario (i.e. total time span, spatial area, or community size).For easier visualisation, results in panels (d-i) are standardised by the area or the number of species sampled (a and m, respectively) A similar relationship arises for resistance.If we replace var(y) in Equation 1b with σ 2 , that is, the observed variance in disturbance impacts on abundances, we find where <σ b > is the mean standard deviation of disturbance impacts observed at scale b, and <ρ(ξ b )> is the average pairwise correlation between observed disturbance impacts.For example, in a plant community, σ b might represent the average distance that extreme weather events push species away from their average abundances, whereas ρ(ξ b ) would indicate whether species respond similarly (ρ(ξ b ) > 0) or differently (ρ(ξ b ) < 0) to those events.
Finally, resilience scaling relationships are somewhat more complicated.Given strong effects of dispersal or species interactions, r can vary greatly across space and time, even within a single measurement scale (Arnoldi et al., 2018).While the precise relationship varies with context (see Appendix B), the median expectation for an aggregate spatial or ecological scale B (e.g. total landscape or total community abundance) can be approximated as where λ is average disturbance frequency.Note that Equation 2c is effectively the ratio of the variance vs. resistance scaling relationships in Equations 2a-2b.Unlike Equations 2a-2b, the relationship in Equation 2c formally applies only to systems with linear dynamics (see Appendix F).However, in practice, the relationship tends to hold across more complex systems (Arnoldi et al., 2019;Zelnik et al., 2018), as we will demonstrate in the modelling tests below.

Modelling tests
To illustrate the statistical scaling relationships presented above, we consider three model structures, describing: (1) a single species occupying a single patch; (2) a single species occupying multiple patches; and (3) multiple species occupying a single patch.These roughly correspond to the simplest structures needed to explore temporal, spatial, and ecological scaling, respectively.Although less general than the statistical scaling functions themselves, we choose these basic models because they allow us to analytically relate expected scaling relationships to simple functions of model parameters, against which the predictions of the statistical scaling functions can be clearly compared.Detailed model derivations, along with comparisons to the classic logistic and Lotka-Volterra models, are available in Appendices B and C.
Our first model of a single species and patch describes a linear stochastic dynamical system near equilibrium (Arnoldi et al., 2018(Arnoldi et al., , 2019;;Lande et al., 1999;Lee et al., 2020).We consider dynamics of standardised abundance x(t) = N(t) − K, where N(t) is species abundance at time t, and K is carrying capacity.In the absence of disturbances, the system is drawn towards the equilibrium x(t) = 0 (i.e.N(t) = K).Thus, x(t) describes the difference between current and equilibrium abundance (Figure 1).In this model, dynamics follow where r is the intrinsic growth rate, and ξ(t) is a stochastic function representing the effects of disturbance events.Disturbances occur at discrete moments in time, t i , with time between disturbances drawn from an exponential distribution with frequency λ (Gillespie, 1976(Gillespie, , 1992)).The effects of disturbances on species abundances are drawn from a normal distribution with mean μ = 0, and standard deviation σ.
Our second model extends Equation 3to consider dynamics of a single species in multiple patches.Here, dynamics follow where x k = N k -K k describes standardised abundance in patch k, D k is dispersal rate out of patch k, a is patch area, A is the area of the total focal region, and r k and K k are patchlevel growth rate and carrying capacity.Disturbances occur simultaneously across patches, but disturbance impacts on each patch, ξ k (t), are drawn from a multivariate normal distribution with inter-patch covariance cov(ξ k , ξ l ).Thus, correlation in stochastic forcing describes, for example, multi-patch responses to a shared environmental disturbance, such as a flood or drought.
Finally, our third model considers multiple interacting species within a single patch.Here, for a community of M total species, dynamics for species i follow where α ij describes the effect of species j on species i, and ξ i (t) describes disturbance impacts on species i.In this model, we linearise dynamics around the multispecies equilibrium achieved by each of the M coexisting species, such that x i = N i − N i * (i.e.distance from the multi-species equilibrium abundance).Under these circumstances, parameters r, K, and α can be interpreted in the same way as in the classic Lotka-Volterra equations.Effects of disturbances in this model also occur simultaneously across species, but their strengths differ, and are drawn from a multivariate normal distribution with covariances cov(ξ i , ξ j ).

(2b) 𝜎
Following the definitions above, stability for the model with one species and patch in Equation 3 relates directly to model parameters, with r describing resilience, resistance inversely related to σ, and invariance inversely related to var(x).As we will discuss, stability in the multi-patch and multi-species models in Equations 4 and 5 varies across spatial and ecological scales, following Equations 2a-2c.Additionally, as discussed in the Measuring stability section below, measurement biases arise in all three models as a function of temporal scale, especially for slow sampling rates.

Simulating dynamics
For each model, we simulated 1000 independent iterations of 100 time-units each, and, after allowing a 20 time-unit "burn-in" period, measured resilience, resistance, and invariance following the empirical methods described below.Simulations were conducted in R version 3.6.1 (R Development Core Team, 2019).Source code is archived on Zenodo (https://doi.org/10.5281/zenodo.4626668).
For the single species, single patch model, we varied temporal scale by "sampling" the simulated timeseries at frequencies ranging from every 0.01 to every 10 timeunits, with growth rate r = 1, standard deviation of disturbance effects σ = √0.1, and average frequency of disturbances λ = 1.Note that because both this model and the multi-species model in Equation 5 can be expressed fully in terms of standardised abundances x or x i , dynamics can be simulated without directly including effects of carrying capacity K.
For spatial scale, we used Equation 4 to simulate a single species in 30 separate patches, and "sampled" regions that included anywhere from 1 to all 30 patches, with r k and K k drawn from normal distributions with mean 1 and standard deviation 0.1, and inter-patch disturbance covariance cov(ξ k , ξ l ) = σ 2 /2 = 0.05 for all patches (i.e.positive correlations among disturbances, with <ρ(ξ b )> = 1/2).For the results presented in the main text, we set dispersal rate D = D k = 1 for all patches (i.e.global dispersal).In the supplement, we also present results for global dispersal rates ranging from D = 0 to D = 2, directional dispersal where D k differs among patches, and dispersal in an open system, where we included an additional loss rate, −LN k , on the right-hand side of Equation 4 to represent net loss of propagules dispersing out of the focal region.
For ecological scale, we used Equation 5 to simulate a single patch with a community of 30 species, and "sampled" groupings containing anywhere from a single species to all 30, with between-species disturbance covariance cov(ξ i , ξ j ) = 0 for species pairs (i.e.no correlation in effects of disturbances).Growth rates r i were again drawn from a normal distribution with mean zero and standard deviation 0.1, and interspecific interaction strengths α ij were drawn from a censored normal distribution with mean −1/2, standard deviation 0.1, and upper limit zero (n.b. we avoid positive interactions in our model because these tend to cause runaway growth).For both the spatial and ecological scale simulations, we set disturbance frequency λ = 1 and used a "fast" sampling interval of one sample per 0.1 time-units in order to avoid measurement biases.

Measuring stability
If sampling rates are slow relative to system dynamics, empirical estimates of resilience and resistance can be biased.This bias occurs because repeated disturbances between observations drive the system away from equilibrium, therefore leading to underestimates of resilience (i.e.r), whereas lags between disturbance events and observations lead to underestimation of disturbance effects (i.e.σ).
In contrast, because sample-size corrected variance is an unbiased estimator, different sampling frequencies do not bias measurements of invariance (i.e.var(x) −1 ).
Importantly, we can leverage information from unbiased estimates of variance to correct biased estimates of resilience and resistance.Following the properties of point processes (Arnoldi et al 2019;Zelnik et al., 2018), average temporal variance of x in the single species, single patch model in Equation 3 is Thus, variance decreases, and invariance increases, with higher resistance (lower σ), higher resilience (larger r), and slower disturbance frequency (lower λ).Building on this formula, we can derive an estimate for how the distance of x from equilibrium changes over time, yielding where t represents the time of an observation, t + τ is the time of the observation made immediately after t, p(t + τ) is the number of disturbances that occurred between t and t + τ, and 1/Λ is the average waiting time between disturbances.For our simulations, we applied this function to observations of individual species and patches to obtain unbiased estimates of r and σ.Because this expression is non-separable, we estimated parameters using a nonlinear optimiser (either nls in the stats package, or gnls in the nlme package in R; Pinheiro & Bates, 2000).More details, and automated functions for applying these corrections, are available in Appendix E and in the ecostatscale R package.

R E SU LT S
Matching theoretical expectations for our singlespecies single-patch model, simulations showed no temporal scaling relationships for invariance, nor for biascorrected estimates of resilience and resistance based on Equation 7 (Figure 3).Raw estimates of both were biased (Figure 3a,b), and for temporal scales above ~3, raw resistance could not be computed, as no disturbance-free reference time steps were available.For all metrics, uncertainty increased at slower sampling rates.
For the spatial model with global dispersal D k = 1 for all patches, both σ B and var(x B ) increased (and thus, resistance and invariance decreased) with spatial scale, matching expected relationships from Equations 2a-2b (Figure 4b,c).Following the expected scaling relationship in Equation 2c, resilience declined slightly with spatial scale, eventually settling on the average inter-patch growth rate r k = 1 (Figure 4a).There was no effect of changing global dispersal rates, nor of allowing variable dispersal rates among patches, largely because of compensatory effects at the patch vs. landscape levels (see Figure S1a,b, and more detailed discussion in Zelnik et al., 2018).In contrast, for the model with an additional F I G U R E 3 Effects of temporal scale on observed stability, following the working definitions in Figures 1 and 2. Recall that resilience is directly proportional to r, whereas resistance and invariance are inversely related to σ and var(x), respectively (signified with "−1" superscript).Shaded regions in (a-b) show ±one standard deviation from the mean for corrected parameter estimates, following Equation 7. Hatched regions show raw estimates of r and σ, calculated by comparing abundances in sequential time windows with differing numbers of disturbances-note resulting bias (see main text for details).Shaded region in (c) shows temporal variance calculated directly from the raw timeseries-note that this estimate is unbiased.Intervals are calculated based on 1000 model iterations simulated over 100 time-units, with sampling interval specified on the horizontal axis, and parameter values r = 1, D = 0, μ = 0, σ 2 = 0.1, and λ = 1.Dashed lines show true values for r and σ, and analytical expectation for var(x) based on Equation 6F I G U R E 4 Effects of spatial scale on observed stability, following the working definitions in Figures 1 and 2. All results are for a system with 30 patches, but measured at different spatial scales.(a) Resilience declines somewhat as a function of spatial scale, due to stronger buffering effects among plots at smaller scales.(b-c) Resistance and invariance differ across scales due to changes in total abundance and between-patch covariance.In all panels, dashed lines show analytical expectation for scaling relationships, following Equations 2a-2c.Sampling interval is 0.1 time-units for all simulations, with global dispersal (D k = D = 1 for all patches) and disturbances with positive between-patch covariance (cov(ξ) = σ 2 /2).Patch-level growth rate r k and carrying capacity K k are drawn from normal distributions with mean 1 and standard deviation 0.1.See Figures S1a-c  loss rate from dispersal outside of the focal area, mean abundance declined as a function of loss rate, leading to large changes in stability estimates (Figure S1c).
Finally, for the multi-species model (Figure 5), both r and σ B increased (and thus, resilience increased, and resistance decreased) with ecological scale following the expected scaling relationships in Equations 2b-2c.Invariance also matched the expected analytical relationship Equation 2a, resulting in a concave-down hump-shaped relationship for var(x B ).

DI SC US SION
Our results demonstrate that resilience, resistance, and invariance can be related to simple statistical properties of dynamical systems, which allow general and robust scaling of these stability metrics across time, space, and level of ecological organisation.Recall that following our definitions, resistance describes the immediate effect of disturbances on abundance, resilience describes the rate at which abundances from disturbances, and invariance describes the joint effects of these two processes on dynamics over time.The scaling functions in Equations 2a-2c demonstrate how to extrapolate all three of these metrics across spatial and ecological scales.Similarly, the relationships among these variables that we summarise in Equation 6, and the correction introduced in Equation 7, demonstrate how to overcome measurement biases that arise as a function of temporal scale.
Jointly, these methods allow empirical estimates of resilience, resistance, and invariance measured at one set of scales to be extrapolated out to other, unobserved scales.These results are therefore important both because they highlight the fact that stability measurements are scale dependent, and because they demonstrate how to correct for this scale-dependence when comparing and interpreting empirical results.The general implications and limitations of these methods are discussed below, and a guide for practitioners showing how to apply our methods using the ecostatscale R package is presented in Box 1.

General scaling relationships
Because they arise from basic mathematical definitions of variance and covariance, there are several important cases for which the scaling relationships in Equations 2a-2c hold exactly, regardless of the underlying processes governing system dynamics.Most obviously, stability measurements from smaller scales can always be extrapolated to larger scales provided that both include the same total extent (e.g. using the full set of small plots in Figure 1d to estimate stability of the single large plot in Figure 1f).However, even if only partial sampling is available, the scaling relationships provide unbiased extrapolations to larger scales so long as smaller scale sampling is representative (i.e. the average conditions at the F I G U R E 5 Effects of ecological scale on observed stability, following the working definitions in Figures 1 and 2. All results are for a system with 30 interacting species, but measured at different ecological scales.(a) Resilience varies with ecological scale in this model because of compensatory dynamics among species.(b) Resistance depends only on instantaneous responses to the disturbance regime, and therefore varies with ecological scale, but not with species interactions (n.b. dashed and dotted lines overlap perfectly).(c) Invariance varies quadratically with ecological scale as a function of average species-level variance in abundance, and between-species correlation-negative correlation implies a concavedown relationship (as shown here), whereas positive correlation implies a concave-up relationship (see Figure S2 in the appendix for other examples from empirical datasets).Note that the vertical axis is log-transformed for easier visualisation.For all panels, dashed lines show analytical expectation for scaling relationships, following Equations 2a-2c.For reference, dotted lines show hypothetical scaling relationships for a system without interactions.Sampling interval is 0.1 time-units for all simulations, with species interactions (α ij ) drawn from a truncated normal distribution (mean = −1/2; standard deviation = 0.1; upper limit = 0), and species growth rates r i drawn from a normal distribution with mean 1 and standard deviation 0.1.Otherwise, intervals and parameters are as described in the legend for Figure 3

Harmonising data
Differences in observational scales greatly influence stability measurements, which presents a problem for synthesis and metanalysis because available data are often measured at many different scales.For example, consider the effects of ecological scales (Figure 5): because all three metrics of stability vary as a function of the number of species included in samples, studies with different sampling designs could reach wildly different conclusions about stability, even if they were all conducted in the same system.A potential solution is to "harmonise" records across studies, by extrapolating estimates to a unified set of scales.Doing so removes the statistical effects of scale differences, meaning that remaining differences should be ecologically meaningful.To improve interpretability, we suggest that harmonisation should generally proceed by extrapolating smaller-scale studies up to larger scales (see General relationships).For example, studies are available with plot sizes of 1, 5, and 10 m 2 , then stability estimates for all three studies should be harmonised to a scale of 10 m 2 .Unbiased estimates of resilience, resistance, invariance are independent of temporal scale, meaning that harmonisation is not necessary to account for differential sampling rates (however, a bias correction such as that in Equation 7, and implemented in the xt2fun function, may still be necessary if sampling is slow relative to system dynamics).In contrast, stability measurements do vary with spatial and ecological scale, meaning that harmonisation is required to account for differences in observational scale.Harmonisation is accomplished by taking available estimates of invariance, resistance, or resilience, plugging them into Equations 2a-2c, respectively, and substituting b with the scales of observations, and B with the larger scale to which observations are to be harmonised.These computations can be carried out using the var_scale, sd_scale, and res_scale functions, respectively.See the help documentation in the ecostatscale R package for more details.

Extrapolating estimates
Equations 2a-2c and the corresponding R functions also allow extrapolation of expected stability beyond the scale of observations.These extrapolations can be useful for, for example, estimating stability at the landscape level based on plot-level observations, or at the community level given observations of a subset of species.As with harmonisation, extrapolation is accomplished by substituting b with the scale of observations, and B with the desired scale at which extrapolations are to be made.Importantly, these extrapolations only describe "average" conditions across observations.Thus, in order for extrapolations to be accurate, observations must be representative of the extrapolated region or species-that is, average conditions across the samples must be the same as the average conditions at the extrapolated scale.
In general, random samples of plots or species should achieve this criterion.However, if heterogeneity is structured across space or time, then this must be accounted for in the sampling design.For example, if a landscape is known to be 20% mesic and 80% xeric, then sampling should be stratified to include this ratio of site types, even if it differs from local conditions near survey plots.

Detecting heterogeneity
Although unobserved heterogeneity can bias extrapolations of stability, data measured across several scales can be leveraged to identify underlying sources of heterogeneity and to validate extrapolations.To do so, we suggest parameterising the statistical scaling relationships in Equations 2a-2c, or the corresponding R functions, with data from the smallest observational scale available, and comparing these to independent estimates of stability observed at larger scales (Pimm et al., 2019;Wang et al., 2017).If empirically observed relationships diverge from statistical expectations, then the associated scaling relationship is likely to be biased, and should not be used to extrapolate stability to other, unobserved scales.Even if observations are only available from a single scale, it may still be possible to identify heterogeneity by aggregating specific combinations of plots or species through "non-random cross-validation" (Wenger & Olden, 2012).For example, if stability varies among plots in different regions of a site (e.g.northwest vs. southeast corner of a field), or among different species groups (e.g.legumes vs. grasses), then these deviations indicate spatial or ecological heterogeneity, respectively.Similarly, if stability varies between estimates derived from the full dataset vs. for example, a subset of the dataset including every second sampling event, then this indicates that temporal sampling was potentially too infrequent to produce unbiased estimates of resilience and resistance, and a correction such as that in Equation 7, or the accompanying xt2fun R function, should be applied.See function help documentation for example implementations.
smaller scale match those at the larger scale).Thus, the scaling definitions should be applicable for most studies that randomly sample species from a larger community, or plots from a larger landscape-which includes a substantial fraction of existing ecological timeseries datasets.For similar reasons, our methods can also be used to project stability measurements from larger scales to smaller scales, although these estimates summarise average conditions at smaller scales-not conditions within individual patches or for individual species.
Another general result is that if sampling is slow relative to dynamics, raw estimates of resilience and resistance will be biased (leading to under-estimates of r and σ).In theory, the bias correction that we present in Equation 7 holds only for simple linear systemshowever, in practice, it will often hold for more complex systems so long as they are undergoing bounded fluctuations around a static mean (Arnoldi et al 2019;Zelnik et al., 2018).For example, estimates from Equation 7converged to the correct parameter values for both our multi-patch and multi-species models.
Our study also highlights the similar effects of spatial vs. ecological scale on stability.For invariance, as predicted in Equation 2a, positive average covariance among species-level abundances or site-level abundances leads to a concave-up scaling relationship for var(x) (e.g. Figure 4c, as might be driven by species in nearby plots all responding similarly to a local disturbance).Negative average covariance leads to a concave-down relationship (e.g. Figure 5c, as might be driven by strong effects of competition).Recall, however, that for the simulation results presented here, we assume that species responses to disturbances are independent (i.e.cov(ξ i , ξ j ) = 0), which allows competition effects to dominate.In real-world systems, these responses are generally thought to covary positively over time (e.g.due to similar effects of droughts or floods across species), leading to net positive average covariance in species abundances, and thus a concave-up relationship for var(x) as a function of ecological scale (Houlahan et al., 2007;Loreau & de Mazancourt, 2013).Examples of concave-up, concave-down, and linear scaling relationships derived from empirical data are shown in Figure S2 in the supplement.
The spatial and ecological scaling relationships for resistance in Equation 2b are similar to those for invariance.However, because we define resistance as the instantaneous effect of disturbances on abundances, it is largely independent of internal system processes, such as inter-patch dispersal or species interactions.Scaling relationships are therefore the same for systems with and without species interactions (Figure 5b; n.b. dashed and dotted lines overlap perfectly).This property means that extrapolating resistance across spatial and ecological scales should be relatively straightforward even in complex systems.
The scaling relationship for resilience is the most variable.For temporal scaling, average per-capita growth rates in our model are constant across space and time, leading to no change in resilience with scale (Figures 3a  and 4a).For multi-patch and multi-species communities, the median scaling relationship for resilience can be approximated following Equation 2c, although for multispecies communities, this relationship only describes the instantaneous recovery rate (i.e. as would be observed immediately after a disturbance).This is because resilience also varies as a function of temporal scale in these systems (see Figure S3)-that is, different return rates occur over short vs. long timespans, because interactions among patches or species can lead to complex, non-monotonic dynamics (see Neubert & Caswell, 1997).In theory, this temporal scaling relationship can be expressed analytically (Eq.6b in Arnoldi et al., 2018), but in practice, doing so requires substantial a priori information, suggesting that empiricists should generally be careful to choose sampling intervals that match the scale of the phenomena being tested when studying highly interconnected or diverse systems.

Interpreting deviations from scaling rules
An important limitation of the statistical scaling relationships in Equations 2a-2c is that they are largely phenomenological.Although they predict how patterns change across scales, they do not explain why these patterns exist.Indeed, there are usually many potential mechanisms that could explain patterns equally well.For example, although a concave-up relationship between var(x) and ecological scale implies negative average correlations in interspecific abundances, these correlations could result from many different mechanisms, including competition, species differences in their response to disturbances, or observation error (Clark et al., 2019;Lee et al., 2020).
Critically, cases where observed scaling relationships diverge from statistical expectations can be especially informative.Because Equations 2a-2c hold only if samples are representative of the full extent under consideration, divergence from these expectations implies that there is unobserved heterogeneity in observations across time, space, or species.For example, if scaling relationships predicted from patch-level or species-level data do not match larger-scale observations (e.g. from remote sensing data, or at the community level), then this might indicate that samples are not representative of conditions at the larger-scale (König et al., 2017;Wang et al., 2017).Similarly, if observed invariance changes as a function of temporal scale, then this divergence may be indicative of regime shifts (Scheffer et al., 2009), or that equilibrium abundances are not static (Pimm et al., 2019).
Note, however, that without additional information, the statistical scaling relationships cannot predict how unobserved heterogeneity outside of the sampled region will influence predictions.In diverse communities, for example, stability varies with community structure, meaning that adding or removing species, or changes in how species interact with one another, will alter the scaling relationships (Arnoldi et al., 2018;Hillebrand et al., 2018;Jacquet et al., 2016).Similarly, if the structure of heterogeneity changes across space or time (e.g.clustered disturbances in one region, vs. disaggregated disturbances in another), then resulting biases will be, for all practical purposes, impossible to correct without a priori knowledge of conditions across the entire landscape (Arnoldi et al., 2019;Hallett et al., 2019;Lee et al., 2020).

Nonstationary systems
In this study, we largely neglect effects of non-stationary behaviour on stability (e.g.dynamics with no fixed equilibrium point).We do so because these kinds of dynamics are generally difficult to summarise through simple statistical relationships.However, it is possible to incorporate some kinds of non-stationary behaviour into our methods.If mean abundances change slowly over time, then our scaling functions can still be applied within subsets of the data.For example, in systems undergoing regime shifts, our scaling functions could be fitted separately for the time periods before vs. after the shift (Carpenter et al., 2001;Scheffer et al., 2009).Similarly, even for systems undergoing rapid changes due to anthropogenic pressures, our methods likely still apply for historical data preceding current disturbance regimes (Coulson, 2021).More broadly, if trends in average abundances are predictable, then we can redefine system dynamics around this moving average.We demonstrate this property in Figure S1c, where losses due to dispersal out of the system lead to declines in average abundance below the expected equilibrium at N k (t) = K k .Without any additional corrections, this shift in equilibrium causes stability scaling relationships to diverge greatly from statistical expectations.However, if we redefine x around the new mean abundance value (i.e.x k,new (t) = N k (t) − <N k >) then we find that the statistical scaling relationships still hold.Using more sophisticated approaches for predicting changes in mean abundance values over time (e.g.Cenci & Saavedra, 2019;Chesson, 2017;Deyle et al., 2016;Hamilton et al., 2017;Karakoç et al., 2020), it may be possible to similarly "detrend" even very complex nonstationary data to make them suitable for our methods.

Conclusions
The strong dependence of ecological stability on scale and context has represented a major challenge for attempts to synthesise and extrapolate measurements from real-world systems.The methods that we present here greatly alleviate these challenges, and provide an opportunity to harmonise measurements from across different sites and systems, and to extrapolate measurements out to scales which are difficult to measure.We are therefore hopeful that these results will be especially helpful in facilitating cross-system comparisons of stability, providing estimates at scales that are more relevant for conservation and management, and identifying underlying drivers of heterogeneity across time, space, and species.
in the supplement for examples of other types of dispersal.Otherwise, intervals and parameters are as described in the legend for Figure3