A unified framework for species spatial patterns: Linking the occupancy area curve, Taylor's Law, the neighborhood density function and two‐plot species turnover

Abstract The description of spatial patterns in species distributions is central to research throughout ecology. In this manuscript, we demonstrate that five of the most widely used species‐level spatial patterns are not only related, but can in fact be quantitatively derived from each other under minimal assumptions: the occupancy area curve, Taylor's Law, the neighborhood density function, a two‐plot variant of Taylor's Law and two‐plot single‐species turnover. We present an overarching mathematical framework and derivations for several theoretical example cases, along with a simulation study and empirical analysis that applies the framework to data from the Barro Colorado Island tropical forest plot. We discuss how knowledge of this mathematical relationship can support the testing of ecological theory, suggest efficient field sampling schemes, highlight the relative importance of plot area and abundance in driving turnover patterns and lay the groundwork for future unified theories of community‐level spatial metrics and multi‐patch spatial patterns.


Defining units of length and area
For notational convenience, and without loss of generality, we define units of length in our landscape such that the species of interest has a density of 1 individual per unit area. Under this definition, the first-order intensity D 1, and the expected number of individuals in a plot of area A is given by .A/ D A D A. Aside from notation, this approach also highlights that every species having a metric with the same parameters (for example, the a and b of Taylor's Law) is expected to show the same general form of the five spatial metrics, regardless of its abundance.
A point pattern with any original first-order intensity 0 , as measured in a landscape in which length is measured in any original units, can be converted to a point pattern in which D 1 by multiplying all measures of length by the conversion factor p 0 . To verify this conversion, consider a square plot that contains an expected number of individuals EOEn. This expectation is invariant to the units in which lengths or distances are measured. In particular EOEn D h 02 0 D h 2 , where the superscript 0 refers to side length and density measured in original units and the non-superscript terms refer to side length and density measured in scaled units. Since in scaled units D 1, it is clear that multiplying the original arbitrary length h 0 by p 0 gives distances measured in scaled units h that ensure that D 1.
For a concrete example, consider an original landscape in which 0 D 4. A square plot of side length h 0 D 2 has EOEn D h 02 0 D 16. Multiplying all distances in the landscape by the factor p 0 D 2 gives a square plot of side length h D 4 in this scaled landscape. In this same plot with distances measured in scaled units, EOEn D 16 D h 2 D 16 and thus D 1, as per the goals of the unit conversion.
Although this conversion is convenient for notational purposes, the validity of the framework described below does not depend on it. To construct or apply the equations below in original, unscaled units, simply replace all occurrences of the variable A with 0 A 0 and all occurrences of the variable h with h 0 p 0 , where the superscript 0 refers to areas and lengths measured in original units.
An additional notational convenience of this unit conversion is that for a landscape where D 1, the neighborhood density function .r/ is identical to the second-order intensity function, 2 .r/. Below, we derive results for the second-order intensity function and report these as the results for both .r/ and 2 .r/ in the main manuscript.

Relationship between plot mean and variance and second-order intensity function
One of the key expressions in this manuscript is the relationship between plot variance and the integral of the secondorder intensity function (Illian et al. 2008;Diggle 2013). This relationship can be understood intuitively by considering the expected number of pairs of points found in a plot.
Consider a single plot of any area and geometry containing some number of individuals of a species, with the number of individuals represented by the random variable n. Define the expected number of unordered pairs of points in this plot as EOER. From basic rules of expectations, we have where and 2 are the mean and variance, respectively, of n.
Next, within this same plot, consider two infinitesimally small discs of area dx 1 and dx 2 , which are small enough that they can only contain zero or one individual. Imagine moving each of these two discs to every possible location within the plot. For each possible pair of locations for the two discs, the discs are separated by a distance r. The probability of encountering an individual in both discs simultaneously, which is equivalent to the expected number of pairs of individuals in the pair of discs, is then equal by definition to 2 .r/. By integrating over all possible locations for each disc, we can thus also calculate the expected number of unordered pairs of individuals in the plot as where L refers to the region of the plot. In a two-dimensional case, for example, the integral above can be written as a quadruple integral over the x and y coordinates of the two discs.
Equating the two equations above leads to the relationship between the mean and variance of the quadrat count distribution and the integral over the second-order intensity function that is described in the main manuscript. Similar logic leads to the expression for the covariance in abundance across two plots as an integral of the second-order intensity function, with the only difference being that the discs dx 1 and dx 2 are constrained to be located in two different plots.
In practice, the most substantial challenge in applying this framework is evaluating the integrals, or in some cases solving the integral equations, resulting from the above. One approach is direct integration, in which the integral is taken over the coordinate locations of the discs dx 2 and dx 2 . As this leads to a double integral for one-dimensional landscapes or a quadruple integral for two-dimensional landscapes, analytic solutions do not appear to be available for many cases. A second approach is based on the line-line picking distribution, which gives the probability that a randomly selected pair of locations (not individuals or points of a point pattern) within a geometrical shape will be separated by a distance r (Mathai 1999). By multiplying these probabilities by the total number of possible locations for the two discs, the number of pairs of locations at each distance r can then be calculated. The multiple integrals resulting from the above can then be expressed as a single integral over the set of allowable distances r between locations at which a point might appear. Both approaches are used in the derivations below.
3 The univariate negative binomial distribution and the occupancy area curve For a single plot of area A, we assume that the quadrat count distribution follows a negative binomial distribution. The negative binomial distribution is given by where n is the abundance of a plot and k (often written as r in a non-ecological context) is interpretable as an aggregation parameter.
For describing occupancy, we are specifically interested in the case of n D 0, which occurs with the probability The mean and variance of the negative binomial distribution are given by D pk=.1 p/ and variance 2 D pk=.1 p/ 2 , leading to the alternative parameterization The occupancy area curve, ‰.A/, is then defined as (He & Gaston 2003) ‰.A/ D 1 P .n D 0jA/ D 1 Since .A/ D A when D 1, as described above, this simplifies to This general equation can be used to calculate the occupancy area curve for all theoretical examples in the main text, given knowledge of Taylor's Law for that particular example.

The bivariate negative binomial distribution
As described in the main text, the traditional form of the bivariate negative binomial distribution (Johnson et al. 1997) does not allow for arbitrary specification of positive correlation. We thus define an alternative bivariate negative binomial distribution based on an additive combination of three independent underlying negative binomial distributions, as described below.
Presume that we have known values for the mean , variance 2 , and covariance COV in abundance for a pair of plots of area A and distance between plots D. Based on these known values, our specific goal is to determine P .n > 0; n 0 > 0/, the joint probability that both plots in a pair are occupied. This expression, along with a known expression for ‰.A/, determines two-plot turnover.
Define the random variable n as the abundances of a species in the first of a pair of plots (all results below apply symmetrically to the second plot in the pair, with abundance n 0 ). When considering only this single plot in isolation, n represents a draw from a negative binomial quadrat count distribution with mean and variance 2 . We consider the random variable n to represent a sum of two underlying random variables, n S and n L . Both n S and n L are negative binomial distributed and are independent of each other. We then have We interpret the random variable n S as a set of individuals in a plot that result from a "small scale" population process (birth, death, and dispersal) that applies only to that plot, independent of the second plot of the pair. This includes, for example, population dynamics internal to the plot as well as short-distance dispersal. The random variable n L , in contrast, represents "large scale" population processes that result in an identical number of individuals being placed in each plot in the pair. The total abundance of a plot represents the contributions of these small-scale and large-scale processes, which operate independently of each other. Note that aside from this mechanistic basis for Eq. 8, this equation can also be interpreted purely phenomenologically as a means of specifying varying degrees of correlation in abundance between two plots.
Equation 8 provides a simple and convenient means of specifying arbitrary correlations between the abundances in two plots based on the relative contributions of n S and n L to n. In the limit EOEn S ! EOEn and EOEn L ! 0, the abundance of the two plots n and n 0 are independent of each other, as the abundance of each plot reflects independent draws of n S . In the limit EOEn L ! EOEn and EOEn S ! 0, both plots will have an identical abundance and hence a correlation of 1.
The probability that one plot is unoccupied is once again given by The probability that both plots are jointly unoccupied, P .n D 0; n 0 D 0/, is given by the product of the probabilities that n S D 0 in the first plot, n 0 S D 0 in the second plot, and n L D 0 (since these three random variables are all independent of each other). This leads to the expression The next step is to determine the values of the three unknown parameters k S , k L and p in terms of the mean , variance 2 , and covariance COV in plot abundance. Considering first the single plot abundance n, the parameters p and k are given by The next task is to solve for k S and k L . The covariance between the abundance in the first plot, n D n S C n L , and the abundance in the second plot, n 0 D n 0 S C n L (note that the same variable n L appears in both expressions), can be expressed as Because n S , n 0 S and n L are independent random variables, Eq. 13 reduces to COV.n; n 0 / D COV.n L ; n L / D 2 L (14) As the variance of a negative binomial distribution is 2 L D pk L =.1 p/ 2 , we can use the relationship COV D 2 L to solve for k L as The probability that both plots are jointly unoccupied can then be expressed as Finally, considering the equation above across multiple plot scales, we have D A under the unit definition above and replace COV.n; n 0 / with the function notation C.A; D/, which expresses that the covariance between n and n 0 will depend on the values of A and D for the pair of plots. This gives This general equation can be used to calculate two-plot turnover for all examples in the main text, given knowledge of Taylor's Law and the two-plot variant of Taylor's Law for that particular example.

The two-plot turnover metric
The Sorensen index of community turnover or distance decay is a widely-used metric for describing two-plot turnover in communities. This index is calculated as twice the number of species shared in common between two plots divided by the sum of the number of species in each plot. When predicted using theory, the index can be defined in terms of expectations as where EOEU is the expected number of species shared between the two plots and EOEX and EOEY are the expected number of species found in each of the two plots. When the two plots are of the same area and EOEX D EOEY , this equation reduces to S D EOEU =EOEX , which is simply the expected number of species shared between the two plots divided by the expected number of species in a single plot. Note that S will likely be dependent both on the area of the two plots and the distance between them.
Analogously for a single species, we can define a turnover metric T as the joint probability of the species occurring in both plots divided by the probability of the species occurring in a single plot.
For the four possible combinations of species occupancy across the two plots, we have Because of the symmetry between n and n 0 when both plots are the same area, this simplifies to P .n > 0; n 0 > 0/ D 1 2P .n > 0; n 0 D 0/ P .n D 0; n 0 D 0/ From the marginal probability for a single plot, we also have Substituting Eq. 22 and Eq. 23 into Eq. 20 gives Recognizing that ‰.A/ D P .n > 0jA/ D 1 P .n D 0jA/, the above can be rearranged and simplified to 6 Example: Power law Taylor's Law in one dimension

Taylor's Law
In this section, we consider the common power law form of Taylor's Law given by applied in a one-dimensional landscape. Such a one-dimensional landscape can describe either a spatial pattern, such as a species population arranged along a river, or a temporal pattern, in which case the "landscape" describes the appearence of individuals in a sample over time. This case provides a uniquely analytically tractable example, with solutions that can be compared to the approximate solutions for the two-dimensional case below. Note that in this one-dimensional case, A is defined as the length of a one-dimensional "plot," which is a line segment.

Neighborhood density function
Using the line-line picking approach (Mathai 1999), we have the general equation where N.r/ is the number of pairs of locations separated by distance r in the plot. Using Taylor's Law from Eq. 26 gives The probability distribution describing the distance r between two randomly chosen locations (distinct from the distribution of distances between known individuals or points in a point pattern) on a line segment is known as the line-line picking distribution (Mathai 1999) and is given by The total number of ordered pairs of locations in a line segment of length A is A 2 , and thus the number of pairs of locations at each distance r is given by Substituting N.r/ into Eq. 28 gives the expression This is a Volterra integral equation of the first kind with a difference kernel, which can be solved for 2 .r/. The solution is well-known to be the second derivative, with respect to A in this example, of the left hand side. The second order intensity is thus where ı.r/ is the Dirac delta function, which is equal to 1 when r D 0 and equal to zero elsewhere. We assume the convention in which the Heaviside function H.Â /, which is the integral of the delta function, has the value H.0/ D 1=2.

Two-plot Taylor's Law
As described in the main text, the second-order intensity function is also related to the covariance in abundance between two non-overlapping plots, which we define as the two-plot variant of Taylor's Law. Define D as the distance between the closest edges of two square plots, where D > 0 to avoid plot overlap. In one dimension, the minimum possible distance between a point occurring in the first plot and a point occurring in the second plot is thus D, and the maximum possible distance is D C 2A. We can then write The number of pairs of points located a distance of D or D C 2A apart approaches zero as r approaches these limits. At r D D C A, there are A such pairs of points. This leads to a two-part expression for N(r), which leads to the equation for C.A; D/ Note that we have dropped the ı.r/ term from Eq. 32 in the above, as the limits of integration above do not include r D 0 so long as D > 0. This integrates to 7 Example: Power law Taylor's Law in two-dimensions

Taylor's Law
In this section, we begin with the common power law form of Taylor's Law, once again given by Note that in this example, A represents a two-dimensional area. All plots are assumed to be square, with side length h D p A.

Second-order intensity function
Using the same line-line picking approach as in Example 1, adjusted for a square plot, we begin with the expression where h D p A is the length of a side of a square plot of area A and N.r/ is the number of pairs of locations at distance r in the plot. The distribution P .r/ for a square is known Weisstein 2020), and multiplying this by h 2 h 2 D h 4 possible pairs of locations gives There does not appear to be a simple, exact solution for 2 .r/ in the integral equation above. In the Supporting Information 2, we show that a very good approximation for 2 .r/ in this two-dimensional case is 2 .r/ cr 2.b 2/ C 1 ı.r/ r where the constant c is given by c D ab.b 1/.2b 1/ .2b 1/ 4.b 1/ As in the previous one-dimensional example, ı.r/ is the Dirac delta function and we assume the Heaviside function has the value H.0/ D 1=2.

