Network inference with LIMITS

We start by loading the seqtime library.

library(seqtime)

Next, we generate a “known” interaction matrix randomly.

N=20
A=generateA(N, c=0.1)

## [1] "Adjusting connectance to 0.1"
## [1] "Initial edge number 400"
## [1] "Initial connectance 1"
## [1] "Number of edges removed 342"
## [1] "Final connectance 0.1"
## [1] "Final connectance: 0.1"

rownames(A)=c(1:N)
colnames(A)=rownames(A)

Below, we plot our known interaction matrix. Green entries stand for a positive interaction of species i with species j, whereas red entries stand for a negative interaction of species i with species j (where i and j are row and column indices, respectively). For instance, a competition consists of two red entries for the negative impact of species i on j and species j on i. The diagonal is entirely red, since the interaction matrix is generated by default with negative values on the diagonal, which represent intra-species competition.

plotA(A, header="Known interaction matrix")

## [1] "Largest value: 0.499967047944665"
## [1] "Smallest value: -0.5"

We can also visualize the interaction matrix as a network. In contrast to correlation networks, the network is directed. Thus, it can represent asymmetric ecological relationships. For example, consider a parasitic relationship, where species i benefits at the cost of species j. It would be represented in the matrix by a positive entry in row i and column j and a negative entry in row j and column i. In the directed network, this parasitic relationship would correspond to two arcs, one from species i to species j and a reverse one from species j to species i. The negative diagonal entries are represented as self-loops.

network=plotA(A,method="network")

## [1] "Largest value: 0.499967047944665"
## [1] "Smallest value: -0.5"
## [1] "Initial edge number 58"
## [1] "Initial connectance 0.1"
## [1] "Final connectance: 0.1"

# use igraph's function tkplot for manual layout of the network

Given our interaction matrix, we can now simulate a test time series with the Ricker community model.

out.ricker=ricker(N,A=A)
tsplot(out.ricker,main="Ricker")

Now, we run LIMITS to test how well it can infer the known interaction matrix from the time series. LIMITS was developed by Fisher and Mehta PLoS ONE 2014. LIMITS needs a minute to run and then returns the estimated interaction matrix.

Aest=limits(out.ricker)$Aest

## [1] "Time series has 20 taxa"

To compare the known and inferred interaction matrix, we can plot them next to each other:

par(mfrow=c(1,2))
plotA(A,header="known")

## [1] "Largest value: 0.499967047944665"
## [1] "Smallest value: -0.5"

plotA(Aest,header="inferred")

## [1] "Largest value: 0.763672072136036"
## [1] "Smallest value: -1.77829739290894"

par(mfrow=c(1,1))

We can also compute the mean cross-correlation of the known and inferred interaction matrix:

crossCor=cor(A,Aest)
mean(diag(crossCor), na.rm=TRUE)

## [1] 0.7347679

Finally, we can plot a few quality estimators for the interaction matrix inference. The quality plot displays the correlation between current and future time points one to five steps ahead (autocor) and the correlation between the original time series and a time series generated from the inferred interaction matrix step by step (cor). We can see that the predicted time series barely outperforms lag-one auto-correlation.

limitsqual=limitsQuality(out.ricker,A=Aest,plot=TRUE)

## [1] "Applying Schur decomposition"
## [1] "Initial edge number 54"
## [1] "Initial connectance 0.0894736842105263"
## [1] "Final connectance: 0.421052631578947"

For comparison, we now generate a neutral time series with the Hubbell model. This model does not take species interactions into account explicitly (though implicitly there is a competition between all species) and hence does not take the interaction matrix A as an input.

out.hubbell=simHubbell(N=N, M=N,I=1500,d=N, m=0.1, tskip=500, tend=1000)
tsplot(out.hubbell,main="Hubbell")

We infer an interaction matrix from the neutral time series:

Aesth=limits(out.hubbell)$Aest

## [1] "Time series has 20 taxa"

The quality plot shows that the inferred interaction matrix leads to a high cross-correlation between predicted and original time series, even though the dynamics in the neutral time series is not determined by species-specific interactions (all interactions are negative and of equal strength, since no species has an advantage over another in the competition for a space in the local community). Thus, a high correlation between predicted and original time series can be misleading as an indicator for the quality of the inferred interaction matrix.

limitsqualh=limitsQuality(out.hubbell,A=Aesth, plot=TRUE)

## [1] "Applying Schur decomposition"
## [1] "Initial edge number 48"
## [1] "Initial connectance 0.0736842105263158"
## [1] "Final connectance: 0.944736842105263"

Karoline Faust

2019-12-15