Probabilistic approaches for investigating species co-occurrence from presence-absence maps

Example datasets

We have analysed two real-life datasets using several pairwise approaches to determine the species-by-species co-occurrence patterns of tree species in Panama and the USA. The first dataset consists of the tree locations, recorded in 2015, from Barro Colorado Island (BCI) in Panama; it is downloaded from “https://forestgeo.si.edu”. The dataset contains 248,835 tree locations capturing 296 tree species of Panama inside a study area of 50 hectare (500 × 1000 m²). In this article, we have selected eight most abundant species to investigate their co-occurrences. See Condit (1998), Hubbell et al. (1999), Condit et al. (2019b), and Condit et al. (2019a) for further details about this dataset.

The second dataset is the famous Lansing dataset from the R package spatstat.data. This dataset provides the locations of 2,251 trees and their botanical classification into six types: black oak, hickory, maples, red oak, white oak, and miscellaneous. These data were collected by Gerrard (1969) from a plot in Lansing Woods, Clinton County, Michigan, USA, that measures 924 ft × 924 ft (19.6 acre). The original plot size has been rescaled to the unit square.

Statistical tests

Let W ⊂ R² denote the study area. We divide the study area into a regular grid, and let N denote the total number of grid cells. We are interested in investigating the association between species A and species B. Let N_a be the number of cells in which species A is present and N_b be the number of cells in which species B is present. Under the randomness assumption, the expected number of cells in which both species A and B are present is E₁ = N_aN_b/N, in which neither of them is present is E₂ = (N − N_a)(N − N_b)/N, in which species A is present but species B is absent is E₃ = N_a(N − N_b)/N, and in which species A is absent but species B is present is E₄ = (N − N_a)N_b/N. For each E_i, i = 1, …, 4, let O_i denote the corresponding random variable. Therefore, the random variables O₁, O₂, O₃, and O₄ represent the number of grid cells (out of N cells) containing both species A and B, containing neither of them, containing only species A, and containing only species B, respectively. Furthermore, let o_i, i = 1, …, 4, denote the observed value corresponding to the random variable O_i. A Chi-squared test is used to test the null hypothesis that the occurrences of species A and B are independent, and the corresponding test statistic can be computed as (1)${\displaystyle {\chi }^2=\sum _{i=1}^4\frac{{\left(o_i-E_i\right)}^2}{E_i}.}$

Under the null hypothesis of independence, the test statistic Eq. (1) follows the Chi-squared distribution with 1 degree of freedom (d.f.). If at least one expected frequency E_i is smaller than ten, then the above test statistic is replaced by the Yates’ corrected version (Yates, 1934): (2)${\displaystyle {\chi }_{\text{Yates}}^2=\sum _{i=1}^4\frac{{\left(|o_i-E_i|-0.5\right)}^2}{E_i}.}$

When (|o_i − E_i| − 0.5) is negative in Eq. (2), it is set equal to zero. The main disadvantage of the Chi-squared test approach is the lack of information. When the null hypothesis is rejected, we need some other method to investigate whether the alternative hypothesis of positive (or negative) association between the two species can be supported.

To develop new approaches for testing associations between species, we utilise the following key results. Under the assumption that every individual is randomly distributed in each cell with equal probability, the random variable O_i follows the binomial distribution with N “trials” and “success” probability E_i/N, which is expressed as O_i ∼ Bin(N, E_i/N). If N is large and the “success” probability E_i/N is small, then O_i is approximately distributed as Poisson distribution with expected value E_i, and this is expressed as $O_i\dot{\sim }\text{Poi}\left(E_i\right)$ for i = 1, …, 4. Similarly, if the randomness assumption holds true, for large N and small E_i (i = 1, …, 4), we also have $O_i+O_j\dot{\sim }\text{Poi}\left(E_i+E_j\right)$ for i ≠ j.