Two-plot Taylor's Law
As it is difficult to directly generalize the line picking distribution given above for a single square to the case of two squares (Chu 2006), we take a two-part line picking approach to deriving the two-plot Taylor's Law.
Consider first two vertical line segments, each of height h and separated by distance d , standing with their bottom ends on the x-axis. The distribution of distances P .r/ between two randomly chosen points, one on each line segment, is given by . The total number of pairs of points located on these line segments is A square plot can then be defined by "sweeping" a vertical line segment of height h across a segment of the x-axis of length h. The number of points shared between two square plots can thus be determined by integrating the number of points shared in the two vertical line segments across allowable locations on the x-axis of the two vertical segments, defined such that the left-most line segment must fall within the region of the left-most plot and the right-most line segment must fall within the region of the right-most plot. The integral over the allowable distances d between the vertical line segments is identical to the integral over pairs of locations in two line segments that is used to calculate the two-plot Taylor's Law in Example 1. By similar logic, we thus have the equation Using the approximation for 2 .r/ in Eq. 40, the value of C.A; D/ can be calculated for any given plot area A and distance between plots D by numerical integration.
In Supporting Information 2, we derive an approximation for the above integral for cases in which D > h, such that plots are separated by a gap equal to at least the length of the side of a plot. This approximation is given by 8 Example: Gaussian second-order intensity

