Tree-based Inference of Species Interaction Network from Abundance Data

TL;DR

This paper introduces a tree-based statistical model for inferring ecological species interaction networks from abundance data, accounting for environmental effects and distinguishing direct from indirect interactions efficiently.

Contribution

It presents a novel, computationally efficient method that incorporates environmental covariates and models species interactions through averaging over tree-shaped networks.

Findings

Method performs well compared to state-of-the-art approaches.

Accounting for covariates reduces spurious edges.

Application to datasets reveals environmental influences on networks.

Abstract

The behavior of ecological systems mainly relies on the interactions between the species it involves. We consider the problem of inferring the species interaction network from abundance data. To be relevant, any network inference methodology needs to handle count data and to account for possible environmental effects. It also needs to distinguish between direct interactions and indirect associations and graphical models provide a convenient framework for this purpose. We introduce a generic statistical model for network inference based on abundance data. The model includes fixed effects to account for environmental covariates and sampling efforts, and correlated random effects to encode species interactions. The inferred network is obtained by averaging over all possible tree-shaped (and therefore sparse) networks, in a computationally efficient manner. An output of the procedure is the probability for each edge to be part of the underlying network. A simulation study shows that the proposed methodology compares well with state-of-the-art approaches, even when the underlying graph strongly differs from a tree. The analysis of two datasets highlights the influence of covariates on the inferred network. Accounting for covariates is critical to avoid spurious edges. The proposed approach could be extended to perform network comparison or to look for missing species.

Tables

Table 1. Table 1: Median and standard-deviation running-time values (in seconds) for Cluster and Erdös structures, including resampling with S = 20 𝑆 20 S=20 for SpiecEasi and EMtree.

	SpiecEasi	gCoda	ecoCopula	MRFcov	MInt	EMtree
Easy	25.45 (1.87)	0.11 (0.06)	5.55 (0.64)	34.51 (3.68)	43.04 (19.76)	11.72 (1.89)
Hard	28.43 (1.30)	0.53 (0.25)	9.6 (0.65)	8.29 (0.36)	33.77 (18.20)	8.17 (0.50)

Table 2. Table 2: Median and standard-deviation of running times for each method in seconds, for n 𝑛 n and p 𝑝 p parameters.

		$n\ge 50$		$p\ge 20$
EMtree	0.44 (0.14)	0.60 (0.17)	0.41 (0.13)	0.76 (0.21)
gCoda	0.11 (26.8)	0.05 (0.05)	0.05 (0.04)	0.09 (0.54)
SpiecEasi	2.09 (0.26)	2.37 (0.28)	2.42 (0.27)	2.42 (0.26)

Table 3. Table 3: Medians and standard-deviation of FDR computed on 100 graphs of each type ( easy : n = 100 , p = 20 formulae-sequence 𝑛 100 𝑝 20 n=100,p=20 , hard : n = 50 , p = 30 formulae-sequence 𝑛 50 𝑝 30 n=50,p=30 )

		SpiecEasi	gCoda	ecoCopula	MRFcov	MInt	EMtree
Easy	Cluster	0.86 (0.20)	0 (0.08)	0.33 (0.14)	0.74 (0.06)	0.38 (0.17)	0.12 (0.09)
	Erdös	0.86 (0.21)	0 (0.15)	0.29 (0.15)	0.73 (0.05)	0.38 (0.15)	0.12 (0.08)
	Scale-free	0.92 (0.04)	0 (0.04)	0.33 (0.11)	0.88 (0.02)	0.73 (0.09)	0.34 (0.08)
Hard	Cluster	0.87 (0.12)	0 (0.20)	0.15 (0.18)	0.78 (0.05)	0.77 (0.09)	0.36 (0.09)
	Erdös.	0.88 (0.11)	0 (0.24)	0 (0.15)	0.78 (0.05)	0.77 (0.10)	0.35 (0.09)
	Scale-free	0.94 (0.05)	0 (0.13)	0 (0.16)	0.89 (0.03)	0.94 (0.03)	0.56 (0.07)

Table 4. Table 4: Medians and standard-deviation of density ratio computed on 100 graphs of each type ( easy : n = 100 , p = 20 formulae-sequence 𝑛 100 𝑝 20 n=100,p=20 , hard : n = 50 , p = 30 formulae-sequence 𝑛 50 𝑝 30 n=50,p=30 )

		SpiecEasi	gCoda	ecoCopula	MRFcov	MInt	EMtree