We use the random variables O_i, i = 1, …, 4, to construct six new test statistics for testing non-random associations between two species. If the occurrences of species A and species B are positively associated, O₁ and O₂ tend to be large, but O₃ and O₄ tend to be small. In contrast, if they are negatively associated, O₁ and O₂ tend to be small, but O₃ and O₄ tend to be large. For testing positive association, we use the following five rejection rules: (i) O₁ is large, (ii) both O₁ and O₂ are large (iii) O₁ + O₂ is large, (iv) both O₃ and O₄ are small, and (v) O₃ + O₄ is small. Note that the test statistics for (i), (iii), and (v) are scaler-valued, whereas the test statistics for (ii) and (iv) are random vectors, given by (O₁, O₂)^⊤ and (O₃, O₄)^⊤, respectively. We define the test statistics for assessing negative association by replacing the term large with the term small and vice versa in (i)–(v). Because the random variables O₁, …, O₄ follow either binomial or Poisson distributions under the null hypothesis, we show below that it is straightforward to compute the p-values using our proposed test statistics.

For testing the alternative hypothesis of positive association, we propose the following six methods to compute the p-values. Let $C_r^N$ (read N choose r) denote the number of combinations of N things taken r at a time. The quantity $C_r^N$ is computed as N!/(r!(N − r)!), where N! is 1 × 2 × ⋯ × N.

• Under binomial distribution, the test statistic O₁ is used to compute the p-value as follows: $P_1=\mathrm{Pr}\left(O_1\ge o_1|H_0\right)={\sum }_{j\ge o_1}{C}_j^N{\left(\frac{E_1}{N}\right)}^j{\left(1-\frac{E_1}{N}\right)}^{N-j}$.

• Under binomial distribution, the test statistic (O₁, O₂)^⊤ is used to compute the p-value as follows: $P_2=\mathrm{Pr}\left(O_1\ge o_1\text{and}{O}_2\ge o_2|H_0\right)\approx {\prod }_{i=1}^2\mathrm{Pr}\left(O_i\ge o_i|H_0\right)={\prod }_{i=1}^2\left[{\sum }_{j\ge o_i}{C}_j^N{\left(\frac{E_i}{N}\right)}^j{\left(1-\frac{E_i}{N}\right)}^{N-j}\right]$.

• Under binomial distribution, the test statistic (O₃, O₄)^⊤ is used to compute the p-value as follows: $P_3=\mathrm{Pr}\left(O_3\le o_3\text{and}{O}_4\le o_4|H_0\right)\approx {\prod }_{i=3}^4\mathrm{Pr}\left(O_i\le o_i|H_0\right)={\prod }_{i=3}^4\left[{\sum }_{j\le o_i}{C}_j^N{\left(\frac{E_i}{N}\right)}^j{\left(1-\frac{E_i}{N}\right)}^{N-j}\right]$.

• Under Poisson distribution, O₁ is used to compute the p-value as follows: $P_4=\mathrm{Pr}\left(O_1\ge o_1\mid H_0\right)={\sum }_{j\ge o_1}\frac{exp\left(-E_1\right)E_1^j}{j!}$.

• Under Poisson distribution, the test statistic O₁ + O₂ is used to compute the p-value as follows: $P_5=P\left(O_1+O_2\ge o_1+o_2\mid H_0\right)={\sum }_{j\ge o_1+o_2}\frac{exp\left[-\left(E_1+E_2\right)\right]{\left(E_1+E_2\right)}^j}{j!}$.

• Under Poisson distribution, the test statistic O₃ + O₄ is used to compute the p-value as follows: $P_6=P\left(O_3+O_4\le o_3+o_4\mid H_0\right)={\sum }_{j\le o_3+o_4}\frac{exp\left[-\left(E_3+E_4\right)\right]{\left(E_3+E_4\right)}^j}{j!}$.

