Using predict()
At its most basic, predict() requires two things: the
nicheR_ellipsoid object produced by
build_ellipsoid(), and newdata, the
environmental data over which predictions will be made. The
newdata must contain columns matching the variable names
used to build the ellipsoid. Only those variables are used, so
subsetting to ref_ellipse$var_names is good practice,
though the function will handle extra columns automatically.
Basic predictions to a data frame
When newdata is a data.frame,
predict() returns a data.frame with class
"nicheR_prediction". By default it includes the predictor
columns followed by two prediction columns: Mahalanobis and
suitability.
pred_df <- predict(ref_ellipse,
newdata = back_data[, ref_ellipse$var_names])
#> Starting: suitability prediction using newdata of class: data.frame...
#> Step: Using 2 predictor variables: bio_1, bio_12
#> Done: Prediction completed successfully. Returned columns: bio_1, bio_12, Mahalanobis, suitability
class(pred_df)
#> [1] "nicheR_prediction" "data.frame"
colnames(pred_df)
#> [1] "bio_1" "bio_12" "Mahalanobis" "suitability"
head(pred_df)
#> bio_1 bio_12 Mahalanobis suitability
#> 1 18.16097 680 59.71094 1.081269e-13
#> 2 18.06556 703 59.62282 1.129974e-13
#> 3 17.95946 725 59.73708 1.067230e-13
#> 4 18.01018 734 58.66810 1.821314e-13
#> 5 18.14458 748 56.37870 5.721652e-13
#> 6 18.36623 771 52.71068 3.581139e-12Basic predictions to a SpatRaster
When newdata is a SpatRaster,
predict() returns a SpatRaster where each
requested output is a named layer. The coordinate reference system,
extent, and resolution always match those of newdata. Cells
with NA in any predictor layer receive NA in
all prediction layers.
pred_rast <- predict(ref_ellipse,
newdata = ma_bios[[ref_ellipse$var_names]],
include_suitability = TRUE,
suitability_truncated = TRUE,
include_mahalanobis = TRUE,
mahalanobis_truncated = TRUE,
keep_data = TRUE)
#> Starting: suitability prediction using newdata of class: SpatRaster...
#> Step: Using 2 predictor variables: bio_1, bio_12
#> Done: Prediction completed successfully. Returned raster layers: bio_1, bio_12, Mahalanobis, suitability, Mahalanobis_trunc, suitability_trunc
pred_rast
#> class : SpatRaster
#> size : 150, 240, 6 (nrow, ncol, nlyr)
#> resolution : 0.1666667, 0.1666667 (x, y)
#> extent : -100, -60, 5, 30 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> sources : ma_bios.tif (2 layers)
#> memory (4 layers)
#> names : bio_1, bio_12, Mahalanobis, suitability, Mahal~trunc, suita~trunc
#> min values : 3.91325, 291, 0.002992, 0, 0.002992, 0
#> max values : 29.390553, 7150, 580.954687, 0.998505, 9.206902, 0.998505
names(pred_rast)
#> [1] "bio_1" "bio_12" "Mahalanobis"
#> [4] "suitability" "Mahalanobis_trunc" "suitability_trunc"Understanding the output
Looking at both outputs above, the default call always returns a Mahalanobis distance and a suitability, regardless of whether the input is a data frame or a raster. These two quantities are mathematically connected, and understanding that connection is essential for interpreting what the model is telling you.
Mahalanobis distance and the ellipsoid
The ellipsoidal niche model in nicheR represents a
species’ niche as a region in multivariate environmental space. The
ellipsoid’s shape comes from the covariance structure of the
environmental conditions the user provides, and its center, the
centroid, is the most favorable environmental
combination.
Every prediction starts with the Mahalanobis distance (\(D^2\)), which measures how far any environmental point \(\mathbf{x}\) is from the niche centroid \(\boldsymbol{\mu}\):
\[D^2 = (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu})\]
where \(\boldsymbol{\Sigma}^{-1}\) is the inverse of the covariance matrix. In simpler terms, \(D^2\) measures how unusual a set of environmental conditions is relative to the niche center, accounting for differences in variable scale and for correlations among variables. Unlike Euclidean distance, it stretches and rotates with the shape of the niche. The ellipsoid surface in E-space is exactly the set of all points at a constant Mahalanobis distance from the centroid.
A smaller \(D^2\) means environmental conditions are similar to those at the niche center. A larger \(D^2\) means they are increasingly different. Critically, \(D^2\) is not bounded between 0 and 1: it is a distance, not a probability.
From distance to suitability: the chi-square and MVN connection
nicheR converts \(D^2\)
to suitability using the multivariate normal (MVN)
kernel:
\[S = \exp\!\left(-\frac{1}{2} D^2\right)\]
This transformation is not arbitrary. The exponent in the MVN probability density function is exactly \(-\frac{1}{2} D^2\), so suitability is proportional to the MVN density at \(\mathbf{x}\). In other words, we are treating the niche as a multivariate Gaussian centered at \(\boldsymbol{\mu}\) with covariance \(\boldsymbol{\Sigma}\), and suitability is the relative likelihood of a species occurring under the conditions at \(\mathbf{x}\). This rescales the distance into a value between 0 and 1 that is more intuitive ecologically: 1 at the centroid, declining smoothly toward 0 as conditions move away from the optimum.
This connection to the MVN also determines the ellipsoid boundary.
Under the MVN assumption, \(D^2\)
follows a chi-square distribution with \(p\) degrees of freedom, where \(p\) is the number of environmental
variables. The confidence level cl stored in the ellipsoid
corresponds to a \(\chi^2_{p,\,\text{cl}}\) quantile: the
ellipsoid surface is the set of points where \(D^2 = \chi^2_{p,\,\text{cl}}\), enclosing a
fraction cl of the MVN probability mass. Setting
cl = 0.95 means the ellipsoid contains 95% of the
theoretical MVN density.
The figure below brings these ideas together. The top row shows the E-space scatter colored by \(D^2\) and the histogram of \(D^2\) values across the background. The bottom row shows the same for suitability. The vertical dashed lines in the histograms mark chi-square quantiles at three probability levels, showing exactly where the ellipsoid boundary falls relative to the distribution of background distances.
The color gradients in the scatter plots are reversed between the two metrics, reflecting their mathematical relationship: \(D^2\) increases as you move away from the centroid, so darker colors near the center indicate lower distance. Suitability decreases with distance, so brighter colors near the center indicate higher suitability. Distance measures how far a point is from the niche center; suitability transforms that distance into an ecologically interpretable value between 0 and 1.
The histograms connect this to the MVN. Under the MVN assumption,
\(D^2\) follows a chi-square
distribution, and the vertical dashed lines mark chi-square quantiles at
50%, 90%, and 99%, which correspond to ellipsoid boundaries at those
confidence levels. The same thresholds appear in the suitability
histogram, showing the equivalent suitability value at each boundary.
This makes it clear that the choice of cl when building the
ellipsoid has a direct and interpretable effect on what suitability
value marks the edge of the niche.
Additional function arguments
Beyond the basic call, predict() lets you control which
outputs are returned independently. The four output types are:
| Argument | Default | Description |
|---|---|---|
include_mahalanobis |
TRUE |
Squared Mahalanobis distance \(D^2\) for every point. Not truncated. Ranges from 0 at the centroid to infinity. |
include_suitability |
TRUE |
Suitability \(S = \exp(-0.5\, D^2)\). Not truncated. Ranges from near 0 to 1. |
mahalanobis_truncated |
FALSE |
\(D^2\) with values outside the
ellipsoid set to NA. |
suitability_truncated |
FALSE |
\(S\) with values outside the ellipsoid set to 0. |
Two additional arguments control what is returned alongside predictions:
keep_data: whether predictor columns are included in the returned object. Default isTRUEfordata.frameinput andFALSEforSpatRaster. Coordinate columns (x,y,lon,lat, etc.) are always retained whenkeep_data = TRUE.adjust_truncation_level: a number strictly between 0 and 1 that overrides theclstored in the ellipsoid for truncated outputs, without refitting the ellipsoid. See Effect of the confidence level on predictions.
All four outputs at once
Any combination of the four flags can be set to TRUE in
a single call. All requested outputs are returned together.
pred_df_all <- predict(ref_ellipse,
newdata = back_data[, ref_ellipse$var_names],
include_mahalanobis = TRUE,
include_suitability = TRUE,
mahalanobis_truncated = TRUE,
suitability_truncated = TRUE)
#> Starting: suitability prediction using newdata of class: data.frame...
#> Step: Using 2 predictor variables: bio_1, bio_12
#> Done: Prediction completed successfully. Returned columns: bio_1, bio_12, Mahalanobis, suitability, Mahalanobis_trunc, suitability_trunc
colnames(pred_df_all)
#> [1] "bio_1" "bio_12" "Mahalanobis"
#> [4] "suitability" "Mahalanobis_trunc" "suitability_trunc"To see what truncation does concretely, compare one row from outside the ellipsoid across all four columns. The non-truncated outputs always have a value; the truncated ones apply the boundary.
# A point outside the ellipsoid
outside_idx <- which(!is.na(pred_df_all$Mahalanobis) &
is.na(pred_df_all$Mahalanobis_trunc))[1]
pred_df_all[outside_idx, ]
#> bio_1 bio_12 Mahalanobis suitability Mahalanobis_trunc suitability_trunc
#> 1 18.16097 680 59.71094 1.081269e-13 NA 0The same applies to raster output:
pred_rast_all <- predict(ref_ellipse,
newdata = ma_bios[[ref_ellipse$var_names]],
include_mahalanobis = TRUE,
include_suitability = TRUE,
mahalanobis_truncated = TRUE,
suitability_truncated = TRUE)
#> Starting: suitability prediction using newdata of class: SpatRaster...
#> Step: Using 2 predictor variables: bio_1, bio_12
#> Done: Prediction completed successfully. Returned raster layers: Mahalanobis, suitability, Mahalanobis_trunc, suitability_trunc
pred_rast_all
#> class : SpatRaster
#> size : 150, 240, 4 (nrow, ncol, nlyr)
#> resolution : 0.1666667, 0.1666667 (x, y)
#> extent : -100, -60, 5, 30 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> source(s) : memory
#> names : Mahalanobis, suitability, Mahalanobis_trunc, suitability_trunc
#> min values : 0.002992, 0, 0.002992, 0
#> max values : 580.954687, 0.998505, 9.206902, 0.998505
names(pred_rast_all)
#> [1] "Mahalanobis" "suitability" "Mahalanobis_trunc"
#> [4] "suitability_trunc"The role of the confidence level and truncation
The ellipsoid boundary is defined by the chi-square cutoff \(\chi^2_{p,\,\text{cl}}\). Without
truncation, \(D^2\) and \(S\) are computed for every point in
newdata regardless of whether it falls inside or outside
the ellipsoid, and every point on Earth receives a positive suitability
value. This can be misleading when the goal is to identify where the
species could potentially occur.
Truncation operationalizes the ellipsoid boundary ecologically:
mahalanobis_truncated = TRUE: points with \(D^2 > \chi^2_{p,\,\text{cl}}\) receiveNA, meaning those environments are outside the niche and their distance is not reported.suitability_truncated = TRUE: points with \(D^2 > \chi^2_{p,\,\text{cl}}\) receive \(S = 0\), meaning unsuitable environments are explicitly scored as zero rather than receiving a small positive value.
The plots in Truncated predictions in E-space and Truncated suitability map below show truncated outputs. Outside points are rendered in grey so the full background extent is preserved in the view.
Effect of the confidence level on predictions
The choice of cl directly controls how much of the
environmental space is classified as suitable. A lower cl
produces a smaller, more conservative ellipsoid boundary; a higher
cl produces a larger, more permissive one. The ellipsoid
object itself remains unchanged, only the boundary used for truncation
shifts.
The adjust_truncation_level argument lets you explore
different thresholds without refitting the ellipsoid:
# More conservative: only the innermost 90% of the MVN distribution
pred_conservative <- predict(ref_ellipse,
newdata = back_data[, ref_ellipse$var_names],
suitability_truncated = TRUE,
adjust_truncation_level = 0.90)
# More permissive: the innermost 99%
pred_permissive <- predict(ref_ellipse,
newdata = back_data[, ref_ellipse$var_names],
suitability_truncated = TRUE,
adjust_truncation_level = 0.99)adjust_truncation_level must be a single number strictly
between 0 and 1 and only affects truncated outputs.
The figure below shows truncated suitability in E-space at three confidence levels side by side.
As cl increases from 0.90 to 0.99 the suitable area
grows larger while the ellipsoid object itself remains unchanged. The
colored interior represents the gradient of suitability within each
boundary and grey points are those classified as outside. The choice of
cl is a meaningful modeling decision that should be
justified in the context of the research question.
