When do we have the power to detect biological interactions in spatial point patterns?

Abstract Uncovering the roles of biotic interactions in assembling and maintaining species‐rich communities remains a major challenge in ecology. In plant communities, interactions between individuals of different species are expected to generate positive or negative spatial interspecific associations over short distances. Recent studies using individual‐based point pattern datasets have concluded that (a) detectable interspecific interactions are generally rare, but (b) are most common in communities with fewer species; and (c) the most abundant species tend to have the highest frequency of interactions. However, it is unclear how the detection of spatial interactions may change with the abundances of each species, or the scale and intensity of interactions. We ask if statistical power is sufficient to explain all three key results. We use a simple two‐species model, assuming no habitat associations, and where the abundances, scale and intensity of interactions are controlled to simulate point pattern data. In combination with an approximation to the variance of the spatial summary statistics that we sample, we investigate the power of current spatial point pattern methods to correctly reject the null model of pairwise species independence. We show the power to detect interactions is positively related to both the abundances of the species tested, and the intensity and scale of interactions, but negatively related to imbalance in abundances. Differences in detection power in combination with the abundance distributions found in natural communities are sufficient to explain all the three key empirical results, even if all pairwise interactions are identical. Critically, many hundreds of individuals of both species may be required to detect even intense interactions, implying current abundance thresholds for including species in the analyses are too low. Sy n thesis. The widespread failure to reject the null model of spatial interspecific independence could be due to low power of the tests rather than any key biological process. Since we do not model habitat associations, our results represent a first step in quantifying sample sizes required to make strong statements about the role of biotic interactions in diverse plant communities. However, power should be factored into analyses and considered when designing empirical studies.


A.2 Summary statistics for bivariate interaction
First assumption is that the expected point count in any set B can be written as an integral where λ i (u) ≥ 0 is called the intensity. For stationary processes λ i (u) ≡ λ i is a constant, and we assume that λ i > 0, so that EN i (B) = λ i |B|.
We define the cross-K function as a function of a distance parameter r > 0 K 12 (r) := λ −1 2 E o1 N 2 (b(o, r)), where b(o, r) is a ball of radius r > 0 centred at the origin, the expectation E o1 is conditional on the joint process having a point of type 1 at the origin o (for stationary processes the exact location does not matter). Heuristically, λ 2 K 12 (r) is the mean abundance of species 2 within distance r of a typical point of species 1. Equivalently we can define K 21 (r), but due to symmetry K 12 = K 21 . The cross-K index is a powerful statistic for testing purposes, but for a more detailed description of spatial interactions we often study the derivative of K 12 , g 12 (r) := K ′ 12 (r) 2πr , known as the cross (or partial) pair correlation function (pcf). The pcf describes the aggregation/segregation of cross species point locations: The probability of having a species 1 point at some small region dx and a species 2 point at some small region dy is given by g 12 (||x − y||)λ 1 λ 2 dxdy. If the processes are independent, g 12 (r) ≡ 1 and K 12 (r) = πr 2 . We say the processes are aggregated if g 12 > 1, and segregated if g 12 < 1, at any particular distance r > 0.
To estimate these quantities several estimators have been proposed, differing in how the observation bias near the borders of W is corrected. We will look at bivariate versions of the globally corrected "Ohser"-type estimators (Illian et al. 2008, p. 230, Ward andFerrandino (1999) and Wiegand et al. 2016) of the form where c(r) is some constant, mainly responsible for scaling and edge correction and possibly depending on r, and f r is some function, for example an indicator function for K 12 and a kernel function for g 12 . Note that while the theoretical K 12 and g 12 are symmetric in the species indices, their estimators are not necessarily so (see e.g. Lotwick and Silverman 1982).

