Jaccard/Tanimoto similarity test and estimation methods

TL;DR

This paper introduces new statistical hypothesis testing methods for the Jaccard/Tanimoto similarity coefficient applied to binary presence-absence data, enabling more accurate and faster analysis of species co-occurrences and other binary datasets.

Contribution

The paper presents unbiased estimation, exact and asymptotic solutions, and efficient algorithms for significance testing of the Jaccard/Tanimoto coefficient, with implementation in an open-source R package.

Findings

Methods produce accurate p-values and false discovery rates

Estimation algorithms are significantly faster than exact solutions

Applicable to diverse binary data in biology and other sciences

Abstract

Binary data are used in a broad area of biological sciences. Using binary presence-absence data, we can evaluate species co-occurrences that help elucidate relationships among organisms and environments. To summarize similarity between occurrences of species, we routinely use the Jaccard/Tanimoto coefficient, which is the ratio of their intersection to their union. It is natural, then, to identify statistically significant Jaccard/Tanimoto coefficients, which suggest non-random co-occurrences of species. However, statistical hypothesis testing using this similarity coefficient has been seldom used or studied. We introduce a hypothesis test for similarity for biological presence-absence data, using the Jaccard/Tanimoto coefficient. Several key improvements are presented including unbiased estimation of expectation and centered Jaccard/Tanimoto coefficients, that account for occurrence probabilities. We derived the exact and asymptotic solutions and developed the bootstrap and measurement concentration algorithms to compute statistical significance of binary similarity. Comprehensive simulation studies demonstrate that our proposed methods produce accurate p-values and false discovery rates. The proposed estimation methods are orders of magnitude faster than the exact solution. The proposed methods are implemented in an open source R package called jaccard (https://cran.r-project.org/package=jaccard). We introduce a suite of statistical methods for the Jaccard/Tanimoto similarity coefficient, that enable straightforward incorporation of probabilistic measures in analysis for species co-occurrences. Due to their generality, the proposed methods and implementations are applicable to a wide range of binary data arising from genomics, biochemistry, and other areas of science.

Equations

$$T(\text{mbox}{\mathbf{y}}_i,\text{mbox}{\mathbf{y}}_j)=\{\begin{array}{ll}\frac{\text{mbox}{\mathbf{y}}_i\cap \text{mbox}{\mathbf{y}}_j}{\text{mbox}{\mathbf{y}}_i\cup \text{mbox}{\mathbf{y}}_j}& \text{if }\text{mbox}{\mathbf{y}}_i\cup \text{mbox}{\mathbf{y}}_j\ne 0\\ \frac{p_i{p}_j}{p_i+p_j-p_i{p}_j}& \text{otherwise.}\end{array}$$

$$\mathbb{E}(T(\text{mbox}{\mathbf{y}}_i,\text{mbox}{\mathbf{y}}_j))=\frac{p_i{p}_j}{p_i+p_j-p_i{p}_j}.$$

$$\mathbb{E}(T(\text{mbox}{\mathbf{y}}_i,\text{mbox}{\mathbf{y}}_j)|N_1+N_2+N_3)$$

$$=\frac{\mathbb{E}(N_1|N_1+N_2+N_3)}{N_1+N_2+N_3}$$

$$=\frac{\frac{p_i{p}_j}{p_i+p_j-p_i{p}_j}(N_1+N_2+N_3)}{N_1+N_2+N_3}$$

$$=\frac{p_i{p}_j}{p_i+p_j-p_i{p}_j}$$

$$\mathbb{E}(T(\text{mbox}{\mathbf{y}}_i,\text{mbox}{\mathbf{y}}_j)|N_1+N_2+N_3)=\frac{p_i{p}_j}{p_i+p_j-p_i{p}_j}$$

$$\mathbb{E}(T(\text{mbox}{\mathbf{y}}_i,\text{mbox}{\mathbf{y}}_j))$$

$$=\frac{p_i{p}_j}{p_i+p_j-p_i{p}_j}\mathbb{P}(N_1+N_2+N_3=0)$$

$$+\frac{p_i{p}_j}{p_i+p_j-p_i{p}_j}\mathbb{P}(N_1+N_2+N_3>0)$$

$$=\frac{p_i{p}_j}{p_i+p_j-p_i{p}_j}.$$

$$T^c(\text{mbox}{\mathbf{y}}_i,\text{mbox}{\mathbf{y}}_j)=T(\text{mbox}{\mathbf{y}}_i,\text{mbox}{\mathbf{y}}_j)-\mathbb{E}[T(\text{mbox}{\mathbf{y}}_i,\text{mbox}{\mathbf{y}}_j)]$$

$$\begin{array}{rl}{H}_0& :T^c(\text{mbox}{\mathbf{y}}_i,\text{mbox}{\mathbf{y}}_j)=0\\ H_1& :T^c(\text{mbox}{\mathbf{y}}_i,\text{mbox}{\mathbf{y}}_j)\ne 0.\end{array}$$

$$T(\text{mbox}{\mathbf{y}}_i,\text{mbox}{\mathbf{y}}_j)=\{\begin{array}{ll}\frac{N_1}{N_1+N_2+N_3}& \text{if }{N}_1+N_2+N_3>0\\ \frac{p_i{p}_j}{p_i+p_j-p_i{p}_j}& \text{otherwise.}\end{array}$$

$$K_{T^c}=\left\{(N_1,N_2,N_3,N_4):\left|\frac{N_1}{N_1+N_2+N_3}-\mathbb{E}[T(y_i,y_j)]\right|\ge |T^c|\right\}.$$

$$K_{{\hat{T}}^c}=\left\{(N_1,N_2,N_3,N_4):\left|\frac{N_1}{N_1+N_2+N_3}-\frac{{\tilde{p}}_i{\tilde{p}}_j}{{\tilde{p}}_i+{\tilde{p}}_j-{\tilde{p}}_i{\tilde{p}}_j}\right|\ge |{\hat{T}}^c|\right\},$$

$$\sqrt{m}{T}^c(\text{mbox}{\mathbf{y}}_i,\text{mbox}{\mathbf{y}}_j)\to \mathcal{N}(0,{\sigma }^2)$$

$${\sigma }^2=\frac{q_1{q}_2(1-q_2)}{(q_1+q_2{)}^3}.$$

$$\sqrt{m}\left((N_1,N_2+N_3)/m-(q_1,q_2)\right)\to \mathcal{N}(0,\Sigma )$$

$$\Sigma =\left[\begin{array}{cc}{q}_1(1-q_1)& -q_1{q}_2\\ -q_1{q}_2& q_2(1-q_2)\end{array}\right].$$

$$\sqrt{m}\left(T(\text{mbox}{\mathbf{y}}_i,\text{mbox}{\mathbf{y}}_j)-\frac{q_1}{q_1+q_2}\right)\to \mathcal{N}(0,\nabla g(q_1,q_2)\Sigma \nabla g(q_1,q_2{)}^T).$$

$$\nabla g(x_1,x_2)=\left[\frac{x_2}{(x_1+x_2{)}^2},\frac{-x_1}{(x_1+x_2{)}^2}\right].$$

$$\nabla g(q_1,q_2)\Sigma \nabla g(q_1,q_2{)}^T=\frac{q_1{q}_2(1-q_2)}{(q_1+q_2{)}^3}.$$

$$2\varphi \left(\frac{\sqrt{m}}{\sigma }\left(T(\text{mbox}{\mathbf{y}}_i,\text{mbox}{\mathbf{y}}_j)-\frac{{\hat{q}}_1}{{\hat{q}}_1+{\hat{q}}_2}\right)\right)-1,$$

$$p_{\epsilon }^L({\hat{T}}^c)=\mathbb{P}\left(K_{{\hat{T}}^c}\cap I_{\epsilon }\right)\le \mathbb{P}\left(K_{{\hat{T}}^c}\right)\le \mathbb{P}\left(K_{{\hat{T}}^c}\cap I_{\epsilon }\right)+\epsilon =p_{\epsilon }^U({\hat{T}}^c).$$

$$p^L({\hat{T}}^c)=$$

$$p\text{-}value=\frac{1\{|t_b^{\ast }|\ge |t|;b=1,\dots ,B\}}{B}.$$

Full text

\startlocaldefs\endlocaldefs

Jaccard/Tanimoto similarity test

and estimation methods

Neo Christopher Chung${}^{\textnormal{1,*}}$, Błażej Miasojedow${}^{\textnormal{2}}$, Michał Startek${}^{\textnormal{1}}$, Anna Gambin${}^{\textnormal{1}}$

${}^{\textnormal{1}}$Institute of Informatics, University of Warsaw

${}^{\textnormal{2}}$Institute of Mathematics, Polish Academy of Sciences

${}^{\textnormal{*}}$nchchung@gmail.com

A survey of presences and absences of specific species across multiple biogeographic units (or bioregions) are used in a broad area of biological studies from ecology to microbiology. Using binary presence-absence data, we evaluate species co-occurrences that help elucidate relationships among organisms and environments. To summarize similarity between occurrences of species, we routinely use the Jaccard/Tanimoto coefficient, which is the ratio of their intersection to their union. It is natural, then, to identify statistically significant Jaccard/Tanimoto coefficients, which suggest non-random co-occurrences of species. However, statistical hypothesis testing using this similarity coefficient has been seldom used or studied.

We introduce a hypothesis test for similarity for biological presence-absence data, using the Jaccard/Tanimoto coefficient. Several key improvements are presented including unbiased estimation of expectation and centered Jaccard/Tanimoto coefficients, that account for occurrence probabilities. The exact and asymptotic solutions are derived. To overcome a computational burden due to high-dimensionality, we propose the bootstrap and measurement concentration algorithms to efficiently estimate statistical significance of binary similarity. Simulation studies demonstrate that our proposed methods produce accurate p-values and false discovery rates. The proposed estimation methods are orders of magnitude faster than the exact solution, particularly with an increasing dimensionality. We showcase their application in evaluating co-occurrences of bird species in 28 islands of Vanuatu and fish species in 3347 freshwater habitats in France. The proposed methods are implemented in a R package called jaccard (https://cran.r-project.org/package=jaccard).

We introduce a suite of statistical methods for the Jaccard/Tanimoto similarity coefficient, that enable straightforward incorporation of probabilistic measures in analysis for species co-occurrences. Due to their generality, the proposed methods and implementations are applicable to a wide range of binary data arising from genomics, biochemistry, and other areas of science.

Keyword: Jaccard, Tanimoto, binary similarity, hypothesis test, co-occurrences, p-value

Funding: Narodowe Centrum Nauki 2014/12/W/ST5/00592 and 2016/23/D/ST6/03613

Background

Analysis of species co-occurrences helps us understand ecological and biological relationships among species. Essentially, the presence (1) and absence (0) of species are surveyed in multiple biogeographic units (or bioregions) using fieldwork, imaging, sequencing, and other techniques. Then, the Jaccard/Tanimoto coefficient is one of the most fundamental and popular similarity measures to compare such biological presence-absence data. Given two presence-absence vectors $\mbox{\bf{y}}_{i}$ and $\mbox{\bf{y}}_{j}$ of length $m$ that represent two different species, the Jaccard/Tanimoto similarity coefficient is the ratio of their intersection to their union, $T(\mbox{\bf{y}}_{i},\mbox{\bf{y}}_{j})=\mbox{\bf{y}}_{i}\cap\mbox{\bf{y}}_{j}/\mbox{\bf{y}}_{i}\cup\mbox{\bf{y}}_{j}$ [1, 2]. This quantification of overlaps allows us to quantify co-existence of species [3, 4, 5, 6]. However, the Jaccard/Tanimoto coefficient lacks probabilistic interpretations or statistical error controls. Surprisingly, its statistical properties, hypothesis testing, and estimation methods for p-values have been inadequately studied. Here, we present a rigorous statistical test evaluating the similarity in presence-absence data, derive exact and asymptotic solutions, and introduce efficient estimation methods for significance of the Jaccard/Tanimoto similarity coefficient.

Generally, analysis of co-occurrences enables us to distinguish generalist species that survive in a broad range of environments from specialists that only thrive in a few localities [7, 8]. Alternatively, similarity between two localities – how two biogeographic units share an overlapping set of species – sheds light on the beta diversity that may arise from ecological processes over time [9, 10, 11]. There has been a long standing discussion on how to conduct association analysis for occurrences of species, including appropriate null models and evaluation techniques [12, 13, 14, 15, 16, 17]. There are also specialized probabilistic approaches, including metrics related to the Jaccard/Tanimoto coefficient [18, 19, 20, 21]. Yet, these studies rarely utilized statistical significance. Therefore, we investigated a hypothesis test using the Jaccard/Tanimoto coefficient that underlies or accompanies most of such association analyses.

The Jaccard/Tanimoto coefficient measuring similarity between two species has long been used to evaluate co-occurrences between species or between biogeographic units [22, 23, 3, 4, 5, 24]. Pioneering early works on probabilistic treatment of the Jaccard/Tanimoto coefficient assume that the probability of species occurrences is 0.5 [22, 23, 5]. These can be seen as special cases of our methods where both probabilities of $\mbox{\bf{y}}_{i}$ and $\mbox{\bf{y}}_{j}$ are set to 0.5. Recently, [24] and [25] proposed estimating p-values with combinatorics and hypergeometric distributions, respectively. We found that they are inaccurate. To provide a comprehensive statistical treatment, we have developed a suite of methods and estimation techniques for rigorously testing similarity between presence-absence data.

We derive a hypothesis test from the first principles using the Jaccard/Tanimoto coefficient. In the process, we propose an unbiased estimation of expectation and a centered Jaccard/Tanimoto coefficient that accounts for different probabilities of species occurrences. The negative and positive values of the centered Jaccard/Tanimoto coefficient naturally correspond to negative and positive association. We introduce an exact distribution of Jaccard/Tanimoto similarity coefficients under independence that is shown to provide accurate p-values. Because the exact solution for a large $m$ is computationally expensive, we have developed two efficient and accurate estimation algorithms. We demonstrate their remarkable accuracy and computational efficiency in comprehensive simulation studies, where p-values and false discovery rates (FDRs) are evaluated. As applications, we evaluated co-occurrences of bird species from $m=28$ islands of Vanuatu and of fish species from $m=3347$ freshwater habitats in France.

All proposed methods are implemented in a statistical programming language R [26], available on the Comprehensive R Archive Network (https://cran.r-project.org/package=jaccard). We additionally provide an interactive web app (https://nnnn.shinyapps.io/jaccard). The implementations are efficient and general, such that the jaccard package can rigorously test similarity between binary data arising from genomics, biochemistry, and others.

Methods

Statistical Model and Test

Quantitative comparison of presence-absence data in ecology and biology plays a crucial role in evaluating species co-existences, biodiversities, and ecosystems. In particular, one may be interested in comparing how species are co-occurring in biogeographic units or how biogeographic units are occupied by certain species. Note that species are used generally to indicate groups of organisms under investigations, such as operational taxonomic units (OTUs); similarly, biogeographic units or bioregions could be distinct survey areas, islands, or habitats. We are interested in statistically testing similarity between a pair of presence-absence data.

Given two presence-absence vectors $\mbox{\bf{y}}_{i}$ and $\mbox{\bf{y}}_{j}$ of length $m$, we are interested in inferring whether they are significantly related. Consider presence (1) and absence (0) of two species are recorded at $m$ biogeographic units. We measure their similarity by the ratio of their intersection to their union, $T(\mbox{\bf{y}}_{i},\mbox{\bf{y}}_{j})=\mbox{\bf{y}}_{i}\cap\mbox{\bf{y}}_{j}/\mbox{\bf{y}}_{i}\cup\mbox{\bf{y}}_{j}$. This is well known as the Jaccard/Tanimoto index or similarity coefficient [1, 2]. In order to utilize the Jaccard/Tanimoto similarity coefficient in a statistically rigorous manner, we propose a family of methods and algorithms (Figure 1).

Under the null model of independence, $\mbox{\bf{y}}_{i}$ and $\mbox{\bf{y}}_{j}$ are assumed to be independent and identically distributed (i.i.d.). They are modeled by a Bernoulli distribution, with corresponding occurrence (i.e., success) probabilities $p_{i}$ and $p_{j}$ $\in[0,1]$. Specifically, for $k=1,\ldots,m$, $\mbox{\bf{y}}_{i,k}\sim_{i.i.d.}\textnormal{Bernoulli}(p_{i})$ and $\mbox{\bf{y}}_{j,k}\sim_{i.i.d.}\textnormal{Bernoulli}(p_{j})$. Because this conventional definition is undefined if both binary vectors contain only zeros such that $\mbox{\bf{y}}_{i}\cup\mbox{\bf{y}}_{j}=0$, we refine the definition of Jaccard/Tanimoto coefficient

[TABLE]

Following the definition of Jaccard/Tanimoto similarity coefficient in Eq. (1), we derive its expected value $\mathbb{E}[T(\mbox{\bf{y}}_{i},\mbox{\bf{y}}_{j})]=\frac{p_{i}p_{j}}{p_{i}+p_{j}-p_{i}p_{j}}$. Substantial deviation from the expected value signifies similarity. Note that the Jaccard/Tanimoto coefficient can also be defined in terms of a multinomial distribution with four categories and $m$ trials (for example, representing $m$ biogeographic units). Four categories arising from presence-absence data are $N_{1}=\mbox{\bf{y}}_{i}\cap\mbox{\bf{y}}_{j}$, $N_{2}=\mbox{\bf{y}}_{i}\cap(1-\mbox{\bf{y}}_{j})$, $N_{3}=(1-\mbox{\bf{y}}_{i})\cap\mbox{\bf{y}}_{j}$ and $N_{4}=m-N_{1}-N_{2}-N_{3}$. From $p_{i}$ and $p_{j}$, probabilities of those four categories are $p_{i}p_{j}$, $p_{i}(1-p_{j})$, $(1-p_{i})p_{j}$ and $(1-p_{i})(1-p_{j})$, respectively. Putting them together, $\boldsymbol{N}=(N_{1},N_{2},N_{3},N_{4})$ is distributed according to a multinomial distribution, $\textnormal{Multi}(m,p_{i}p_{j},p_{i}(1-p_{j}),(1-p_{i})p_{j},(1-p_{i})(1-p_{j}))$.

Proposition 1**.**

If $\mbox{\bf{y}}_{i}$ and $\mbox{\bf{y}}_{j}$ are independent, then

[TABLE]

Proof.

First, we compute conditional expectation given $N_{1}+N_{2}+N_{3}$. We observe that $N_{1}|N_{1}+N_{2}+N_{3}$ follows $\textnormal{Bernoulli}(N_{1}+N_{2}+N_{3},\frac{p_{i}p_{j}}{p_{i}+p_{j}-p_{i}p_{j}})$. Hence, on set $N_{1}+N_{2}+N_{3}>0$, we have

[TABLE]

and on set $N_{1}+N_{2}+N_{3}=0$, we have

[TABLE]

Therefore,

[TABLE]

∎

This allows us to define the centered Jaccard/Tanimoto coefficient as

[TABLE]

This accounts for expected values, naturally distinguishing negative and positive associations. Generally, we would like to measure the deviation of an observed coefficient from an expected value, instead of simply looking at a magnitude of an observed statistics. Furthermore, this centered coefficient may be scaled by variance in order to span a pre-defined range.

To evaluate whether $\mbox{\bf{y}}_{i}$ and $\mbox{\bf{y}}_{j}$ are independent, a following statistical hypothesis testing is performed:

[TABLE]

The null hypothesis $H_{0}$ is that the centered Jaccard/Tanimoto coefficient equals zero. Note that this is equivalent to that the conventional (uncentered) Jaccard/Tanimoto coefficient equals an expected value under independence. Therefore, although we propose and use the centered coefficient, this hypothesis testing is attributed to both uncentered and centered versions. Then, a p-value indicates a probability of observing a coefficient equal to or more extreme than an observed coefficient under the null hypothesis.

Distribution of the Jaccard/Tanimoto Coefficient

To obtain its p-value, we derive the distribution of Jaccard/Tanimoto coefficient under the null hypothesis. In terms of $\boldsymbol{N}=(N_{1},N_{2},N_{3},N_{4})$, the Jaccard/Tanimoto coefficient can be expressed as

[TABLE]

When $p_{i}$ and $p_{j}$ are known, the p-value is given by $\mathbb{P}(K_{T^{c}})$ where

[TABLE]

However, in practice, probabilities $p_{i}$ and $p_{j}$ are usually unknown. Therefore, we define the centered Jaccard/Tanimoto coefficient by $\hat{T}^{c}=T-\frac{\hat{p}_{i}\hat{p}_{j}}{\hat{p}_{i}+\hat{p}_{j}-\hat{p}_{i}\hat{p}_{j}}$, where $\hat{p}_{i}=\frac{\sum{\mbox{\bf{y}}_{i}}}{m}$, $\hat{p}_{j}=\frac{\sum{\mbox{\bf{y}}_{j}}}{m}$ are standard estimators of $p_{i}$ and $p_{j}$ respectively. Plug-in estimates of $\mathbb{E}[T(y_{i},y_{j})]$ into Eq. (4) will result in conservative behaviors, since we estimate the probabilities on the same sample that we want to perform the test. Then, the estimates of expectation are biased toward the observed value of Jaccard/Tanimoto coefficient. To overcome this bias, we estimate probabilities $p_{i}$ and $p_{j}$ for each configuration $(N_{1},N_{2},N_{3},N_{4})$ separately.

So in this case, the critical region is defined as follows

[TABLE]

where $\tilde{p}_{i}=\frac{N_{1}+N_{2}}{m}$ and $\tilde{p}_{j}=\frac{N_{1}+N_{3}}{m}$.

Because the exact distribution is computationally expensive (see Results for comparison), we introduce an asymptotic approximation when $m\to\infty$. It may be useful when dealing with very large binary data, where computational power is a bottleneck. Denote by $q_{1}=p_{i}p_{j}$ the probability that both $\mbox{\bf{y}}_{i}$ and $\mbox{\bf{y}}_{j}$ have ones, and by $q_{2}=p_{i}+p_{j}-2p_{i}p_{j}$ the probability that only one of two vectors has one. Similarly, $\hat{q}_{1}$ and $\hat{q}_{2}$ are defined with the plug-in estimators. As $m\to\infty$, we can estimate the variance:

Proposition 2**.**

If $\mbox{\bf{y}}_{i}$ and $\mbox{\bf{y}}_{j}$ are independent then

[TABLE]

as $m\to\infty$, where

[TABLE]

Proof.

Theorem 14.6 of [27] states that

[TABLE]

where

[TABLE]

Then, we define function $g(x_{1},x_{2})=\frac{x_{1}}{x_{1}+x_{2}}$ and apply the delta method. So, we get

[TABLE]

The gradient of $g$ is

[TABLE]

Finally, after simplification, we obtain

[TABLE]

∎

In practice, probabilities $p_{i}$ and $p_{j}$ are unknown and need to be estimated. Recall that $\hat{p}_{i}=\frac{\#\{y_{ik}=1\}}{m}$ and $\hat{p}_{j}=\frac{\#\{y_{jk}=1\}}{m}$. We define $\hat{q}_{1}$ and $\hat{q}_{2}$ by replacing in definition of $q_{1}$ and $q_{2}$ true probabilities $p_{i}$ and $p_{j}$ by its estimators. So based on Proposition 2 we are able to approximate p-values as follow:

[TABLE]

where $\phi=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}e^{-x^{2}/2}dx$ is a standard Gaussian cumulative distribution function (CDF).

Measure Concentration Algorithm

The distribution of the centered Jaccard/Tanimoto coefficient can be expressed in terms of the multinomial distribution. However, evaluating a significance test based on this representation requires exhaustive computations. It needs summation over all possible states of the multinomial distribution. For the centered Jaccard/Tanimoto coefficient between $y_{i}$ and $y_{j}$, we need to compute probability of event $K_{\hat{T}^{c}}$ defined by Eq. (5).

This can be quickly and accurately estimated by the measure concentration algorithm (MCA) with a known error bound [28]. For every $\varepsilon>0$, we will construct $I_{\varepsilon}$, a set of $(N_{1},N_{2},N_{3},N_{4})$ with $N_{1}+N_{2}+N_{3}+N_{4}=m$, such that $\mathbb{P}(N_{1},N_{2},N_{3},N_{4})\in I_{\varepsilon}\geq 1-\varepsilon$. Given the set $I_{\varepsilon}$, we have following bounds

[TABLE]

In addition, $p^{U}_{\varepsilon}(\hat{T}^{c})-p^{L}_{\varepsilon}(\hat{T}^{c})=\varepsilon$.

The idea behind the algorithm is that a multinomial distribution concentrates around its mode. Two possible states $\boldsymbol{N}=(N_{1},N_{2},N_{3},N_{4})$ and $\boldsymbol{N}^{\prime}=(N_{1}^{\prime},N_{2}^{\prime},N_{3}^{\prime},N_{4}^{\prime})$ are neighbors, $\boldsymbol{N}\sim\boldsymbol{N^{\prime}}$, if $\sum_{i=1}^{4}|N_{i}-N^{\prime}_{i}|=2$. This means that $\boldsymbol{N}^{\prime}$ can be obtained from $\boldsymbol{N}$ by moving one element to a different class. We construct the set $I_{\varepsilon}$ as follows.

At the onset, $I_{\varepsilon}$ contains only the mode of multinomial distribution. We find the mode by a simple hill climbing algorithm, which starts with a state close to the mean of the multinomial distribution and follows the direction of increasing probability until the maximum is reached. Because of unimodality, it is indeed a global maximum. In the next steps, we add the neighbors of states which were previously visited. The procedure is repeated until the total probability of set $I_{\varepsilon}$ reaches the desired value $1-\varepsilon$. The details of the above method can be found in [28]. We construct the set $I_{\varepsilon}$ and we estimate the p-value by

[TABLE]

Bootstrap Procedure

The bootstrap procedure has gained mainstream popularity for its wide applicability and statistical treatments [29]. Creating an empirical distribution of null statistics allows for a flexible and robust estimation of p-values and related statistics. We show how to use the resampling with replacement to obtain statistical significance of $T^{c}(\mbox{\bf{y}}_{i},\mbox{\bf{y}}_{j})$. Particularly, resampling with replacement $\mbox{\bf{y}}_{i}$ and $\mbox{\bf{y}}_{j}$, separately, breaks any potential dependency. This allows us to calculate an empirical distribution of Jaccard/Tanimoto coefficients under the null hypothesis:

The expectation of Jaccard/Tanimoto coefficients is estimated directly from resampled vectors $\mbox{\bf{y}}^{*}_{i}$ and $\mbox{\bf{y}}^{*}_{j}$, that are effectively independent. Therefore, each iteration provides randomness, which helps avoid a bias related to using an estimated expectation based only on observation. Previously, there are early works in Monte Carlo procedures [30, 14] and published statistical tables for assessing randomness in species co-occurances [22, 23]. However, earlier works have assumed that a probability of occurrences is 0.5 regardless of species or biogeographic units. Permutation methods based on conventional uncentered coefficients are available in R packages, whose operating characteristics are not described in details [31, 32].

The resolution of the empirical null distribution depends on $B$, where the larger $B$ will result in more precise estimation of p-values. Although the choice of $B$ would likely be dictated by $n$ and $m$, as well as available computational power, we recommend setting $B$ to at least 5-10 times of $m$. In our simulation studies, the total bootstrap iterations is set to $B=5\times m$, which are shown to be both accurate and fast. When comparing a very large set of species or OTUs, it may be helpful to pool null statistics to increase the p-value resolution and speed up the computation.

Results and discussion

Simulation Studies

We have developed statistical methods and algorithms to obtain statistical significance of Jaccard/Tanimoto similarity coefficients for biological presence-absence data. Beyond deriving the exact solution, we introduce the measurement concentration algorithm (MCA) and bootstrap method. We characterize their operating characteristics by comprehensive simulation studies where a wide range of parameters for presence-absence datasets are considered. Our goal is to maintain theoretically correct behaviors of p-values. Null p-values corresponding to $H_{0}$ are evaluated against a Uniform(0,1) distribution. False discovery rates (FDRs) are directly estimated from p-values produced by our methods to demonstrate an overall error control.

First, we conducted 5 simulation scenarios using different underlying occurrence probabilities $p=0.1,0.3,0.5,0.7,0.9$ to generate independent presence-absence datasets. In essence, they are two species of length $m=100$ that exhibit unrelated co-occurrence patterns, where a proportion of presence (1’s) ranges from $10\%$ to $90\%$. For each of simulation scenarios, a total of $2000$ comparisons were made using a length $m=100$. Without any information about simulation parameters, our proposed methods are applied on an identically simulated dataset (Figure 2). Theoretically correct p-values under the null hypothesis (null p-values) should form a Uniform distribution between 0 and 1, which are denoted by dashed diagonal lines in QQ plots. An upward deviation from diagonals shows an anti-conservative bias, as shown among some asymptotic p-values. In all scenarios, p-values from the exact solution, bootstrap ($B=500$), and measure concentration (accuracy $=1\times 10^{-5}$) algorithms follow a theoretically correct Uniform(0,1) distribution (Figure 2). Asymptotic approximation is inconsistent; its behavior is anti-conservative with $p=0.3,0.5$ and slightly conservative with $p=0.7,0.9$. Asymptotic approximation should only be used when computational time is a critical bottleneck.

Second, we generated a mixture of independent and dependent datasets out of $n=2000$ presence-absence vectors (of $n=2000$ species observed in $m=200$ biogeographic units) to evaluate false discovery rates. In three separate scenarios, we simulated $25\%$, $50\%$, and $75\%$ of $n=2000$ species to be independent, resulting in null proportions of $\pi_{0}=.25,.50,.75$ respectively. For example, a scenario with $\pi_{0}=.75$ produces 500 out of $n=2000$ presence-absence variables that are truly associated with the query variable. Then, our proposed asymptotic approximation, bootstrap method, and measure concentration algorithm (MCA) are used to automatically compute p-values. To account for variation in simulation, we repeated each scenario 20 times. False discovery rates (FDRs) and $\pi_{0}$ are estimated by the q-value methodology [33]. Q-values are evaluated against FDR thresholds, so that we can evaluate accuracy of observed FDRs (Figure 3). Twenty simulation replications are shown in semi-transparent shades, whereas their group averages for 3 methods are shown as solid lines. An upward deviation as shown by asymptotic approximation indicates an overall anti-conservative behavior, likely due to $m\not\to\infty$. The bootstrap and MCA maintain the overall error rates, where the bootstrap exhibits slightly conservative characteristics (Figure 3).

Third, we compared the computational efficiency of our proposed methods using our jaccard package on RStudio Cloud (Intel Xeon 2.90GHz and 1GB RAM), with R 3.5.0. We measured the runtime for a range of lengths $m=50,\ldots,500$. For each $m$, we applied the proposed methods 10 times, with the bootstrap iteration $B=5\times m$ and MCA accuracy of $1\times 10^{-5}$. The average runtimes are shown Figure 4. Our proposed computational methods show drastic improvement over the exact solution as $m$ increases. The asymptotic approximation is mostly instantaneous. When the similarity between two presence-absence vectors of length $m=500$ were tested using the jaccard package, the exact solution was prohibitively slow, taking 41.5s on average. The bootstrap method was 449.8 times (0.09s) faster, whereas MCA was 92.5 times (0.45s) faster than the exact solution. Furthermore, we compared the runtimes of estimation methods for $m=1000,\ldots,10000$ (Figure S1). The gain in computational efficiency is more pronounced as the dimension (i.e., a length of presence-absence vectors) grows in size.

Last, a simulation study with $p=0.5$ and $m=200$ was used to evaluate two recent methods of species co-occurrences analysis. We generated independent presence-absence data where two species are truly unrelated. Then, methods of combinatorics [24] and hypergeometric distributions [25] are applied to obtain p-values. We followed the recommendations given in each paper, displaying four possible p-values from [24] (Figure S2) and two one-sided p-values from [25] (Figure S3). We observe these p-values under the null hypothesis to substantially deviate from theoretically correct Uniform(0,1) distributions.

Applications in Species Co-occurrences

To show applications in statistically testing biological presence-absence data, the proposed methods are applied to species co-occurrence data. We investigated bird species on 28 islands in the Republic of Vanuatu, that are available in [6] and analyzed in several pioneering studies in non-random co-occurrences of species [12, 14, 15, 16]. The data is consisted of presence and absence of bird species in 28 islands of Vanuatu, which used to be known as the New Hebrides. Three generalist species that existed in all 28 islands were removed from our analysis. We are interested in identifying what pairs of species exhibit statistically significantly co-occurrences.

For $n=53$ bird species in $m=28$ islands, we obtained 1378 pair-wise Jaccard/Tanimoto similarity coefficients. The conventional Jaccard/Tanimoto coefficients depends strongly on their expected values under independence (Figure 5). Similarly, the conventional Jaccard/Tanimoto coefficients are substantially correlated with the proportion of occurrences, with a Pearson correlation of $0.43$ (p-value $<2.2\times 10^{-16}$). Relying only on similarity coefficients would miss non-random co-occurrences among bird species that live in a few islands (Figure S4). Our proposed methods account for co-occurrences that would be expected under independence. Histograms of the uncentered and centered Jaccard/Tanimoto coefficients are compared in Figure S5.

We computed statistical significance by applying the bootstrap method with $B=5000$ and MCA with accuracy of $1\times 10^{-5}$. Our two computational approaches estimated p-values that are almost identical with their mean squared deviation of $1.15\times 10^{-4}$ (Figure S6). Significant results that are substantially deviating from random samples indicate non-random co-occurrences of species (Figure 6). Out of 1378 pairs of species that were tested, the proportion of independent specie pairs was estimated to be $24\%$ using q-value methodology [33]. Then, we calculated false discovery rates (FDRs) from 1378 pair-wise p-values. We discovered that 374 ($27\%$) pairs are deemed significant at a q-value threshold of 0.10.

Additionally, we applied the Jaccard/Tanimoto similarity tests among fish species in French freshwater streams, surveyed over a long period of time [34]. Briefly, the presence and absence data of the $n=32$ most common fish species in $m=3347$ sites across French rivers are obtained during 1980 - 1991 [34]. Our analysis estimates that about $84.3\%$ of 496 pairs are estimated to be non-randomly co-occurring. As surveyed for over a decade across Fresh rivers and surrounding habitats, it is reasonable that many fish species are interacting or influenced by related climate conditions. There are 21 pairs of species with q-values $>0.1$ (corresponding p-values ranging from 0.637 and 0.969). For example, the centered statistics between Pungitius pungitius and Cyprinus carpio is $3.31\times 10^{-4}$, whereas that between Pungitius pungitius and Lota lota is $-4.40\times 10^{-4}$. P. pungitius is a small fish species typically riding in thick submerged vegetation with the breeding season falling in April - July. C. carpio and L. lota are much bigger species and generally prefers a large body of water.

Conclusion

From biogeography to microbiology, evaluating similarity among species and biogeographic units is fundamental to assessing co-existence and biodiversity. Having observed occurrences of species in multiple biogeographic units, one of the primary goals in analyzing presence-absence data is to identify non-random co-occurrences. Even if two species would be present independently of each other, they may occur together by chance. For the last 30 years, the Jaccard/Tanimoto coefficient has been shown to be highly useful for quantitative analysis of co-occurrences that help inform systematic relationship among species [3, 4, 5]. We have developed a rigorous statistical framework and methods to efficiently calculate statistical significance of such similarity and to identify non-random co-occurrences.

For testing co-occurrences using the Jaccard/Tanimoto coefficient, we introduce exact and asymptotic solutions, as well as bootstrap and measure concentration algorithm. The proposed suite of statistical methods can provide a rigorous guideline to identify related species. Through comprehensive simulation studies, we characterized their operating characteristics using p-values and FDRs. The proposed bootstrap and measure concentration algorithms are highly accurate and efficient, providing orders of magnitude improvement in a computational speed. We have implemented the proposed methods in an open source R package and a Shiny web app. A user can upload a dataset to be analyzed, and create histograms and heat maps automatically. This will facilitate adaptation of p-values, FDRs, and related quantities in analyzing species co-occurrences.

Beyond species co-occurrences, the Jaccard/Tanimoto coefficient is used in diverse areas of biological science where binary data are observed and compared. When molecules and reactions are represented as hashed fingerprints, it is used for quantitative comparisons and classifications [35, 36, 37]. Similarity between biochemical reactions can be tested by applying our methods on their corresponding fingerprints. In genomics, the standard tools such as BEDTools [38] evaluate genomic intervals using the Jaccard/Tanimoto coefficients. Given genomic intervals from two samples or groups, one could test whether their overlap is statistically significant, providing evidences for shared genomic variations. Due to the popularity of Jaccard/Tanimoto coefficients, the proposed suite of methods would be useful in a broad range of scientific applications.

Figures

Supplementary Figures