Food web prediction is trivial

if you do it wrong

I currently spend a lot of time figuring out ways to predict ecological networks, and food webs and other antagonistic interactions in particular. But as it turns out, it is very easy to game this problem and obtain fantastic performance For your model. All you have to do is be spectacularly wrong. Don’t try this at home, and please don’t try this at work.

We will only need to use two packages, which if anything, goes to show that even with the simplest of setups, we can do a lot of damage. We will also grab Statistics from the standard library, and this should be enough.

using EcologicalNetworks
using DecisionTree
using Statistics

I apologize in advance for what follows.

Getting the data we have

We will attempt to forecast the host-parasite metaweb from Eurasia, which is my usual benchmark dataset because it is spatially replicated, has many overlapping species, and still displays some interesting variation. It comes with EcologicalNetworks through the web_of_life function, so we can get it, and transform it to a bipartite metaweb.

good_dataset(ds) = contains("Hadfield")(ds.Reference)
N = map(i -> web_of_life(i.ID), filter(good_dataset, web_of_life()))
B = convert.(BipartiteNetwork, N)
M = reduce(union, B)
206×121 bipartite  ecological network (Bool, String) (L: 2131)

The M network is a metaweb with all species and interactions across the 51 locations, and this is what we will use for the prediction. Because the metaweb stores boolean interactions (yes/no), this is a classification problem. And for this, we need features.

Faking the data we don’t have

What I would really like to do is predict interactions from species traits. I cannot do that because I do not have species traits. This looks like an unsolvable situation.

Wrong! What I do have is a series of ways to generate pseudo-random numbers. Let’s assume that both the hosts and the parasites have three traits each, and let’s also assume that globally, these traits follow a uniform distribution between 0 and 1, because this has, as you will see, literally no effect on the performance of the prediction.

host_traits = rand(Float64, (3, richness(M; dims=2)))
parasite_traits = rand(Float64, (3, richness(M; dims=1)))

Well, that was easy. Now, we want to make predictions, so we will unfold the interaction matrix (this is our target for prediction), and generate a feature vector for every entry in the matrix, which will have 6 entries (the traits of both species). If the interactions are determined by traits, whatever predictive tool we use should figure out associations between both.

Generating features and labels

The next step is to build a matrix of features:

matrix_entries = prod(size(M))
labels = zeros(Bool, matrix_entries)
features = zeros(Float64, (6, matrix_entries))

cursor = 0
for i in 1:richness(M; dims=1)
  for j in 1:richness(M; dims=2)
    global cursor += 1
    labels[cursor] = M[i,j]
    features[1:3, cursor] = host_traits[:,j]
    features[4:6, cursor] = parasite_traits[:,i]
  end
end

And we’re good to go. Notice how we haven’t made any effort to make the code efficient or nice? It’s because the problem is, as you will see, so trivial, that we don’t need to do this. Anything goes, my dudes.

Nailing the prediction

Let’s just drop all of this into a random forest. We will need to transpose the features vectors, because DecisionTree and Flux have different opinions about how these data should be laid out, and I like Flux’s convention more.

Anyways, we can get a sense of the performance of the model with n-fold cross-validation, using 20 folds, and 3 sub-features every time:

n_folds=20; n_subfeatures=3;
accuracy = nfoldCV_forest(labels, features', n_folds, n_subfeatures);

With this done, we can look at the accuracy:

mean(accuracy)
0.9132022471910112

Wow! This is very good! Mission accomplished, we have estimated about 90% of the interactions correctly.

Except we didn’t.

Wait, what?

Here’s a fun fact about networks - they are extremely sparse. They tend not to have a lot of interactions at all. We can actually measure the proportion of the matrix which is interactions, and in the case of our metaweb, this is

connectance(M)
0.08549305945598973

This means that a model that would predict all zeros would get an accuracy of about

1.0 - connectance(M)
0.9145069405440103

You might notice that this is just about the same as the accuracy of our “model” above. I like this example, because it shows how easy it is to get fooled by the numbers that come out of the packages. An accuracy of 90% seems good (it’s really not that fantastic but ecological data are small and complicated, so it’s nothing to be ashamed of). But in this case, all that it reflects is that most of the dataset is zeroes. A direct consequence of this is: whatever happens to the ones doesn’t matter. By chance, we should get a lot of predicted zeroes, and we will also have a lot of true zeroes, and so the True Negative compartment is going to be so overwhelmingly large that accuracy is bound to look good.

For an in-progress paper, we are working on a worked version of this example where we show that it is possible to produce a model that reaches 95% accuracy on this dataset and that actually means something (I will post this here when the paper is more advanced).

Now is also a good time to remember that we are using random traits. This does not mean that interactions are random - it only means that models are really good at learning something from the data, even if this something is a bunch of nonsense. In this case, the nonsense that the trees learned is “pretty much any combination of values will give you a zero”. I like this example as a cautionary tale - defining a “good” fit for the model demands careful consideration.

There is a very easy way to see that the model is not doing really well - Cohen’s κ. To do it a little more properly, we can get 15000 samples in our training set, and the rest in the testing set.

using StatsBase
train = sample(1:length(labels), 15000, replace=false)
test = filter(i -> !(i in train), 1:length(labels))

model = build_forest(labels[train], features[:,train]')
preds = apply_forest(model, features[:,test]')

Let’s look at the confusion matrix (which DecisionTree kindly gives us with the accurracy and κ scores):

confusion_matrix(labels[test], preds)

The value of κ is around 0.05, which usually indicates “very low accord” between the prediction and reality. Notice that the accuracy on this model is still very high, but the confusion matrix makes it very clear that only 31 interactions were correctly predicted. The model is bad, and we should not use it.

In the next post (or series of posts?), I will show a few ways to get some actual predictions on this dataset, using a mix of good practices and insights about ecology. Stay tuned.