If the p-value is below the significance level, we conclude that the occurrences of species A and species B are positively associated. These six methods of computing p-values can be used for testing negative association by changing the directions of the inequalities in P₁ − P₆. For example, if $P_1=P_r\left(O_1\le o_1|H_0\right)={\sum }_{j\le o_1}{C}_j^N{\left(\frac{E_1}{N}\right)}^j{\left(1-\frac{E_1}{N}\right)}^{N-j}$ is lower than the significance level then we have enough evidence to support that the occurrences of the two species are negatively associated.

Besides these six approaches to calculating p-values, we also consider the probabilistic approach of Veech (2013), which in our simulation study has been used as a benchmark for comparison with the other approaches. The p-value corresponding to this probabilistic approach is given by ${\displaystyle P_7=P_r\left(O_1\ge o_1|H_0\right)=\sum _{j\ge o_1}\frac{C_j^N\times C_{N_b-j}^{N-j}\times C_{N_a-j}^{N-N_b}}{C_{N_b}^N\times C_{N_a}^N},}$ where max{0, N_a + N_b − N} ≤ j ≤ min{N_a, N_b}, is used as p-value for testing positive association. It is notable that this method is similar to binomial test but not equivalent. The numerator of the equation is the total number of ways that species A and B could be distributed among N cells for a given N_a, N_b and j. The denominator represents the total number of ways that species A and B can be arranged among N cells without regard for j (Veech, 2013). This quantity P₇ can be served as p-value for testing negative association through the same modification used for P₁.

Griffith, Veech & Marsh (2016) developed an R package, cooccur, to implement Veech (2013)’s method. It is highly accessible and handles large datasets with high performance. All computational procedures in this research were implemented in R software. The R source code is provided in the Supplemental Files and the repository at https://github.com/Yamei628/Prob-Cooccur.

Simulation scenarios

We conduct a simulation study to compare the performance of our methods with the method proposed by Veech (2013). Functions runifpoint and rpoispp in spatstat, an R software package developed by Baddeley & Turner (2005), are used to simulate spatial point patterns. The spatial point patterns are simulated in a 1 × 1 plot. The plot is divided into a regular grid where the cell size is 0.05 × 0.05 and the number of cells is 400. Function lets.presab.points in letsR, a R software package developed by Vilela & Villalobos (2015), is used to convert spatial point patterns into presence-absence maps. The test for investigating the association between two species is based on the presence-absence map. Figure 1 illustrates an example of a uniform point pattern and its presence-absence map.

The spatial point patterns are simulated under the following five scenarios.

• Scenario 1: The occurrences of species A and species B are independent. The abundances of these two species are the same and given at various values, n_a = n_b = 20, 40, 100, 200, 300, …, 1000.

• Scenario 2: The occurrences of species A and species B are positively associated. The intensities of these two species are the same and given in various settings. The average presence cells of these two species are similar.

• Scenario 3: The occurrences of species A and species B are positively associated. The intensity of species A is fixed, but the intensity of species B is given in various settings. The average presence cells of these two species are dissimilar.

• Scenario 4: The occurrences of species A and species B are negatively associated. The intensities of these two species are both given in various settings. The average presence cells of these two species are similar.

• Scenario 5: The occurrences of species A and species B are negatively associated. The intensity of species A is fixed, but the intensity of species B is given in various settings. The average presence cells of these two species are dissimilar.

