---
title: "3. Niche modeling"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{3. Niche modeling}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  out.width = "90%"
)
has_terra <- requireNamespace("terra", quietly = TRUE)
has_nicheR <- requireNamespace("nicheR", quietly = TRUE)
```

After thinning, the environmental niche is summarised by a multivariate
ellipsoid via `fit_ellipsoid()`. Two estimators are available:

- `"covmat"` — classical sample mean and covariance.
- `"mve"` — robust Minimum Volume Ellipsoid (Rousseeuw, 1985).

The ellipsoid boundary is defined by a chi-square cutoff on the squared
Mahalanobis distance, controlled by `level` (default 95 %).

```{r setup}
library(bean)
data(origin_dat_prepared, package = "bean")
data(thinned_stochastic, package = "bean")
data(thinned_deterministic, package = "bean")
env_vars <- c("bio_1", "bio_4", "bio_12", "bio_15")
```

## Fit three ellipsoids

We fit one ellipsoid to the raw prepared data and one to each of the
thinned datasets, so we can compare what bias correction does to the
inferred niche.

```{r}
origin_ellipse <- fit_ellipsoid(
  data = origin_dat_prepared,
  env_vars = env_vars,
  method = "covmat",
  level = 0.95
)

stochastic_ellipse <- fit_ellipsoid(
  data = thinned_stochastic$thinned_data,
  env_vars = env_vars,
  method = "covmat",
  level = 0.95
)

deterministic_ellipse <- fit_ellipsoid(
  data = thinned_deterministic$thinned_points,
  env_vars = env_vars,
  method = "covmat",
  level = 0.95
)

origin_ellipse
stochastic_ellipse
deterministic_ellipse
```

## Visualise the ellipsoids (2-D slices)

```{r, fig.width = 6, fig.height = 6}
plot(origin_ellipse, dims = c("bio_1", "bio_12"))
plot(stochastic_ellipse, dims = c("bio_1", "bio_12"))
plot(deterministic_ellipse, dims = c("bio_1", "bio_12"))
```

For an interactive 3-D view, install the optional package **rgl** and
call `plot(origin_ellipse, dims = c(1, 2, 3))`. If `rgl` is not
installed, the plot falls back to the 2-D view of the first two
requested variables.

## Inspecting the inferred niche

`fit_ellipsoid()` returns an object that exposes the centroid, the
covariance matrix, the precomputed inverse, the variable names, and the
subset of points classified as inside or outside the ellipsoid. These
are the building blocks downstream species-distribution-modelling
pipelines need.

```{r}
origin_ellipse$centroid
origin_ellipse$cov_matrix
origin_ellipse$dimensions
origin_ellipse$var_names
head(origin_ellipse$points_in_ellipse)
```

## Projecting suitability with nicheR

If you intend to project a `bean_ellipsoid` into geographic space,
please install the **nicheR** package and use its `predict()` method.
`bean_ellipsoid` objects carry `"nicheR_ellipsoid"` as a second S3
class, so once nicheR is attached its `predict.nicheR_ellipsoid()`
method dispatches on them automatically — no conversion step is
required. If you use the prediction step in published work, please
cite nicheR (see the References section at the bottom of this
vignette).

### Installing and loading nicheR

`nicheR` is available on CRAN:

```{r, eval = FALSE}
install.packages("nicheR")
library(nicheR)
```

### Predicting suitability

When both **nicheR** and **terra** are available, the chunk below runs
the prediction for all three ellipsoids and plots the suitability
layers side by side. When either package is missing, the chunk is
skipped automatically.

The ellipsoids were fitted on the scaled environmental variables in
`origin_dat_prepared` (the output of `prepare_bean(transform = "scale")`).
The raster must be scaled the same way before being passed to
`predict()`, otherwise the Mahalanobis distances would be computed in a
different coordinate system from the one in which the ellipsoid was
defined. We use `terra::scale()` to apply standardisation layer by
layer.

```{r, eval = has_nicheR && has_terra, fig.width = 9, fig.height = 4}
library(nicheR)
library(terra)

# Load the raster and scale each layer to mean 0, SD 1 — matching how
# the ellipsoids were trained.
env <- terra::rast(system.file("extdata", "thai_env.tif",
                               package = "bean"))
env_scaled <- terra::scale(env)

# Project each ellipsoid back into geographic space.
origin_suit <- predict(
  origin_ellipse,
  newdata             = env_scaled,
  include_suitability = TRUE,
  include_mahalanobis = FALSE
)

stochastic_suit <- predict(
  stochastic_ellipse,
  newdata             = env_scaled,
  include_suitability = TRUE,
  include_mahalanobis = FALSE
)

deterministic_suit <- predict(
  deterministic_ellipse,
  newdata             = env_scaled,
  include_suitability = TRUE,
  include_mahalanobis = FALSE
)

# Three-panel comparison: Original / Stochastic / Deterministic.
# A shared colour scale (range = [0, 1]) makes the panels directly
# comparable; the same legend applies to all three.
op <- par(mfrow = c(1, 3), mar = c(2.5, 2.5, 3, 4))
terra::plot(origin_suit[["suitability"]],
            main  = "Original (unthinned)",
            range = c(0, 1))
terra::plot(stochastic_suit[["suitability"]],
            main  = "Stochastic thinning",
            range = c(0, 1))
terra::plot(deterministic_suit[["suitability"]],
            main  = "Deterministic thinning",
            range = c(0, 1))
par(op)
```

Because the raw data is biased toward heavily sampled environments, the
unthinned model typically over-predicts those conditions. Both thinned
models produce a narrower, less inflated suitability surface — the
expected effect of bias correction.

## References

- Castaneda-Guzman, M., Hughes, C., Paansri, P., & Cobos, M. E. (2026).
  *nicheR: Ellipsoid-Based Virtual Niches and Visualization.* R package
  version 0.1.0. <https://github.com/castanedaM/nicheR>
- Rousseeuw, P. J. (1985). Multivariate estimation with high breakdown
  point. *Mathematical Statistics and Applications, Vol. B*, 283–297.
- Van Aelst, S. & Rousseeuw, P. (2009). Minimum volume ellipsoid.
  *Wiley Interdisciplinary Reviews: Computational Statistics*, 1(1),
  71–82.