Second-order intensity function
In this section, we begin with the assumption of a Gaussian form of the second order intensity function, 2 .r/, which is given by 2 .r/ D˛e ˇr 2 C 1 (45)

Taylor's Law
Unlike the two examples above, we here use a direct integration approach to derive the remaining spatial metrics from the neighborhood density function. We begin by considering a single square plot of area A, with sides of length h, placed with its bottom left corner at the origin. The coordinates .x 1 ; y 1 / and .x 2 ; y 2 / represent the location of two small discs within the square. Defining the distance between the two discs as r D p .
x 2 x 1 / 2 C .y 2 y 1 / 2 leads to the expression Substituting the Gaussian second order intensity function leads to Rearranging leads to where erf.x/ is the error function.

Two-plot Taylor's Law
The covariance in abundance between two plots can be evaluated in a similar fashion to the variance in a single plot. We have the general expression where the right hand side differs from the previous calculation of single-plot variance only in the limits of integration for x 2 . The variable D here represents the distance along the x-axis between the square in which the first disc and second disc are located. Using a Gaussian intensity function and following the same procedure as for single-plot variance leads to 9 Example: Complete spatial randomness For comparison to the relationships above, we consider the form of the five spatial metrics in the case of complete spatial randomness (CSR), which is defined by the two equivalent statements that (a) the number of points in a plot of area A follows a Poisson distribution with mean A, and (b) each point is randomly and independently located within the larger landscape (Diggle 2013). To demonstrate the quantitative relationship between the five spatial metrics for the case of CSR, we use definition (a) above, which defines both the occupancy area curve and Taylor's Law, and use this to derive the other three spatial metrics.
The quadrat count distribution under CSR takes the form of a Poisson distribution, which is the limiting form of the negative binomial distribution as k ! 1 for a constant mean. We thus use equations for the Poisson distribution in this example where needed, rather than those for the negative binomial distribution given above.