Easy	Cluster	0.16 (0.11)	0.05 (0.07)	1.04 (0.48)	2.26 (0.58)	0.30 (0.13)	0.62 (0.14)
	Erdös	0.15 (0.09)	0.06 (0.08)	0.95 (0.50)	2.23 (0.42)	0.30 (0.14)	0.58 (0.12)
	Scale-free	0.63 (0.13)	0.08 (0.07)	0.92 (0.30)	4.86 (0.44)	0.81 (0.25)	0.96 (0.08)
Hard	Cluster	0.21 (0.08)	0.02 (0.03)	0.02 (0.17)	1.65 (0.33)	0.68 (0.30)	0.43 (0.10)
	Erdös	0.21 (0.08)	0.02 (0.02)	0.00 (0.18)	1.56 (0.32)	0.66 (0.25)	0.42 (0.10)
	Scale-free	0.61 (0.12)	0.04 (0.03)	0.08 (0.24)	3.29 (0.40)	3.63 (1.08)	0.94 (0.09)

Table 5. Table 5: Percentage of empty networks computed on 100 graphs of each type ( easy : n = 100 , p = 20 formulae-sequence 𝑛 100 𝑝 20 n=100,p=20 , hard : n = 50 , p = 30 formulae-sequence 𝑛 50 𝑝 30 n=50,p=30 )

		SpiecEasi	gCoda	ecoCopula	MRFcov	MInt	EMtree
Easy	Cluster	1.77	13.89	1.74	0.00	0.00	0.00
	Erdös	0.68	11.95	0.99	0.00	0.83	0.00
	Scale-free	0.00	1.88	0.00	0.00	0.00	0.00
Hard	Cluster	0.00	14.05	23.40	0.00	0.00	0.00
	Erdös	0.00	20.85	27.28	0.00	0.00	0.00
	Scale-free	0.00	5.97	15.46	0.00	0.00	0.00

Table 6. Table 6: Median and standard-deviation of running times in seconds for scale-free structures, for two dataset dimensions ( easy : n = 100 𝑛 100 n=100 , p = 20 𝑝 20 p=20 , hard : n = 50 𝑛 50 n=50 , p = 30 𝑝 30 p=30 ).

	SpiecEasi	gCoda	ecoCopula	MRFcov	MInt	EMtree
Easy	29.73 (2.00)	1.29 ( 0.30)	28.14 (1.46)	48.7 (2.32)	138.13 (39.60)	50.19 (7.81)
Hard	29.38 (1.31)	40.73 (20.94)	27.16 (1.30)	14.1 (0.36)	95.27 (46.34)	23.59 (2.09)

Table 7. Table 7: Median and standard-deviation running-time values in seconds for inference of Erdös structure with EMtree and different values of the number of sub-samples S 𝑆 S .

$S$	1	2	10	20	50	150
Easy	0.66 (0.15)	1.86 (0.23)	7.00 (0.81)	12.29 (1.27)	29.50 (3.39)	87.30 (10.36)
Hard	0.45 (0.12)	1.44 (0.14)	5.06 (0.78)	8.97 (0.87)	23.35 (2.40)	69.29 (10.83)

Table 8. Table 8: Median and standard-deviation of running times for each method in seconds, for n 𝑛 n and p 𝑝 p parameters. corresponding to Erdös and cluster structures with 5 / p 5 𝑝 5/p densities.

		$n\ge 50$		$p\ge 20$
EMtree	0.41 (0.11)	0.60 (0.15)	0.38 (0.12)	0.71 (0.21)
gCoda	0.12 (0.47)	0.07 (0.03)	0.05 (0.03)	0.09 (0.06)
SpiecEasi	2.41 (0.25)	2.41 (0.25)	2.39 (0.25)	2.42 (0.25)

Equations

$$Y_{ij}\sim \mathcal{P}\left(\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\theta }}_j+o_{ij}+Z_{ij})\right).$$

$$p(T)=\prod _{(j,k)\in T}{\beta }_{jk}/B,$$

$$p({\mathbf{Z}}_i\mid T)=\prod _{j=1}^{p}p(Z_{ij}|T)\prod _{(j,k)\in T}{\psi }_{jk}({\mathbf{Z}}_i),$$

