Subtle and obvious differences between SDM and trophic SDM

Compare trophic and classical SDMs: difference and equalities

Here, we present the differences between a classic SDM and our trophic SDM, starting from the simplest case when these two approaches coincide to more realistic cases where they do not.

The simplest case

Let’s consider a very simple case with one environmental covariate x, and two species y_A and y_B. We assume species y_A to be basal and species y_B to feed on y_A. This leads to the following directed acyclic graph that resumes the dependencies across variables, assuming bottom-up control.

library(igraph)
G=graph_from_edgelist(as.matrix(data.frame(from = c("x","x","A"),
                                           to =  c("A","B","B"))))
plot(G)

plot of chunk igraph

We now assume that data on species distributions are continuous observable entities (e.g. basal area, biomass), that we can model with a linear regression. We name y_iA and y_iB the observations of species y_A and y_B at site i, respectively, and we assume that the environment has a quadratic effect on both species. Therefore, using classical SDMs, model equations are:
$$ \begin{aligned} y_{iA} &= \beta^{SDM}_{0A} + \beta^{SDM}_{1A} x_i + \beta^{SDM}_{2A}x_i^2 + \epsilon_{iA}, \ \ \epsilon_{iA} \sim N(0, \sigma_A^{SDM}) \\ y_{iB} &= \beta^{SDM}_{0B} + \beta^{SDM}_{1B} x_i + \beta^{SDM}_{2B}x_i^2 + \epsilon_{iB}, \ \ \epsilon_{iB} \sim N(0, \sigma_B^{SDM}) \end{aligned} $$

Instead, our trophic model includes a dependence of y_B on its prey y_A:
$$ \begin{aligned} y_{iA} &= \beta^{tSDM}_{0A} + \beta^{tSDM}_{1A} x_i + \beta^{tSDM}_{2A}x_i^2 + \epsilon_{iA}, \ \ \epsilon_{iA} \sim N(0, \sigma_A^{tSDM}) \\ y_{iB} &= \beta^{tSDM}_{0B} + \beta^{tSDM}_{1B} x_i + \beta^{tSDM}_{2B}x_i^2 + \alpha y_{iA} + \epsilon_{iB}, \ \ \epsilon_{iB} \sim N(0, \sigma_B^{tSDM}) \end{aligned} $$

We hereafter focus on the expectation of the models (corresponding to the predicted values), focusing in particular on species y_B for which the two models are different. For SDMs:
E[y_iB| x_i] = β_0B^SDM + β_1B^SDMx_i + β_2B^SDMx_i² For trophic SDM, we need to condition on a value of y_iA. We here choose to condition on the expected value of y_iA, which corresponds to using predictions from y_A to predict y_B. In a linear model, with Gaussian data and a linear interaction between prey and predator, we have that E_A[E_B[y_B| y_A]] = E_B[y_B| E[y_A]] (which is not necessarily true in other cases).
We thus have:
E[y_iB| x_i, E[y_iA]] = β_0B^tSDM + β_1B^tSDMx_i + β_2B^tSDMx_i² + αE[y_iA] We can rewrite the expected value of y_B only as a function of the environment x only. $$ \begin{aligned} E[y_{iB} | \ x_i, E[y_{iA}] ] &= \beta^{tSDM}_{0B} + \beta^{tSDM}_{1B} x_i + \beta^{tSDM}_{2B}x_i^2 + \alpha (\beta^{tSDM}_{0A} + \beta^{tSDM}_{1A} x_i + \beta^{tSDM}_{2A}x_i^2 ) \\ & = \beta^{tSDM}_{0B} + \alpha \beta^{tSDM}_{0A} + ( \beta^{tSDM}_{1B} + \alpha \beta^{tSDM}_{1A}) x_i + (\beta^{tSDM}_{2B} + \alpha \beta^{tSDM}_{2A}) x_i^2 \\ &= \beta^{SDM}_{0B} + \beta^{SDM}_{1B} x_i + \beta^{SDM}_{2B}x_i^2 \end{aligned} $$