For Scenario 1, the function runifpoint is used to simulate two random point patterns of given abundances. The Chi-squared statistic Eq. (1) (or Eq. (2) when the Yates’ correction is applicable is used to determine whether these two species are independent. For the rest four scenarios, the intensities are given and the occurrences of two species are either positively or negatively associated. Function rpoispp is used to generate point patterns under these scenarios. Let λ_a(x, y) and λ_b(x, y) be the intensity function of species A and species B, respectively, where (x, y) are the spatial coordinates of any location in (0, 1)² ⊂ R². Figure 2 shows the presence-absence maps of species A and species B where λ_a(x, y) = 400x and λ_b(x, y) = 800x. This figure demonstrates an example that the occurrences of two species are positively associated. In contrast, Fig. 3 shows an example that the occurrences of two species are negatively associated, where λ_a(x, y) = 400x and λ_b(x, y) = 400 − 400x.

The level of significance is set at 0.05. For evaluating the proposed methods, the proportion of rejecting the null hypothesis is evaluated based on 500 simulations generated using the true underlying intensity functions.

Simulation results

Table 1 presents the simulation results of the Chi-squared test of random association (i.e., the null hypothesis states that the two point patterns are independent of each other) corresponding to Scenario 1 based on the test statistic defined in Eq. (2). The first two columns of the table present the abundance values n_a and n_b considered for two species A and B, respectively. The next two columns present the average number of cells (${\bar{N}}_a$ and ${\bar{N}}_b$) that are associated with species presence, and the last column presents the estimated type-I error rate based on 500 simulations.

When the abundance is between 200 and 700, the presence rate (the proportion of the number of presence cells to the number of the overall cells) is between 39.3% and 82.5%. It is an ideal range to use the Chi-squared statistic for testing independence since the estimated type-I error rates are close to the nominal significance level of 5%. When the abundances are outside the range of 200 to 700, the estimated type-I error rates are lower than the nominal significance level. These low estimated type-I error rates indicate that the Chi-squared test of independence may produce high type-II error rates (thus, low statistical powers) when the two species are not independent and have presence rates outside the range of 40% to 80%.

The simulation results of Scenario 2–5 are listed in Tables 2, 3, 4 and 5, respectively. The intensities λ_a(x, y) and λ_b(x, y) are listed in the first two columns of the tables when the two intensities are different; otherwise, the common intensity is presented in the first column (e.g., in Table 2). In the last seven columns of Tables 2–5, we report the proportions of null hypothesis rejection corresponding to the seven approaches described in the previous section. Note that for Scenario 2–5, these proportions are estimated statistical powers corresponding to the seven approaches. For the ith approach (i = 1, …, 7), we have used the notation 1 − β_i to denote the corresponding statistical power, where β_i is the estimated type-II error rate.

Tables 2–5 show that the estimated powers of the approaches that use binomial distributions to compute p-values are usually higher than those that use Poisson distributions. In particular, the results show that 1 − β₁ is higher than 1 − β₄ for all the scenarios reported in Tables 2–5. The two approaches use the same random variable O₁ to compute the p-values; however, 1 − β₁ is computed assuming a binomial distribution, whereas 1 − β₄ is computed assuming a Poisson distribution. The second and third approaches, both of which used two variables to compute p-values, performed better than the first approach. Overall, the second approach performed best and produced the highest estimated power in most scenarios.

In the scenarios with positive associations (i.e., Scenario 2 and 3), the approaches based on P₂, P₃, and P₇ performed equally well when the average abundances were large for both species. For example, in Scenario 2, the estimated powers of the three approaches are unity for all the cases with large abundances. However, the second approach performed better than the other two in the cases with low abundances. When the average abundance is 50, Table 2 shows that 1 − β₂ is 0.254, which is higher than 1 − β₃ (0.064) and 1 − β₇ (0.154). Similarly, when the average abundance ${\bar{n}}_b$ is 50 in Table 3, the estimated power of the second approach is 0.452, which is higher than 0.372 and 0.362, estimated powers of the third and seventh approaches, respectively.

In the scenarios with negative associations (i.e., Scenario 4 and 5), the top three performers with high estimated powers are the approaches based on P₂, P₃, and P₇. These three approaches consistently outperformed the other approaches for the cases reported in Tables 4–5. For Scenario 4, the second approach performed best, especially in the extreme cases where there are very low abundances or very high abundances. For example, when the average abundance is considered to be 50, the estimated power based on P₂ is 0.224, which is higher than the estimated powers 0.02 and 0.114 based on P₃ and P₇, respectively. Moreover, when the average abundance is high (approximately 2,500), the estimated power based on P₂ equal to 0.784 (in the last row of Table 4), which is significantly higher than the estimated powers 0.000 and 0.222 based on P₃ and P₇, respectively. For Scenario 5, when large abundance values are considered for species B, the third approach performed slightly better than the second and seventh approaches. For example, in Table 5, the estimated power 1 − β₃ is the highest when the species abundance ${\bar{n}}_b$ is between 500 and 2500. As with the other scenarios, the second approach performed best for Scenario 5 with low abundance values.

We investigated the type-I error rates of the seven approaches using the nominal significance level 0.05. Table 6 shows the estimated type-I error rates, each calculated based on 500 null samples simulated using Scenario 1. The first two rows and the last row of the table shows that the estimated type-I error rates of the second approach (based on P₂) for testing against positive associations are close to the nominal level of 5%. Thus, the second approach will be the most powerful for testing positive association between two species when both presence rates are either below 22.5% or above 90%. When the presence rates are between 22.5% and 90% (i.e., average presence cells are between 90 and 360), the use of the second approach may lead to inflated type-I error rates. In this scenario, we recommend to use either the third or the seventh approach. In the case of assessing negative associations, the estimated type-I error rates of the second approach are close to the nominal level when the presence rates of both species exceed 85% (i.e., presence cells are more than 340). Because estimated powers of the second approach for these presence rates are highest among the approaches considered in this article, we recommend using the second approach to assess negative associations in the case of large presence rates. Although the estimated powers are high for the second approaches when the presence rates are between 22.5% and 85%, the type-I error estimates for this approach are higher than the nominal level of 5%. In contrast, the type-I error estimates of the seventh approach are below the nominal level. Thus, we recommend using the seventh approach for assessing negative associations when both presence rates are between 22.5% and 85%.

Data examples

BCI data analysis

We analyse the popular tropical tree data from Barro Colorado Island (BCI), Panama. The eight most abundant tropical tree species are selected. Our objective is to investigate all 28 pairwise co-occurrence patterns among these eight species. The dataset corresponding to each species is originally stored as a point pattern where each point corresponds to a tree location within the study region. For investigating the co-occurrence pattern of any two given species, we convert the corresponding point patterns into two 5 m × 5 m gridded presence-absence maps, which allowed us to compute important presence/absence rates for performing statistical tests based on the p-values P₁ − P₇. Figure 4 shows the 5m × 5m gridded presence-absence maps of these eight tree species.

Table 7 provides the names and identifiers of the eight species from BCI, along with some key summary measures of them, namely, the total number of trees within the study region (abundance), the total number of cells that are presences (presence cells), and the presence rate (computed as the percentage of total number of cells). For performing pairwise comparisons, we first use the Chi-squared test to test the null hypothesis of independence, and then, if the null hypothesis is rejected, we use the p-values P_i, i =1 , …, 7, to determine whether there is any significant association between two species. Table 8 shows the results of the statistical tests for alternative hypotheses of positive and negative associations. This table also reported the p-values corresponding to the Chi-squared test of independence.

Table 7:

Eight most abundant species in BCI tree data.

The data contains tree locations and the study area is 50-ha (500 × 1000 m²). The data are converted into 5 × 5 m gridded presence-absence maps. The presence cells and the presence rates are computed from the presence-absence maps.

No.	Species	Abundance	Presence cells	Presence rate
1	Hybanthus prunifolius	37123	12053	60.3%
2	Faramea occidentalis	26420	12060	60.3%
3	Desmopsis panamensis	13048	7901	39.5%
4	Trichilia tuberculata	11779	6521	32.6%
5	Alseis blackiana	8700	5737	28.7%
6	Oenocarpus mapora	7947	1623	8.1%
7	Mouriri myrtilloides	7538	5305	26.5%
8	Swartzia simplex	6103	4991	25.0%

DOI: 10.7717/peerj.15907/table-7

Table 8:

Using eight most abundant species in BCI tree data for testing positive association.

The first two columns are the identity numbers (given in Table 7) of two species. The p-values by using the Chi-squared statistic to test the null hypothesis of independence are listed in the third column. The rest columns are the p-values of P_i for testing positive or negative association, i = 1, …, 7.

Species No.		Chi-sq. test p-value	p-values for testing positive association							p-values for testing negative association
A	B		P1	P2	P3	P4	P5	P6	P7	P1	P2	P3	P4	P5	P6	P7
1	2	0.000	0.000	0.000	0.000	0.003	0.000	0.000	0.000	1.000	1.000	1.000	0.997	1.000	1.000	1.000
1	3	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000
1	4	0.000	0.002	0.000	0.000	0.005	0.000	0.001	0.000	0.998	0.994	0.991	0.995	1.000	0.999	1.000
1	5	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000
1	6	0.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000
1	7	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000
1	8	0.002	0.037	0.003	0.002	0.050	0.028	0.041	0.001	0.965	0.889	0.890	0.952	0.973	0.959	0.999
2	3	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000
2	4	0.000	0.001	0.000	0.000	0.003	0.000	0.000	0.000	0.999	0.996	0.993	0.997	1.000	1.000	1.000
2	5	0.242	0.250	0.071	0.064	0.269	0.223	0.243	0.124	0.756	0.545	0.558	0.736	0.780	0.760	0.882
2	6	0.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000
2	7	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000
2	8	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	1.000	1.000	0.999	1.000	1.000	1.000	1.000
3	4	0.000	0.000	0.000	0.000	0.001	0.001	0.001	0.000	1.000	0.988	0.992	0.999	0.999	1.000	1.000
3	5	0.000	0.002	0.000	0.000	0.003	0.006	0.003	0.000	0.998	0.970	0.976	0.997	0.994	0.997	1.000
3	6	0.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000
3	7	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	1.000	0.996	0.998	1.000	1.000	1.000	1.000
3	8	0.000	0.000	0.000	0.000	0.000	0.002	0.001	0.000	1.000	0.985	0.990	1.000	0.998	0.999	1.000
4	5	0.376	0.744	0.483	0.463	0.733	0.691	0.720	0.816	0.264	0.094	0.106	0.274	0.312	0.284	0.193
4	6	0.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000
4	7	0.992	0.501	0.251	0.252	0.500	0.499	0.500	0.502	0.509	0.257	0.255	0.509	0.505	0.504	0.511
4	8	0.024	0.049	0.009	0.016	0.056	0.117	0.078	0.013	0.953	0.783	0.765	0.946	0.885	0.924	0.988
5	6	0.000	1.000	0.997	0.998	1.000	0.999	1.000	1.000	0.000	0.000	0.000	0.000	0.001	0.000	0.000
5	7	0.000	0.000	0.000	0.000	0.000	0.006	0.001	0.000	1.000	0.974	0.986	1.000	0.994	0.999	1.000
5	8	0.000	0.004	0.000	0.002	0.006	0.039	0.014	0.000	0.996	0.914	0.921	0.995	0.962	0.986	1.000
6	7	0.000	1.000	0.992	0.996	1.000	0.997	1.000	1.000	0.000	0.000	0.000	0.000	0.003	0.000	0.000
6	8	0.000	1.000	0.996	0.998	1.000	0.998	1.000	1.000	0.000	0.000	0.000	0.000	0.002	0.000	0.000
7	8	0.000	0.003	0.000	0.001	0.004	0.038	0.012	0.000	0.997	0.919	0.929	0.997	0.963	0.989	1.000

DOI: 10.7717/peerj.15907/table-8

Notes:

The bold font indicates the p-value (by using Chi-squared statistic) which is lower than the 5% significance level.

Table 8 shows that the Chi-squared test of independence is rejected for all but three pairs: (i) Faramea occidentalis and Alseis blackiana (No. 2 and 5), (ii) Trichilia tuberculata and Alseis blackiana (No. 4 and 5), and (iii) Trichilia tuberculata and Mouriri myrtilloides (No. 4 and 7). Moreover, the results indicate that the six out of eight species (No. 1, 2, 3, 4, 7, and 8) are positively associated; all the pairwise comparisons except the one between the species No. 4 and 7 have produced statistically significant p-values at the 5% significance level. The positive association between the species Hybanthus prunifolius (No. 1) and Swartzia simplex (No. 8) is visible in Fig. 4, as they are nearly absent in the middle-left region. The species Alseis blackiana (No. 5) also exhibits positive associations with the four (No. 1, 3, 7, and 8) out of six species in the above mentioned group.

From Fig. 4, we observe that species No. 1 distributes densely on both left-hand and right-hand sides. Species No. 4 and 7 only aggregates on either of the two sides. Consequently, the testing results show No. 1 is positively associated with species No. 4 and 7, but the latter two species are not positively associated. This example demonstrates that we can not guarantee the positive association between two species if they are positively associated with the same species.

Table 7 reveals that the species Oenocarpus mapora (No. 6) possesses a unique feature of a high abundance but low presence. This may be evident from the spatial distribution displayed in Fig. 4F. While it exhibits clustering at the local scale, there appear to be large distances between clusters at the global scale (i.e., within the entire study region). The results in Table 8 show that Oenocarpus mapora is negatively associated with all the other species.

Lansing data analysis

Table 9 shows the identification number (No.), species name, total number of tree locations, total number of presence cells, and presence rate for the six tree species in the Lansing dataset. Figure 5 shows the 0.05 × 0.05 gridded presence-absence maps. We examine all 15 pairwise co-occurrence patterns that are possible among the six species.

Table 9:

Lansing data contains locations of six tree species. The study plot is 19.6-acre (924 × 924 ft²).

The plot has been re-scaled to the unit square. We converted it into 0.05 × 0.05 gridded presence-absence maps. The presence cells and the presence rates are computed from the presence-absence maps.

No.	Species	Abundance	Presence cells	Presence rate
1	black oak	135	93	23.3%
2	hickory	703	285	71.3%
3	maple	514	212	53.0%
4	miscellaneous	105	72	18.0%
5	red oak	346	198	49.5%
6	white oak	448	247	61.8%

DOI: 10.7717/peerj.15907/table-9

The results of the Chi-squared tests, the tests of positive associations, and the tests of negative associations are presented in Table 10. The results of the Chi-squared test in the third column show that the null hypothesis of random association can be rejected for six pairs at the 5% significance level. These pairs are: (i) black oak and maple (p-value: 0.015); (ii) black oak and miscellaneous (p-value <0.001); (iii) hickory and maple (p-value <0.001); (iv) hickory and miscellaneous (p-value: 0.017); (v) hickory and white oak(p-value: 0.006); and (vi) maple and miscellaneous (p-value <0.001).

Based on most of the p-values reported in Table 10, we conclude that, except for the pair maple and miscellaneous, all other pairs mentioned above are negatively associated. Note that all the p-values (shown in columns 4–10 in Table 10) obtained from testing the positive co-occurrence between maple and miscellaneous are smaller than the standard significance level of 0.05, indicating that this pair has a strong positive co-occurrence pattern. The presence-absence maps in Fig. 5 can also be used to informally verify that the pair maple and miscellaneous exhibit similar presence-absence patterns within the study window. Both categories are rare in the top-left part of the study window, while they are relatively abundant in the bottom and top-right parts.

Both the species black oak and hickory are negatively associated with the same two species, namely, maple and miscellaneous, one may suspect that they form a positive co-occurrence pattern. However, based on the Chi-squared test result, we can conclude that the species black oak and hickory are randomly co-occurring (p-value: 0.22).

Table 10:

Using Lansing data for testing positive association.

The first two columns are the identity numbers (given in Table 9) of two species. The p-values by using the Chi-squared statistic to test the null hypothesis of independence are listed in the third column. The rest columns are the p-values of P_i for testing positive or negative association, i = 1, …, 7.

Species No.		Chi-sq. test p-value	p-values for testing positive association							p-values for testing negative association
A	B		P1	P2	P3	P4	P5	P6	P7	P1	P2	P3	P4	P5	P6	P7
1	2	0.215	0.281	0.085	0.067	0.296	0.233	0.286	0.133	0.762	0.562	0.598	0.744	0.790	0.736	0.916
1	3	0.015	0.953	0.829	0.825	0.942	0.937	0.927	0.995	0.065	0.010	0.010	0.078	0.073	0.083	0.010
1	4	0.000	1.000	0.929	0.942	1.000	0.958	0.991	1.000	0.000	0.000	0.001	0.000	0.048	0.011	0.000
1	5	0.472	0.705	0.452	0.459	0.693	0.675	0.683	0.799	0.352	0.140	0.135	0.362	0.351	0.343	0.274
1	6	0.186	0.800	0.593	0.616	0.780	0.804	0.777	0.925	0.244	0.072	0.060	0.262	0.218	0.243	0.115
2	3	0.000	0.984	0.983	0.982	0.955	0.998	0.998	1.000	0.021	0.000	0.000	0.053	0.002	0.003	0.000
2	4	0.017	0.908	0.772	0.791	0.893	0.924	0.858	0.993	0.120	0.022	0.010	0.137	0.089	0.156	0.014
2	5	0.190	0.284	0.062	0.063	0.320	0.209	0.213	0.115	0.750	0.615	0.613	0.709	0.810	0.808	0.922
2	6	0.006	0.896	0.879	0.869	0.826	0.953	0.964	0.998	0.123	0.004	0.007	0.194	0.055	0.042	0.004
3	4	0.000	0.014	0.001	0.000	0.019	0.027	0.028	0.000	0.990	0.920	0.924	0.987	0.977	0.977	1.000
3	5	0.556	0.649	0.424	0.432	0.626	0.671	0.677	0.755	0.394	0.154	0.149	0.412	0.356	0.348	0.312
3	6	0.985	0.515	0.266	0.274	0.508	0.504	0.514	0.533	0.527	0.283	0.274	0.526	0.524	0.515	0.549
4	5	0.868	0.570	0.311	0.322	0.565	0.545	0.555	0.616	0.500	0.248	0.237	0.502	0.483	0.473	0.486
4	6	0.885	0.489	0.242	0.253	0.488	0.477	0.489	0.498	0.574	0.314	0.309	0.572	0.553	0.537	0.607
5	6	0.502	0.656	0.446	0.454	0.628	0.688	0.693	0.781	0.385	0.141	0.136	0.407	0.338	0.332	0.285

DOI: 10.7717/peerj.15907/table-10

Notes:

The bold font indicates the p-value (by using Chi-squared statistic) which is lower than the 5% significance level.

Many of the aforementioned negative associations can also be verified by visually comparing the point patterns or presence-absence maps between two species. For example, Figure 5 illustrates that, while both maple and miscellaneous are sparse in the upper-left and upper-right corners of the study window, the species black oak and hickory are dense in these two corners. The statistical tests become indispensable when other simpler methods, such as visual inspection of point patterns, do not provide a clear picture of the co-occurrence pattern for any given species pair. For example, the patterns of hickory and white oak do not provide any clear indication about their co-occurrence pattern; however, the p-values obtained from testing the negative co-occurrence reveal that hickory and white oak might be negatively associated with each other.