$$p_{\mathbf{\beta },\mathbf{\Sigma },\mathbf{\theta }}(T,\mathbf{Z},\mathbf{Y})=p_{\mathbf{\beta }}(T)\times p_{\mathbf{\Sigma }}(\mathbf{Z}|T)\times p_{\mathbf{\theta }}(\mathbf{Y}|\mathbf{Z}),$$

$$E\left(\log p_{\mathbf{\beta },\mathbf{\Sigma },\mathbf{\theta }}(T,\mathbf{Z},\mathbf{Y})|\mathbf{Y}\right)\simeq \sum _{1\le j<k\le p}{P}_{jk}\log \left({\beta }_{jk}{\hat{\psi }}_{jk}\right)-\log B+\text{cst}$$

$$\mathcal{J}(\mathbf{Y};\mathbf{\theta },\mathbf{\Sigma },\tilde{p}):=\log p_{\mathbf{\theta },\mathbf{\Sigma }}(\mathbf{Y})-KL\left(\tilde{p}(\mathbf{Z})||p_{\mathbf{\theta },\mathbf{\Sigma }}(\mathbf{Z}\mid \mathbf{Y})\right),$$

$$\log p_{\mathbf{\Sigma }}(\mathbf{Z}\mid T)$$

$$=-\frac{n}{2}\log {\sigma }_j^2-\frac{1}{2}\sum _{j=1}^p\sum _{i=1}^n\frac{Z_{ij}^2}{{\sigma }_j^2}-\frac{n}{2}\sum _{(j,k)\in T}\log (1-{\rho }_{jk}^2)$$

$$-\frac{1}{2}\sum _{(j,k)\in T}\frac{1}{1-{\rho }_{jk}^2}\sum _{i=1}^n\left({\rho }_{jk}^2\frac{Z_{ij}^2}{{\sigma }_j^2}+{\rho }_{jk}^2\frac{Z_{ik}^2}{{\sigma }_k^2}-2{\rho }_{jk}\frac{Z_{ij}{Z}_{ik}}{{\sigma }_j{\sigma }_k}\right)+\text{cst}$$

$${\hat{\sigma }}_j^2=\frac{1}{n}\sum _{i}E(Z_{ij}^2\mid \mathbf{Y}),{\hat{\rho }}_{jk}=\frac{1}{n}\sum _{i}E(Z_{ij}{Z}_{ik}\mid \mathbf{Y})/({\hat{\sigma }}_j{\hat{\sigma }}_k).$$

$$E\left(\log p_{\hat{\mathbf{\Sigma }}}(\mathbf{Z}\mid T)\mid \mathbf{Y}\right)=-\frac{n}{2}\log {\hat{\sigma }}_j^2-\frac{n}{2}\sum _{(j,k)\in T}\log (1-{\hat{\rho }}_{jk}^2)+\text{cst}.$$

$$E\left(\log p_{\mathbf{\beta }}(T)+\log p_{\hat{\mathbf{\Sigma }}}(\mathbf{Z}\mid T)\mid \mathbf{Y}\right)$$

$$=-\log B+\sum _{T\in \mathcal{T}}p(T\mid \mathbf{Y})\sum _{(j,k)\in T}\left(\log {\beta }_{jk}+\log {\hat{\psi }}_{jk}\right)+\text{cst}$$

$$=-\log B+\sum _{(j,k)}\text{mathds}\mathrm{P}\{(j,k)\in T\mid \mathbf{Y}\}\left(\log {\beta }_{jk}+\log {\hat{\psi }}_{jk}\right)+\text{cst},$$