In this very simple case of a linear regression model with a linear interactions between the predator and the preys, the estimation of the regression coefficients is such that β_.B^SDM = β_.B^tSDM + αβ_.A^tSDM As a consequence the models estimate the same realised niche, and predictions from trophic and classical SDMs coincide.
However, our trophic model allows a separation of the direct effect of the environment on y_B (i.e., the potential niche of y_B) from its indirect effect due to the effect of y_A on y_B. In other words, using the jargon of structural equation models, β_1B^tSDM and β_2B^tSDM are the direct effects of the environment on species y_B (i.e., its potential niche), the terms αβ_1A^tSDM and αβ_2A^tSDM are the indirect effect of the environment on species y_B. The sum of these two effects, β_1B^tSDM + αβ_1A^tSDM and β_2B^tSDM + αβ_2A^tSDM corresponds to the total effect of the environment on y_B, i.e., the realised niche of y_B. In this very particular case, SDM can correctly estimate the total effect (the realised niche), but cannot separate the direct and indirect effect of the environment on y_B and thus they do not correctly estimate the potential niche.

Simulation

We hereafter present a short example confirming the analytical derivation of these results.

First, we simulate species distribution data accordingly to the directed acyclic graph presented above. Therefore, the environment has a quadratic effect on species y_A and species y_B, with a linear effect of y_A on y_B.

set.seed(1712)
# Create the environmental gradient
X = seq(0,5,length.out=100)
# Simulate data for the basal species A 
A = 3*X + 2*X^2 + rnorm(100, sd = 5)
# Simulate data for the predator species B
B = 2*X -3*X^2 + 3*A  + rnorm(100, sd = 5)

oldpar = par()
par(mfrow=c(1,3))
plot(X,A, main = "Species A realised \n and potential niche")
plot(X,B, main = "Species B realised niche")
plot(X, 2*X -3*X^2, main = "Species B potential niche")

plot of chunk simGauss


par(mfrow = oldpar$mfrow)

We now run both a classic SDM and a trophic SDM

# For species A, the two models coincide
m_A = lm(A ~ X + I(X^2))

# For species B, the two models are different:

## trophic model
m_trophic = lm(B ~ X + I(X^2) + A)
## classic SDM
m_SDM = lm(B ~ X + I(X^2))

We can now test the relationships between the coefficients of SDM and trophic SDM.

# Show the equality of coefficients
all.equal(coef(m_SDM), coef(m_trophic)[-4] + coef(m_trophic)[4]*coef(m_A))
#> [1] TRUE

# Prediction coincide
pred_SDM = m_SDM$fitted.values
pred_trophic = predict(m_trophic,
                       newdata = data.frame(X = X, A = m_A$fitted.values))

all.equal(pred_SDM, pred_trophic)
#> [1] TRUE

The coefficient of classic SDM for species y_B are equal to the environmental coefficients of trophic SDM for species y_B plus the indirect effect of the environment on species y_B (given by the product of the biotic term and the environmental coefficients of species y_A). As a consequence, the predictions of both models coincide when the predicted value of species y_A is used to predict y_B. However, unlike the traditional SDM, our trophic model allows to retrieve the correct potential niche of species y_B

# Potential niche infered by SDM 
plot(X, m_SDM$fitted.values, ylim = c(-65, 120), ylab = "B",
     main = "Potential niche of yB")

points(X, cbind(1, X, X^2) %*% coef(m_trophic)[-4], col ="red")
points(X,  2*X -3*X^2, col ="blue")
legend(-0.1, 120, legend=c("Potential niche inferred by SDM",
                           "Potential niche inferred by trophic SDM",
                           "True fundamenta niche"), pch =1,
       col=c("black","red", "blue"), cex = 0.5)

plot of chunk potential

Does this always hold?

We showed that SDM and trophic SDM estimate the same realized niche, and thus lead to the same predictions, when all the following conditions hold:

The relationship between predator and prey is linear
Species distribution are continuous abundances
The statistical model employed is a linear regression
The selected environmental variables are the same for prey and predators

However, as soon as at least one of these conditions does not hold, the estimations of the realized niche (and so model prediction) differ. We hereafter show that releasing just one of these constraints at a time is enough to obtain predictions from trophic SDM that are different, and better, than classical SDMs. Therefore, the example presented above is a very particular case that is useful to understand model similarities, but that rarely occurs in application.

Non-linear relationship between prey and predator