A.3 Covariance of the summary statistics
The key ideas for this Section follow those given by Lotwick and Silverman 1982. Consider two point processes X 1 and X 2 with fixed counts n 1 , n 2 in a bounded window W . We are interested in the covariance between different distances r > 0, s > 0 of estimators of type where f r is symmetric in x − y (this can be extended to non-symmetric functions). Then, if X 1 and X 2 are independent, Cov[T (r), T (s)] = (n 1 n 2 ) −2 [ n 1 (n 1 − 1)n 2 (n 2 − 1)a 1 (r, s) + n 1 n 2 (n 1 − 1)a 2 (r, s) + n 1 n 2 (n 2 − 1)a 3 (r, s) + a 4 (r, s) − (n 1 n 2 ) 2 a 5 (r, s)] where with expectations over random locations on W such that the pair x, x ′ follow the distribution of X 1 and the pair y, y ′ follow the distribution of X 2 . Heuristically, for the K 12 which effectively counts point pairs, the terms reflect the probabilities of all different types of point pair occurrences, and then get multiplied by corresponding number of possible combinations. For example, a 1 counts the occurrences of two separate pairs which both have a point from each species, and a 2 and a 3 count the times an individual has two neighbours both of the opposite species. These terms are all influenced by the edge of the observation window W , and are therefore sensitive to the size of the neighbourhood radii r, s. In the case of i.i.d. uniform locations (that is, assuming X 1 and X 2 follow Poisson processes and we condition on fixed n 1 and n 2 ), we have a 1 = a 5 and a 2 = a 3 , and the covariance, say Cov 0 (r, s), simplifies to Cov 0 (r, s) = (n 1 n 2 ) −1 [ (n 1 + n 2 − 2)a 2 (r, s) + a 4 (r, s) − (n 1 + n 2 − 1)a 1 (r, s)] = (n 1 n 2 ) −1 [ (n 1 + n 2 )c 2 (r, s) + c 3 (r, s)] where we have further simplified using c 2 = a 2 (r, s) − a 1 (r, s) c 3 = a 4 (r, s) + a 1 (r, s) − a 2 (r, s) to highlight the effect of sample sizes. The additional constant for the statistic M is c 1 (r, s) = c(r)c(s), which mainly scales away biases due to edge effects. This is the form we provide in Equation 2 and use for our examples.
For the general situation where the processes have some pair correlation functions g 1 and g 2 , additional terms appear such that From these expression we can see that if the marginal processes tend towards internal clustering, i.e. g 1 ≥ 1 and g 2 ≥ 1 the extra terms will be positive and the covariances will be higher than with just uniformly random marginals. This means that in our examples where the species 2 is slightly clustered (cf. Section B), our plug-in variances are slightly lower than in truth, leading to over-estimation of power. This explains the "optimistic" bias we see in e.g. Figures 1 and S8. Also noteworthy is that for processes tending towards internal segregation (g 1 ≤ 1, g 2 ≤ 1) the extra terms are negative and the covariances and variances are lower than for uniformly random marginals. From a statistical point of view the best patterns to analyse are strongly internally segregated. The evaluation of Cov 0 (r, s) is not trivial, as all the terms are geometrical integrals determined by the weighting function f and the observation window W . For estimator of the pair correlation function g 12 , consider the box-kernel The quantities I 3 and I 4 can be approximated numerically using Monte Carlo integration, and I 1 is up to a constant the isotropised set covariance of W which has a closed form solution for some elementary shapes of W (rectangle, disc; Illian et al., 2008, p. 485).

A.4 Gaussian approximation ofK 12
The mathematics of the limiting behaviour of the "Ohser"-type estimators are beyond this study, and for progress in this regard we refer to Heinrich (2015). We resort to the same argument as Wiegand et al. (2016): the empirical plots do not show signs against normality apart from very short distances due to the positivity constraint. Fig. S5 illustrates this (compare to Wiegand et al. 2016, Fig. S3). Note the accuracy of the analytical formula derived in Appendix A.3 for the variance.

B Model generated data
The process is inspired by the shot-noise product Cox processes (Jalilian et al., 2015), and is constructed hierarchically. First, let X 1 be a stationary Poisson process with intensity λ 1 . Then conditional on a realisation x 1 of X 1 , let X 2 be an inhomogeneous Poisson process with intensity function λ 2 (u; x 1 ) = e a x∈x 1 (1 + bh(u − x)) ∈ u ∈ W where h(v) = k τ (v)/k τ (0) with k τ a 2D kernel function (probability density) with standard deviation τ > 0, and a ∈ R, b > −1 are parameters controlling the intensity and the interaction, respectively. The joint model is stationary, and isotropic if k is isotropic, and has λ 2 = exp (a + λ 1 b/k τ (0)) , g 11 (u) ≡ 1, g 12 (u) = 1+bh(u), g 22 (u) = exp λ 1 b 2 (h * h)(u) where * denotes convolution. From these properties we see that if −1 < b < 0 the two species exhibit segregation (g 12 < 1), and if b > 0 the two species exhibit aggregation or clustering (g 12 > 1), and when b = 0 the two species are independent. We also see that the both types of interactions result in clustering of species 2 (g 22 ≥ 1). The range of interaction (if defined via the pair correlation) is controlled by the parameter τ . In our examples we use a Gaussian kernel, for which the range, i.e. h is non-zero, is approximately 2τ . Fig. S6 shows two examples of the process with identical type 1 patterns, together with their K 12 and g 12 estimates.  Figure S7: Power of K 12 -based pointwise cross-species independence tests with varying degrees of imbalance n 1 = 30 ≤ n 2 . Range of interaction 2τ = 10. The true power is estimated using 5000 repeated tests with 199 random shifts each. Fig. S8 provides evidence that the analytical power formula works also for the crosspair correlation function g 12 . Compare Fig. S8 to Fig. 1. The optimistic bias in the power is again due to the slight downwards bias in the variance as species 2 is clustered.