Occupancy area curve
Under CSR, the quadrat count distribution is a Poisson distribution. In a landscape with scaled units such that D 1 and .A/ D A, we thus have

Taylor's Law
Since the variance and mean of the Poisson distribution are equal, Taylor's Law is given by

Second-order intensity function
Using Eq. 46 with 2 .A/ D A and A 2 D h 4 , we have We can see that the solution to this equation is which is consistent with the observation that points are randomly and independently distributed. In the more general case in which ¤ 1, 2 .r/ D 2 under CSR.

Two-plot turnover
Under random placement, the univariate quadrat count distribution is Poisson and the bivariate quadrat count distribution is the product of two independent Poisson random variables. As a result, P .n D 0; n 0 D 0jA; D/ D P .n D 0jA; D/P .n D 0jA; D/, with P .n D 0jA; D/ D e A as noted above. The two-plot turnover metric from Eq. 24 thus becomes This is once again consistent with the principle of CSR, as the conditional occupancy of a second plot given the occupancy of a first plot is equal to the unconditional occupancy of the second plot.

Application to Maximum Entropy Theory
The Maximum Entropy Theory of Ecology (METE) (Harte et al. 2008;Harte 2011;Harte & Newman 2014) is an ecological theory that predicts many well-known ecological patterns based on knowledge of a small number of state variables and the principle of maximizing information entropy (MaxEnt) (Jaynes 1982). Details of the theory can be found in the references above.
Here we focus on one specific prediction of the theory, which is that the quadrat count distribution, denoted ….n/ in the theory, will take the form of a geometric distribution across all spatial scales. This result arises from a relatively straightforward application of MaxEnt to the problem of determining the probability that an urn will contain a given number of balls when only the mean number of balls per urn is known. The MaxEnt prediction for this problem, which is analogous to the quadrat count distribution in the spatial framework derived here, is the geometric distribution (Jaynes 1982;Harte et al. 2008). When applied in a spatial context in METE, this distribution is not scale dependent and is thus predicted to apply at all spatial scales, for all possible mean numbers of individuals per plot.
The relationship between the variance and mean of a univariate negative binomial quadrat count distribution is 2 D C 2 =k. Using k D 1 to retrieve a geometric distribution and scaling units such that D 1 and hence D A, the METE prediction for Taylor's Law is thus .A/ 2 D A C A 2 .
Using the approach of direct integration, the relationship between Taylor's Law and the second-order intensity function is given by Eq. 46. Substituting Taylor's Law under METE into the left hand side and recognizing that A D h 2 for a square plot gives This equation takes the same form as that for the case of complete spatial randomness (Eq. 54). The solution to this integral equation is We then have the equation for covariance or two-plot Taylor's Law, which is leading to METE thus predicts that the covariance in plot abundance is a positive value that depends on plot area but is constant with distance. This theory thus makes the empirically implausible predictions that the correlation in abundances in two plots of fixed area (a) does not decrease with increasing distance, and (b) does not asymptote to zero at large distances.
Within the context of maximum entropy theories, we note that both Banavar et al. (2010) and Haegeman and Etienne (2010) apply the maximum entropy principle in a manner that instead leads to a Poisson quadrat count distribution. Under these theories, the predicted spatial distribution of individuals is thus consistent with complete spatial randomness (see Table 1 and main text).   Figure SI1.2: Influence of plot area A and intraplot distance D on values of the two-plot turnover metric. Panels show turnover metric calculated for Taylor's power law in two-dimensions and a Gaussian second-order intensity. Plot area, as opposed to intraplot distance, can be seen to be a primary influence on values of the turnover metric in both examples. For a species with a power law Taylor's Law, for example, a plot of A D 1 has a turnover metric of 0.70 at D D 1, decaying to 0.66 at large D. This turnover at large D is approximately the same as the turnover at D D 1 for a plot of A D 0:85. For these parameters, a 15% difference in plot area, or equivalently expected species abundance, gives approximately the same change in the turnover metric as a several order of magnitude increase in interplot distance. Parameters are the same as those used in Figure 3 in the main text. T Two-plot turnover Figure SI1.3: Change in shape of metrics calculated from Gaussian second-order intensity with change in˛andˇparameters. Black lines in both subplots show predicted metrics for the case of a Gaussian-second order intensity function with parameters D 2:7 andˇD 0:8, equal to the median parameters across the Barro Colorado Island data set. In (a), the parameterˇis held constant at this value, and colored lines show predicted metrics for a one order of magnitude increase (red dashed), a two order of magnitude increase (red dotted), a one order of magnitude decrease (blue dashed), and two order of magnitude decrease (blue dotted) in the parameter˛relative to the black line. In (b), similar calculations are performed but holding˛constant and allowingˇto vary. In general, large variation inˇhas less of an effect on the predicted metrics than variation in˛. Figure 5 shows that disagreements in the estimate ofˇacross empirical metrics are generally larger than disagreements in the estimate of˛. T Two-plot turnover Figure SI1.4: Predicted shapes of five spatial metrics based on fits to empirical shapes of all five metrics, for species Casearia arborea. Casearia arborea has an abundance of 108 individuals in the Barro Colorado Island forest plot and is the species with median abundance in the lowest abundance tertile. Each panel shows five predictions for the named spatial metric, with each prediction based on parameters fitted from the empirical shape of one spatial metric (neighborhood density function: black, Taylor's Law: blue, occupancy area curve: green, two-plot Taylor's Law: red, turnover: orange). For example, in the Taylor's Law panel, the black line shows the predicted Taylor's Law based on fitted parameters from the empirical neighborhood density function, the green line shows the predicted Taylor's Law based on fitted parameters from the empirical occupancy area curve, etc. All fitting and predictions use the parametric forms of the five spatial metrics associated with the theoretical example of a Gaussian second-order intensity function (see main text). Overall, the five predictions for each metric can be seen to be relatively similar to each other compared to variation in best fitting metric shapes across all species (light gray lines). Similar plots for all 171 species analyzed in the Barro Colorado Island dataset are available in Supporting Information 4. T Two-plot turnover Figure SI1.5: Predicted shapes of five spatial metrics based on fits to empirical shapes of all five metrics, for species Chrysochlamys eclipes. Chrysochlamys eclipes has an abundance of 432 individuals in the Barro Colorado Island forest plot and is the species with median abundance in the middle abundance tertile. The interpretation of panels and lines in this figure is identical to that of Figure SI1.4. Overall, the five predictions for each metric can be seen to be relatively similar to each other compared to variation in best fitting metric shapes across all species. Predictions based on empirical fits to Taylor's Law and the two-plot Taylor's Law, however, can be seen to systematically deviate by a small amount from predictions based on fits to the other three metrics. This deviation is consistent with the overall patterns in Figure 5 in the main text, including both the low correlation between fitted parameters from Taylor's Law and several other metrics and the high correlation between fitted parameters from Taylor's Law and the two-plot Taylor's Law. T Two-plot turnover Figure SI1.6: Predicted shapes of five spatial metrics based on fits to empirical shapes of all five metrics, for species Psychotria horizontalis. Psychotria horizontalis has an abundance of 1,856 individuals in the Barro Colorado Island forest plot and is the species with median abundance in the highest abundance tertile. The interpretation of panels and lines in this figure is identical to that of Figure SI1.4. Similar to Figure SI1.4, predictions group into two sets, with those based on empirical fits to Taylor's Law and the two-plot Taylor's Law deviating from predictions based on fits to the other three metrics.