$$[\mathbf{Q}{]}_{jk}=\{\begin{array}{ll}-w_{jk}& 1\le j<k\le p\\ {\sum }_{u=1}^p{w}_{ju}& 1\le j=k\le p.\end{array}$$

$$W:=\sum _{T\in \mathcal{T}}\prod _{(j,k)\in T}{w}_{jk}=|{\mathbf{Q}}^{uv}|.$$

$$[\mathbf{M}{]}_{jk}=\{\begin{array}{ll}{\left[({\mathbf{Q}}^{pp}{)}^{-1}\right]}_{jj}+{\left[({\mathbf{Q}}^{pp}{)}^{-1}\right]}_{kk}-2{\left[({\mathbf{Q}}^{pp}{)}^{-1}\right]}_{jk}& 1\le j<k<p\\ {\left[({\mathbf{Q}}^{pp}{)}^{-1}\right]}_{jj}& k=p,1\le j\le p\\ 0& 1\le j=k\le p.\end{array}$$

Full text

1. The behavior of ecological systems mainly relies on the interactions between the species it involves. We consider the problem of inferring the species interaction network from abundance data. To be relevant, any network inference methodology needs to handle count data and to account for possible environmental effects. It also needs to distinguish between direct interactions and indirect associations and graphical models provide a convenient framework for this purpose.

2. We introduce a generic statistical model for network inference based on abundance data. The model includes fixed effects to account for environmental covariates and sampling efforts, and correlated random effects to encode species interactions. The inferred network is obtained by averaging over all possible tree-shaped (and therefore sparse) networks, in a computationally efficient manner. An output of the procedure is the probability for each edge to be part of the underlying network.

3. A simulation study shows that the proposed methodology compares well with state-of-the-art approaches, even when the underlying graph strongly differs from a tree. The analysis of two datasets highlights the influence of covariates on the inferred network.

4. Accounting for covariates is critical to avoid spurious edges. The proposed approach could be extended to perform network comparison or to look for missing species.

Key-words:

abundance data, covariates adjustment, EM algorithm, graphical models, matrix tree theorem, Poisson log-Normal model, species interaction network

1 Introduction

There is a growing awareness of biotic interactions being crucial components of biodiversity and relevant descriptors of ecosystems (Valiente-Banuet et al., 2015; Jordano, 2016). Such interactions can be conveniently represented by networks, which have been increasingly studied and used in recent years for describing and understanding living systems in ecology (Poisot et al., 2016), microbiology (Faust and Raes, 2012) or genomics (Evans et al., 2016). Observing species interactions is a laborious task which restricts them to certain categories (e.g. trophic, pollination), while many other mutualistic and/or antagonistic interactions may be hard to observe and key in the system organization (e.g. communication, shelter sharing). Many efforts have been devoted in the last decade to get a more complete picture of the biotic interactions existing between species living in the same niche.

Network reconstruction.

A first attempt consists in using observed interactions to predict other possible links based on species traits matching (see e.g. Olito and Fox, 2015; Bartomeus et al., 2016; Weinstein and Graham, 2017; Graham and Weinstein, 2018). The interaction strength can also be predicted (Wells and O’Hara, 2013). This can be viewed as a prediction task, and modern approaches arising from signal processing and machine learning have been also proposed (Desjardins-Proulx et al., 2017; Stock et al., 2017; Dallas et al., 2017). We name these approaches network reconstruction to distinguish them from network inference, which is the problem we consider in this article.

Network inference.

Network inference approaches also aim at retrieving the interactions among species, but do not rely on observed interactions and therefore, remain agnostic as for their type. Such approaches have been developed in many domains ranging from cell biology (Friedman, 2004, to infer gene regulatory networks) to neurosciences (Zhu and Cribben, 2018, to deciphere brain connectivity structures). In ecology, it will typically aim at inferring the set of biotic interactions linking species from the same guild. As summarized in Fig. 1, network inference takes as input measures on species at similar sites, and returns a network of species direct interactions. The importance of distinguishing between direct interaction and indirect association between species is explained in Popovic et al. (2019). To be accurate, network inference must account for environmental covariates to prevent the inference of spurious interactions resulting from abiotic effects. Fig. 1 illustrates this phenomenon: ($c$) corresponds to the case where two species (1 and 4) are not in direct interaction, but are affected by the variations of the same environmental covariate $x$. ($d$) displays the network when $x$ is not accounted for: a spurious edge appears between these two species.

Joint species distribution models.

The rational between network inference is that interactions between species must affect their joint distribution in a series of similar sites. Such approaches necessarily rely on a joint species distribution model (JSDM), as opposed to species distribution models (Elith and Leathwick, 2009) where species are traditionally considered as disconnected entities. A JSDM is a probabilistic model describing the species simultaneous presence/absence (Harris, 2015; Ovaskainen et al., 2017) or joint abundances (Popovic et al., 2018, 2019). An important feature of JSDMs is to include environmental covariates to account for abiotic interactions.

Recently, latent variable models have received attention in community ecology as they provide a convenient way to model the dependence structure between species (Warton et al., 2015). The JSDM proposed by Popovic et al. (2018, 2019) involves a latent layer. So does the Poisson log-Normal model (PLN, Aitchison and Ho, 1989), which combines generalized linear models to account for covariates and offsets, and a Gaussian latent structure to describe the species interactions. It can be seen as a multivariate mixed model, in which correlated random effects encode the dependency between the species abundances.