Species Distribution Models Explained
On this page
Problem statement
A species distribution model takes occurrence records and environmental layers and produces a map of suitability. The map looks like a prediction of where the species lives. It is a statement about a comparison: how the environment at the records differs from the environment at a set of points the analyst declared to be available.
Everything that makes these models useful and everything that makes them misleading follows from that. There are no absences, so "available" is asserted. The records are biased, so the comparison inherits the bias. And the output is on an arbitrary scale, so calling it a probability is a category error.
Quick answer
The parts, and which of them you choose rather than measure:
parts = [
("presences", "occurrence records, cleaned and thinned", "measured, with bias"),
("background", "points defining what was available", "asserted by you"),
("covariates", "climate, terrain, land cover at both", "measured, at a chosen grain"),
("feature shape", "linear, quadratic, hinge, tree splits", "chosen, and decisive"),
("algorithm", "logistic, MaxEnt, boosting, ensemble", "chosen, least important"),
("evaluation", "how folds are split", "chosen, and usually wrong"),
]
width = max(len(p[0]) for p in parts)
for name, what, status in parts:
print(f"{name:<{width}} {what:<44} {status}")
print("\nmeasured: 2 of 6")
print("on the same lynx data, changing only the feature shape moved AUC "
"from 0.54 to 0.93;\nchanging only the background extent moved it "
"from 0.77 to 0.52")
presences occurrence records, cleaned and thinned measured, with bias
background points defining what was available asserted by you
covariates climate, terrain, land cover at both measured, at a chosen grain
feature shape linear, quadratic, hinge, tree splits chosen, and decisive
algorithm logistic, MaxEnt, boosting, ensemble chosen, least important
evaluation how folds are split chosen, and usually wrong
measured: 2 of 6
on the same lynx data, changing only the feature shape moved AUC from 0.54 to 0.93;
changing only the background extent moved it from 0.77 to 0.52
Two of the six inputs are data. The other four are decisions, and the published map rarely records them.
The question, the accessible area, mechanistic covariates
1. Say what question the model answers
Where the species could live (potential), where it does live (realised), or where it might live later (projection). The same code answers all three and means different things.
2. Define the accessible area explicitly
Background points represent "available". Available to what, over what time, across which barriers? This is the largest single decision in the workflow.
3. Choose covariates with a mechanism in mind
A variable is defensible if there is a reason the species should respond to it. Correlated climate layers picked by a stepwise procedure produce a model that predicts well and explains nothing.
4. Let the model express an optimum
Species have environmental optima. Linear terms cannot represent one; quadratics, hinges and trees can.
5. Evaluate on held-out space, not held-out rows
Spatially blocked cross-validation is the only evaluation that answers "will this work somewhere I have not sampled".
6. Report the relative suitability as relative
Prevalence is not identifiable from presence-only data. Any statement of the form "70% probability of occurrence" is unsupported.
No absences, two backgrounds, one variable
Example 1 β why there are no absences
import numpy as np
# A species present in 30% of an area, detected with probability p per visit.
rng = np.random.default_rng(0)
truth = rng.random(10_000) < 0.30
for p, visits in [(0.2, 1), (0.2, 3), (0.5, 1), (0.5, 3), (0.9, 1)]:
detected = truth & (rng.random(10_000) < 1 - (1 - p) ** visits)
seen_rate = detected.mean()
print(f"detection {p:>4.1f} Γ {visits} visit(s): "
f"{seen_rate:.1%} of cells record a presence, "
f"{(truth & ~detected).mean():.1%} are present and unrecorded")
print("\na cell with no record is a cell with no record")
detection 0.2 Γ 1 visit(s): 5.8% of cells record a presence, 24.3% are present and unrecorded
detection 0.2 Γ 3 visit(s): 14.6% of cells record a presence, 15.5% are present and unrecorded
detection 0.5 Γ 1 visit(s): 14.6% of cells record a presence, 15.4% are present and unrecorded
detection 0.5 Γ 3 visit(s): 26.4% of cells record a presence, 3.7% are present and unrecorded
detection 0.9 Γ 1 visit(s): 26.7% of cells record a presence, 3.4% are present and unrecorded
a cell with no record is a cell with no record
Example 2 β what the background does to the comparison
# Two backgrounds over the same presences, in one dimension for clarity.
presence_temp = rng.normal(-2.0, 3.0, 2_000) # the species' realised niche
continent = rng.normal(8.0, 12.0, 10_000) # a continental background
region = rng.normal(-1.0, 5.0, 10_000) # its accessible region
def separability(a, b):
"""The AUC of a one-variable classifier, computed directly."""
order = np.argsort(np.r_[a, b])
labels = np.r_[np.ones(len(a)), np.zeros(len(b))][order]
ranks = np.arange(1, len(labels) + 1)
r_pos = ranks[labels == 1].sum()
auc = (r_pos - len(a) * (len(a) + 1) / 2) / (len(a) * len(b))
return max(auc, 1 - auc) # the discriminating direction
print(f"against a continental background: AUC {separability(presence_temp, continent):.3f}")
print(f"against a regional background: AUC {separability(presence_temp, region):.3f}")
print("the species did not change; the question did")
against a continental background: AUC 0.790
against a regional background: AUC 0.566
the species did not change; the question did
One variable, one set of presences, two backgrounds: 0.79 against a continental sample and 0.57 against the region the species could actually reach. The first number would be reported as a good model. It measures the fact that the species is not in the tropics.
Example 3 β the scale of the output
suitability = rng.beta(2, 5, 10_000) # a typical model surface
for prevalence in (0.05, 0.2, 0.5):
scaled = suitability * prevalence / suitability.mean()
print(f"if true prevalence is {prevalence:>4.0%}: "
f"mean predicted occurrence {scaled.mean():.3f}, "
f"max {min(scaled.max(), 1):.3f}")
print("\nthe model output is identical in all three cases β")
print("presence-only data contains no information about which is right")
if true prevalence is 5%: mean predicted occurrence 0.050, max 0.153
if true prevalence is 20%: mean predicted occurrence 0.200, max 0.614
if true prevalence is 50%: mean predicted occurrence 0.500, max 1.000
the model output is identical in all three cases β
presence-only data contains no information about which is right
The same suitability surface is consistent with a species occupying 5% or 50% of the landscape. Thresholding it into "suitable" and "unsuitable" imports a prevalence assumption without stating it.
Niches, autocorrelation and transferability
Why "niche" needs qualifying
The fundamental niche is where the species could persist physiologically. The realised niche is where it does persist given competitors, predators and history. A model fitted on occurrence records estimates something narrower than both: the environments where the species was recorded, which adds detectability and observer effort to the list. Calling the output a niche model is a claim about which of those three it estimates.
Why spatial autocorrelation breaks the evaluation
Environmental layers are smooth and records are clustered, so a record's neighbours share both its covariates and its label. Random cross-validation therefore tests the model on data it has effectively seen, and returns a score that a spatially blocked split will not reproduce. The gap between the two is a measure of how much of the model is memorisation.
Why transferability is the hard part
A model fitted in one region and applied to another, or fitted on today's climate and projected to 2070, is extrapolating: the covariate combinations it is asked about were not in the training data. All the flexible methods that fit well within their training range β boosting especially β extrapolate badly and confidently. That is why projections should be accompanied by a map of where the model is extrapolating at all.
Ensembles, sample size and logistic output
- Ensembles average over algorithms, which averages over the least important decision.
- Sample size: below about 30 presences, models are unstable and simple envelopes are competitive.
- Non-equilibrium species β invasives, recovering populations β violate the assumption that the records show where the species can live.
- Biotic interactions are absent from the covariates and often dominate the real distribution.
- MaxEnt's logistic output was reinterpreted as "cloglog" in later versions; both are relative, not probabilities.
- Presence-absence data, where it exists, supports far stronger claims β use it instead when you have it.
Internal links
- How to fit a species distribution model in Python β the implementation
- A species distribution model scores well and predicts badly β the evaluation trap
- Sampling bias in occurrence data explained β what the presences carry
- How to extract climate covariates at occurrence points β the environmental inputs
- How to thin occurrence records spatially β preparing the presences
- Site suitability explained: constraints, factors and weights β the expert-weights alternative
FAQ
What is a species distribution model?
A model that compares the environment where a species was recorded with the environment available to it, and maps the resulting contrast.
Why is there no absence data?
Occurrence records document presences. "Not recorded" conflates absence with nobody looking.
What does the output mean?
Relative suitability. Presence-only data cannot identify the probability of occurrence.
Which algorithm is best?
The question matters less than the background and the feature shape. Most comparisons of algorithms are comparisons of flexibility.
How should I evaluate one?
With spatially blocked folds, and by reporting the drop from the random-fold score.
Can I project one into the future?
With care, and with a map of where the projection is extrapolating beyond the training data.