C Additional power estimates
The estimation was carried out in this example with the bandwidths 0.15 |W | n 1 as the samples sizes were balanced.  Figure S8: Power of g 12 -based pointwise cross-species independence tests. Range of interaction 2τ = 10. The true power is estimated using 5000 repeated tests with 199 random shifts each.

C.1 Testing without simulation
Note that the approximate Gaussianity and the approximate variance formula lead directly to a χ 2 -test of independence without random shift simulations, much like in the work by Wiegand et al. (2016). Procedure: 1. EstimateK 12 (r) for one r

C.2 Pointwise test vs testing over a range
In the simulation experiments we control all factors, so we can choose the distance of the pointwise test to be optimal, i.e. the distance which we know the power is highest. Table  S1 compares this optimal pointwise power to the power of a test where instead of a single distance an interval of distances is tested simultaneously using a deviation test (see e.g. Myllymäki et al., 2017).  Table S1: Power comparison of the Studentised L 2 deviation test over two distance intervals ("r=1-10" and "r=1-20") with the pointwise power formula at the known optimal distance ("pw.o."). The K 12 statistic, and the deviation test powers estimated using 1000 simulations per model and/or setting as indicated, with 199 random shifts each.
We tried using the covariance formula to combine several distances to a χ 2 -test, but the very short distance asymmetry and the non-central χ 2 did not immediately lead to a useful power approximation of the Studentised L 2 test.

C.3 Improving power by combining summaries
A simple way to improve power is to combine several summaries in the test statistic. As an example, we combined the K 12 with the nearest neighbour distance distribution function D 12 (Van Lieshout and Baddeley, 1999) by using the pointwise test statistic Fig. S9 depicts the pointwise powers for K 12 , D 12 and the combination when the data was generated by our bivariate model with b = 0.5, 2τ = 10. The nearest neighbour summary operates only at short distances as it saturates to 1 quickly, and for distances > 5 is inferior to K 12 in this scenario. But as it captures different information than the K 12 , combining it with K 12 increases the power, at least when r < 10. After r > 10 the combined pointwise power is diminished as the nearest neighbour summary provides no help yet is weighted equally with K 12 in making the decision. Weighting the statistics by their useful ranges is therefore recommended.  Figure S9: Power of K 12 , D 12 and K 12 + D 12 -based pointwise cross-species independence tests when species are significantly segregated. Table S2 gives the powers of a test with T = 20 r=1 T 2 KD (r), over the distances 1 − 20. The power with the combined statistic is higher than with either of the components alone for small samples.

C.4 Sample size and observation window
Under most circumstances samples sizes can only really be increased by increasing the area of observation, and the connection λ i ≈ n i /Area can be used to get a rough idea of the requirements. First, a pilot study needs to be conducted to estimate λ i (see Illian et al. 2008 for estimation techniques). Then we need to determine the minimum Area, accounting for imbalance between species if that is needed. Table S3 gives some example calculations when a square area is used (note that the window geometry might affect the power; see C.5).

C.5 Additional factors
The geometry of the area has an effect on the estimator's variance and hence the power of the test, but according to the analytical formula the effect is relatively small. For example, if we change from a square shape to an elongated rectangle shape with equal area but width-to-height -ratio 3, and consider interaction b = 0.25 and type I error level α = 5%, the power drops from 33.2% to 32.9% with 2τ = 10 and 78.0% to 76.8% with 2τ = 20 for sample size (80, 80), and from 15.4% to 15.3% with 2τ = 10 and 41.4% to 40.5% with 2τ = 20 for sample size (30, 80). Increasing the type I error level α increases the power as illustrated in Table S4. From the table we can see that a 5% increase in α can reduce the type II error β = 1 − power by more than 10%. So in scenarios where we can tolerate some extra false positive discoveries with the simultaneous decrease in false negatives, for example when pre-screening a large data set for more involved downstream analysis on found interacting pairs, adjustments to α should be considered.