We hereafter release the conditions on the linearity of the relationship between y_B and y_A. Indeed, we now consider that species y_B has a quadratic response to y_A. The classic SDM does not change, as we model species as a function of the environment only. Instead, the trophic model for species y_B becomes:
y_iB = β_0B^tSDM + β_1B^tSDMx_i + β_2B^tSDMx_i² + αy_iA² + ϵ_iB, ϵ_iB ∼ N(0, σ_B^tSDM) Predictions between trophic SDM and SDM are now different, since: $$ \begin{aligned} E[y_{iB} | \ x_i, E[y_{iA}] ] &= \beta^{tSDM}_{0B} + \beta^{tSDM}_{1B} x_i + \beta^{tSDM}_{2B}x_i^2 + \alpha (\beta^{tSDM}_{0A} + \beta^{tSDM}_{1A} x_i + \beta^{tSDM}_{2A}x_i^2 )^2 \\ & = \beta^{tSDM}_{0B} + \alpha(\beta^{tSDM}_{0A})^2+ (\beta^{tSDM}_{1B} + 2\alpha\beta^{tSDM}_{0A}\beta^{tSDM}_{1A}) x_i \\ & \ \ \ \ + (\beta^{tSDM}_{2B} + \alpha (2\beta^{tSDM}_{0A}\beta^{tSDM}_{2A} + (\beta^{tSDM}_{1A})^2))x_i^2 + 2\alpha\beta^{tSDM}_{1A}\beta^{tSDM}_{2A} x_i^3 + \alpha(\beta^{tSDM}_{2A})^2x_i^4 \\ &\neq \beta^{SDM}_{0B} + \beta^{SDM}_{1B} x_i + \beta^{SDM}_{2B}x_i^2 \end{aligned} $$

In order to correctly infer the realised niche, classic SDMs should be specified with a polynomial of degree 4. As soon as we add multiple environmental variables and link between species, classic SDM should be specified with extremely complex formula, which is unfeasible for real application.

Simulation

We provide a simple example to demonstrate the analytical computations with a simulation. We use the same example of above, but we now simulate a quadratic relationship of y_B with respect to y_A.

B = 2*X -3*X^2 + 3*A^2  + rnorm(100, sd = 5)

oldpar = par()
par(mfrow=c(1,3))
plot(X,A, main = "Species A realised \n and potential niche")
plot(X,B, main = "Species B realised niche")
plot(X, 2*X -3*X^2, main = "Species B potential niche")

plot of chunk simQuadr

par(mfrow = oldpar$mfrow)

# Fit a trophic SDM
m_trophic = lm(B ~ X + I(X^2) + I(A^2))

We now compare the regression coefficients and predictions from the two models

# Show the equality of coefficients
all.equal(coef(m_SDM), coef(m_trophic)[-4] + coef(m_trophic)[4]*coef(m_A))
#> [1] "Mean relative difference: 0.185587"

# Prediction are now different
pred_SDM = m_SDM$fitted.values
pred_trophic = predict(m_trophic,
                       newdata = data.frame(X = X, A = m_A$fitted.values))

all.equal(pred_SDM, pred_trophic)
#> [1] "Mean relative difference: 53.61947"

#Compare R2
cor(B, pred_SDM)^2
#> [1] 0.8339402
cor(B, pred_trophic)^2
#> [1] 0.9605349

The relationship between the coefficients of the two models does not hold anymore more, and predictions are now different, with trophic SDM reaching a far better predictive accuracy than classic SDM for species y_B (R2 of 0.83 and 0.96 respectively).

Non-continuous observable entities (presence-absence, counts)

Let’s now consider here the case of presence-absences (i.e., binary data), without modifying the other two conditions. Therefore, we use a (generalized) linear model like logistic regression, and we model a linear relationship (at the link level) between predator and prey. Thus, for SDM we have:
$$ \begin{aligned} p(y_{iA}= 1) &= \text{logit}^{-1}(\beta^{SDM}_{0A} + \beta^{SDM}_{1A} x_i + \beta^{SDM}_{2A}x_i^2) \\ p(y_{iB} = 1) &= \text{logit}^{-1}(\beta^{SDM}_{0B} + \beta^{SDM}_{1B} x_i + \beta^{SDM}_{2B}x_i^2) \end{aligned} $$
Instead, our trophic model includes a dependence of y_B on its prey y_A:
p(y_iB = 1) = logit⁻¹(β_0B^tSDM + β_1B^tSDMx_i + β_2B^tSDMx_i² + αy_A) The equality between the models does not hold anymore due to the non linearity of the logit function. Indeed: $$ \begin{aligned} E[y_{iB} | \ x_i, E[y_{iA}] ] &= p(y_{iB} = 1 | \ x_i, p(y_{iA} = 1)) \\ &= \text{logit}^{-1}(\beta^{tSDM}_{0B} + \beta^{tSDM}_{1B} x_i + \beta^{tSDM}_{2B}x_i^2 + \alpha \text{logit}^{-1}(\beta^{tSDM}_{0A} + \beta^{tSDM}_{1A} x_i + \beta^{tSDM}_{2A}x_i^2))\\ &\neq \text{logit}^{-1}(\beta^{SDM}_{0B} + \beta^{SDM}_{1B} x_i + \beta^{SDM}_{2B}x_i^2) \end{aligned} $$

Simulations

We now confirm this analytical derivation with a simple example. The code is very similar to the previous one, but we now generate binary data from a logistic model.

# Generate the link layer basal species
V_A= 3*X - 2*X^2
# Transform to probability
prob_A = 1/(1 + exp(-V_A))
# Sample the binary data
A = rbinom(n = 100, size = 1, prob = prob_A)
# Same for species B
V_B = -3*X +1*X^2 + 3*A
prob_B = 1/(1 + exp(-V_B))
B = rbinom(n = 100, size = 1, prob = prob_B)

oldpar = par()
par(mfrow=c(1,3))
plot(X,prob_A, main = "Species A realised \n and potential niche")
plot(X,prob_B, main = "Species B realised niche")
plot(X, 1/(1+exp(-(-3*X +1*X^2))), main = "Species B potential niche")

plot of chunk simlogi


par(mfrow = oldpar$mfrow)

We then fit a classic SDM and a trophic SDM

# For species A, the two models coincide
m_A = glm(A ~ X + I(X^2), family = binomial(link = "logit"))

# For species B, the two models are different:
m_SDM = glm(B ~ X + I(X^2), family = binomial(link = "logit"))
m_trophic = glm(B ~ X + I(X^2) + A, family = binomial(link = "logit"))

We now test the relationships between model coefficients and whether the equality of model predictions still holds

# Equality of the coefficient does not hold anymore
all.equal(coef(m_SDM), coef(m_trophic)[-4] + coef(m_trophic)[4]*coef(m_A))
#> [1] "Mean relative difference: 2.305022"

# Predictions do not coincide anymore
pred_SDM = m_SDM$fitted.values
pred_trophic = predict(m_trophic,
                       newdata = data.frame(X = X, A = m_A$fitted.values))

all.equal(pred_SDM, pred_trophic)
#> [1] "Mean relative difference: 2.915665"

# Predictions are improved by trophic SDM

## AUC of species B for SDM
dismo::evaluate(p = pred_SDM[which(B==1)],
                a = pred_SDM[which(B==0)] )@auc
#> [1] 0.8408333

## AUC of species B for trophic SDM
dismo::evaluate(p = pred_trophic[which(B==1)],
                a = pred_trophic[which(B==0)] )@auc
#> [1] 0.8529167

The relationship between the coefficients of the two models does not hold anymore more, and predictions are now different, with trophic SDM that slightly improve classic SDM (AUC of species y_B goes from 0.84 to 0.85).

Non-linear statistical model

We finally show how using a non-linear statistical model leads to different predictions between classical SDMs and a trophic SDM with the example of a random forest. Since we cannot provide any analytical derivation for random forest, we simply provide an example. Data are simulated exactly as in the very first examples (continuous observable entities).

# Simulate data for the basal species A 
A = 3*X + 2*X^2 + rnorm(100, sd = 5)
# Simulate data for the predator species B
B = 2*X -3*X^2 + 3*A  + rnorm(100, sd = 5)

oldpar = par()
par(mfrow=c(1,3))
plot(X,A, main = "Species A realised \n and potential niche")
plot(X,B, main = "Species B realised niche")
plot(X, 2*X -3*X^2, main = "Species B potential niche")

plot of chunk simRF


par(mfrow = oldpar$mfrow)

We infer statistical models with random forest

library(randomForest)
#> randomForest 4.7-1.1
#> Type rfNews() to see new features/changes/bug fixes.

# Species A
m_A  = randomForest(A ~ X + I(X^2))

# Species B trophic SDM
m_trophic = randomForest(B ~ X + I(X^2) + A)
# Species B classic SDM
m_SDM <- randomForest(B ~ X + I(X^2))

We now compare the predictions from the two models.

# Prediction are now different
pred_SDM = m_SDM$predicted
pred_trophic = predict(m_trophic,
                       newdata = data.frame(X = X, A = m_A$predicted))

all.equal(pred_SDM, pred_trophic)
#> [1] "Mean relative difference: 0.08021014"

#Compare R2
cor(B, pred_SDM)^2
#> [1] 0.8075415
cor(B, pred_trophic)^2
#> [1] 0.8670482

Predictions from the model are different, and trophic SDM has a larger R2 than classical SDMs (0.86 versus 0.80).

- Compare trophic and classical SDMs: difference and equalities
- Does this always hold?